BOPES Sampler algorithm for Chemistry

.pylintdict

@@ -41,6 +41,7 @@ backend
 backends
 baritompa
 barkoutsos
+basedriver
 bb
 benchmarking
 bergholm
@@ -167,6 +168,7 @@ ee
 egger
 eigen
 eigenenergies
+eigensolution
 eigensolver
 eigensolvers
 eigenstate
@@ -203,10 +205,14 @@ exponentials
 exponentiated
 exponentiations
 expr
+extrapolator
+extrapolators
+extrapolator's
 factorizers
 factr
 fcidump
 fcompiler
+fd
 fermionic
 fermionicoperator
 fermions
@@ -249,6 +255,7 @@ grinko
 groundenergy
 groundstates
 grover
+gsc
 gset
 gsls
 gto
@@ -259,6 +266,7 @@ hamiltonian
 hamiltonians
 hardcoded
 hartree
+hartrees
 hassidim
 hb
 hhl
@@ -301,6 +309,7 @@ instantiation
 integercomparator
 intelvem
 interatomic
+internuclear
 ints
 invertible
 ip
@@ -362,6 +371,7 @@ lnot
 loglik
 loglikelihood
 logn
+lookback
 lor
 lowerbound
 lp
@@ -408,6 +418,7 @@ moller
 momentums
 monomial
 monomials
+morse
 mosca
 motived
 mprev
@@ -680,6 +691,7 @@ terra
 terra's
 tgt
 th
+thermo
 tiajb
 tibja
 timelimit
@@ -768,6 +780,7 @@ xuxv
 xyz
 yamamoto
 yao
+ydata
 YFinance
 yu
 yy

qiskit/chemistry/algorithms/pes_samplers/__init__.py

@@ -0,0 +1,18 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 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.
+
+"""Potential energy surface samplers"""
+
+from .bopes_sampler import BOPESSampler
+
+__all__ = ['BOPESSampler'
+           ]

qiskit/chemistry/algorithms/pes_samplers/bopes_sampler.py

