Hide keyboard shortcuts

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 

2 

3"""Module for setting up velocity distributions such as Maxwell–Boltzmann. 

4 

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. 

9 

10""" 

11 

12import numpy as np 

13from ase.parallel import world 

14from ase import units 

15from ase.md.md import process_temperature 

16 

17# define a ``zero'' temperature to avoid divisions by zero 

18eps_temp = 1e-12 

19 

20 

21class UnitError(Exception): 

22 """Exception raised when wrong units are specified""" 

23 

24 

25def force_temperature(atoms, temperature, unit="K"): 

26 """ force (nucl.) temperature to have a precise value 

27 

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 """ 

36 

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)) 

43 

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)) 

50 

51 

52def _maxwellboltzmanndistribution(masses, temp, communicator=None, rng=None): 

53 """Return a Maxwell-Boltzmann distribution with a given temperature. 

54 

55 Paremeters: 

56 

57 masses: float 

58 The atomic masses. 

59 

60 temp: float 

61 The temperature in electron volt. 

62 

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 

67 

68 rng: numpy RNG (optional) 

69 The random number generator. Default: np.random 

70 

71 Returns: 

72 

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 

84 

85 

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. 

92 

93 Parameters: 

94 

95 atoms: Atoms object 

96 The atoms. Their momenta will be modified. 

97 

98 temp: float (deprecated) 

99 The temperature in eV. Deprecated, used temperature_K instead. 

100 

101 temperature_K: float 

102 The temperature in Kelvin. 

103 

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 

108 

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. 

113 

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") 

123 

124 

125def Stationary(atoms, preserve_temperature=True): 

126 "Sets the center-of-mass momentum to zero." 

127 

128 # Save initial temperature 

129 temp0 = atoms.get_temperature() 

130 

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) 

139 

140 if preserve_temperature: 

141 force_temperature(atoms, temp0) 

142 

143 

144def ZeroRotation(atoms, preserve_temperature=True): 

145 "Sets the total angular momentum to zero by counteracting rigid rotations." 

146 

147 # Save initial temperature 

148 temp0 = atoms.get_temperature() 

149 

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)) 

163 

164 if preserve_temperature: 

165 force_temperature(atoms, temp0) 

166 

167 

168def n_BE(temp, omega): 

169 """Bose-Einstein distribution function. 

170 

171 Args: 

172 temp: temperature converted to eV (*units.kB) 

173 omega: sequence of frequencies converted to eV 

174 

175 Returns: 

176 Value of Bose-Einstein distribution function for each energy 

177 

178 """ 

179 

180 omega = np.asarray(omega) 

181 

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 

188 

189 

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. 

203 

204 Parameters: 

205 

206 force_constants: array of size 3N x 3N 

207 force constants (Hessian) of the system in eV/Ų 

208 

209 masses: array of length N 

210 masses of the structure in amu 

211 

212 temp: float (deprecated) 

213 Temperature converted to eV (T * units.kB). Deprecated, use 

214 ``temperature_K``. 

215 

216 temperature_K: float 

217 Temperature in Kelvin. 

218 

219 rng: function 

220 Random number generator function, e.g., np.random.rand 

221 

222 quantum: bool 

223 True for Bose-Einstein distribution, False for Maxwell-Boltzmann 

224 (classical limit) 

225 

226 plus_minus: bool 

227 Displace atoms with +/- the amplitude accoding to PRB 94, 075125 

228 

229 return_eigensolution: bool 

230 return eigenvalues and eigenvectors of the dynamical matrix 

231 

232 failfast: bool 

233 True for sanity checking the phonon spectrum for negative 

234 frequencies at Gamma 

235 

236 Returns: 

237 

238 Displacements, velocities generated from the eigenmodes, 

239 (optional: eigenvalues, eigenvectors of dynamical matrix) 

240 

241 Purpose: 

242 

243 Excite phonon modes to specified temperature. 

244 

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. 

249 

250 In other words the system will be equilibrated for an MD run at 

251 that temperature. 

252 

253 force_constants should be the matrix as force constants, e.g., 

254 as computed by the ase.phonons module. 

255 

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 

259 

260 .. code-block:: none 

261 

262 

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 

268 

269 

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 

275 

276 Reference: [West, Estreicher; PRL 96, 22 (2006)] 

277 """ 

278 

279 # Handle the temperature units 

280 temp = units.kB * process_temperature(temp, temperature_K, 'eV') 

281 

282 # Build dynamical matrix 

283 rminv = (masses ** -0.5).repeat(3) 

284 dynamical_matrix = force_constants * rminv[:, None] * rminv[None, :] 

285 

286 # Solve eigenvalue problem to compute phonon spectrum and eigenvectors 

287 w2_s, X_is = np.linalg.eigh(dynamical_matrix) 

288 

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])) 

296 

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)) 

301 

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) 

307 

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 

315 

316 if plus_minus: 

317 # create samples by multiplying the amplitude with +/- 

318 # according to Eq. 5 in PRB 94, 075125 

319 

320 spread = (-1) ** np.arange(nw) 

321 

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]) 

327 

328 # Create velocities und displacements from the amplitudes and 

329 # eigenvectors 

330 A_s *= spread 

331 phi_s = 2.0 * np.pi * rng(nw) 

332 

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] 

337 

338 # Assign displacements 

339 d_ac = (A_s * X_acs).sum(axis=2) 

340 d_ac /= np.sqrt(masses)[:, None] 

341 

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. 

346 

347 # Box Muller [en.wikipedia.org/wiki/Box–Muller_transform]: 

348 spread = np.sqrt(-2.0 * np.log(1.0 - rng(nw))) 

349 

350 # assign amplitudes and phases 

351 A_s *= spread 

352 phi_s = 2.0 * np.pi * rng(nw) 

353 

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] 

357 

358 d_ac = (A_s * np.sin(phi_s) * X_acs).sum(axis=2) 

359 d_ac /= np.sqrt(masses)[:, None] 

360 

361 if return_eigensolution: 

362 return d_ac, v_ac, w2_s, X_is 

363 # else 

364 return d_ac, v_ac 

365 

366 

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. 

380 

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. 

384 

385 Parameters: 

386 

387 atoms: ase.atoms.Atoms() object 

388 Positions and momenta of this object are perturbed. 

389 

390 force_constants: ndarray of size 3N x 3N 

391 Force constants for the the structure represented by atoms in eV/Ų 

392 

393 temp: float (deprecated). 

394 Temperature in eV. Deprecated, use ``temperature_K`` instead. 

395 

396 temperature_K: float 

397 Temperature in Kelvin. 

398 

399 rng: Random number generator 

400 RandomState or other random number generator, e.g., np.random.rand 

401 

402 quantum: bool 

403 True for Bose-Einstein distribution, False for Maxwell-Boltzmann 

404 (classical limit) 

405 

406 failfast: bool 

407 True for sanity checking the phonon spectrum for negative frequencies 

408 at Gamma. 

409 

410 """ 

411 

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 ) 

424 

425 # Assign new positions (with displacements) and velocities 

426 atoms.positions += d_ac 

427 atoms.set_velocities(v_ac)