VQE Application in Quantum Chemistry Computing

   

Overview

Quantum chemistry refers to solving the numerical values of the time-dependent or time-independent Schrödinger equations by using the basic theory and method of quantum mechanics. Quantum chemical simulation on high-performance computers has become an important method to study the physical and chemical properties of materials. However, the exact solution of the Schrödinger equation has exponential complexity, which severely constrains the scale of the chemical system that can be simulated. The development of quantum computing in recent years provides a feasible way to solve this problem. It is expected that the Schrödinger equation can be solved with high accuracy on quantum computers under the complexity of polynomials.

Peruzzo et al. first applied the VQE and unitary coupled-cluster theory to quantum chemistry simulation in 2014 to solve the ground state energy of He-H+. The VQE is a hybrid quantum-classical algorithm and is widely used in chemical simulation based on quantum algorithms. This tutorial describes how to use the VQE to solve the ground-state energy of a molecular system.

This tutorial consists of the following parts:

1. Introduction to the quantum chemistry

2. VQE application

3. Using MindSpore Quantum to perform VQE simulation with efficient and automatic derivation

This document applies to the CPU environment.

Environment Preparation

In this tutorial, the following environments need to be installed:

• NumPy

• SciPy

• MindSpore Quantum

• MindSpore

• PySCF

• OpenFermion

• OpenFermion-PySCF

The preceding dependencies can be installed by running the pip command.

The chemical simulation packages PySCF and OpenFermion are only supported on Linux or Mac systems. On Windows 10, they need to be installed within the Ubuntu subsystem.

Importing Dependencies

Import the modules on which this tutorial depends.

from openfermion.chem import MolecularData
from openfermionpyscf import run_pyscf
from mindquantum.core.gates import X
from mindquantum.core.circuit import Circuit
from mindquantum.core.operators import Hamiltonian
from mindquantum.simulator import Simulator
from mindquantum.algorithm.nisq import generate_uccsd
import mindspore as ms

ms.set_context(mode=ms.PYNATIVE_MODE, device_target="CPU")

Quantum Chemistry Computing Method

The core of quantum chemistry is to solve the Schrödinger equation. In general, the solution of time-dependent Schrödinger equation is more complex, so Born-Oppenheimer approximation (BO approximation) is introduced. In BO approximation, the mass of the nucleus is far greater than that of electrons, and the velocity of the nucleus is far lower than that of electrons. Therefore, the nucleus and electrons can be separated from each other, and the time-independent electron motion equation (also called the time-independent Schrödinger equation) can be obtained as follows:

$$\hat{H}|\mathrm{\Psi }⟩=E|\mathrm{\Psi }⟩$$

$\hat{H}$ contains the following three items:

$$\hat{H}={\hat{K}}_e+{\hat{V}}_{ee}+{\hat{V}}_{Ne}$$

They are electron kinetic energy, electron-electron potential energy and electron-nuclear potential energy.

There are many numerical algorithms that can be used to solve the time-independent Schrödinger equation. This tutorial introduces one of these methods: the wave function. Wave function directly solves the eigenfunction and eigenenergy of a given molecular Hamiltonian. At present, there are a large number of open-source software packages, such as PySCF, which can be implemented. Here is a simple example: lithium hydride molecules, using the OpenFermion and OpenFermion-PySCF plug-ins. First, define the molecular structure:

dist = 1.5
geometry = [
    ["Li", [0.0, 0.0, 0.0 * dist]],
    ["H", [0.0, 0.0, 1.0 * dist]],
]
basis = "sto3g"
spin = 0
print("Geometry: \n", geometry)

Geometry:
 [['Li', [0.0, 0.0, 0.0]], ['H', [0.0, 0.0, 1.5]]]

The code above defines a Li-H key with a length of 1.5Å molecules. The STO-3G basis set is used for computing. Then, OpenFermion-PySCF is used to call PySCF to perform Hartree-Fock (HF), coupled-cluster with singles and doubles (CCSD), and full configuration interaction (FCI) computing. These three methods belong to the wave function. Before starting the computing, first make a brief introduction to these methods.

Wave Function

