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
|
"""Berendsen NVT dynamics class."""
import numpy as np
from ase.md.md import MolecularDynamics
from ase.parallel import world
class NVTBerendsen(MolecularDynamics):
def __init__(self, atoms, timestep, temperature=None, taut=None,
fixcm=True, *, temperature_K=None,
trajectory=None, logfile=None, loginterval=1,
communicator=world, append_trajectory=False):
"""Berendsen (constant N, V, T) molecular dynamics.
Parameters:
atoms: Atoms object
The list of atoms.
timestep: float
The time step in ASE time units.
temperature: float
The desired temperature, in Kelvin.
temperature_K: float
Alias for *temperature*
taut: float
Time constant for Berendsen temperature coupling in ASE
time units.
fixcm: bool (optional)
If True, the position and momentum of the center of mass is
kept unperturbed. Default: True.
trajectory: Trajectory object or str (optional)
Attach trajectory object. If *trajectory* is a string a
Trajectory will be constructed. Use *None* for no
trajectory.
logfile: file object or str (optional)
If *logfile* is a string, a file with that name will be opened.
Use '-' for stdout.
loginterval: int (optional)
Only write a log line for every *loginterval* time steps.
Default: 1
append_trajectory: boolean (optional)
Defaults to False, which causes the trajectory file to be
overwriten each time the dynamics is restarted from scratch.
If True, the new structures are appended to the trajectory
file instead.
"""
MolecularDynamics.__init__(self, atoms, timestep, trajectory,
logfile, loginterval,
append_trajectory=append_trajectory)
if taut is None:
raise TypeError("Missing 'taut' argument.")
self.taut = taut
self.temperature = self._process_temperature(temperature,
temperature_K, 'K')
self.fix_com = 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=None, *, temperature_K=None):
self.temperature = self._process_temperature(temperature,
temperature_K, 'K')
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
p = self.atoms.get_momenta()
p = scl_temperature * p
self.atoms.set_momenta(p)
return
def step(self, forces=None):
"""Move one timestep forward using Berenden NVT molecular dynamics."""
self.scale_velocities()
# one step velocity verlet
atoms = self.atoms
if forces is None:
forces = atoms.get_forces(md=True)
p = self.atoms.get_momenta()
p += 0.5 * self.dt * forces
if self.fix_com:
# 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)
forces = self.atoms.get_forces(md=True)
atoms.set_momenta(self.atoms.get_momenta() + 0.5 * self.dt * forces)
return forces
|
