1"""Berendsen NPT dynamics class.""" 

2 

3import numpy as np 

4import warnings 

5 

6from ase.md.nvtberendsen import NVTBerendsen 

7import ase.units as units 

8 

9 

10class NPTBerendsen(NVTBerendsen): 

11 def __init__(self, atoms, timestep, temperature=None, 

12 *, temperature_K=None, pressure=None, pressure_au=None, 

13 taut=0.5e3 * units.fs, taup=1e3 * units.fs, 

14 compressibility=None, compressibility_au=None, fixcm=True, 

15 trajectory=None, 

16 logfile=None, loginterval=1, append_trajectory=False): 

17 """Berendsen (constant N, P, T) molecular dynamics. 

18 

19 This dynamics scale the velocities and volumes to maintain a constant 

20 pressure and temperature. The shape of the simulation cell is not 

21 altered, if that is desired use Inhomogenous_NPTBerendsen. 

22 

23 Parameters: 

24 

25 atoms: Atoms object 

26 The list of atoms. 

27 

28 timestep: float 

29 The time step in ASE time units. 

30 

31 temperature: float 

32 The desired temperature, in Kelvin. 

33 

34 temperature_K: float 

35 Alias for ``temperature``. 

36 

37 pressure: float (deprecated) 

38 The desired pressure, in bar (1 bar = 1e5 Pa). Deprecated, 

39 use ``pressure_au`` instead. 

40 

41 pressure: float 

42 The desired pressure, in atomic units (eV/Å^3). 

43 

44 taut: float 

45 Time constant for Berendsen temperature coupling in ASE 

46 time units. Default: 0.5 ps. 

47 

48 taup: float 

49 Time constant for Berendsen pressure coupling. Default: 1 ps. 

50 

51 compressibility: float (deprecated) 

52 The compressibility of the material, in bar-1. Deprecated, 

53 use ``compressibility_au`` instead. 

54 

55 compressibility_au: float 

56 The compressibility of the material, in atomic units (Å^3/eV). 

57 

58 fixcm: bool (optional) 

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

60 kept unperturbed. Default: True. 

61 

62 trajectory: Trajectory object or str (optional) 

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

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

65 trajectory. 

66 

67 logfile: file object or str (optional) 

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

69 Use '-' for stdout. 

70 

71 loginterval: int (optional) 

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

73 Default: 1 

74 

75 append_trajectory: boolean (optional) 

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

77 overwriten each time the dynamics is restarted from scratch. 

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

79 file instead. 

80 

81 

82 """ 

83 

84 NVTBerendsen.__init__(self, atoms, timestep, temperature=temperature, 

85 temperature_K=temperature_K, 

86 taut=taut, fixcm=fixcm, trajectory=trajectory, 

87 logfile=logfile, loginterval=loginterval, 

88 append_trajectory=append_trajectory) 

89 self.taup = taup 

90 self.pressure = self._process_pressure(pressure, pressure_au) 

91 if compressibility is not None and compressibility_au is not None: 

92 raise TypeError( 

93 "Do not give both 'compressibility' and 'compressibility_au'") 

94 if compressibility is not None: 

95 # Specified in bar, convert to atomic units 

96 warnings.warn(FutureWarning( 

97 "Specify the compressibility in atomic units.")) 

98 self.set_compressibility( 

99 compressibility_au=compressibility / (1e5 * units.Pascal)) 

100 else: 

101 self.set_compressibility(compressibility_au=compressibility_au) 

102 

103 def set_taup(self, taup): 

104 self.taup = taup 

105 

106 def get_taup(self): 

107 return self.taup 

108 

109 def set_pressure(self, pressure=None, *, pressure_au=None, 

110 pressure_bar=None): 

111 self.pressure = self._process_pressure(pressure, pressure_bar, 

112 pressure_au) 

113 

114 def get_pressure(self): 

115 return self.pressure 

116 

117 def set_compressibility(self, *, compressibility_au): 

118 self.compressibility = compressibility_au 

119 

120 def get_compressibility(self): 

121 return self.compressibility 

122 

123 def set_timestep(self, timestep): 

124 self.dt = timestep 

125 

126 def get_timestep(self): 

127 return self.dt 

128 

129 def scale_positions_and_cell(self): 

130 """ Do the Berendsen pressure coupling, 

131 scale the atom position and the simulation cell.""" 

132 

133 taupscl = self.dt / self.taup 

134 stress = self.atoms.get_stress(voigt=False, include_ideal_gas=True) 

135 old_pressure = -stress.trace() / 3 

136 scl_pressure = (1.0 - taupscl * self.compressibility / 3.0 * 

137 (self.pressure - old_pressure)) 

138 

139 cell = self.atoms.get_cell() 

140 cell = scl_pressure * cell 

141 self.atoms.set_cell(cell, scale_atoms=True) 

142 

143 def step(self, forces=None): 

144 """ move one timestep forward using Berenden NPT molecular dynamics.""" 

145 

146 NVTBerendsen.scale_velocities(self) 

147 self.scale_positions_and_cell() 

148 

149 # one step velocity verlet 

150 atoms = self.atoms 

151 

152 if forces is None: 

153 forces = atoms.get_forces(md=True) 

