# Ground state compute issue by qiskit

I use the following code to calculate the ground state for the LiH molecule in an active space. I come across two problems. The first is I found that the Hartree Fock state gave energy that is far from the ground state energy. But the qc state is the true Hartree Fock state. The second problem is using the ground_state, it just outputs 0. Can anybody help me with that?

from qiskit import *
import numpy as np

#Circuit imports
from qiskit_nature.drivers import UnitsType
from qiskit_nature.drivers.second_quantization import PySCFDriver
from qiskit_nature.problems.second_quantization.electronic import ElectronicStructureProblem
from qiskit_nature.circuit.library import HartreeFock
from qiskit_nature.transformers.second_quantization.electronic import FreezeCoreTransformer, ActiveSpaceTransformer
from qiskit_nature.algorithms import GroundStateEigensolver
from qiskit_nature.results import EigenstateResult
from qiskit import Aer
from qiskit_nature.mappers.second_quantization import ParityMapper, JordanWignerMapper
from qiskit_nature.converters.second_quantization import QubitConverter
from qiskit.algorithms import VQE, NumPyMinimumEigensolver
from qiskit.opflow.state_fns import StateFn, CircuitStateFn
import matplotlib
from qiskit.tools.visualization import circuit_drawer
from qiskit.quantum_info import state_fidelity, Statevector
matplotlib.use('Agg')
np.set_printoptions(threshold=np.inf)

driver = PySCFDriver(atom='H .0, .0, .0; Li .0, .0, 0.5',
unit=UnitsType.ANGSTROM,
basis='sto3g')

# set parameters
molecule = "LiH"
set_as = True # set the active space or not
as_particle = 2 # number of particles in the active space
as_mol_orbital = 3 # number of molecular orbitals in the active space
at = ActiveSpaceTransformer(as_particle, as_mol_orbital, active_orbitals=[1,2,5])
ft = FreezeCoreTransformer()
as_problem = ElectronicStructureProblem(driver, transformers=[at])

# generate the second-quantized operators for active space
as_second_q_ops = as_problem.second_q_ops()
as_main_op = as_second_q_ops[0]

as_particle_number = as_problem.grouped_property_transformed.get_property("ParticleNumber")
as_num_particles = (as_particle_number.num_alpha, as_particle_number.num_beta)
as_num_spin_orbitals = as_particle_number.num_spin_orbitals

as_qubit_op = converter.convert(main_op, num_particles=num_particles)
as_init_state = HartreeFock(as_num_spin_orbitals, as_num_particles, converter)
qc = QuantumCircuit(as_num_spin_orbitals)
qc.x(1)
qc.x(3)

# calculate the energy
numpy_solver = NumPyMinimumEigensolver()
calc = GroundStateEigensolver(converter, numpy_solver)
calc.returns_groundstate()
as_res_ref = calc.solve(as_problem)
ground_state = as_res_ref.eigenstates[0]
print('G.S. energy is', as_res_ref.eigenenergies.min())
gse = (~ground_state @ as_qubit_op @ ground_state).eval()
print('Second G.S. energy is',gse)

e = ~StateFn(as_qubit_op) @ CircuitStateFn(primitive=as_init_state, coeff=1.)
e = e.eval()
print('HF state energy is',e)
e2 = ~StateFn(as_qubit_op) @ CircuitStateFn(primitive=qc, coeff=1.)
e2 = e2.eval()
print('qc state energy is',e2)
$$`$$
• Your code snippet is incomplete. You don't have the converter specified. Working with the assumption that you are using the JordanWignerMapper, then qc is not the correct HF groundstate. What makes you believe that it should be? Nov 9 at 7:22
• In this case, the molecule is weakly correlated, which means the HF state has a fairly large overlap with the ground state. For the HF state generate by qiskit. The energy is totally off. Nov 9 at 8:00
• Your reply still does not explain why you believe that the circuit qc.x(1); qc.x(3) should be the correct HF state compared to qc.x(0); qc.x(3) (which is what Qiskit Nature's HF state does), matching the expectation that you have one alpha and one beta electron in the lowest alpha and beta spin orbitals, respectively (assuming JW mapping once again). Nov 9 at 9:16
• For this image ,nature.com/articles/s41467-019-10988-2/figures/2, it can be seen that the energy difference of the HF state with ground state is around 0.02 when bond length is 1.5A. However, the program output a difference, which is 0.8. Nov 9 at 12:28