Coverage for /builds/kinetik161/ase/ase/md/langevin.py: 94.44%
72 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1"""Langevin dynamics class."""
2from typing import IO, Optional, Union
4import numpy as np
6from ase import Atoms, units
7from ase.md.md import MolecularDynamics
8from ase.parallel import DummyMPI, world
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: Optional[float] = None,
22 friction: Optional[float] = None,
23 fixcm: bool = True,
24 *,
25 temperature_K: Optional[float] = None,
26 trajectory: Optional[str] = None,
27 logfile: Optional[Union[IO, str]] = None,
28 loginterval: int = 1,
29 communicator=world,
30 rng=None,
31 append_trajectory: bool = False,
32 ):
33 """
34 Parameters:
36 atoms: Atoms object
37 The list of atoms.
39 timestep: float
40 The time step in ASE time units.
42 temperature: float (deprecated)
43 The desired temperature, in electron volt.
45 temperature_K: float
46 The desired temperature, in Kelvin.
48 friction: float
49 A friction coefficient in inverse ASE time units.
50 For example, set ``0.01 / ase.units.fs`` to provide
51 0.01 fs\\ :sup:`−1` (10 ps\\ :sup:`−1`).
53 fixcm: bool (optional)
54 If True, the position and momentum of the center of mass is
55 kept unperturbed. Default: True.
57 rng: RNG object (optional)
58 Random number generator, by default numpy.random. Must have a
59 standard_normal method matching the signature of
60 numpy.random.standard_normal.
62 logfile: file object or str (optional)
63 If *logfile* is a string, a file with that name will be opened.
64 Use '-' for stdout.
66 trajectory: Trajectory object or str (optional)
67 Attach trajectory object. If *trajectory* is a string a
68 Trajectory will be constructed. Use *None* (the default) for no
69 trajectory.
71 communicator: MPI communicator (optional)
72 Communicator used to distribute random numbers to all tasks.
73 Default: ase.parallel.world. Set to None to disable communication.
75 append_trajectory: bool (optional)
76 Defaults to False, which causes the trajectory file to be
77 overwritten each time the dynamics is restarted from scratch.
78 If True, the new structures are appended to the trajectory
79 file instead.
81 The temperature and friction are normally scalars, but in principle one
82 quantity per atom could be specified by giving an array.
84 RATTLE constraints can be used with these propagators, see:
85 E. V.-Eijnden, and G. Ciccotti, Chem. Phys. Lett. 429, 310 (2006)
87 The propagator is Equation 23 (Eq. 39 if RATTLE constraints are used)
88 of the above reference. That reference also contains another
89 propagator in Eq. 21/34; but that propagator is not quasi-symplectic
90 and gives a systematic offset in the temperature at large time steps.
91 """
92 if friction is None:
93 raise TypeError("Missing 'friction' argument.")
94 self.fr = friction
95 self.temp = units.kB * self._process_temperature(temperature,
96 temperature_K, 'eV')
97 self.fix_com = fixcm
98 if communicator is None:
99 communicator = DummyMPI()
100 self.communicator = communicator
101 if rng is None:
102 self.rng = np.random
103 else:
104 self.rng = rng
105 MolecularDynamics.__init__(self, atoms, timestep, trajectory,
106 logfile, loginterval,
107 append_trajectory=append_trajectory)
108 self.updatevars()
110 def todict(self):
111 d = MolecularDynamics.todict(self)
112 d.update({'temperature_K': self.temp / units.kB,
113 'friction': self.fr,
114 'fixcm': self.fix_com})
115 return d
117 def set_temperature(self, temperature=None, temperature_K=None):
118 self.temp = units.kB * self._process_temperature(temperature,
119 temperature_K, 'eV')
120 self.updatevars()
122 def set_friction(self, friction):
123 self.fr = friction
124 self.updatevars()
126 def set_timestep(self, timestep):
127 self.dt = timestep
128 self.updatevars()
130 def updatevars(self):
131 dt = self.dt
132 T = self.temp
133 fr = self.fr
134 masses = self.masses
135 sigma = np.sqrt(2 * T * fr / masses)
137 self.c1 = dt / 2. - dt * dt * fr / 8.
138 self.c2 = dt * fr / 2 - dt * dt * fr * fr / 8.
139 self.c3 = np.sqrt(dt) * sigma / 2. - dt**1.5 * fr * sigma / 8.
140 self.c5 = dt**1.5 * sigma / (2 * np.sqrt(3))
141 self.c4 = fr / 2. * self.c5
143 def step(self, forces=None):
144 atoms = self.atoms
145 natoms = len(atoms)
147 if forces is None:
148 forces = atoms.get_forces(md=True)
150 # This velocity as well as rnd_pos, rnd_mom and a few other
151 # variables are stored as attributes, so Asap can do its magic
152 # when atoms migrate between processors.
153 self.v = atoms.get_velocities()
155 xi = self.rng.standard_normal(size=(natoms, 3))
156 eta = self.rng.standard_normal(size=(natoms, 3))
158 # When holonomic constraints for rigid linear triatomic molecules are
159 # present, ask the constraints to redistribute xi and eta within each
160 # triple defined in the constraints. This is needed to achieve the
161 # correct target temperature.
162 for constraint in self.atoms.constraints:
163 if hasattr(constraint, 'redistribute_forces_md'):
164 constraint.redistribute_forces_md(atoms, xi, rand=True)
165 constraint.redistribute_forces_md(atoms, eta, rand=True)
167 self.communicator.broadcast(xi, 0)
168 self.communicator.broadcast(eta, 0)
170 # To keep the center of mass stationary, we have to calculate
171 # the random perturbations to the positions and the momenta,
172 # and make sure that they sum to zero.
173 self.rnd_pos = self.c5 * eta
174 self.rnd_vel = self.c3 * xi - self.c4 * eta
175 if self.fix_com:
176 self.rnd_pos -= self.rnd_pos.sum(axis=0) / natoms
177 self.rnd_vel -= (self.rnd_vel *
178 self.masses).sum(axis=0) / (self.masses * natoms)
180 # First halfstep in the velocity.
181 self.v += (self.c1 * forces / self.masses - self.c2 * self.v +
182 self.rnd_vel)
184 # Full step in positions
185 x = atoms.get_positions()
187 # Step: x^n -> x^(n+1) - this applies constraints if any.
188 atoms.set_positions(x + self.dt * self.v + self.rnd_pos)
190 # recalc velocities after RATTLE constraints are applied
191 self.v = (self.atoms.get_positions() - x - self.rnd_pos) / self.dt
192 forces = atoms.get_forces(md=True)
194 # Update the velocities
195 self.v += (self.c1 * forces / self.masses - self.c2 * self.v +
196 self.rnd_vel)
198 # Second part of RATTLE taken care of here
199 atoms.set_momenta(self.v * self.masses)
201 return forces