Coverage for /builds/kinetik161/ase/ase/md/nptberendsen.py: 63.86%
83 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1"""Berendsen NPT dynamics class."""
2import warnings
3from typing import IO, Optional, Union
5import numpy as np
7from ase import Atoms, units
8from ase.md.nvtberendsen import NVTBerendsen
11class NPTBerendsen(NVTBerendsen):
12 def __init__(
13 self,
14 atoms: Atoms,
15 timestep: float,
16 temperature: Optional[float] = None,
17 *,
18 temperature_K: Optional[float] = None,
19 pressure: Optional[float] = None,
20 pressure_au: Optional[float] = None,
21 taut: float = 0.5e3 * units.fs,
22 taup: float = 1e3 * units.fs,
23 compressibility: Optional[float] = None,
24 compressibility_au: Optional[float] = None,
25 fixcm: bool = True,
26 trajectory: Optional[str] = None,
27 logfile: Optional[Union[IO, str]] = None,
28 loginterval: int = 1,
29 append_trajectory: bool = False,
30 ):
31 """Berendsen (constant N, P, T) molecular dynamics.
33 This dynamics scale the velocities and volumes to maintain a constant
34 pressure and temperature. The shape of the simulation cell is not
35 altered, if that is desired use Inhomogenous_NPTBerendsen.
37 Parameters:
39 atoms: Atoms object
40 The list of atoms.
42 timestep: float
43 The time step in ASE time units.
45 temperature: float
46 The desired temperature, in Kelvin.
48 temperature_K: float
49 Alias for ``temperature``.
51 pressure: float (deprecated)
52 The desired pressure, in bar (1 bar = 1e5 Pa). Deprecated,
53 use ``pressure_au`` instead.
55 pressure: float
56 The desired pressure, in atomic units (eV/Å^3).
58 taut: float
59 Time constant for Berendsen temperature coupling in ASE
60 time units. Default: 0.5 ps.
62 taup: float
63 Time constant for Berendsen pressure coupling. Default: 1 ps.
65 compressibility: float (deprecated)
66 The compressibility of the material, in bar-1. Deprecated,
67 use ``compressibility_au`` instead.
69 compressibility_au: float
70 The compressibility of the material, in atomic units (Å^3/eV).
72 fixcm: bool (optional)
73 If True, the position and momentum of the center of mass is
74 kept unperturbed. Default: True.
76 trajectory: Trajectory object or str (optional)
77 Attach trajectory object. If *trajectory* is a string a
78 Trajectory will be constructed. Use *None* for no
79 trajectory.
81 logfile: file object or str (optional)
82 If *logfile* is a string, a file with that name will be opened.
83 Use '-' for stdout.
85 loginterval: int (optional)
86 Only write a log line for every *loginterval* time steps.
87 Default: 1
89 append_trajectory: boolean (optional)
90 Defaults to False, which causes the trajectory file to be
91 overwriten each time the dynamics is restarted from scratch.
92 If True, the new structures are appended to the trajectory
93 file instead.
96 """
98 NVTBerendsen.__init__(self, atoms, timestep, temperature=temperature,
99 temperature_K=temperature_K,
100 taut=taut, fixcm=fixcm, trajectory=trajectory,
101 logfile=logfile, loginterval=loginterval,
102 append_trajectory=append_trajectory)
103 self.taup = taup
104 self.pressure = self._process_pressure(pressure, pressure_au)
105 if compressibility is not None and compressibility_au is not None:
106 raise TypeError(
107 "Do not give both 'compressibility' and 'compressibility_au'")
108 if compressibility is not None:
109 # Specified in bar, convert to atomic units
110 warnings.warn(FutureWarning(
111 "Specify the compressibility in atomic units."))
112 self.set_compressibility(
113 compressibility_au=compressibility / (1e5 * units.Pascal))
114 else:
115 self.set_compressibility(compressibility_au=compressibility_au)
117 def set_taup(self, taup):
118 self.taup = taup
120 def get_taup(self):
121 return self.taup
123 def set_pressure(self, pressure=None, *, pressure_au=None,
124 pressure_bar=None):
125 self.pressure = self._process_pressure(pressure, pressure_bar,
126 pressure_au)
128 def get_pressure(self):
129 return self.pressure
131 def set_compressibility(self, *, compressibility_au):
132 self.compressibility = compressibility_au
134 def get_compressibility(self):
135 return self.compressibility
137 def set_timestep(self, timestep):
138 self.dt = timestep
140 def get_timestep(self):
141 return self.dt
143 def scale_positions_and_cell(self):
144 """ Do the Berendsen pressure coupling,
145 scale the atom position and the simulation cell."""
147 taupscl = self.dt / self.taup
148 stress = self.atoms.get_stress(voigt=False, include_ideal_gas=True)
149 old_pressure = -stress.trace() / 3
150 scl_pressure = (1.0 - taupscl * self.compressibility / 3.0 *
151 (self.pressure - old_pressure))
153 cell = self.atoms.get_cell()
154 cell = scl_pressure * cell
155 self.atoms.set_cell(cell, scale_atoms=True)
157 def step(self, forces=None):
158 """ move one timestep forward using Berenden NPT molecular dynamics."""
160 NVTBerendsen.scale_velocities(self)
161 self.scale_positions_and_cell()
163 # one step velocity verlet
164 atoms = self.atoms
166 if forces is None:
167 forces = atoms.get_forces(md=True)
169 p = self.atoms.get_momenta()
170 p += 0.5 * self.dt * forces
172 if self.fix_com:
173 # calculate the center of mass
174 # momentum and subtract it
175 psum = p.sum(axis=0) / float(len(p))
176 p = p - psum
178 self.atoms.set_positions(
179 self.atoms.get_positions() +
180 self.dt * p / self.atoms.get_masses()[:, np.newaxis])
182 # We need to store the momenta on the atoms before calculating
183 # the forces, as in a parallel Asap calculation atoms may
184 # migrate during force calculations, and the momenta need to
185 # migrate along with the atoms. For the same reason, we
186 # cannot use self.masses in the line above.
188 self.atoms.set_momenta(p)
189 forces = self.atoms.get_forces(md=True)
190 atoms.set_momenta(self.atoms.get_momenta() + 0.5 * self.dt * forces)
192 return forces
194 def _process_pressure(self, pressure, pressure_au):
195 """Handle that pressure can be specified in multiple units.
197 For at least a transition period, Berendsen NPT dynamics in ASE can
198 have the pressure specified in either bar or atomic units (eV/Å^3).
200 Two parameters:
202 pressure: None or float
203 The original pressure specification in bar.
204 A warning is issued if this is not None.
206 pressure_au: None or float
207 Pressure in ev/Å^3.
209 Exactly one of the two pressure parameters must be different from
210 None, otherwise an error is issued.
212 Return value: Pressure in eV/Å^3.
213 """
214 if (pressure is not None) + (pressure_au is not None) != 1:
215 raise TypeError("Exactly one of the parameters 'pressure',"
216 + " and 'pressure_au' must"
217 + " be given")
219 if pressure is not None:
220 w = ("The 'pressure' parameter is deprecated, please"
221 + " specify the pressure in atomic units (eV/Å^3)"
222 + " using the 'pressure_au' parameter.")
223 warnings.warn(FutureWarning(w))
224 return pressure * (1e5 * units.Pascal)
225 else:
226 return pressure_au
229class Inhomogeneous_NPTBerendsen(NPTBerendsen):
230 """Berendsen (constant N, P, T) molecular dynamics.
232 This dynamics scale the velocities and volumes to maintain a constant
233 pressure and temperature. The size of the unit cell is allowed to change
234 independently in the three directions, but the angles remain constant.
236 Usage: NPTBerendsen(atoms, timestep, temperature, taut, pressure, taup)
238 Parameters
239 ----------
240 mask : tuple[int]
241 Specifies which axes participate in the barostat. Default (1, 1, 1)
242 means that all axes participate, set any of them to zero to disable
243 the barostat in that direction.
244 """
246 def __init__(self, *args, mask=(1, 1, 1), **kwargs):
247 NPTBerendsen.__init__(self, *args, **kwargs)
248 self.mask = mask
250 def scale_positions_and_cell(self):
251 """ Do the Berendsen pressure coupling,
252 scale the atom position and the simulation cell."""
254 taupscl = self.dt * self.compressibility / self.taup / 3.0
255 stress = - self.atoms.get_stress(include_ideal_gas=True)
256 if stress.shape == (6,):
257 stress = stress[:3]
258 elif stress.shape == (3, 3):
259 stress = [stress[i][i] for i in range(3)]
260 else:
261 raise ValueError('Cannot use a stress tensor of shape ' +
262 str(stress.shape))
263 pbc = self.atoms.get_pbc()
264 scl_pressurex = 1.0 - taupscl * (self.pressure - stress[0]) \
265 * pbc[0] * self.mask[0]
266 scl_pressurey = 1.0 - taupscl * (self.pressure - stress[1]) \
267 * pbc[1] * self.mask[1]
268 scl_pressurez = 1.0 - taupscl * (self.pressure - stress[2]) \
269 * pbc[2] * self.mask[2]
270 cell = self.atoms.get_cell()
271 cell = np.array([scl_pressurex * cell[0],
272 scl_pressurey * cell[1],
273 scl_pressurez * cell[2]])
274 self.atoms.set_cell(cell, scale_atoms=True)