MDALGO

MDALGO = 0 | 1 | 2 | 3 | 11 | 21 | 13
Default: MDALGO = 0

Description: MDALGO specifies the molecular dynamics simulation protocol (in case IBRION=0 and VASP was compiled with -Dtbdyn).

MDALGO=0: Standard molecular dynamics

Standard molecular dynamics (IBRION=0), the same behavior as if VASP were compiled without -Dtbdyn.

MDALGO=1: Andersen thermostat

For the description of the Andersen thermostat see: Andersen thermostat.

The Andersen thermostat is only available for the NVT case.

Standard molecular dynamics in

For a standard molecular dynamics run with Anderson thermostat, one has to:

Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW.
Set MDALGO=1, and choose an appropriate setting for ANDERSEN_PROB

Constrained molecular dynamics

For a description of constrained molecular dynamics see Constrained molecular dynamics.

For a constrained molecular dynamics run with Andersen thermostat, one has to:

Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW
Set MDALGO=1, and choose an appropriate setting for ANDERSEN_PROB
Define geometric constraints in the ICONST-file, and set the STATUS parameter for the constrained coordinates to 0
When the free-energy gradient is to be computed, set LBLUEOUT=.TRUE.

Slow-growth approach

For a description of slow-growth approach see Slow-growth approach.

For a slow-growth simulation, one has to perform a calcualtion very similar to Constrained molecular dynamics but additionally the transformation velocity-related INCREM-tag for each geometric parameter with STATUS=0 has to be specified:

Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW
Set MDALGO=1, and choose an appropriate setting for ANDERSEN_PROB
Define geometric constraints in the ICONST-file, and set the STATUS parameter for the constrained coordinates to 0
When the free-energy gradient is to be computed, set LBLUEOUT=.TRUE.

Specify the transformation velocity-related INCREM-tag for each geometric parameter with STATUS=0.

Monitoring geometric parameters

Geometric parameters with STATUS = 7 in the ICONST-file are monitored during the MD simulation. The corresponding values are written onto the REPORT-file, for each MD step, after the lines following the string Monit_coord.

Sometimes it is desirable to terminate the simulation if all values of monitored parameters get larger that some predefined upper and/or lower limits. These limits can be set by the user by means of the VALUE_MAX and VALUE_MIN-tags.

To monitor geometric parameters during an MD run:

Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW
Set MDALGO=1, and choose an appropriate setting for ANDERSEN_PROB
Define geometric constraints in the ICONST-file, and set the STATUS parameter for the constrained coordinates to 7
Optionally, set the upper and/or lower limits for the coordinates, by means of the VALUE_MAX and VALUE_MIN-tags, respectively.

Special case: NVE ensemble

NVE ensemble calculations can be also run by selecting the Anderson thermostat and setting the update collision probability (ANDERSEN_PROB) to zero.

To run an NVE ensemble:

Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW.
Set MDALGO=1 and ANDERSEN_PROB=0.0.

MDALGO=2: Nose-Hoover thermostat

For the description of the Nose-Hoover thermostat see: Nose-Hoover thermostat.

The Nose-Hoover thermostat is currently only available for the NVT ensemble.

Standard molecular dynamics

For a standard molecular dynamics run with Nose-Hoover thermostat, one has to:

Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW.
Set MDALGO=2, and choose an appropriate setting for SMASS.

Constrained molecular dynamics

For a description of constrained molecular dynamics see Constrained molecular dynamics.

For a constrained molecular dynamics run with Nose-Hoover thermostat, one has to:

Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW.
Set MDALGO=2, and choose an appropriate setting for SMASS.
Define geometric constraints in the ICONST-file, and set the STATUS parameter for the constrained coordinates to 0.
When the free-energy gradient is to be computed, set LBLUEOUT=.TRUE.

Slow-growth approach

For a description of slow-growth approach see Slow-growth approach.

For a slow-growth approach run with Nose-Hoover thermostat, one has to:

Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW
Set MDALGO=2, and choose an appropriate setting for SMASS
Define geometric constraints in the ICONST-file, and set the STATUS parameter for the constrained coordinates to 0
When the free-energy gradient is to be computed, set LBLUEOUT=.TRUE.

Specify the transformation velocity-related INCREM-tag for each geometric parameter with STATUS=0

Monitoring geometric parameters

