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 #print("old_pressure", old_pressure, self.pressure) 140 #print("volume scaling by:", scl_pressure) 141 142 cell = self.atoms.get_cell() 143 cell = scl_pressure * cell 144 self.atoms.set_cell(cell, scale_atoms=True) 145 146 def step(self, forces=None): 147 """ move one timestep forward using Berenden NPT molecular dynamics.""" 148 149 NVTBerendsen.scale_velocities(self) 150 self.scale_positions_and_cell() 151 152 # one step velocity verlet 153 atoms = self.atoms 154 155 if forces is None: 156 forces = atoms.get_forces(md=True) 157 158 p = self.atoms.get_momenta() 159 p += 0.5 * self.dt * forces 160 161 if self.fix_com: 162 # calculate the center of mass 163 # momentum and subtract it 164 psum = p.sum(axis=0) / float(len(p)) 165 p = p - psum 166 167 self.atoms.set_positions( 168 self.atoms.get_positions() + 169 self.dt * p / self.atoms.get_masses()[:, np.newaxis]) 170 171 # We need to store the momenta on the atoms before calculating 172 # the forces, as in a parallel Asap calculation atoms may 173 # migrate during force calculations, and the momenta need to 174 # migrate along with the atoms. For the same reason, we 175 # cannot use self.masses in the line above. 176 177 self.atoms.set_momenta(p) 178 forces = self.atoms.get_forces(md=True) 179 atoms.set_momenta(self.atoms.get_momenta() + 0.5 * self.dt * forces) 180 181 return forces 182 183 def _process_pressure(self, pressure, pressure_au): 184 """Handle that pressure can be specified in multiple units. 185 186 For at least a transition period, Berendsen NPT dynamics in ASE can 187 have the pressure specified in either bar or atomic units (eV/Å^3). 188 189 Two parameters: 190 191 pressure: None or float 192 The original pressure specification in bar. 193 A warning is issued if this is not None. 194 195 pressure_au: None or float 196 Pressure in ev/Å^3. 197 198 Exactly one of the two pressure parameters must be different from 199 None, otherwise an error is issued. 200 201 Return value: Pressure in eV/Å^3. 202 """ 203 if (pressure is not None) + (pressure_au is not None) != 1: 204 raise TypeError("Exactly one of the parameters 'pressure'," 205 + " and 'pressure_au' must" 206 + " be given") 207 208 if pressure is not None: 209 w = ("The 'pressure' parameter is deprecated, please" 210 + " specify the pressure in atomic units (eV/Å^3)" 211 + " using the 'pressure_au' parameter.") 212 warnings.warn(FutureWarning(w)) 213 return pressure * (1e5 * units.Pascal) 214 else: 215 return pressure_au 216 217 218class Inhomogeneous_NPTBerendsen(NPTBerendsen): 219 """Berendsen (constant N, P, T) molecular dynamics. 220 221 This dynamics scale the velocities and volumes to maintain a constant 222 pressure and temperature. The size of the unit cell is allowed to change 223 independently in the three directions, but the angles remain constant. 224 225 Usage: NPTBerendsen(atoms, timestep, temperature, taut, pressure, taup) 226 227 atoms 228 The list of atoms. 229 230 timestep 231 The time step. 232 233 temperature 234 The desired temperature, in Kelvin. 235 236 taut 237 Time constant for Berendsen temperature coupling. 238 239 fixcm 240 If True, the position and momentum of the center of mass is 241 kept unperturbed. Default: True. 242 243 pressure 244 The desired pressure, in bar (1 bar = 1e5 Pa). 245 246 taup 247 Time constant for Berendsen pressure coupling. 248 249 compressibility 250 The compressibility of the material, water 4.57E-5 bar-1, in bar-1 251 252 mask 253 Specifies which axes participate in the barostat. Default (1, 1, 1) 254 means that all axes participate, set any of them to zero to disable 255 the barostat in that direction. 256 """ 257 258 def __init__(self, atoms, timestep, temperature=None, 259 *, temperature_K=None, 260 taut=0.5e3 * units.fs, pressure=None, 261 pressure_au=None, taup=1e3 * units.fs, 262 compressibility=None, compressibility_au=None, 263 mask=(1, 1, 1), fixcm=True, trajectory=None, 264 logfile=None, loginterval=1): 265 266 NPTBerendsen.__init__(self, atoms, timestep, temperature=temperature, 267 temperature_K=temperature_K, 268 taut=taut, taup=taup, pressure=pressure, 269 pressure_au=pressure_au, 270 compressibility=compressibility, 271 compressibility_au=compressibility_au, 272 fixcm=fixcm, trajectory=trajectory, 273 logfile=logfile, loginterval=loginterval) 274 self.mask = mask 275 276 def scale_positions_and_cell(self): 277 """ Do the Berendsen pressure coupling, 278 scale the atom position and the simulation cell.""" 279 280 taupscl = self.dt * self.compressibility / self.taup / 3.0 281 stress = - self.atoms.get_stress(include_ideal_gas=True) 282 if stress.shape == (6,): 283 stress = stress[:3] 284 elif stress.shape == (3, 3): 285 stress = [stress[i][i] for i in range(3)] 286 else: 287 raise ValueError('Cannot use a stress tensor of shape ' + 288 str(stress.shape)) 289 pbc = self.atoms.get_pbc() 290 scl_pressurex = 1.0 - taupscl * (self.pressure - stress[0]) \ 291 * pbc[0] * self.mask[0] 292 scl_pressurey = 1.0 - taupscl * (self.pressure - stress[1]) \ 293 * pbc[1] * self.mask[1] 294 scl_pressurez = 1.0 - taupscl * (self.pressure - stress[2]) \ 295 * pbc[2] * self.mask[2] 296 cell = self.atoms.get_cell() 297 cell = np.array([scl_pressurex * cell[0], 298 scl_pressurey * cell[1], 299 scl_pressurez * cell[2]]) 300 self.atoms.set_cell(cell, scale_atoms=True) 301