cp2k/docs/methods/sampling/molecular_dynamics.md

# Molecular Dynamics

Molecular dynamics (MD) simulates how atoms and molecules move in time. It is widely used to sample
finite-temperature behaviour, compute thermodynamic averages, and analyse dynamical properties such
as diffusion, vibrations, structural fluctuations, and reaction events. MD trajectories are also a
common starting point for more advanced simulation and analysis techniques.

In CP2K, the forces used for MD are provided by a selected `FORCE_EVAL` method. The same MD driver
can therefore be used with classical force fields, machine-learning potentials, QM/MM,
semi-empirical methods, and Born-Oppenheimer ab-initio molecular dynamics (AIMD).

```{youtube} cGDiVSWpXLc
---
url_parameters: ?start=2
align: center
privacy_mode:
---
```

## Basic setup

Select an MD calculation with [`RUN_TYPE MD`](#CP2K_INPUT.GLOBAL.RUN_TYPE). The main MD controls are
in [MOTION/MD](#CP2K_INPUT.MOTION.MD): [ENSEMBLE](#CP2K_INPUT.MOTION.MD.ENSEMBLE) selects the
propagator, [STEPS](#CP2K_INPUT.MOTION.MD.STEPS) or [MAX_STEPS](#CP2K_INPUT.MOTION.MD.MAX_STEPS)
sets the length of the run, [TIMESTEP](#CP2K_INPUT.MOTION.MD.TIMESTEP) sets the integration step,
and [TEMPERATURE](#CP2K_INPUT.MOTION.MD.TEMPERATURE) is used for velocity initialisation and as the
constant-temperature target.

A minimal fixed-energy MD block is:

```none
&GLOBAL
  RUN_TYPE MD
  PROJECT_NAME my_md
&END GLOBAL

&MOTION
  &MD
    ENSEMBLE NVE
    STEPS 10000
    TIMESTEP 0.5
    TEMPERATURE 300
  &END MD
&END MOTION
```

The force method is defined independently in [FORCE_EVAL](#CP2K_INPUT.FORCE_EVAL).

## Preliminary considerations

Universal to every atomistic simulation and every program, a number of crucial points must be
carefully considered before setting up a molecular dynamics task and investing copious amounts of
computational resources. As a starter, this section discusses preliminary factors in the choice of
the spatiotemporal scale of the trajectory, periodic boundary condition. force method, and the
initial structure.

### Space and time scale

Even when taking the ergodic hypothesis as granted, statistical analysis on target properties at
thermodynamic equilibrium has to be performed on top of sufficient sampling. The (auto)correlation
length and (auto)correlation time provide estimates for the space and time scale where measurements
of some property retain a memory of past states and are not fully statistically independent. It is
important to learn about their typical values in the same or similar systems, and take adequate
multiples as the size of the system and the length of the trajectory.

On the other hand, non-equilibrium irreversible processes such as phase transitions and chemical
reactions are also observable at some time and space resolution. Experimentally it is determined by
the sensitivity of instrumental characterization and analysis methods, with certain precise
definitions of the states before and after the process. Unfortunately, many experiments take place
in real life in a macroscopic time interval that exceeds the capability of simulation; for instance,
a chemical reaction taking days, hours, or even just seconds to happen may not be replicatable in
AIMD simulations with trajectory on the picosecond or nanosecond scale starting with only a bunch of
reactants randomly packed together, however thermodynamically favorable it is. Again, it is better
to survey literatures with relevant kinetic data and have a clear understanding of the typical free
energy barriers and reaction rates before proceeding to simulations. Even static thermochemistry
calculations for the profile of energy landscape via optimization of minima and transition states,
vibrational analysis and evaluation of single-point energy would be helpful. To simulate difficult
reactions by molecular dynamics in a short trajectory, advanced _enhanced sampling_ methods such as
[](./metadynamics) may be necessary.

### Periodic boundary condition

**Periodic boundary condition** ({term}`PBC`) may or may not be used in the molecular dynamics
simulation, and there are caveats for both ends.

On one hand, very small periodic cells for production MD should be avoided unless the finite-size
effects are known to be acceptable. Small cells produce artificial correlations through periodic
images, exaggerate statistical temperature and pressure fluctuations, and can make NPT cell dynamics
unstable or unphysical. Sometimes using a supercell of a crystalline structure is preferred if its
primitive cell is too small.

On the other hand, without PBC in MD, a common risk is that strong repulsion or very weak attraction
causes the system to expand and effectively diminishes desired interactions between molecules. A
possible remedy is to define an external potential or a constraint that act as a soft wall to
reflect or push stray molecules back to the center, but this is artificial and arbitrary.

### Force method

To propagate the trajectory in a molecular dynamics simulation, the energy and force is evaluated on
each timestep to provide acceleration via Newton's second law. The choice of the force method
determines the description of chemical bonding and/or weak interactions, the configuration space
that can be sampled, and the eventual computational cost. Careful selection of function forms and
parameters as well as thorough validation against the target property is vital to any application;
in particular, the part of "cross" interactions between two or more components described with more
than one unified force method requires very well-founded justification.

The strength and limitation of each force method is decided by design and parametrization. For
instance, classical force fields in molecular mechanics usually model chemical bonds by harmonic
potentials on the stretching, bending and rigid dihedral terms (and periodic torsion potentials on
the flexible dihedral terms). The simplicity of function forms allows for speedy evaluation and the
conformational fluctuation or conversion can be captured well, but the bonding pattern is fixed to
near equilibrium and bond breaking and forming in chemical reactions cannot be captured. To describe
the bond rearrangement and the underlying electronic structure, AIMD based on quantum chemical
methods becomes requisite even though it is slower and more costly. The middle grounds have recently
becoming filled with emerging methods like the machine-learning potentials.

### Initial structure

Because most of the points apply here as well, please first refer to the documentation about the
[starting structure and cell](../optimization/geometry_and_cell_opt.md#starting-structure-and-cell)
of a geometry optimization.

In general, start from a physically reasonable structure with sensible bond lengths, intermolecular
distances, density, and cell parameters. A majority of the force methods demonstrate extremely
strong repulsion at very small nuclei distances, and thus molecular dynamics simulation employing
these methods is not a reliable way to fix severe close contacts, unrealistic densities, or badly
prepared cells, because such problems can cause unstable forces, SCF failures, or immediate heating
and atom ejection.

In fact, when starting from scratch, a geometry optimization preceding molecular dynamics simulation
is highly recommended. During geometry optimization, the structure is relaxed and the maximum force
on atoms is significantly lowered, which puts an upper limit to the initial acceleration in the
subsequent molecular dynamics simulation and prevents rapid atomic motion and extreme temperature
rise. With electronic-structure methods, the geometry optimization also offers a chance to confirm
SCF convergence and estimate the time consumption on each step. The only "downside" is that the
temperature may start from or drop to a relatively low value at early simulation due to the small
forces and slow motion, and may take a little longer to heat up to the desired temperature; see the
section on equilibration below for how to address this.

## Initial velocities

If velocities are not provided by a restart file or an external trajectory, CP2K can initialise
atomic velocities from [TEMPERATURE](#CP2K_INPUT.MOTION.MD.TEMPERATURE). The
[INITIALIZATION_METHOD](#CP2K_INPUT.MOTION.MD.INITIALIZATION_METHOD) keyword controls how the
initial positions and velocities are generated.

For isolated molecules, clusters, slabs, or other systems where global translation or rotation is
not part of the intended dynamics, it can be useful to remove centre-of-mass translation and, when
meaningful, overall angular motion; see [COMVEL_TOL](#CP2K_INPUT.MOTION.MD.COMVEL_TOL),
[ANGVEL_TOL](#CP2K_INPUT.MOTION.MD.ANGVEL_TOL), [ANGVEL_ZERO](#CP2K_INPUT.MOTION.MD.ANGVEL_ZERO),
and related keywords. For periodic bulk systems, small total-momentum drift is usually less
important.

## Choosing an ensemble

The most common choices are `NVE`, `NVT`, `NPT_I`, `NPT_F`, and `LANGEVIN`.

- `NVE` is microcanonical dynamics. Use it to test the time step, force convergence, grid settings,
  and extrapolation through energy conservation. It is also suitable for production when isolated
  Hamiltonian dynamics is desired.
- `NVT` is canonical sampling at fixed cell. It is usually the most convenient ensemble for
  equilibration and for production when the volume is known or fixed.
- `NPT_I` uses isotropic cell fluctuations. It is appropriate when only the cell volume should
  fluctuate, for example for liquids, gases, or approximately isotropic solids.
- `NPT_F` uses a flexible cell. Use it when the cell shape and anisotropic stress relaxation are
  part of the physical question, for example for solids or interfaces.
- `LANGEVIN` is stochastic constant-temperature dynamics. It is useful for robust equilibration,
  dissipative dynamics, or models where a stochastic heat bath is intentional.

Other supported values of [ENSEMBLE](#CP2K_INPUT.MOTION.MD.ENSEMBLE) are more specialised: `NPE_I`,
`NPE_F`, `ISOKIN`, `REFTRAJ`, `MSST`, `MSST_DAMPED`, `HYDROSTATICSHOCK`, `NVT_ADIABATIC`, and
`NPT_IA`. See the input reference for details.

Note that NPT simulations are sensitive to system size and barostat settings; small cells can show
large pressure fluctuations and artificial cell oscillations. For small cells, consider fixed-cell
`NVT` after a separate cell optimisation, or use a larger supercell before production `NPT`.

## Thermostats

Constant-temperature ensembles use [MOTION/MD/THERMOSTAT](#CP2K_INPUT.MOTION.MD.THERMOSTAT). The
input keyword default for [TYPE](#CP2K_INPUT.MOTION.MD.THERMOSTAT.TYPE) is currently `NOSE`, but for
routine `NVT` equilibration and production `CSVR` is often the safer practical default.

- `CSVR`: recommended general-purpose choice. It gives canonical sampling through velocity rescaling
  and is robust for both equilibration and production.
- `NOSE`: Nose-Hoover chain thermostat. It can be appropriate for well-equilibrated systems, but it
  is less robust when the system is far from equilibrium or very small, and should not be the
  default recommendation for new users.
- `GLE`: generalised Langevin thermostat. Use it for specialised coloured-noise applications or
  workflows that need a fitted GLE kernel.
- `AD_LANGEVIN`: adaptive Langevin thermostat. Use it when adaptive stochastic thermostatting is
  specifically desired.

A practical `NVT` setup is:

```none
&MOTION
  &MD
    ENSEMBLE NVT
    TEMPERATURE 300
    TIMESTEP 0.5
    &THERMOSTAT
      TYPE CSVR
      &CSVR
        TIMECON 100.0
      &END CSVR
    &END THERMOSTAT
  &END MD
&END MOTION
```

The [CSVR/TIMECON](#CP2K_INPUT.MOTION.MD.THERMOSTAT.CSVR.TIMECON) (unit: fs) value controls the
coupling strength. Useful starting values are `50` for rapid heating or cooling, `100` as a simple
robust default, and `200` or larger for gentler production sampling after equilibration.

Do not use crude velocity rescaling as a production thermostat.
[TEMP_TOL](#CP2K_INPUT.MOTION.MD.TEMP_TOL) can trigger velocity rescaling when the instantaneous
temperature deviates from the target, but this keyword is obsolescent; a CSVR thermostat with a
short time constant is usually the better early- equilibration replacement.

## Barostats and NPT simulations

NPT simulations additionally use [MOTION/MD/BAROSTAT](#CP2K_INPUT.MOTION.MD.BAROSTAT). The pressure
target is set by [PRESSURE](#CP2K_INPUT.MOTION.MD.BAROSTAT.PRESSURE), and the barostat time scale by
[TIMECON](#CP2K_INPUT.MOTION.MD.BAROSTAT.TIMECON). For `NPT_F`,
[VIRIAL](#CP2K_INPUT.MOTION.MD.BAROSTAT.VIRIAL) can restrict which Cartesian components of the
virial are used for cell relaxation.

Use `NPT_I` when isotropic volume fluctuations are sufficient. Use `NPT_F` only when the cell shape
should respond to anisotropic stress. Full-cell NPT is powerful but can amplify problems from small
cells, poor initial structures, noisy stresses, or too aggressive barostat coupling. A robust
workflow is to equilibrate first in `NVT`, then run a short NPT test, inspect the cell trajectory
and pressure oscillations, and only then start production sampling. If the cell volume or shape
oscillates strongly, use a larger barostat time constant, a larger supercell, or fixed-cell `NVT` if
pressure sampling is not essential.

## Time step and validation

The time step must resolve the fastest relevant nuclear motion. For AIMD with light atoms, 0.5 ~ 1.0
fs is a common starting range. Use a smaller value for high temperature, weak bonds, reactive
events, poor starting structures, or hydrogen-rich systems. Constraints, deuteration, classical
force fields, or multiple-time-step schemes may allow larger time steps, but always validate the
choice.

A sensitive validation test is a short `NVE` trajectory. The total-energy oscillation and long-time
drift should be small relative to the kinetic energy or target temperature. Energy conservation
alone is not always sufficient for heterogeneous systems; also inspect the region or property of
interest when only part of the system matters.

## Electronic-structure settings for AIMD

Each AIMD step requires an electronic-structure calculation. Important Quickstep settings include
[EPS_SCF](#CP2K_INPUT.FORCE_EVAL.DFT.SCF.EPS_SCF),
[EPS_DEFAULT](#CP2K_INPUT.FORCE_EVAL.DFT.QS.EPS_DEFAULT), and the plane-wave grid parameters in
[MGRID](#CP2K_INPUT.FORCE_EVAL.DFT.MGRID), especially
[CUTOFF](#CP2K_INPUT.FORCE_EVAL.DFT.MGRID.CUTOFF) and
[REL_CUTOFF](#CP2K_INPUT.FORCE_EVAL.DFT.MGRID.REL_CUTOFF).

AIMD settings do not always need to match high-precision static energies or final property
calculations. For long trajectories, stable forces and negligible drift are often more important
than over-converging every SCF cycle. For example, if one element in your system use `CUTOFF 400` in
optimization and/or energy tasks, then in AIMD you can set to `CUTOFF 300`. However, too loose
settings can create noisy forces, poor energy conservation, artificial heating, or wrong dynamics.

## Extrapolation of the electronic initial guess

In a normal MD trajectory, consecutive geometries are close to each other. CP2K can use previous
electronic solutions to build a better initial guess for the next step through
[EXTRAPOLATION](#CP2K_INPUT.FORCE_EVAL.DFT.QS.EXTRAPOLATION) and
[EXTRAPOLATION_ORDER](#CP2K_INPUT.FORCE_EVAL.DFT.QS.EXTRAPOLATION_ORDER). A good extrapolation keeps
the electronic state continuous and usually reduces the number of SCF iterations at each MD step,
which is often crucial for efficiency.

For standard Born-Oppenheimer MD, `ASPC` is the default and generally the best first choice. It is
similar in spirit to `PS`, but is designed for MD stability rather than only initial-guess accuracy.
`EXTRAPOLATION_ORDER 3` is a good starting point for `ASPC` and `PS`.

`GEXT_PROJ` and `GEXT_PROJ_QTR` are also worth considering. They can use higher extrapolation
orders, typically 4 ~ 10. `GEXT_PROJ_QTR` is better suited for long-term dynamical stability in
principle.

## Equilibration

Equilibration should remove artefacts of the initial structure and velocity distribution before any
production statistics are collected. It does not have to be performed in `NVT`. A good equilibration
setup is usually the same as, or close to, the intended production ensemble: use `NVT` for
fixed-cell canonical sampling, `NPT_I` or `NPT_F` when pressure and cell relaxation are part of the
target workflow, and `LANGEVIN` when stochastic damping is desired. A typical protocol is:

- Prepare a reasonable starting structure;
- Set the target temperature directly for ordinary heating, or use staged heating/cooling when an
  abrupt temperature change is likely to destabilise the trajectory or when the thermal path is
  important;
- Equilibrate in the chosen target-like ensemble with suitable thermostat and, when needed, barostat
  time constants;
- Continue with adjusted thermostat or barostat parameters for gentler production sampling;
- Discard the equilibration part before computing averages.

For heating with `CSVR`, a practical two-stage strategy is to use a relatively short thermostat time
constant first, for example `TIMECON 50`, to bring the system close to the target temperature
quickly, and then switch to a gentler production setting such as `TIMECON 200`. This is useful when
the initial and target temperatures are not too far apart and the structure is already reasonable,
but a short equilibration stage is still desired before collecting statistics.

If a separate equilibration stage is not necessary, or if a simple robust setup is preferred, one
can also run the whole trajectory with an intermediate coupling strength, for example `TIMECON 100`.
This is often sufficient for ordinary room-temperature or moderately high-temperature simulations,
and the equilibration and production parts may then differ only in which initial frames are
discarded from the analysis.

If the target temperature is very high, directly heating to that temperature can make the simulation
likely to crash. In such cases, staged heating is recommended: increase
[TEMPERATURE](#CP2K_INPUT.MOTION.MD.TEMPERATURE) through several short equilibration stages until
the final target temperature is reached. This is especially useful for AIMD at several thousand
kelvin, reactive systems, poor initial structures, or systems with close contacts. The thermostat
`TIMECON` still controls the coupling strength within each stage, but it does not replace the need
for a gradual temperature schedule when an abrupt jump would destabilise the trajectory.

For cooling or annealing, it is usually best to lower the target temperature in stages, especially
when the final structure depends on the cooling path. One-step cooling can also be used, but the
thermostat coupling should be weak enough, with a sufficiently large `TIMECON` so that the
temperature decreases slowly rather than being quenched abruptly. The appropriate schedule depends
on the system and on whether the goal is sampling, relaxation, or deliberate annealing.

For stochastic equilibration, use the `LANGEVIN` ensemble and adjust
[GAMMA](#CP2K_INPUT.MOTION.MD.LANGEVIN.GAMMA): smaller values are closer to NVE-like dynamics and
equilibrate more slowly; larger values thermostat more aggressively.

Temperature is an instantaneous kinetic estimator and is expected to fluctuate. A single snapshot at
the target temperature is not evidence of equilibration. Check structural, energetic, and
property-specific observables. For heterogeneous systems, do not rely only on total energy or global
temperature; also inspect the region of interest.

## Trajectories, energies, and restarts

The trajectory is controlled by [MOTION/PRINT/TRAJECTORY](#CP2K_INPUT.MOTION.PRINT.TRAJECTORY). By
default, positions are written to files such as `project-pos-1.xyz`, where `project` is the value of
[PROJECT_NAME](#CP2K_INPUT.GLOBAL.PROJECT_NAME). Velocities, forces, cell vectors, stresses, and
restart files are controlled by other print keys in [MOTION/PRINT](#CP2K_INPUT.MOTION.PRINT).

The file `project-1.ener` is often the first file to inspect. It contains columns such as

```none
md_step, time[fs], e_kin[hartree], temp[K], e_pot[hartree], e_tot[hartree], elapsed_time_per_step[s]
```

For thermostatted or barostatted simulations, also inspect the corresponding conserved quantity or
cell output when available. For production runs, write restart files frequently enough that the
trajectory can be continued without losing significant sampling time.

## Validation checklist

Before running a long production trajectory, it is suggested to check that:

- The planned length of trajectory and reserved computational time is long enough to overcome
  (auto)correlations and/or observe phenomena;
- The force method and key parameters is validated on the particular system and is suitable for
  describing the target property;
- The structure is physically reasonable and the shortest interatomic distances are sensible;
- The periodic cell is large enough for the target property;
- The SCF procedure converges reliably on representative snapshots;
- The extrapolation method reduces SCF iterations without causing instability;
- The time step gives acceptable energy conservation in a short test;
- Thermostat and barostat time constants do not dominate the physical dynamics of interest;
- NPT cell fluctuations are reasonable and not dominated by small-cell oscillations;
- Output frequencies are sufficient for the target analysis but not so high that I/O dominates.