One of the methods to solve the time-independent Schrödinger equation is the Hartree-Fock (HF) method, which was proposed by Hartree et al. in the 1930s and is the basic method in quantum chemistry computing. The HF method introduces a single determinant approximation, that is, a wave function of the $N$-electronic system is represented by a wave function in a determinant form:

$$|\mathrm{\Psi }⟩=|{\psi }_1{\psi }_2{\psi }_3\dots {\psi }_N⟩$$

Where $|{\psi }_1{\psi }_2{\psi }_3\dots ⟩$ represents the Nth-order determinants formed by a set of spin-orbit wave functions $\{{\pi }_i\}$. The spin-orbit wave function ${\psi }_i$ may be further expanded with a set of basis functions in known forms:

$${\psi }_i={\varphi }_i{\eta }_i$$

$${\varphi }_i=\sum _{\mu }{C}_{\mu i}{\chi }_{\mu }$$

$\{{\chi }_{\mu }\}$ is referred to as a basis function, and may be a Gaussian function or the like. This approximation considers the exchange between electrons, but neglects the correlation between electrons, so it cannot correctly compute the properties such as dissociation energy.

The improvement of the HF method can be derived from the wave function expansion theorem. The wave function expansion theorem can be expressed as follows: if $\{{\psi }_i\}$ is a complete set of spin-orbit wave functions, the $N$-electronic system wave function may be accurately expanded by a determinant wave function formed by $\{{\psi }_i\}$:

$$|\mathrm{\Psi }⟩=\sum _{i_1<i_2<\cdots <i_N}^{\mathrm{\infty }}{C}_{i_1{i}_2\dots i_N}|{\psi }_{i_1}{\psi }_{i_2}\dots {\psi }_{i_N}⟩$$

You can obtain the configuration interaction (CI) method:

$$|{\mathrm{\Psi }}_{CI}⟩=C_0|{\mathrm{\Psi }}_{HF}⟩+\sum _{i\in occ,a\notin occ}^{a\to \mathrm{\infty }}{C}_i^a|{\mathrm{\Psi }}_i^a⟩+\sum _{ij\in occ,ab\notin occ}^{ab\to \mathrm{\infty }}{C}_{ij}^{ab}|{\mathrm{\Psi }}_{ij}^{ab}⟩$$

$|{\mathrm{\Psi }}_i^a⟩+\dots$ in the preceding formula represents a single excitation wave function of the electron from the orbit $i$ to the orbit $a$, and so on. A CI that considers only single excitation and double excitation is called a configuration interaction with singles and doubles (CISD). The CI that takes into account all the ground-state HF wave functions to N excitation wave functions is called full configuration interaction (FCI). The FCI wave function is the exact solution of the time-independent Schrödinger equation under the given basis function.

Second Quantization

Under the second quantization expression, the Hamiltonian of the system has the following form:

$$\hat{H}=\sum _{p,q}{h}_q^p{E}_q^p+\sum _{p,q,r,s}\frac{1}{2}{g}_{rs}^{pq}{E}_{rs}^{pq}$$

$E_q^p$ and $E_{rs}^{pq}$ are as follows:

$$\begin{array}{r}\begin{array}{rl}{E}_{rs}^{pq}& =a_p^{†}{a}_q^{†}{a}_r{a}_s\\ E_q^p& =a_p^{†}{a}_q\end{array}\end{array}$$

$a_p^{†}$ and $a_q$ are creation operator and annihilation operator, respectively.

The excited-state wave function can be expressed conveniently by using a second quantization expression method:

$$|{\mathrm{\Psi }}_{ijk\dots }^{abc\dots }⟩=a_a^{†}{a}_b^{†}{a}_c^{†}\dots a_i{a}_j{a}_k\dots |\mathrm{\Psi }⟩$$

An improvement to the CI method is the coupled-cluster theory (CC). Exponential operators are introduced to CC:

$$|{\mathrm{\Psi }}_{CC}⟩=\exp (\hat{T})|{\mathrm{\Psi }}_{HF}⟩$$

The coupled-cluster operator $\hat{T}$ is the sum of excitation operators.