Geometric parameters with STATUS = 7 in the ICONST-file are monitored during the MD simulation. The corresponding values are written onto the REPORT-file, for each MD step, after the lines following the string Monit_coord.

Sometimes it is desirable to terminate the simulation if all values of monitored parameters get larger that some predefined upper and/or lower limits. These limits can be set by the user by means of the VALUE_MAX and VALUE_MIN-tags.

To monitor geometric parameters during an MD run:

Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW
Set MDALGO=2, and choose an appropriate setting for SMASS
Define geometric constraints in the ICONST-file, and set the STATUS parameter for the constrained coordinates to 7
Optionally, set the upper and/or lower limits for the coordinates, by means of the VALUE_MAX and VALUE_MIN-tags, respectively.

MDALGO=3: Langevin thermostat

Note: Geometric constraints and metadynamics are not available for Langevin dynamics.

NVT-simulation with Langevin thermostat

The Langevin thermostat^[1] maintains the temperature through a modification of Newton's equations of motion

{\displaystyle \dot{r_i} = p_i/m_i \qquad \dot{p_i} = F_i - {\gamma}_i\,p_i + f_i, }

where F_i is the force acting on atom i due to the interaction potential, γ_i is a friction coefficient, and f_i is a random force with dispersion σ_i related to γ_i through:

{\displaystyle \sigma_i^2 = 2\,m_i\,{\gamma}_i\,k_B\,T/{\Delta}t }

with Δt being the time-step used in the MD to integrate the equations of motion. Obviously, Langevin dynamics is identical to the classical Hamiltonian in the limit of vanishing γ.

To run an NVT-simulation with a Langevin thermostat, one has to:

Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW
Set ISIF=2
Set MDALGO=3 to invoke the Langevin thermostat
Specify friction coefficients for all species in the POSCAR file, by means of the LANGEVIN_GAMMA-tag.

NpT-simulation with Langevin thermostat

In the method of Parrinello and Rahman^[2]^[3], the equations of motion for ionic and lattice degrees-of-freedom are derived from the following Lagrangian:

{\displaystyle {\mathcal{L}}(s,h,\dot{s},\dot{h})=\frac{1}{2}\sum_i^N m_i \dot{s_i}^t\,G \dot{s_i}-V(s,h) +\frac{1}{2}W\,\mathrm{Tr}(\dot{h}^t \dot{h}) - p_\mathrm{ext}\Omega, }

where s_i is a position vector in fractional coordinates for atom i, h is the matrix formed by lattice vectors, the tensor G is defined as G=h^th, p_ext is the external pressure, Ω is the cell volume (Ω=det h), and W is a constant with the dimensionality of mass. Integrating equations of motion based on the Parrinello-Rahman Lagrangian generates an NpH ensemble, where the enthalpy $H=E+p_\mathrm{ext}\Omega$ is the constant of motion. The Parrinello-Rahman method can be combined with a Langevin thermostat^[1] to generate an NpT-ensemble.

To run an NpT-simulation (Parinello-Rahman dynamics) with a Langevin thermostat, one has to:

Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW
Set ISIF=3 to allow for relaxation of the cell volume and shape. At the moment, dynamics with fixed volume+variable shape (ISIF=4) or fixed shape+variable volume (ISIF=7) are not available.
Set MDALGO=3 to invoke the Langevin thermostat
Specify friction coefficients for all species in the POSCAR file, by means of the LANGEVIN_GAMMA-tag.
Specify a separate set of friction coefficient for the lattice degrees-of-freedom, using the LANGEVIN_GAMMA_L-tag.
Set a mass for the lattice degrees-of-freedom, using the PMASS-tag.
Optionally, one may define an external pressure (in kB), by means of the PSTRESS-tag.

The temperatures listed in the OSZICAR are computed using the kinetic energy including contribution from both atomic and lattice degrees of freedom. The external pressure for a simulation can be computed as one third of the trace of the stress-tensor corrected for kinetic contributions, listed in the OUTCAR file. This can be achieved, e.g. using the following command: grep "Total+kin" OUTCAR| awk 'BEGIN {p=0.} {p+=($2+$3+$4)/3.} END {print "pressure (kB):",p}'

