Biased molecular dynamics

Biased molecular dynamics' (MD) refers to advanced MD-simulation methods that introduce a bias potential. One of the most important purposes of using bias potentials is to enhance the sampling of phase space with low probability density (e.g., transition regions of chemical reactions). Depending on the type of sampling and in combination with the corresponding statistical methods one then has access to important thermodynamic quantities like, e.g., free energies. Biased molecular dynamics comes in very different flavors such as, e.g., umbrella sampling and umbrella integration, to name a few. For a comprehensive description (especially about umbrella sampling), we refer the interested user to Ref. written by D. Frenkel and B. Smit.

The probability density for a geometric parameter ξ of the system driven by a Hamiltonian

[math]\displaystyle{ H(q,p) = T(p) + V(q), \; }[/math]

with T(p), and V(q) being kinetic, and potential energies, respectively, can be written as

[math]\displaystyle{ P(\xi_i)=\frac{\int\delta\Big(\xi(q)-\xi_i\Big) \exp\left\{-H(q,p)/k_B\,T\right\} dq\,dp}{\int \exp\left\{-H(q,p)/k_B\,T\right\}dq\,dp} = \langle\delta\Big(\xi(q)-\xi_i\Big)\rangle_{H}. }[/math]

The term [math]\displaystyle{ \langle X \rangle_H }[/math] stands for a thermal average of quantity X evaluated for the system driven by the Hamiltonian H.

If the system is modified by adding a bias potential [math]\displaystyle{ \tilde{V}(\xi) }[/math] acting on one or multiple selected internal coordinates of the system ξ=ξ(q), the Hamiltonian takes the form

[math]\displaystyle{ \tilde{H}(q,p) = H(q,p) + \tilde{V}(\xi), }[/math]

and the probability density of ξ in the biased ensemble is

[math]\displaystyle{ \tilde{P}(\xi_i)= \frac{\int \delta\Big(\xi(q)-\xi_i\Big) \exp\left\{-\tilde{H}(q,p)/k_B\,T\right\} dq\,dp}{\int \exp\left\{-\tilde{H}(q,p)/k_B\,T\right\}dq\,dp} = \langle\delta\Big(\xi(q)-\xi_i\Big)\rangle_{\tilde{H}}. }[/math]

It can be shown that the biased and unbiased averages are related via

[math]\displaystyle{ P(\xi_i)=\tilde{P}(\xi_i) \frac{\exp\left\{\tilde{V}(\xi)/k_B\,T\right\}}{\langle \exp\left\{\tilde{V}(\xi)/k_B\,T\right\} \rangle_{\tilde{H}}}. }[/math]

More generally, an observable

[math]\displaystyle{ \langle A \rangle_{H} = \frac{\int A(q) \exp\left\{-H(q,p)/k_B\,T\right\} dq\,dp}{\int \exp\left\{-H(q,p)/k_B\,T\right\}dq\,dp} }[/math]

can be expressed in terms of thermal averages within the biased ensemble as

[math]\displaystyle{ \langle A \rangle_{H} =\frac{\langle A(q) \,\exp\left\{\tilde{V}(\xi)/k_B\,T\right\} \rangle_{\tilde{H}}}{\langle \exp\left\{\tilde{V}(\xi)/k_B\,T\right\} \rangle_{\tilde{H}}}. }[/math]

One of the most popular methods using bias potentials is umbrella sampling. This method uses a bias potential to enhance sampling of ξ in regions with low P(ξ_i), e.g., transition regions of chemical reactions. The correct distributions are recovered afterward using the equation for [math]\displaystyle{ \langle A \rangle_{H} }[/math] above.

How to

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

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

1. Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW.
2. Choose thermostat:
   1. Set MDALGO=1 (or MDALGO=11 in VASP 5.x), and choose an appropriate setting for ANDERSEN_PROB.
   2. Set MDALGO=2 (or MDALGO=21 in VASP 5.x), and choose an appropriate setting for SMASS.
3. In order to avoid updating of the bias potential, set HILLS_BIN=NSW.
4. Define collective variables in the ICONST-file, and set the STATUS parameter for the collective variables to 5.
5. Define the bias potential in the PENALTYPOT file if necessary.

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