Notebook

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).

In [1]:

# 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)

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

In [2]:

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

In [3]:

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)

In [4]:

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¶

In [5]:

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.8999999999999999])]
[('H', [0.0, 0.0, 0.0]), ('H', [0.0, 0.0, 0.5])]

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 boostrapping is set to True, all previous optimal parameters are used as initial points for the next runs.

In [81]:

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)

# Specify degree of freedom (points of interest)
points = np.linspace(0.1, 2, 20)
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

Compare to classical eigensolver¶

In [82]:

# 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

Plot results¶

In [83]:

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')

Out[83]:

Text(0, 0.5, 'Energy')

Compare number of evaluations¶

In [84]:

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:
1714
bootstrapping
Total evaluations for bootstrapping:
1266
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 specifies an integer value to num_bootstrap - which 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.

Window Extrapolator: An extrapolator which wraps another extrapolator, limiting the internal extrapolator's ground truth parameter set to a fixed window size
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.
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.
Polynomial Extrapolator: An extrapolator based on fitting each parameter to a polynomial function of a user-specified degree.
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¶

In [85]:

# 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¶

In [86]:

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

Out[86]:

Text(0, 0.5, 'Energy')

Compare number of evaluations¶

In [87]:

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:
1714
bootstrapping
Total evaluations for bootstrapping:
1266
np
Total evaluations for np:
0
poly
Total evaluations for poly:
715
diff_model
Total evaluations for diff_model:
928

In [88]:

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

Version Information

Qiskit Software	Version
Qiskit	0.23.0
Terra	0.16.0
Aer	0.7.0
Ignis	0.5.0
Aqua	0.8.0
IBM Q Provider	0.11.0
System information
Python	3.6.12 \|Anaconda, Inc.\| (default, Sep 8 2020, 17:50:39) [GCC Clang 10.0.0 ]
OS	Darwin
CPUs	4
Memory (Gb)	16.0
Mon Oct 19 16:12:19 2020 CEST

This code is a part of Qiskit

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.