Important: In Parinello-Rahman^[2]^[3] dynamics (NpT), the stress tensor is used to define forces on lattice degrees-of-freedom. In order to achieve a reasonable quality of sampling (and to avoid numerical problems), it is essential to eliminate Pulay stress. Unfortunately, this usually requires a rather large value of ENCUT. PREC=low, frequently used in NVT-MD is not recommended for molecular dynamics with variable cell volume.

Stochastic boundary conditions

In some cases it is desirable to study approach of initially non-equilibrium system to equilibrium. Examples of such simulations include the impact problems when a particle with large kinetic energy hits a surface or calculation of friction force between two surfaces sliding with respect to each other. As shown by Toton et al.^[4], this type of problems can be studied using the stochastic boundary conditions (SBC) derived from the generalized Langevin equation by Kantorovich and Rompotis.^[5] In this approach, the system of interest is divided into three regions: (a) fixed atoms, (b) the internal (Newtonian) atoms moving according to Newtonian dynamics, and (c) a buffer region of Langevin atoms (i.e., atoms governed by Langevin equations of motion) located between (a) and (b).

The role of the Langevin atoms is to dissipate heat, while the fixed atoms are needed solely to create the correct potential well for the Langevin atoms to move in. The Newtonian region should include all atoms relevant to the process under study: in the case of the impact problem, for instance, the Newtonian region should contain atoms of the molecule hitting the surface and several uppermost layers of the material forming the surface. Performing molecular dynamics with such a setup guarantees that the system (possibly out of equilibrium initially) arrives at the appropriate canonical distribution.

To run a simulation with stochastic boundary conditions, one has to:

Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW
Set ISIF=2
Set MDALGO=3 to invoke the Langevin thermostat
Prepare the POSCAR file in such a way that the Newtonian and Langevin atoms are treated as different species (even if they are chemically identical). In your POSCAR, use Selective Dynamics and the corresponding logical flags to define the frozen and moveable atoms.
Specify friction coefficients γ, for all species in the POSCAR file, by means of the LANGEVIN_GAMMA-tag: set the friction coefficients to 0 for all fixed and Newtonian atoms, and choose a proper γ for the Langevin atoms.

Practical example

Consider a system consisting of 16 hydrogen and 48 silicon atoms. Suppose that eight silicon atoms are considered to be Langevin atoms and the remaining 32 Si atoms are either fixed or Newtonian atoms. The Langevin and Newtonian (or fixed) atoms should be considered as different species, i.e., the POSCAR-file should contain the line like this:

Si H Si
40 16 8

As only the final eight Si atoms are considered to be Langevin atoms, the INCAR-file should contain the following line defining the friction coefficients:

LANGEVIN_GAMMA = 0.0   0.0   10.0

i.e., for all non-Langevin atoms, γ should be set to zero.

MDALGO=11: Biased-MD and metadynamics with Andersen thermostat

Andersen thermostat

For a short description of the Andersen thermostat see its section under MDALGO=1.

Metadynamics

For a description of metadynamics see Metadynamics.

For a metadynamics run with Andersen thermostat, one has to:

Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW
Set MDALGO=11, and choose an appropriate setting for ANDERSEN_PROB
Set the parameters HILLS_H, HILLS_W, and HILLS_BIN
Define collective variables in the ICONST-file, and set the STATUS parameter for the collective variables to 5
If needed, define the bias potential in the PENALTYPOT-file

The actual time-dependent bias potential is written to the HILLSPOT-file, which is updated after adding a new Gaussian. At the beginning of the simulation, VASP attempts to read the initial bias potential from the PENALTYPOT-file. For the continuation of a metadynamics run, copy HILLSPOT to PENALTYPOT. The values of all collective variables for each MD step are listed in REPORT-file, check the lines after the string Metadynamics.

Biased molecular dynamics

For a description of biased molecular dynamics see Biased molecular dynamics.

For a biased molecular dynamics run with Andersen thermostat, one has to:

Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW
Set MDALGO=11, and choose an appropriate setting for ANDERSEN_PROB
In order to avoid updating of the bias potential, set HILLS_BIN=NSW
Define collective variables in the ICONST-file, and set the STATUS parameter for the collective variables to 5
Define the bias potential in the PENALTYPOT-file

The values of all collective variables for each MD step are listed in the REPORT-file, check the lines after the string Metadynamics.

MDALGO=21: Biased-MD and metadynamics with Nose-Hoover Thermostat

Biased-molecular- or Metadynamics with Nose-Hoover Thermostat (SMASS needs to be specified in the INCAR file).

