Zum Haubtinhald schringsn

Sample-basierte Krylov-Quantediagonalisieerung von'm fermionischn Giddermodäll

Schätzung vom Zeitaufwand: Neun Segunden uff'm Heron-r2-Brozässor (ACHTUNG: Des is nur 'ne Schätzung. Dei Laufzeit gann abwäschn.)

Hintergrund

Deer Tutorial zeechd, wie de sample-basierte Quantediagonalisieerung (SQD) nutzen kannst, um die Grundsuschstandsenergie von'm fermionischn Giddermodäll ze schätzn. Genauer gschaut studie'rn mer 's eindimensionale Einfach-Störschdeeln-Anderson-Modäll (SIAM), das ze'm Beschrääbn von magnetischn Störschdeeln in Metallen gebracht wird.

Deer Tutorial folcd 'nm ähnlichn Abloof wie de verwoandde Tutorial Sample-basierte Quantediagonalisieerung von'm Chemie-Hamiltonian. A a wichtschn Unterschied gibbts ba de Art, wie die Quanteschaltgreese gebaut wärn. De ander Tutorial nutzd'n heuristischn Variationsansatz, was bein Chemie-Hamiltonians mit möglicherweise Millionen von Wächselwirkungstermen schee is. Deer Tutorial hingeen nutzd Schaltgreese, die die Zeitentwiggelung durschn Hamiltonian annähern. Solche Schaltgreese kenn recht dief sei, was denn Oosotz fier Anwendungen uff Giddermodälle bässer macht. Die Suschstandsvektorn, was diese Schaltgreese vorberaadn, billn die Grundlache fier 'n Krylov-Unterroum, un doorum konvergierd der Algorithmus nachwäslich un effizient zum Grundsuschstand, wenn geeignete Voraussetzungen erfüllt sinn.

Dr Oosotz, der in deem Tutorial verwendet wird, kann als Gombination von den Techniggen von SQD un Krylov-Quantediagonalisieerung (KQD) oangschaut wärn. Die kombinierte Methode wird manchmool als sample-basierte Krylov-Quantediagonalisieerung (SQKD) bezeichned. Schau mer beim Tutorial Krylov-Quantediagonalisieerung von Gidder-Hamiltonians fier 'ne Eiführung in die KQD-Methode.

Deer Tutorial basiert uff der Arweet "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization", wo de fier mehr Eizelheitn nachschaun kannst.

Einfach-Störschdeeln-Anderson-Modäll (SIAM)

Dr eindimensionale SIAM-Hamiltonian is 'ne Summe von drei Tärmen:

H=Himp+Hbath+Hhyb,H = H_{\textrm{imp}}+ H_\textrm{bath} + H_\textrm{hyb},

wo

Himp=ε(n^d+n^d)+Un^dn^d,Hbath=tj=0σ{,}L1(c^j,σc^j+1,σ+c^j+1,σc^j,σ),Hhyb=Vσ{,}(d^σc^0,σ+c^0,σd^σ).\begin{align*} H_\textrm{imp} &= \varepsilon \left( \hat{n}_{d\uparrow} + \hat{n}_{d\downarrow} \right) + U \hat{n}_{d\uparrow}\hat{n}_{d\downarrow}, \\ H_\textrm{bath} &= -t \sum_{\substack{\mathbf{j} = 0\\ \sigma\in \{\uparrow, \downarrow\}}}^{L-1} \left(\hat{c}^\dagger_{\mathbf{j}, \sigma}\hat{c}_{\mathbf{j}+1, \sigma} + \hat{c}^\dagger_{\mathbf{j}+1, \sigma}\hat{c}_{\mathbf{j}, \sigma} \right), \\ H_\textrm{hyb} &= V\sum_{\sigma \in \{\uparrow, \downarrow \}} \left(\hat{d}^\dagger_\sigma \hat{c}_{0, \sigma} + \hat{c}^\dagger_{0, \sigma} \hat{d}_{\sigma} \right). \end{align*}

