<jupyter_start><jupyter_text>Unitary ansatz entering the VQEThe final energy output of a VQE calculation will crucially depend on the ansatz/form of the parameterized unitary $\hat U(\boldsymbol{\theta})$ employed in state preparation. Here we review two popular approaches, the unitary coupled cluster and qubit coupled cluster methodologies, and benchmark them for energy calculations of small molecules.<jupyter_code>import numpy as np
import tequila as tq
from utility import *
threshold = 1e-6 #Cutoff for UCC MP2 amplitudes and QCC ranking gradients<jupyter_output><empty_output><jupyter_text>Unitary Coupled Cluster (UCC)The UCC ansatz is obtained by 'unitarizing' the traditional coupled cluster ansatz,$$ e^{\hat T} \rightarrow e^{\hat T - \hat T^\dagger} \equiv \hat U_{\text{UCC}}$$Due to non-commutativity of terms in $\hat T - \hat T^\dagger$, the UCC ansatz does not have a straightforward decomposition in terms of circuit primitives implementable on the quantum computer. Therefore, to obtain a form which can be compiled, we employ the Trotter approximation. The accuracy of the circuit ansatz relative to the exact UCC operator will be dependent on how many Trotter steps are employed. The number of Trotter steps is commonly set to its minimal value of one to avoid excessive circuit depth.<jupyter_code>trotter_steps = 1<jupyter_output><empty_output><jupyter_text>H2 in STO-3G basisBelow is a sample VQE simulation using the UCCSD ansatz compiled using a single trotter step for H$_2$ in minimal basis at $R=2.5$ (Angstrom). For comparison, we can run FCI to obtain the true ground state energy.<jupyter_code>xyz_data = get_molecular_data('h2', geometry=2.5, xyz_format=True)
basis='sto-3g'

h2 = tq.quantumchemistry.Molecule(geometry=xyz_data, basis_set=basis)

print('Number of spin-orbitals (qubits): {} \n'.format(2*h2.n_orbitals))

E_FCI = h2.compute_energy(method='fci')

print('FCI energy: {}'.format(E_FCI))<jupyter_output>Number of spin-orbitals (qubits): 4 

FCI energy: -0.9360549199436616<jupyter_text>The UCCSD VQE optimization is of the form$$E = \min_{\boldsymbol{\theta}} \langle \text{HF} | \hat U_{\text{UCC}}^\dagger(\boldsymbol{\theta}) \hat H  \hat U_{\text{UCC}} (\boldsymbol{\theta}) | \text{HF} \rangle $$To expedite the optimization process, we can set the initial guess of the amplitudes to zero, i.e. the optimization will begin at the Hartree Fock state. This heuristic is best suited for when Hartree Fock is believed to be in qualitative agreement with the true ground state. To further alleviate quantum resources, we can estimate the amplitudes using classical electronic structure methods (here, MP2 perturbation theory), and only include the unitaries with non-zero estimated amplitudes.<jupyter_code>H = h2.make_hamiltonian()

print("\nHamiltonian has {} terms\n".format(len(H)))

U_UCCSD = h2.make_uccsd_ansatz(initial_amplitudes='MP2',threshold=threshold, trotter_steps=trotter_steps)

E = tq.ExpectationValue(H=H, U=U_UCCSD)

print('\nNumber of UCCSD amplitudes: {} \n'.format(len(E.extract_variables())))

print('\nStarting optimization:\n')

result = tq.minimize(objective=E, method="BFGS", initial_values={k:0.0 for k in E.extract_variables()}, tol=1e-6)

print('\nObtained UCCSD energy: {}'.format(result.energy))<jupyter_output>Hamiltonian has 15 terms


Number of UCCSD amplitudes: 1 


Starting optimization:

Optimizer: <class 'tequila.optimizers.optimizer_scipy.OptimizerSciPy'> 
backend         : qulacs
samples         : None
save_history    : True
noise           : None

Method          : BFGS
Objective       : 1 expectationvalues
gradient        : 32 expectationvalues

active variables : 1

