Quantum Chemistry: Simulation of Lithium Cations

Published by strangeworks (01/22/2021)

Written by Bikash's Quantum for the Strangeworks Platform

Molecular simulation, one of the promising applications of a quantum computer, is performed mostly on neutral species to calculate their ground state energies. The available algorithms are not suitable for the charged species as they usually tend to converge the energy to one of the neutral systems. Here, we show the quantum simulation of both dilithium and lithium cations on the quantum simulators. Li2+\text{Li}_{2}^+​ can be simulated using 12, 14 and 20 qubits and Li+\text{Li}^+​ using 10, 8 and 4 qubits. Simulations are performed in Qiskit using the variational quantum eigensolver algorithm. Different variational forms can be used including UCCSD, RYRZ and Swap RZ optimization techniques. Considering the active space of the system, we can reduce the number of molecular orbitals and simulate the system using 12 and 14 qubits, thus obtaining a good approximation to the energy value.


Following Feynman's idea, a number of works have been proposed in the field of quantum computation. Quantum simulation has become a field of great interest catching the eyes of researchers working in the field of condensed matter, nuclear physics, chemistry, biology, cosmology and so on. Quantum chemistry being an area for the rigorous application of quantum computation, algorithms for the efficient simulation of molecular systems have been proposed. These are limited to the neutral molecules as the optimization routines tend to converge the energy to the ground state of the nearest neutral molecular configuration. Charged atomic or molecular configurations are yet to be addressed. It brings a great zeal of interest to see how we can approximate the ground and excited state energies for the charged systems with the available algorithms, different encoding methods, and quantum chemistry programs. Variational Quantum Eigensolver (VQE) is best suited for creating wavefunctions for molecules. For a large number of qubit requirement and a charged molecule or atom, the algorithm does not provide good approximations to the energy. We attempt simulating a lesser-known charged molecular system dilithium cation and the well known lithium cation. Li2+\text{Li}_2^{+}​being a 20 qubit system which can be further reduced to 12 and 14 qubits by considering only the active space of the molecule. The 20-qubit circuit optimization is found to give poor approximations and thus the qubits are reduced. Li+\text{Li}^+ finds its applications particularly in medicine and batteries, and thus remains an important system.

Structure of the Li2+\text{Li}_{2}^{+}​ and Li+\text{Li}^{+}​. Red dots in the center are nucleus, holes depict the vacant positions and black dots are electrons.

The dilithium cation (Li2+\text{Li}_{2}^{+}​) is a seven-body system with 5 electrons and two nuclei and lithium ion (Li+\text{Li}^{+}​) is a three-body system with 2 electrons and 1 nucleus. The structure of Li2+\text{Li}_{2}^{+}​ and Li+\text{Li}^{+}​ ion is shown in Fig. 1.

The trial guess wavefunction for these systems is generated by VQE algorithm. To get a good approximation for Li2+\text{Li}_{2}^{+}​, we reduce the active space. The importance of VQE lies on the fact that it is able to prepare the ansatz by varying the parameters of the circuit. Due to it's high fidelity, robustness to errors and less resource requirement, it is best suited for molecular simulations. The backbone of VQE is the variational principle,

∣⟨ψ∣H∣ψ⟩∣∣⟨ψ∣ψ⟩∣≥E\frac{|\langle \psi| H | \psi \rangle |}{|\langle \psi | \psi \rangle|} \geq E \hspace{1cm}

where E is the minimum possible energy eigenvalue. The variational principle ensures that this expectation value is always greater than the smallest eigenvalue of H. Here, we show the calculation of the ground state energy of Li2+\text{Li}_{2}^{+}​ and Li+\text{Li}^{+}​ ion using VQE.


The Second-Quantization

We know the Schrodinger equation, 

H^∣ψ⟩=E∣ψ⟩\hat{H}| \psi \rangle = E | \psi \rangle

Where H^\hat{H}​ is Hamiltonian of the system and EE​ is the energy eigenvalue. This Hamiltonian consists of two parts, the kinetic energy term T^\hat{T}​ and the potential energy term U^\hat{U}​.

H^=T^+U^\hat{H} = \hat{T}+ \hat{U}​​

For the calculation of ground state or excited state energy of any system, we need to solve the Eq. 2. Here we deal with ions. In this case, the Hamiltonian describes a system of electrons and nuclei in an amalgamation of potential comprised by the Coulombic potential of each particle. To reduce the complexity of the calculation, we consider the Born-Oppenheimer (BO) approximation. As the mass of nucleus is much greater than the electron, the nucleus is treated as a stationary particle in BO approximation. Our reduced Hamiltonian is given by 

