This is a Python program written to calculate the ground state energy of the helium hydride ion with focus on classically simulating & investigating the variational eigensolver of Peruzzo, McClean et al. arXiv:1304.3061v1 (October 2017). There are two purposes for this program. One is to be able to quickly test ideas that may improve the efficiency of the variational eigensolver (esp. with view towards hybrid classical-quantum methods - e.g. optimisation under noise - or merging ideas from the phase estimation algorithm). Two is to plot out the quantum variational surface and the classically reachable points on it (e.g. stabiliser states) to better understand the quantum advantage.

Technical details:

The program is based upon 2 STO spatial minimal basis functions as in Szabo & Ostlund (S.O.) which is then used to carry out a self-consistent field calculation to solve the Fock eigen-equation. As in S.O. the 2 STO basis functions are each approximated by contracting 3 Gaussian functions (hence basis name is in fact STO-3G), this can be regarded as a technical trick which makes it much easier to evaluate integrals involving the STO functions by giving these integrals analytical expressions. The code up to this point (i.e. pre_scf.py & scf.py) is simply a Python translation of the FORTRAN found in the Appendix B of S.O.

In post_scf.py, I recalculate the (electronic) Hamiltonian matrix elements between the 2 Fock eigen-spatial-orbitals which give the coefficients in the second quantisation formulation of the Hamiltonian; post_scf.py then goes on to express this Hamiltonian in the 4 - qubit basis under the J-W transformation (ai† -> Z.Z...Z.sigma-.I.I...I, ai -> Z.Z...Z.sigma+.I.I...I) which then gives a 16-by-16 matrix via the standard matrix representation of Pauli operators. Note that this matrix representation corresponds to the state-qubit transformation |X1 X2> -> |1 1 0 0>, |X1 X4> -> |1 0 0 1> etc. which then gives a 16-by-1 vector via the canonical ordering of binary strings.

In ve.py, I implement the unitary coupled cluster trial state exp(T-T†)|HF> where T = T1 + T2 in the 2 electron case and T1 = sum(t_ij.ai†aj†), T2 = sum(t_ijkl.ai†aj†akal), the sum in each case is finite because of the introduction of a finite # of basis functions). The energy associated with this trial state can be computed in two ways by ve.energy depending on the mode. In mode = 0, the energy is calculated directly via matrix multiplication. In mode = 1, it samples from the output distribution of quantum circuits whose output expectation is a summand of this energy, the number of samples taken can be varied which affects the noisiness of this summand. The overall energy is then calculated classically by adding up the summands. With the energy function defined, ve.py finally uses the Nelder-Mead method to minimise this energy with respect to trial state parameters t.

For example outputs, try running ve.optimise(), ve.plot_2d(), ve.plot_3d() or ve.plot_disa().