Coverage for /builds/debichem-team/python-ase/ase/md/velocitydistribution.py : 82.30%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# VelocityDistributions.py -- set up a velocity distribution
3"""Module for setting up velocity distributions such as Maxwell–Boltzmann.
5Currently, only a few functions are defined, such as
6MaxwellBoltzmannDistribution, which sets the momenta of a list of
7atoms according to a Maxwell-Boltzmann distribution at a given
8temperature.
10"""
12import numpy as np
13from ase.parallel import world
14from ase import units
15from ase.md.md import process_temperature
17# define a ``zero'' temperature to avoid divisions by zero
18eps_temp = 1e-12
21class UnitError(Exception):
22 """Exception raised when wrong units are specified"""
25def force_temperature(atoms, temperature, unit="K"):
26 """ force (nucl.) temperature to have a precise value
28 Parameters:
29 atoms: ase.Atoms
30 the structure
31 temperature: float
32 nuclear temperature to set
33 unit: str
34 'K' or 'eV' as unit for the temperature
35 """
37 if unit == "K":
38 E_temp = temperature * units.kB
39 elif unit == "eV":
40 E_temp = temperature
41 else:
42 raise UnitError("'{}' is not supported, use 'K' or 'eV'.".format(unit))
44 if temperature > eps_temp:
45 E_kin0 = atoms.get_kinetic_energy() / len(atoms) / 1.5
46 gamma = E_temp / E_kin0
47 else:
48 gamma = 0.0
49 atoms.set_momenta(atoms.get_momenta() * np.sqrt(gamma))
52def _maxwellboltzmanndistribution(masses, temp, communicator=None, rng=None):
53 """Return a Maxwell-Boltzmann distribution with a given temperature.
55 Paremeters:
57 masses: float
58 The atomic masses.
60 temp: float
61 The temperature in electron volt.
63 communicator: MPI communicator (optional)
64 Communicator used to distribute an identical distribution to
65 all tasks. Set to 'serial' to disable communication (setting to None
66 gives the default). Default: ase.parallel.world
68 rng: numpy RNG (optional)
69 The random number generator. Default: np.random
71 Returns:
73 A numpy array with Maxwell-Boltzmann distributed momenta.
74 """
75 if rng is None:
76 rng = np.random
77 if communicator is None:
78 communicator = world
79 xi = rng.standard_normal((len(masses), 3))
80 if communicator != 'serial':
81 communicator.broadcast(xi, 0)
82 momenta = xi * np.sqrt(masses * temp)[:, np.newaxis]
83 return momenta
86def MaxwellBoltzmannDistribution(
87 atoms, temp=None, *, temperature_K=None,
88 communicator=None, force_temp=False,
89 rng=None
90):
91 """Sets the momenta to a Maxwell-Boltzmann distribution.
93 Parameters:
95 atoms: Atoms object
96 The atoms. Their momenta will be modified.
98 temp: float (deprecated)
99 The temperature in eV. Deprecated, used temperature_K instead.
101 temperature_K: float
102 The temperature in Kelvin.
104 communicator: MPI communicator (optional)
105 Communicator used to distribute an identical distribution to
106 all tasks. Set to 'serial' to disable communication. Leave as None to
107 get the default: ase.parallel.world
109 force_temp: bool (optinal, default: False)
110 If True, random the momenta are rescaled so the kinetic energy is
111 exactly 3/2 N k T. This is a slight deviation from the correct
112 Maxwell-Boltzmann distribution.
114 rng: Numpy RNG (optional)
115 Random number generator. Default: numpy.random
116 """
117 temp = units.kB * process_temperature(temp, temperature_K, 'eV')
118 masses = atoms.get_masses()
119 momenta = _maxwellboltzmanndistribution(masses, temp, communicator, rng)
120 atoms.set_momenta(momenta)
121 if force_temp:
122 force_temperature(atoms, temperature=temp, unit="eV")
125def Stationary(atoms, preserve_temperature=True):
126 "Sets the center-of-mass momentum to zero."
128 # Save initial temperature
129 temp0 = atoms.get_temperature()
131 p = atoms.get_momenta()
132 p0 = np.sum(p, 0)
133 # We should add a constant velocity, not momentum, to the atoms
134 m = atoms.get_masses()
135 mtot = np.sum(m)
136 v0 = p0 / mtot
137 p -= v0 * m[:, np.newaxis]
138 atoms.set_momenta(p)
140 if preserve_temperature:
141 force_temperature(atoms, temp0)
144def ZeroRotation(atoms, preserve_temperature=True):
145 "Sets the total angular momentum to zero by counteracting rigid rotations."
147 # Save initial temperature
148 temp0 = atoms.get_temperature()
150 # Find the principal moments of inertia and principal axes basis vectors
151 Ip, basis = atoms.get_moments_of_inertia(vectors=True)
152 # Calculate the total angular momentum and transform to principal basis
153 Lp = np.dot(basis, atoms.get_angular_momentum())
154 # Calculate the rotation velocity vector in the principal basis, avoiding
155 # zero division, and transform it back to the cartesian coordinate system
156 omega = np.dot(np.linalg.inv(basis), np.select([Ip > 0], [Lp / Ip]))
157 # We subtract a rigid rotation corresponding to this rotation vector
158 com = atoms.get_center_of_mass()
159 positions = atoms.get_positions()
160 positions -= com # translate center of mass to origin
161 velocities = atoms.get_velocities()
162 atoms.set_velocities(velocities - np.cross(omega, positions))
164 if preserve_temperature:
165 force_temperature(atoms, temp0)
168def n_BE(temp, omega):
169 """Bose-Einstein distribution function.
171 Args:
172 temp: temperature converted to eV (*units.kB)
173 omega: sequence of frequencies converted to eV
175 Returns:
176 Value of Bose-Einstein distribution function for each energy
178 """
180 omega = np.asarray(omega)
182 # 0K limit
183 if temp < eps_temp:
184 n = np.zeros_like(omega)
185 else:
186 n = 1 / (np.exp(omega / (temp)) - 1)
187 return n
190def phonon_harmonics(
191 force_constants,
192 masses,
193 temp=None,
194 *,
195 temperature_K=None,
196 rng=np.random.rand,
197 quantum=False,
198 plus_minus=False,
199 return_eigensolution=False,
200 failfast=True,
201):
202 r"""Return displacements and velocities that produce a given temperature.
204 Parameters:
206 force_constants: array of size 3N x 3N
207 force constants (Hessian) of the system in eV/Ų
209 masses: array of length N
210 masses of the structure in amu
212 temp: float (deprecated)
213 Temperature converted to eV (T * units.kB). Deprecated, use
214 ``temperature_K``.
216 temperature_K: float
217 Temperature in Kelvin.
219 rng: function
220 Random number generator function, e.g., np.random.rand
222 quantum: bool
223 True for Bose-Einstein distribution, False for Maxwell-Boltzmann
224 (classical limit)
226 plus_minus: bool
227 Displace atoms with +/- the amplitude accoding to PRB 94, 075125
229 return_eigensolution: bool
230 return eigenvalues and eigenvectors of the dynamical matrix
232 failfast: bool
233 True for sanity checking the phonon spectrum for negative
234 frequencies at Gamma
236 Returns:
238 Displacements, velocities generated from the eigenmodes,
239 (optional: eigenvalues, eigenvectors of dynamical matrix)
241 Purpose:
243 Excite phonon modes to specified temperature.
245 This excites all phonon modes randomly so that each contributes,
246 on average, equally to the given temperature. Both potential
247 energy and kinetic energy will be consistent with the phononic
248 vibrations characteristic of the specified temperature.
250 In other words the system will be equilibrated for an MD run at
251 that temperature.
253 force_constants should be the matrix as force constants, e.g.,
254 as computed by the ase.phonons module.
256 Let X_ai be the phonon modes indexed by atom and mode, w_i the
257 phonon frequencies, and let 0 < Q_i <= 1 and 0 <= R_i < 1 be
258 uniformly random numbers. Then
260 .. code-block:: none
263 1/2
264 _ / k T \ --- 1 _ 1/2
265 R += | --- | > --- X (-2 ln Q ) cos (2 pi R )
266 a \ m / --- w ai i i
267 a i i
270 1/2
271 _ / k T \ --- _ 1/2
272 v = | --- | > X (-2 ln Q ) sin (2 pi R )
273 a \ m / --- ai i i
274 a i
276 Reference: [West, Estreicher; PRL 96, 22 (2006)]
277 """
279 # Handle the temperature units
280 temp = units.kB * process_temperature(temp, temperature_K, 'eV')
282 # Build dynamical matrix
283 rminv = (masses ** -0.5).repeat(3)
284 dynamical_matrix = force_constants * rminv[:, None] * rminv[None, :]
286 # Solve eigenvalue problem to compute phonon spectrum and eigenvectors
287 w2_s, X_is = np.linalg.eigh(dynamical_matrix)
289 # Check for soft modes
290 if failfast:
291 zeros = w2_s[:3]
292 worst_zero = np.abs(zeros).max()
293 if worst_zero > 1e-3:
294 msg = "Translational deviate from 0 significantly: "
295 raise ValueError(msg + "{}".format(w2_s[:3]))
297 w2min = w2_s[3:].min()
298 if w2min < 0:
299 msg = "Dynamical matrix has negative eigenvalues such as "
300 raise ValueError(msg + "{}".format(w2min))
302 # First three modes are translational so ignore:
303 nw = len(w2_s) - 3
304 n_atoms = len(masses)
305 w_s = np.sqrt(w2_s[3:])
306 X_acs = X_is[:, 3:].reshape(n_atoms, 3, nw)
308 # Assign the amplitudes according to Bose-Einstein distribution
309 # or high temperature (== classical) limit
310 if quantum:
311 hbar = units._hbar * units.J * units.s
312 A_s = np.sqrt(hbar * (2 * n_BE(temp, hbar * w_s) + 1) / (2 * w_s))
313 else:
314 A_s = np.sqrt(temp) / w_s
316 if plus_minus:
317 # create samples by multiplying the amplitude with +/-
318 # according to Eq. 5 in PRB 94, 075125
320 spread = (-1) ** np.arange(nw)
322 # gauge eigenvectors: largest value always positive
323 for ii in range(X_acs.shape[-1]):
324 vec = X_acs[:, :, ii]
325 max_arg = np.argmax(abs(vec))
326 X_acs[:, :, ii] *= np.sign(vec.flat[max_arg])
328 # Create velocities und displacements from the amplitudes and
329 # eigenvectors
330 A_s *= spread
331 phi_s = 2.0 * np.pi * rng(nw)
333 # Assign velocities, sqrt(2) compensates for missing sin(phi) in
334 # amplitude for displacement
335 v_ac = (w_s * A_s * np.sqrt(2) * np.cos(phi_s) * X_acs).sum(axis=2)
336 v_ac /= np.sqrt(masses)[:, None]
338 # Assign displacements
339 d_ac = (A_s * X_acs).sum(axis=2)
340 d_ac /= np.sqrt(masses)[:, None]
342 else:
343 # compute the gaussian distribution for the amplitudes
344 # We need 0 < P <= 1.0 and not 0 0 <= P < 1.0 for the logarithm
345 # to avoid (highly improbable) NaN.
347 # Box Muller [en.wikipedia.org/wiki/Box–Muller_transform]:
348 spread = np.sqrt(-2.0 * np.log(1.0 - rng(nw)))
350 # assign amplitudes and phases
351 A_s *= spread
352 phi_s = 2.0 * np.pi * rng(nw)
354 # Assign velocities and displacements
355 v_ac = (w_s * A_s * np.cos(phi_s) * X_acs).sum(axis=2)
356 v_ac /= np.sqrt(masses)[:, None]
358 d_ac = (A_s * np.sin(phi_s) * X_acs).sum(axis=2)
359 d_ac /= np.sqrt(masses)[:, None]
361 if return_eigensolution:
362 return d_ac, v_ac, w2_s, X_is
363 # else
364 return d_ac, v_ac
367def PhononHarmonics(
368 atoms,
369 force_constants,
370 temp=None,
371 *,
372 temperature_K=None,
373 rng=np.random,
374 quantum=False,
375 plus_minus=False,
376 return_eigensolution=False,
377 failfast=True,
378):
379 r"""Excite phonon modes to specified temperature.
381 This will displace atomic positions and set the velocities so as
382 to produce a random, phononically correct state with the requested
383 temperature.
385 Parameters:
387 atoms: ase.atoms.Atoms() object
388 Positions and momenta of this object are perturbed.
390 force_constants: ndarray of size 3N x 3N
391 Force constants for the the structure represented by atoms in eV/Ų
393 temp: float (deprecated).
394 Temperature in eV. Deprecated, use ``temperature_K`` instead.
396 temperature_K: float
397 Temperature in Kelvin.
399 rng: Random number generator
400 RandomState or other random number generator, e.g., np.random.rand
402 quantum: bool
403 True for Bose-Einstein distribution, False for Maxwell-Boltzmann
404 (classical limit)
406 failfast: bool
407 True for sanity checking the phonon spectrum for negative frequencies
408 at Gamma.
410 """
412 # Receive displacements and velocities from phonon_harmonics()
413 d_ac, v_ac = phonon_harmonics(
414 force_constants=force_constants,
415 masses=atoms.get_masses(),
416 temp=temp,
417 temperature_K=temperature_K,
418 rng=rng.rand,
419 plus_minus=plus_minus,
420 quantum=quantum,
421 failfast=failfast,
422 return_eigensolution=False,
423 )
425 # Assign new positions (with displacements) and velocities
426 atoms.positions += d_ac
427 atoms.set_velocities(v_ac)