H^=−∑i∇i22−∑i,IZI∣ri−RI∣+∑i≠I12∣ri−RI∣\hat{H} = - \sum_{i} \frac{\nabla _{i}^{2}}{2} -\sum_{i, I} \frac{Z _{I}}{|r_{i} - R_{I}|} + \sum_{i \neq I} \frac{1}{2|r_{i} - R_{I}|}

Where indices ii​ and jj​ are used for electrons, index II​ is used for nucleus, RIR_{I}​ is the location of IthI_{th}​ nucleus,ZIZ_{I}​ is the atomic number of IthI_{th}​ nucleus, rir_{i}​  denotes the position of ithi_{th}​ electron. 

The first term in Eq. 4 denotes the total kinetic energy of the electrons, the second term denotes the Coulomb interaction between electrons and fixed nuclei and the third term represents the Coulomb interaction between the two different electrons.

The mapping of the Hamiltonian to Pauli matrices is done using the second quantization, introduced by Paul Dirac in 1927. The concept of the second quantization is useful to describe and analyz quantum many-body system. The most import thing about this formalism is the introduction of anti-commuting creation (a†a^{\dagger}​) and annihilation (aa​) operators. Here, ai†a_{i}^{\dagger}​ and aja_{j}​ follow the anti-commutation relations given by,

{ai†,aj}=δij,{ai,aj}=0,{ai†,aj†}=0\{a_i^{\dagger},a_j\} = \delta_{ij}, \{a_i,a_j\} = 0, \{a_i^{\dagger},a_j^{\dagger}\} = 0

The fermionic creation operator increases the occupational number of an orbital by one and the annihilation operator decreases it by one. Using the Jordan-Wigner or Bravyi-Kitaev or parity transform, the fermionic Hamiltonian can be mapped onto the spin operators. The second quantized form of the Hamiltonian for fermions is given by,

H=∑i,jhijai†aj+12∑i,j,k,lhijklai†aj†akalH= \sum_{i,j}h_{ij} a_i^{\dagger}a_j +\frac{1}{2}\sum_{i,j,k,l}h_{ijkl} a_i^{\dagger}a_j^{\dagger}a_k a_l

where hijh_{ij}​ and hijklh_{ijkl} are one and two electron integrals are given by

hij=dr1χi(r1)(122σZr1 Ro)χj(r1)h_{ij} = \int d\vec{r_1} \chi_i^*(\vec{r_1}) (\frac{- \nabla_1^{2}}{2} -{\sum_ \sigma\frac{Z}{|\vec{r_1}\ -R_o|}})\chi_j(\vec{r_1})


hijkl=dr1dr2χi(r1)χj(r2)(r1r2)χk(r2)χl(r1)2h_{ijkl} = \int \frac{ d\vec{r_1} d\vec{r_2} \chi_i^*(\vec{r_1}) \chi_j^*(\vec{r_2})}{|(r_1-r_2)|} \frac {{\chi_{k}(\vec{r_2})}{\chi_l(\vec{r_1})}}{2}

where χi(r1⃗){\chi_i(\vec{r_1})}​ is the ii​th spin orbital of the electron 1, ri⃗{\vec{r_i}} is the position of the ithi^{\text {th}}​ electron, r12r_{12}​ is the distance between the two points r1r_{1}​ and r2r_{2}​, Z is the total charge on the nucleus, and RoR_{o}​ is the position of the nucleus. 

The Hamiltonians of the Li+\text{Li}^{+}​ and Li2+Li_{2}^{+}​ are determined using the PyQuante on QISKit. Then it is converted into Pauli's spin operators form using the Jordan-Wigner transformation. The calculated Hamiltonian in Pauli form of Li+\text{Li}^+​ for four qubits is 

HLi+=−4.391IIII+1.588IIIZ+0.1648IIZI+1.588IZII+0.165ZIII+0.097IIZZ+0.4155IZIZ+0.00243XXYY+0.00243YYYY+0.00243XXXX+0.00243YYXX+0.099ZIIZ+0.099IZZI+0.078ZIZI+0.0967ZZIIH_{Li^+} = -4.391 IIII + 1.588 IIIZ + 0.1648 IIZI+1.588 IZII +0.165 ZIII + 0.097 IIZZ + 0.4155 IZIZ + 0.00243 XXYY + 0.00243 YYYY + 0.00243 XXXX+ 0.00243 YYXX + 0.099 ZIIZ + 0.099 IZZI + 0.078 ZIZI + 0.0967 ZZII

