Ab Initio Molecular Dynamics

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)

Further, it can be desirable to change the smearing of your system (ISMEAR) to the Fermi Smearing, since this explicitly considered the temperature applied to the system also for the electrons. The occupation numbers at the Fermi edge are then modeled with a Fermi-Dirac distribution at the MD temperature (Phyiscally correct), instead of a somewhat arbitrary Gaussian function. For Fermi smearing, the smearing width SIGMA needs to be set to resemble the temperature, which can be done with the formula:

 SIGMA = 8.61733333E-5 * T(K)

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.

To monitor the progress and stability of a AIMD simulation, the check_geoopt.py script can be used.