integrators

class ufedmm.integrators.AbstractMiddleRespaIntegrator(temperature, step_size, num_rattles=0, scheme='VV-Middle', respa_loops=[1], bath_loops=1, intertwine=True, embodied_force_groups=[], unroll_loops=True)

Bases: CustomIntegrator

An abstract class for middle-type, multiple time-scale integrators.

Warning

This class is meant for inheritance only and does not actually include thermostatting. Concrete subclasses are available, such as MiddleMassiveNHCIntegrator and GeodesicLangevinIntegrator, for instance.

Child classes will differ by the thermostat algorithm, which must be implemented by overriding the _bath method (see the example below). Temperature is treated as a per-dof parameter so as to allow adiabatic simulations.

The following ODE system is solved for every degree of freedom in the system, with possibly \(n_c\) holonomic constraints and with forces possibly split into \(m\) parts according to their characteristic time scales:

\[\begin{split}& \dot{r}_i = v_i \\ & \dot{v}_i = \frac{\sum_{k=1}^m F_i^{[k]}}{m_i} + \sum_{k=1}^{n_c} \lambda_k \nabla_{r_i} \sigma_k + \mathrm{bath}(T_i, v_i) \\ & \sigma_k(\mathbf{r}) = 0\end{split}\]

An approximate solution is obtained by applying the Trotter-Suzuki splitting formula. In the particular case of two time scales, the default splitting scheme goes as follows:

\[e^{\Delta t\mathcal{L}} = e^{\frac{\Delta t}{2}\mathcal{L}^{[1]}_v} \left[ e^{\frac{\Delta t}{2 n_0}\mathcal{L}^{[0]}_v} \left( e^{\frac{\Delta t}{2 n_0 n_b}\mathcal{L}_r} e^{\frac{\Delta t}{n_0 n_b}\mathcal{L}_\mathrm{bath}} e^{\frac{\Delta t}{2 n_0 n_b}\mathcal{L}_r} \right)^{n_b} e^{\frac{\Delta t}{2 n_0}\mathcal{L}^{[0]}_v} \right]^{n_0} e^{\frac{\Delta t}{2}\mathcal{L}^{[1]}_v}\]

Each exponential operator is the solution of a particular subsystem of equations. If \(\mathrm{bath}(T_i, v_i) = 0\), the scheme above is time-reversible, measure-preserving, and symplectic. It is referred to as the VV-Middle scheme [2], where VV stands for Velocity Verlet. An alternative approach is also available, which is:

\[e^{\Delta t\mathcal{L}} = \left[ \left( e^{\frac{\Delta t}{2 n_0 n_b}\mathcal{L}_r} e^{\frac{\Delta t}{n_0 n_b}\mathcal{L}_\mathrm{bath}} e^{\frac{\Delta t}{2 n_0 n_b}\mathcal{L}_r} \right)^{n_b} e^{\frac{\Delta t}{n_0}\mathcal{L}^{[0]}_v} \right]^{n_0} e^{\Delta t \mathcal{L}^{[1]}_v}\]

This is referred to as the LF-Middle scheme [2], where LF stands for Leap-Frog. In contrast to the previous scheme, it is not time-reversible. However, in single time-scale simulations, the two approaches result in equivalent coordinate trajectories, while the latter provides a velocity trajectory more consistent with the Maxwell-Boltzmann distribution at the specified temperature [2].

Parameters:

• temperature (unit.Quantity) – The temperature of the heat bath.

• step_size (float or unit.Quantity) – The outer step size with which to integrate the equations of motion.

Keyword Arguments:

• num_rattles (int, default=0) – The number of RATTLE computations for geodesic integration [3]. If num_rattles=0, then no constraints are considered at all.

• scheme (str, default='VV-Middle') – Which splitting scheme will be used. Valid options are ‘VV-Middle’ and ‘LF-Middle’.

• respa_loops (list(int), default=[1]) – A list of m integers, where respa_loops[k] determines how many substeps with force group k are internally executed for every step with force group k+1.