Hierbei sinn cj,σ/cj,σc^\dagger_{\mathbf{j},\sigma}/c_{\mathbf{j},\sigma} die fermionischn Erzeuchungs-/Vernichtungsoperatorn fier de jth\mathbf{j}^{\textrm{th}} Badschdeelle mit Schdrehimpuls σ\sigma, d^σ/d^σ\hat{d}^\dagger_{\sigma}/\hat{d}_{\sigma} sinn Erzeuchungs-/Vernichtungsoperatorn fier de Störschdeellenmode, un n^dσ=d^σd^σ\hat{n}_{d\sigma} = \hat{d}^\dagger_{\sigma} \hat{d}_{\sigma}. tt, UU un VV sinn reelle Zahln, die Hüpfen, Ortswächselwirkung un Hybridisieerung beschrääbn, un ε\varepsilon is 'ne reelle Zahl, die das chemische Potenzial feschdlächt.

Beachde, dass der Hamiltonian 'ne spezifische Instanz vom allgemäänen Wächselwirkungselektron-Hamiltonian is,

H=p,qσhpqa^pσa^qσ+p,q,r,sστhpqrs2a^pσa^qτa^sτa^rσ=H1+H2,\begin{align*} H &= \sum_{\substack{p, q \\ \sigma}} h_{pq} \hat{a}^\dagger_{p\sigma} \hat{a}_{q\sigma} + \sum_{\substack{p, q, r, s \\ \sigma \tau}} \frac{h_{pqrs}}{2} \hat{a}^\dagger_{p\sigma} \hat{a}^\dagger_{q\tau} \hat{a}_{s\tau} \hat{a}_{r\sigma} \\ &= H_1 + H_2, \end{align*}

wo H1H_1 aus Eeinkörpertärmen beschdäht, die quadratisch in dn fermionischn Erzeuchungs- un Vernichtungsoperatoorn sinn, un H2H_2 aus Zweikörpertärmen, die quartisch sinn. Fier's SIAM gillt

H2=Un^dn^dH_2 = U \hat{n}_{d\uparrow}\hat{n}_{d\downarrow}

un H1H_1 enthält de reschtlichn Tärme im Hamiltonian. Um'n Hamiltonian brogrammatiisch darstellen ze kenne, späächrn mer die Matrix hpqh_{pq} un'n Tensoor hpqrsh_{pqrs}.

Orts- un Impulsbasis

Wechn der näherungsweesn Translatschienssiemmetrie in HbathH_\textrm{bath} erwaardn mer nedd, dass der Grundsuschstand in der Ortsbasis (der Orbitalbasis, in der der Hamiltonian oangegääbn is) schütterschd is. Die Leistung von SQD is nur garantiert, wenn der Grundsuschstand schüttersch is, also nur uff wenige Berechnungsbasissuschstände erhebliches Gewicht hat. Um die Schütterschgeet vom Grundsuschstand ze värbässern, führn mer die Simulatschn in der Orbitalbasis dursch, in der HbathH_\textrm{bath} diagonal is. Mir nenn diese Basis die Impulsbasis. Weil HbathH_\textrm{bath} 'n quadratischr fermionischr Hamiltonian is, gann er effizient dursch 'ne Orbitalrotatschn diagonalisierd wärn.

Näherungsweise Zeitentwiggelung durschn Hamiltonian

Um die Zeitentwiggelung durschn Hamiltonian ze näherungswäse ze beschrääbn, verwändn mer 'ne Trotter-Suzuki-Zerlegung zweiter Ordnung,

eiΔtHeiΔt2H2eiΔtH1eiΔt2H2. e^{-i \Delta t H} \approx e^{-i\frac{\Delta t}{2} H_2} e^{-i\Delta t H_1} e^{-i\frac{\Delta t}{2} H_2}.

Untr der Jordan-Wigner-Transformatschn entspricht die Zeitentwiggelung dursch H2H_2 'm eenzelnen CPhase-Gate zwischen dn Spin-rauf- un Spin-runner-Orbitalen uff der Störschdeellenstelle. Weil H1H_1 'n quadratischr fermionischr Hamiltonian is, entspricht die Zeitentwiggelung dursch H1H_1 'ner Orbitalrotatschn.

