Note

This page was generated from tutorials/chemistry/05_Sampling_potential_energy_surfaces.ipynb.

Run interactively in the IBM Quantum lab.

Sampling the potential energy surface

Introduction

This interactive notebook demonstrates how to utilize the Potential Energy Surface (PES) samplers algorithm of qiskit chemistry to generate the dissociation profile of a molecule. We use the Born-Oppenhemier Potential Energy Surface (BOPES)and demonstrate how to exploit bootstrapping and extrapolation to reduce the total number of function evaluations in computing the PES using the Variational Quantum Eigensolver (VQE).

# import common packages
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from functools import partial

# qiskit
from qiskit.aqua import QuantumInstance
from qiskit import BasicAer
from qiskit.aqua.algorithms import NumPyMinimumEigensolver, VQE, IQPE
from qiskit.aqua.components.optimizers import SLSQP
from qiskit.circuit.library import ExcitationPreserving
from qiskit import BasicAer
from qiskit.aqua.algorithms import NumPyMinimumEigensolver, VQE
from qiskit.aqua.components.optimizers import SLSQP

# chemistry components
from qiskit.chemistry.components.initial_states import HartreeFock
from qiskit.chemistry.components.variational_forms import UCCSD
from qiskit.chemistry.drivers import PySCFDriver, UnitsType
from qiskit.chemistry.core import Hamiltonian, TransformationType, QubitMappingType
from qiskit.aqua.algorithms import VQAlgorithm, VQE, MinimumEigensolver
from qiskit.chemistry.transformations import FermionicTransformation
from qiskit.chemistry.drivers import PySCFDriver
from qiskit.chemistry.algorithms.ground_state_solvers import GroundStateEigensolver
from qiskit.chemistry.algorithms.pes_samplers.bopes_sampler import BOPESSampler
from qiskit.chemistry.drivers.molecule import Molecule
from qiskit.chemistry.algorithms.pes_samplers.extrapolator import *

import warnings
warnings.simplefilter('ignore', np.RankWarning)

/home/computertreker/git/qiskit/qiskit-dup/.tox/docs/lib/python3.7/site-packages/qiskit/chemistry/__init__.py:170: DeprecationWarning: The package qiskit.chemistry is deprecated. It was moved/refactored to qiskit_nature (pip install qiskit-nature). For more information see <https://github.com/Qiskit/qiskit-aqua/blob/master/README.md#migration-guide>
  warn_package('chemistry', 'qiskit_nature', 'qiskit-nature')

Here, we use the H2 molecule as a model system for testing.

ft = FermionicTransformation()
driver = PySCFDriver()
solver = VQE(quantum_instance=
             QuantumInstance(backend=BasicAer.get_backend('statevector_simulator')))
me_gsc = GroundStateEigensolver(ft, solver)

/home/computertreker/git/qiskit/qiskit-dup/.tox/docs/lib/python3.7/site-packages/qiskit/aqua/quantum_instance.py:137: DeprecationWarning: The class qiskit.aqua.QuantumInstance is deprecated. It was moved/refactored to qiskit.utils.QuantumInstance (pip install qiskit-terra). For more information see <https://github.com/Qiskit/qiskit-aqua/blob/master/README.md#migration-guide>
  'qiskit-terra')
/home/computertreker/git/qiskit/qiskit-dup/.tox/docs/lib/python3.7/site-packages/qiskit/aqua/components/optimizers/optimizer.py:50: DeprecationWarning: The package qiskit.aqua.components.optimizers is deprecated. It was moved/refactored to qiskit.algorithms.optimizers (pip install qiskit-terra). For more information see <https://github.com/Qiskit/qiskit-aqua/blob/master/README.md#migration-guide>
  'qiskit.algorithms.optimizers', 'qiskit-terra')
/home/computertreker/git/qiskit/qiskit-dup/.tox/docs/lib/python3.7/site-packages/qiskit/aqua/algorithms/vq_algorithm.py:72: DeprecationWarning: The class qiskit.aqua.algorithms.VQAlgorithm is deprecated. It was moved/refactored to qiskit.algorithms.VariationalAlgorithm (pip install qiskit-terra). For more information see <https://github.com/Qiskit/qiskit-aqua/blob/master/README.md#migration-guide>
  'qiskit-terra')

stretch1 = partial(Molecule.absolute_stretching, atom_pair=(1, 0))
mol = Molecule(geometry=[('H', [0., 0., 0.]),
                        ('H', [0., 0., 0.3])],
                       degrees_of_freedom=[stretch1],
                       )

# pass molecule to PSYCF driver
driver = PySCFDriver(molecule=mol)

print(mol.geometry)

[('H', [0.0, 0.0, 0.0]), ('H', [0.0, 0.0, 0.3])]

Make a perturbation to the molecule along the absolute_stretching dof