• bath_loops (int, default=1) – The number of iterations of the bath operator per each step at time scale 0. This is useful when the bath operator is not exact, but derived from a splitting solution.

• embodied_force_groups (list(int), default=[]) – A list of indices of force groups. The presence of an index i is this list means that the contribution of force group i is embodied in force group i+1. Therefore, such contribution must be properly subtracted during the integration at time scale i+1. This feature requires OpenMM 7.5 or a newer version.

• unroll_loops (bool, default=True) – Whether the integrator loops should be unrolled for improving efficiency. Using unroll_loops=False can be useful for printing the integrator steps.

Example

>>> from ufedmm import integrators
>>> from openmm import unit
>>> class MiddleNoseHooverIntegrator(integrators.AbstractMiddleRespaIntegrator):
...     def __init__(self, ndof, tau, temperature, step_size, num_rattles=1):
...         super().__init__(
...             temperature, step_size, num_rattles, 'VV-Middle', [1], 1
...         )
...         kB = 8.3144626E-3*unit.kilojoules_per_mole/unit.kelvin
...         gkT = ndof*unit.MOLAR_GAS_CONSTANT_R*temperature
...         self.addGlobalVariable('gkT', gkT)
...         self.addGlobalVariable('Q', gkT*tau**2)
...         self.addGlobalVariable('v_eta', 0)
...         self.addGlobalVariable('twoK', 0)
...         self.addGlobalVariable('scaling', 1)
...     def _bath(self, fraction):
...         self.addComputeSum('twoK', 'm*v*v')
...         self.addComputeGlobal(
...             'v_eta', f'v_eta + {0.5*fraction}*dt*(twoK - gkT)/Q'
...         )
...         self.addComputeGlobal('scaling', f'exp(-{fraction}*dt*v_eta)')
...         self.addComputePerDof('v', f'v*scaling')
...         self.addComputeGlobal(
...             'v_eta', f'v_eta + {0.5*fraction}*dt*(scaling^2*twoK - gkT)/Q'
...         )
>>> integrator = MiddleNoseHooverIntegrator(
...     500, 10*unit.femtoseconds, 300*unit.kelvin,
...     1*unit.femtoseconds, num_rattles=0
... )
>>> print(integrator)
Per-dof variables:
  kT
Global variables:
  gkT = 1247.169392722986
  Q = 0.1247169392722986
  v_eta = 0.0
  twoK = 0.0
  scaling = 1.0
Computation steps:
   0: allow forces to update the context state
   1: v <- v + 0.5*dt*f/m
   2: x <- x + 0.5*dt*v
   3: twoK <- sum(m*v*v)
   4: v_eta <- v_eta + 0.5*dt*(twoK - gkT)/Q
   5: scaling <- exp(-1.0*dt*v_eta)
   6: v <- v*scaling
   7: v_eta <- v_eta + 0.5*dt*(scaling^2*twoK - gkT)/Q
   8: x <- x + 0.5*dt*v
   9: v <- v + 0.5*dt*f/m

class ufedmm.integrators.CustomIntegrator(temperature, step_size)

Bases: CustomIntegrator

An extension of openmm.CustomIntegrator class with an extra per-dof variable named temperature, whose content is the temperature of the heat bath associated to each degree of freedom. A per-dof temperature is necessary if the extended-space variables and the physical system are coupled adiabatically to thermostats at different temperatures. Otherwise, any other OpenMM integrator can be used.

Parameters:

• temperature (float or unit.Quantity) – The temperature.

• step_size (float or unit.Quantity) – The step size with which to integrate the equations of motion.

step(self, steps)

Advance a simulation through time by taking a series of time steps.

Parameters:

steps (int) – the number of time steps to take

class ufedmm.integrators.GeodesicLangevinIntegrator(temperature, friction_coefficient, step_size, num_rattles=1, scheme='LF-Middle', **kwargs)

Bases: AbstractMiddleRespaIntegrator

A geodesic Langevin integrator [3], which can be integrated by using either the LF-Middle or the VV-Middle scheme [2].

Parameters:

• temperature (float or unit.Quantity) – The temperature.

