/usr/share/python-ase/doc/ase/md.rst is in python-ase-doc 3.12.0-2.
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 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 | ==================
Molecular dynamics
==================
.. module:: ase.md
:synopsis: Molecular Dynamics
Typical computer simulations involve moving the atoms around, either
to optimize a structure (energy minimization) or to do molecular
dynamics. This chapter discusses molecular dynamics, energy
minimization algorithms will be discussed in the :mod:`ase.optimize`
section.
A molecular dynamics object will operate on the atoms by moving them
according to their forces - it integrates Newton's second law
numerically. A typical molecular dynamics simulation will use the
`Velocity Verlet dynamics`_. You create the
:class:`ase.md.verlet.VelocityVerlet` object, giving it the atoms and a time
step, and then you perform dynamics by calling its :meth:`run` method::
dyn = VelocityVerlet(atoms, dt=5.0 * units.fs,
trajectory='md.traj', logfile='md.log')
dyn.run(1000) # take 1000 steps
A number of different algorithms can be used to perform molecular
dynamics, with slightly different results.
Choosing the time step
======================
All the dynamics objects need a time step. Choosing it too small will
waste computer time, choosing it too large will make the dynamics
unstable, typically the energy increases dramatically (the system
"blows up"). If the time step is only a little to large, the lack of
energy conservation is most obvious in `Velocity Verlet dynamics`_,
where energy should otherwise be conserved.
Experience has shown that 5 femtoseconds is a good choice for most metallic
systems. Systems with light atoms (e.g. hydrogen) and/or with strong
bonds (carbon) will need a smaller time step.
All the dynamics objects documented here are sufficiently related to
have the same optimal time step.
File output
===========
The time evolution of the system can be saved in a trajectory file,
by creating a trajectory object, and attaching it to the dynamics
object. This is documented in the module :mod:`ase.io.trajectory`.
You can attach the trajectory explicitly to the dynamics object, and
you may want to use the optional ``interval`` argument, so every
time step is not written to the file.
Alternatively, you can just use the ``trajectory`` keyword when
instantiating the dynamics object as in the example above. In this
case, a ``loginterval`` keyword may also be supplied to specify the
frequency of writing to the trajectory. The loginterval keyword will
apply to both the trajectory and the logfile.
Logging
=======
A logging mechanism is provided, printing time; total, potential and
kinetic energy; and temperature (calculated from the kinetic energy).
It is enabled by giving the ``logfile`` argument when the dynamics
object is created, ``logfile`` may be an open file, a filename or the
string '-' meaning standard output. Per default, a line is printed
for each timestep, specifying the ``loginterval`` argument will chance
this to a more reasonable frequency.
The logging can be customized by explicitly attaching a
:class:`MDLogger` object to the dynamics::
from ase.md import MDLogger
dyn = VelocityVerlet(atoms, dt=2*ase.units.fs)
dyn.attach(MDLogger(dyn, atoms, 'md.log', header=False, stress=False,
peratom=True, mode="a"), interval=1000)
This example will skip the header line and write energies per atom
instead of total energies. The parameters are
``header``: Print a header line defining the columns.
``stress``: Print the six components of the stress tensor.
``peratom``: Print energy per atom instead of total energy.
``mode``: If 'a', append to existing file, if 'w' overwrite
existing file.
Despite appearances, attaching a logger like this does *not* create a
cyclic reference to the dynamics.
.. note::
If building your own logging class, be sure not to attach the dynamics
object directly to the logging object. Instead, create a weak reference
using the ``proxy`` method of the ``weakref`` package. See the
*ase.md.MDLogger* source code for an example. (If this is not done, a
cyclic reference may be created which can cause certain calculators,
such as Jacapo, to not terminate correctly.)
Constant NVE simulations (the microcanonical ensemble)
======================================================
Newton's second law preserves the total energy of the system, and a
straightforward integration of Newton's second law therefore leads to
simulations preserving the total energy of the system (E), the number
of atoms (N) and the volume of the system (V). The most appropriate
algorithm for doing this is velocity Verlet dynamics, since it gives
very good long-term stability of the total energy even with quite
large time steps. Fancier algorithms such as Runge-Kutta may give
very good short-term energy preservation, but at the price of a slow
drift in energy over longer timescales, causing trouble for long
simulations.
In a typical NVE simulation, the temperature will remain approximately
constant, but if significant structural changes occurs they may result
in temperature changes. If external work is done on the system, the
temperature is likely to rise significantly.
Velocity Verlet dynamics
------------------------
.. module:: ase.md.verlet
.. class:: VelocityVerlet(atoms, timestep)
``VelocityVerlet`` is the only dynamics implementing the NVE ensemble.
It requires two arguments, the atoms and the time step. Choosing
a too large time step will immediately be obvious, as the energy will
increase with time, often very rapidly.
Example: See the tutorial :ref:`md_tutorial`.
Constant NVT simulations (the canonical ensemble)
=================================================
Since Newton's second law conserves energy and not temperature,
simulations at constant temperature will somehow involve coupling the
system to a heat bath. This cannot help being somewhat artificial.
Two different approaches are possible within ASE. In Langevin
dynamics, each atom is coupled to a heat bath through a fluctuating
force and a friction term. In Nosé-Hoover dynamics, a term
representing the heat bath through a single degree of freedom is
introduced into the Hamiltonian.
Langevin dynamics
-----------------
.. module:: ase.md.langevin
.. class:: Langevin(atoms, timestep, temperature, friction)
The Langevin class implements Langevin dynamics, where a (small)
friction term and a fluctuating force are added to Newton's second law
which is then integrated numerically. The temperature of the heat
bath and magnitude of the friction is specified by the user, the
amplitude of the fluctuating force is then calculated to give that
temperature. This procedure has some physical justification: in a
real metal the atoms are (weakly) coupled to the electron gas, and the
electron gas therefore acts like a heat bath for the atoms. If heat
is produced locally, the atoms locally get a temperature that is
higher than the temperature of the electrons, heat is transferred to
the electrons and then rapidly transported away by them. A Langevin
equation is probably a reasonable model for this process.
A disadvantage of using Langevin dynamics is that if significant heat
is produced in the simulation, then the temperature will stabilize at
a value higher than the specified temperature of the heat bath, since
a temperature difference between the system and the heat bath is
necessary to get a finite heat flow. Another disadvantage is that the
fluctuating force is stochastic in nature, so repeating the simulation
will not give exactly the same trajectory.
When the ``Langevin`` object is created, you must specify a time step,
a temperature (in energy units) and a friction. Typical values for
the friction are 0.01-0.02 atomic units.
::
# Room temperature simulation
dyn = Langevin(atoms, 5 * units.fs, units.kB * 300, 0.002)
Both the friction and the temperature can be replaced with arrays
giving per-atom values. This is mostly useful for the friction, where
one can choose a rather high friction near the boundaries, and set it
to zero in the part of the system where the phenomenon being studied
is located.
Nosé-Hoover dynamics
--------------------
In Nosé-Hoover dynamics, an extra term is added to the Hamiltonian
representing the coupling to the heat bath. From a pragmatic point of
view one can regard Nosé-Hoover dynamics as adding a friction term to
Newton's second law, but dynamically changing the friction coefficient
to move the system towards the desired temperature. Typically the
"friction coefficient" will fluctuate around zero.
Nosé-Hoover dynamics is not implemented as a separate class, but is a
special case of NPT dynamics.
Berendsen NVT dynamics
-----------------------
.. module:: ase.md.nvtberendsen
.. class:: NVTBerendsen(atoms, timestep, temperature, taut, fixcm)
In Berendsen NVT simulations the velocities are scaled to achieve the desired
temperature. The speed of the scaling is determined by the parameter taut.
This method does not result proper NVT sampling but it usually is
sufficiently good in practise (with large taut). For discussion see
the gromacs manual at www.gromacs.org.
*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.
::
# Room temperature simulation (300K, 0.1 fs time step)
dyn = NVTBerendsen(atoms, 0.1 * units.fs, 300, taut=0.5*1000*units.fs)
Constant NPT simulations (the isothermal-isobaric ensemble)
===========================================================
.. module:: ase.md.npt
.. class:: NPT(atoms, timestep, temperature, externalstress, ttime, pfactor, mask=None)
Dynamics with constant pressure (or optionally, constant stress) and
constant temperature (NPT or N,stress,T ensemble). It uses the
combination of Nosé-Hoover and Parrinello-Rahman dynamics proposed by
Melchionna et al. [1] and later modified by Melchionna [2]. The
differential equations are integrated using a centered difference
method [3]. Details of the implementation are available in the
document XXX NPTdynamics.tex, distributed with the module.
The dynamics object is called with the following parameters:
*atoms*:
The atoms object.
*timestep*:
The timestep in units matching eV, Å, u. Use the *units.fs* constant.
*temperature*:
The desired temperature in eV.
*externalstress*:
The external stress in eV/Å^3. Either a symmetric
3x3 tensor, a 6-vector representing the same, or a scalar
representing the pressure. Note that the stress is positive in
tension whereas the pressure is positive in compression: giving a
scalar p is equivalent to giving the tensor (-p. -p, -p, 0, 0, 0).
*ttime*:
Characteristic timescale of the thermostat. Set to None to
disable the thermostat.
*pfactor*:
A constant in the barostat differential equation. If a
characteristic barostat timescale of ptime is desired, set pfactor
to ptime^2 * B (where B is the Bulk Modulus). Set to None to
disable the barostat. Typical metallic bulk moduli are of the order
of 100 GPa or 0.6 eV/Å^3.
*mask=None*:
Optional argument. A tuple of three integers (0 or 1),
indicating if the system can change size along the three Cartesian
axes. Set to (1,1,1) or None to allow a fully flexible
computational box. Set to (1,1,0) to disallow elongations along the
z-axis etc.
Useful parameter values:
* The same *timestep* can be used as in Verlet dynamics, i.e. 5 fs is fine
for bulk copper.
* The *ttime* and *pfactor* are quite critical[4], too small values may
cause instabilites and/or wrong fluctuations in T / p. Too
large values cause an oscillation which is slow to die. Good
values for the characteristic times seem to be 25 fs for *ttime*,
and 75 fs for *ptime* (used to calculate pfactor), at least for
bulk copper with 15000-200000 atoms. But this is not well
tested, it is IMPORTANT to monitor the temperature and
stress/pressure fluctuations.
It has the following methods:
.. method:: NPT.run(n):
Perform n timesteps.
.. method:: NPT.initialize():
Estimates the dynamic variables for time=-1 to start the
algorithm. This is automatically called before the first timestep.
.. method:: NPT.set_stress():
Set the external stress. Use with care. It is
preferable to set the right value when creating the object.
.. method:: NPT.set_mask():
Change the mask. Use with care, as you may "freeze" a
fluctuation in the strain rate.
.. method:: NPT.set_strain_rate(eps):
Set the strain rate. ``eps`` must be an upper-triangular matrix.
If you set a strain rate along a direction that is "masked out"
(see ``set_mask``), the strain rate along that direction will be
maintained constantly.
.. method:: NPT.get_strain_rate():
Set the instantaneous strain rate (due to the fluctuations in the
shape of the computational box).
.. method:: NPT.get_gibbs_free_energy():
Gibbs free energy is supposed to be
preserved by this dynamics. This is mainly intended as a diagnostic
tool.
References:
[1] S. Melchionna, G. Ciccotti and B. L. Holian, Molecular Physics
78, p. 533 (1993).
[2] S. Melchionna, Physical Review E 61, p. 6165 (2000).
[3] B. L. Holian, A. J. De Groot, W. G. Hoover, and C. G. Hoover,
Physical Review A 41, p. 4552 (1990).
[4] F. D. Di Tolla and M. Ronchetti, Physical Review E 48, p. 1726 (1993).
Berendsen NPT dynamics
-----------------------
.. module:: ase.md.nptberendsen
.. class:: NPTBerendsen(atoms, timestep, temperature, taut, pressure, taup, compressibility, fixcm)
In Berendsen NPT simulations the velocities are scaled to achieve the desired
temperature. The speed of the scaling is determined by the parameter taut.
The atom positions and the simulation cell are scaled in order to achieve
the desired pressure.
This method does not result proper NPT sampling but it usually is
sufficiently good in practise (with large taut and taup). For discussion see
the gromacs manual at www.gromacs.org. or amber at ambermd.org
*atoms*:
The list of atoms.
*timestep*:
The time step.
*temperature*:
The desired temperature, in Kelvin.
*taut*:
Time constant for Berendsen temperature coupling.
*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
*fixcm*:
If True, the position and momentum of the center of mass is
kept unperturbed. Default: True.
::
# Room temperature simulation (300K, 0.1 fs time step, atmospheric pressure)
dyn = NPTBerendsen(atoms, timestep=0.1*units.fs, temperature=300,
taut=0.1*1000*units.fs, pressure = 1.01325,
taup=1.0*1000*units.fs, compressibility=4.57e-5)
|