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
|
"""Molecular Dynamics."""
import warnings
import numpy as np
from ase.optimize.optimize import Dynamics
from ase.md.logger import MDLogger
from ase.io.trajectory import Trajectory
from ase import units
def process_temperature(temperature, temperature_K, orig_unit):
"""Handle that temperature can be specified in multiple units.
For at least a transition period, molecular dynamics in ASE can
have the temperature specified in either Kelvin or Electron
Volt. The different MD algorithms had different defaults, by
forcing the user to explicitly choose a unit we can resolve
this. Using the original method then will issue a
FutureWarning.
Four parameters:
temperature: None or float
The original temperature specification in whatever unit was
historically used. A warning is issued if this is not None and
the historical unit was eV.
temperature_K: None or float
Temperature in Kelvin.
orig_unit: str
Unit used for the `temperature`` parameter. Must be 'K' or 'eV'.
Exactly one of the two temperature parameters must be different from
None, otherwise an error is issued.
Return value: Temperature in Kelvin.
"""
if (temperature is not None) + (temperature_K is not None) != 1:
raise TypeError("Exactly one of the parameters 'temperature',"
+ " and 'temperature_K', must be given")
if temperature is not None:
w = "Specify the temperature in K using the 'temperature_K' argument"
if orig_unit == 'K':
return temperature
elif orig_unit == 'eV':
warnings.warn(FutureWarning(w))
return temperature / units.kB
else:
raise ValueError("Unknown temperature unit " + orig_unit)
assert temperature_K is not None
return temperature_K
class MolecularDynamics(Dynamics):
"""Base-class for all MD classes."""
def __init__(self, atoms, timestep, trajectory, logfile=None,
loginterval=1, append_trajectory=False):
"""Molecular Dynamics object.
Parameters:
atoms: Atoms object
The Atoms object to operate on.
timestep: float
The time step in ASE time units.
trajectory: Trajectory object or str
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.
"""
# dt as to be attached _before_ parent class is initialized
self.dt = timestep
Dynamics.__init__(self, atoms, logfile=None, trajectory=None)
self.masses = self.atoms.get_masses()
self.max_steps = None
if 0 in self.masses:
warnings.warn('Zero mass encountered in atoms; this will '
'likely lead to errors if the massless atoms '
'are unconstrained.')
self.masses.shape = (-1, 1)
if not self.atoms.has('momenta'):
self.atoms.set_momenta(np.zeros([len(self.atoms), 3]))
# Trajectory is attached here instead of in Dynamics.__init__
# to respect the loginterval argument.
try:
if trajectory is not None:
if isinstance(trajectory, str):
mode = "a" if append_trajectory else "w"
trajectory = self.ensureclose(
Trajectory(trajectory, mode=mode, atoms=atoms)
)
self.attach(trajectory, interval=loginterval)
if logfile:
logger = self.ensureclose(
MDLogger(dyn=self, atoms=atoms, logfile=logfile))
self.attach(logger, loginterval)
except BaseException:
self._closefiles()
raise
def todict(self):
return {'type': 'molecular-dynamics',
'md-type': self.__class__.__name__,
'timestep': self.dt}
def irun(self, steps=50):
""" Call Dynamics.irun and adjust max_steps """
self.max_steps = steps + self.nsteps
return Dynamics.irun(self)
def run(self, steps=50):
""" Call Dynamics.run and adjust max_steps """
self.max_steps = steps + self.nsteps
return Dynamics.run(self)
def get_time(self):
return self.nsteps * self.dt
def converged(self):
""" MD is 'converged' when number of maximum steps is reached. """
return self.nsteps >= self.max_steps
def _get_com_velocity(self, velocity):
"""Return the center of mass velocity.
Internal use only. This function can be reimplemented by Asap.
"""
return np.dot(self.masses.ravel(), velocity) / self.masses.sum()
# Make the process_temperature function available to subclasses
# as a static method. This makes it easy for MD objects to use
# it, while functions in md.velocitydistribution have access to it
# as a function.
_process_temperature = staticmethod(process_temperature)
|