$$\hat{T}=\sum _{p\notin occ,q\in occ}{\theta }_q^p{E}_q^p+\sum _{pq\notin occ,rs\in occ}{\theta }_{rs}^{pq}{E}_{rs}^{pq}+\dots$$

$\theta$ is similar to $C$ in the CI method, and is the parameter to be solved. It is easy to know from the Taylor’s expansion of the exponent that even if the coupled-cluster operator $\hat{T}$ includes only a low-order excitation term, the $\exp (\hat{T})$ may also implicitly include a high-order excitation term. This also makes a convergence speed of the CC method to the FCI wave function much faster than that of the CI. For truncation to K excitation, for example, K=2, the accuracy of the CCSD exceeds that of the CISD.

The effect of electron correlation is to reduce the total energy, so the ground state energy of HF is slightly higher than that of CCSD and FCI. In addition, it is easy to find that the computing volume of FCI is much greater than that of CCSD and HF. The MolecularData function encapsulated by OpenFermion and the run_pyscf function encapsulated by OpenFermion-PySCF are used for demonstration.

molecule_of = MolecularData(
    geometry,
    basis,
    multiplicity=2 * spin + 1,
    data_directory='./'
)
molecule_of = run_pyscf(
    molecule_of,
    run_scf=1,
    run_ccsd=1,
    run_fci=1
)

print("Hartree-Fock energy: %20.16f Ha" % (molecule_of.hf_energy))
print("CCSD energy: %20.16f Ha" % (molecule_of.ccsd_energy))
print("FCI energy: %20.16f Ha" % (molecule_of.fci_energy))

Hartree-Fock energy:  -7.8633576215351182 Ha
CCSD energy:  -7.8823529091526794 Ha
FCI energy:  -7.8823622867987284 Ha

In the preceding example, HF, CCSD, and FCI are used to compute the total energy. If you collect statistics on the runtime, you will find that $T_{HF}<T_{CCSD}\ll T_{FCI}$. It is more obvious if you use the system with larger calculation amount, such as ethylene molecule. In addition, for the total computed energy, you will find that $E_{HF}>E_{CCSD}>E_{FCI}$. After the computing is complete, save the result to the molecule_file file (molecule_of.filename).

molecule_of.save()
molecule_file = molecule_of.filename
print(molecule_file.split('/')[-1])

H1-Li1_sto3g_singlet

One of the major obstacles to quantum chemistry is the volume of computation. As the system size (electron number and atomic number) increases, the time required for solving the FCI wave function and ground state energy increases by about $2^N$. Even for small molecules such as ethylene molecules, FCI computing is not easy. Quantum computers provide a possible solution to this problem. Research shows that quantum computers can simulate the time-dependent evolution of Hamiltonian in terms of polynomial time complexity. Compared with classical computers, quantum computers exponentially accelerate the chemical simulation on quantum processors. This tutorial introduces one of the quantum algorithms: VQE.

Variational Quantum Eigensolver (VQE)

The VQE is a hybrid quantum-classical algorithm. It uses the variational principle to solve the ground state wave function. The optimization of variational parameters is carried out on the classical computer.

Variational Principle

The variational principle may be expressed in the following form:

$$E_0\le \frac{⟨{\mathrm{\Psi }}_t|\hat{H}|{\mathrm{\Psi }}_t⟩}{⟨{\mathrm{\Psi }}_t|{\mathrm{\Psi }}_t⟩}$$

In the preceding formula, $|{\mathrm{\Psi }}_t⟩$ indicates the probe wave function. The variational principle shows that the ground state energy obtained by any probe wave function is always greater than or equal to the real ground state energy under certain conditions. The variational principle provides a method for solving the molecular ground state Schrödinger equation. A parameterized function $f(\theta )$ is used as an approximation of the accurate ground state wave function, and the accurate ground state energy is approximated by optimizing the parameter $\theta$.

Initial State Preparation

The $N$-electron HF wave function also has a very concise form under the quadratic quantization expression:

$$|{\mathrm{\Psi }}_{HF}⟩=\prod _{N-1}^{i=0}{a}_i^{†}|0⟩$$

