Qiskit Aqua: Chemistry advanced how to

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

In the aqua_chemistry_howto example, we show how to configure different parameters in an input dictionary for different experiments in Aqua Chemistry. Nonetheless, for the users who intersted in experimenting with their new algorithm or new components in the algorithm, which is outside of the bounds of the high-level APIs, the users may want to break down the computation into smaller steps.

In this notebook, we decompose the steps to compute the ground state energy of a molecule into 4 steps:

define a molecule and get integrals from a driver (PySCF)
define how to map Fermionic Hamiltonian into qubit Hamiltonian
initiate and config dynamically-loaded instances, such as algorithm, optimizer, variational_form, and initial_state
run the algorithm and retrieve the results



In [1]:

    
# import common packages
import numpy as np
from collections import OrderedDict

# lib from Qiskit Aqua Chemistry
from qiskit_aqua_chemistry import FermionicOperator

# lib from Qiskit Aqua
from qiskit_aqua import Operator
from qiskit_aqua import (get_algorithm_instance, get_optimizer_instance, 
                          get_variational_form_instance, get_initial_state_instance)

# lib for driver
from qiskit_aqua_chemistry.drivers import ConfigurationManager

Step 1: define a molecule

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



In [2]:

    
# using driver to get fermionic Hamiltonian
# PySCF example
cfg_mgr = ConfigurationManager()
pyscf_cfg = OrderedDict([('atom', 'Li .0 .0 .0; H .0 .0 1.6'), 
                         ('unit', 'Angstrom'), 
                         ('charge', 0), 
                         ('spin', 0), 
                         ('basis', 'sto3g')])
section = {}
section['properties'] = pyscf_cfg
driver = cfg_mgr.get_driver_instance('PYSCF')
molecule = driver.run(section)

Step 2: Prepare qubit Hamiltonian

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



In [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))

# 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 = qubitOp.two_qubit_reduced_operator(num_particles) if qubit_reduction else qubitOp
qubitOp.chop(10**-10)









    



HF energy: -8.854072040283643
# of electrons: 4
# of spin orbitals: 12

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



In [4]:

    
# Using exact eigensolver to get the smallest eigenvalue
exact_eigensolver = get_algorithm_instance('ExactEigensolver')
exact_eigensolver.init_args(qubitOp, k=1)
ret = exact_eigensolver.run()
print('The computed energy is: {:.12f}'.format(ret['eigvals'][0].real))
print('The total ground state energy is: {:.12f}'.format(ret['eigvals'][0].real + energy_shift + nuclear_repulsion_energy))









    



The computed energy is: -1.077059745735
The total ground state energy is: -7.881072044031

Step 3: Initiate and config dynamically-loaded instances

To run VQE with UCCSD variational form, we require

VQE algorithm
Classical Optimizer
UCCSD variational form
Prepare the initial state into HartreeFock state

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

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

Note: If you do not store your token yet, use IBMQ.save_accounts() to store it first.



In [5]:

    
from qiskit import IBMQ
IBMQ.load_accounts()



In [6]:

    
# setup COBYLA optimizer
max_eval = 200
cobyla = get_optimizer_instance('COBYLA')
cobyla.set_options(maxiter=max_eval)

# setup HartreeFock state
HF_state = get_initial_state_instance('HartreeFock')
HF_state.init_args(qubitOp.num_qubits, num_spin_orbitals, map_type, 
                   qubit_reduction, num_particles)

# setup UCCSD variational form
var_form = get_variational_form_instance('UCCSD')
var_form.init_args(qubitOp.num_qubits, depth=1, 
                   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_algorithm = get_algorithm_instance('VQE')
vqe_algorithm.setup_quantum_backend(backend='statevector_simulator')
vqe_algorithm.init_args(qubitOp, 'matrix', var_form, cobyla)

Step 4: Run algorithm and retrieve the results

The smallest eigenvalue is stored in the first entry of the eigvals key.



In [7]:

    
results = vqe_algorithm.run()
print('The computed ground state energy is: {:.12f}'.format(results['eigvals'][0]))
print('The total ground state energy is: {:.12f}'.format(results['eigvals'][0] + energy_shift + nuclear_repulsion_energy))
print("Parameters: {}".format(results['opt_params']))









    



The computed ground state energy is: -1.077059715882
The total ground state energy is: -7.881072014178
Parameters: [ 0.03406673  0.00540282  0.0345916   0.00525172 -0.03877438  0.06030781
  0.06050618 -0.11645761]