Die Krylov-Basissuschstände {ψk}k=0D1\{ |\psi_k\rangle \}_{k=0}^{D-1}, wo DD die Dimensschn vom Krylov-Unterroum is, wärn dursch wiederhooltes Anwändn von 'm eenzelnen Trotter-Schredd gebilded, also

ψk[eiΔt2H2eiΔtH1eiΔt2H2]kψ0. |\psi_k\rangle \approx \left[e^{-i\frac{\Delta t}{2} H_2} e^{-i\Delta t H_1} e^{-i\frac{\Delta t}{2} H_2} \right]^k\ket{\psi_0}.

Im nachfolgenden SQD-basieerdn Workflow werden mer von deem Schdroum von Schaltgreesn Samples zieen un die gombinierte Menge von Bitfolgen mit SQD noochberarbeiten. Deer Oosotz unterschäädet sich von'm oodern Tutorial Sample-basierte Quantediagonalisieerung von'm Chemie-Hamiltonian, wo Samples von 'nem eenzischn heuristischn Variationsschaltgrees gezochn worn sinn.

Vorraussetzungen

Bevor de mit deem Tutorial oofoongst, schau druff, dass folchndes installierd is:

  • Qiskit SDK v1.0 oder neeschr, mit Visualisierungs-Unterstützung
  • Qiskit Runtime v0.22 oder neeschr (pip install qiskit-ibm-runtime)
  • SQD Qiskit Addon v0.11 oder neeschr (pip install qiskit-addon-sqd)
  • ffsim (pip install ffsim)

Schdregg 1: Probläm uff 'n Quanteschaltgrees abbildn

Erschd generieern mer'n SIAM-Hamiltonian in der Ortsbasis. Dr Hamiltonian wird dursch die Matrix hpqh_{pq} un'n Tensoor hpqrsh_{pqrs} dargestellt. Dann rotieern mer ihn zur Impulsbasis. In der Ortsbasis lächn mer die Störschdeelle uff die erschde Schdelle. Wenn mer aber zur Impulsbasis wechseln, versiehbn mer die Störschdeelle uff 'ne zendraaln Schdelle, um Wächselwirkungen mit oodrn Orbitalen ze erleichtern.

# Added by doQumentation — installs packages not in the Binder environment
!pip install -q ffsim qiskit-addon-sqd
import numpy as np

def siam_hamiltonian(
norb: int,
hopping: float,
onsite: float,
hybridization: float,
chemical_potential: float,
) -> tuple[np.ndarray, np.ndarray]:
"""Hamiltonian for the single-impurity Anderson model."""
# Place the impurity on the first site
impurity_orb = 0

# One body matrix elements in the "position" basis
h1e = np.zeros((norb, norb))
np.fill_diagonal(h1e[:, 1:], -hopping)
np.fill_diagonal(h1e[1:, :], -hopping)
h1e[impurity_orb, impurity_orb + 1] = -hybridization
h1e[impurity_orb + 1, impurity_orb] = -hybridization
h1e[impurity_orb, impurity_orb] = chemical_potential

# Two body matrix elements in the "position" basis
h2e = np.zeros((norb, norb, norb, norb))
h2e[impurity_orb, impurity_orb, impurity_orb, impurity_orb] = onsite

return h1e, h2e

def momentum_basis(norb: int) -> np.ndarray:
"""Get the orbital rotation to change from the position to the momentum basis."""
n_bath = norb - 1

# Orbital rotation that diagonalizes the bath (non-interacting system)
hopping_matrix = np.zeros((n_bath, n_bath))
np.fill_diagonal(hopping_matrix[:, 1:], -1)
np.fill_diagonal(hopping_matrix[1:, :], -1)
_, vecs = np.linalg.eigh(hopping_matrix)

# Expand to include impurity
orbital_rotation = np.zeros((norb, norb))
# Impurity is on the first site
orbital_rotation[0, 0] = 1
orbital_rotation[1:, 1:] = vecs

# Move the impurity to the center
new_index = n_bath // 2
perm = np.r_[1 : (new_index + 1), 0, (new_index + 1) : norb]
orbital_rotation = orbital_rotation[:, perm]

