13. Unitary Coupled Cluster (UCC) for pairing Hamiltonians#
Qiskit Natureの中に、UCC系のansatz及びVQEの実装がある:
これらを利用してpairing HamiltonianでUCC計算を行ってみよう。 QiskitのUCC実装は、軌道をup/down spinに分けるなど、電子系/量子化学での計算を念頭においている。 したがって、核子系や、pairingのpair自由度を畳んでhard-core boson化した場合は、適切なオプションの指定が必要となる。
その他、qiskit_algorithms.minimum_eigensolvers.VQEが、古い(v2ではduplicated)Estimatorを使うので、将来的にはコードの修正が必要となる点に注意。
Nq = 6, Nocc = 3の場合の厳密解を求めておこう。
Show code cell source
import numpy as np
from itertools import combinations
from qiskit.quantum_info import SparsePauliOp
from qiskit_nature.second_q.mappers import JordanWignerMapper
from qiskit_nature.second_q.circuit.library.ansatzes import UCC
from qiskit_nature.second_q.circuit.library import UCCSD, HartreeFock
from qiskit_algorithms.optimizers import SLSQP
from qiskit_algorithms.minimum_eigensolvers import AdaptVQE, VQE
from qiskit.primitives import Estimator as oEstimator #duplicate!
from qiskit_aer.primitives import Estimator, EstimatorV2
class PairingHamiltonian:
def __init__(self, Norb, Nocc, gval, delta_eps=1.0):
self.Norb = Norb
self.Nocc = Nocc
self.delta_eps = delta_eps
self.gval = gval
self.basis = self.make_basis()
self.epsilon = self.eval_epsilon()
self.Hmat = self.eval_Hmat()
def make_basis(self):
self.basis = []
for occ in combinations(range(self.Norb), self.Nocc):
self.basis.append(occ)
return self.basis
def eval_epsilon(self):
self.epsilon = [ 2 * i * self.delta_eps for i in range(self.Norb) ]
return self.epsilon
def eval_Hmat(self):
dim = len(self.basis)
self.Hmat = np.zeros((dim, dim))
for bra_idx, bra in enumerate(self.basis):
for ket_idx, ket in enumerate(self.basis):
# Hamming distance
diff = [ i for i in bra if i not in ket ]
same = [ i for i in bra if i in ket ]
# for SPE term
if bra_idx == ket_idx:
self.Hmat[bra_idx, ket_idx] += np.sum( [self.epsilon[i] for i in same])
self.Hmat[bra_idx, ket_idx] += - self.gval * len(same)
# for pairing term
if len(diff) == 1:
self.Hmat[bra_idx, ket_idx] = - self.gval
return self.Hmat
def tuple_to_bitstring(tup, Norb, rev=True):
bitint = 0
for i in tup:
bitint += 2**i
if rev:
bitstring = "|"+format(bitint, f'0{Norb}b')[::-1]+">"
else:
bitstring = "|"+format(bitint, f'0{Norb}b')+">"
return bitstring
params_exact = np.array([-0.48104276, -1.03976498, -0.98963981, -1.18481738, -0.54832984])
Norb = 6
Nocc = 3
gval = 0.33
Hamil = PairingHamiltonian(Norb, Nocc, gval)
evals, evecs = np.linalg.eigh(Hamil.Hmat)
evals = np.linalg.eigvalsh(Hamil.Hmat)
Egs_exact = evals[0]
E_HF = Hamil.Hmat[0,0]
print("Egs_exact: ", Egs_exact, " E_HF", E_HF)
SPEs = Hamil.epsilon
pauli_list = [ ]
obs = [ ]
coeffs = [ ]
# I term
coeff = 0.0
op = "I" * Norb
for i in range(Norb):
coeff += 0.5 * ( SPEs[i] - Hamil.gval )
obs += [op]
coeffs += [coeff]
# -Zp term
for i in range(Norb):
op = "I" * Norb
op = op[:i] + "Z" + op[i+1:]
coeff = -0.5 * ( SPEs[i] - Hamil.gval )
op = op[::-1]
obs += [op]
coeffs += [coeff]
# XX+YY term
for i in range(Hamil.Norb):
for j in range(i+1, Hamil.Norb):
factor = - Hamil.gval / 2
op = "I" * Norb
op = op[:i] + "X" + op[i+1:j] + "X" + op[j+1:]
op = op[::-1]
obs += [op]
coeffs += [ factor ]
op = "I" * Norb
op = op[::-1]
op = op[:i] + "Y" + op[i+1:j] + "Y" + op[j+1:]
obs += [op]
coeffs += [ factor ]
hamiltonian_op = SparsePauliOp(obs, coeffs)
Egs_exact: 4.741264131588444 E_HF 5.01
13.1. UCCによるVQE計算#
UCCやUCCSDクラスを使う場合は、
Norb: 軌道数(up+down spinそれぞれを同じ数持つことが想定されている)num_particles: up/down spinのparticle数generalized:Trueにすると、hole-holeのような項も考慮される。preserve_spin:Trueにすると、up/down spinの数を保持するようなUCCが実装される。(proton/neutronに応用できる)
などのオプションがある。詳しくはdocumentationを参照のこと。
いま、pairの自由度を畳んでいるので、UCCS的な使い方で、pUCCDに対応する。
num_particles = (Nocc//2, Nocc-Nocc//2) #←粒子数が多いとき、α/β側の軌道を全部埋めると励起が考えられないようなので、少し変な実装にした。
mapper = JordanWignerMapper()
init_state=HartreeFock(
Norb//2,
num_particles=num_particles,
qubit_mapper=mapper,
)
display(init_state.draw())
ucc_s = UCC(
Norb//2,
num_particles=num_particles,
excitations="s",
qubit_mapper = mapper,
generalized=False,
initial_state=init_state,
preserve_spin=False,
reps=1
)
print("UCCS: ")
print("Excitation list: ", ucc_s.excitation_list)
print("Number of parameters: ", ucc_s.num_parameters)
┌───┐
q_0: ┤ X ├
└───┘
q_1: ─────
q_2: ─────
┌───┐
q_3: ┤ X ├
├───┤
q_4: ┤ X ├
└───┘
q_5: ─────
UCCS:
Excitation list: [((0,), (1,)), ((0,), (2,)), ((0,), (5,)), ((3,), (1,)), ((3,), (2,)), ((3,), (5,)), ((4,), (1,)), ((4,), (2,)), ((4,), (5,))]
Number of parameters: 9
回路ができたので、VQEを実行してみよう。 適当にiteration回数を設定して...
niter = 50
optimizer = SLSQP(maxiter=niter, ftol=1e-6)
vqe = VQE(oEstimator(), ucc_s, optimizer)
vqe.initial_point = np.zeros(ucc_s.num_parameters, dtype=float)
vqe_result = vqe.compute_minimum_eigenvalue(hamiltonian_op)
binding_energy = vqe_result.optimal_value
print("E(pUCCD): ", binding_energy, "Egs_exact: ", Egs_exact, " E_HF", E_HF)
/var/folders/w4/qfsp76x90732_cjtryhq_19w0000gn/T/ipykernel_19095/3267568472.py:4: DeprecationWarning: The class ``qiskit.primitives.estimator.Estimator`` is deprecated as of qiskit 1.2. It will be removed no earlier than 3 months after the release date. All implementations of the `BaseEstimatorV1` interface have been deprecated in favor of their V2 counterparts. The V2 alternative for the `Estimator` class is `StatevectorEstimator`.
vqe = VQE(oEstimator(), ucc_s, optimizer)
E(pUCCD): 4.750238826425197 Egs_exact: 4.741264131588444 E_HF 5.01
DOCIの場合の配位が20個、pUCCDでは9個の自由度でcorrelation energy(ここではE(HF)-E(FCI))の95%以上を説明していることがわかる。