154 

155 p = self.atoms.get_momenta() 

156 p += 0.5 * self.dt * forces 

157 

158 if self.fix_com: 

159 # calculate the center of mass 

160 # momentum and subtract it 

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

162 p = p - psum 

163 

164 self.atoms.set_positions( 

165 self.atoms.get_positions() + 

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

167 

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

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

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

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

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

173 

174 self.atoms.set_momenta(p) 

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

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

177 

178 return forces 

179 

180 def _process_pressure(self, pressure, pressure_au): 

181 """Handle that pressure can be specified in multiple units. 

182 

183 For at least a transition period, Berendsen NPT dynamics in ASE can 

184 have the pressure specified in either bar or atomic units (eV/Å^3). 

185 

186 Two parameters: 

187 

188 pressure: None or float 

189 The original pressure specification in bar. 

190 A warning is issued if this is not None. 

191 

192 pressure_au: None or float 

193 Pressure in ev/Å^3. 

194 

195 Exactly one of the two pressure parameters must be different from 

196 None, otherwise an error is issued. 

197 

198 Return value: Pressure in eV/Å^3. 

199 """ 

200 if (pressure is not None) + (pressure_au is not None) != 1: 

201 raise TypeError("Exactly one of the parameters 'pressure'," 

202 + " and 'pressure_au' must" 

203 + " be given") 

204 

205 if pressure is not None: 

206 w = ("The 'pressure' parameter is deprecated, please" 

207 + " specify the pressure in atomic units (eV/Å^3)" 

208 + " using the 'pressure_au' parameter.") 

209 warnings.warn(FutureWarning(w)) 

210 return pressure * (1e5 * units.Pascal) 

211 else: 

212 return pressure_au 

213 

214 

215class Inhomogeneous_NPTBerendsen(NPTBerendsen): 

216 """Berendsen (constant N, P, T) molecular dynamics. 

217 

218 This dynamics scale the velocities and volumes to maintain a constant 

219 pressure and temperature. The size of the unit cell is allowed to change 

220 independently in the three directions, but the angles remain constant. 

221 

222 Usage: NPTBerendsen(atoms, timestep, temperature, taut, pressure, taup) 

223 

224 atoms 

225 The list of atoms. 

226 

227 timestep 

228 The time step. 

229 

230 temperature 

231 The desired temperature, in Kelvin. 

232 

233 taut 

234 Time constant for Berendsen temperature coupling. 

235 

236 fixcm 

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

238 kept unperturbed. Default: True. 

239 

240 pressure 

241 The desired pressure, in bar (1 bar = 1e5 Pa). 

242 

243 taup 

244 Time constant for Berendsen pressure coupling. 

245 

246 compressibility 

247 The compressibility of the material, water 4.57E-5 bar-1, in bar-1 

248 

249 mask 

250 Specifies which axes participate in the barostat. Default (1, 1, 1) 

251 means that all axes participate, set any of them to zero to disable 

252 the barostat in that direction. 

253 """ 

254 

255 def __init__(self, atoms, timestep, temperature=None, 

256 *, temperature_K=None, 

257 taut=0.5e3 * units.fs, pressure=None, 

258 pressure_au=None, taup=1e3 * units.fs, 

259 compressibility=None, compressibility_au=None, 

260 mask=(1, 1, 1), fixcm=True, trajectory=None, 

261 logfile=None, loginterval=1): 

262 

263 NPTBerendsen.__init__(self, atoms, timestep, temperature=temperature, 

264 temperature_K=temperature_K, 

265 taut=taut, taup=taup, pressure=pressure, 

266 pressure_au=pressure_au, 

267 compressibility=compressibility, 

268 compressibility_au=compressibility_au, 

269 fixcm=fixcm, trajectory=trajectory, 

270 logfile=logfile, loginterval=loginterval) 

271 self.mask = mask 

272 

273 def scale_positions_and_cell(self): 

274 """ Do the Berendsen pressure coupling, 

275 scale the atom position and the simulation cell.""" 

276 

277 taupscl = self.dt * self.compressibility / self.taup / 3.0 

278 stress = - self.atoms.get_stress(include_ideal_gas=True) 

279 if stress.shape == (6,): 

280 stress = stress[:3] 

281 elif stress.shape == (3, 3): 

282 stress = [stress[i][i] for i in range(3)] 

283 else: 

284 raise ValueError('Cannot use a stress tensor of shape ' + 

285 str(stress.shape)) 

286 pbc = self.atoms.get_pbc() 

287 scl_pressurex = 1.0 - taupscl * (self.pressure - stress[0]) \ 

288 * pbc[0] * self.mask[0] 

289 scl_pressurey = 1.0 - taupscl * (self.pressure - stress[1]) \ 

290 * pbc[1] * self.mask[1] 

291 scl_pressurez = 1.0 - taupscl * (self.pressure - stress[2]) \ 

292 * pbc[2] * self.mask[2] 

293 cell = self.atoms.get_cell() 

294 cell = np.array([scl_pressurex * cell[0], 

295 scl_pressurey * cell[1], 

296 scl_pressurez * cell[2]]) 

297 self.atoms.set_cell(cell, scale_atoms=True)