• friction_coefficient (float or unit.Quantity) – The friction coefficient.

• step_size (float or unit.Quantity) – The time-step size.

Keyword Arguments:

• num_rattles (int, default=1) – The number of RATTLE computations for geodesic integration [3]. If num_rattles=0, then no constraints are considered at all.

• scheme (str, default='LF-Middle') – Which splitting scheme will be used. Valid options are ‘VV-Middle’ and ‘LF-Middle’.

• **kwargs – All other keyword arguments in AbstractMiddleRespaIntegrator.

Example

>>> import ufedmm
>>> dt = 2*unit.femtoseconds
>>> temp = 300*unit.kelvin
>>> gamma = 10/unit.picoseconds
>>> ufedmm.GeodesicLangevinIntegrator(
...     temp, gamma, dt, num_rattles=1, scheme='VV-Middle'
... )
Per-dof variables:
  kT, x0
Global variables:
  friction = 10.0
Computation steps:
   0: allow forces to update the context state
   1: v <- v + 0.5*dt*f/m
   2: constrain velocities
   3: x <- x + 0.5*dt*v
   4: x0 <- x
   5: constrain positions
   6: v <- v + (x - x0)/(0.5*dt)
   7: constrain velocities
   8: v <- z*v + sqrt((1 - z*z)*kT/m)*gaussian; z = exp(-friction*1.0*dt)
   9: constrain velocities
  10: x <- x + 0.5*dt*v
  11: x0 <- x
  12: constrain positions
  13: v <- v + (x - x0)/(0.5*dt)
  14: constrain velocities
  15: v <- v + 0.5*dt*f/m
  16: constrain velocities

class ufedmm.integrators.MiddleMassiveGGMTIntegrator(temperature, time_constant, step_size, **kwargs)

Bases: AbstractMiddleRespaIntegrator

A massive, middle-type Generalized Gaussian Moment Thermostat [4] solver with optional multiple time-scale integration via RESPA.

To enable RESPA, the forces in OpenMM system must be split into distinct force groups and the keyword respa_loop (see below) must be a list with multiple entries.

Parameters:

• temperature (float or unit.Quantity) – The temperature.

• time_constant (float or unit.Quantity) – The characteristic time constant.

• step_size (float or unit.Quantity) – The time-step size.

Keyword Arguments:

**kwargs – All keyword arguments in AbstractMiddleRespaIntegrator, except num_rattles.

Example

>>> import ufedmm
>>> temp, tau, dt = 300*unit.kelvin, 10*unit.femtoseconds, 2*unit.femtoseconds
>>> integrator = ufedmm.MiddleMassiveGGMTIntegrator(temp, tau, dt)
>>> print(integrator)
Per-dof variables:
  kT, Q1, Q2, v1, v2
Computation steps:
   0: allow forces to update the context state
   1: v <- v + 0.5*dt*f/m
   2: x <- x + 0.5*dt*v
   3: v1 <- v1 + 0.5*dt*(m*v^2 - kT)/Q1
   4: v2 <- v2 + 0.5*dt*((m*v^2)^2/3 - kT^2)/Q2
   5: v <- v*exp(-1.0*dt*(v1 + kT*v2))/sqrt(1 + 2.0*dt*m*v^2*v2/3)
   6: v1 <- v1 + 0.5*dt*(m*v^2 - kT)/Q1
   7: v2 <- v2 + 0.5*dt*((m*v^2)^2/3 - kT^2)/Q2
   8: x <- x + 0.5*dt*v
   9: v <- v + 0.5*dt*f/m

class ufedmm.integrators.MiddleMassiveNHCIntegrator(temperature, time_constant, step_size, nchain=2, track_energy=False, **kwargs)

Bases: AbstractMiddleRespaIntegrator

A massive, middle-type Nose-Hoover Chain Thermostat solver [5] with optional multiple time-scale integration via RESPA.

To enable RESPA, the forces in OpenMM system must be split into distinct force groups and the keyword respa_loop (see below) must be a list with multiple entries.

Parameters:

• temperature (float or unit.Quantity) – The temperature.

• time_constant (float or unit.Quantity) – The characteristic time constant.

• step_size (float or unit.Quantity) – The time-step size.

Keyword Arguments:

• nchain (int, default=2) – The number of thermostats in each Nose-Hoover chain.

• track_energy (bool, default=False) – Whether to track the thermostat energy term.

• **kwargs – All keyword arguments in AbstractMiddleRespaIntegrator, except num_rattles.

Example

>>> import ufedmm
>>> temp, tau, dt = 300*unit.kelvin, 10*unit.femtoseconds, 2*unit.femtoseconds
>>> integrator = ufedmm.MiddleMassiveNHCIntegrator(
...     temp, tau, dt, respa_loops=[4, 1], unroll_loops=False
... )
>>> print(integrator)
Per-dof variables:
  kT, Q, v1, v2
Global variables:
  irespa0 = 0.0
Computation steps:
   0: allow forces to update the context state
   1: v <- v + 0.5*dt*f1/m
   2: irespa0 <- 0
   3: while (irespa0 < 3.5):
   4:    v <- v + 0.125*dt*f0/m
   5:    x <- x + 0.125*dt*v
   6:    v2 <- v2 + 0.125*dt*(Q*v1^2 - kT)/Q
   7:    v1 <- (v1*z + 0.125*dt*(m*v^2 - kT)/Q)*z; z=exp(-0.0625*dt*v2)
   8:    v <- v*exp(-0.25*dt*v1)
   9:    v1 <- (v1*z + 0.125*dt*(m*v^2 - kT)/Q)*z; z=exp(-0.0625*dt*v2)
  10:    v2 <- v2 + 0.125*dt*(Q*v1^2 - kT)/Q
  11:    x <- x + 0.125*dt*v
  12:    v <- v + 0.125*dt*f0/m
  13:    irespa0 <- irespa0 + 1
  14: end
  15: v <- v + 0.5*dt*f1/m

class ufedmm.integrators.RegulatedNHLIntegrator(temperature, time_constant, friction_coefficient, step_size, regulation_parameter, semi_regulated=True, split_ornstein_uhlenbeck=True, **kwargs)

Bases: AbstractMiddleRespaIntegrator

A regulated version of the massive Nose-Hoover-Langevin [6, 7] method. Regulation means that the system Hamiltonian is modified so that velocities remain below a temperature-dependent limit. This is closely related to the SIN(R) method [8] and allows multiple time-scale integration with very large outer time steps, without resonance.

The following SDE system is solved for every degree of freedom in the system:

\[\begin{split}& dr_i = v_i dt \\ & dp_i = F_i dt - v_{\eta_i} m_i v_i dt \\ & dv_{\eta_i} = \frac{1}{Q}\left(\frac{n+1}{n} m_i v_i^2 - k_B T\right) dt - \gamma v_{\eta_i} dt + \sqrt{\frac{2\gamma k_B T}{Q}} dW_i,\end{split}\]

where:

\[v_i = c_i \tanh\left(\frac{p_i}{m_i c_i}\right).\]

Here, \(n\) is the regulation parameter and \(c_i\) is the maximum speed for degree of freedom i, given by \(c_i = \sqrt{\frac{n k T}{m_i}}\). The inertial parameter \(Q\) is defined as \(Q = n k_B T \tau^2\), with \(\tau\) being a relaxation time [9]. An approximate solution is obtained by applying the Trotter-Suzuki splitting formula:

\[e^{\Delta t\mathcal{L}} = e^{\frac{\Delta t}{2}\mathcal{L}^1_p} \left[e^{\frac{\delta t}{2}\mathcal{L}^0_p} e^{\frac{\delta t}{2}\mathcal{L}_r} e^{\delta t \mathcal{L}_\mathrm{bath}} e^{\frac{\delta t}{2}\mathcal{L}_r} e^{\frac{\delta t}{2}\mathcal{L}^0_p}\right]^m e^{\frac{\Delta t}{2}\mathcal{L}^1_p}\]