@@ -0,0 +1,196 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 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.
+
+"""The calculation of points on the Born-Oppenheimer Potential Energy Surface (BOPES)."""
+
+import logging
+from typing import Optional, List, Dict
+
+import numpy as np
+from qiskit.aqua import AquaError
+from qiskit.aqua.algorithms import VQAlgorithm
+from qiskit.chemistry.drivers import BaseDriver
+from qiskit.chemistry.algorithms.ground_state_solvers import GroundStateSolver
+from qiskit.chemistry.results.bopes_sampler_result import BOPESSamplerResult
+from qiskit.chemistry.algorithms.pes_samplers.extrapolator import Extrapolator, WindowExtrapolator
+from qiskit.chemistry.results import EigenstateResult
+
+logger = logging.getLogger(__name__)
+
+
+class BOPESSampler:
+    """Class to evaluate the Born-Oppenheimer Potential Energy Surface (BOPES)."""
+
+    def __init__(self,
+                 gss: GroundStateSolver,
+                 tolerance: float = 1e-3,
+                 bootstrap: bool = True,
+                 num_bootstrap: Optional[int] = None,
+                 extrapolator: Optional[Extrapolator] = None) -> None:
+        """
+        Args:
+            gss: GroundStateSolver
+            tolerance: Tolerance desired for minimum energy.
+            bootstrap: Whether to warm-start the solution of variational minimum eigensolvers.
+            num_bootstrap: Number of previous points for extrapolation
+                and bootstrapping. If None and a list of extrapolators is defined,
+                the first two points will be used for bootstrapping.
+                If no extrapolator is defined and bootstrap is True,
+                all previous points will be used for bootstrapping.
+            extrapolator: Extrapolator objects that define space/window
+                           and method to extrapolate variational parameters.
+
+        Raises:
+            AquaError: If ``num_boostrap`` is an integer smaller than 2, or
+                if ``num_boostrap`` is larger than 2 and the extrapolator is not an instance of
+                ``WindowExtrapolator``.
+        """
+
+        self._gss = gss
+        self._tolerance = tolerance
+        self._bootstrap = bootstrap
+        self._driver = None     # type: Optional[BaseDriver]
+        self._points = None     # type: Optional[List[float]]
+        self._energies = None   # type: Optional[List[float]]
+        self._raw_results = None        # type: Optional[Dict[float, EigenstateResult]]
+        self._points_optparams = None   # type: Optional[Dict[float, List[float]]]
+        self._num_bootstrap = num_bootstrap
+        self._extrapolator = extrapolator
+
+        if self._extrapolator:
+            if num_bootstrap is None:
+                # set default number of bootstrapping points to 2
+                self._num_bootstrap = 2
+            elif num_bootstrap >= 2:
+                if not isinstance(self._extrapolator, WindowExtrapolator):
+                    raise AquaError(
+                        'If num_bootstrap >= 2 then the extrapolator must be an instance '
+                        'of WindowExtrapolator, got {} instead'.format(self._extrapolator))
+                self._num_bootstrap = num_bootstrap
+                self._extrapolator.window = num_bootstrap  # window for extrapolator
+            else:
+                raise AquaError(
+                    'num_bootstrap must be None or an integer greater than or equal to 2')
+
+        if isinstance(self._gss.solver, VQAlgorithm):
+            # Save initial point passed to min_eigensolver;
+            # this will be used when NOT bootstrapping
+            self._initial_point = self._gss.solver.initial_point
+
+    def sample(self, driver: BaseDriver, points: List[float]) -> BOPESSamplerResult:
+        """Run the sampler at the given points, potentially with repetitions.
+
+        Args:
+            driver: BaseDriver specific for the problem. The driver should be based on
+                    a Molecule object that has perturbations to be varied.
+            points: The points along the degrees of freedom to evaluate.
+
+        Returns:
+            BOPES Sampler Result
+
+        Raises:
+            AquaError: if the driver does not have a molecule specified.
+        """
+
+        self._driver = driver
+
+        if self._driver.molecule is None:
+            raise AquaError('Driver MUST be configured with a Molecule.')
+
+        # full dictionary of points
+        self._raw_results = self._run_points(points)
+        # create results dictionary with (point, energy)
+        self._points = list(self._raw_results.keys())
+        self._energies = []
+        for key in self._raw_results:
+            energy = self._raw_results[key].computed_electronic_energy + \
+                     self._raw_results[key].nuclear_repulsion_energy
+            self._energies.append(energy)
+
+        result = BOPESSamplerResult(self._points, self._energies, self._raw_results)
+
+        return result
+
+    def _run_points(self, points: List[float]) -> Dict[float, EigenstateResult]:
+        """Run the sampler at the given points.
+
+        Args:
+            points: the points along the single degree of freedom to evaluate
+
+        Returns:
+            The results for all points.
+        """
+        raw_results = dict()   # type: Dict[float, EigenstateResult]
+        if isinstance(self._gss.solver, VQAlgorithm):
+            self._points_optparams = dict()
+            self._gss.solver.initial_point = self._initial_point
+
+        # Iterate over the points
+        for i, point in enumerate(points):
+            logger.info('Point %s of %s', i + 1, len(points))
+            raw_result = self._run_single_point(point)  # dict of results
+            raw_results[point] = raw_result
+
+        return raw_results
+
+    def _run_single_point(self, point: float) -> EigenstateResult:
+        """Run the sampler at the given single point
+
+        Args:
+            point: The value of the degree of freedom to evaluate.
+
+        Returns:
+            Results for a single point.
+        """
+
+        # update molecule geometry and  thus resulting Hamiltonian based on specified point
+        self._driver.molecule.perturbations = [point]
+
+        # find closest previously run point and take optimal parameters
+        if isinstance(self._gss.solver, VQAlgorithm) and self._bootstrap:
+            prev_points = list(self._points_optparams.keys())
+            prev_params = list(self._points_optparams.values())
+            n_pp = len(prev_points)
+
+            # set number of points to bootstrap
+            if self._extrapolator is None:
+                n_boot = len(prev_points)  # bootstrap all points
+            else:
+                n_boot = self._num_bootstrap
+
+            # Set initial params # if prev_points not empty
+            if prev_points:
+                if n_pp <= n_boot:
+                    distances = np.array(point) - \
+                                np.array(prev_points).reshape(n_pp, -1)
+                    # find min 'distance' from point to previous points
+                    min_index = np.argmin(np.linalg.norm(distances, axis=1))
+                    # update initial point
+                    self._gss.solver.initial_point = prev_params[min_index]
+                else:  # extrapolate using saved parameters
+                    opt_params = self._points_optparams
+                    param_sets = self._extrapolator.extrapolate(points=[point],
+                                                                param_dict=opt_params)
+                    # update initial point, note param_set is a list
+                    # param set is a dictionary
+                    self._gss.solver.initial_point = param_sets.get(point)
+
+        # the output is an instance of EigenstateResult
+        result = self._gss.solve(self._driver)
+
+        # Save optimal point to bootstrap
+        if isinstance(self._gss.solver, VQAlgorithm):
+            # at every point evaluation, the optimal params are updated
+            optimal_params = self._gss.solver.optimal_params
+            self._points_optparams[point] = optimal_params
+
+        return result

