<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <img src="../../../images/qiskit-heading.gif" alt="Note
Step1:  Step 1
Step2:  Step 2
Step3:  We use the classical eigen decomposition to get the smallest eigenvalue as a reference.
Step4:  Step 3
Step5:  Step 4
<ASSISTANT_TASK:>
Python Code:
# 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
Explanation: <img src="../../../images/qiskit-heading.gif" alt="Note: In order for images to show up in this jupyter notebook you need to select File => Trusted Notebook" width="500 px" align="left">
Qiskit Aqua: Chemistry advanced how to
The latest version of this notebook is available on https://github.com/Qiskit/qiskit-tutorial.
Contributors
Richard Chen<sup>[1]</sup>, Antonio Mezzacapo<sup>[1]</sup>, Marco Pistoia<sup>[1]</sup>, Stephen Wood<sup>[1]</sup>
Affiliation
<sup>[1]</sup>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:
  1. define a molecule and get integrals from a driver (PySCF)
  2. define how to map Fermionic Hamiltonian into qubit Hamiltonian
  3. initiate and config dynamically-loaded instances, such as algorithm, optimizer, variational_form, and initial_state
  4. run the algorithm and retrieve the results
End of explanation
# 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)
Explanation: 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.
End of explanation
# 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)
Explanation: 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.
End of explanation
# 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))
Explanation: We use the classical eigen decomposition to get the smallest eigenvalue as a reference.
End of explanation
from qiskit import IBMQ
IBMQ.load_accounts()
# 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)
Explanation: 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.
End of explanation
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']))
Explanation: Step 4: Run algorithm and retrieve the results
The smallest eigenvalue is stored in the first entry of the eigvals key.
End of explanation
<END_TASK>