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':

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

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

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)