انتقل إلى المحتوى الرئيسي

التقطير الكمي لكريلوف القائم على العينات لنموذج شبكة فيرميوني

تقدير الاستخدام: تسع ثوانٍ على معالج Heron r2 (ملاحظة: هذا تقدير فحسب. قد يختلف وقت التشغيل الفعلي.)

الخلفية النظرية​

يوضح هذا البرنامج التعليمي كيفية استخدام التقطير الكمي القائم على العينات (SQD) لتقدير طاقة الحالة الأساسية لنموذج شبكة فيرميوني. تحديداً، ندرس نموذج Anderson أحادي الشائبة أحادي البعد (SIAM)، الذي يُستخدم لوصف الشوائب المغناطيسية المضمّنة في المعادن.

يتبع هذا البرنامج التعليمي سير عمل مشابهاً للبرنامج التعليمي المرتبط التقطير الكمي القائم على العينات لـ Hamiltonian الكيمياء. غير أن الفارق الجوهري يكمن في طريقة بناء الدوائر الكمية. يستخدم البرنامج التعليمي الآخر نموذجاً تباينياً استدلالياً، وهو جذاب لـ Hamiltonian الكيمياء الذي قد يضم ملايين حدود التفاعل. في المقابل، يستخدم هذا البرنامج التعليمي دوائر تُقرّب التطور الزمني بواسطة الـ Hamiltonian. يمكن أن تكون هذه الدوائر عميقة، مما يجعل هذا النهج أكثر ملاءمةً لتطبيقات نماذج الشبكات. تشكّل متجهات الحالة التي تُعدّها هذه الدوائر أساساً لـ Krylov subspace، وبالتالي تتقارب الخوارزمية بصورة مثبتة وفعّالة نحو الحالة الأساسية في ظل افتراضات مناسبة.

يمكن النظر إلى النهج المستخدم في هذا البرنامج التعليمي باعتباره مزيجاً من التقنيات المستخدمة في SQD وKrylov quantum diagonalization (KQD). يُشار إلى هذا النهج المدمج أحياناً بالتقطير الكمي لكريلوف القائم على العينات (SQKD). راجع Krylov quantum diagonalization of lattice Hamiltonians للاطلاع على برنامج تعليمي حول طريقة KQD.

يستند هذا البرنامج التعليمي إلى العمل البحثي "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization"، الذي يمكن الرجوع إليه للاطلاع على مزيد من التفاصيل.

نموذج Anderson أحادي الشائبة (SIAM)​

Hamiltonian نموذج SIAM أحادي البعد هو مجموع ثلاثة حدود:

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

حيث

Himp=ε(n^d↑+n^d↓)+Un^d↑n^d↓,Hbath=−t∑j=0σ∈{↑,↓}L−1(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*}

هنا، cj,σ†/cj,σc^\dagger_{\mathbf{j},\sigma}/c_{\mathbf{j},\sigma} هي عوامل الإنشاء/الإفناء الفيرميونية لموقع الحمّام jth\mathbf{j}^{\textrm{th}} بالسبين σ\sigma، وd^σ†/d^σ\hat{d}^\dagger_{\sigma}/\hat{d}_{\sigma} هي عوامل الإنشاء/الإفناء لنمط الشائبة، وn^dσ=d^σ†d^σ\hat{n}_{d\sigma} = \hat{d}^\dagger_{\sigma} \hat{d}_{\sigma}. أما tt وUU وVV فهي أعداد حقيقية تصف تفاعلات القفز والموقع والتهجين، وε\varepsilon عدد حقيقي يحدد الجهد الكيميائي.

لاحظ أن الـ Hamiltonian هو حالة خاصة من Hamiltonian الإلكترون التفاعلي العام،

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*}

حيث يتضمن H1H_1 الحدود أحادية الجسيم التربيعية في عوامل الإنشاء والإفناء الفيرميونية، ويتضمن H2H_2 الحدود ثنائية الجسيم الرباعية. بالنسبة لـ SIAM،

H2=Un^d↑n^d↓H_2 = U \hat{n}_{d\uparrow}\hat{n}_{d\downarrow}

ويحتوي H1H_1 على بقية الحدود في الـ Hamiltonian. لتمثيل الـ Hamiltonian برمجياً، نخزّن المصفوفة hpqh_{pq} والموتر hpqrsh_{pqrs}.

أساسا الموضع والزخم​

نظراً للتماثل الانتقالي التقريبي في HbathH_\textrm{bath}، لا نتوقع أن تكون الحالة الأساسية متفرقة في أساس الموضع (الأساس المداري الذي يُحدَّد فيه الـ Hamiltonian أعلاه). لا يُضمن أداء SQD إلا إذا كانت الحالة الأساسية متفرقة، أي أنها تمتلك وزناً ملحوظاً على عدد صغير فحسب من حالات الأساس الحسابية. لتحسين تفرق الحالة الأساسية، نُجري المحاكاة في الأساس المداري الذي يكون فيه HbathH_\textrm{bath} قطرياً. نسمي هذا الأساس أساس الزخم. بما أن HbathH_\textrm{bath} هو Hamiltonian فيرميوني تربيعي، يمكن تقطيره بكفاءة عبر دوران مداري.

التطور الزمني التقريبي بواسطة الـ Hamiltonian​

لتقريب التطور الزمني بواسطة الـ Hamiltonian، نستخدم تحليل Trotter-Suzuki من الدرجة الثانية،

e−iΔtH≈e−iΔt2H2e−iΔtH1e−iΔ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}.

