Unitary Selective Coupled Cluster
The Unitary Coupled-Cluster with Singles and Doubles (UCCSD) method parameterizes coupled-cluster variational parameters for each fermionic excitation from a reference state, whether single or double excitations. This ansatz has been successfully applied to compute accurate electronic energies for small molecules. However, for large molecules, this approach can introduce many redundant, unimportant excitations, resulting in a large number of parameters to optimize and excessively long circuits
The UsCC method systematically reduces the energy error with an increasing ansatz size for a set of test molecules. A significant advantage of the UsCC method is that expanding the ansatz does not necessitate additional measurements on a quantum computer. The implementation of the UsCC method begins with the reference Hartree–Fock state. One- and two-body excitations (from occupied to unoccupied orbitals) with corresponding Hamiltonian matrix elements above a given threshold \(\epsilon\) are included in the ansatz. A Variational Quantum Eigensolver (VQE) is then performed on this ansatz until convergence, completing the first iteration. In subsequent iterations, the threshold is updated such that \(\epsilon_i = \frac{\epsilon_{i-1}}{2}\). New, higher-order excitations are generated by applying one- and two-body excitations to those already included in the ansatz. The metric to compare to the threshold for inclusion of these excitations is obtained by multiplying the parameters obtained after optimization at the previous step with the Hamiltonian matrix element corresponding to the single or double excitation added. As in previous iterations, coefficients that are above the threshold are included in the ansatz and then optimized through VQE. The optimization continues until a convergence criterion is met or sufficient excitations have been included.
The ansatz generation starts with the reference Hartree–Fock state, with no initial excitation amplitudes available. First, all possible single- and double-electronic excitations are generated. These excitations are expressed in spin-block notation: the occupied orbitals, from which excitations are performed, are labeled $i$, $j$, $k$, and $l$, while the virtual orbitals, to which electrons are excited, are labeled $p$, $q$, $r$, and $s$. For instance, single and double excitations are denoted as $[i, p]$ and $[i, j, p, q]$, respectively.
UsCCSDTQ-VQE Algorithm
- Generate single and double excitations for a given molecule.
For all single (S) and double (D) excitations \([i, p]\) and \([i, j, p, q]\) in UCCSD:
- If \(h_1[i, p]\) and \(h_2[i, j, p, q]\) are larger than \(\epsilon_1\) then:
- Add to ansatz.
- end if
end for
Repeat.
Run VQE with the current ansatz to compute energy, update amplitudes for each excitation present in ansatz.
For each single $[i, p]$ or double $[i, j, p, q]$ excitation present in ansatz using $t_1$ and $t_2$ values from the previous iteration and additional excitations $[k, r]$ or $[k, l, r, s]$, generate triple and quadruple excitations with the following coefficients:
- $t_1[i, p] \cdot h_2[j, k, q, r]$
- $h_1[i, p] \cdot t_2[j, k, q, r]$
- $t_2[i, j, p, q] \cdot h_1[k, r]$
- $h_2[i, j, p, q] \cdot t_1[k, r]$
- $t_2[i, j, p, q] \cdot h_2[k, l, r, s]$
For each excitation, if the absolute value of the largest coefficient computed in step 3 is larger than $\epsilon_n$ on iteration $n$, add this excitation to ansatz.
Code Source Implementation
Step 1
Import necessary library and define the fucntions
1from openvqe.common_files.qubit_pool import QubitPool
2from openvqe.common_files.molecule_factory import MoleculeFactory
3from openvqe.ucc_family.get_energy_ucc import EnergyUCC
4from qat.fermion.transforms import (get_jw_code, recode_integer)
5from qat.fermion.chemistry.pyscf_tools import perform_pyscf_computation
6from qat.fermion.chemistry.pyscf_tools import perform_pyscf_computation
7from qat.fermion.chemistry import MolecularHamiltonian, MoleculeInfo
8from qat.fermion.chemistry.ucc import guess_init_params, get_hf_ket, get_cluster_ops
9from qat.fermion.transforms import transform_to_jw_basis # , transform_to_bk_basis, transform_to_parity_basis
10from qat.fermion.transforms import recode_integer, get_jw_code # , get_bk_code, get_parity_code
11from qat.lang.AQASM import Program, X
12from qat.fermion.trotterisation import make_trotterisation_routine
13from qat.fermion import ElectronicStructureHamiltonian
14import scipy.optimize
15from qat.fermion.chemistry.ucc_deprecated import build_ucc_ansatz
16from qat.lang.AQASM import Program
17from qat.qpus import get_default_qpu
18from openvqe.common_files.circuit import count
19import numpy as np
20
21class EnergyUCC:
22 def ucc_action(self, theta_current, hamiltonian_sp, cluster_ops_sp, hf_init_sp, theta_thresh=1e-7):
23 """
24 It maps the exponential of cluster operators ("cluster_ops_sp") associated by their parameters ("theta_current")
25 using the CNOTS-staircase method, which is done by "build_ucc_ansatz" which creates the circuit on the top of
26 the HF-state ("hf_init_sp"). Then, this function also calculates the expected value of the hamiltonian ("hamiltonian_sp").
27
28 Parameters
29 ----------
30 theta_current: List<float>
31 the Parameters of the cluster operators
32
33 hamiltonian_sp: Hamiltonian
34 Hamiltonian in the spin representation
35
36 cluster_ops_sp: list[Hamiltonian]
37 list of spin cluster operators
38
39 hf_init_sp: int
40 the integer corresponds to the hf_init (The Hartree-Fock state in integer representation) obtained by using
41 "qat.fermion.transforms.record_integer".
42
43 Returns
44 --------
45 res.value: float
46 the resulted energy
47
48 """
49 qpu = 0
50 prog = 0
51 reg = 0
52 qpu = get_default_qpu()
53 prog = Program()
54 reg = prog.qalloc(hamiltonian_sp.nbqbits)
55 qrout = 0
56 for n_term, (term, theta_term) in enumerate(zip(cluster_ops_sp, theta_current)):
57 init = hf_init_sp if n_term == 0 else 0
58 qprog = build_ucc_ansatz([term], init, n_steps=1)
59 if abs(theta_term) > theta_thresh:
60 prog.apply(qprog([theta_term]), reg)
61 circ = prog.to_circ()
62
63 return circ
64
65 def get_optimization_func(self, circ, qpu, hamiltonian_sp, method, nqbits, psi0, energy_list, fid_list):
66 # below, I show an example of minimization of the energy where I store the energy and fidelity for each parameter
67 # I use a gradient-free procedure.
68 def my_func(x):
69 """returns energy given parameter x, and stores it + the fidelity of state"""
70 circ1 = circ.bind_variables(
71 {k: v for k, v in zip(sorted(circ.get_variables()), x)}
72 )
73 res0 = qpu.submit(circ1.to_job(observable=hamiltonian_sp))
74 energy = res0.value
75 energy_list[method].append(energy)
76
77 # additional computation to compute fidelity (just for my own information)
78 res = qpu.submit(circ1.to_job())
79 psi = np.zeros((2**nqbits,), complex)
80 for sample in res:
81 psi[sample.state.int] = sample.amplitude
82 fid = abs(psi.conj().dot(psi0)) ** 2
83 fid_list[method].append(fid)
84
85 return energy
86
87 return my_func
88
89
90 def get_grad_func(self, circ, qpu, hamiltonian_sp):
91 # here I show a gradient-based minimization strategy
92 def my_grad(x):
93 grads = circ.to_job(observable=hamiltonian_sp).gradient()
94 grad_list = []
95 for var_name in sorted(circ.get_variables()):
96 list_jobs = grads[var_name]
97 # list_jobs contains jobs to compute E(theta+pi/2) and E(theta-pi/2)
98 # the gradient w.r.t theta is then 0.5 (E(theta+pi/2) - E(theta-pi/2))
99 grad = 0.0
100 for ind in range(len(list_jobs)):
101 circ1 = list_jobs[ind].circuit.bind_variables(
102 {k: v for k, v in zip(circ.get_variables(), x)}
103 )
104 job = circ1.to_job(observable=list_jobs[ind].observable)
105 res = qpu.submit(job)
106 grad += 0.5 * res.value
107 grad_list.append(grad)
108 return grad_list
109
110
111 def prepare_state_ansatz(
112 self, hamiltonian_sp, cluster_ops_sp, hf_init_sp, parameters
113 ):
114 """
115 It constructs the trial wave function (ansatz)
116
117 Parameters
118 ----------
119 hamiltonian_sp: Hamiltonian
120 Hamiltonian in the spin representation
121
122 cluster_ops_sp: list[Hamiltonian]
123 list of spin cluster operators
124
125 hf_init_sp: int
126 the integer corresponds to the hf_init (The Hartree-Fock state in integer representation) obtained by using
127 "qat.fermion.transforms.record_integer".
128
129 parameters: List<float>
130 the Parameters for the trial wave function to be constructed
131
132
133
134 Returns
135 --------
136 curr_state: qat.core.Circuit
137 the circuit that represent the trial wave function
138
139 """
140 qpu = get_default_qpu()
141 prog = Program()
142 reg = prog.qalloc(hamiltonian_sp.nbqbits)
143 for n_term, (term, theta_term) in enumerate(zip(cluster_ops_sp, parameters)):
144 init = hf_init_sp if n_term == 0 else 0
145 qprog = build_ucc_ansatz([term], init, n_steps=1)
146 prog.apply(qprog([theta_term]), reg)
147 circ = prog.to_circ()
148 curr_state = circ
149 return curr_state
150
151 def get_energies(
152 self,
153 hamiltonian_sp,
154 cluster_ops_sp,
155 hf_init_sp,
156 theta_current,
157 fci,
158 method="BFGS"
159 ):
160 """
161 It calls internally the functions "ucc_action" and "prepare_state_ansatz", and uses scipy.optimize to
162 return the properties of the ucc energy and wave function.
163
164 Parameters
165 ----------
166 hamiltonian_sp: Hamiltonian
167 Hamiltonian in the spin representation
168
169 cluster_ops_sp: list[Hamiltonian]
170 list of spin cluster operators
171
172 hf_init_sp: int
173 the integer corresponds to the hf_init (The Hartree-Fock state in integer representation) obtained by using
174 "qat.fermion.transforms.record_integer".
175
176 theta_current: List<float>
177 the Parameters of the cluster operators of "cluster_ops_sp"
178
179
180 theta_current2: List<float>
181 the Parameters of the cluster operators of "pool_generator"
182
183
184 theta_current2: List<float>
185 the Parameters of the cluster operators of "pool_generator"
186
187 fci: float
188 the full configuration interaction energy (for any basis set)
189
190
191 Returns
192 --------
193 iterations: Dict
194 the minimum energy and the optimized parameters
195
196 result: Dict
197 the number of CNOT gates, the number of operators/parameters, and the substraction of the optimized energy from fci.
198
199 """
200
201 iterations = {
202 "minimum_energy_result_guess": [],
203 "theta_optimized_result": [],
204 }
205 result = {}
206 tolerance = 10 ** (-4)
207 print("tolerance= ", tolerance)
208 print("method= ", method)
209
210 theta_optimized_result = []
211
212
213
214 opt_result = scipy.optimize.minimize(
215 lambda theta: self.ucc_action(
216 theta, hamiltonian_sp, cluster_ops_sp, hf_init_sp
217 ),
218 x0=theta_current,
219 method=method,
220 tol=tolerance,
221 options={"maxiter": 50000, "disp": True},
222 )
223
224 xlist = opt_result.x
225
226 for si in range(len(theta_current)):
227 theta_optimized_result.append(xlist[si])
228
229 curr_state_result = self.prepare_state_ansatz(
230 hamiltonian_sp, cluster_ops_sp, hf_init_sp, theta_optimized_result
231 )
232
233 gates = curr_state_result.ops
234 cnot = count("CNOT", gates)
235 iterations["minimum_energy_result_guess"].append(opt_result.fun)
236 iterations["theta_optimized_result"].append(theta_optimized_result)
237 result["CNOT"] = cnot
238 result["len_op"] = len(theta_optimized_result)
239 result["energies_substracted_from_FCI"] = abs(opt_result.fun - fci)
240 return iterations, result
Step 2
Callout the function and excecute the program
1from openvqe.common_files.molecule_factory import MoleculeFactory
2#from openvqe.ucc_family.get_energy_ucc import EnergyUCC
3from openvqe.common_files.generator_excitations import _apply_transforms
4from qat.fermion.transforms import get_jw_code, recode_integer
5from qat.core import Term
6from qat.fermion import FermionHamiltonian
7
8
9def get_h(hamiltonian, exc):
10 hamiltonian_string = str(hamiltonian)
11 lines = hamiltonian_string.strip().split("\n")
12 if len(exc) == 2:
13 id_str = f"(Cc|[{exc[0]}, {exc[1]}])"
14 else:
15 id_str = f"(CCcc|[{exc[0]}, {exc[1]}, {exc[2]}, {exc[3]}])"
16
17 lines = [line for line in lines if id_str in line]
18 if len(lines) == 0:
19 print("zero h:", exc)
20 return complex(0.0)
21 assert len(lines) == 1
22
23 line = lines[0]
24 start_ind = line.find("(") + 1
25 end_ind = line.find(")")
26 value = line[start_ind:end_ind]
27 return complex(value)
28
29
30def init(molecule_symbol):
31 (
32 hamiltonian,
33 hamiltonian_sp,
34 _,
35 active_noons,
36 _,
37 info,
38 ) = MoleculeFactory().generate_hamiltonian(
39 molecule_symbol, active=False, transform="JW"
40 )
41
42 (
43 _,
44 cluster_ops,
45 _,
46 theta_MP2,
47 hf_init,
48 ) = MoleculeFactory().generate_cluster_ops(
49 molecule_symbol, type_of_generator="UCCSD", transform="JW", active=False
50 )
51
52 hf_init_sp = recode_integer(hf_init, get_jw_code(hamiltonian_sp.nbqbits))
53 fci = info["FCI"]
54 print("FCI", fci)
55 nqbits = len(active_noons)
56
57 excitations = [tuple(op.terms[0].qbits) for op in cluster_ops]
58 current_coeffs = {}
59
60 # TODO:
61 # theta_MP2 = [x if x != 0 else 0.01 for x in theta_MP2]
62 theta_MP2 = [0.01 for x in theta_MP2]
63 for i, exc in enumerate(excitations):
64 current_coeffs[exc] = {"h": get_h(hamiltonian, exc)}
65
66 return hamiltonian_sp, hf_init_sp, fci, current_coeffs, nqbits
67
68
69def build_cluster_ops(exc_list, nqbits: int):
70 """Adapted from myqlm"""
71 t_opti = []
72
73 for op_index in exc_list:
74 current_excitation_op = []
75
76 op_description, indices, indices_conj = None, None, None
77 if len(op_index) == 2:
78 op_description = "Cc"
79 indices, indices_conj = list(op_index), [op_index[1], op_index[0]]
80
81 elif len(op_index) == 4:
82 op_description = "CCcc"
83 indices, indices_conj = list(op_index), [
84 op_index[2],
85 op_index[3],
86 op_index[0],
87 op_index[1],
88 ]
89 elif len(op_index) == 6:
90 op_description = "CCCccc"
91 indices, indices_conj = list(op_index), [
92 op_index[3],
93 op_index[4],
94 op_index[5],
95 op_index[0],
96 op_index[1],
97 op_index[2],
98 ]
99
100 elif len(op_index) == 8:
101 op_description = "CCCCcccc"
102 indices, indices_conj = list(op_index), [
103 op_index[4],
104 op_index[5],
105 op_index[6],
106 op_index[7],
107 op_index[0],
108 op_index[1],
109 op_index[2],
110 op_index[3],
111 ]
112
113 current_excitation_op.append(Term(1j, op_description, indices))
114 current_excitation_op.append(Term(-1j, op_description, indices_conj))
115 t_opti.append(FermionHamiltonian(nqbits=nqbits, terms=current_excitation_op))
116
117 return t_opti
118
119
120def update_ansatz(current_excitations, current_coeffs, thresh, nqbits):
121 all_excitations = list(current_coeffs.keys())
122
123 for optimized_exc in current_excitations:
124 for exc in all_excitations:
125 exc_len = len(exc)
126 opt_exc_len = len(optimized_exc)
127
128 # single-single
129 if exc_len == 2 and opt_exc_len == 2:
130 new_exc = sorted([exc[0], optimized_exc[0]]) + sorted(
131 [exc[1], optimized_exc[1]]
132 )
133 elif exc_len == 2 and opt_exc_len == 4:
134 new_exc = sorted([exc[0], optimized_exc[0], optimized_exc[1]]) + sorted(
135 [exc[1], optimized_exc[2], optimized_exc[3]]
136 )
137 elif exc_len == 4 and opt_exc_len == 2:
138 new_exc = sorted([exc[0], exc[1], optimized_exc[0]]) + sorted(
139 [exc[2], exc[3], optimized_exc[1]]
140 )
141 elif exc_len == 4 and opt_exc_len == 4:
142 new_exc = sorted(
143 [exc[0], exc[1], optimized_exc[0], optimized_exc[1]]
144 ) + sorted([exc[2], exc[3], optimized_exc[2], optimized_exc[3]])
145 else:
146 continue
147
148 if len(set(new_exc)) != len(new_exc):
149 continue
150 th = 0
151 ht = 0
152 if "t" in current_coeffs[optimized_exc] and "h" in current_coeffs[exc]:
153 # print(
154 # thresh,
155 # f"t{opt_exc_len//2}",
156 # f"h{exc_len//2}",
157 # len(new_exc) // 2,
158 # current_coeffs[optimized_exc]["t"],
159 # current_coeffs[exc]["h"],
160 # )
161 th = current_coeffs[optimized_exc]["t"] * current_coeffs[exc]["h"]
162 if "t" in current_coeffs[exc] and "h" in current_coeffs[optimized_exc]:
163 # print(
164 # thresh,
165 # f"h{opt_exc_len//2}",
166 # f"t{exc_len//2}",
167 # len(new_exc) // 2,
168 # current_coeffs[optimized_exc]["h"],
169 # current_coeffs[exc]["t"],
170 # )
171 ht = current_coeffs[optimized_exc]["h"] * current_coeffs[exc]["t"]
172
173 new_exc = tuple(new_exc)
174 if new_exc not in current_coeffs:
175 current_coeffs[new_exc] = {}
176
177 current_coeffs[new_exc]["th"] = max(
178 current_coeffs[new_exc].get("th", 0), abs(th)
179 )
180 current_coeffs[new_exc]["ht"] = max(
181 current_coeffs[new_exc].get("ht", 0), abs(ht)
182 )
183
184 for exc in current_coeffs.keys():
185 if exc in current_excitations:
186 continue
187
188 coeffs = list(map(abs, current_coeffs[exc].values()))
189 max_coeff = max(coeffs)
190 if max_coeff > thresh:
191 print("ADDED", exc)
192 current_excitations.append(exc)
193
194 current_theta = [
195 current_coeffs[exc]["t"] if "t" in current_coeffs[exc] else 0.01
196 for exc in current_excitations
197 ]
198
199 cluster_ops = build_cluster_ops(current_excitations, nqbits)
200 _, _, cluster_ops_sp = _apply_transforms(cluster_ops, "JW")
201
202 return current_excitations, current_theta, cluster_ops_sp
203
204
205def iterate(hamiltonian_sp, ansatz_ops, hf_init_sp, theta_current, fci):
206 iterations, result = EnergyUCC().get_energies(
207 hamiltonian_sp, ansatz_ops, hf_init_sp, theta_current, fci, "BFGS"
208 )
209 new_theta = iterations["theta_optimized_result"][0]
210 energy = result["energies_substracted_from_FCI"]
211 print("ENERGY ERROR", energy)
212
213 return new_theta, energy
214
215
216def energy_stop_condition(energy):
217 return False
218
219
220thresholds = iter(
221 [0.04, 0.02, 0.01, 0.005, 0.002, 0.001, 0.0005, 0.0002, 1e-4, 1e-5, 1e-6, 1e-7]
222)
223hamiltonian_sp, hf_init_sp, fci, current_coeffs, nqbits = init("LIH")
224current_excitations = []
225thresh = next(thresholds)
226
227just_count = 0
228while True:
229 prev_num_excitations = len(current_excitations)
230 current_excitations, current_theta, cluster_ops_sp = update_ansatz(
231 current_excitations, current_coeffs, thresh, nqbits
232 )
233 if len(current_excitations) == prev_num_excitations:
234 try:
235 thresh = next(thresholds)
236 except StopIteration:
237 break
238 continue
239 print("Optimization:", just_count)
240 just_count += 1
241
242
243
244 new_theta, energy = iterate(
245 hamiltonian_sp, cluster_ops_sp, hf_init_sp, current_theta, fci
246 )
247 if energy_stop_condition(energy):
248 break
249
250 for i, exc in enumerate(current_excitations):
251 current_coeffs[exc]["t"] = new_theta[i]
252
253
254for exc in current_excitations:
255 print(exc)
Analytical Results
We test UsCCSDTQ-VQE on H$_4$ and LiH molecules using PMRS method (described in the previous section). We emulate 8 qubits in STO-3G basis set. Numerical results are shown in Figure 1 for H$_4$ and 2 LiH for as below:
Figure 1: PMRS application on UsCCSDTQ-VQE algorithm using H4 molecule (8 qubits) at 0.85Å, using BFGS optimizer: (left figure) shows the energy obtained from UsCC simulations with respect to the number of selective parameters. (Right figure) shows the number of function evaluations per each iteration. Error is the difference between UsCCDTQ energies and FCI. STO-3G basis set is used.
Figure 2: PMRS application on UsCCSDTQ-VQE algorithm using LiH molecule (12 qubits) at 1.45Å, using BFGS optimizer: (left figure) shows the energy obtained from UsCC simulations with respect to the number of selective parameters. (Right figure) shows the number of function evaluations per each iteration. Error is the difference between UsCCDTQ energies and FCI. STO-3G basis set is used.
We observe that only three iterations are needed to reach $10^{-9}$ Hartee for H$_4$ (with overall 19 parameters to optimize, with function evaluations 225 (first iteration), 198 (second iteration) and 180( third iteration)) and 13 iterations are needed for LiH reach to 1.06 $\times 10^{-5}$ Hartree, with about 45 function evaluations. We observe that PMRS helps reduce the number of function evaluations in both systems, which would exceed 1000 if PMRS were not used. This is demonstrated in the below table
| $\bm{r_{\text{Li} - \text{H}}}$ | Selective Parameters | Functions Evaluation | Gradient Evaluation | Error (Ha) |
|---|---|---|---|---|
| $0.5\mathring{\text{A}}$ | 54 | 10938 | 226 | $3.7\times 10^{-13}$ |
| $1.0\mathring{\text{A}}$ | 68 | 10795 | 247 | $4.4\times 10^{-13}$ |
| $1.5\mathring{\text{A}}$ | 72 | 10640 | 242 | $7.84\times 10^{-13}$ |
| $2.0\mathring{\text{A}}$ | 68 | 11406 | 263 | $1.96\times 10^{-12}$ |
Table 1: Test UsCCSDTQ-VQE on LiH 12 qubits. BFGS optimizer, tolerance $10^{-6}$ (Hartree), STO-3G basis set is used. (Molecule LiH at $r_{\text{Li} -\text{H}} = \left[0.5 - 2.0 \right] \mathring{\text{A}}$)
References
Fedorov, Dmitry A., et al. "Unitary selective coupled-cluster method." Quantum 6 (2022): 703.Haidar, Mohammad, et al. "Open source variational quantum eigensolver extension of the quantum learning machine for quantum chemistry." Wiley Interdisciplinary Reviews: Computational Molecular Science 13.5 (2023): e1664.
