1"""Molecular Dynamics.""" 

2 

3import warnings 

4import numpy as np 

5 

6from ase.optimize.optimize import Dynamics 

7from ase.md.logger import MDLogger 

8from ase.io.trajectory import Trajectory 

9from ase import units 

10 

11 

12def process_temperature(temperature, temperature_K, orig_unit): 

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

14 

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

16 have the temperature specified in either Kelvin or Electron 

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

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

19 this. Using the original method then will issue a 

20 FutureWarning. 

21 

22 Four parameters: 

23 

24 temperature: None or float 

25 The original temperature specification in whatever unit was 

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

27 the historical unit was eV. 

28 

29 temperature_K: None or float 

30 Temperature in Kelvin. 

31 

32 orig_unit: str 

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

34 

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

36 None, otherwise an error is issued. 

37 

38 Return value: Temperature in Kelvin. 

39 """ 

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

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

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

43 if temperature is not None: 

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

45 if orig_unit == 'K': 

46 return temperature 

47 elif orig_unit == 'eV': 

48 warnings.warn(FutureWarning(w)) 

49 return temperature / units.kB 

50 else: 

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

52 

53 assert temperature_K is not None 

54 return temperature_K 

55 

56 

57class MolecularDynamics(Dynamics): 

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

59 

60 def __init__(self, atoms, timestep, trajectory, logfile=None, 

61 loginterval=1, append_trajectory=False): 

62 """Molecular Dynamics object. 

63 

64 Parameters: 

65 

66 atoms: Atoms object 

67 The Atoms object to operate on. 

68 

69 timestep: float 

70 The time step in ASE time units. 

71 

72 trajectory: Trajectory object or str 

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

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

75 trajectory. 

76 

77 logfile: file object or str (optional) 

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

79 Use '-' for stdout. 

80 

81 loginterval: int (optional) 

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

83 Default: 1 

84 

85 append_trajectory: boolean (optional) 

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

87 overwriten each time the dynamics is restarted from scratch. 

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

89 file instead. 

90 """ 

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

92 self.dt = timestep 

93 

94 Dynamics.__init__(self, atoms, logfile=None, trajectory=None) 

95 

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

97 self.max_steps = None 

98 

99 if 0 in self.masses: 

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

101 'likely lead to errors if the massless atoms ' 

102 'are unconstrained.') 

103 

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

105 

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

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

108 

109 # Trajectory is attached here instead of in Dynamics.__init__ 

110 # to respect the loginterval argument. 

111 if trajectory is not None: 

112 if isinstance(trajectory, str): 

113 mode = "a" if append_trajectory else "w" 

114 trajectory = self.closelater( 

115 Trajectory(trajectory, mode=mode, atoms=atoms) 

116 ) 

117 self.attach(trajectory, interval=loginterval) 

118 

119 if logfile: 

120 logger = self.closelater( 

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

122 self.attach(logger, loginterval) 

123 

124 def todict(self): 

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

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

127 'timestep': self.dt} 

128 

129 def irun(self, steps=50): 

130 """ Call Dynamics.irun and adjust max_steps """ 

131 self.max_steps = steps + self.nsteps 

132 return Dynamics.irun(self) 

133 

134 def run(self, steps=50): 

135 """ Call Dynamics.run and adjust max_steps """ 

136 self.max_steps = steps + self.nsteps 

137 return Dynamics.run(self) 

138 

139 def get_time(self): 

140 return self.nsteps * self.dt 

141 

142 def converged(self): 

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

144 return self.nsteps >= self.max_steps 

145 

146 def _get_com_velocity(self, velocity): 

147 """Return the center of mass velocity. 

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

149 """ 

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

151 

152 # Make the process_temperature function available to subclasses 

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

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

155 # as a function. 

156 _process_temperature = staticmethod(process_temperature)