Skip to content

Ab Initio Molecular Dynamics

Julien Steffen edited this page Mar 11, 2024 · 2 revisions

Direct calculation of time-dependent properties with DFT become significantly expensive if the system contains more than 50-100 atoms and if more than some 1000 of steps shall be calculated. If you do the calculations on a cluster with a walltime limit, the script md_long.sh can be used to manage arbitrary long trajectories.

Build on the settings in the general input section, the following additions and changes should be made in the INCAR file:

  • IBRION = 0 Activates the MD mode
  • MDALGO = [number] If a NVT ensemble with constant volume shall be simulated, two options are possible. 1 activates the stochastic Andersen thermostat, 2 activates the deterministic Nose-Hoover thermostat. For the simulation of a NpT ensemble with constant pressure and variable volume, 3 activates the Langevin thermostat/barostat. In all cases, the initial velocities are generated by chance (accorting to the Maxwell-Boltzmann velocity distribution). If the POSCAR has been obtained by copying a CONTCAR of a previous calculation, the velocities are read in from this file (below the coordinate section).
  • NSW = [number] The number of time steps.
  • POTIM = [value] The MD time step in fs. If only heavy elements are in the system, time steps between 2 and 10 fs might be ok (depending on how heavy the lightest element is). If hydrogen is present as well, the time step should not be larger than 1 fs. Alternatively, the mass of the lightest atoms can be altered, if pure conformational sampling is desired (in the POTCAR, see
    here.)
  • TEBEG = [value] The initial temperature in K.
  • TEEND = [value] The final temperature in K. If TEBEG is set to a different value, the temperature will vary linearly during the dynamics from the start to the end.
  • SMASS = 0 If the Nose-Hoover thermostat is used, this sets the mass of the extra degree of freedom, regulating the stiffness around which the actual temperature varies around the desired temperature. The default value 0 sets the mass to 40 time steps, which should be good in most cases.
  • IWAVPR = 12 How the new orbitals are predicted from the old ones. For MD calculations, VASP strongly recommend 12. Might lead to problems for machine learning calculations, see [section](Machine Learning).
  • LMAXMIX = 4 PAW charge densities from up to d-electrons are passed through the Broyden charge-density mixer. Should be OK for most situations.
  • MAXMIX = 100 Number of previous steps stored in the Broyden charge mixer for subsequent MD steps.
  • LANGEVIN_GAMMA = [values] Friction coefficients for the Langevin thermostat in NpT calculations, one value per element. Values around 5.0 are good in most cases, e.g., LANGEVIN_GAMMA = 5.0 5.0.
  • LANGEVIN_GAMMA_L = [value] Frictition coefficient for lattice degrees of freedom, 5.0 should be OK.
  • PMASS = 1000 Fictitious ass of the lattice degree of freedom (1000 should be a reasonable value)

For NpT dynamics it is usually desirable to keep the shape of the simulation cell. This can be achieved by setting ISIF = 8 or with the following ICONST file:

LA 1 2 0
LA 1 3 0
LA 2 3 0
LR 1 0
LR 2 0
LR 3 0
S  1  0  0  0  0  0 0
S  0  1  0  0  0  0 0
S  0  0  1  0  0  0 0
S  0  0  0  1 -1  0 0
S  0  0  0  1  0 -1 0
S  0  0  0  0  1 -1 0

For NpT simulations of large systems, it might occur that the application of the simulation box contraints leads to problems, resulting in errors like: Error: RATTLE_vel algorithm did not converge!. In this case, vary LANGEVIN_GAMMA, LANGEVIN_GAMMA_L and PMASS until it works.