1"""Berendsen NVT dynamics class.""" 

2 

3import numpy as np 

4from ase.md.md import MolecularDynamics 

5from ase.parallel import world 

6 

7 

8class NVTBerendsen(MolecularDynamics): 

9 def __init__(self, atoms, timestep, temperature=None, taut=None, 

10 fixcm=True, *, temperature_K=None, 

11 trajectory=None, logfile=None, loginterval=1, 

12 communicator=world, append_trajectory=False): 

13 """Berendsen (constant N, V, T) molecular dynamics. 

14 

15 Parameters: 

16 

17 atoms: Atoms object 

18 The list of atoms. 

19 

20 timestep: float 

21 The time step in ASE time units. 

22 

23 temperature: float 

24 The desired temperature, in Kelvin. 

25 

26 temperature_K: float 

27 Alias for *temperature* 

28 

29 taut: float 

30 Time constant for Berendsen temperature coupling in ASE 

31 time units. 

32 

33 fixcm: bool (optional) 

34 If True, the position and momentum of the center of mass is 

35 kept unperturbed. Default: True. 

36 

37 trajectory: Trajectory object or str (optional) 

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

39 Trajectory will be constructed. Use *None* for no 

40 trajectory. 

41 

42 logfile: file object or str (optional) 

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

44 Use '-' for stdout. 

45 

46 loginterval: int (optional) 

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

48 Default: 1 

49 

50 append_trajectory: boolean (optional) 

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

52 overwriten each time the dynamics is restarted from scratch. 

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

54 file instead. 

55 

56 """ 

57 

58 MolecularDynamics.__init__(self, atoms, timestep, trajectory, 

59 logfile, loginterval, 

60 append_trajectory=append_trajectory) 

61 if taut is None: 

62 raise TypeError("Missing 'taut' argument.") 

63 self.taut = taut 

64 self.temperature = self._process_temperature(temperature, 

65 temperature_K, 'K') 

66 

67 self.fix_com = fixcm # will the center of mass be held fixed? 

68 self.communicator = communicator 

69 

70 def set_taut(self, taut): 

71 self.taut = taut 

72 

73 def get_taut(self): 

74 return self.taut 

75 

76 def set_temperature(self, temperature=None, *, temperature_K=None): 

77 self.temperature = self._process_temperature(temperature, 

78 temperature_K, 'K') 

79 

80 def get_temperature(self): 

81 return self.temperature 

82 

83 def set_timestep(self, timestep): 

84 self.dt = timestep 

85 

86 def get_timestep(self): 

87 return self.dt 

88 

89 def scale_velocities(self): 

90 """ Do the NVT Berendsen velocity scaling """ 

91 tautscl = self.dt / self.taut 

92 old_temperature = self.atoms.get_temperature() 

93 

94 scl_temperature = np.sqrt(1.0 + 

95 (self.temperature / old_temperature - 1.0) * 

96 tautscl) 

97 # Limit the velocity scaling to reasonable values 

98 if scl_temperature > 1.1: 

99 scl_temperature = 1.1 

100 if scl_temperature < 0.9: 

101 scl_temperature = 0.9 

102 

103 p = self.atoms.get_momenta() 

104 p = scl_temperature * p 

105 self.atoms.set_momenta(p) 

106 return 

107 

108 def step(self, forces=None): 

109 """Move one timestep forward using Berenden NVT molecular dynamics.""" 

110 self.scale_velocities() 

111 

112 # one step velocity verlet 

113 atoms = self.atoms 

114 

115 if forces is None: 

116 forces = atoms.get_forces(md=True) 

117 

118 p = self.atoms.get_momenta() 

119 p += 0.5 * self.dt * forces 

120 

121 if self.fix_com: 

122 # calculate the center of mass 

123 # momentum and subtract it 

124 psum = p.sum(axis=0) / float(len(p)) 

125 p = p - psum 

126 

127 self.atoms.set_positions( 

128 self.atoms.get_positions() + 

129 self.dt * p / self.atoms.get_masses()[:, np.newaxis]) 

130 

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

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

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

134 # migrate along with the atoms. For the same reason, we 

135 # cannot use self.masses in the line above. 

136 

137 self.atoms.set_momenta(p) 

138 forces = self.atoms.get_forces(md=True) 

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

140 

141 return forces