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