Note

This page was generated from tutorials/chemistry/1_programmatic_approach.ipynb.

Qiskit Chemistry, Programmatic Approach

The latest version of this notebook is available on https://github.com/Qiskit/qiskit-tutorial.


Contributors

Richard Chen[1], Antonio Mezzacapo[1], Marco Pistoia[1], Stephen Wood[1] ### Affiliation - [1]IBMQ

Introduction

This notebook illustrates how to use Qiskit Chemistry’s programmatic APIs.

In this notebook, we decompose the computation of the ground state energy of a molecule into 4 steps: 1. Define a molecule and get integrals from a computational chemistry driver (PySCF in this case) 2. Construct a Fermionic Hamiltonian and map it onto a qubit Hamiltonian 3. Instantiate and initialize dynamically-loaded algorithmic components, such as the quantum algorithm VQE, the optimizer and variational form it will use, and the initial_state to initialize the variational form 4. Run the algorithm on a quantum backend and retrieve the results

[1]:
# import common packages
import numpy as np

from qiskit import Aer

# lib from Qiskit Aqua
from qiskit.aqua import QuantumInstance
from qiskit.aqua.algorithms import VQE, NumPyMinimumEigensolver
from qiskit.aqua.operators import Z2Symmetries
from qiskit.aqua.components.optimizers import COBYLA

# lib from Qiskit Aqua Chemistry
from qiskit.chemistry import FermionicOperator
from qiskit.chemistry.drivers import PySCFDriver, UnitsType
from qiskit.chemistry.components.variational_forms import UCCSD
from qiskit.chemistry.components.initial_states import HartreeFock

Step 1: Define a molecule

Here, we use LiH in the sto3g basis with the PySCF driver as an example. The molecule object records the information from the PySCF driver.

[2]:
# using driver to get fermionic Hamiltonian
# PySCF example
driver = PySCFDriver(atom='Li .0 .0 .0; H .0 .0 1.6', unit=UnitsType.ANGSTROM,
                     charge=0, spin=0, basis='sto3g')
molecule = driver.run()

Step 2: Prepare qubit Hamiltonian

Here, we setup the to-be-frozen and to-be-removed orbitals to reduce the problem size when we map to the qubit Hamiltonian. Furthermore, we define the mapping type for the qubit Hamiltonian. For the particular parity mapping, we can further reduce the problem size.

[3]:
# please be aware that the idx here with respective to original idx
freeze_list = [0]
remove_list = [-3, -2] # negative number denotes the reverse order
map_type = 'parity'

h1 = molecule.one_body_integrals
h2 = molecule.two_body_integrals
nuclear_repulsion_energy = molecule.nuclear_repulsion_energy

num_particles = molecule.num_alpha + molecule.num_beta
num_spin_orbitals = molecule.num_orbitals * 2
print("HF energy: {}".format(molecule.hf_energy - molecule.nuclear_repulsion_energy))
print("# of electrons: {}".format(num_particles))
print("# of spin orbitals: {}".format(num_spin_orbitals))
HF energy: -8.854072040283643
# of electrons: 4
# of spin orbitals: 12
[4]:
# prepare full idx of freeze_list and remove_list
# convert all negative idx to positive
remove_list = [x % molecule.num_orbitals for x in remove_list]
freeze_list = [x % molecule.num_orbitals for x in freeze_list]
# update the idx in remove_list of the idx after frozen, since the idx of orbitals are changed after freezing
remove_list = [x - len(freeze_list) for x in remove_list]
remove_list += [x + molecule.num_orbitals - len(freeze_list)  for x in remove_list]
freeze_list += [x + molecule.num_orbitals for x in freeze_list]

# prepare fermionic hamiltonian with orbital freezing and eliminating, and then map to qubit hamiltonian
# and if PARITY mapping is selected, reduction qubits
energy_shift = 0.0
qubit_reduction = True if map_type == 'parity' else False

ferOp = FermionicOperator(h1=h1, h2=h2)
if len(freeze_list) > 0:
    ferOp, energy_shift = ferOp.fermion_mode_freezing(freeze_list)
    num_spin_orbitals -= len(freeze_list)
    num_particles -= len(freeze_list)
if len(remove_list) > 0:
    ferOp = ferOp.fermion_mode_elimination(remove_list)
    num_spin_orbitals -= len(remove_list)

