Coverage for ase / md / md.py: 84.13%

1# fmt: off 

2 

3"""Molecular Dynamics.""" 

4import warnings 

5from typing import IO, Optional, Union 

6 

7import numpy as np 

8 

9from ase import Atoms, units 

10from ase.md.logger import MDLogger 

11from ase.optimize.optimize import BaseDynamics 

12 

13 

14def process_temperature( 

15 temperature: Optional[float], 

16 temperature_K: Optional[float], 

17 orig_unit: str, 

18) -> float: 

19 """Handle that temperature can be specified in multiple units. 

20 

21 For at least a transition period, molecular dynamics in ASE can 

22 have the temperature specified in either Kelvin or Electron 

23 Volt. The different MD algorithms had different defaults, by 

24 forcing the user to explicitly choose a unit we can resolve 

25 this. Using the original method then will issue a 

26 FutureWarning. 

27 

28 Four parameters: 

29 

30 temperature: None or float 

31 The original temperature specification in whatever unit was 

32 historically used. A warning is issued if this is not None and 

33 the historical unit was eV. 

34 

35 temperature_K: None or float 

36 Temperature in Kelvin. 

37 

38 orig_unit: str 

39 Unit used for the `temperature`` parameter. Must be 'K' or 'eV'. 

40 

41 Exactly one of the two temperature parameters must be different from 

42 None, otherwise an error is issued. 

43 

44 Return value: Temperature in Kelvin. 

45 """ 

46 if (temperature is not None) + (temperature_K is not None) != 1: 

47 raise TypeError("Exactly one of the parameters 'temperature'," 

48 + " and 'temperature_K', must be given") 

49 if temperature is not None: 

50 w = "Specify the temperature in K using the 'temperature_K' argument" 

51 if orig_unit == 'K': 

52 return temperature 

53 elif orig_unit == 'eV': 

54 warnings.warn(FutureWarning(w)) 

55 return temperature / units.kB 

56 else: 

57 raise ValueError("Unknown temperature unit " + orig_unit) 

58 

59 assert temperature_K is not None 

60 return temperature_K 

61 

62 

63class MolecularDynamics(BaseDynamics): 

64 """Base-class for all MD classes.""" 

65 

66 def __init__( 

67 self, 

68 atoms: Atoms, 

69 timestep: float, 

70 trajectory: Optional[str] = None, 

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

72 loginterval: int = 1, 

73 **kwargs, 

74 ): 

75 """Molecular Dynamics object. 

76 

77 Parameters 

78 ---------- 

79 atoms : Atoms object 

80 The Atoms object to operate on. 

81 

82 timestep : float 

83 The time step in ASE time units. 

84 

85 trajectory : Trajectory object or str 

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

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

88 trajectory. 

89 

90 logfile : file object or str (optional) 

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

92 Use '-' for stdout. 

93 

94 loginterval : int, default: 1 

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

96 

97 kwargs : dict, optional 

98 Extra arguments passed to :class:`~ase.optimize.optimize.Dynamics`. 

99 """ 

100 # dt as to be attached _before_ parent class is initialized 

101 self.dt = timestep 

102 

103 super().__init__( 

104 atoms, 

105 logfile=None, 

106 trajectory=trajectory, 

107 loginterval=loginterval, 

108 **kwargs, 

109 ) 

110 

111 # Some codes (e.g. Asap) may be using filters to 

112 # constrain atoms or do other things. Current state of the art 

113 # is that "atoms" must be either Atoms or Filter in order to 

114 # work with dynamics. 

115 # 

116 # In the future, we should either use a special role interface 

117 # for MD, or we should ensure that the input is *always* a Filter. 

118 # That way we won't need to test multiple cases. Currently, 

119 # we do not test /any/ kind of MD with any kind of Filter in ASE. 

120 self.atoms = atoms 

121 self.masses = self.atoms.get_masses() 

122 

123 if 0 in self.masses: 

124 warnings.warn('Zero mass encountered in atoms; this will ' 

125 'likely lead to errors if the massless atoms ' 

126 'are unconstrained.') 

127 

128 self.masses.shape = (-1, 1) 

129 

130 if not self.atoms.has('momenta'): 

131 self.atoms.set_momenta(np.zeros([len(self.atoms), 3])) 

132 

133 if logfile: 

134 logger = self.closelater( 

135 MDLogger(dyn=self, atoms=atoms, logfile=logfile, 

136 comm=self.comm)) 

137 self.attach(logger, loginterval) 

138 

139 def todict(self): 

140 return {'type': 'molecular-dynamics', 

141 'md-type': self.__class__.__name__, 

142 'timestep': self.dt} 

143 

144 def irun(self, steps=50): 

145 """Run molecular dynamics algorithm as a generator. 

146 

147 Parameters 

148 ---------- 

149 steps : int, default=DEFAULT_MAX_STEPS 

150 Number of molecular dynamics steps to be run. 

151 

152 Yields 

153 ------ 

154 converged : bool 

155 True if the maximum number of steps are reached. 

156 """ 

157 # update the maximum number of steps 

158 self.max_steps = self.nsteps + steps 

159 

160 if self.nsteps == 0: 

161 # For historical reasons we do a magical incantation 

162 # here with forces, log, observers. 

163 self.atoms.get_forces() 

164 self.log() 

165 self.call_observers() 

166 

167 yield self.nsteps == self.max_steps 

168 

169 # run the algorithm until converged or max_steps reached 

170 while self.nsteps < self.max_steps: 

171 self.step() 

172 self.nsteps += 1 

173 self.atoms.get_forces() 

174 self.log() 

175 self.call_observers() 

176 yield self.nsteps == self.max_steps 

177 

178 def log(self): 

179 # Subclasses can override this. 

180 return None 

181 

182 def run(self, steps=50): 

183 """Run molecular dynamics algorithm. 

184 

185 Parameters 

186 ---------- 

187 steps : int, default=DEFAULT_MAX_STEPS 

188 Number of molecular dynamics steps to be run. 

189 

190 Returns 

191 ------- 

192 converged : bool 

193 True if the maximum number of steps are reached. 

194 """ 

195 converged = steps == 0 

196 for converged in self.irun(steps=steps): 

197 pass 

198 return converged 

199 

200 def get_time(self): 

201 return self.nsteps * self.dt 

202 

203 def converged(self, gradient=None): 

204 """ MD is 'converged' when number of maximum steps is reached. """ 

205 # We take gradient now (due to optimizers). Should refactor. 

206 return self.nsteps >= self.max_steps 

207 

208 def _get_com_velocity(self, velocity): 

209 """Return the center of mass velocity. 

210 Internal use only. This function can be reimplemented by Asap. 

211 """ 

212 return np.dot(self.masses.ravel(), velocity) / self.masses.sum() 

213 

214 # Make the process_temperature function available to subclasses 

215 # as a static method. This makes it easy for MD objects to use 

216 # it, while functions in md.velocitydistribution have access to it 

217 # as a function. 

218 _process_temperature = staticmethod(process_temperature)