오비탈 최적화로 SQD 추정값 개선하기

Sample-based quantum diagonalization(SQD)은 전자 구성의 고정된 부분 공간에서 Hamiltonian을 대각화하여 바닥 상태 에너지를 근사합니다. 이 추정값은 Hamiltonian을 표현하는 오비탈 기저에 따라 달라지며, 오비탈 최적화(OO)는 이러한 자유도를 활용하여 부분 공간을 확대하지 않고도 에너지를 낮춥니다.

이 안내서는 $N_2$ 분자에 SQD를 실행한 다음 오비탈 최적화로 결과를 개선합니다. Hamiltonian을 표현하고 에너지를 최소화하는 오비탈 회전을 찾는 데 ffsim을 사용합니다.

SQD 실행하기

분자 오비탈(MO) 기저에서 $N_2$ 의 분자 적분을 구성하고, 균일한 무작위 샘플을 생성한 뒤, SQD를 실행하여 바닥 상태 근사를 얻습니다.

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"ffsim": "ffsim", "numpy": "numpy", "pyscf": "pyscf", "qiskit_addon_sqd": "qiskit-addon-sqd"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

import numpy as np
import pyscf
import pyscf.cc
import pyscf.mcscf
from qiskit_addon_sqd.counts import generate_bit_array_uniform
from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian

# Specify molecule properties
num_orbitals = 16
num_elec_a = num_elec_b = 5
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
    basis="6-31g",
    symmetry="Dooh",
)

# Define active space
n_frozen = 2
active_space = range(n_frozen, mol.nao_nr())

# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
num_orbitals = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
num_elec_a = (n_electrons + mol.spin) // 2
num_elec_b = (n_electrons - mol.spin) // 2
cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)

# Compute exact energy
exact_energy = cas.run().e_tot

# Create a seed to control randomness throughout this workflow
rng = np.random.default_rng(24)

# Generate random samples
bit_array = generate_bit_array_uniform(
    10_000, num_orbitals * 2, rand_seed=rng
)

# Run SQD
result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=100,
    norb=num_orbitals,
    nelec=(num_elec_a, num_elec_b),
    num_batches=1,
    max_iterations=5,
    symmetrize_spin=True,
    seed=rng,
)

converged SCF energy = -108.835236570775
CASCI E = -109.046671778080  E(CI) = -32.8155692383187  S^2 = 0.0000000

sqd_energy = result.energy + nuclear_repulsion_energy
print(f"Exact energy:  {exact_energy:.8f}")
print(f"SQD energy:    {sqd_energy:.8f}")

Exact energy:  -109.04667178
SQD energy:    -108.98469255

오비탈 최적화하기

오비탈 최적화는 변분 에너지를 낮추는 오비탈 회전을 찾습니다.

E = \langle \psi | \mathcal{U}^\dagger\, H\, \mathcal{U} | \psi \rangle

이는 SQD 바닥 상태 근사 $|\psi\rangle$ 의 에너지입니다. 오비탈 회전은 $N \times N$ 유니터리 행렬 $\mathbf{U}$ ( $N$ 은 공간 오비탈의 수)로 지정되며, 이 행렬은 다음 연산자를 통해 다체 상태에 작용합니다.

\mathcal{U} = \exp\left[\sum_{pq, \sigma} \log(\mathbf{U})_{pq}\, a^\dagger_{p\sigma} a_{q\sigma}\right].

ffsim.optimize_orbitals는 행렬 $\mathbf{U}$ 를 반환하며, 이를 오비탈 기저에 적용하는 것(hamiltonian.rotated 사용)은 상태에 $\mathcal{U}$ 를 적용하는 것과 동등합니다. 자세한 내용은 ffsim 오비탈 회전 설명을 참고하세요.

오비탈을 회전하면 부분 공간에서 보는 Hamiltonian이 바뀌므로, 에너지가 더 이상 개선되지 않을 때까지 다음 두 단계를 번갈아 수행합니다:

고정된 구성 집합에 대해 현재 기저에서 Hamiltonian을 대각화합니다.
결과 상태의 에너지를 최소화하는 회전을 찾아 오비탈을 최적화한 다음, 적분을 새 기저로 회전합니다.

오비탈 회전 단계는 상태의 1체 및 2체 축약 밀도 행렬(RDM)로부터 에너지를 최소화하는 회전을 찾는 ffsim.optimize_orbitals에 맡깁니다. 자세한 내용은 Sec. II A 4를 참고하세요.

여기서 오비탈 최적화가 도움이 되는 이유

SCF 분자 오비탈(MO) 기저는 완전 CI 문제에 대해서는 오비탈 회전에 대해 정류(stationary) 상태입니다. 그러나 SQD는 작게 잘린 부분 공간(여기서는 약 1,900만 개의 완전 CI 행렬식 중 수백 개의 CI 문자열)에서 작동하며, 이 경우 MO 기저는 일반적으로 최적이 아니므로 오비탈을 회전하면 부분 공간이 표현할 수 있는 에너지가 낮아집니다.

