/usr/share/pyshared/ase/md/velocitydistribution.py is in python-ase 3.6.0.2515-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 | # VelocityDistributions.py -- set up a velocity distribution
"""Module for setting up e.g. Maxwell-Boltzmann velocity distributions.
Currently, only one function is defined, MaxwellBoltzmannDistribution,
which sets the momenta of a list of atoms according to a
Maxwell-Boltzmann distribution at a given temperature.
"""
import sys
import numpy as np
from ase.parallel import world
def _maxwellboltzmanndistribution(masses, temp, communicator=world):
# For parallel GPAW simulations, the random velocities should be
# distributed. Uses gpaw world communicator as default, but allow
# option of specifying other communicator (for ensemble runs)
xi = np.random.standard_normal((len(masses), 3))
communicator.broadcast(xi, 0)
momenta = xi * np.sqrt(masses * temp)[:, np.newaxis]
return momenta
def MaxwellBoltzmannDistribution(atoms, temp, communicator=world):
"""Sets the momenta to a Maxwell-Boltzmann distribution."""
momenta = _maxwellboltzmanndistribution(atoms.get_masses(), temp,
communicator)
atoms.set_momenta(momenta)
def Stationary(atoms):
"Sets the center-of-mass momentum to zero."
p = atoms.get_momenta()
p0 = np.sum(p, 0)
# We should add a constant velocity, not momentum, to the atoms
m = atoms.get_masses()
mtot = np.sum(m)
v0 = p0 / mtot
p -= v0 * m[:, np.newaxis]
atoms.set_momenta(p)
def ZeroRotation(atoms):
"Sets the total angular momentum to zero by counteracting rigid rotations."
# Find the principal moments of inertia and principal axes basis vectors
Ip, basis = atoms.get_moments_of_inertia(vectors=True)
# Calculate the total angular momentum and transform to principal basis
Lp = np.dot(basis, atoms.get_angular_momentum())
# Calculate the rotation velocity vector in the principal basis, avoiding
# zero division, and transform it back to the cartesian coordinate system
omega = np.dot(np.linalg.inv(basis), np.select([Ip > 0], [Lp / Ip]))
# We subtract a rigid rotation corresponding to this rotation vector
com = atoms.get_center_of_mass()
positions = atoms.get_positions()
positions -= com # translate center of mass to origin
velocities = atoms.get_velocities()
atoms.set_velocities(velocities - np.cross(omega, positions))
|