Coverage for /builds/kinetik161/ase/ase/md/md.py: 82.69%
52 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1"""Molecular Dynamics."""
2import warnings
3from typing import IO, Optional, Union
5import numpy as np
7from ase import Atoms, units
8from ase.io.trajectory import Trajectory
9from ase.md.logger import MDLogger
10from ase.optimize.optimize import Dynamics
13def process_temperature(
14 temperature: Optional[float],
15 temperature_K: Optional[float],
16 orig_unit: str,
17) -> float:
18 """Handle that temperature can be specified in multiple units.
20 For at least a transition period, molecular dynamics in ASE can
21 have the temperature specified in either Kelvin or Electron
22 Volt. The different MD algorithms had different defaults, by
23 forcing the user to explicitly choose a unit we can resolve
24 this. Using the original method then will issue a
25 FutureWarning.
27 Four parameters:
29 temperature: None or float
30 The original temperature specification in whatever unit was
31 historically used. A warning is issued if this is not None and
32 the historical unit was eV.
34 temperature_K: None or float
35 Temperature in Kelvin.
37 orig_unit: str
38 Unit used for the `temperature`` parameter. Must be 'K' or 'eV'.
40 Exactly one of the two temperature parameters must be different from
41 None, otherwise an error is issued.
43 Return value: Temperature in Kelvin.
44 """
45 if (temperature is not None) + (temperature_K is not None) != 1:
46 raise TypeError("Exactly one of the parameters 'temperature',"
47 + " and 'temperature_K', must be given")
48 if temperature is not None:
49 w = "Specify the temperature in K using the 'temperature_K' argument"
50 if orig_unit == 'K':
51 return temperature
52 elif orig_unit == 'eV':
53 warnings.warn(FutureWarning(w))
54 return temperature / units.kB
55 else:
56 raise ValueError("Unknown temperature unit " + orig_unit)
58 assert temperature_K is not None
59 return temperature_K
62class MolecularDynamics(Dynamics):
63 """Base-class for all MD classes."""
65 def __init__(
66 self,
67 atoms: Atoms,
68 timestep: float,
69 trajectory: Optional[str] = None,
70 logfile: Optional[Union[IO, str]] = None,
71 loginterval: int = 1,
72 append_trajectory: bool = False,
73 ):
74 """Molecular Dynamics object.
76 Parameters:
78 atoms: Atoms object
79 The Atoms object to operate on.
81 timestep: float
82 The time step in ASE time units.
84 trajectory: Trajectory object or str
85 Attach trajectory object. If *trajectory* is a string a
86 Trajectory will be constructed. Use *None* for no
87 trajectory.
89 logfile: file object or str (optional)
90 If *logfile* is a string, a file with that name will be opened.
91 Use '-' for stdout.
93 loginterval: int (optional)
94 Only write a log line for every *loginterval* time steps.
95 Default: 1
97 append_trajectory: boolean (optional)
98 Defaults to False, which causes the trajectory file to be
99 overwriten each time the dynamics is restarted from scratch.
100 If True, the new structures are appended to the trajectory
101 file instead.
102 """
103 # dt as to be attached _before_ parent class is initialized
104 self.dt = timestep
106 super().__init__(atoms, logfile=None, trajectory=None)
108 # Some codes (e.g. Asap) may be using filters to
109 # constrain atoms or do other things. Current state of the art
110 # is that "atoms" must be either Atoms or Filter in order to
111 # work with dynamics.
112 #
113 # In the future, we should either use a special role interface
114 # for MD, or we should ensure that the input is *always* a Filter.
115 # That way we won't need to test multiple cases. Currently,
116 # we do not test /any/ kind of MD with any kind of Filter in ASE.
117 self.atoms = atoms
118 self.masses = self.atoms.get_masses()
120 if 0 in self.masses:
121 warnings.warn('Zero mass encountered in atoms; this will '
122 'likely lead to errors if the massless atoms '
123 'are unconstrained.')
125 self.masses.shape = (-1, 1)
127 if not self.atoms.has('momenta'):
128 self.atoms.set_momenta(np.zeros([len(self.atoms), 3]))
130 # Trajectory is attached here instead of in Dynamics.__init__
131 # to respect the loginterval argument.
132 if trajectory is not None:
133 if isinstance(trajectory, str):
134 mode = "a" if append_trajectory else "w"
135 trajectory = self.closelater(
136 Trajectory(trajectory, mode=mode, atoms=atoms)
137 )
138 self.attach(trajectory, interval=loginterval)
140 if logfile:
141 logger = self.closelater(
142 MDLogger(dyn=self, atoms=atoms, logfile=logfile))
143 self.attach(logger, loginterval)
145 def todict(self):
146 return {'type': 'molecular-dynamics',
147 'md-type': self.__class__.__name__,
148 'timestep': self.dt}
150 def irun(self, steps=50):
151 """Run molecular dynamics algorithm as a generator.
153 Parameters
154 ----------
155 steps : int, default=DEFAULT_MAX_STEPS
156 Number of molecular dynamics steps to be run.
158 Yields
159 ------
160 converged : bool
161 True if the maximum number of steps are reached.
162 """
163 return Dynamics.irun(self, steps=steps)
165 def run(self, steps=50):
166 """Run molecular dynamics algorithm.
168 Parameters
169 ----------
170 steps : int, default=DEFAULT_MAX_STEPS
171 Number of molecular dynamics steps to be run.
173 Returns
174 -------
175 converged : bool
176 True if the maximum number of steps are reached.
177 """
178 return Dynamics.run(self, steps=steps)
180 def get_time(self):
181 return self.nsteps * self.dt
183 def converged(self):
184 """ MD is 'converged' when number of maximum steps is reached. """
185 return self.nsteps >= self.max_steps
187 def _get_com_velocity(self, velocity):
188 """Return the center of mass velocity.
189 Internal use only. This function can be reimplemented by Asap.
190 """
191 return np.dot(self.masses.ravel(), velocity) / self.masses.sum()
193 # Make the process_temperature function available to subclasses
194 # as a static method. This makes it easy for MD objects to use
195 # it, while functions in md.velocitydistribution have access to it
196 # as a function.
197 _process_temperature = staticmethod(process_temperature)