where \(\delta t = \frac{\Delta t}{m}\). Each exponential operator above is the solution of a differential equation.

The exact solution for the physical-system part is:

\[r_i(t) = r_i^0 + c_i \mathrm{tanh}\left(\frac{p_i}{m c_i}\right) t\]

\[p_i(t) = p_i^0 + F_i t\]

The bath propagator is further split as:

\[e^{\delta t \mathcal{L}_\mathrm{bath}} = e^{\frac{\delta t}{2m}\mathcal{L}_B} e^{\frac{\delta t}{2m}\mathcal{L}_S} e^{\frac{\delta t}{m}\mathcal{L}_O} e^{\frac{\delta t}{2m}\mathcal{L}_S} e^{\frac{\delta t}{2m}\mathcal{L}_B}\]

Part ‘B’ is a boost, whose solution is:

\[v_{\eta_i}(t) = v_{\eta_i}^0 + \frac{1}{Q}\left(\frac{n+1}{n} m_i v_i^2 - k_B T\right) t\]

Part ‘S’ is a scaling, whose solution is:

\[p_i(t) = m_i c_i \mathrm{arcsinh}\left[ \sinh\left(\frac{p_i^0}{m_i c_i}\right) e^{- v_{\eta_i} t} \right]\]

Part ‘O’ is an Ornstein–Uhlenbeck process, whose solution is:

\[v_{\eta_i}(t) = v_{\eta_i}^0 e^{-\gamma t} + \sqrt{\frac{k_B T}{Q}(1-e^{-2\gamma t})} R_N\]

where \(R_N\) is a normally distributed random number.

Parameters:

• step_size (float or unit.Quantity) – The outer step size with which to integrate the equations of motion.

• loops (int) – The number of internal substeps at each time step.

• temperature (unit.Quantity) – The temperature of the heat bath.

• time_scale (unit.Quantity (time)) – The relaxation time (\(\tau\)) of the Nose-Hoover thermostat.

• friction_coefficient (unit.Quantity (1/time)) – The friction coefficient (\(\gamma\)) of the Langevin thermostat.

• regulation_parameter (int or float) – The regulation parameter n.

Keyword Arguments:

• semi_regulated (bool, default=True) – Whether to use the semi-regulated NHL version of the method instead of its fully-regulated version.

• split_ornstein_uhlenbeck (bool, default=True) – Whether to split the drifted Ornstein-Uhlenbeck operator.

• **kwargs – All keyword arguments in AbstractMiddleRespaIntegrator, except num_rattles.

ufedmm.integrators.add_inner_nonbonded_force(system, inner_switch, inner_cutoff, force_group_index)

Adds a new force group to an openmm.System containing a openmm.NonbondedForce object, with the purpose of performing multiple time-scale integration via the RESPA2 splitting scheme of Morrone, Zhou, and Berne [10]. Also, assigns the provided force_group_index to this new group and force_group_index+1 to the original openmm.NonbondedForce. When used in an instance of any subclass of AbstractMiddleRespaIntegrator, the new force group must be set as embodied by the openmm.NonbondedForce as opposed to being complimentary to it.

Parameters:

• system (openmm.System) – The system the inner force will be added to, which must contain a openmm.NonbondedForce.

• inner_switch (float or unit.Quantity) – The inner switching distance, where the interaction of an atom pair begins to switch off to zero.

• inner_cutoff (float or unit.Quantity) – The inner cutoff distance, where the interaction of an atom pairs completely switches off.

• force_group_index (int) – The force group the new interactions will belong to. The old openmm.NonbondedForce will be automatically assigned to force_group_index+1.

Example

>>> import ufedmm
>>> from openmm import unit
>>> dt = 2*unit.femtoseconds
>>> temp = 300*unit.kelvin
>>> tau = 10*unit.femtoseconds
>>> gamma = 10/unit.picoseconds
>>> model = ufedmm.AlanineDipeptideModel()
>>> ufedmm.add_inner_nonbonded_force(
...     model.system, 5*unit.angstroms, 8*unit.angstroms, 1
... )