The above formula builds a bridge from quantum chemical wave function to quantum computing: $|0⟩$ is used to represent a non-occupied orbit, and $|1⟩$ is used to represent an orbit occupied by an electron. Therefore, the $N$-electron HF wave function may be mapped to a string of $M+N$ quantum bits $|00\dots 11\dots ⟩$. $M$ indicates the number of unoccupied tracks.

The following code builds an HF initial state wave function corresponding to the LiH molecule. In Jordan-Wigner transformation, $N$ $\text{X}$ gates are applied to $|000\dots ⟩$.

hartreefock_wfn_circuit = Circuit([X.on(i) for i in range(molecule_of.n_electrons)])
print(hartreefock_wfn_circuit)

      ┏━━━┓
q0: ──┨╺╋╸┠───
      ┗━━━┛
      ┏━━━┓
q1: ──┨╺╋╸┠───
      ┗━━━┛
      ┏━━━┓
q2: ──┨╺╋╸┠───
      ┗━━━┛
      ┏━━━┓
q3: ──┨╺╋╸┠───
      ┗━━━┛

We can build a probe wave function in the following form:

$$|{\mathrm{\Psi }}_t⟩=U(\theta )|{\mathrm{\Psi }}_{HF}⟩$$

$U(\theta )$ represents a unitary transformation that may be simulated by using a quantum circuit. $|{\mathrm{\Psi }}_{HF}⟩$ is used as an initial state, and may be conveniently prepared by using a plurality of single-bit $\text{X}$ gates. A specific form of the $U(\theta )|{\mathrm{\Psi }}_{HF}⟩$ is also referred to as wave function ansatz.

Wave Function Ansatz

The coupled-cluster theory mentioned above is a very efficient wave function ansatz. To use it on a quantum computer, you need to make the following modifications:

$$|{\mathrm{\Psi }}_{UCC}⟩=\exp (\hat{T}-{\hat{T}}^{†})|{\mathrm{\Psi }}_{HF}⟩$$

UCC is short for unitary coupled-cluster theory. ${\hat{T}}^{†}$ represents the Hermite conjugate of $\hat{T}$. In this way, $\exp (\hat{T}-{\hat{T}}^{†})$ is the unitary operator. Peruzzo et al. first performed chemical simulation experiments on quantum computers using VQE and unitary coupled-cluster with singles and doubles (UCCSD) in 2014. It should be noted that, by default, the parameter $\{\theta \}$ in the coupled-cluster operator is a real number. There is no problem with this hypothesis in molecular systems. In periodic systems, the study of Liu Jie et al. suggests that a unitary coupled-cluster can result in errors due to the neglect of the complex numbers. This tutorial does not discuss the application of unitary coupled-cluster in periodic systems.

The generate_uccsd function in the circuit module of MindSpore Quantum can be used to read the computing result saved in molecule_file, build the UCCSD wave function by one click, and obtain the corresponding quantum circuit.

ansatz_circuit, \
init_amplitudes, \
ansatz_parameter_names, \
hamiltonian_QubitOp, \
n_qubits, n_electrons = generate_uccsd(molecule_file, threshold=-1)

ccsd:-7.882352909152679.
fci:-7.882362286798728.

generate_uccsd packs functions related to the unitary coupled-cluster, including multiple steps such as deriving a molecular Hamiltonian, building a unitary coupled-cluster ansatz operator, and extracting a coupled-cluster coefficient computed by CCSD. This function reads the molecule by entering its file path. The parameter th indicates the to-be-updated gradient threshold of a parameter in the quantum circuit. In the section Building a Unitary Coupled-Cluster Ansatz Step by Step, we will demonstrate how to use the related interfaces of MindSpore Quantum to complete the steps. A complete quantum circuit includes an HF initial state and a UCCSD ansatz, as shown in the following code:

total_circuit = hartreefock_wfn_circuit + ansatz_circuit
total_circuit.summary()
print("Number of parameters: %d" % (len(ansatz_parameter_names)))

                            Circuit Summary                            