mol.perturbations = [0.2]
print(mol.geometry)

mol.perturbations = [0.6]
print(mol.geometry)

[('H', [0.0, 0.0, 0.0]), ('H', [0.0, 0.0, 0.5])]
[('H', [0.0, 0.0, 0.0]), ('H', [0.0, 0.0, 0.8999999999999999])]

Calculate bond dissociation profile using BOPES Sampler

Here, we pass the molecular information and the VQE to a built-in type called the BOPES Sampler. The BOPES Sampler allows the computation of the potential energy surface for a specified set of degrees of freedom/points of interest.

First we compare no bootstrapping vs bootstrapping

Bootstrapping the BOPES sampler involves utilizing the optimal variational parameters for a given degree of freedom, say r (ex. interatomic distance) as the initial point for VQE at a later degree of freedom, say r + \(\epsilon\). By default, if bootstrapping is set to True, all previous optimal parameters are used as initial points for the next runs.

distance1 = partial(Molecule.absolute_distance, atom_pair=(1, 0))
mol = Molecule(geometry=[('H', [0., 0., 0.]),
                        ('H', [0., 0., 0.3])],
                       degrees_of_freedom=[distance1],
                       )

# pass molecule to PSYCF driver
driver = PySCFDriver(molecule=mol)

# Specify degree of freedom (points of interest)
points = np.linspace(0.25, 2, 30)
results_full = {} # full dictionary of results for each condition
results = {} # dictionary of (point,energy) results for each condition
conditions = {False: 'no bootstrapping', True: 'bootstrapping'}


for value, bootstrap in conditions.items():
    # define instance to sampler
    bs = BOPESSampler(
        gss=me_gsc
        ,bootstrap=value
        ,num_bootstrap=None
        ,extrapolator=None)
    # execute
    res = bs.sample(driver,points)
    results_full[f'{bootstrap}'] =  res.raw_results
    results[f'points_{bootstrap}'] = res.points
    results[f'energies_{bootstrap}'] = res.energies

/home/computertreker/git/qiskit/qiskit-dup/.tox/docs/lib/python3.7/site-packages/qiskit/chemistry/fermionic_operator.py:386: DeprecationWarning: The package qiskit.aqua.operators is deprecated. It was moved/refactored to qiskit.opflow (pip install qiskit-terra). For more information see <https://github.com/Qiskit/qiskit-aqua/blob/master/README.md#migration-guide>
  pauli_list = WeightedPauliOperator(paulis=[])
/home/computertreker/git/qiskit/qiskit-dup/.tox/docs/lib/python3.7/site-packages/qiskit/chemistry/fermionic_operator.py:394: DeprecationWarning: The variable qiskit.aqua.aqua_globals is deprecated. It was moved/refactored to qiskit.utils.algorithm_globals (pip install qiskit-terra). For more information see <https://github.com/Qiskit/qiskit-aqua/blob/master/README.md#migration-guide>
  task_args=(threshold,), num_processes=aqua_globals.num_processes)

Compare to classical eigensolver

# define numpy solver
solver_numpy = NumPyMinimumEigensolver()
me_gsc_numpy = GroundStateEigensolver(ft, solver_numpy)
bs_classical = BOPESSampler(
               gss=me_gsc_numpy
               ,bootstrap=False
               ,num_bootstrap=None
               ,extrapolator=None)
# execute
res_np = bs_classical.sample(driver, points)
results_full['np'] =  res_np.raw_results
results['points_np'] = res_np.points
results['energies_np'] = res_np.energies

/home/computertreker/git/qiskit/qiskit-dup/.tox/docs/lib/python3.7/site-packages/qiskit/aqua/algorithms/minimum_eigen_solvers/minimum_eigen_solver.py:38: DeprecationWarning: The package qiskit.aqua.algorithms.minimum_eigen_solvers is deprecated. It was moved/refactored to qiskit.algorithms.minimum_eigen_solvers (pip install qiskit-terra). For more information see <https://github.com/Qiskit/qiskit-aqua/blob/master/README.md#migration-guide>
  'qiskit-terra')
/home/computertreker/git/qiskit/qiskit-dup/.tox/docs/lib/python3.7/site-packages/qiskit/aqua/algorithms/eigen_solvers/eigen_solver.py:38: DeprecationWarning: The package qiskit.aqua.algorithms.eigen_solvers is deprecated. It was moved/refactored to qiskit.algorithms.eigen_solvers (pip install qiskit-terra). For more information see <https://github.com/Qiskit/qiskit-aqua/blob/master/README.md#migration-guide>
  'qiskit-terra')

Plot results

fig = plt.figure()
for value, bootstrap in conditions.items():
    plt.plot(results[f'points_{bootstrap}'], results[f'energies_{bootstrap}'], label = f'{bootstrap}')
