3
$\begingroup$

I'm new to QISKIT Aqua and I'm trying to follow the Ground State tutorial to solve the Water molecule ground state, but when I execute the calculations the kernel stops. I'm using IBM Q Experience jupyter notebooks and the only thing I changed from the original tutorial is the molecule. I Changed it by excluding the previous molecule and adding this to the code:

import numpy as np

angle = np.deg2rad(104.45/2)
dist = 0.9584
molecule = Molecule(geometry=[['O', [0., 0., 0.]],
                              ['H', [dist*np.sin(angle), -dist*np.cos(angle), 0.]],
                              ['H', [-dist*np.sin(angle), -dist*np.sin(angle), 0.]]],
                     charge=0, multiplicity=1)
$\endgroup$
4
  • $\begingroup$ What is the error you encountered? If you are using Jupyter notebook in the IBM Quantum Experience then it could means that you are exceeding the allotted memory... but I don't think that is the case for water unless you are using a pretty sophisticated basis set. $\endgroup$ – KAJ226 Dec 21 '20 at 21:07
  • $\begingroup$ While running I noticed that the memory peaks at 6.9GB just before stopping. I think it might be a memory problem so I'll try to run on my own PC. Other doubt I have is the "charge" value for the water molecule, as it's a polar molecule. Is 0 right? What value should I input? $\endgroup$ – Nillmer Dec 21 '20 at 21:16
  • $\begingroup$ That should be fine I think. I put a working code in the answer. I used the Z-matrix representation for ease of use but XYZ format should also work. One possible issue you might encounter is the"failed electronic structure" error from PySCF. I see that this sometime occur... but it would work if you rerun that same point separately... $\endgroup$ – KAJ226 Dec 21 '20 at 21:38
  • 1
    $\begingroup$ It worked perfectly, thank you! $\endgroup$ – Nillmer Dec 21 '20 at 22:46
3
$\begingroup$

Here is the working code. Because of the symmetry in water molecule, I used the Z Matrix representation instead of the usual XYZ format. You can and should be able to execute this within the IBM Quantum Experience.

import os
os.environ['QISKIT_SETTINGS']=''
import matplotlib.pyplot as plt
%matplotlib inline
%config InlineBackend.figure_format = 'svg' # Makes the images look nice
import numpy as np
from qiskit.aqua import QuantumInstance
from qiskit.aqua.algorithms import NumPyMinimumEigensolver, VQE, NumPyEigensolver
from qiskit.aqua.components.optimizers import COBYLA
from qiskit.aqua.components.initial_states import Custom
from qiskit import IBMQ, BasicAer, Aer
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
from qiskit.circuit.library import TwoLocal
import warnings
warnings.filterwarnings("ignore",category=DeprecationWarning)

provider = IBMQ.load_account()
simulator_backend = Aer.get_backend("qasm_simulator")
hardware_backend = provider.get_backend('ibmq_valencia')
quantum_instance = QuantumInstance(backend = simulator_backend, shots= 1, optimization_level= 3)

atom_structure = "H; O 1 0.958; H 2 0.958 1 104.47"  
orbital_reduction=[ ] 
[unit, charge, spin, basis]= [UnitsType.ANGSTROM, 0, 0, 'sto6g'] 
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 = 'auto')
qubitOp, aux_ops = ferOp.run(qmolecule)

energy_shift = ferOp._energy_shift
nuclear_repulsion_energy = ferOp._nuclear_repulsion_energy 
total_shift = energy_shift + nuclear_repulsion_energy 

print('HF Energy is:' , ferOp._hf_energy )
print('\n *********************************Exact Eigensolver*************************************** \n')
exact_result = NumPyMinimumEigensolver(qubitOp).run() 
ee_result = ferOp.process_algorithm_result(exact_result)
print(ee_result)

    
#----------------------- VQE PROCESSED ---------------------------# 
var_form= TwoLocal(qubitOp.num_qubits, ['ry'], ['cx'], entanglement= 'linear',
                    reps= 2,skip_unentangled_qubits=True, skip_final_rotation_layer = False, insert_barriers=True,
                    initial_state = None) 


optimizer = COBYLA(maxiter= 50,tol=0.000000001)   
vqe = VQE(qubitOp, var_form, optimizer, initial_point = None, 
          include_custom = True, aux_operators=aux_ops)
vqe_result = vqe.run(quantum_instance)
vqe_processed_result = ferOp.process_algorithm_result(vqe_result)
vqe_energy = vqe_result['eigenvalue'] + total_shift
total_ang_moment = vqe_processed_result['total_angular_momentum'] 
print('Optimal VQE energy at this step is:', vqe_energy)
print('The corresponding Total Angular Momentum is', total_ang_moment)
vqe.get_optimal_circuit().draw('mpl', style = {'name': 'bw'}, scale = 0.75, plot_barriers = False) 

Ran this gave me the output:

HF Energy is: -75.6787913259016

 *********************************Exact Eigensolver*************************************** 

=== GROUND STATE ENERGY ===
 
* Electronic ground state energy (Hartree): -84.916233233599
  - computed part:      -23.633272241623
  - frozen energy part: -61.282960991976
  - particle hole part: 0.0
~ Nuclear repulsion energy (Hartree): 9.187404411009
> Total ground state energy (Hartree): -75.72882882259
Optimal VQE energy at this step is: (-72.69078459811197+0j)
The corresponding Total Angular Momentum is 1.2057181036792375

and the optimal circuit: enter image description here

As you can see this var_form is too long to be executed with good result for the hardware so maybe you can modify it. Changing the entangle_map might help. So instead of using entanglement= 'linear' you can set it to be something else. Note that you can also specify how the CNOT gates are connected. That is if you specify:

var_form= TwoLocal(qubitOp.num_qubits, ['ry'], ['cx'], entanglement= [ [0,1],[2,3],[4,5],[6,7] ],
                    reps= 1,skip_unentangled_qubits=True, skip_final_rotation_layer = False, insert_barriers=True,
                    initial_state = None) 

Then you would get something like:

     ┌──────────┐ ░       ░  ┌──────────┐
q_0: ┤ RY(θ[0]) ├─░───■───░──┤ RY(θ[8]) ├
     ├──────────┤ ░ ┌─┴─┐ ░  ├──────────┤
q_1: ┤ RY(θ[1]) ├─░─┤ X ├─░──┤ RY(θ[9]) ├
     ├──────────┤ ░ └───┘ ░ ┌┴──────────┤
q_2: ┤ RY(θ[2]) ├─░───■───░─┤ RY(θ[10]) ├
     ├──────────┤ ░ ┌─┴─┐ ░ ├───────────┤
q_3: ┤ RY(θ[3]) ├─░─┤ X ├─░─┤ RY(θ[11]) ├
     ├──────────┤ ░ └───┘ ░ ├───────────┤
q_4: ┤ RY(θ[4]) ├─░───■───░─┤ RY(θ[12]) ├
     ├──────────┤ ░ ┌─┴─┐ ░ ├───────────┤
q_5: ┤ RY(θ[5]) ├─░─┤ X ├─░─┤ RY(θ[13]) ├
     ├──────────┤ ░ └───┘ ░ ├───────────┤
q_6: ┤ RY(θ[6]) ├─░───■───░─┤ RY(θ[14]) ├
     ├──────────┤ ░ ┌─┴─┐ ░ ├───────────┤
q_7: ┤ RY(θ[7]) ├─░─┤ X ├─░─┤ RY(θ[15]) ├
     └──────────┘ ░ └───┘ ░ └───────────┘
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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