Coverage for /builds/debichem-team/python-ase/ase/md/langevin.py : 94.44%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""Langevin dynamics class."""
3import numpy as np
5from ase.md.md import MolecularDynamics
6from ase.parallel import world, DummyMPI
7from ase import units
10class Langevin(MolecularDynamics):
11 """Langevin (constant N, V, T) molecular dynamics."""
13 # Helps Asap doing the right thing. Increment when changing stuff:
14 _lgv_version = 4
16 def __init__(self, atoms, timestep, temperature=None, friction=None,
17 fixcm=True, *, temperature_K=None, trajectory=None,
18 logfile=None, loginterval=1, communicator=world,
19 rng=None, append_trajectory=False):
20 """
21 Parameters:
23 atoms: Atoms object
24 The list of atoms.
26 timestep: float
27 The time step in ASE time units.
29 temperature: float (deprecated)
30 The desired temperature, in electron volt.
32 temperature_K: float
33 The desired temperature, in Kelvin.
35 friction: float
36 A friction coefficient, typically 1e-4 to 1e-2.
38 fixcm: bool (optional)
39 If True, the position and momentum of the center of mass is
40 kept unperturbed. Default: True.
42 rng: RNG object (optional)
43 Random number generator, by default numpy.random. Must have a
44 standard_normal method matching the signature of
45 numpy.random.standard_normal.
47 logfile: file object or str (optional)
48 If *logfile* is a string, a file with that name will be opened.
49 Use '-' for stdout.
51 trajectory: Trajectory object or str (optional)
52 Attach trajectory object. If *trajectory* is a string a
53 Trajectory will be constructed. Use *None* (the default) for no
54 trajectory.
56 communicator: MPI communicator (optional)
57 Communicator used to distribute random numbers to all tasks.
58 Default: ase.parallel.world. Set to None to disable communication.
60 append_trajectory: bool (optional)
61 Defaults to False, which causes the trajectory file to be
62 overwritten each time the dynamics is restarted from scratch.
63 If True, the new structures are appended to the trajectory
64 file instead.
66 The temperature and friction are normally scalars, but in principle one
67 quantity per atom could be specified by giving an array.
69 RATTLE constraints can be used with these propagators, see:
70 E. V.-Eijnden, and G. Ciccotti, Chem. Phys. Lett. 429, 310 (2006)
72 The propagator is Equation 23 (Eq. 39 if RATTLE constraints are used)
73 of the above reference. That reference also contains another
74 propagator in Eq. 21/34; but that propagator is not quasi-symplectic
75 and gives a systematic offset in the temperature at large time steps.
76 """
77 if friction is None:
78 raise TypeError("Missing 'friction' argument.")
79 self.fr = friction
80 self.temp = units.kB * self._process_temperature(temperature,
81 temperature_K, 'eV')
82 self.fix_com = fixcm
83 if communicator is None:
84 communicator = DummyMPI()
85 self.communicator = communicator
86 if rng is None:
87 self.rng = np.random
88 else:
89 self.rng = rng
90 MolecularDynamics.__init__(self, atoms, timestep, trajectory,
91 logfile, loginterval,
92 append_trajectory=append_trajectory)
93 self.updatevars()
95 def todict(self):
96 d = MolecularDynamics.todict(self)
97 d.update({'temperature_K': self.temp / units.kB,
98 'friction': self.fr,
99 'fixcm': self.fix_com})
100 return d
102 def set_temperature(self, temperature=None, temperature_K=None):
103 self.temp = units.kB * self._process_temperature(temperature,
104 temperature_K, 'eV')
105 self.updatevars()
107 def set_friction(self, friction):
108 self.fr = friction
109 self.updatevars()
111 def set_timestep(self, timestep):
112 self.dt = timestep
113 self.updatevars()
115 def updatevars(self):
116 dt = self.dt
117 T = self.temp
118 fr = self.fr
119 masses = self.masses
120 sigma = np.sqrt(2 * T * fr / masses)
122 self.c1 = dt / 2. - dt * dt * fr / 8.
123 self.c2 = dt * fr / 2 - dt * dt * fr * fr / 8.
124 self.c3 = np.sqrt(dt) * sigma / 2. - dt**1.5 * fr * sigma / 8.
125 self.c5 = dt**1.5 * sigma / (2 * np.sqrt(3))
126 self.c4 = fr / 2. * self.c5
128 def step(self, forces=None):
129 atoms = self.atoms
130 natoms = len(atoms)
132 if forces is None:
133 forces = atoms.get_forces(md=True)
135 # This velocity as well as xi, eta and a few other variables are stored
136 # as attributes, so Asap can do its magic when atoms migrate between
137 # processors.
138 self.v = atoms.get_velocities()
140 self.xi = self.rng.standard_normal(size=(natoms, 3))
141 self.eta = self.rng.standard_normal(size=(natoms, 3))
143 # When holonomic constraints for rigid linear triatomic molecules are
144 # present, ask the constraints to redistribute xi and eta within each
145 # triple defined in the constraints. This is needed to achieve the
146 # correct target temperature.
147 for constraint in self.atoms.constraints:
148 if hasattr(constraint, 'redistribute_forces_md'):
149 constraint.redistribute_forces_md(atoms, self.xi, rand=True)
150 constraint.redistribute_forces_md(atoms, self.eta, rand=True)
152 self.communicator.broadcast(self.xi, 0)
153 self.communicator.broadcast(self.eta, 0)
155 # First halfstep in the velocity.
156 self.v += (self.c1 * forces / self.masses - self.c2 * self.v +
157 self.c3 * self.xi - self.c4 * self.eta)
159 # Full step in positions
160 x = atoms.get_positions()
161 if self.fix_com:
162 old_com = atoms.get_center_of_mass()
163 # Step: x^n -> x^(n+1) - this applies constraints if any.
164 atoms.set_positions(x + self.dt * self.v + self.c5 * self.eta)
165 if self.fix_com:
166 atoms.set_center_of_mass(old_com)
168 # recalc velocities after RATTLE constraints are applied
169 self.v = (self.atoms.get_positions() - x -
170 self.c5 * self.eta) / self.dt
171 forces = atoms.get_forces(md=True)
173 # Update the velocities
174 self.v += (self.c1 * forces / self.masses - self.c2 * self.v +
175 self.c3 * self.xi - self.c4 * self.eta)
177 if self.fix_com: # subtract center of mass vel
178 self.v -= self._get_com_velocity(self.v)
180 # Second part of RATTLE taken care of here
181 atoms.set_momenta(self.v * self.masses)
183 return forces