import ffsim
from pyscf import fci

# ffsim's ``MolecularHamiltonian`` uses the same "chemist" ordering for the two-body
# tensor as PySCF's ``eri``, and stores the nuclear repulsion energy as the constant
# term so that expectation values come out as total energies.

대각화와 오비탈 최적화를 번갈아 수행하기

각 반복에서 오비탈 회전의 효과만 분리해서 볼 수 있도록, 대각화 부분 공간을 위에서 SQD가 찾은 구성으로 고정합니다. 각 반복은 다음을 수행합니다:

PySCF의 선택적 CI 솔버를 사용하여 현재 기저에서 고정된 부분 공간에 대해 Hamiltonian을 대각화합니다.
결과 상태의 RDM을 구성합니다. ffsim.optimize_orbitals에 필요한 것은 이것뿐입니다.
오비탈을 최적화합니다. ffsim.optimize_orbitals가 에너지를 최소화하는 회전을 반환하며, 이를 적분에 적용하여 개선된 기저로 이동합니다.

각 최적화 단계 이전의 에너지를 기록합니다. 반복마다 기저가 개선되므로, 이 수열은 고정된 부분 공간에서 얻을 수 있는 최적의 에너지를 향해 단조롭게 감소합니다.

# Fix the diagonalization subspace to the configurations found by SQD.
ci_strings = (result.sci_state.ci_strs_a, result.sci_state.ci_strs_b)
nelec = (num_elec_a, num_elec_b)

# Start from the MO basis in which we ran SQD.
hamiltonian_opt = ffsim.MolecularHamiltonian(
    hcore, eri, constant=nuclear_repulsion_energy
)

num_iters = 10
for i in range(num_iters):
    # Diagonalize over the fixed subspace in the current basis.
    myci = fci.selected_ci.SelectedCI()
    myci = fci.addons.fix_spin_(myci, ss=spin_sq)
    _, amplitudes = fci.selected_ci.kernel_fixed_space(
        myci,
        hamiltonian_opt.one_body_tensor,
        hamiltonian_opt.two_body_tensor,
        num_orbitals,
        nelec,
        ci_strs=ci_strings,
    )

    # Build the RDMs and record the energy before re-optimizing the orbitals.
    dm1, dm2 = myci.make_rdm12(amplitudes, num_orbitals, nelec)
    rdm = ffsim.ReducedDensityMatrix(dm1, dm2)
    energy = rdm.expectation(hamiltonian_opt).real
    print(f"Iteration {i}: energy = {energy:.8f}")

    # Rotate the Hamiltonian into the energy-minimizing basis for the next iteration.
    # optimize_orbitals returns the unitary matrix U minimizing
    # rdm.rotated(U).expectation(hamiltonian), equivalently
    # rdm.expectation(hamiltonian.rotated(U.conj().T)), so we rotate by U^dagger.
    orbital_rotation = ffsim.optimize_orbitals(rdm, hamiltonian_opt)
    hamiltonian_opt = hamiltonian_opt.rotated(orbital_rotation.T.conj())

Iteration 0: energy = -108.98452447
Iteration 1: energy = -108.99981993
Iteration 2: energy = -109.00585329
Iteration 3: energy = -109.00816569
Iteration 4: energy = -109.00936616
Iteration 5: energy = -109.01014322
Iteration 6: energy = -109.01069439
Iteration 7: energy = -109.01109308
Iteration 8: energy = -109.01138928
Iteration 9: energy = -109.01161411

결과 비교하기

오비탈 최적화는 고정 부분 공간 추정값을 개선하여, 정확한 에너지보다 높은 값을 유지하면서 그와의 차이를 상당 부분 줄입니다.

# Diagonalize once more in the final optimized basis to report the improved energy.
myci = fci.selected_ci.SelectedCI()
myci = fci.addons.fix_spin_(myci, ss=spin_sq)
_, amplitudes = fci.selected_ci.kernel_fixed_space(
    myci,
    hamiltonian_opt.one_body_tensor,
    hamiltonian_opt.two_body_tensor,
    num_orbitals,
    nelec,
    ci_strs=ci_strings,
)
dm1, dm2 = myci.make_rdm12(amplitudes, num_orbitals, nelec)
energy_after_oo = (
    ffsim.ReducedDensityMatrix(dm1, dm2).expectation(hamiltonian_opt).real
)

print(f"Exact energy:      {exact_energy:.8f}")
print(f"SQD energy (MO):   {sqd_energy:.8f}")
print(f"Energy after OO:   {energy_after_oo:.8f}")

Exact energy:      -109.04667178
SQD energy (MO):   -108.98469255
Energy after OO:   -109.01178727

SQD 실행하기​

오비탈 최적화하기​

여기서 오비탈 최적화가 도움이 되는 이유​

대각화와 오비탈 최적화를 번갈아 수행하기​

결과 비교하기​

SQD 실행하기

오비탈 최적화하기

여기서 오비탈 최적화가 도움이 되는 이유

대각화와 오비탈 최적화를 번갈아 수행하기

결과 비교하기