The Hamiltonian of Li2+\text{Li}_{2}^{+}​ contains nearly 3000 terms and all the other Hamiltonian used here also contains hundreds of terms.  

Jordan-Wigner transformation

Jordan-Wigner transformation is a method for encoding second quantized Hamiltonian, in the form of fermionic creation and annihilation operator, into spin operators. Then, it can be easily mapped on qubits. This method requires n qubits to store the occupation number of n spin-orbitals. An occupied orbital is equivalent to the qubit in state ∣1⟩|1 \rangle​ and an unoccupied orbital is equivalent to the qubit in state ∣0⟩|0 \rangle​. This new kind of basis is known as the occupation basis. The process of mapping is,

aj=Q⊗Zj−1⊗Zj−2....⊗Z0,aj†=Q†⊗Zj−1⊗Zj−2....⊗Z0a{_j} = Q{\otimes}Z{_{j-1}}{\otimes}Z{_{j-2}}....{\otimes}Z{_0},a_{j}^{\dagger} = Q^{\dagger}{\otimes}Z{_{j-1}}{\otimes}Z{_{j-2}}....{\otimes}Z{_0}


Q=X+iY2Q = \frac{X + iY}{2}

Q†=X−iY2Q{^\dagger} = \frac{X - iY}{2}

Where X, Y and Z are Pauli's spin operators. Each qubit stores the occupation number of the orbital in Fock space. Other encoding methods like Bravyi-Kitaev encoding and parity encoding also could be used for the Hamiltonian transformation. In Bravyi-Kitaev model, the qubit with even index stores the occupation number of orbitals, as in Jordan-Wigner transformation, and the qubits with odd index stores the parity of a particular set of occupation numbers corresponding to that set of orbitals.

Variational Quantum Eigensolver

VQE is a semi Quantum/Classical and the most efficient hybrid quantum computing algorithm till date. We run a quantum subroutine and then optimize the eigenvalues using a classical optimization loop. VQE, first demonstrated in 2014 is widely used for molecular simulations. The VQE algorithm consists of two parts:


We prepare the initial state of the system ψ{\psi}​, also known as the ansatz and measure the expectation value ∣⟨ψ∣H∣ψ⟩∣{|\langle \psi| H | \psi \rangle |}​ using the quantum simulators. The variational principle ensures that the calculated expectation value is always greater than the smallest eigenvalue of H. The ansatz is prepared by entangling all the qubits with various single-qubit rotation gates. This ansatz preparation could be done by using different combination of entangler and rotation gates. The complexity of initial state increases as the size of the system (i.e. the number of qubits required for simulation) increases. For dilithium cation, we use pre-specified QISKit variational forms.


Each time after the calculation of expectation value, the classical optimizer automatically feeds it back to an optimization loop to obtain the minimum possible eigenvalue. Some easily available classical optimizers are Nelder-Mead, Powell, Coybala and gradient descent methods. These optimizers iterate until the energy converges. VQE is based on the Rayleigh-Ritz variational principle,

∣⟨ψ(θi)∣H∣ψ(θi)⟩∣≥E0{|\langle \psi ({\theta_i})| H | \psi ({\theta_i})\rangle |} \geq E_0

where E0E_0​ is the ground state energy. The wave function ∣ψ(θi)⟩|\psi(\theta_i)\rangle is a trial wave function. The quantum circuit consists of 3 parts:

Initial state preparation (ansatz)

This part creates the initial wave-function (ψ\psi​). In QISKit, the state preparation is done using different variational forms like RyRz, UCCSD, and SwapRz with a given initial state like HF (Hartree - Fock) state. These variational forms provide a different combination of entangled qubits with different unitary rotation gates. The circuit depth can be increased to produce more complex state as required. For an n-qubit system, the number of parameters for initial rotations would be 3n(n−1)3n(n-1)​ and the number of unitary gates required is given by n(n−1)n(n-1)​.

The Hamiltonian

As our Hamiltonian is in the form of Pauli X, Y and Z gates, we put each Hamiltonian term after ansatz circuit one by one. 


Finally, we measure the expectation values of each term in the Hamiltonian one by one. 

Schematic of quantum circuit used for VQE. Ansatz is created by using different combination of rotation gates with entanglement.  

Optimization Methods

The main classical part of the VQE algorithm is the optimization of the obtained result and it is also the last step for finding the minimum (ground state) energy. We choose the best-suited optimizer for our system. Different classical optimizers like particle swarm optimization, Nelder-Mead, Powell, L-BFGS-B, and Coybala, etc are used. These optimizers mainly come in two categories. The L-BFGS-B follows the gradient-based method, in which the optimization is done using the gradient of the objective function. Some other optimizers use direct search methods. For example, COYBALA optimizes by using the linear approximation not by using the gradient of the function and so on.