qiskit/chemistry/algorithms/pes_samplers/energy_surface_spline.py

@@ -0,0 +1,96 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 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.
+
+"""
+An spline interpolation method for data fitting
+"""
+
+import scipy.interpolate as interp
+from scipy.optimize import minimize_scalar
+
+from qiskit.chemistry.algorithms.pes_samplers.potentials.potential_base import EnergySurfaceBase
+
+
+class EnergySurface1DSpline(EnergySurfaceBase):
+    """
+    A simple cubic spline interpolation for the potential energy surface.
+    """
+
+    def __init__(self):
+        """
+        Constructor.
+        """
+        super().__init__()
+        self._eval = None
+        self.eval_d = None
+        self.min_x = None
+        self.min_val = None
+        self.x_left = None
+        self.x_right = None
+
+    # Implementing the EnergySurfaceBase interface
+
+    def eval(self, x):
+        """
+        After fitting the data to the fit function, predict the energy
+            at a point x.
+        """
+        assert self._eval is not None
+        result = self._eval(x)
+
+        return result
+
+    def fit(self, xdata, ydata, initial_vals=None, bounds_list=None):
+
+        newx = xdata
+        newy = ydata
+
+        tck = interp.splrep(newx, newy, k=3)
+
+        self._eval = lambda x: interp.splev(x, tck)
+        self.eval_d = lambda x: interp.splev(x, tck, der=1)
+
+        result = minimize_scalar(self._eval)
+        assert result.success
+
+        self.min_x = result.x
+        self.min_val = result.fun
+        self.x_left = min(xdata)
+        self.x_right = max(xdata)
+
+    def get_equilibrium_geometry(self, scaling=1.0):
+        """
+        Returns the geometry for the minimal energy (scaled by 'scaling')
+        Default units (scaling=1.0) are Angstroms. Scale by 1E-10 to get
+        meters.
+        """
+        assert self.min_x is not None
+        return self.min_x * scaling
+
+    def get_minimal_energy(self, scaling=1.0):
+        """
+        Returns the value of the minimal energy (scaled by 'scaling')
+        Default units (scaling=1.0) are J/mol. Scale appropriately for
+        Hartrees.
+        """
+        assert self.min_val is not None
+        return self.min_val * scaling
+
+    def get_trust_region(self):
+        """
+        Returns the bounds of the region (in space) where the energy
+        surface implementation can be trusted. When doing spline
+        interpolation, for example, that would be the region where data
+        is interpolated (vs. extrapolated) from the arguments of
+        fit().
+        """
+        return (self.x_left, self.x_right)