Why is HF Initial point seems to be better than MP2 initial point in qiskit?

Question

I am trying to design a custom VQE algorithm using qiskit. The main customization of the algorithm is in the minimization of expectation value of the Hamiltonian. To accomplish this, I have used the following code. I have chosen the LiH molecule in sto3g basis.

First, I perform SCF calculation using pyscf just to get the number of orbitals, orbital occupation etc.

Code:

from pyscf import gto, scf, fci
import numpy as np

mol = gto.Mole()
mol.build(atom='Li 0 0 0; H 0 0 1', basis='sto3g', symmetry=1)

mf = scf.RHF(mol)
mf.kernel()
print ('Occupation in MOs: ', mf.mo_occ)

Output:

converged SCF energy = -7.76736213574856
Occupation in MOs:  [2. 2. 0. 0. 0. 0.]

Next, I generate the HF reference state and the UCCSD operator.

from qiskit_nature.circuit.library.initial_states.hartree_fock import *
from qiskit_nature.circuit.library.ansatzes import UCC

#Writing the HF state in qubit occupation form 
#The following code gives a circuit to prepare the HF state from the state with all qubits =0
HF = HartreeFock(num_spin_orbitals=len(mf.mo_occ)*2, num_particles=mol.nelec, qubit_converter=QubitConverter(mapper=JordanWignerMapper(), two_qubit_reduction=True))
#print (HF)

UCC_operator = UCC(num_spin_orbitals=len(mf.mo_occ)*2, num_particles=mol.nelec, excitations='sd', qubit_converter=QubitConverter(mapper=JordanWignerMapper(), two_qubit_reduction=True), initial_state = HF)

Then I write down the molecule whose energy needs to be calculated and obtain its Hamiltonian.

from qiskit_nature.drivers import UnitsType, Molecule
from qiskit_nature.drivers.second_quantization import ElectronicStructureDriverType, ElectronicStructureMoleculeDriver
from qiskit_nature.problems.second_quantization import ElectronicStructureProblem

molecule = Molecule(
    # coordinates are given in Angstrom
    geometry=[
        ["Li", [0.0, 0.0, 0.0]],
        ["H", [0.0, 0.0, 1.0]]
    ],
    multiplicity=1,  # = 2*spin + 1
    charge=0,)

driver = ElectronicStructureMoleculeDriver(molecule, basis="sto3g", 
                    driver_type=ElectronicStructureDriverType.PYSCF)

problem = ElectronicStructureProblem(driver)
second_q_op = problem.second_q_ops()
qubit_converter = QubitConverter(mapper = JordanWignerMapper(), two_qubit_reduction = True)
Hamiltonian = qubit_converter.convert(second_q_op[0], num_particles=problem.num_particles)

Now since I am customizing in the minimization part of the VQE algorithm, I will just prepare the VQE instance and then use the get_energy_evaluation method that qiskit provides to use UCCSD-VQE circuit as a parametrized circuit to be evaluated at the desired parameter values.

from qiskit import Aer
from qiskit.algorithms.optimizers import SLSQP
from qiskit.algorithms import VQE

backend = Aer.get_backend('aer_simulator_statevector')
optimizer = SLSQP(maxiter=100)

algorithm = VQE(UCC_operator, optimizer=optimizer, quantum_instance=backend)

energy_eval_func = algorithm.get_energy_evaluation(Hamiltonian, False)

Now that I have obtained the function to evaluate expectation value of the Hamiltonian of the molecule wrt to a desired UCC state (i.e. T1 and T2 amplitudes), I choose 2 kinds of parameters for the UCC operator and observe the results.

1. All the parameters are equal to 0. This is equivalent to the HFInitialpoint. I use this as follows.

t = UCC_operator.parameters
num_params = len(t)
energy_eval_func(np.zeros(len(t))
print (energy_eval_func)

Output:

-9.354893768508214

2. I choose the parameter elements to be MP2 amplitudes. This can be achieved by using MP2InitialPoint as follows.

from qiskit_nature.algorithms.initial_points import *

mp2_initial_point = MP2InitialPoint()
mp2_initial_point.compute(UCC_operator, driver_results)
initial_point = mp2_initial_point.to_numpy_array()

energy_eval_func(initial_point)
print (energy_eval_func)

Output:

-9.324343952128133

This is what I do not understand. Given that the theoretical energy is -9.371991912791, why is the energy lower for HF initial point than the MP2 initial point? I was of the opinion that it should have been the other way around. I do not see where have I gone wrong. It would be really helpful if you could clarify this doubt. Thank you.

Answer 1

The reason that you are not observing the expected MP2 energy after that first energy evaluation is due to the fact that you are combining it with the unitary coupled cluster ansatz. Compared to the classical coupled cluster ansatz, the difference is that you are also including the complex conjugate of the double excitations. This is necessary in quantum computing because we need to perform unitary operations.

However, that being said, while looking into this in a bit more detail, I did find a discrepancy of 0.1mHa in the computed MP2 energy correction. I have not yet been able to debug this 100% but will soon open an issue on the Qiskit Nature repository about this. I did not want to leave this question unanswered that much longer though. I will add a comment here, with a link to that issue (and hopefully corresponding fix) once published.

EDIT: I am finally adding the link to the issue I opened on Qiskit Nature: https://github.com/Qiskit/qiskit-nature/issues/739