╭──────────────────────┬──────────────────────────────────────────────╮
│ Info                 │ value                                        │
├──────────────────────┼──────────────────────────────────────────────┤
│ Number of qubit      │ 12                                           │
├──────────────────────┼──────────────────────────────────────────────┤
│ Total number of gate │ 12612                                        │
│ Barrier              │ 2560                                         │
│ Noise Channel        │ 0                                            │
│ Measurement          │ 0                                            │
├──────────────────────┼──────────────────────────────────────────────┤
│ Parameter gate       │ 640                                          │
│ 44 ansatz parameters │ p0, p8, p1, p9, p2, p10, p3, p11, p4, p12... │
╰──────────────────────┴──────────────────────────────────────────────╯

Number of parameters: 44

For the LiH molecule, the UCCSD wave function ansatz includes 44 variational parameters. The total number of quantum bit gates of the circuit is 12612, and a total of 12 quantum bits are needed for simulation.

VQE Procedure

The procedure for solving the molecular ground state by using the VQE is as follows:

1. Prepare the HF initial state: $|00\dots 11\dots ⟩$.

2. Define the wave function ansatz, such as UCCSD.

3. Convert the wave function into a variational quantum circuit.

4. Initialize the variational parameters, for example, set all parameters to 0.

5. Obtain the energy $E(\theta )$ of the molecular Hamiltonian under the set of variational parameters and the derivative $\{\partial E/\partial {\theta }_i\}$ of the energy about the parameters by means of multiple measurements on the quantum computer.

6. Use optimization algorithms, such as gradient descent and BFGS, to update variational parameters on classical computers.

7. Transfer the new variational parameters to the quantum circuit for updating.

8. Repeat steps 5 to 7 until the convergence criteria are met.

9. End.

In step 5, the derivative $\{\partial E/\partial {\theta }_i\}$ of the energy about the parameter may be computed by using a parameter-shift rule on a quantum computer, or may be computed by simulating a parameter-shift rule or a finite difference method in a simulator. This is a relatively time-consuming process. Based on the MindSpore framework, MindSpore Quantum provides the automatic derivation function similar to machine learning, which can efficiently compute the derivatives of variational quantum circuits in simulation. The following uses MindSpore Quantum to build a parameterized UCCSD quantum circuit with an automatic derivation function:

sim = Simulator('mqvector', total_circuit.n_qubits)
molecule_pqc = sim.get_expectation_with_grad(Hamiltonian(hamiltonian_QubitOp), total_circuit)

You can obtain the energy $E(\theta )=⟨{\mathrm{\Psi }}_{UCC}(\theta )|\hat{H}|{\mathrm{\Psi }}_{UCC}(\theta )⟩$ corresponding to the variational parameter and the derivative of each variational parameter by transferring a specific value of the parameter to molecule_pqc.

For example, we can use the following code to get the expectation of hamiltonian and the corresponding gradient when initial parameters of variational quantum circuit is zero.

import numpy as np

n_params = len(total_circuit.params_name)
p0 = np.zeros(n_params)
f, g = molecule_pqc(p0)
print("Energy: ", f, "\nshape: ", f.shape, '\n')
print("Gradient: ", g, "\nshape: ", g.shape)

Energy:  [[-7.86335762+0.j]]
shape:  (1, 1)

Gradient:  [[[-5.76761758e-11+0.j -8.60518134e-02+0.j  1.19275904e-08+0.j
   -4.85545124e-02+0.j -1.24342308e-14+0.j -3.92769093e-02+0.j
   -3.28108247e-15+0.j -9.59481745e-02+0.j  1.24424913e-14+0.j
   -3.92769093e-02+0.j -3.82408561e-15+0.j -9.59481745e-02+0.j
   -8.60105773e-11+0.j -2.89649669e-02+0.j  2.97699096e-10+0.j
   -4.91813233e-01+0.j -9.35293242e-04+0.j -1.63755737e-16+0.j
   -3.35024598e-17+0.j  1.55629243e-17+0.j  1.59462515e-16+0.j
    2.67655514e-17+0.j -2.16577662e-17+0.j  5.06813117e-03+0.j
    1.08542346e-02+0.j -1.28614265e-02+0.j -4.90653908e-17+0.j
    1.71728860e-17+0.j  1.33973787e-01+0.j -3.03063680e-02+0.j
   -3.37324552e-18+0.j  9.19976050e-19+0.j -1.05649122e-28+0.j
    5.33586109e-17+0.j -3.68357159e-17+0.j  3.33490478e-17+0.j
    1.09202635e-28+0.j  1.15391965e-16+0.j -3.03063680e-02+0.j
   -5.27452935e-17+0.j  4.81071347e-17+0.j -4.22508288e-17+0.j
   -2.44510199e-17+0.j -1.68035039e-03+0.j]]]
