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
|
"""Andersen dynamics class."""
from numpy import random, cos, pi, log, ones, repeat
from ase.md.md import MolecularDynamics
from ase.parallel import world, DummyMPI
from ase import units
class Andersen(MolecularDynamics):
"""Andersen (constant N, V, T) molecular dynamics."""
def __init__(self, atoms, timestep, temperature_K, andersen_prob,
fixcm=True, trajectory=None, logfile=None, loginterval=1,
communicator=world, rng=random, append_trajectory=False):
""""
Parameters:
atoms: Atoms object
The list of atoms.
timestep: float
The time step in ASE time units.
temperature_K: float
The desired temperature, in Kelvin.
andersen_prob: float
A random collision probability, typically 1e-4 to 1e-1.
With this probability atoms get assigned random velocity components.
fixcm: bool (optional)
If True, the position and momentum of the center of mass is
kept unperturbed. Default: True.
rng: RNG object (optional)
Random number generator, by default numpy.random. Must have a
random_sample method matching the signature of
numpy.random.random_sample.
logfile: file object or str (optional)
If *logfile* is a string, a file with that name will be opened.
Use '-' for stdout.
trajectory: Trajectory object or str (optional)
Attach trajectory object. If *trajectory* is a string a
Trajectory will be constructed. Use *None* (the default) for no
trajectory.
communicator: MPI communicator (optional)
Communicator used to distribute random numbers to all tasks.
Default: ase.parallel.world. Set to None to disable communication.
append_trajectory: bool (optional)
Defaults to False, which causes the trajectory file to be
overwritten each time the dynamics is restarted from scratch.
If True, the new structures are appended to the trajectory
file instead.
The temperature is imposed by stochastic collisions with a heat bath
that acts on velocity components of randomly chosen particles.
The algorithm randomly decorrelates velocities, so dynamical properties
like diffusion or viscosity cannot be properly measured.
H. C. Andersen, J. Chem. Phys. 72 (4), 2384–2393 (1980)
"""
self.temp = units.kB * temperature_K
self.andersen_prob = andersen_prob
self.fix_com = fixcm
self.rng = rng
if communicator is None:
communicator = DummyMPI()
self.communicator = communicator
MolecularDynamics.__init__(self, atoms, timestep, trajectory,
logfile, loginterval,
append_trajectory=append_trajectory)
def set_temperature(self, temperature_K):
self.temp = units.kB * temperature_K
def set_andersen_prob(self, andersen_prob):
self.andersen_prob = andersen_prob
def set_timestep(self, timestep):
self.dt = timestep
def boltzmann_random(self, width, size):
x = self.rng.random_sample(size=size)
y = self.rng.random_sample(size=size)
z = width * cos(2 * pi * x) * (-2 * log(1 - y))**0.5
return z
def get_maxwell_boltzmann_velocities(self):
natoms = len(self.atoms)
masses = repeat(self.masses, 3).reshape(natoms, 3)
width = (self.temp / masses)**0.5
velos = self.boltzmann_random(width, size=(natoms, 3))
return velos # [[x, y, z],] components for each atom
def step(self, forces=None):
atoms = self.atoms
if forces is None:
forces = atoms.get_forces()
self.v = atoms.get_velocities()
# Random atom-wise variables are stored as attributes and broadcasted:
# - self.random_com_velocity # added to all atoms if self.fix_com
# - self.random_velocity # added to some atoms if the per-atom
# - self.andersen_chance # andersen_chance <= andersen_prob
# a dummy communicator will be used for serial runs
if self.fix_com:
# add random velocity to center of mass to prepare Andersen
width = (self.temp / sum(self.masses))**0.5
self.random_com_velocity = (ones(self.v.shape)
* self.boltzmann_random(width, (3)))
self.communicator.broadcast(self.random_com_velocity, 0)
self.v += self.random_com_velocity
self.v += 0.5 * forces / self.masses * self.dt
# apply Andersen thermostat
self.random_velocity = self.get_maxwell_boltzmann_velocities()
self.andersen_chance = self.rng.random_sample(size=self.v.shape)
self.communicator.broadcast(self.random_velocity, 0)
self.communicator.broadcast(self.andersen_chance, 0)
self.v[self.andersen_chance <= self.andersen_prob] \
= self.random_velocity[self.andersen_chance <= self.andersen_prob]
x = atoms.get_positions()
if self.fix_com:
old_com = atoms.get_center_of_mass()
self.v -= self._get_com_velocity(self.v)
# Step: x^n -> x^(n+1) - this applies constraints if any
atoms.set_positions(x + self.v * self.dt)
if self.fix_com:
atoms.set_center_of_mass(old_com)
# recalc velocities after RATTLE constraints are applied
self.v = (atoms.get_positions() - x) / self.dt
forces = atoms.get_forces(md=True)
# Update the velocities
self.v += 0.5 * forces / self.masses * self.dt
if self.fix_com:
self.v -= self._get_com_velocity(self.v)
# Second part of RATTLE taken care of here
atoms.set_momenta(self.v * self.masses)
return forces
|
