Skip to main content
`Ctrl`+`K`

# 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 is discussed in the Structure optimization section.

A molecular dynamics object will operate on the atoms by moving them according to their forces - it integrates Newton’s second law numerically. In the most basic form, this preserves the particle number, volume (actually, the shape of the unit cell), and the energy, generating an NVE or microcanonical ensemble. Other ensembles (NVT and NpT) can be generated by coupling the motion of the atoms to a simple model of the surroundings.

Simple example: A molecular dynamics simulation in the NVE ensemble will use the Velocity Verlet dynamics. You create the `ase.md.verlet.VelocityVerlet` object, giving it the atoms and a time step, and then you perform dynamics by calling its `run()` method:

    dyn = VelocityVerlet(atoms, dt=5.0 * units.fs,
                         trajectory='md.traj', logfile='md.log')
    dyn.run(1000)  # take 1000 steps

Note

Prior to ASE version 3.21.0, inconsistent units were used to specify temperature. Some modules expected kT (in eV), others T (in Kelvin). From ASE 3.21.0, all molecular dynamics modules expecting a temperature take a parameter `temperature_K` which is the temperature in Kelvin. For compatibility, they still accept the `temperature` parameter in the same unit as previous versions (eV or K), but using the old parameter will issue a warning.

## 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, maybe as small as 1 or 2 fs.

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 `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 `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 to not terminate correctly.)

class ase.md.MDLogger(dyn: ~typing.Any, atoms: ~ase.atoms.Atoms, logfile: ~typing.IO | str, header: bool = True, stress: bool = False, peratom: bool = False, mode: str = 'a', comm=<ase.parallel.MPI object>)

  Class for logging molecular dynamics simulations.

  Parameters: dyn: The dynamics. Only a weak reference is kept.

  atoms: The atoms.

  logfile: File name or open file, “-” meaning standard output.

  stress=False: Include stress in log.

  peratom=False: Write energies per atom.

  mode=”a”: How the file is opened if logfile is a filename.

## 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

class ase.md.verlet.VelocityVerlet(atoms: Atoms, timestep: float, trajectory: str | None = None, logfile: IO | str | None = None, loginterval: int = 1, **kwargs)

  MD with NVE ensemble and velocity Verlet time integration.

  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, default: 1) – Only write a log line for every loginterval time steps.

    - kwargs (dict, optional) – Extra arguments passed to `Dynamics`.

