Source code for openpathsampling.engines.toy.pes

import numpy as np

from openpathsampling.netcdfplus import StorableObject


# The decorator @restores_ allows us to restore the object from a JSON
# string completely and can thus be stored automatically

[docs]class PES(StorableObject): """Abstract base class for toy potential energy surfaces. """ # For now, we only support additive combinations; maybe someday that can # include multiplication, too
[docs] def __init__(self): super(PES, self).__init__()
def __add__(self, other): return PES_Add(self, other) def __sub__(self, other): return PES_Sub(self, other) def kinetic_energy(self, sys): """Default kinetic energy implementation. Parameters ---------- sys : :class:`.ToyEngine` engine contains its state, including velocities and masses """ v = sys.velocities m = sys.mass return 0.5*np.dot(m, np.multiply(v, v))
class PES_Combination(PES): """Mathematical combination of two potential energy surfaces. Abstract base class. Parameters ---------- pes1 : :class:`.PES` first potential energy surface of the combination pes2 : :class:`.PES` second potential energy surface of the combination fcn : function of two variables function to combine the PES energies dfdx_fcn : function of two variables function to combine the PES (first) derivatives """ def __init__(self, pes1, pes2, fcn, dfdx_fcn): super(PES_Combination, self).__init__() self.pes1 = pes1 self.pes2 = pes2 self._fcn = fcn self._dfdx_fcn = dfdx_fcn def V(self, sys): """Potential energy Parameters ---------- sys : :class:`.ToyEngine` engine contains its state, including velocities and masses Returns ------- float the potential energy """ return self._fcn(self.pes1.V(sys), self.pes2.V(sys)) def dVdx(self, sys): """Derivative of potential energy (-force) Parameters ---------- sys : :class:`.ToyEngine` engine contains its state, including velocities and masses Returns ------- np.array the derivatives of the potential at this point """ return self._dfdx_fcn(self.pes1.dVdx(sys), self.pes2.dVdx(sys))
[docs]class PES_Sub(PES_Combination): """Difference of two potential energy surfaces; pes1 - pes2 Parameters ---------- pes1 : :class:`.PES` first potential energy surface of the combination pes2 : :class:`.PES` second potential energy surface of the combination """
[docs] def __init__(self, pes1, pes2): super(PES_Sub, self).__init__( pes1, pes2, lambda a, b: a - b, lambda a, b: a - b )
[docs]class PES_Add(PES_Combination): """Sum of two potential energy surfaces; pes1 + pes 2 Parameters ---------- pes1 : :class:`.PES` first potential energy surface of the combination pes2 : :class:`.PES` second potential energy surface of the combination """
[docs] def __init__(self, pes1, pes2): super(PES_Add, self).__init__( pes1, pes2, lambda a, b: a + b, lambda a, b: a + b )
[docs]class HarmonicOscillator(PES): r"""Simple harmonic oscillator. Independent in each degree of freedom. :math:`V(x) = \sum_i A_i * mass_i * omega_i^2 * (x_i - x0_i)^2` Parameters ---------- A : list of float omega : list of float x0 : list of float """
[docs] def __init__(self, A, omega, x0): super(HarmonicOscillator, self).__init__() self.A = np.array(A) self.omega = np.array(omega) self.x0 = np.array(x0)
def __repr__(self): # pragma: no cover repr_str = "HarmonicOscillator({obj.A}, {obj.omega}, {obj.x0})" return repr_str.format(obj=self) def to_dict(self): dct = super(HarmonicOscillator, self).to_dict() dct['A'] = dct['A'].tolist() dct['omega'] = dct['omega'].tolist() dct['x0'] = dct['x0'].tolist() return dct def V(self, sys): """Potential energy Parameters ---------- sys : :class:`.ToyEngine` engine contains its state, including velocities and masses Returns ------- float the potential energy """ dx = sys.positions - self.x0 k = self.omega*self.omega*sys.mass return 0.5*np.dot(self.A * k, dx * dx) def dVdx(self, sys): """Derivative of potential energy (-force) Parameters ---------- sys : :class:`.ToyEngine` engine contains its state, including velocities and masses Returns ------- np.array the derivatives of the potential at this point """ dx = sys.positions - self.x0 k = self.omega*self.omega*sys.mass return self.A*k*dx
[docs]class Gaussian(PES): r"""Gaussian given by: :math:`A*exp(-\sum_i alpha[i]*(x[i]-x0[i])^2)` Parameters ---------- A : float amplitude of the Gaussian alpha : list of float Gaussian width parameter x0 : list of float center of the Gaussian """
[docs] def __init__(self, A, alpha, x0): super(Gaussian, self).__init__() self.A = A self.alpha = np.array(alpha) self.x0 = np.array(x0) self._local_dVdx = np.zeros(self.x0.size)
def to_dict(self): dct = super(Gaussian, self).to_dict() dct['alpha'] = dct['alpha'].tolist() dct['x0'] = dct['x0'].tolist() return dct def __repr__(self): # pragma: no cover return "Gaussian({o.A}, {o.alpha}, {o.x0})".format(o=self) def V(self, sys): """Potential energy Parameters ---------- sys : :class:`.ToyEngine` engine contains its state, including velocities and masses Returns ------- float the potential energy """ dx = sys.positions - self.x0 return self.A*np.exp(-np.dot(self.alpha, np.multiply(dx, dx))) def dVdx(self, sys): """Derivative of potential energy (-force) Parameters ---------- sys : :class:`.ToyEngine` engine contains its state, including velocities and masses Returns ------- np.array the derivatives of the potential at this point """ dx = sys.positions - self.x0 exp_part = self.A*np.exp(-np.dot(self.alpha, np.multiply(dx, dx))) for i in range(len(dx)): self._local_dVdx[i] = -2*self.alpha[i]*dx[i]*exp_part return self._local_dVdx
[docs]class OuterWalls(PES): r"""Creates an x**6 barrier around the system. :math:`V(x) = \sum_i sigma_i * (x_i - x0_i)^6` Parameters ---------- sigma : list of float linear scaling of the potential wall x0 : list of float center of the potential """
[docs] def __init__(self, sigma, x0): super(OuterWalls, self).__init__() self.sigma = np.array(sigma) self.x0 = np.array(x0) self._local_dVdx = np.zeros(self.x0.size)
def to_dict(self): dct = super(OuterWalls, self).to_dict() dct['x0'] = dct['x0'].tolist() dct['sigma'] = dct['sigma'].tolist() return dct def __repr__(self): # pragma: no cover return "OuterWalls({o.sigma}, {o.x0})".format(o=self) def V(self, sys): """Potential energy Parameters ---------- sys : :class:`.ToyEngine` engine contains its state, including velocities and masses Returns ------- float the potential energy """ dx = sys.positions - self.x0 myV = 0.0 for i in range(len(dx)): myV += self.sigma[i]*dx[i]**6 return myV def dVdx(self, sys): """Derivative of potential energy (-force) Parameters ---------- sys : :class:`.ToyEngine` engine contains its state, including velocities and masses Returns ------- np.array the derivatives of the potential at this point """ dx = sys.positions - self.x0 for i in range(len(dx)): self._local_dVdx[i] = 6.0*self.sigma[i]*dx[i]**5 return self._local_dVdx
[docs]class LinearSlope(PES): r"""Linear potential energy surface. :math:`V(x) = \sum_i m_i * x_i + c` Parameters ---------- m : list of float slope c : float energy offset """
[docs] def __init__(self, m, c): super(LinearSlope, self).__init__() self.m = m self.c = c self._local_dVdx = self.m self.dim = len(self.m)
def __repr__(self): # pragma: no cover return "LinearSlope({o.m}, {o.c})".format(o=self) def V(self, sys): """Potential energy Parameters ---------- sys : :class:`.ToyEngine` engine contains its state, including velocities and masses Returns ------- float the potential energy """ return np.dot(self.m, sys.positions) + self.c def dVdx(self, sys): """Derivative of potential energy (-force) Parameters ---------- sys : :class:`.ToyEngine` engine contains its state, including velocities and masses Returns ------- np.array the derivatives of the potential at this point """ # this is independent of the position return self._local_dVdx
class DoubleWell(PES): r"""Simple double-well potential. Independent in each degree of freedom. V(x) = \sum_i A_i * (x_i**2 - x0_i**2)**2 WARNING: Two minima only in one dimension, otherwise there are more! Parameters ---------- A : list of float potential prefactor for in each degree of freedom. x0 : list of float minimum position in each degree of freedom. """ def __init__(self, A, x0): super(DoubleWell, self).__init__() self.A = np.array(A) self.x0 = np.array(x0) def __repr__(self): # pragma: no cover repr_str = "DoubleWell({obj.A}, {obj.x0})" return repr_str.format(obj=self) def to_dict(self): dct = super(DoubleWell, self).to_dict() dct['A'] = dct['A'].tolist() dct['x0'] = dct['x0'].tolist() return dct def V(self, sys): """Potential energy Parameters ---------- sys : :class:`.ToyEngine` engine contains its state, including velocities and masses Returns ------- float the potential energy """ dx2 = sys.positions * sys.positions - self.x0 * self.x0 return np.dot(self.A, dx2 * dx2) def dVdx(self, sys): """Derivative of potential energy (-force) Parameters ---------- sys : :class:`.ToyEngine` engine contains its state, including velocities and masses Returns ------- np.array the derivatives of the potential at this point """ dx2 = sys.positions * sys.positions - self.x0 * self.x0 return 4 * self.A * sys.positions * dx2