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) mf.kernel() 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 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, 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.
- 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)
- I choose the parameter elements to be MP2 amplitudes. This can be achieved by using
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)
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.