qubitOp = ferOp.mapping(map_type=map_type, threshold=0.00000001)
qubitOp = Z2Symmetries.two_qubit_reduction(qubitOp, num_particles) if qubit_reduction else qubitOp
qubitOp.chop(10**-10)

print(qubitOp.print_details())
print(qubitOp)
IIII    (-0.2076593350197073+0j)
IIIZ    (-0.09376337484626507+0j)
IIZX    (-0.0031775814548800733+0j)
IIIX    (0.0031775814548800733+0j)
IIXX    (-0.0012513965999528314+0j)
IIYY    (0.0012513965999528314+0j)
IIZZ    (-0.21162509515109656+0j)
IIXZ    (0.019200533863101478+0j)
IIXI    (0.019200533863101478+0j)
IIZI    (0.35810269945770484+0j)
IZII    (0.09376337484626497+0j)
ZXII    (0.0031775814548800686+0j)
IXII    (0.0031775814548800686+0j)
XXII    (-0.0012513965999528375+0j)
YYII    (0.0012513965999528375+0j)
ZZII    (-0.21162509515109656+0j)
XZII    (-0.019200533863101464+0j)
XIII    (0.019200533863101464+0j)
ZIII    (-0.3581026994577048+0j)
IZIZ    (-0.12182774215820326+0j)
IZZX    (0.01214489722808736+0j)
IZIX    (-0.01214489722808736+0j)
IZXX    (0.03169874598733563+0j)
IZYY    (-0.03169874598733563+0j)
IXIZ    (0.01214489722808736+0j)
ZXIZ    (0.01214489722808736+0j)
IXZX    (-0.0032659954996687823+0j)
ZXZX    (-0.0032659954996687823+0j)
IXIX    (0.0032659954996687823+0j)
ZXIX    (0.0032659954996687823+0j)
IXXX    (-0.008650156860622697+0j)
ZXXX    (-0.008650156860622697+0j)
IXYY    (0.008650156860622697+0j)
ZXYY    (0.008650156860622697+0j)
YYIZ    (0.03169874598733563+0j)
XXIZ    (-0.03169874598733563+0j)
YYZX    (-0.008650156860622697+0j)
XXZX    (0.008650156860622697+0j)
YYIX    (0.008650156860622697+0j)
XXIX    (-0.008650156860622697+0j)
YYXX    (-0.030981613344623248+0j)
XXXX    (0.030981613344623248+0j)
YYYY    (0.030981613344623248+0j)
XXYY    (-0.030981613344623248+0j)
ZZIZ    (0.05590251078516957+0j)
ZZZX    (0.001871042751423773+0j)
ZZIX    (-0.001871042751423773+0j)
ZZXX    (0.003104004116063206+0j)
ZZYY    (-0.003104004116063206+0j)
XIIZ    (0.01284172318076957+0j)
XZIZ    (-0.01284172318076957+0j)
XIZX    (-0.002352152173256113+0j)
XZZX    (0.002352152173256113+0j)
XIIX    (0.002352152173256113+0j)
XZIX    (-0.002352152173256113+0j)
XIXX    (-0.007975908750574185+0j)
XZXX    (0.007975908750574185+0j)
XIYY    (0.007975908750574185+0j)
XZYY    (-0.007975908750574185+0j)
ZIIZ    (0.11346110712684569+0j)
ZIZX    (-0.010838363828759677+0j)
ZIIX    (0.010838363828759677+0j)
ZIXX    (-0.03355135311123067+0j)
ZIYY    (0.03355135311123067+0j)
IZZZ    (-0.05590251078516957+0j)
IZXZ    (-0.01284172318076957+0j)
IZXI    (-0.01284172318076957+0j)
IXZZ    (-0.001871042751423773+0j)
ZXZZ    (-0.001871042751423773+0j)
IXXZ    (0.002352152173256113+0j)
ZXXZ    (0.002352152173256113+0j)
IXXI    (0.002352152173256113+0j)
ZXXI    (0.002352152173256113+0j)
YYZZ    (-0.003104004116063206+0j)
XXZZ    (0.003104004116063206+0j)
YYXZ    (0.007975908750574185+0j)
XXXZ    (-0.007975908750574185+0j)
YYXI    (0.007975908750574185+0j)
XXXI    (-0.007975908750574185+0j)
ZZZZ    (0.08447056807294112+0j)
ZZXZ    (-0.008994911953942263+0j)
ZZXI    (-0.008994911953942263+0j)
XIZZ    (-0.008994911953942263+0j)
XZZZ    (0.008994911953942263+0j)
XIXZ    (0.006612047066159558+0j)
XZXZ    (-0.006612047066159558+0j)
XIXI    (0.006612047066159558+0j)
XZXI    (-0.006612047066159558+0j)
ZIZZ    (0.060358912810790984+0j)
ZIXZ    (0.01101923164472494+0j)
ZIXI    (0.01101923164472494+0j)
IZZI    (0.11346110712684569+0j)
IXZI    (-0.010838363828759677+0j)
ZXZI    (-0.010838363828759677+0j)
YYZI    (-0.03355135311123067+0j)
XXZI    (0.03355135311123067+0j)
ZZZI    (-0.060358912810790984+0j)
XIZI    (-0.01101923164472494+0j)
XZZI    (0.01101923164472494+0j)
ZIZI    (-0.1134468030036645+0j)

