/usr/share/pyshared/ase/md/nptberendsen.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 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 | """Berendsen NPT dynamics class."""
import numpy as np
#from ase.md import MolecularDynamics
from ase.md.nvtberendsen import NVTBerendsen
import ase.units as units
#import math
class NPTBerendsen(NVTBerendsen):
"""Berendsen (constant N, P, T) molecular dynamics.
This dynamics scale the velocities and volumes to maintain a constant
pressure and temperature. The shape of the simulation cell is not
altered, if that is desired use Inhomogenous_NPTBerendsen.
Usage: NPTBerendsen(atoms, timestep, temperature, taut, pressure, taup)
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.
pressure
The desired pressure, in bar (1 bar = 1e5 Pa).
taup
Time constant for Berendsen pressure coupling.
compressibility
The compressibility of the material, water 4.57E-5 bar-1, in bar-1
"""
def __init__(self, atoms, timestep, temperature, taut=0.5e3*units.fs,
pressure = 1.01325, taup=1e3*units.fs,
compressibility=4.57e-5, fixcm=True,
trajectory=None, logfile=None, loginterval=1):
NVTBerendsen.__init__(self, atoms, timestep, temperature, taut, fixcm,
trajectory,
logfile, loginterval)
self.taup = taup
self.pressure = pressure
self.compressibility = compressibility
def set_taup(self, taut):
self.taut = taut
def get_taup(self):
return self.taut
def set_pressure(self, pressure):
self.pressure = pressure
def get_pressure(self):
return self.pressure
def set_compressibility(self, compressibility):
self.compressibility = compressibility
def get_compressibility(self):
return self.compressibility
def set_timestep(self, timestep):
self.dt = timestep
def get_timestep(self):
return self.dt
def scale_positions_and_cell(self):
""" Do the Berendsen pressure coupling,
scale the atom positon and the simulation cell."""
taupscl = self.dt / self.taup
stress = self.atoms.get_stress()
old_pressure = self.atoms.get_isotropic_pressure(stress)
scl_pressure = 1.0 - taupscl * self.compressibility / 3.0 * \
(self.pressure - old_pressure)
#print "old_pressure", old_pressure
#print "volume scaling by:", scl_pressure
cell = self.atoms.get_cell()
cell = scl_pressure * cell
self.atoms.set_cell(cell, scale_atoms=True)
def step(self, f):
""" move one timestep forward using Berenden NPT molecular dynamics."""
NVTBerendsen.scale_velocities(self)
self.scale_positions_and_cell()
#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
class Inhomogenous_NPTBerendsen(NPTBerendsen):
"""Berendsen (constant N, P, T) molecular dynamics.
This dynamics scale the velocities and volumes to maintain a constant
pressure and temperature. The size of the unit cell is allowed to change
independently in the three directions, but the angles remain constant.
Usage: NPTBerendsen(atoms, timestep, temperature, taut, pressure, taup)
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.
pressure
The desired pressure, in bar (1 bar = 1e5 Pa).
taup
Time constant for Berendsen pressure coupling.
compressibility
The compressibility of the material, water 4.57E-5 bar-1, in bar-1
"""
def scale_positions_and_cell(self):
""" Do the Berendsen pressure coupling,
scale the atom positon and the simulation cell."""
taupscl = self.dt * self.compressibility / self.taup / 3.0
stress = - self.atoms.get_stress() * 1e-5 / units.Pascal
if stress.shape == (6,):
stress = stress[:3]
elif stress.shape == (3,3):
stress = [stress[i][i] for i in range(3)]
else:
raise ValueError("Cannot use a stress tensor of shape " + str(stress.shape))
scl_pressurex = 1.0 - taupscl * (self.pressure - stress[0])
scl_pressurey = 1.0 - taupscl * (self.pressure - stress[1])
scl_pressurez = 1.0 - taupscl * (self.pressure - stress[2])
cell = self.atoms.get_cell()
cell = np.array([scl_pressurex * cell[0],scl_pressurey * cell[1],scl_pressurez * cell[2]])
self.atoms.set_cell(cell, scale_atoms=True)
|