في إطار تحويل Jordan-Wigner، يُعادل التطور الزمني بواسطة H2H_2 بوابة CPhase واحدة بين المدارات ذات السبين المتجه لأعلى والسبين المتجه لأسفل عند موقع الشائبة. ولأن H1H_1 هو Hamiltonian فيرميوني تربيعي، يُعادل التطور الزمني بواسطة H1H_1 دوراناً مداريا.

تتشكّل حالات أساس Krylov {∣ψk⟩}k=0D−1\{ |\psi_k\rangle \}_{k=0}^{D-1}، حيث DD هو بُعد الفضاء الجزئي لـ Krylov، عبر تطبيق متكرر لخطوة Trotter واحدة، وبذلك

∣ψk⟩≈[e−iΔt2H2e−iΔtH1e−iΔ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}.

في سير عمل SQD التالي، سنأخذ عينات من هذه المجموعة من الدوائر ونعالج المجموعة المجمّعة من البتسلاسل باستخدام SQD. يتناقض هذا النهج مع النهج المستخدم في البرنامج التعليمي المرتبط التقطير الكمي القائم على العينات لـ Hamiltonian الكيمياء، حيث كانت العينات تُسحب من دائرة تباينية استدلالية واحدة.

المتطلبات​

قبل البدء في هذا البرنامج التعليمي، تأكد من تثبيت ما يلي:

  • Qiskit SDK الإصدار 1.0 أو أحدث، مع دعم التمثيل البصري
  • Qiskit Runtime الإصدار 0.22 أو أحدث (pip install qiskit-ibm-runtime)
  • SQD Qiskit addon الإصدار 0.11 أو أحدث (pip install qiskit-addon-sqd)
  • ffsim (pip install ffsim)

الخطوة 1: تعيين المسألة إلى دائرة كمية​

أولاً، نُولّد Hamiltonian نموذج SIAM في أساس الموضع. يُمثَّل الـ Hamiltonian بالمصفوفة hpqh_{pq} والموتر hpqrsh_{pqrs}. ثم نُحوّله إلى أساس الزخم. في أساس الموضع، نضع الشائبة في الموقع الأول. غير أننا عند التحويل إلى أساس الزخم، ننقل الشائبة إلى موقع مركزي لتيسير التفاعلات مع المدارات الأخرى.

# 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

بعد ذلك، نُولّد الدوائر لإنتاج حالات أساس Krylov. بالنسبة لكل نوع من أنواع السبين، تُعطى الحالة الابتدائية ∣ψ0⟩\ket{\psi_0} بوصفها تراكباً لجميع الإثارات الممكنة للإلكترونات الثلاثة الأقرب إلى مستوى Fermi في أقرب 4 أنماط فارغة بدءاً من الحالة ∣00⋯0011⋯11⟩|00\cdots 0011 \cdots 11\rangle، ويتحقق ذلك بتطبيق سبع بوابات XXPlusYYGate. تُنتج الحالات المتطورة زمنياً عبر تطبيقات متتالية لخطوة Trotter من الدرجة الثانية.

للاطلاع على وصف أكثر تفصيلاً لهذا النموذج وكيفية تصميم الدوائر، يُرجى الرجوع إلى "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

الخطوة 2: تحسين المسألة للتنفيذ الكمي​

بعد إنشاء الدوائر، يمكننا تحسينها لجهاز مستهدف. نختار وحدة المعالجة الكمية (QPU) الأقل ازدحاماً والتي تحتوي على 127 كيوبت على الأقل. اطلع على وثائق Qiskit IBM® Runtime للمزيد من المعلومات.

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

الآن، نستخدم Qiskit لتحويل الدوائر إلى الصيغة الملائمة للواجهة الخلفية المستهدفة.

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)

الخطوة 3: التنفيذ باستخدام Qiskit primitives​

بعد تحسين الدوائر للتنفيذ على الأجهزة، أصبحنا مستعدين لتشغيلها على الجهاز المستهدف وجمع العينات لتقدير طاقة الحالة الأساسية. بعد استخدام الأداة الأولية Sampler لأخذ عينات من البتسلاسل في كل دائرة، نجمع جميع النتائج في قاموس عدّ واحد ونرسم رسماً بيانياً لأكثر 20 بتسلسلة عينةً شيوعاً.

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

الخطوة 4: المعالجة اللاحقة وإعادة النتيجة إلى الصيغة الكلاسيكية المطلوبة​

الآن، نُشغّل خوارزمية SQD باستخدام الدالة diagonalize_fermionic_hamiltonian. راجع توثيق API للاطلاع على شرح وسيطات هذه الدالة.

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

تُرسم نتائج الكود التالي في مخططين. يُظهر المخطط الأول الطاقة المحسوبة كدالة لعدد تكرارات استرداد التهيئة، ويُظهر المخطط الثاني متوسط إشغال كل مدار مكاني بعد التكرار الأخير. كطاقة مرجعية، نستخدم نتائج حساب DMRG أُجري بصورة منفصلة.

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

التحقق من الطاقة​

الطاقة التي تُعيدها SQD مضمونة بوصفها حداً أعلى لطاقة الحالة الأساسية الحقيقية. يمكن التحقق من قيمة الطاقة لأن SQD تُعيد أيضاً معاملات متجه الحالة الذي يُقرّب الحالة الأساسية. يمكنك حساب الطاقة من متجه الحالة باستخدام مصفوفتي كثافة الاختزال أحادية وثنائية الجسيم، كما هو موضح في كتلة الكود التالية.

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

المراجع​