E=-0.70294360  angles= {(1, 0, 1, 0): 0.0}  samples= None
E=-0.92698008  angles= {(1, 0, 1, 0): -0.5644201636314392}  samples= None
E=-0.93354955  angles= {(1, 0, 1, 0): -0.756477668453395}  samples= None
E=-0.93605486  angles= {(1, 0, 1, 0): -0.6900732739717247}  samples= None
E=-0.93605492  angles= {(1, 0, 1, 0): -0.6904080632337539}  samples= None
E=-0.93605492  angles= {(1, 0, 1, 0): -0.6904070264134937}  samples= None
Optimization terminated successfully.
         Current function value: -0.936055
         Iterations: 5
         Function evaluations: 6
         Gradient evaluations: 6

Obtained UCCSD energy: -0.9[...]<jupyter_text>We see that the converged UCCSD energy is in exact agreement with the FCI energy, as expected for a $2$-electron system.  H2O in 6-31G basisNow let us try a larger problem of H$_2$O in 6-31G basis. However, we will restrict the active space. The unrestricted problem leads to a $14$-qubit Hamiltonian, and $34$ UCCSD amplitudes to optimize even after neglecting the zero MP2 amplitudes. Therefore, we will remove some orbital degrees of freedom which are less important in accurately describing the electronic structure. By freezing all orbitals other than $0b_1$, $1b_1$, $2a_1$, and $3a_1$, we reduce the problem to an $8$-qubit Hamiltonian with $8$ UCCSD variational amplitudes.<jupyter_code>xyz_data = get_molecular_data('h2o', geometry=1, xyz_format=True)

basis = '6-31g'
active = {'B1':[0,1], 'A1':[2,3]}
h2o = tq.quantumchemistry.Molecule(geometry=xyz_data, basis_set = basis, active_orbitals = active)

print('Number of spin-orbitals (qubits): {} \n'.format(2*h2o.n_orbitals))

E_FCI = h2o.compute_energy(method='fci')

print('FCI energy: {}'.format(E_FCI))<jupyter_output>There are known issues with some psi4 methods and frozen virtual orbitals. Proceed with fingers crossed for hf.
Number of spin-orbitals (qubits): 8 

There are known issues with some psi4 methods and frozen virtual orbitals. Proceed with fingers crossed for fci.
FCI energy: -75.98980140961186<jupyter_text>We will then run the UCCSD VQE simulation (***warning: tq.minimize will take several minutes - 1 hour + to converge for a VQE instance of this size.*** Smaller active spaces can be employed to lower VQE simulation runtimes).<jupyter_code>H = h2o.make_hamiltonian()

print("\nHamiltonian has {} terms\n".format(len(H)))

U_UCCSD = h2o.make_uccsd_ansatz(initial_amplitudes='MP2',threshold=threshold, trotter_steps=trotter_steps)

E = tq.ExpectationValue(H=H, U=U_UCCSD)

print('\nNumber of UCCSD amplitudes: {} \n'.format(len(E.extract_variables())))

print('\nStarting optimization:\n')

result = tq.minimize(objective=E, method="BFGS", initial_values={k:0.0 for k in E.extract_variables()}, tol=1e-4)

print('\nObtained UCCSD energy: {}'.format(result.energy))<jupyter_output>Hamiltonian has 185 terms


Number of UCCSD amplitudes: 8 


Starting optimization:

Optimizer: <class 'tequila.optimizers.optimizer_scipy.OptimizerSciPy'> 
backend         : qulacs
samples         : None
save_history    : True
noise           : None

Method          : BFGS
Objective       : 1 expectationvalues
gradient        : 512 expectationvalues

active variables : 8

