/doc/ase/md.rst

https://gitlab.com/vote539/ase · ReStructuredText · 409 lines · 289 code · 120 blank · 0 comment · 0 complexity · 96eb9a013b0f4e5d71a240e22e754af9 MD5 · raw file

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