Coverage for /builds/ase/ase/ase/md/md.py: 80.43%

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 Dynamics 

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(Dynamics): 

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 self.attach(logger, loginterval) 

137 

138 def todict(self): 

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

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

141 'timestep': self.dt} 

142 

143 def irun(self, steps=50): 

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

145 

146 Parameters 

147 ---------- 

148 steps : int, default=DEFAULT_MAX_STEPS 

149 Number of molecular dynamics steps to be run. 

150 

151 Yields 

152 ------ 

153 converged : bool 

154 True if the maximum number of steps are reached. 

155 """ 

156 return Dynamics.irun(self, steps=steps) 

157 

158 def run(self, steps=50): 

159 """Run molecular dynamics algorithm. 

160 

161 Parameters 

162 ---------- 

163 steps : int, default=DEFAULT_MAX_STEPS 

164 Number of molecular dynamics steps to be run. 

165 

166 Returns 

167 ------- 

168 converged : bool 

169 True if the maximum number of steps are reached. 

170 """ 

171 return Dynamics.run(self, steps=steps) 

172 

173 def get_time(self): 

174 return self.nsteps * self.dt 

175 

176 def converged(self, gradient=None): 

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

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

179 return self.nsteps >= self.max_steps 

180 

181 def _get_com_velocity(self, velocity): 

182 """Return the center of mass velocity. 

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

184 """ 

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

186 

187 # Make the process_temperature function available to subclasses 

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

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

190 # as a function. 

191 _process_temperature = staticmethod(process_temperature)