/doc/ase/md.rst

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