Simulation of a helium molecule using Qiskit

Question

I am trying to compute the ground state energy of He-He using VQE. For this purpose I have utilized Qiskit and have written the following script:

from qiskit.chemistry.drivers import PySCFDriver, UnitsType, Molecule
from qiskit.chemistry.transformations import FermionicTransformation, FermionicQubitMappingType
from qiskit.aqua.algorithms import NumPyMinimumEigensolver
from qiskit import BasicAer
from qiskit.aqua import QuantumInstance
from qiskit.chemistry.algorithms.ground_state_solvers.minimum_eigensolver_factories import VQEUCCSDFactory
from qiskit.aqua.algorithms import VQE
from qiskit.circuit.library import TwoLocal
from qiskit.chemistry.algorithms.ground_state_solvers import GroundStateEigensolver

molecule = Molecule(geometry=[['He', [0., 0., 0.]],
                              ['He', [0., 0., 3.1]]],
                     charge=0, multiplicity=1)

driver = PySCFDriver(molecule = molecule, unit=UnitsType.ANGSTROM, basis='sto3g')
transformation = FermionicTransformation(qubit_mapping=FermionicQubitMappingType.JORDAN_WIGNER)

numpy_solver = NumPyMinimumEigensolver()
vqe_solver = VQEUCCSDFactory(QuantumInstance(BasicAer.get_backend('statevector_simulator')))


num_qubits = 4
tl_circuit = TwoLocal(num_qubits, ['h', 'rx'], 'cz',
                      entanglement='full', reps=3, parameter_prefix = 'y')

tl_circuit.draw(output = 'mpl')

another_solver = VQE(var_form = tl_circuit,
                     quantum_instance = QuantumInstance(BasicAer.get_backend('statevector_simulator')))


calc = GroundStateEigensolver(transformation, vqe_solver)
res = calc.solve(driver)

print(transformation)
print(res)

When I run it I get a following error message:

Traceback (most recent call last):
  File "test_remove.py", line 37, in <module>
    res = calc.solve(driver)
  File "/Users/monica/miniconda3/lib/python3.7/site-packages/qiskit/chemistry/algorithms/ground_state_solvers/ground_state_eigensolver.py", line 94, in solve
    solver = self._solver.get_solver(self.transformation)
  File "/Users/monica/miniconda3/lib/python3.7/site-packages/qiskit/chemistry/algorithms/ground_state_solvers/minimum_eigensolver_factories/vqe_uccsd_factory.py", line 212, in get_solver
    same_spin_doubles=self._same_spin_doubles)
  File "/Users/monica/miniconda3/lib/python3.7/site-packages/qiskit/chemistry/components/variational_forms/uccsd.py", line 153, in __init__
    excitation_type=self._excitation_type,)
  File "/Users/monica/miniconda3/lib/python3.7/site-packages/qiskit/chemistry/components/variational_forms/uccsd.py", line 509, in compute_excitation_lists
    raise ValueError('No unoccupied orbitals')
ValueError: No unoccupied orbitals

Can anyone explain to me what is going on? What am I doing wrong?

Answer 1

Try this:

from qiskit import IBMQ
from qiskit.aqua import QuantumInstance
from qiskit.aqua.components.optimizers import COBYLA
from qiskit.aqua.algorithms import VQE
from qiskit.circuit.library import TwoLocal
from qiskit.chemistry.qmolecule import QMolecule
from qiskit.chemistry.drivers import PySCFDriver, UnitsType
from qiskit.chemistry import FermionicOperator
from qiskit.chemistry.core import Hamiltonian, TransformationType, QubitMappingType

import warnings
warnings.filterwarnings("ignore",category=DeprecationWarning)


##### Defining Molecular Structure #######
atom_structure = "He 0. 0. 0.; He 3.1 0. 0"  
[unit, charge, spin, basis]= [UnitsType.ANGSTROM, 0, 0, 'sto3g'] 
driver = PySCFDriver(atom= atom_structure, unit=unit,charge=charge, spin=spin, basis=basis)
qmolecule = driver.run()
ferOp = Hamiltonian(transformation = TransformationType.FULL , 
                    qubit_mapping = QubitMappingType.PARITY,
                    two_qubit_reduction = True,
                    freeze_core = True, 
                    orbital_reduction = [ ], 
                    z2symmetry_reduction = None)
qubitOp, aux_ops = ferOp.run(qmolecule)
print( qubitOp.print_details() )

##### Exact Diagonalization Result ########
from qiskit.aqua.algorithms import NumPyMinimumEigensolver
exact_result = NumPyMinimumEigensolver(qubitOp).run() 
print('THis is the exact result:', exact_result)


##### VQE Result #######
tl_circuit = TwoLocal(qubitOp.num_qubits, ['h', 'rx'], 'cz',
                      entanglement='full', reps=3, parameter_prefix = 'y')
provider = IBMQ.load_account()
backend = provider.get_backend('ibmq_qasm_simulator')
quantum_instance = QuantumInstance(backend,
                                   shots = 8192,
                                   initial_layout = None,
                                   optimization_level = 3)  

optimizer = COBYLA(maxiter= 100,tol=0.000000001)   
vqe = VQE(qubitOp, tl_circuit, optimizer,initial_point= None)
vqe_result = vqe.run(quantum_instance)
vqe_result

Note that by specifying freeze_core = True and z2symmetry_reduction = 'auto' will help with reducing the number of qubits needed to solve your problem.