Coverage for openxps/integrators/regulated_nhl.py: 100%
59 statements
« prev ^ index » next coverage.py v7.11.3, created at 2025-11-13 22:08 +0000
1"""
2.. module:: openxps.integrators.regulated_nhl
3 :platform: Linux, MacOS, Windows
4 :synopsis: Regulated Nosé-Hoover-Langevin (NHL) integrator.
6.. classauthor:: Charlles Abreu <craabreu@gmail.com>
8"""
10import numpy as np
11import openmm as mm
12from openmm import unit as mmunit
14from openxps.utils import preprocess_args
16from .utils import IntegratorMixin, add_property
19@add_property("temperature", mmunit.kelvin)
20@add_property("time constant", mmunit.picosecond)
21@add_property("friction coefficient", 1 / mmunit.picosecond)
22class RegulatedNHLIntegrator(IntegratorMixin, mm.CustomIntegrator):
23 r"""A Regulated Nosé-Hoover-Langevin (NHL) integrator :cite:`Abreu2021`.
25 The Regulated Nosé-Hoover-Langevin (NHL) integrator is effective at maintaining
26 temperature control under steady-state conditions, such as in d-AFED/TAMD
27 :cite:`Maragliano2006,Abrams2008` or UFED :cite:`Chen2012` simulations.
29 .. note::
30 The friction coefficient is initially set to the inverse of the time constant,
31 but this value can be customized via :meth:`setFrictionCoefficient`.
33 .. warning::
34 This integrator is incompatible with systems subject to constraints.
36 Parameters
37 ----------
38 temperature
39 The temperature of the heat bath.
40 timeConstant
41 The time constant (a.k.a. coupling/damping/relaxation time) of the thermostat.
42 stepSize
43 The integration step size.
45 Keyword Arguments
46 -----------------
47 regulationParameter
48 The regulation parameter (see :cite:`Abreu2021` for details).
49 forceFirst
50 If True, the integrator will apply a force-first scheme. Otherwise, it will
51 apply a symmetric operator splitting scheme.
53 Example
54 -------
55 >>> import openxps as xps
56 >>> from openmm import unit
57 >>> from openmmtools import testsystems
59 Symmetric scheme (default)
61 >>> integrator = xps.integrators.RegulatedNHLIntegrator(
62 ... temperature=300 * unit.kelvin,
63 ... timeConstant=100 * unit.femtoseconds,
64 ... stepSize=2 * unit.femtoseconds,
65 ... )
66 >>> integrator
67 Per-dof variables:
68 v1
69 Global variables:
70 kT = 2.494...
71 invQ = 40.090...
72 a = 0.9801...
73 b = 1.9801...
74 Computation steps:
75 0: allow forces to update the context state
76 1: v <- v + 0.5*dt*f/m
77 2: x <- x + 0.5*dt*c*tanh(v/c); c=sqrt(1*kT/m)
78 3: v1 <- v1 + 0.5*dt*(m*v*c*tanh(v/c) - kT)*invQ; c=sqrt(1*kT/m)
79 4: v <- v*exp(-0.5*dt*v1)
80 5: v1 <- a*v1 + b*gaussian
81 6: v <- v*exp(-0.5*dt*v1)
82 7: v1 <- v1 + 0.5*dt*(m*v*c*tanh(v/c) - kT)*invQ; c=sqrt(1*kT/m)
83 8: x <- x + 0.5*dt*c*tanh(v/c); c=sqrt(1*kT/m)
84 9: v <- v + 0.5*dt*f/m
86 Force-first scheme
88 >>> integrator_ff = xps.integrators.RegulatedNHLIntegrator(
89 ... temperature=300 * unit.kelvin,
90 ... timeConstant=100 * unit.femtoseconds,
91 ... stepSize=2 * unit.femtoseconds,
92 ... forceFirst=True
93 ... )
94 >>> model = testsystems.AlanineDipeptideVacuum()
95 >>> try:
96 ... integrator_ff.registerWithSystem(model.system, False)
97 ... except ValueError as e:
98 ... print(e)
99 Regulated NHL integrators do not support constraints.
100 """
102 @preprocess_args
103 def __init__( # noqa: PLR0913
104 self,
105 temperature: mmunit.Quantity,
106 timeConstant: mmunit.Quantity,
107 stepSize: mmunit.Quantity,
108 *,
109 regulationParameter: float = 1,
110 forceFirst: bool = False,
111 ) -> None:
112 super().__init__(stepSize)
113 self._forceFirst = forceFirst
114 self._regulation_parameter = regulationParameter
115 self._init_temperature(temperature)
116 self._init_time_constant(timeConstant)
117 self._init_friction_coefficient(1 / timeConstant)
118 self._add_variables()
119 self.addUpdateContextState()
120 self.setKineticEnergyExpression(
121 f"0.5*m*v*c*tanh(v/c); c=sqrt({self._regulation_parameter}*kT/m)"
122 )
123 self._add_boost(1 if forceFirst else 0.5)
124 self._add_translation(0.5)
125 self._add_thermostat(1)
126 self._add_translation(0.5)
127 if not forceFirst:
128 self._add_boost(0.5)
130 def _add_variables(self) -> None:
131 self.addGlobalVariable("kT", 0)
132 self.addGlobalVariable("invQ", 0)
133 self.addGlobalVariable("a", 0)
134 self.addGlobalVariable("b", 0)
135 self.addPerDofVariable("v1", 0)
136 self._update_global_variables()
138 def _update_global_variables(self) -> None:
139 tau = self.getTimeConstant()
140 dt = self.getStepSize()
141 friction = self.getFrictionCoefficient()
142 kt = mmunit.MOLAR_GAS_CONSTANT_R * self.getTemperature()
143 self.setGlobalVariableByName("kT", kt)
144 self.setGlobalVariableByName("invQ", 1 / (kt * tau**2))
145 self.setGlobalVariableByName("a", np.exp(-friction * dt))
146 self.setGlobalVariableByName("b", np.sqrt(1 - np.exp(-2 * friction * dt)) / tau)
148 def _add_translation(self, fraction: float) -> None:
149 self.addComputePerDof(
150 "x",
151 f"x + {fraction}*dt*c*tanh(v/c); c=sqrt({self._regulation_parameter}*kT/m)",
152 )
154 def _add_boost(self, fraction: float) -> None:
155 self.addComputePerDof("v", f"v + {fraction}*dt*f/m")
157 def _add_v1_boost(self, fraction: float) -> str:
158 self.addComputePerDof(
159 "v1",
160 f"v1 + {fraction}*dt*(m*v*c*tanh(v/c) - kT)*invQ"
161 f"; c=sqrt({self._regulation_parameter}*kT/m)",
162 )
164 def _add_v_scaling(self, fraction: float) -> str:
165 return self.addComputePerDof("v", f"v*exp(-{fraction}*dt*v1)")
167 def _add_thermostat(self, fraction: float) -> None:
168 self._add_v1_boost(0.5 * fraction)
169 self._add_v_scaling(0.5 * fraction)
170 self.addComputePerDof("v1", "a*v1 + b*gaussian")
171 self._add_v_scaling(0.5 * fraction)
172 self._add_v1_boost(0.5 * fraction)
174 def registerWithSystem(self, system: mm.System, isExtension: bool) -> None:
175 if system.getNumConstraints() > 0:
176 raise ValueError("Regulated NHL integrators do not support constraints.")