Coverage for /builds/ase/ase/ase/md/langevin.py: 92.11%
76 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1"""Langevin dynamics class."""
3import warnings
4from typing import Optional
6import numpy as np
8from ase import Atoms, units
9from ase.md.md import MolecularDynamics
12class Langevin(MolecularDynamics):
13 """Langevin (constant N, V, T) molecular dynamics."""
15 # Helps Asap doing the right thing. Increment when changing stuff:
16 _lgv_version = 5
18 def __init__(
19 self,
20 atoms: Atoms,
21 timestep: float,
22 temperature: Optional[float] = None,
23 friction: Optional[float] = None,
24 fixcm: bool = True,
25 *,
26 temperature_K: Optional[float] = None,
27 rng=None,
28 **kwargs,
29 ):
30 """
31 Parameters
32 ----------
33 atoms: Atoms object
34 The list of atoms.
36 timestep: float
37 The time step in ASE time units.
39 temperature: float (deprecated)
40 The desired temperature, in electron volt.
42 temperature_K: float
43 The desired temperature, in Kelvin.
45 friction: float
46 A friction coefficient in inverse ASE time units.
47 For example, set ``0.01 / ase.units.fs`` to provide
48 0.01 fs\\ :sup:`−1` (10 ps\\ :sup:`−1`).
50 fixcm: bool (optional)
51 If True, the position and momentum of the center of mass is
52 kept unperturbed. Default: True.
54 rng: RNG object (optional)
55 Random number generator, by default numpy.random. Must have a
56 standard_normal method matching the signature of
57 numpy.random.standard_normal.
59 **kwargs : dict, optional
60 Additional arguments passed to :class:~ase.md.md.MolecularDynamics
61 base class.
63 The temperature and friction are normally scalars, but in principle one
64 quantity per atom could be specified by giving an array.
66 RATTLE constraints can be used with these propagators, see:
67 E. V.-Eijnden, and G. Ciccotti, Chem. Phys. Lett. 429, 310 (2006)
69 The propagator is Equation 23 (Eq. 39 if RATTLE constraints are used)
70 of the above reference. That reference also contains another
71 propagator in Eq. 21/34; but that propagator is not quasi-symplectic
72 and gives a systematic offset in the temperature at large time steps.
73 """
74 if 'communicator' in kwargs:
75 msg = (
76 '`communicator` has been deprecated since ASE 3.25.0 '
77 'and will be removed in ASE 3.26.0. Use `comm` instead.'
78 )
79 warnings.warn(msg, FutureWarning)
80 kwargs['comm'] = kwargs.pop('communicator')
82 if friction is None:
83 raise TypeError("Missing 'friction' argument.")
84 self.fr = friction
85 self.temp = units.kB * self._process_temperature(
86 temperature, temperature_K, 'eV'
87 )
88 self.fix_com = fixcm
90 if rng is None:
91 self.rng = np.random
92 else:
93 self.rng = rng
94 MolecularDynamics.__init__(self, atoms, timestep, **kwargs)
95 self.updatevars()
97 def todict(self):
98 d = MolecularDynamics.todict(self)
99 d.update(
100 {
101 'temperature_K': self.temp / units.kB,
102 'friction': self.fr,
103 'fixcm': self.fix_com,
104 }
105 )
106 return d
108 def set_temperature(self, temperature=None, temperature_K=None):
109 self.temp = units.kB * self._process_temperature(
110 temperature, temperature_K, 'eV'
111 )
112 self.updatevars()
114 def set_friction(self, friction):
115 self.fr = friction
116 self.updatevars()
118 def set_timestep(self, timestep):
119 self.dt = timestep
120 self.updatevars()
122 def updatevars(self):
123 dt = self.dt
124 T = self.temp
125 fr = self.fr
126 masses = self.masses
127 sigma = np.sqrt(2 * T * fr / masses)
129 self.c1 = dt / 2.0 - dt * dt * fr / 8.0
130 self.c2 = dt * fr / 2 - dt * dt * fr * fr / 8.0
131 self.c3 = np.sqrt(dt) * sigma / 2.0 - dt**1.5 * fr * sigma / 8.0
132 self.c5 = dt**1.5 * sigma / (2 * np.sqrt(3))
133 self.c4 = fr / 2.0 * self.c5
135 def step(self, forces=None):
136 atoms = self.atoms
137 natoms = len(atoms)
139 if forces is None:
140 forces = atoms.get_forces(md=True)
142 # This velocity as well as rnd_pos, rnd_mom and a few other
143 # variables are stored as attributes, so Asap can do its magic
144 # when atoms migrate between processors.
145 self.v = atoms.get_velocities()
147 xi = self.rng.standard_normal(size=(natoms, 3))
148 eta = self.rng.standard_normal(size=(natoms, 3))
150 # When holonomic constraints for rigid linear triatomic molecules are
151 # present, ask the constraints to redistribute xi and eta within each
152 # triple defined in the constraints. This is needed to achieve the
153 # correct target temperature.
154 for constraint in self.atoms.constraints:
155 if hasattr(constraint, 'redistribute_forces_md'):
156 constraint.redistribute_forces_md(atoms, xi, rand=True)
157 constraint.redistribute_forces_md(atoms, eta, rand=True)
159 self.comm.broadcast(xi, 0)
160 self.comm.broadcast(eta, 0)
162 # To keep the center of mass stationary, we have to calculate
163 # the random perturbations to the positions and the momenta,
164 # and make sure that they sum to zero. This perturbs the
165 # temperature slightly, and we have to correct.
166 self.rnd_pos = self.c5 * eta
167 self.rnd_vel = self.c3 * xi - self.c4 * eta
168 if self.fix_com:
169 factor = np.sqrt(natoms / (natoms - 1.0))
170 self.rnd_pos -= self.rnd_pos.sum(axis=0) / natoms
171 self.rnd_vel -= (self.rnd_vel * self.masses).sum(axis=0) / (
172 self.masses * natoms
173 )
174 self.rnd_pos *= factor
175 self.rnd_vel *= factor
177 # First halfstep in the velocity.
178 self.v += (
179 self.c1 * forces / self.masses - self.c2 * self.v + self.rnd_vel
180 )
182 # Full step in positions
183 x = atoms.get_positions()
185 # Step: x^n -> x^(n+1) - this applies constraints if any.
186 atoms.set_positions(x + self.dt * self.v + self.rnd_pos)
188 # recalc velocities after RATTLE constraints are applied
189 self.v = (self.atoms.get_positions() - x - self.rnd_pos) / self.dt
190 forces = atoms.get_forces(md=True)
192 # Update the velocities
193 self.v += (
194 self.c1 * forces / self.masses - self.c2 * self.v + self.rnd_vel
195 )
197 # Second part of RATTLE taken care of here
198 atoms.set_momenta(self.v * self.masses)
200 return forces