Coverage for /builds/kinetik161/ase/ase/md/verlet.py: 89.29%

1import warnings 

2from typing import IO, Optional, Union 

3 

4import numpy as np 

5 

6from ase import Atoms 

7from ase.md.md import MolecularDynamics 

8 

9 

10class VelocityVerlet(MolecularDynamics): 

11 def __init__( 

12 self, 

13 atoms: Atoms, 

14 timestep: Optional[float] = None, 

15 trajectory: Optional[str] = None, 

16 logfile: Optional[Union[IO, str]] = None, 

17 loginterval: int = 1, 

18 dt: Optional[float] = None, 

19 append_trajectory: bool = False, 

20 ): 

21 """Molecular Dynamics object. 

22 

23 Parameters: 

24 

25 atoms: Atoms object 

26 The Atoms object to operate on. 

27 

28 timestep: float 

29 The time step in ASE time units. 

30 

31 trajectory: Trajectory object or str (optional) 

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

33 Trajectory will be constructed. Default: None. 

34 

35 logfile: file object or str (optional) 

36 If *logfile* is a string, a file with that name will be opened. 

37 Use '-' for stdout. Default: None. 

38 

39 loginterval: int (optional) 

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

41 Default: 1 

42 

43 append_trajectory: boolean 

44 Defaults to False, which causes the trajectory file to be 

45 overwriten each time the dynamics is restarted from scratch. 

46 If True, the new structures are appended to the trajectory 

47 file instead. 

48 

49 dt: float (deprecated) 

50 Alias for timestep. 

51 """ 

52 if dt is not None: 

53 warnings.warn( 

54 FutureWarning( 

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

56 timestep = dt 

57 if timestep is None: 

58 raise TypeError('Missing timestep argument') 

59 

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

61 loginterval, 

62 append_trajectory=append_trajectory) 

63 

64 def step(self, forces=None): 

65 

66 atoms = self.atoms 

67 

68 if forces is None: 

69 forces = atoms.get_forces(md=True) 

70 

71 p = atoms.get_momenta() 

72 p += 0.5 * self.dt * forces 

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

74 r = atoms.get_positions() 

75 

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

77 # RATTLE algorithm: 

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

79 if atoms.constraints: 

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

81 

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

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

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

85 # migrate along with the atoms. 

86 atoms.set_momenta(p, apply_constraint=False) 

87 

88 forces = atoms.get_forces(md=True) 

89 

90 # Second part of RATTLE will be done here: 

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

92 return forces