Hide keyboard shortcuts

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 

2 

3from ase.md.md import MolecularDynamics 

4import warnings 

5 

6 

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. 

11 

12 Parameters: 

13 

14 atoms: Atoms object 

15 The Atoms object to operate on. 

16 

17 timestep: float 

18 The time step in ASE time units. 

19 

20 trajectory: Trajectory object or str (optional) 

21 Attach trajectory object. If *trajectory* is a string a 

22 Trajectory will be constructed. Default: None. 

23 

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. 

27 

28 loginterval: int (optional) 

29 Only write a log line for every *loginterval* time steps. 

30 Default: 1 

31 

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. 

37 

38 dt: float (deprecated) 

39 Alias for timestep. 

40 """ 

41 if dt is not None: 

42 warnings.warn( 

43 FutureWarning( 

44 'dt variable is deprecated; please use timestep.')) 

45 timestep = dt 

46 if timestep is None: 

47 raise TypeError('Missing timestep argument') 

48 

49 MolecularDynamics.__init__(self, atoms, timestep, trajectory, logfile, 

50 loginterval, 

51 append_trajectory=append_trajectory) 

52 

53 def step(self, forces=None): 

54 

55 atoms = self.atoms 

56 

57 if forces is None: 

58 forces = atoms.get_forces(md=True) 

59 

60 p = atoms.get_momenta() 

61 p += 0.5 * self.dt * forces 

62 masses = atoms.get_masses()[:, np.newaxis] 

63 r = atoms.get_positions() 

64 

65 # if we have constraints then this will do the first part of the 

66 # RATTLE algorithm: 

67 atoms.set_positions(r + self.dt * p / masses) 

68 if atoms.constraints: 

69 p = (atoms.get_positions() - r) * masses / self.dt 

70 

71 # We need to store the momenta on the atoms before calculating 

72 # the forces, as in a parallel Asap calculation atoms may 

73 # migrate during force calculations, and the momenta need to 

74 # migrate along with the atoms. 

75 atoms.set_momenta(p, apply_constraint=False) 

76 

77 forces = atoms.get_forces(md=True) 

78 

79 # Second part of RATTLE will be done here: 

80 atoms.set_momenta(atoms.get_momenta() + 0.5 * self.dt * forces) 

81 return forces