The gradient-based methods converge the function at a local minima, which is no use for us as we need global minima of the function. The direct search method is derivative free and good for the systems where we use noisy and discontinuous functions. However, these optimizers are system dependent, means the accuracy of the result will vary with the different number of qubits used (i.e. complexity of the system). Here we use COYBALA (Constrained Optimization BY Linear Approximation) optimizer for our energy calculation. Powell or Nelder-Mead can also be used.


The Variational Quantum Eigensolver is used for the simulation, which is a semi-quantum algorithm, with the quantum chemistry package PyQuante. All the calculations are done using the QISKit package.

To determine the ground state energy, we need the h values, which is calculated by PyQuante. Manually the h values can be calculated by solving the one and two-electron integrals. The expectation value of the Hamiltonian is therefore,

⟨ψ∣H∣ψ⟩=Σ0imaxhi⟨ψ∣Hi∣ψ⟩\langle \psi| H | \psi \rangle= \Sigma_{0} ^{i_{max}} h_{i} \langle \psi|{H_{i}}|\psi \rangle

Where H is the Hamiltonian of the system, ii​ is the index running from 0 to the maximum number of Hamiltonian terms, HiH_{i}​ denotes each Hamiltonian terms, hih_{i}​ is the hh​ values for each HiH_{i} respectively.

For initial geometry optimization of the ions, we perform Restricted Open-Shell Hartree-Fock (ROHF) calculation in STO-3G basis using the PyQuante package and obtain the Hamiltonian of the system. We use the Jordan-Wigner transformation to convert the Hamiltonian in spin operators form.  Then we use the optimized initial HF state as the ansatz (∣ψ(θ)⟩| \psi (\theta) \rangle​) for the system. 

We reduce the active space of the ions by manually selecting the active orbitals and calculate the Hamiltonian for the different number of qubits. Which we further fed into the quantum system. 

We use different variational forms with the initial HF state and Hamiltonian of the system and start our quantum machine.

We calculate the expectation value of the Hamiltonian (i.e. ⟨ψ∣H^∣ψ⟩\langle \psi | \hat{H} | \psi \rangle​) by using statevector simulator.

Now we use a classical optimizer (like Coybala, Powell, etc.) to find the minimum ground state energy of the system by iterating the process multiple times. 

Manually, the expectation value of each Hamiltonian terms (in Pauli form) can be calculated by placing each term one by one in the circuit after ansatz circuit and measuring all the qubits. All the gates are to be placed at their respected qubits. But for obtaining the ground state energy, we have to run the circuit multiple times by changing the parameters, which we can do using QISKit and Python programming language. For further calculations, the Ref. can be followed. The initial parameters may be chosen at random. It can be adjusted according to the measured expectation value.

Quantum circuit for calculating average of XXXX term in the Hamiltonian

from qiskit import QuantumCircuit, execute
from strangeworks.qiskit import get_backend
from math import pi

#taking 4 quantum and 4 classical registers
qc = QuantumCircuit(4, 4)

#preparing the ansatz



#applying XXXX terms of the Hamiltonian

#measuring the qubits
qc.measure([0,1,2,3], [0,1,2,3])

# get a Strangeworks backend
backend = get_backend("qasm_simulator")

# execute circuit
execute(qc, backend, shots=1024)

In the circuit, u3 operators, are used where the parameters θ\theta​, ϕ\phi​ and λ\lambda can be taken in a range, and average of XXXX term can be calculated from the probabilities. Similarly, the total Hamiltonian can be calculated in a range of parameters and using the above mentioned optimized techniques, it can be used and thus ground state energy values can be evaluated.


To conclude, we have shown that the VQE is good for small ionic systems but as the complexity of the system increases the error margin increases. We have demonstrated that the UCCSD variational forms are much better than any other forms for chemistry related calculations. It is observed that the energy calculation is affected by limiting our system by considering Slater type basis set and by reducing the active space. In future, more complex molecules/ions can be simulated on the Noisy Intermediate-Scale Quantum (NISQ) era quantum systems. Still, we need to develop more advanced quantum machines to simulate them.


Utkarsh Singh, Shubham Kumar, Bikash K. Behera, and Prasanta K. Panigrahi, Efficient Simulation of Charged Systems Using Variational Quantum Eigensolver, [DOI: 10.13140/RG.2.2.29794.48328]