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.


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)
print ('Occupation in MOs: ', mf.mo_occ)


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
        ["Li", [0.0, 0.0, 0.0]],
        ["H", [0.0, 0.0, 1.0]]
    multiplicity=1,  # = 2*spin + 1

driver = ElectronicStructureMoleculeDriver(molecule, basis="sto3g", 

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)
print (energy_eval_func)


  1. 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()

print (energy_eval_func)



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.


1 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

  • $\begingroup$ Thank you for the answer. So are you saying in the first paragraph that the results I have obtained are consistent with the theory? In the sense MP2 initial point is not giving a better result because "It need not" as I am using UCC and not CC. Now for the second paragraph, are you saying, the exact value of the MP2 energy I obtained is incorrect by about 0.1mH and its due to the qiskit code itself and has nothing to do with UCC vs CC explanation you gave above? $\endgroup$
    – Kurious
    Jul 14, 2022 at 10:00
  • $\begingroup$ For the first paragraph: indeed, it does not need to match the MP2 energy because the ansatz is UCC rather than CC. In that sense everything is working as expected. For the second paragraph: This is not for the value that you printed above but the value reported by mp2_initial_point.get_energy(). Comparing this value with the one computed by PySCF shows a discrepancy of 0.1mHa (see also the link I finally added to the answer above) $\endgroup$
    – mrossinek
    Jul 15, 2022 at 12:28
  • $\begingroup$ Excellent! Thank you. Actually what you have explained in the link, solved another of my doubt. I was getting incorrect values in MP2 initial point calculated from qiskit at few places when compared with PySCF. They were precisely corresponding to the excitations you have mentioned in the link and 2 more corresponding to the identical spin excitations (alpha, alpha goes to alpha, alpha and same for beta) which can be derived from the former amplitudes using this identity. $\endgroup$
    – Kurious
    Jul 15, 2022 at 15:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.