E=-75.98071144  angles= {(3, 0, 3, 0): 0.0, (3, 1, 2, 0): 0.0, (2, 1, 2, 1): 0.0, (2, 0, 3, 1): 0.0, (3, 1, 3, 1): 0.0, (3, 0, 2, 1): 0.0, (2, 0, 2, 0): 0.0, (2, 1, 3, 0): 0.0}  samples= None
E=-75.57991345  angles= {(3, 0, 3, 0): -0.1516876220703125, (3, 1, 2, 0): -0.04468536376953125, (2, 1, 2, 1): -0.06198883056640625, (2, 0, 3, 1): -0.04468536376953125, (3, 1, 3, 1): -0.0497894287109375, (3, 0, 2, 1): 0.15348052978515625, (2, 0, 2, 0): -0.06015777587890625, (2, 1, 3, 0): 0.15348052978515625}  samples= None
E=-75.98420833  angles= {(3, 0, 3, 0): -0.012046029637663564, (3, 1, 2, 0): -0.0035486166174326273, (2, 1, [...]<jupyter_text>We obtain a UCCSD energy within chemical accuracy (< 1.6 miliHartree from FCI energy) for the H$_2$O molecule near equilibrium geometry. As the bond distance approaches the dissociation limit, the energy deviation from FCI is typically expected to increase as electronic correlations increase during covalent bond-breaking. Furthermore, as HF becomes energetically more distant from the FCI wavefunction, the initial guess of all amplitudes being zero may lead to a local minimum. One can then instead initialize the amplitudes using random guesses, and repeat for $n$ samples to attempt to find the global minimum. To facilitate this, one can run the commented lines below. This will execute $n$ instances of tq.minimize for uniformly random initial parameters in interval $[0, 2\pi]$ and return the lowest optimization result.<jupyter_code>#n=10
#result = minimize_E_random_guesses(objective=E, method='BFGS', tol=1e-4, n=n)<jupyter_output><empty_output><jupyter_text>Qubit Coupled Cluster (QCC)In contrast to UCC, the QCC methodology makes no direct reference to fermionic algebra and seeks to construct an efficient ansatz directly in qubit-space by finding multi-qubit Pauli strings (entanglers) which lower energy. This is done through an energy-lowering heuristic employing the energy gradient with respect to a Pauli strings variational amplitude. As opposed to UCCSD, the circuit depth and number of parameter is chosen to meet hardware limitations, i.e. one must choose how many exponentiated Pauli strings will be entering the QCC ansatz. H2 in STO-3G basisBelow we perform the entangler screening protocol for H2 in minimal basis, and obtain one grouping of entanglers with non-zero energy gradient. We then select one of them to be used in the QCC VQE simulation.<jupyter_code>xyz_data = get_molecular_data('h2', geometry=2.5, xyz_format=True)
basis='sto-3g'

h2 = tq.quantumchemistry.Molecule(geometry=xyz_data, basis_set='sto-3g')

hf_reference = hf_occ(2*h2.n_orbitals, h2.n_electrons)

H = h2.make_hamiltonian()

print("\nHamiltonian has {} terms\n".format(len(H)))

#Define number of entanglers to enter ansatz
n_ents = 1

#Rank entanglers using energy gradient criterion
ranked_entangler_groupings = generate_QCC_gradient_groupings(H.to_openfermion(), 
                                                             2*h2.n_orbitals, 
                                                             hf_reference, 
                                                             cutoff=threshold)

print('Grouping gradient magnitudes (Grouping : Gradient magnitude):')
for i in range(len(ranked_entangler_groupings)):
    print('{} : {}'.format(i+1,ranked_entangler_groupings[i][1]))


entanglers = get_QCC_entanglers(ranked_entangler_groupings, n_ents, 2*h2.n_orbitals)

print('\nSelected entanglers:')
for ent in entanglers:
    print(ent)<jupyter_output>Hamiltonian has 15 terms

Grouping gradient magnitudes (Grouping : Gradient magnitude):
1 : 0.2822

Selected entanglers:
1.0 [X0 Y1 X2 X3]<jupyter_text>Once the QCC ranking procedure has been ran, we can simulate the QCC VQE optimization with the generated entanglers. The VQE optimization for the QCC ansatz is of the form$$E = \min_{\boldsymbol{\Omega}, \boldsymbol{\tau}} \langle \boldsymbol{\Omega} | U_{\text{ENT}}^\dagger (\boldsymbol{\tau}) \hat H  U_{\text{ENT}} (\boldsymbol{\tau}) | \boldsymbol{\Omega} \rangle $$where $\boldsymbol{\Omega}$ denote collective Euler angles parameterizing single-qubit rotations, and $\boldsymbol{\tau}$ are entangler amplitudes.<jupyter_code>#Mean-field part of U (Omega):    
U_MF = construct_QMF_ansatz(n_qubits = 2*h2.n_orbitals)
#Entangling part of U:
U_ENT = construct_QCC_ansatz(entanglers)

U_QCC = U_MF + U_ENT

E = tq.ExpectationValue(H=H, U=U_QCC)

initial_vals = init_qcc_params(hf_reference, E.extract_variables())

#Minimize wrt the entangler amplitude and MF angles:
result = tq.minimize(objective=E, method="BFGS", initial_values=initial_vals, tol=1.e-6)

print('\nObtained QCC energy ({} entanglers): {}'.format(len(entanglers), result.energy))<jupyter_output>Optimizer: <class 'tequila.optimizers.optimizer_scipy.OptimizerSciPy'> 
backend         : qulacs
samples         : None
save_history    : True
noise           : None

Method          : BFGS
Objective       : 1 expectationvalues
gradient        : 18 expectationvalues

active variables : 9

E=-0.70294360  angles= {gamma_3: 0.0, beta_1: 3.141592653589793, gamma_0: 0.0, tau_0: 0.0, gamma_2: 0.0, gamma_1: 0.0, beta_2: 0.0, beta_0: 3.141592653589793, beta_3: 0.0}  samples= None
E=-0.77938639  angles= {gamma_3: 0.0, beta_1: 3.141592653589793, gamma_0: 0.0, tau_0: 0.2822100520133972, gamma_2: 0.0, gamma_1: 0.0, beta_2: 0.0, beta_0: 3.141592653589793, beta_3: 0.0}  samples= None
E=-0.93592356  angles= {gamma_3: 0.0, beta_1: 3.141592653589793, gamma_0: 0.0, tau_0: 1.411050260066986, gamma_2: 0.0, gamma_1: 0.0, beta_2: 0.0, beta_0: 3.141592653589793, beta_3: 0.0}  samples= None
E=-0.93603455  angles= {gamma_3: 0.0, beta_1: 3.141592653589793, gamma_0: 0.0, tau_0: 1.3689076316049542, gamma_2: 0.0, [...]<jupyter_text>We see that the QCC energy converged to the FCI energy with only a single entangler!  H2O in 6-31G basisLet us move on to the problem of H$_2$O in 6-31G basis using the same active space as the UCCSD example.<jupyter_code>xyz_data = get_molecular_data('h2o', geometry=1, xyz_format=True)

basis = '6-31g'
active = {'B1':[0,1], 'A1':[2,3]}
h2o = tq.quantumchemistry.Molecule(geometry=xyz_data, basis_set = basis, active_orbitals = active)
hf_reference = hf_occ(2*h2o.n_orbitals, h2o.n_electrons)


H = h2o.make_hamiltonian()

print("\nHamiltonian has {} terms\n".format(len(H)))

#Define number of entanglers to enter ansatz
n_ents = 6

#Rank entanglers using energy gradient criterion
ranked_entangler_groupings = generate_QCC_gradient_groupings(H.to_openfermion(), 
                                                             2*h2o.n_orbitals, 
                                                             hf_reference, 
                                                             cutoff=threshold)

print('Grouping gradient magnitudes (Grouping : Gradient magnitude):')
for i in range(len(ranked_entangler_groupings)):
    print('{} : {}'.format(i+1,ranked_entangler_groupings[i][1]))

entanglers = get_QCC_entanglers(ranked_entangler_groupings, n_ents, 2*h2o.n_orbitals)

print('\nSelected entanglers:')
for ent in entanglers:
    print(ent)<jupyter_output>There are known issues with some psi4 methods and frozen virtual orbitals. Proceed with fingers crossed for hf.

Hamiltonian has 185 terms

Grouping gradient magnitudes (Grouping : Gradient magnitude):
1 : 0.0758
2 : 0.0371
3 : 0.0371
4 : 0.031
5 : 0.0301
6 : 0.0249
7 : 0.0198
8 : 0.0198
9 : 0.0173
10 : 0.0173
11 : 0.0
12 : 0.0

Selected entanglers:
1.0 [X0 Y1 X6 X7]
1.0 [X0 Y3 X5 X6]
1.0 [X1 Y2 X4 X7]
1.0 [X2 Y3 X4 X5]
1.0 [X0 Y1 X4 X5]
1.0 [X2 Y3 X6 X7]<jupyter_text>With the 6 chosen entanglers, we can then run the QCC VQE simulation.<jupyter_code>#Mean-field part of U (Omega):    
U_MF = construct_QMF_ansatz(n_qubits = 2*h2o.n_orbitals)
#Entangling part of U:
U_ENT = construct_QCC_ansatz(entanglers)

U_QCC = U_MF + U_ENT

E = tq.ExpectationValue(H=H, U=U_QCC)

initial_vals = init_qcc_params(hf_reference, E.extract_variables())

#Minimize wrt the entangler amplitude and MF angles:
result = tq.minimize(objective=E, method="BFGS", initial_values=initial_vals, tol=1.e-4)


print('\nObtained QCC energy ({} entanglers): {}'.format(len(entanglers), result.energy))<jupyter_output>Optimizer: <class 'tequila.optimizers.optimizer_scipy.OptimizerSciPy'> 
backend         : qulacs
samples         : None
save_history    : True
noise           : None

Method          : BFGS
Objective       : 1 expectationvalues
gradient        : 44 expectationvalues

active variables : 22

E=-75.98071144  angles= {gamma_7: 0.0, tau_0: 0.0, gamma_5: 0.0, beta_1: 3.141592653589793, tau_3: 0.0, tau_5: 0.0, tau_4: 0.0, gamma_1: 0.0, beta_0: 3.141592653589793, gamma_3: 0.0, gamma_4: 0.0, beta_5: 0.0, beta_2: 3.141592653589793, beta_4: 0.0, gamma_2: 0.0, beta_7: 0.0, tau_1: 0.0, tau_2: 0.0, gamma_6: 0.0, beta_3: 3.141592653589793, beta_6: 0.0, gamma_0: 0.0}  samples= None
E=-75.98783364  angles= {gamma_7: 0.0, tau_0: 0.07584381103515625, gamma_5: 0.0, beta_1: 3.141592653589793, tau_3: 0.030994415283203125, tau_5: 0.02489471435546875, tau_4: 0.030078887939453125, gamma_1: 0.0, beta_0: 3.141592653589793, gamma_3: 0.0, gamma_4: 0.0, beta_5: 0.0, beta_2: 3.141592653589793, beta_4: 0.0, gamma_2: [...]<jupyter_text>We obtain chemical accuracy for water near equilibrium geometry with only 6 entanglers. The obtained energy is not as accurate as that of UCCSD for this problem, however the QCC optimization may be performed at a fraction of the UCCSD circuit depth. One can also increase the number of entanglers entering the QCC ansatz to increase accuracy. As a final check, one can always run $n$ VQE trials with random initial guesses to test if the optimization fell into a local minimum. ***(Warning: Completing n=10 trials may take a few minutes for this VQE instance).***<jupyter_code>n = 10
result = minimize_E_random_guesses(objective=E, method='BFGS', tol=1e-4, n=n)

print('\nObtained QCC energy ({} entanglers): {}'.format(len(entanglers), result))<jupyter_output>Optimizer: <class 'tequila.optimizers.optimizer_scipy.OptimizerSciPy'> 
backend         : qulacs
samples         : None
save_history    : True
noise           : None

Method          : BFGS
Objective       : 1 expectationvalues
gradient        : 44 expectationvalues

active variables : 22

E=-74.72863483  angles= {gamma_7: 1.2943697738877236, tau_0: 4.988480926540945, gamma_5: 0.9187555063348196, beta_1: 1.829015364772612, tau_3: 5.66076752282859, tau_5: 6.188245020859796, tau_4: 3.271842760348092, gamma_1: 5.386789573793364, beta_0: 2.16863319397154, gamma_3: 6.229607259068195, gamma_4: 3.8082917717557754, beta_5: 5.928799869627327, beta_2: 1.5977383419153086, beta_4: 1.4267087149043718, gamma_2: 1.5785345140689353, beta_7: 4.925899849934009, tau_1: 3.037487940287173, tau_2: 0.4091435985766853, gamma_6: 5.485195124612168, beta_3: 5.740555378623027, beta_6: 5.966738647629863, gamma_0: 3.9960084393454562}  samples= None
E=-74.85002750  angles= {gamma_7: 1.2934428024521767, tau_0: 5.0634[...]
{'fraction_non_alphanumeric': 0.107995962767747, 'fraction_numerical': 0.10121116967589996, 'mean_word_length': 4.383338364020525, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '220702', 'license': 'MIT', 'n_tokens_mistral': 7386, 'n_tokens_neox': 6166, 'n_words': 1893, 'path': 'Project_2_VQE_Molecules/S3_Unitary_Ansatz.ipynb', 'repo_name': 'CDLQuantum2020Week2Team1/CohortProject_2020', 'size': '17834'}