Coverage for ase / md / langevin.py: 83.33%
78 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 08:22 +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 'communicator' in kwargs:
86 msg = (
87 '`communicator` has been deprecated since ASE 3.25.0 '
88 'and will be removed in ASE 3.26.0. Use `comm` instead.'
89 )
90 warnings.warn(msg, FutureWarning)
91 kwargs['comm'] = kwargs.pop('communicator')
93 if friction is None:
94 raise TypeError("Missing 'friction' argument.")
95 self.fr = friction
96 self.temp = units.kB * self._process_temperature(
97 temperature, temperature_K, 'eV'
98 )
100 if fixcm:
101 msg = (
102 'The implementation of `fixcm=True` in `Langevin` does not '
103 'strictly sample the correct NVT distributions. '
104 'The deviations are typically small for large systems but can '
105 'be more pronounced for small systems. '
106 'Use `fixcm=False` together with `ase.constraints.FixCom`. '
107 '`fixcm` is deprecated since ASE 3.28.0 and will be removed in '
108 'a future release.'
109 )
110 warnings.warn(msg, FutureWarning)
111 self.fix_com = fixcm
113 if rng is None:
114 self.rng = np.random
115 else:
116 self.rng = rng
117 MolecularDynamics.__init__(self, atoms, timestep, **kwargs)
118 self.updatevars()
120 def todict(self):
121 d = MolecularDynamics.todict(self)
122 d.update(
123 {
124 'temperature_K': self.temp / units.kB,
125 'friction': self.fr,
126 'fixcm': self.fix_com,
127 }
128 )
129 return d
131 def set_temperature(self, temperature=None, temperature_K=None):
132 self.temp = units.kB * self._process_temperature(
133 temperature, temperature_K, 'eV'
134 )
135 self.updatevars()
137 def set_friction(self, friction):
138 self.fr = friction
139 self.updatevars()
141 def set_timestep(self, timestep):
142 self.dt = timestep
143 self.updatevars()
145 def updatevars(self):
146 dt = self.dt
147 T = self.temp
148 fr = self.fr
149 masses = self.masses
150 sigma = np.sqrt(2 * T * fr / masses)
152 self.c1 = dt / 2.0 - dt * dt * fr / 8.0
153 self.c2 = dt * fr / 2 - dt * dt * fr * fr / 8.0
154 self.c3 = np.sqrt(dt) * sigma / 2.0 - dt**1.5 * fr * sigma / 8.0
155 self.c5 = dt**1.5 * sigma / (2 * np.sqrt(3))
156 self.c4 = fr / 2.0 * self.c5
158 def step(self, forces=None):
159 atoms = self.atoms
160 natoms = len(atoms)
162 if forces is None:
163 forces = atoms.get_forces(md=True)
165 # This velocity as well as rnd_pos, rnd_mom and a few other
166 # variables are stored as attributes, so Asap can do its magic
167 # when atoms migrate between processors.
168 self.v = atoms.get_velocities()
170 xi = self.rng.standard_normal(size=(natoms, 3))
171 eta = self.rng.standard_normal(size=(natoms, 3))
173 # When holonomic constraints for rigid linear triatomic molecules are
174 # present, ask the constraints to redistribute xi and eta within each
175 # triple defined in the constraints. This is needed to achieve the
176 # correct target temperature.
177 for constraint in self.atoms.constraints:
178 if hasattr(constraint, 'redistribute_forces_md'):
179 constraint.redistribute_forces_md(atoms, xi, rand=True)
180 constraint.redistribute_forces_md(atoms, eta, rand=True)
182 self.comm.broadcast(xi, 0)
183 self.comm.broadcast(eta, 0)
185 # To keep the center of mass stationary, we have to calculate
186 # the random perturbations to the positions and the momenta,
187 # and make sure that they sum to zero. This perturbs the
188 # temperature slightly, and we have to correct.
189 self.rnd_pos = self.c5 * eta
190 self.rnd_vel = self.c3 * xi - self.c4 * eta
191 # https://gitlab.com/ase/ase/-/merge_requests/3986
192 # The implementation of `fixcm=True` does not strictly sample the
193 # correct NVT distributions for positions and momenta.
194 # It is deprecated and will be removed in a future ASE release.
195 if self.fix_com:
196 factor = np.sqrt(natoms / (natoms - 1.0))
197 self.rnd_pos -= self.rnd_pos.sum(axis=0) / natoms
198 self.rnd_vel -= (self.rnd_vel * self.masses).sum(axis=0) / (
199 self.masses * natoms
200 )
201 self.rnd_pos *= factor
202 self.rnd_vel *= factor
204 # First halfstep in the velocity.
205 self.v += (
206 self.c1 * forces / self.masses - self.c2 * self.v + self.rnd_vel
207 )
209 # Full step in positions
210 x = atoms.get_positions()
212 # Step: x^n -> x^(n+1) - this applies constraints if any.
213 atoms.set_positions(x + self.dt * self.v + self.rnd_pos)
215 # recalc velocities after RATTLE constraints are applied
216 self.v = (self.atoms.get_positions() - x - self.rnd_pos) / self.dt
217 forces = atoms.get_forces(md=True)
219 # Update the velocities
220 self.v += (
221 self.c1 * forces / self.masses - self.c2 * self.v + self.rnd_vel
222 )
224 # Second part of RATTLE taken care of here
225 atoms.set_momenta(self.v * self.masses)
227 return forces