Metadynamics

For a description of metadynamics see Metadynamics.

For a metadynamics run with Nose-Hoover thermostat, one has to:

Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW
Set MDALGO=21, and choose an appropriate setting for SMASS
Set the parameters HILLS_H, HILLS_W, and HILLS_BIN
Define collective variables in the ICONST-file, and set the STATUS parameter for the collective variables to 5
If needed, define the bias potential in the PENALTYPOT-file

The actual time-dependent bias potential is written to the HILLSPOT-file, which is updated after adding a new Gaussian. At the beginning of the simulation, VASP attempts to read the initial bias potential from the PENALTYPOT-file. For the continuation of a metadynamics run, copy HILLSPOT to PENALTYPOT. The values of all collective variables for each MD step are listed in REPORT-file, check the lines after the string Metadynamics.

Biased molecular dynamics

For a description of biased molecular dynamics see Biased molecular dynamics.

For a biased molecular dynamics run with Nose-Hoover thermostat, one has to:

Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW
Set MDALGO=21, and choose an appropriate setting for SMASS
In order to avoid updating of the bias potential, set HILLS_BIN=NSW
Define collective variables in the ICONST-file, and set the STATUS parameter for the collective variables to 5
Define the bias potential in the PENALTYPOT-file

The values of all collective variables for each MD step are listed in the REPORT-file, check the lines after the string Metadynamics.

MDALGO=13: Multiple Anderson thermostats

Up to three user-defined atomic subsystems may be coupled with independent Andersen thermostats^[6] (see remarks under MDALGO=1 as well). The POSCAR file must be organized such that the positions of atoms of subsystem i+1 are defined after those for the subsystem i, and the following flags must be set by the user:

NSUBSYS=[int array]

Define the last atom for each subsystem (two or three values must be supplied). For instance, if total of 20 atoms is defined in the POSCAR file, and the initial 10 atoms belong to the subsystem 1, the next 7 atoms to the subsystem 2, and the last 3 atoms to the subsystem 3, NSUBSYS should be defined as follows:

NSUBSYS= 10 17 20

Note that the last number in the previous example is actually redundant (clearly the last three atoms belong to the last subsystem) and does not have to be user-supplied.

TSUBSYS=[real array]

Simulation temperature for each subsystem

PSUBSYS=[real array]

Collision probability for atoms in each subsystem. Only the values 0≤PSUBSYS≤1 are allowed.

Related Tags and Sections

IBRION, ISIF, SMASS, ANDERSEN_PROB, RANDOM_SEED, LBLUEOUT, SHAKETOL, SHAKEMAXITER, HILLS_H, HILLS_W, HILLS_BIN, INCREM, VALUE_MIN, VALUE_MAX, LANGEVIN_GAMMA, LANGEVIN_GAMMA_L, PMASS, NSUBSYS, TSUBSYS, PSUBSYS, ICONST, PENALTYPOT, HILLSPOT, REPORT

Examples that use this tag

References

↑ ^a ^b M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Oxford university press: New York, 1991.
↑ ^a ^b M. Parrinello and A. Rahman, Phys. Rev. Lett. 45, 1196 (1980).
↑ ^a ^b M. Parrinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981).
↑ D. Toton, C. D. Lorenz, N. Rompotis, N. Martsinovich, and L. Kantorovich, J. Phys.: Condens. Matter 22, 074205 (2010).
↑ L. Kantorovich and N. Rompotis, Phys. Rev. B 78, 094305 (2008).
↑ H. C. Andersen, J. Chem. Phys. 72, 2384 (1980).

[Allen91-1] M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Oxford university press: New York, 1991.

[Parrinello80-2] M. Parrinello and A. Rahman, Phys. Rev. Lett. 45, 1196 (1980).

[Parrinello81-3] M. Parrinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981).

[Toton10-4] D. Toton, C. D. Lorenz, N. Rompotis, N. Martsinovich, and L. Kantorovich, J. Phys.: Condens. Matter 22, 074205 (2010).

[Kantorovich08-5] L. Kantorovich and N. Rompotis, Phys. Rev. B 78, 094305 (2008).

[Andersen80-6] H. C. Andersen, J. Chem. Phys. 72, 2384 (1980).

[1]

[2]

[3]

[4]

[5]

[6]