Biased molecular dynamics: Difference between revisions

From VASP Wiki
VisualWikitext
No edit summary
(Removed redirect to Category:Biased molecular dynamics)
Tag: Removed redirect
 
(9 intermediate revisions by 2 users not shown)
Line 1: Line 1:
''Biased molecular dynamics''' (MD) refers to advanced [[:Category:Molecular dynamics|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{{cite|torrie:jcp:1977}} and umbrella integration{{cite|kaestner:jcp:2005}}, to name a few. For a comprehensive description (especially about umbrella sampling), we refer the interested user to Ref. {{cite|frenkel:ap-book:2002}} written by D. Frenkel and B. Smit.


The probability density for a geometric parameter ξ of the system driven by a Hamiltonian:
The probability density for a geometric parameter ξ of the system driven by a Hamiltonian
:<math>
:<math>
H(q,p) = T(p) + V(q), \;
H(q,p) = T(p) + V(q), \;
</math>
</math>
with ''T''(''p''), and ''V''(''q'') being kinetic, and potential energies, respectively, can be written as:
with ''T''(''p''), and ''V''(''q'') being kinetic, and potential energies, respectively, can be written as
:<math>
:<math>
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} =
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} =
Line 11: Line 12:
The term <math>\langle X \rangle_H</math> stands for a thermal average of quantity ''X'' evaluated for the system driven by the Hamiltonian ''H''.
The term <math>\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>\tilde{V}(\xi)</math> acting on one or multiple selected internal coordinates of the system &xi;=&xi;(''q''), the Hamiltonian takes a form:
If the system is modified by adding a bias potential <math>\tilde{V}(\xi)</math> acting on one or multiple selected internal coordinates of the system &xi;=&xi;(''q''), the Hamiltonian takes the form
:<math>
:<math>
\tilde{H}(q,p) = H(q,p) + \tilde{V}(\xi),
\tilde{H}(q,p) = H(q,p) + \tilde{V}(\xi),
</math>
</math>
and the probability density of &xi; in the biased ensemble is:
and the probability density of &xi; in the biased ensemble is
:<math>
:<math>
\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}}
\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>
</math>
It can be shown that the biased and unbiased averages are related via a simple formula:
It can be shown that the biased and unbiased averages are related via
:<math>
:<math>
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}}}.
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>
</math>
More generally, an observable <math>\langle A \rangle_{H}</math>:
More generally, an observable
:<math>
:<math>
\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}
\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>
</math>
can be expressed in terms of thermal averages within the biased ensemble:
can be expressed in terms of thermal averages within the biased ensemble as
:<math>
:<math>
\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}}}.
\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>
</math>
Simulation methods such as umbrella sampling<ref name="Torrie77"/> use a bias potential to enhance sampling of &xi; in regions with low ''P''(&xi;<sub>''i''</sub>) such as transition regions of chemical
reactions.
The correct distributions are recovered afterwards using the equation for <math>\langle A \rangle_{H}</math> above.


A more detailed description of the method can be found in Ref.<ref name="FrenkelSmit"/>.
One of the most popular methods using bias potentials is umbrella sampling{{cite|torrie:jcp:1977}}. This method uses a bias potential to enhance sampling of &xi; in regions with low ''P''(&xi;<sub>''i''</sub>), e.g., transition regions of chemical reactions. The correct distributions are recovered afterward using the equation for <math>\langle A \rangle_{H}</math> above.
Biased molecular dynamics can be used also to introduce soft geometric constraints in which the controlled geometric parameter is not strictly constant, instead it oscillates in a narrow interval
of values.
The bias potentials are supported in both the NVT and NpT MD simulations regardless of the particular thermostat and/or barostat setting.
Different types of potentials can be freely combined.
[[File:Bias potentials.png|thumb|Graphical representation of (a) harmonic, (b) Fermi function-shaped, and (c) and Gauss function-shaped bias potentials.  ]]
 
== Supported types of bias potentials ==
Presently, the following types of bias potential are supported:
 
 
*sum of Harmonic potentials (see curve (a) on the plot above)
:<math>
\tilde{V}(\xi_1,\dots,\xi_{M_8}) = \sum_{\mu=1}^{M}\frac{1}{2}\kappa_{\mu} (\xi_{\mu}(q)-\xi_{0\mu})^2 \;
</math>
where the sum runs all (<math>M_8</math>) coordinates the potential acts upon, which are defined in the {{FILE|ICONST}}-file by setting the status to 8.
The parameters of the potential, <math>\kappa_{\mu}</math> (force constant) and <math>\xi_{0\mu}</math> (minimum of potential), are defined, respectively, via the keywords {{TAG|SPRING_K}} and {{TAG|SPRING_R0}} in {{FILE|INCAR}}.
Optionally, it is also possible to change the value of <math>\xi_{0\mu}</math> every MD step at a constant rate defined via parameter {{TAG|SPRING_V}}.  
The number of items defined via {{TAG|SPRING_K}}, {{TAG|SPRING_R0}}, and {{TAG|SPRING_V}} must be equal to <math>M_8</math>, otherwise the calculation terminates with an error message.
This form of bias potential is employed in several simulation protocols, such as the umbrella sampling<ref name="Torrie77"/>, umbrella integration, or steered MD, and is useful also in cases where the <math>\xi_{\mu}</math> values need restrained.
 
*sum of Fermi-like step functions (see curve (b) on the plot above)
:<math>
\tilde{V}(\xi_1,\dots,\xi_{M_4}) = \sum_{\mu=1}^{M_4}\frac{A_{\mu}}{1+\text{exp}\left [-D_{\mu}\frac{\xi(q)}{\xi_{0\mu}} -1 \right ]} \;
</math>
where the sum runs all (<math>M_4</math>) coordinates the potential acts upon, which are defined in the {{FILE|ICONST}}-file by setting the status to 4.
The parameters of the potential, <math>A_{\mu}</math> (the height of step), <math>D_{\mu}</math> (controlling the slope around the point <math>\xi_{0\mu}</math>), and  <math>\xi_{0\mu}</math> (position of the step), are defined, respectively,
via the keywords {{TAG|FBIAS_A}}, {{TAG|FBIAS_D}}, and {{TAG|FBIAS_R0}} in {{FILE|INCAR}}.
The number of items defined via {{TAG|FBIAS_A}}, {{TAG|FBIAS_D}}, and {{TAG|FBIAS_R0}} must be equal to <math>M_4</math>, otherwise the calculation terminates with an error message.
This form of potential is suitable especially for imposing restriction on the upper (or lower) limit of value of <math>\xi</math>.
 
*sum of Gauss functions (see curve (b) on the plot above)
:<math>
\tilde{V}(\xi_1,\dots,\xi_{M})  = \sum_{\nu=1}^{N_5}h_{\nu}\text{exp}\left [-\frac{\sum_{\mu=1}^{M_5}(\xi_{\mu}(q)-\xi_{0\nu,\mu})^2}{2w_{\nu}^2}  \right ], \;
</math>
where <math>N_5</math> is the number of Gaussian functions and
<math>M_5</math> is the number of coordinates the potential acts upon.
The latter defined in the {{FILE|ICONST}}-file by setting the status to 5.
The parameters of the potentials, <math>h_{\nu}</math> and <math>w_{\nu}</math> are defined in the {{FILE|PENALTYPOT}}-file.
The this type of bias potential is primarily intended for the use in metadynamics but since Gaussians can be used as basis functions for more general shapes, they can be used also to prepare various atypically-shaped bias potentials. 
 
== Examples of usage ==
 
 
== Andersen thermostat ==
 
* For a biased molecular dynamics run with Andersen thermostat, one has to:
#Set the standard MD-related tags: {{TAG|IBRION}}=0, {{TAG|TEBEG}}, {{TAG|POTIM}}, and {{TAG|NSW}}
#Set {{TAG|MDALGO}}=1 ({{TAG|MDALGO}}=11 in VASP 5.x), and choose an appropriate setting for {{TAG|ANDERSEN_PROB}}
#In order to avoid updating of the bias potential, set {{TAG|HILLS_BIN}}={{TAG|NSW}}
#Define collective variables in the {{FILE|ICONST}}-file, and set the <tt>STATUS</tt> parameter for the collective variables to 5
#Define the bias potential in the {{FILE|PENALTYPOT}}-file
 
== Nose-Hoover thermostat ==
 
* For a biased molecular dynamics run with Nose-Hoover thermostat, one has to:
#Set the standard MD-related tags: {{TAG|IBRION}}=0, {{TAG|TEBEG}}, {{TAG|POTIM}}, and {{TAG|NSW}}
#Set {{TAG|MDALGO}}=2 ({{TAG|MDALGO}}=21 in VASP 5.x), and choose an appropriate setting for {{TAG|SMASS}}
#In order to avoid updating of the bias potential, set {{TAG|HILLS_BIN}}={{TAG|NSW}}
#Define collective variables in the {{FILE|ICONST}}-file, and set the <tt>STATUS</tt> parameter for the collective variables to 5
#Define the bias potential in the {{FILE|PENALTYPOT}}-file
 
The values of all collective variables for each MD step are listed in the {{FILE|REPORT}}-file, check the lines after the string <tt>Metadynamics</tt>.


== References ==
== References ==
<references>
<references/>
<ref name="Torrie77">[http://dx.doi.org/10.1016/0021-9991(77)90121-8 G. M. Torrie and J. P. Valleau, J. Comp. Phys. 23, 187 (1977).]</ref>
<ref name="FrenkelSmit">D. Frenkel and B. Smit, ''Understanding molecular simulations: from algorithms to applications'', Academic Press: San Diego, 2002.</ref>
</references>
----


[[Category:Molecular dynamics]][[Category:Biased molecular dynamics]][[Category:Theory]][[Category:Howto]]
[[Category:Advanced molecular-dynamics sampling]][[Category:Theory]]

Latest revision as of 11:40, 16 October 2024

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[1] and umbrella integration[2], to name a few. For a comprehensive description (especially about umbrella sampling), we refer the interested user to Ref. [3] written by D. Frenkel and B. Smit.

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

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

The term 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 acting on one or multiple selected internal coordinates of the system ξ=ξ(q), the Hamiltonian takes the form

and the probability density of ξ in the biased ensemble is

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

More generally, an observable

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

One of the most popular methods using bias potentials is umbrella sampling[1]. 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 above.

References

  1. ↑ a b G. M. Torrie and J. P. Valleau, J. Comp. Phys. 23, 187 (1977).
  2. ↑ J. Kästner, and W. Thiel, Bridging the gap between thermodynamic integration and umbrella sampling provides a novel analysis method: “Umbrella integration”, J. Chem. Phys. 123, 144104 (2005).
  3. ↑ D. Frenkel and B. Smit, Understanding molecular simulations: from algorithms to applications, Academic Press: San Diego, 2002.
Retrieved from "https://www.vasp.at/wiki/index.php?title=Biased_molecular_dynamics&oldid=26636"