/usr/share/pyshared/ase/md/nvtberendsen.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 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 | """Berendsen NVT dynamics class."""
import sys
import numpy as np
from ase.md.md import MolecularDynamics
from ase.parallel import world
class NVTBerendsen(MolecularDynamics):
"""Berendsen (constant N, V, T) molecular dynamics.
Usage: NVTBerendsen(atoms, timestep, temperature, taut, fixcm)
atoms
The list of atoms.
timestep
The time step.
temperature
The desired temperature, in Kelvin.
taut
Time constant for Berendsen temperature coupling.
fixcm
If True, the position and momentum of the center of mass is
kept unperturbed. Default: True.
"""
def __init__(self, atoms, timestep, temperature, taut, fixcm=True,
trajectory=None, logfile=None, loginterval=1,
communicator=world):
MolecularDynamics.__init__(self, atoms, timestep, trajectory,
logfile, loginterval)
self.taut = taut
self.temperature = temperature
self.fixcm = fixcm # will the center of mass be held fixed?
self.communicator = communicator
def set_taut(self, taut):
self.taut = taut
def get_taut(self):
return self.taut
def set_temperature(self, temperature):
self.temperature = temperature
def get_temperature(self):
return self.temperature
def set_timestep(self, timestep):
self.dt = timestep
def get_timestep(self):
return self.dt
def scale_velocities(self):
""" Do the NVT Berendsen velocity scaling """
tautscl = self.dt / self.taut
old_temperature = self.atoms.get_temperature()
scl_temperature = np.sqrt(1.0+ (self.temperature/ old_temperature- 1.0)
*tautscl)
#limit the velocity scaling to reasonable values
if scl_temperature > 1.1:
scl_temperature = 1.1
if scl_temperature < 0.9:
scl_temperature = 0.9
atoms = self.atoms
p = self.atoms.get_momenta()
p = scl_temperature * p
self.atoms.set_momenta(p)
return
def step(self, f):
""" move one timestep forward using Berenden NVT molecular dynamics."""
self.scale_velocities()
#one step velocity verlet
atoms = self.atoms
p = self.atoms.get_momenta()
p += 0.5 * self.dt * f
if self.fixcm:
# calculate the center of mass
# momentum and subtract it
psum = p.sum(axis=0) / float(len(p))
p = p - psum
self.atoms.set_positions(self.atoms.get_positions() +
self.dt * p / self.atoms.get_masses()[:,np.newaxis])
# We need to store the momenta on the atoms before calculating
# the forces, as in a parallel Asap calculation atoms may
# migrate during force calculations, and the momenta need to
# migrate along with the atoms. For the same reason, we
# cannot use self.masses in the line above.
self.atoms.set_momenta(p)
f = self.atoms.get_forces()
atoms.set_momenta(self.atoms.get_momenta() + 0.5 * self.dt * f)
return f
|