shape:  (1, 1, 44)

Throw the above calculation, we get the energy and gradient value and the user can use these data according their practical needs. Now we following the step of (5)~(7) of optimization of VQE, to optimize of variational quantum circuit. Here we use the optimizer in scipy to optimize our quantum circuit. First, we need to define the optimization function that scipy required:

def fun(p0, molecule_pqc, energy_list=None):
    f, g = molecule_pqc(p0)
    f = np.real(f)[0, 0]
    g = np.real(g)[0, 0]
    if energy_list is not None:
        energy_list.append(f)
        if len(energy_list) % 5 == 0:
            print(f"Step: {len(energy_list)},\tenergy: {f}")
    return f, g

fun(p0, molecule_pqc)

(-7.863357621536964,
 array([-5.76761758e-11, -8.60518134e-02,  1.19275904e-08, -4.85545124e-02,
        -1.24342308e-14, -3.92769093e-02, -3.28108247e-15, -9.59481745e-02,
         1.24424913e-14, -3.92769093e-02, -3.82408561e-15, -9.59481745e-02,
        -8.60105773e-11, -2.89649669e-02,  2.97699096e-10, -4.91813233e-01,
        -9.35293242e-04, -1.63755737e-16, -3.35024598e-17,  1.55629243e-17,
         1.59462515e-16,  2.67655514e-17, -2.16577662e-17,  5.06813117e-03,
         1.08542346e-02, -1.28614265e-02, -4.90653908e-17,  1.71728860e-17,
         1.33973787e-01, -3.03063680e-02, -3.37324552e-18,  9.19976050e-19,
        -1.05649122e-28,  5.33586109e-17, -3.68357159e-17,  3.33490478e-17,
         1.09202635e-28,  1.15391965e-16, -3.03063680e-02, -5.27452935e-17,
         4.81071347e-17, -4.22508288e-17, -2.44510199e-17, -1.68035039e-03]))

Here, the fun that we define can correctly to return the data that we need: a real energy value, and a array of gradient value with the same size of parameters. Now, we use bfgs optimizer in scipy to finish the optimization.

from scipy.optimize import minimize

energy_list = []
res = minimize(fun, p0, args=(molecule_pqc, energy_list), method='bfgs', jac=True)

Step: 5,        energy: -7.880227726111826
Step: 10,       energy: -7.8818171240648365
Step: 15,       energy: -7.882213242905283
Step: 20,       energy: -7.8823453369936445
Step: 25,       energy: -7.8823524949903465
Step: 30,       energy: -7.8823526912721755
Step: 35,       energy: -7.882352703403696
Step: 40,       energy: -7.882352708341612

So, we finished the gradient optimization of variational quantum circuit. Here, energy_list is going to store the energy during optimization. Here, we briefly introduce the usage of minimize:

• fun: The first arg is the function you want to optimize.

• p0: The second arg is the initial value of variables.

• args: The other argument of fun except the first argument. According the definition of fun, we choose args=(molecule_pqc, energy_list).

• method: The optimization algorithm. Here we use a second order optimization algorithm bfgs. For more optimization algorithm, please refer: scipy tutorial.

• jac: To info that whether fun return gradient. Here we use True, because MindSpore Quantum can calculate the accuracy gradient value of variational quantum circuit. If use False, minimize framework will calculate the approximated gradient value base on difference method.

res is the optimization result of scipy, including optimized parameters, the optimized value and evolution steps.

print(f"Ground state: \n{res.fun}\n")
print(f"FCI: \n-7.882362286798721\n")
print(f"Optimized amplitudes: \n{res.x}")

Ground state:
-7.882352708341612

FCI:
-7.882362286798721