plt.plot(results['points_np'], results['energies_np'], label = 'numpy')
plt.legend()
plt.title('Dissociation profile')
plt.xlabel('Interatomic distance')
plt.ylabel('Energy')

Text(0, 0.5, 'Energy')

Compare number of evaluations

for condition, result_full in results_full.items():
    print(condition)
    print('Total evaluations for ' + condition + ':')
    sum = 0
    for key in result_full:
        if condition is not 'np':
            sum += result_full[key]['raw_result']['cost_function_evals']
        else:
            sum = 0
    print(sum)

no bootstrapping
Total evaluations for no bootstrapping:
2741
bootstrapping
Total evaluations for bootstrapping:
952
np
Total evaluations for np:
0

Extrapolation

Here, an extrapolator added that will try to fit each (param,point) set to some specified function and suggest an initial parameter set for the next point (degree of freedom).

• Extrapolator is based on an external extrapolator that sets the ‘window’, that is, the number of previous points to use for extrapolation, while the internal extrapolator proceeds with the actual extrapolation.

• In practice, the user sets the window by specifying an integer value to num_bootstrap - which is also the number of previous points to use for bootstrapping. Additionally, the external extrapolator defines the space within how to extrapolate - here PCA, clustering and the standard window approach.

In practice, if no extrapolator is defined and bootstrapping is True, then all previous points will be bootstrapped. If an extrapolator list is defined and no points are specified for bootstrapping, then the extrapolation will be done based on all previous points.

1. Window Extrapolator: An extrapolator which wraps another extrapolator, limiting the internal extrapolator’s ground truth parameter set to a fixed window size

2. PCA Extrapolator: A wrapper extrapolator which reduces the points’ dimensionality with PCA, performs extrapolation in the transformed pca space, and untransforms the results before returning.

3. Sieve Extrapolator: A wrapper extrapolator which performs an extrapolation, then clusters the extrapolated parameter values into two large and small clusters, and sets the small clusters’ parameters to zero.

4. Polynomial Extrapolator: An extrapolator based on fitting each parameter to a polynomial function of a user-specified degree.

5. Differential Extrapolator: An extrapolator based on treating each param set as a point in space, and performing regression to predict the param set for the next point. A user-specified degree also adds derivatives to the values in the point vectors which serve as features in the training data for the linear regression.

Here we test two different extrapolation techniques

# define different extrapolators
degree = 1
extrap_poly = Extrapolator.factory("poly", degree = degree)
extrap_diff = Extrapolator.factory("diff_model", degree = degree)
extrapolators = {'poly': extrap_poly, 'diff_model': extrap_diff}

for key in extrapolators:
    extrap_internal = extrapolators[key]
    extrap = Extrapolator.factory("window", extrapolator = extrap_internal)
    # define extrapolator
    # BOPES sampler
    bs_extr = BOPESSampler(
        gss=me_gsc
        ,bootstrap=True
        ,num_bootstrap=None
        ,extrapolator=extrap)
    # execute
    res = bs_extr.sample(driver, points)

    results_full[f'{key}']= res.raw_results
    results[f'points_{key}'] = res.points
    results[f'energies_{key}'] = res.energies

Plot results

fig = plt.figure()
for value, bootstrap in conditions.items():
    plt.plot(results[f'points_{bootstrap}'], results[f'energies_{bootstrap}'], label = f'{bootstrap}')
plt.plot(results['points_np'], results['energies_np'], label = 'numpy')
for condition in extrapolators.keys():
    print(condition)
    plt.plot(results[f'points_{condition}'], results[f'energies_{condition}'], label = condition)
plt.legend()
plt.title('Dissociation profile')
plt.xlabel('Interatomic distance')
plt.ylabel('Energy')

poly
diff_model

Text(0, 0.5, 'Energy')

Compare number of evaluations

for condition, result_full in results_full.items():
    print(condition)
    print('Total evaluations for ' + condition + ':')
    sum = 0
    for key in results_full[condition].keys():
        if condition is not 'np':
                sum += results_full[condition][key]['raw_result']['cost_function_evals']
        else:
                sum = 0
    print(sum)

no bootstrapping
Total evaluations for no bootstrapping:
2741
bootstrapping
Total evaluations for bootstrapping:
952
np
Total evaluations for np:
0
poly
Total evaluations for poly:
542
diff_model
Total evaluations for diff_model:
940

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

Version Information

Qiskit Software	Version
Qiskit	0.25.4
Terra	0.17.2
Aer	0.8.2
Ignis	0.6.0
Aqua	0.9.1
IBM Q Provider	0.12.3
System information
Python	3.7.7 (default, Apr 22 2020, 19:15:10) [GCC 9.3.0]
OS	Linux
CPUs	32
Memory (Gb)	125.71903228759766
Tue May 25 16:18:41 2021 EDT

This code is a part of Qiskit

© Copyright IBM 2017, 2021.

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.