return orbital_rotation

def rotated(
h1e: np.ndarray, h2e: np.ndarray, orbital_rotation: np.ndarray
) -> tuple[np.ndarray, np.ndarray]:
"""Rotate the orbital basis of a Hamiltonian."""
h1e_rotated = np.einsum(
"ab,Aa,Bb->AB",
h1e,
orbital_rotation,
orbital_rotation.conj(),
optimize="greedy",
)
h2e_rotated = np.einsum(
"abcd,Aa,Bb,Cc,Dd->ABCD",
h2e,
orbital_rotation,
orbital_rotation.conj(),
orbital_rotation,
orbital_rotation.conj(),
optimize="greedy",
)
return h1e_rotated, h2e_rotated

# Total number of spatial orbitals, including the bath sites and the impurity
# This should be an even number
norb = 20

# System is half-filled
nelec = (norb // 2, norb // 2)
# One orbital is the impurity, the rest are bath sites
n_bath = norb - 1

# Hamiltonian parameters
hybridization = 1.0
hopping = 1.0
onsite = 10.0
chemical_potential = -0.5 * onsite

# Generate Hamiltonian in position basis
h1e, h2e = siam_hamiltonian(
norb=norb,
hopping=hopping,
onsite=onsite,
hybridization=hybridization,
chemical_potential=chemical_potential,
)

# Rotate to momentum basis
orbital_rotation = momentum_basis(norb)
h1e_momentum, h2e_momentum = rotated(h1e, h2e, orbital_rotation.T.conj())
# In the momentum basis, the impurity is placed in the center
impurity_index = n_bath // 2

Als nächstes generieern mer die Schaltgreese, um die Krylov-Basissuschstände herzustellen. Fier jede Schbezies von Schbins is der Ausgangsuschstand ψ0\ket{\psi_0} als Überlagerung aller möglichn Anregungen von den drei Elektronen, die 'm Fermi-Niveau am nächstn sinn, in die 4 nächstn leern Modes ab dem Zuschstand 00001111|00\cdots 0011 \cdots 11\rangle gääbn, un mit siebm XXPlusYYGates realisiert. Die zeidentwickeltn Suschstände wärn dursch sukzessive Anwendung von 'm Trotter-Schredd zweiter Ordnung erzeugt.

Fier 'ne ausführlichere Beschräbung von deem Modäll un wie die Schaltgreese entworfen sinn, schau noo "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization".

from typing import Sequence

import ffsim
import scipy
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit import CircuitInstruction, Qubit
from qiskit.circuit.library import CPhaseGate, XGate, XXPlusYYGate

def prepare_initial_state(qubits: Sequence[Qubit], norb: int, nocc: int):
"""Prepare initial state."""
x_gate = XGate()
rot = XXPlusYYGate(0.5 * np.pi, -0.5 * np.pi)
for i in range(nocc):
yield CircuitInstruction(x_gate, [qubits[i]])
yield CircuitInstruction(x_gate, [qubits[norb + i]])
for i in range(3):
for j in range(nocc - i - 1, nocc + i, 2):
yield CircuitInstruction(rot, [qubits[j], qubits[j + 1]])
yield CircuitInstruction(
rot, [qubits[norb + j], qubits[norb + j + 1]]
)
yield CircuitInstruction(rot, [qubits[j + 1], qubits[j + 2]])
yield CircuitInstruction(
rot, [qubits[norb + j + 1], qubits[norb + j + 2]]
)

def trotter_step(
qubits: Sequence[Qubit],
time_step: float,
one_body_evolution: np.ndarray,
h2e: np.ndarray,
impurity_index: int,
norb: int,
):
"""A Trotter step."""
# Assume the two-body interaction is just the on-site interaction of the impurity
onsite = h2e[
impurity_index, impurity_index, impurity_index, impurity_index
]
# Two-body evolution for half the time
yield CircuitInstruction(
CPhaseGate(-0.5 * time_step * onsite),
[qubits[impurity_index], qubits[norb + impurity_index]],
)
# One-body evolution for the full time
yield CircuitInstruction(
ffsim.qiskit.OrbitalRotationJW(norb, one_body_evolution), qubits
)
# Two-body evolution for half the time
yield CircuitInstruction(
CPhaseGate(-0.5 * time_step * onsite),
[qubits[impurity_index], qubits[norb + impurity_index]],
)

# Time step
time_step = 0.2
# Number of Krylov basis states
krylov_dim = 8

# Initialize circuit
qubits = QuantumRegister(2 * norb, name="q")
circuit = QuantumCircuit(qubits)

# Generate initial state
for instruction in prepare_initial_state(qubits, norb=norb, nocc=norb // 2):
circuit.append(instruction)
circuit.measure_all()

# Create list of circuits, starting with the initial state circuit
circuits = [circuit.copy()]

# Add time evolution circuits to the list
one_body_evolution = scipy.linalg.expm(-1j * time_step * h1e_momentum)
for i in range(krylov_dim - 1):
# Remove measurements
circuit.remove_final_measurements()
# Append another Trotter step
for instruction in trotter_step(
qubits,
time_step,
one_body_evolution,
h2e_momentum,
impurity_index,
norb,
):
circuit.append(instruction)
# Measure qubits
circuit.measure_all()
# Add a copy of the circuit to the list
circuits.append(circuit.copy())
circuits[0].draw("mpl", scale=0.4, fold=-1)

Output of the previous code cell

circuits[-1].draw("mpl", scale=0.4, fold=-1)

Output of the previous code cell

Schdregg 2: Probläm fier die Quanteausführung optimiern

Jetzt, wo mer die Schaltgreese erstellt ham, kenne mer se fier 'n Zielhardware optimiern. Mir nähmn'n am wenischstn ausgelastndn QPU mit mindeschdens 127 Qubits. Schau in die Qiskit IBM® Runtime-Doks fier mehr Infos.

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
backend = service.least_busy(
operational=True, simulator=False, min_num_qubits=127
)
print(f"Using backend {backend.name}")
Using backend ibm_fez

Jetzt nutzn mer Qiskit, um die Schaltgreese fier'n Ziel-Backend ze transpiließrn.

from qiskit.transpiler import generate_preset_pass_manager

pass_manager = generate_preset_pass_manager(
optimization_level=3, backend=backend
)
isa_circuits = pass_manager.run(circuits)

Schdregg 3: Mit Qiskit-Primitives ausführn

Nooch'm Optimiern der Schaltgreese fier die Hardwareausführung sinn mer bereet, se uff der Zielhardware ze laafen ze lassn un Samples fier die Grundsuschstandsenergischätzung ze sammelln. Nooch'm Nutzen vom Sampler-Primitive zum Samplern von Bitfolgen aus jedn Schaltgrees kombiniern mer alle Ärgebnisse in 'n eenzschn Counts-Dictionary un plottn die 20 am häufichsdn gesampletn Bitfolgen.

from qiskit.visualization import plot_histogram
from qiskit_ibm_runtime import SamplerV2 as Sampler

# Sample from the circuits
sampler = Sampler(backend)
job = sampler.run(isa_circuits, shots=500)
from qiskit.primitives import BitArray

# Combine the counts from the individual Trotter circuits
bit_array = BitArray.concatenate_shots(
[result.data.meas for result in job.result()]
)

plot_histogram(bit_array.get_counts(), number_to_keep=20)

Output of the previous code cell

Schdregg 4: Noochbearbeitung un Rückgaab von Ärgebnissn im gwünschtn glassischn Format

Jetzt laafen mer'n SQD-Algorithmus mit der diagonalize_fermionic_hamiltonian-Funktschen. Schau in die API-Dokumentatschn fier Erklärungen zu de Argumände von deeser Funktschen.

from qiskit_addon_sqd.fermion import (
SCIResult,
diagonalize_fermionic_hamiltonian,
)

# List to capture intermediate results
result_history = []

def callback(results: list[SCIResult]):
result_history.append(results)
iteration = len(result_history)
print(f"Iteration {iteration}")
for i, result in enumerate(results):
print(f"\tSubsample {i}")
print(f"\t\tEnergy: {result.energy}")
print(
f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
)

rng = np.random.default_rng(24)
result = diagonalize_fermionic_hamiltonian(
h1e_momentum,
h2e_momentum,
bit_array,
samples_per_batch=100,
norb=norb,
nelec=nelec,
num_batches=3,
max_iterations=5,
symmetrize_spin=True,
callback=callback,
seed=rng,
)
Iteration 1
Subsample 0
Energy: -28.61321893815165
Subspace dimension: 10609
Subsample 1
Energy: -28.628985564542244
Subspace dimension: 13924
Subsample 2
Energy: -28.620151775558114
Subspace dimension: 10404
Iteration 2
Subsample 0
Energy: -28.656893066053115
Subspace dimension: 34225
Subsample 1
Energy: -28.65277622004119
Subspace dimension: 38416
Subsample 2
Energy: -28.670856034959165
Subspace dimension: 39601
Iteration 3
Subsample 0
Energy: -28.684787675404362
Subspace dimension: 42436
Subsample 1
Energy: -28.676984757118426
Subspace dimension: 50176
Subsample 2
Energy: -28.671581704249885
Subspace dimension: 40804
Iteration 4
Subsample 0
Energy: -28.6859683054753
Subspace dimension: 47961
Subsample 1
Energy: -28.69418206537316
Subspace dimension: 51529
Subsample 2
Energy: -28.686083516445752
Subspace dimension: 51529
Iteration 5
Subsample 0
Energy: -28.694665630711178
Subspace dimension: 50625
Subsample 1
Energy: -28.69505984237118
Subspace dimension: 47524
Subsample 2
Energy: -28.6942873883992
Subspace dimension: 48841

De nachfolgend Codezelle plotd die Ärgebnisse. Dr erschde Bläd zeechd die berechnete Energie als Funktschen von der Anzool von Konfiguratschensweederhärstellungsieteratschen, un dr zweede Bläd zeechd die durchschnittliche Besetzung von jedn räumlichn Orbital nooch der letztn Iteratschn. Als Referenzenergie nutzen mer die Ärgebnisse von 'ner DMRG-Berechnungch, die getrenndt durchgführd worn is.

import matplotlib.pyplot as plt

dmrg_energy = -28.70659686

min_es = [
min(result, key=lambda res: res.energy).energy
for result in result_history
]
min_id, min_e = min(enumerate(min_es), key=lambda x: x[1])

# Data for energies plot
x1 = range(len(result_history))

# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, min_es, label="energy", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].axhline(
y=dmrg_energy, color="#BF5700", linestyle="--", label="DMRG energy"
)
axs[0].set_title("Approximated Ground State Energy vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy", fontdict={"fontsize": 12})
axs[0].legend()

# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})

print(f"Reference (DMRG) energy: {dmrg_energy:.5f}")
print(f"SQD energy: {min_e:.5f}")
print(f"Absolute error: {abs(min_e - dmrg_energy):.5f}")
plt.tight_layout()
plt.show()
Reference (DMRG) energy: -28.70660
SQD energy: -28.69506
Absolute error: 0.01154

Output of the previous code cell

Überprüfung von der Energie

Die von SQD zurückgegääbne Energie is garantiert 'ne obere Schranke fier die waare Grundsuschstandsenergie. Dr Weärt von der Energie gann überprüfd wärn, weil SQD ooch die Koeffiziendn vom Suschstandsvektoor, der'n Grundsuschstand näherungsweese beschreibt, zurückgibd. De kannst die Energie aus'm Suschstandsvektoor mit seinen 1- un 2-Deilchn-Dichtematrizen berechnern, wie in der nachfolgenden Codezelle gezeechd.

rdm1 = result.sci_state.rdm(rank=1, spin_summed=True)
rdm2 = result.sci_state.rdm(rank=2, spin_summed=True)

energy = np.sum(h1e_momentum * rdm1) + 0.5 * np.sum(h2e_momentum * rdm2)

print(f"Recomputed energy: {energy:.5f}")
Recomputed energy: -28.69506

Quelln