`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 Molecular dynamics.

## 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. Such algorithms can be stochastic or deterministic, and the coupling to the atoms typically occurs through a rescaling of the velocities. In the Langevin algorithm, a friction term and a fluctuating force is used instead.

In the NVT ensemble both the kinetic energy and the total energy will fluctuate. Some algorithms do not correctly reproduce these fluctuations, this may be an issue in particular for small systems, and is noted in the overview below.

Recommended algorithms:

Langevin dynamics

  A friction and fluctuating force is added to the equation of motion. Advantages: Simple. Correctly samples the ensemble. Disadvantage: Stochastic.

Nosé-Hoover chain

  In Nosé-Hoover dynamics, the velocities are rescaled at each time step. The scaling factor is itself a dynamic variable. For stable operation, a chain of thermostats is used. Advantages: Deterministic. Well-studied by the Stat. Mek. community. Correctly samples the ensemble. Disadvantages: The fluctuations tend to show a period given by the thermostat time scale. Slow to reach the correct temperature if started from a significantly wrong temperature.

Bussi dynamics

  Rescales the velocities at each time step, using a stochastic algorithm to ensure the correct fluctuations, unlike the closely related Berendsen algorithm. Advantages: Simple. Correctly samples the ensemble. Disadvantage: Stochastic.

Not recommended algorithms:

Andersen dynamics

  At each time step a fraction of the atoms have their velocity replaced with one drawn from the Maxwell-Boltzmann distribution. Disadvantage: Dramatically alters a few atoms instead of gently perturbing them all. Velocities are artificially decorrelated over a time corresponding to the thermalization time.

Berendsen NVT dynamics

  Rescale the velocities at each time step, so the kinetic energy exponentially approaches the correct one. Advantage: Very efficient for initializing a system to the correct initial temperature. Disadvantage: Almost completely supresses fluctuations in the total energy.

### Langevin dynamics

class ase.md.langevin.Langevin(atoms: Atoms, timestep: float, temperature: float | None = None, friction: float | None = None, fixcm: bool = True, *, temperature_K: float | None = None, rng=None, **kwargs)

  Langevin (constant N, V, T) molecular dynamics.

  Parameters:

    - atoms (Atoms object) – The list of atoms.

    - timestep (float) – The time step in ASE time units.

    - temperature (float (deprecated)) – The desired temperature, in electron volt.

    - temperature_K (float) – The desired temperature, in Kelvin.

    - friction (float) – A friction coefficient in inverse ASE time units. For example, set `0.01 / ase.units.fs` to provide 0.01 fs−1 (10 ps−1).

    - fixcm (bool (optional)) – If True, the position and momentum of the center of mass is kept unperturbed. Default: True.

    - rng (RNG object (optional)) – Random number generator, by default numpy.random. Must have a standard_normal method matching the signature of numpy.random.standard_normal.

    - **kwargs (dict, optional) – Additional arguments passed to :class:~ase.md.md.MolecularDynamics base class.

    - scalars (The temperature and friction are normally)

    - one (but in principle)

    - array. (quantity per atom could be specified by giving an)

    - propagators (RATTLE constraints can be used with these)

    - see

    - V.-Eijnden (E.)

    - Ciccotti (and G.)

    - 429 (Chem. Phys. Lett.)

    - (2006) (310)

    - used) (The propagator is Equation 23 (Eq. 39 if RATTLE constraints are)

    - another (of the above reference. That reference also contains)

    - quasi-symplectic (propagator in Eq. 21/34; but that propagator is not)

    - steps. (and gives a systematic offset in the temperature at large time)

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 Kelvin) and a friction coefficient. A typical range for the friction coefficient may be 0.001–0.1 fs−1 (1–100 ps−1). For example, you can give a friction coefficient of 0.01 fs−1 (10 ps−1) as

    dyn = Langevin(
        atoms,
        timestep=5.0 * units.fs,
        temperature_K=300.0,  # temperature in K
        friction=0.01 / units.fs,
    )

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

class ase.md.nose_hoover_chain.NoseHooverChainNVT(atoms: Atoms, timestep: float, temperature_K: float, tdamp: float, tchain: int = 3, tloop: int = 1, **kwargs)

  Isothermal molecular dynamics with Nose-Hoover chain.

  This implementation is based on the Nose-Hoover chain equations and the Liouville-operator derived integrator for non-Hamiltonian systems [1-3].

  - [1] G. J. Martyna, M. L. Klein, and M. E. Tuckerman, J. Chem. Phys. 97,

      2635 (1992). https://doi.org/10.1063/1.463940

  - [2] M. E. Tuckerman, J. Alejandre, R. López-Rendón, A. L. Jochim,

      and G. J. Martyna, J. Phys. A: Math. Gen. 39, 5629 (2006). https://doi.org/10.1088/0305-4470/39/19/S18

  - [3] M. E. Tuckerman, Statistical Mechanics: Theory and Molecular

      Simulation, Oxford University Press (2010).

  While the algorithm and notation for the thermostat are largely adapted from Ref. [4], the core equations are detailed in the implementation note available at lan496/lan496.github.io.

  - [4] M. E. Tuckerman, Statistical Mechanics: Theory and Molecular

      Simulation, 2nd ed. (Oxford University Press, 2009).

  Parameters:

    - atoms (ase.Atoms) – The atoms object.

    - timestep (float) – The time step in ASE time units.

    - temperature_K (float) – The target temperature in K.

    - tdamp (float) – The characteristic time scale for the thermostat in ASE time units. Typically, it is set to 100 times of \(timestep\).

    - tchain (int) – The number of thermostat variables in the Nose-Hoover chain.

    - tloop (int) – The number of sub-steps in thermostat integration.

    - **kwargs (dict, optional) – Additional arguments passed to :class:~ase.md.md.MolecularDynamics base class.

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. To give a more stable dynamics, the friction coefficient is itself thermalized in the same way, leading to a chain of thermostats (a Nosé-Hoover chain, in ASE per default the chain length is 3).

### Bussi dynamics

class ase.md.bussi.Bussi(atoms, timestep, temperature_K, taut, rng=None, **kwargs)

  Bussi stochastic velocity rescaling (NVT) molecular dynamics.

  Based on the paper from Bussi et al., J. Chem. Phys. 126, 014101 (2007) (also available from https://arxiv.org/abs/0803.4060).

  Parameters:

    - atoms (Atoms) – The atoms object.

    - timestep (float) – The time step in ASE time units.

    - temperature_K (float) – The desired temperature, in Kelvin.

    - taut (float) – Time constant for Bussi temperature coupling in ASE time units.

    - rng (RNG object, optional) – Random number generator, by default numpy.random.

    - **kwargs (dict, optional) – Additional arguments are passed to :class:~ase.md.md.MolecularDynamics base class.

The Bussi class implements the Bussi dynamics, where constant temperature is imposed by a stochastic velocity rescaling algorithm. The thermostat is conceptually very close to the Berendsen thermostat, but does sample the canonical ensemble correctly. Given that the thermostat is both correct and global, it is advised to use it for production runs.

### Andersen dynamics

class ase.md.andersen.Andersen(atoms: ~ase.atoms.Atoms, timestep: float, temperature_K: float, andersen_prob: float, fixcm: bool = True, rng=<module 'numpy.random' from '/usr/lib/python3/dist-packages/numpy/random/__init__.py'>, **kwargs)

  Andersen (constant N, V, T) molecular dynamics.

  Parameters:

    - atoms (Atoms object) – The list of atoms.

    - timestep (float) – The time step in ASE time units.

    - temperature_K (float) – The desired temperature, in Kelvin.

    - andersen_prob (float) – A random collision probability, typically 1e-4 to 1e-1. With this probability atoms get assigned random velocity components.

    - fixcm (bool (optional)) – If True, the position and momentum of the center of mass is kept unperturbed. Default: True.

    - rng (RNG object, default: `numpy.random`) – Random number generator. This must have the `random` method with the same signature as `numpy.random.random`.

    - **kwargs (dict, optional) – Additional arguments passed to :class:~ase.md.md.MolecularDynamics base class.

    - bath (The temperature is imposed by stochastic collisions with a heat)

    - particles. (that acts on velocity components of randomly chosen)

    - velocities (The algorithm randomly decorrelates)

    - properties (so dynamical)

    - measured. (like diffusion or viscosity cannot be properly)

    - Andersen (H. C.)

    - (4) (J. Chem. Phys. 72)

    - (1980) (2384–2393)

The Andersen class implements Andersen dynamics, where constant temperature is imposed by stochastic collisions with a heat bath. With a (small) probability (`andersen_prob`) the collisions act occasionally on velocity components of randomly selected particles. Upon a collision the new velocity is drawn from the Maxwell-Boltzmann distribution at the corresponding temperature. The system is then integrated numerically at constant energy according to the Newtonian laws of motion. The collision probability is defined as the average number of collisions per atom and timestep. The algorithm generates a canonical distribution. [1] However, due to the random decorrelation of velocities, the dynamics are unphysical and cannot represent dynamical properties like e.g. diffusion or viscosity. Another disadvantage is that the collisions are stochastic in nature, so repeating the simulation will not give exactly the same trajectory.

When the `Andersen` object is created, you must specify a time step, a temperature (in Kelvin) and a collision probability. Typical values for this probability are in the order of 1e-4 to 1e-1.

    # Room temperature simulation (300 Kelvin, Andersen probability: 0.002)
    dyn = Andersen(atoms, 5 * units.fs, 300, 0.002)

References:

[1] D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic Press, London, 1996)

### Berendsen NVT dynamics

class ase.md.nvtberendsen.NVTBerendsen(atoms: Atoms, timestep: float, temperature: float | None = None, taut: float | None = None, fixcm: bool = True, *, temperature_K: float | None = None, **kwargs)

  Berendsen (constant N, V, T) molecular dynamics.

  Parameters:

    - atoms (Atoms object) – The list of atoms.

    - timestep (float) – The time step in ASE time units.

    - temperature (float) – The desired temperature, in Kelvin.

    - temperature_K (float) – Alias for temperature

    - taut (float) – Time constant for Berendsen temperature coupling in ASE time units.

    - fixcm (bool (optional)) – If True, the position and momentum of the center of mass is kept unperturbed. Default: True.

    - **kwargs (dict, optional) – Additional arguments passed to :class:~ase.md.md.MolecularDynamics base class.

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 practice (with large taut). For discussion see the gromacs manual at www.gromacs.org.

    # 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)

Constant pressure (or for solids, constant stress) is usually obtained by adding a barostat to one of the NVT algorithms above. ASE currently lacks a good NPT algorithm. The following two are available.

Algorithms:

Berendsen NPT dynamics

  This is a variation of Berendsen NVT dynamics with a barostat added. The size of the unit cell is rescaled after each time step, so the pressure / stress approaches the desired pressure. It exists in two variations, one where the shape of the unit cell is preserved and one where it is allowed to vary. Disadvantage: Fluctuations in both total energy and pressure are suppressed compared to the correct NPT ensemble. For large systems, this is not expected to be serious.

NPT

  An implementation of NPT dynamics combining a Nosé-Hoover thermostat with a Parinello-Rahman barostat, according to Melchionna et al., see below. Not recommended! The dynamics tend to be unstable, especially if started with a temperature or pressure that is different from the desired. The fluctuations seem to often be wrong.

### Berendsen NPT dynamics

class ase.md.nptberendsen.NPTBerendsen(atoms: Atoms, timestep: float, temperature: float | None = None, *, temperature_K: float | None = None, pressure: float | None = None, pressure_au: float | None = None, taut: float = 49.11347394232032, taup: float = 98.22694788464064, compressibility: float | None = None, compressibility_au: float | None = None, fixcm: bool = True, trajectory: str | None = None, logfile: IO | str | None = None, loginterval: int = 1, append_trajectory: bool = False)

  Berendsen (constant N, P, T) molecular dynamics.

  This dynamics scale the velocities and volumes to maintain a constant pressure and temperature. The shape of the simulation cell is not altered, if that is desired use Inhomogenous_NPTBerendsen.

  Parameters:

  atoms: Atoms object

    The list of atoms.

  timestep: float

    The time step in ASE time units.

  temperature: float

    The desired temperature, in Kelvin.

  temperature_K: float

    Alias for `temperature`.

  pressure: float (deprecated)

    The desired pressure, in bar (1 bar = 1e5 Pa). Deprecated, use `pressure_au` instead.

  pressure_au: float

    The desired pressure, in atomic units (eV/Å^3).

  taut: float

    Time constant for Berendsen temperature coupling in ASE time units. Default: 0.5 ps.

  taup: float

    Time constant for Berendsen pressure coupling. Default: 1 ps.

  compressibility: float (deprecated)

    The compressibility of the material, in bar-1. Deprecated, use `compressibility_au` instead.

  compressibility_au: float

    The compressibility of the material, in atomic units (Å^3/eV).

  fixcm: bool (optional)

    If True, the position and momentum of the center of mass is kept unperturbed. Default: True.

  trajectory: Trajectory object or str (optional)

    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.

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. The time scale of this barostat is given by the parameter taup.

This method does not result proper NPT sampling but it usually is sufficiently good in practice (with large taut and taup). For discussion see the gromacs manual at www.gromacs.org. or amber at ambermd.org

    # Room temperature simulation (300K, 0.1 fs time step, atmospheric pressure)
    dyn = NPTBerendsen(atoms, timestep=0.1 * units.fs, temperature_K=300,
                       taut=100 * units.fs, pressure_au=1.01325 * units.bar,
                       taup=1000 * units.fs, compressibility_au=4.57e-5 / units.bar)

### Nosé-Hoover-Parinello-Rahman NPT dynamics

class ase.md.npt.NPT(atoms: Atoms, timestep: float, temperature: float | None = None, externalstress: float | None = None, ttime: float | None = None, pfactor: float | None = None, *, temperature_K: float | None = None, mask: Tuple[int] | ndarray | None = None, trajectory: str | None = None, logfile: IO | str | None = None, loginterval: int = 1, append_trajectory: bool = False)

  Constant pressure/stress and temperature dynamics.

  Combined Nose-Hoover and Parrinello-Rahman dynamics, creating an NPT (or N,stress,T) ensemble.

  The method is the one proposed by Melchionna et al. [1] and later modified by Melchionna [2]. The differential equations are integrated using a centered difference method [3]. See also NPTdynamics.tex

  The dynamics object is called with the following parameters:

  atoms: Atoms object

    The list of atoms.

  timestep: float

    The timestep in units matching eV, Å, u.

  temperature: float (deprecated)

    The desired temperature in eV.

  temperature_K: float

    The desired temperature in K.

  externalstress: float or nparray

    The external stress in eV/A^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: float

    Characteristic timescale of the thermostat, in ASE internal units Set to None to disable the thermostat.

    WARNING: Not specifying ttime sets it to None, disabling the thermostat.

  pfactor: float

    A constant in the barostat differential equation. If a characteristic barostat timescale of ptime is desired, set pfactor to ptime^2 * B (where ptime is in units matching eV, Å, u; and B is the Bulk Modulus, given in eV/Å^3). Set to None to disable the barostat. Typical metallic bulk moduli are of the order of 100 GPa or 0.6 eV/A^3.

    WARNING: Not specifying pfactor sets it to None, disabling the barostat.

  mask: None or 3-tuple or 3x3 nparray (optional)

    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. mask may also be specified as a symmetric 3x3 array indicating which strain values may change.

  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.

  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).

  run(steps)

    Perform a number of time steps.

    Parameters:

      steps (int) – Number of dynamics steps to be run.

    Yields:

      complete (bool) – True if the maximum number of steps are reached.

  set_stress(stress)

    Set the applied stress.

    Must be a symmetric 3x3 tensor, a 6-vector representing a symmetric 3x3 tensor, or a number representing the pressure.

    Use with care, it is better to set the correct stress when creating the object.

  set_temperature(temperature=None, *, temperature_K=None)

    Set the temperature.

    Parameters:

    temperature: float (deprecated)

      The new temperature in eV. Deprecated, use `temperature_K`.

    temperature_K: float (keyword-only argument)

      The new temperature, in K.

  set_mask(mask)

    Set the mask indicating dynamic elements of the computational box.

    If set to None, all elements may change. If set to a 3-vector of ones and zeros, elements which are zero specify directions along which the size of the computational box cannot change. For example, if mask = (1,1,0) the length of the system along the z-axis cannot change, although xz and yz shear is still possible. May also be specified as a symmetric 3x3 array indicating which strain values may change.

    Use with care, as you may “freeze in” a fluctuation in the strain rate.

  set_fraction_traceless(fracTraceless)

    set what fraction of the traceless part of the force on eta is kept.

    By setting this to zero, the volume may change but the shape may not.

  get_strain_rate()

    Get the strain rate as a triangular 3x3 matrix.

    This includes the fluctuations in the shape of the computational box.

  set_strain_rate(rate)

    Set the strain rate. Must be a triangular 3x3 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.

  get_time()

    Get the elapsed time.

  initialize()

    Initialize the dynamics.

    The dynamics requires positions etc for the two last times to do a timestep, so the algorithm is not self-starting. This method performs a ‘backwards’ timestep to generate a configuration before the current.

    This is called automatically the first time `run()` is called.

  get_gibbs_free_energy()

    Return the Gibb’s free energy, which is supposed to be conserved.

    Requires that the energies of the atoms are up to date.

    This is mainly intended as a diagnostic tool. If called before the first timestep, Initialize will be called.

  zero_center_of_mass_momentum(verbose=0)

    Set the center of mass momentum to zero.

  attach(function, interval=1, *args, **kwargs)

    Attach callback function or trajectory.

    At every interval steps, call function with arguments args and keyword arguments kwargs.

    If function is a trajectory object, its write() method is attached, but if function is a BundleTrajectory (or another trajectory supporting set_extra_data(), said method is first used to instruct the trajectory to also save internal data from the NPT dynamics object.

This module is not recommended! There is a strong tendency for oscillations in the temperature and/or pressure, unless the starting configuration is chosen with great care.

## Contour Exploration

class ase.md.contour_exploration.ContourExploration(atoms: ~ase.atoms.Atoms, maxstep: float = 0.5, parallel_drift: float = 0.1, energy_target: float | None = None, angle_limit: float | None = 20.0, potentiostat_step_scale: float | None = None, remove_translation: bool = False, use_frenet_serret: bool = True, initialization_step_scale: float = 0.01, use_target_shift: bool = True, target_shift_previous_steps: int = 10, use_tangent_curvature: bool = False, rng=<module 'numpy.random' from '/usr/lib/python3/dist-packages/numpy/random/__init__.py'>, force_consistent: bool | None = None, trajectory: str | None = None, logfile: ~typing.IO | str | None = None, append_trajectory: bool = False, loginterval: int = 1)

  Contour Exploration object.

  Parameters:

  atoms: Atoms object

    The Atoms object to operate on. Atomic velocities are required for the method. If the atoms object does not contain velocities, random ones will be applied.

  maxstep: float

    Used to set the maximum distance an atom can move per iteration (default value is 0.5 Å).

  parallel_drift: float

    The fraction of the update step that is parallel to the contour but in a random direction. Used to break symmetries.

  energy_target: float

    The total system potential energy for that the potentiostat attepts to maintain. (defaults the initial potential energy)

  angle_limit: float or None

    Limits the stepsize to a maximum change of direction angle using the curvature. Gives a scale-free means of tuning the stepsize on the fly. Typically less than 30 degrees gives reasonable results but lower angle limits result in higher potentiostatic accuracy. Units of degrees. (default 20°)

  potentiostat_step_scale: float or None

    Scales the size of the potentiostat step. The potentiostat step is determined by linear extrapolation from the current potential energy to the target_energy with the current forces. A potentiostat_step_scale > 1.0 overcorrects and < 1.0 undercorrects. By default, a simple heuristic is used to selected the valued based on the parallel_drift. (default None)

  remove_translation: boolean

    When True, the net momentum is removed at each step. Improves potentiostatic accuracy slightly for bulk systems but should not be used with constraints. (default False)

  use_frenet_serret: Bool

    Controls whether or not the Taylor expansion of the Frenet-Serret formulas for curved path extrapolation are used. Required for using angle_limit based step scalling. (default True)

  initialization_step_scale: float

    Controls the scale of the initial step as a multiple of maxstep. (default 1e-2)

  use_target_shift: boolean

    Enables shifting of the potentiostat target to compensate for systematic undercorrection or overcorrection by the potentiostat. Uses the average of the target_shift_previous_steps to prevent coupled occilations. (default True)

  target_shift_previous_steps: int

    The number of pevious steps to average when using use_target_shift. (default 10)

  use_tangent_curvature: boolean

    Use the velocity unit tangent rather than the contour normals from forces to compute the curvature. Usually not as accurate. (default False)

  rng: a random number generator

    > Lets users control the random number generator for the parallel_drift vector. (default numpy.random)

    force_consistent: boolean

      (default None)

  trajectory: Trajectory object or str (optional)

    Attach trajectory object. If trajectory is a string a Trajectory will be constructed. Default: None.

  logfile: file object or str (optional)

    If logfile is a string, a file with that name will be opened. Use ‘-’ for stdout. Default: None.

  loginterval: int (optional)

    Only write a log line for every loginterval time steps. Default: 1

  append_trajectory: boolean

    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.

Contour Exploration evolves the system along constant potentials energy contours on the potential energy surface. The method uses curvature based extrapolation and a potentiostat to correct for potential energy errors. It is similar to molecular dynamics but with a potentiostat rather than a thermostat. Without changes in kinetic energy, it is more useful to automatically scale step sizes to the curvature of the potential energy contour via an `angle_limit` while enforcing a `maxstep` to ensure potentiostatic accuracy. [1] To escape loops on the pontential energy surface or to break symmetries, a random drift vector parallel to the contour can be applied as a fraction of the step size via `parallel_drift`. Contour exploration cannot be used at minima since the contour is a single point.

    # Contour exploration at the current potential energy
    dyn = ContourExploration(atoms)

References:

[1] M. J. Waters and J. M. Rondinelli, Contour Exploration with Potentiostatic Kinematics arXiv:2103.08054.

## Velocity distributions

A selection of functions are provided to initialize atomic velocities to the correct temperature.

ase.md.velocitydistribution.MaxwellBoltzmannDistribution(atoms: ~ase.atoms.Atoms, temp: float | None = None, *, temperature_K: float | None = None, comm=<ase.parallel.MPI object>, communicator=None, force_temp: bool = False, rng=None)

  Set the atomic momenta to a Maxwell-Boltzmann distribution.

  Parameters:

    - atoms (Atoms object) – The atoms. Their momenta will be modified.

    - temp (float (deprecated)) – The temperature in eV. Deprecated, use `temperature_K` instead.

    - temperature_K (float) – The temperature in Kelvin.

    - comm (MPI communicator, default: `ase.parallel.world`) –

      Communicator used to distribute an identical distribution to all tasks.

      Added in version 3.26.0.

    - communicator –

      Deprecated since version 3.26.0: To be removed in ASE 3.27.0 in favor of `comm`.

    - force_temp (bool (optional, default: False)) – If True, the random momenta are rescaled so the kinetic energy is exactly 3/2 N k T. This is a slight deviation from the correct Maxwell-Boltzmann distribution.

    - rng (Numpy RNG (optional)) – Random number generator. Default: numpy.random If you would like to always get the identical distribution, you can supply a random seed like `rng=numpy.random.RandomState(seed)`, where seed is an integer.

ase.md.velocitydistribution.Stationary(atoms: Atoms, preserve_temperature: bool = True)

  Sets the center-of-mass momentum to zero.

ase.md.velocitydistribution.ZeroRotation(atoms: Atoms, preserve_temperature: float = True)

  Sets the total angular momentum to zero by counteracting rigid rotations.

ase.md.velocitydistribution.PhononHarmonics(atoms, force_constants, temp=None, *, temperature_K=None, rng=<module 'numpy.random' from '/usr/lib/python3/dist-packages/numpy/random/__init__.py'>, quantum=False, plus_minus=False, return_eigensolution=False, failfast=True)

  Excite phonon modes to specified temperature.

  This will displace atomic positions and set the velocities so as to produce a random, phononically correct state with the requested temperature.

  Parameters:

  atoms: ase.atoms.Atoms() object

    Positions and momenta of this object are perturbed.

  force_constants: ndarray of size 3N x 3N

    Force constants for the the structure represented by atoms in eV/Å²

  temp: float (deprecated).

    Temperature in eV. Deprecated, use `temperature_K` instead.

  temperature_K: float

    Temperature in Kelvin.

  rng: Random number generator

    RandomState or other random number generator, e.g., np.random.rand

  quantum: bool

    True for Bose-Einstein distribution, False for Maxwell-Boltzmann (classical limit)

  failfast: bool

    True for sanity checking the phonon spectrum for negative frequencies at Gamma.

ase.md.velocitydistribution.phonon_harmonics(force_constants, masses, temp=None, *, temperature_K=None, rng=<bound method RandomState.rand of RandomState(MT19937)>, quantum=False, plus_minus=False, return_eigensolution=False, failfast=True)

  Return displacements and velocities that produce a given temperature.

  Parameters:

  force_constants: array of size 3N x 3N

    force constants (Hessian) of the system in eV/Å²

  masses: array of length N

    masses of the structure in amu

  temp: float (deprecated)

    Temperature converted to eV (T * units.kB). Deprecated, use `temperature_K`.

  temperature_K: float

    Temperature in Kelvin.

  rng: function

    Random number generator function, e.g., np.random.rand

  quantum: bool

    True for Bose-Einstein distribution, False for Maxwell-Boltzmann (classical limit)

  plus_minus: bool

    Displace atoms with +/- the amplitude accoding to PRB 94, 075125

  return_eigensolution: bool

    return eigenvalues and eigenvectors of the dynamical matrix

  failfast: bool

    True for sanity checking the phonon spectrum for negative frequencies at Gamma

  Returns:

  Displacements, velocities generated from the eigenmodes, (optional: eigenvalues, eigenvectors of dynamical matrix)

  Purpose:

  Excite phonon modes to specified temperature.

  This excites all phonon modes randomly so that each contributes, on average, equally to the given temperature. Both potential energy and kinetic energy will be consistent with the phononic vibrations characteristic of the specified temperature.

  In other words the system will be equilibrated for an MD run at that temperature.

  force_constants should be the matrix as force constants, e.g., as computed by the ase.phonons module.

  Let X_ai be the phonon modes indexed by atom and mode, w_i the phonon frequencies, and let 0 < Q_i <= 1 and 0 <= R_i < 1 be uniformly random numbers. Then

                   1/2
      _     / k T \     ---  1  _             1/2
      R  += | --- |      >  --- X   (-2 ln Q )    cos (2 pi R )
       a    \  m  /     ---  w   ai         i                i
                a        i    i

                   1/2
      _     / k T \     --- _            1/2
      v   = | --- |      >  X  (-2 ln Q )    sin (2 pi R )
       a    \  m  /     ---  ai        i                i
                a        i

  Reference: [West, Estreicher; PRL 96, 22 (2006)]

## Post-simulation Analysis

Functionality is provided to perform analysis of atomic/molecular behaviour as calculation in a molecular dynamics simulation. Currently, this is presented as a class to address the Einstein equation for diffusivity.

class ase.md.analysis.DiffusionCoefficient(traj, timestep, atom_indices=None, molecule=False)

  This class calculates the Diffusion Coefficient for the given Trajectory using the Einstein Equation:

  ..math:: left langle left | r(t) - r(0) right | ^{2} right rangle = 2nDt

  where r(t) is the position of atom at time t, n is the degrees of freedom and D is the Diffusion Coefficient

  Solved herein by fitting with \(y = mx + c\), i.e. \(\frac{1}{2n} \left \langle \left | r(t) - r(0) \right | ^{2} \right \rangle = Dt\), with m = D and c = 0

  wiki : https://en.wikibooks.org/wiki/Molecular_Simulation/Diffusion_Coefficients

  Parameters:

    - traj (Trajectory) – Trajectory of atoms objects (images)

    - timestep (Float) – Timestep between each image in the trajectory, in ASE timestep units (For an MD simulation with timestep of N, and images written every M iterations, our timestep here is N * M)

    - atom_indices (List of Int) – The indices of atoms whose Diffusion Coefficient is to be calculated explicitly

    - molecule (Boolean) – Indicate if we are studying a molecule instead of atoms, therefore use centre of mass in calculations