Coverage for ase / md / nvtberendsen.py: 80.00%

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

2 

3import warnings 

4 

5import numpy as np 

6 

7from ase import Atoms 

8from ase.md.md import MolecularDynamics 

9 

10 

11class NVTBerendsen(MolecularDynamics): 

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

13 

14 Berendsen *et al.*, J. Chem. Phys. 81, 3684–3690 (1984). 

15 https://doi.org/10.1063/1.448118 

16 """ 

17 

18 def __init__( 

19 self, 

20 atoms: Atoms, 

21 timestep: float, 

22 temperature: float | None = None, 

23 taut: float | None = None, 

24 fixcm: bool = True, 

25 *, 

26 temperature_K: float | None = None, 

27 **kwargs, 

28 ): 

29 """ 

30 Parameters 

31 ---------- 

32 atoms: Atoms object 

33 The list of atoms. 

34 

35 timestep: float 

36 The time step in ASE time units. 

37 

38 temperature: float 

39 The desired temperature, in Kelvin. 

40 

41 temperature_K: float 

42 Alias for *temperature* 

43 

44 taut: float 

45 Time constant for Berendsen temperature coupling in ASE 

46 time units. 

47 

48 fixcm: bool (optional) 

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

50 kept unperturbed. Default: True. 

51 

52 **kwargs : dict, optional 

53 Additional arguments passed to :class:~ase.md.md.MolecularDynamics 

54 base class. 

55 

56 """ 

57 if 'communicator' in kwargs: 

58 msg = ( 

59 '`communicator` has been deprecated since ASE 3.25.0 ' 

60 'and will be removed in ASE 3.26.0. Use `comm` instead.' 

61 ) 

62 warnings.warn(msg, FutureWarning) 

63 kwargs['comm'] = kwargs.pop('communicator') 

64 

65 MolecularDynamics.__init__(self, atoms, timestep, **kwargs) 

66 

67 if taut is None: 

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

69 self.taut = taut 

70 self.temperature = self._process_temperature( 

71 temperature, temperature_K, 'K' 

72 ) 

73 

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

75 

76 def set_taut(self, taut): 

77 self.taut = taut 

78 

79 def get_taut(self): 

80 return self.taut 

81 

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

83 self.temperature = self._process_temperature( 

84 temperature, temperature_K, 'K' 

85 ) 

86 

87 def get_temperature(self): 

88 return self.temperature 

89 

90 def set_timestep(self, timestep): 

91 self.dt = timestep 

92 

93 def get_timestep(self): 

94 return self.dt 

95 

96 def scale_velocities(self): 

97 """Do the NVT Berendsen velocity scaling""" 

98 tautscl = self.dt / self.taut 

99 old_temperature = self.atoms.get_temperature() 

100 

101 scl_temperature = np.sqrt( 

102 1.0 + (self.temperature / old_temperature - 1.0) * tautscl 

103 ) 

104 # Limit the velocity scaling to reasonable values 

105 if scl_temperature > 1.1: 

106 scl_temperature = 1.1 

107 if scl_temperature < 0.9: 

108 scl_temperature = 0.9 

109 

110 p = self.atoms.get_momenta() 

111 p = scl_temperature * p 

112 self.atoms.set_momenta(p) 

113 return 

114 

115 def step(self, forces=None): 

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

117 self.scale_velocities() 

118 

119 # one step velocity verlet 

120 atoms = self.atoms 

121 

122 if forces is None: 

123 forces = atoms.get_forces(md=True) 

124 

125 p = self.atoms.get_momenta() 

126 p += 0.5 * self.dt * forces 

127 

128 if self.fix_com: 

129 # calculate the center of mass 

130 # momentum and subtract it 

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

132 p = p - psum 

133 

134 self.atoms.set_positions( 

135 self.atoms.get_positions() 

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

137 ) 

138 

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

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

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

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

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

144 

145 self.atoms.set_momenta(p) 

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

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

148 

149 return forces