Coverage for /builds/debichem-team/python-ase/ase/md/md.py: 80.43%

1"""Molecular Dynamics.""" 

2import warnings 

3from typing import IO, Optional, Union 

4 

5import numpy as np 

6 

7from ase import Atoms, units 

8from ase.md.logger import MDLogger 

9from ase.optimize.optimize import Dynamics 

10 

11 

12def process_temperature( 

13 temperature: Optional[float], 

14 temperature_K: Optional[float], 

15 orig_unit: str, 

16) -> float: 

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

18 

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

20 have the temperature specified in either Kelvin or Electron 

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

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

23 this. Using the original method then will issue a 

24 FutureWarning. 

25 

26 Four parameters: 

27 

28 temperature: None or float 

29 The original temperature specification in whatever unit was 

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

31 the historical unit was eV. 

32 

33 temperature_K: None or float 

34 Temperature in Kelvin. 

35 

36 orig_unit: str 

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

38 

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

40 None, otherwise an error is issued. 

41 

42 Return value: Temperature in Kelvin. 

43 """ 

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

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

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

47 if temperature is not None: 

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

49 if orig_unit == 'K': 

50 return temperature 

51 elif orig_unit == 'eV': 

52 warnings.warn(FutureWarning(w)) 

53 return temperature / units.kB 

54 else: 

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

56 

57 assert temperature_K is not None 

58 return temperature_K 

59 

60 

61class MolecularDynamics(Dynamics): 

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

63 

64 def __init__( 

65 self, 

66 atoms: Atoms, 

67 timestep: float, 

68 trajectory: Optional[str] = None, 

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

70 loginterval: int = 1, 

71 **kwargs, 

72 ): 

73 """Molecular Dynamics object. 

74 

75 Parameters 

76 ---------- 

77 atoms : Atoms object 

78 The Atoms object to operate on. 

79 

80 timestep : float 

81 The time step in ASE time units. 

82 

83 trajectory : Trajectory object or str 

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

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

86 trajectory. 

87 

88 logfile : file object or str (optional) 

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

90 Use '-' for stdout. 

91 

92 loginterval : int, default: 1 

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

94 

95 kwargs : dict, optional 

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

97 """ 

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

99 self.dt = timestep 

100 

101 super().__init__( 

102 atoms, 

103 logfile=None, 

104 trajectory=trajectory, 

105 loginterval=loginterval, 

106 **kwargs, 

107 ) 

108 

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

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

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

112 # work with dynamics. 

113 # 

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

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

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

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

118 self.atoms = atoms 

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

120 

121 if 0 in self.masses: 

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

123 'likely lead to errors if the massless atoms ' 

124 'are unconstrained.') 

125 

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

127 

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

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

130 

131 if logfile: 

132 logger = self.closelater( 

133 MDLogger(dyn=self, atoms=atoms, logfile=logfile)) 

134 self.attach(logger, loginterval) 

135 

136 def todict(self): 

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

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

139 'timestep': self.dt} 

140 

141 def irun(self, steps=50): 

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

143 

144 Parameters 

145 ---------- 

146 steps : int, default=DEFAULT_MAX_STEPS 

147 Number of molecular dynamics steps to be run. 

148 

149 Yields 

150 ------ 

151 converged : bool 

152 True if the maximum number of steps are reached. 

153 """ 

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

155 

156 def run(self, steps=50): 

157 """Run molecular dynamics algorithm. 

158 

159 Parameters 

160 ---------- 

161 steps : int, default=DEFAULT_MAX_STEPS 

162 Number of molecular dynamics steps to be run. 

163 

164 Returns 

165 ------- 

166 converged : bool 

167 True if the maximum number of steps are reached. 

168 """ 

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

170 

171 def get_time(self): 

172 return self.nsteps * self.dt 

173 

174 def converged(self): 

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

176 return self.nsteps >= self.max_steps 

177 

178 def _get_com_velocity(self, velocity): 

179 """Return the center of mass velocity. 

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

181 """ 

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

183 

184 # Make the process_temperature function available to subclasses 

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

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

187 # as a function. 

188 _process_temperature = staticmethod(process_temperature)