Coverage for /builds/debichem-team/python-ase/ase/md/verlet.py : 88.46%
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
1import numpy as np
3from ase.md.md import MolecularDynamics
4import warnings
7class VelocityVerlet(MolecularDynamics):
8 def __init__(self, atoms, timestep=None, trajectory=None, logfile=None,
9 loginterval=1, dt=None, append_trajectory=False):
10 """Molecular Dynamics object.
12 Parameters:
14 atoms: Atoms object
15 The Atoms object to operate on.
17 timestep: float
18 The time step in ASE time units.
20 trajectory: Trajectory object or str (optional)
21 Attach trajectory object. If *trajectory* is a string a
22 Trajectory will be constructed. Default: None.
24 logfile: file object or str (optional)
25 If *logfile* is a string, a file with that name will be opened.
26 Use '-' for stdout. Default: None.
28 loginterval: int (optional)
29 Only write a log line for every *loginterval* time steps.
30 Default: 1
32 append_trajectory: boolean
33 Defaults to False, which causes the trajectory file to be
34 overwriten each time the dynamics is restarted from scratch.
35 If True, the new structures are appended to the trajectory
36 file instead.
38 dt: float (deprecated)
39 Alias for timestep.
40 """
41 if dt is not None:
42 warnings.warn(FutureWarning('dt variable is deprecated; please use timestep.'))
43 timestep = dt
44 if timestep is None:
45 raise TypeError('Missing timestep argument')
47 MolecularDynamics.__init__(self, atoms, timestep, trajectory, logfile,
48 loginterval,
49 append_trajectory=append_trajectory)
51 def step(self, forces=None):
53 atoms = self.atoms
55 if forces is None:
56 forces = atoms.get_forces(md=True)
58 p = atoms.get_momenta()
59 p += 0.5 * self.dt * forces
60 masses = atoms.get_masses()[:, np.newaxis]
61 r = atoms.get_positions()
63 # if we have constraints then this will do the first part of the
64 # RATTLE algorithm:
65 atoms.set_positions(r + self.dt * p / masses)
66 if atoms.constraints:
67 p = (atoms.get_positions() - r) * masses / self.dt
69 # We need to store the momenta on the atoms before calculating
70 # the forces, as in a parallel Asap calculation atoms may
71 # migrate during force calculations, and the momenta need to
72 # migrate along with the atoms.
73 atoms.set_momenta(p, apply_constraint=False)
75 forces = atoms.get_forces(md=True)
77 # Second part of RATTLE will be done here:
78 atoms.set_momenta(atoms.get_momenta() + 0.5 * self.dt * forces)
79 return forces