Representation: paulis, qubits: 4, size: 100

We use the classical eigen decomposition to get the smallest eigenvalue as a reference.

[5]:
# Using exact eigensolver to get the smallest eigenvalue
exact_eigensolver = NumPyMinimumEigensolver(qubitOp)
ret = exact_eigensolver.run()
print('The computed energy is: {:.12f}'.format(ret.eigenvalue.real))
print('The total ground state energy is: {:.12f}'.format(ret.eigenvalue.real + energy_shift + nuclear_repulsion_energy))
The computed energy is: -1.077059745735
The total ground state energy is: -7.881072044031

Step 3: Initiate and configure dynamically-loaded instances

To run VQE with the UCCSD variational form, we require:

  • VQE algorithm

  • Classical Optimizer

  • UCCSD variational form

  • Prepare the initial state in the HartreeFock state

[Optional] Setup token to run the experiment on a real device

If you would like to run the experiment on a real device, you need to setup your account first.

Note: If you did not store your token yet, use IBMQ.save_account('MY_API_TOKEN') to store it first.

[6]:
# from qiskit import IBMQ
# provider = IBMQ.load_account()
[7]:
backend = Aer.get_backend('statevector_simulator')
[8]:
# setup COBYLA optimizer
max_eval = 200
cobyla = COBYLA(maxiter=max_eval)

# setup HartreeFock state
HF_state = HartreeFock(num_spin_orbitals, num_particles, map_type,
                       qubit_reduction)

# setup UCCSD variational form
var_form = UCCSD(num_orbitals=num_spin_orbitals, num_particles=num_particles,
                 active_occupied=[0], active_unoccupied=[0, 1],
                 initial_state=HF_state, qubit_mapping=map_type,
                 two_qubit_reduction=qubit_reduction, num_time_slices=1)

# setup VQE
vqe = VQE(qubitOp, var_form, cobyla)
quantum_instance = QuantumInstance(backend=backend)

Step 4: Run algorithm and retrieve the results

[9]:
results = vqe.run(quantum_instance)
print('The computed ground state energy is: {:.12f}'.format(results.eigenvalue.real))
print('The total ground state energy is: {:.12f}'.format(results.eigenvalue.real + energy_shift + nuclear_repulsion_energy))
print("Parameters: {}".format(results.optimal_point))
The computed ground state energy is: -1.057852461436
The total ground state energy is: -7.861864759732
Parameters: [3.03358125e-06 7.13391799e-05 3.03361004e-06 6.99450928e-05]
[10]:
import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

Version Information

Qiskit SoftwareVersion
Qiskit0.21.0
Terra0.15.2
Aer0.6.1
Ignis0.4.0
Aqua0.7.5
IBM Q Provider0.9.0
System information
Python3.8.5 (default, Jul 21 2020, 12:22:34) [GCC 7.5.0]
OSLinux
CPUs2
Memory (Gb)6.791389465332031
Wed Sep 16 14:26:20 2020 UTC

This code is a part of Qiskit

© Copyright IBM 2017, 2020.

This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.

Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.

[11]: