Coverage for ase / md / langevin.py: 86.49%
74 statements
« prev ^ index » next coverage.py v7.14.0, created at 2026-05-21 15:52 +0000
1"""Langevin dynamics class."""
3import warnings
5import numpy as np
7from ase import Atoms, units
8from ase.md.md import MolecularDynamics
11class Langevin(MolecularDynamics):
12 """Langevin (constant N, V, T) molecular dynamics."""
14 # Helps Asap doing the right thing. Increment when changing stuff:
15 _lgv_version = 5
17 def __init__(
18 self,
19 atoms: Atoms,
20 timestep: float,
21 temperature: float | None = None,
22 friction: float | None = None,
23 fixcm: bool = True,
24 *,
25 temperature_K: float | None = None,
26 rng=None,
27 **kwargs,
28 ):
29 """
30 Parameters
31 ----------
32 atoms: Atoms object
33 The list of atoms.
35 timestep: float
36 The time step in ASE time units.
38 temperature: float (deprecated)
39 The desired temperature, in electron volt.
41 temperature_K: float
42 The desired temperature, in Kelvin.
44 friction: float
45 A friction coefficient in inverse ASE time units.
46 For example, set ``0.01 / ase.units.fs`` to provide
47 0.01 fs\\ :sup:`−1` (10 ps\\ :sup:`−1`).
49 fixcm: bool (optional)
50 If True, the position and momentum of the center of mass is
51 kept unperturbed. Default: True.
53 .. deprecated:: 3.28.0
55 The implementation of ``fixcm=True`` does not strictly sample
56 the correct NVT distributions for positions and momenta.
57 The deviations are typically small for large systems but can be
58 more pronounced for small systems.
60 Use :class:`~ase.constraints.FixCom` instead.
61 ``fixcm`` will be removed from ASE in the near future.
63 rng: RNG object (optional)
64 Random number generator, by default numpy.random. Must have a
65 standard_normal method matching the signature of
66 numpy.random.standard_normal.
68 **kwargs : dict, optional
69 Additional arguments passed to
70 :class:`~ase.md.md.MolecularDynamics` base class.
72 Notes
73 -----
74 The temperature and friction are normally scalars, but in principle one
75 quantity per atom could be specified by giving an array.
77 RATTLE constraints can be used with these propagators, see:
78 E. V.-Eijnden, and G. Ciccotti, Chem. Phys. Lett. 429, 310 (2006)
80 The propagator is Equation 23 (Eq. 39 if RATTLE constraints are used)
81 of the above reference. That reference also contains another
82 propagator in Eq. 21/34; but that propagator is not quasi-symplectic
83 and gives a systematic offset in the temperature at large time steps.
84 """
85 if friction is None:
86 raise TypeError("Missing 'friction' argument.")
87 self.fr = friction
88 self.temp = units.kB * self._process_temperature(
89 temperature, temperature_K, 'eV'
90 )
92 if fixcm:
93 msg = (
94 'The implementation of `fixcm=True` in `Langevin` does not '
95 'strictly sample the correct NVT distributions. '
96 'The deviations are typically small for large systems but can '
97 'be more pronounced for small systems. '
98 'Use `fixcm=False` together with `ase.constraints.FixCom`. '
99 '`fixcm` is deprecated since ASE 3.28.0 and will be removed in '
100 'a future release.'
101 )
102 warnings.warn(msg, FutureWarning)
103 self.fix_com = fixcm
105 if rng is None:
106 self.rng = np.random
107 else:
108 self.rng = rng
109 MolecularDynamics.__init__(self, atoms, timestep, **kwargs)
110 self.updatevars()
112 def todict(self):
113 d = MolecularDynamics.todict(self)
114 d.update(
115 {
116 'temperature_K': self.temp / units.kB,
117 'friction': self.fr,
118 'fixcm': self.fix_com,
119 }
120 )
121 return d
123 def set_temperature(self, temperature=None, temperature_K=None):
124 self.temp = units.kB * self._process_temperature(
125 temperature, temperature_K, 'eV'
126 )
127 self.updatevars()
129 def set_friction(self, friction):
130 self.fr = friction
131 self.updatevars()
133 def set_timestep(self, timestep):
134 self.dt = timestep
135 self.updatevars()
137 def updatevars(self):
138 dt = self.dt
139 T = self.temp
140 fr = self.fr
141 masses = self.masses
142 sigma = np.sqrt(2 * T * fr / masses)
144 self.c1 = dt / 2.0 - dt * dt * fr / 8.0
145 self.c2 = dt * fr / 2 - dt * dt * fr * fr / 8.0
146 self.c3 = np.sqrt(dt) * sigma / 2.0 - dt**1.5 * fr * sigma / 8.0
147 self.c5 = dt**1.5 * sigma / (2 * np.sqrt(3))
148 self.c4 = fr / 2.0 * self.c5
150 def step(self, forces=None):
151 atoms = self.atoms
152 natoms = len(atoms)
154 if forces is None:
155 forces = atoms.get_forces(md=True)
157 # This velocity as well as rnd_pos, rnd_mom and a few other
158 # variables are stored as attributes, so Asap can do its magic
159 # when atoms migrate between processors.
160 self.v = atoms.get_velocities()
162 xi = self.rng.standard_normal(size=(natoms, 3))
163 eta = self.rng.standard_normal(size=(natoms, 3))
165 # When holonomic constraints for rigid linear triatomic molecules are
166 # present, ask the constraints to redistribute xi and eta within each
167 # triple defined in the constraints. This is needed to achieve the
168 # correct target temperature.
169 for constraint in self.atoms.constraints:
170 if hasattr(constraint, 'redistribute_forces_md'):
171 constraint.redistribute_forces_md(atoms, xi, rand=True)
172 constraint.redistribute_forces_md(atoms, eta, rand=True)
174 self.comm.broadcast(xi, 0)
175 self.comm.broadcast(eta, 0)
177 # To keep the center of mass stationary, we have to calculate
178 # the random perturbations to the positions and the momenta,
179 # and make sure that they sum to zero. This perturbs the
180 # temperature slightly, and we have to correct.
181 self.rnd_pos = self.c5 * eta
182 self.rnd_vel = self.c3 * xi - self.c4 * eta
183 # https://gitlab.com/ase/ase/-/merge_requests/3986
184 # The implementation of `fixcm=True` does not strictly sample the
185 # correct NVT distributions for positions and momenta.
186 # It is deprecated and will be removed in a future ASE release.
187 if self.fix_com:
188 factor = np.sqrt(natoms / (natoms - 1.0))
189 self.rnd_pos -= self.rnd_pos.sum(axis=0) / natoms
190 self.rnd_vel -= (self.rnd_vel * self.masses).sum(axis=0) / (
191 self.masses * natoms
192 )
193 self.rnd_pos *= factor
194 self.rnd_vel *= factor
196 # First halfstep in the velocity.
197 self.v += (
198 self.c1 * forces / self.masses - self.c2 * self.v + self.rnd_vel
199 )
201 # Full step in positions
202 x = atoms.get_positions()
204 # Step: x^n -> x^(n+1) - this applies constraints if any.
205 atoms.set_positions(x + self.dt * self.v + self.rnd_pos)
207 # recalc velocities after RATTLE constraints are applied
208 self.v = (self.atoms.get_positions() - x - self.rnd_pos) / self.dt
209 forces = atoms.get_forces(md=True)
211 # Update the velocities
212 self.v += (
213 self.c1 * forces / self.masses - self.c2 * self.v + self.rnd_vel
214 )
216 # Second part of RATTLE taken care of here
217 atoms.set_momenta(self.v * self.masses)
219 return forces