Optimized amplitudes:
[ 2.38815317e-04  1.89073105e-03  3.52372251e-02  1.60368559e-02
 -4.30846145e-09  9.09437670e-04 -5.45444468e-10  1.41627870e-02
  3.62569484e-09  9.08701303e-04 -1.05938178e-10  1.41711693e-02
 -5.47703644e-04  4.26818370e-04  2.87168405e-03  5.38109720e-02
  2.34920831e-04 -9.88498886e-08 -2.51432020e-07  2.54936335e-07
 -2.65409411e-08  4.64583011e-08 -4.42267421e-08  1.32809550e-05
 -1.04167881e-04  7.98983208e-04  5.39567507e-10 -2.28346700e-10
 -5.50005419e-02  3.09112487e-03 -4.88778663e-10 -4.89650670e-09
  4.92879467e-09  1.13389943e-07  4.45780819e-09 -1.85315776e-08
 -1.42394064e-12  4.57319777e-10  3.09100040e-03  5.36042785e-08
 -3.87192856e-09  9.99056234e-10  1.51453640e-10  3.72805204e-04]

We can see here the result of ucc method is very close to FCI method with very good accuracy.

Building a Unitary Coupled-Cluster Ansatz Step by Step

In the preceding part, the generate_uccsd is used to build all the content required for designing the unitary coupled-cluster. In this section, the steps are split, we get the coupled-cluster operator, the corresponding quantum circuit and the initial guess of the variational parameters from the classical CCSD results. First, import some extra dependencies, including the related functions in MindSpore Quantum.

from mindquantum.algorithm.nisq import Transform
from mindquantum.algorithm.nisq import get_qubit_hamiltonian
from mindquantum.algorithm.nisq import uccsd_singlet_generator, uccsd_singlet_get_packed_amplitudes
from mindquantum.core.operators import TimeEvolution

The molecule Hamiltonian uses get_qubit_hamiltonian to read the previous computing result. The result is as follows:

hamiltonian_QubitOp = get_qubit_hamiltonian(molecule_of)

The unitary coupled-cluster operator $\hat{T}-{\hat{T}}^{†}$ can be built using uccsd_singlet_generator. Provide the total number of quantum bits (total number of spin orbits) and the total number of electrons, and set anti_hermitian=True.

ucc_fermion_ops = uccsd_singlet_generator(
    molecule_of.n_qubits, molecule_of.n_electrons, anti_hermitian=True)

The ucc_fermion_ops built in the previous step is parameterized. Use the Jordan-Wigner transformation to map the Fermi excitation operator to the Pauli operator:

ucc_qubit_ops = Transform(ucc_fermion_ops).jordan_wigner()

Next, we need to obtain the quantum circuit corresponding to the unitary operator $\exp (\hat{T}-{\hat{T}}^{†})$. TimeEvolution can generate the circuit corresponding to $\exp (-i\hat{H}t)$, where $\hat{H}$ is a Hermitian operator, and $t$ is a real number. Note that when TimeEvolution is used, ucc_qubit_ops already contains the complex number factor $i$. Therefore, you need to divide ucc_qubit_ops by $i$ or extract the imaginary part of ucc_qubit_ops.

ansatz_circuit = TimeEvolution(ucc_qubit_ops.imag, 1.0).circuit
ansatz_parameter_names = ansatz_circuit.params_name

ansatz_parameter_names is used to record the parameter names in the circuit. So far, we have obtained the contents required by the VQE quantum circuit, including the Hamiltonian hamiltonian_QubitOp and the parameterized wave function ansatz ansatz_circuit. By referring to the preceding steps, we can obtain a complete state preparation circuit. hartreefock_wfn_circuit mentioned above is used as the Hartree-Fock reference state:

total_circuit = hartreefock_wfn_circuit + ansatz_circuit
total_circuit.summary()

                                  Circuit Summary                                  
╭──────────────────────┬──────────────────────────────────────────────────────────╮
│ Info                 │ value                                                    │
├──────────────────────┼──────────────────────────────────────────────────────────┤
│ Number of qubit      │ 12                                                       │
├──────────────────────┼──────────────────────────────────────────────────────────┤
│ Total number of gate │ 12612                                                    │
│ Barrier              │ 2560                                                     │
│ Noise Channel        │ 0                                                        │
│ Measurement          │ 0                                                        │
├──────────────────────┼──────────────────────────────────────────────────────────┤
│ Parameter gate       │ 640                                                      │
│ 44 ansatz parameters │ s_0, d1_0, s_1, d1_1, s_2, d1_2, s_3, d1_3, s_4, d1_4... │
╰──────────────────────┴──────────────────────────────────────────────────────────╯

Next, you need to provide a reasonable initial value for the variational parameter. The PQCNet built in the preceding text uses 0 as the initial guess by default, which is feasible in most cases. However, using CCSD’s computational data as a starting point for UCC may be better. Use the uccsd_singlet_get_packed_amplitudes function to extract CCSD parameters from molecule_of.

init_amplitudes_ccsd = uccsd_singlet_get_packed_amplitudes(
    molecule_of.ccsd_single_amps, molecule_of.ccsd_double_amps, molecule_of.n_qubits, molecule_of.n_electrons)
init_amplitudes_ccsd = [init_amplitudes_ccsd[param_i] for param_i in ansatz_parameter_names]

Just like the previous method, we can get the grad_ops with MindSpore Quantum, and optimize it with scipy.

grad_ops = Simulator('mqvector', total_circuit.n_qubits).get_expectation_with_grad(
    Hamiltonian(hamiltonian_QubitOp.real),
    total_circuit)

init_amplitudes_ccsd (coupled-cluster coefficient computed by CCSD) is used as an initial variational parameter:

energy_list = []
res = minimize(fun, init_amplitudes_ccsd, args=(grad_ops, energy_list), method='bfgs', jac=True)

Step: 5,        energy: -7.878223283110822
Step: 10,       energy: -7.880288481802435
Step: 15,       energy: -7.8820356683191415
Step: 20,       energy: -7.882302370884269
Step: 25,       energy: -7.882349803535784
Step: 30,       energy: -7.882352702053716
Step: 35,       energy: -7.88235270792034
Step: 40,       energy: -7.882352708346457

The final optimized result is shown as below.

print(f"Ground state: \n{res.fun}\n")
print(f"FCI: \n-7.882362286798721\n")
print(f"Optimized amplitudes: \n{res.x}")

Ground state:
-7.882352708346457

FCI:
-7.882362286798721

Optimized amplitudes:
[-2.38636816e-04  1.89071890e-03 -3.52372349e-02  1.60368129e-02
  1.34665623e-08  9.09430837e-04 -4.17115460e-10  1.41641317e-02
 -1.28144111e-08  9.08694024e-04  4.20204892e-10  1.41698036e-02
  5.47710797e-04  4.26820871e-04 -2.87169433e-03  5.38109417e-02
  2.34627619e-04  4.32417990e-07 -1.26325933e-08 -1.39715964e-08
 -4.27213491e-07  1.24814811e-08  1.37984739e-08  1.32946010e-05
  7.99032002e-04 -1.04121599e-04  1.25575236e-09 -1.24558942e-09
 -5.50005457e-02  3.09142613e-03  1.53511697e-08 -3.40030804e-08
  3.29912584e-08  3.67660373e-07 -2.59584687e-08  2.22720149e-10
  4.47254154e-11  1.06718989e-09  3.09086800e-03 -3.54929776e-07
  1.95048358e-08  5.27817957e-09 -1.03614664e-09  3.72830373e-04]

Summary

In this case, the ground state energy of the LiH molecule is obtained by using scipy in two methods. In the first method, we use the generate_uccsd function packaged by MindSpore Quantum to generate a quantum neural network that can solve this problem. In the second method, we build a similar gradient operator step by step. The final results are the same.

from mindquantum.utils.show_info import InfoTable

InfoTable('mindquantum', 'scipy', 'numpy')

Software	Version
mindquantum	0.9.11
scipy	1.10.1
numpy	1.23.5
System	Info
Python	3.9.16
OS	Linux x86_64
Memory	8.3 GB
CPU Max Thread	8
Date	Mon Jan 1 00:25:57 2024