Category:Biased molecular dynamics: Difference between revisions

From VASP Wiki
No edit summary
((move edit by Karsai to main))
Line 1: Line 1:
'''Biased molecular dynamics''' (MD) refers to advanced [[:Category:Molecular dynamics|MD-simulation methods]] that introduce a ''bias potential''.
'''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. In the following, a short theoretical overview, an introduction to the available bias potentials, and a description of how to use bias potentials in VASP is given. 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.
 
For instance, umbrella sampling{{cite|torrie:jcp:1977}} uses a bias potential to enhance the sampling of a geometric parameter in regions with low probability density, e.g., transition regions of chemical reactions.
The correct distributions are recovered afterward using an exact relation between the biased and unbiased systems. D. Frenkel and B. Smit provide a more detailed description of the method in Ref. {{cite|frenkel:ap-book:2002}}.
'''Biased MD''' can also be used to introduce soft geometric constraints in which the controlled geometric parameter is not strictly constant. Instead, it oscillates in a narrow interval of values.
 
In VASP, the bias potentials are supported in both the [[NVT ensemble|NVT]] and [[NpT ensemble|NpT]] MD simulations, regardless of the particular [[:Category:Thermostats|thermostat]] and/or [[PSTRESS|barostat]] setting.
Different types of potentials can be freely combined. The values of all collective variables defined in the {{FILE|ICONST}} file for each MD step are listed in the {{FILE|REPORT}} file. Check the lines after the string <tt>Metadynamics</tt>.
[[File:Bias potentials.png|thumb|Graphical representation of (a) harmonic, (b) Fermi function-shaped, and (c) and Gauss function-shaped bias potentials. ]]


== Theoretical background ==
== Theoretical background ==


The probability density for a geometric parameter &xi; of the system driven by a Hamiltonian:
The probability density for a geometric parameter &xi; 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 22: Line 14:
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, e.g., umbrella sampling{{cite|torrie:jcp:1977}}, use 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. Ref.{{cite|frenkel:ap-book:2002}} provides a more detailed description of the method.


== Supported types of bias potentials ==
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.
Presently, the following types of bias potential are supported:
 
== How to ==
*First one needs to setup a [[:Category:Molecular dynamics|standard molecular dynamics]] run. The bias potentials are supported in both the [[NVT ensemble|NVT]] and [[NpT ensemble|NpT]] MD simulations, regardless of the particular [[:Category:Thermostats|thermostat]] and/or [[PSTRESS|barostat]] setting.
{{NB|mind|Mind that for VASP 5.x the biased molecular dynamics runs have to be chosen by adding 10 to the chosen value of {{TAG|MDALGO}}. E.g. {{TAG|MDALGO}}{{=}}12 instead of {{TAG|MDALGO}}{{=}}2 has to be chosen for Nose-Hoover thermostat.|:}}
*Then one needs to set the geometric parameters and bias potential types. The geometric parameters &xi;, also called collective variables, that are controlled via the potentials are defined in the {{TAG|ICONST}} file. The type of bias potential is also set in this file. The format of this file follows this layout:
flag item(1) ... item(N) status
Here <code>flag</code> specifies the type of geometric parameters (bond lengths, angles, etc.), <code>item(1) ... item(N)</code> the actual geometric items (atom numbers, etc.) and  the <code>status</code> sets the type of used bias potential. [[ICONST|Here]] we advise the user to also look into the description of the full capabilities of the geometric parameters.
*To avoid the update of the bias potential {{TAG|HILLS_BIN}}={{TAG|NSW}} is set by default.
*Next one needs to set the parameters of the potential in the [[INCAR|INCAR]] (harmonic and step function) file or in the [[PENALTYPOT|PENALTYPOT]] file (Gaussian). The following table summarizes the available potentials with their corresponding parameters:
:{| cellpadding="3" cellspacing="0" border="1"
|colspan="1" style="text-align: center;"|Potential type ||colspan="1" style="text-align: center;"| <code>status</code> in [[ICONST|ICONST]] || colspan="1" style="text-align: center;"| [[INCAR|INCAR]] parameters
|-
|colspan="1" style="text-align: center;"|Harmonic potential || colspan="1" style="text-align: center;"| 8 || colspan="1" style="text-align: center;"|{{TAG|SPRING_K}}, {{TAG|SPRING_R0}}, optional {{TAG|SPRING_V0}}
|-
|colspan="1" style="text-align: center;"|Step function || colspan="1" style="text-align: center;"| 4 || {{TAG|FBIAS_A}}, {{TAG|FBIAS_D}}, and {{TAG|FBIAS_R0}}
|-
|colspan="1" style="text-align: center;"|Gaussian potential || colspan="1" style="text-align: center;"| 5 || parameters set in [[PENALTYPOT|PENALTYPOT]] file
|}
 
=== Available potentials ===
[[File:Bias potentials.png|thumb| Fig.1) Graphical representation of (a) harmonic, (b) Fermi function-shaped, and (c) and Gauss function-shaped bias potentials.  ]]
In the following details are given on how to control the available bias potentials in VASP that are plotted in Fig.1.


=== Harmonic potentials ===
==== Harmonic potentials ====


:A sum of Harmonic potentials (see curve (a) on the plot above)
:A sum of Harmonic potentials (curve (a) in Fig.1)
::<math>
::<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 \;
\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>
</math>
:where the sum runs over all (<math>M_8</math>) coordinates the potential acts upon. The coordinates are defined in the {{FILE|ICONST}} file by setting the <code>status</code> to 8. The parameters of the potential are the force constant <math>\kappa_{\mu}</math> ({{TAG|SPRING_K}}) and the minimum or the potential <math>\xi_{0\mu}</math> {{TAG|SPRING_R0}}. These must be set in the {{FILE|INCAR}} file. Optionally, it is possible to change the value of <math>\xi_{0\mu}</math> every MD step at a constant rate defined via the {{FILE|INCAR}} tag {{TAG|SPRING_V0}}.  
:where the sum runs over all (<math>M_8</math>) coordinates the potential acts upon. The potential is chosen in the {{FILE|ICONST}} file by setting the <code>status</code> to 8. The parameters of the potential are the force constant <math>\kappa_{\mu}</math> ({{TAG|SPRING_K}}) and the minimum or the potential <math>\xi_{0\mu}</math> {{TAG|SPRING_R0}}. These must be set in the {{FILE|INCAR}} file. Optionally, it is possible to change the value of <math>\xi_{0\mu}</math> every MD step at a constant rate defined via the {{FILE|INCAR}} tag {{TAG|SPRING_V0}}.  
{{NB|mind|The number of items defined via {{TAG|SPRING_K}}, {{TAG|SPRING_R0}}, and {{TAG|SPRING_V0}} must be equal to <math>M_8</math>. Otherwise, the calculation terminates with an error message.|:}}  
{{NB|mind|The number of items defined via {{TAG|SPRING_K}}, {{TAG|SPRING_R0}}, and {{TAG|SPRING_V0}} 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{{cite|torrie:jcp:1977}}, umbrella integration, or steered MD, and is useful also in cases where the <math>\xi_{\mu}</math> values need to be restrained.  
:This form of bias potential is employed in several simulation protocols, such as the umbrella sampling{{cite|torrie:jcp:1977}}, umbrella integration, or steered MD, and is useful also in cases where the <math>\xi_{\mu}</math> values need to be restrained.


=== Step function ===
==== Step function ====


:A sum of Fermi-like step functions (see curve (b) on the plot above)
:A sum of Fermi-like step functions (curve (b) in Fig.1)
::<math>
::<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 ]} \;
\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>
</math>
:where the sum runs over all (<math>M_4</math>) coordinates the potential acts upon. The coordinates are defined in the {{FILE|ICONST}} file by setting the <code>status</code> to 4. The parameters of the potential are the height of the step (<math>A_{\mu}</math> set by {{TAG|FBIAS_A}}), the slope around the point <math>\xi_{0\mu}</math> (<math>D_{\mu}</math> set by {{TAG|FBIAS_D}}), and the position of the step (<math>\xi_{0\mu}</math> set by {{TAG|FBIAS_R0}}). These must be set in the {{FILE|INCAR}} file.  
:where the sum runs over all (<math>M_4</math>) coordinates the potential acts upon. The potential is chosen in the {{FILE|ICONST}} file by setting the <code>status</code> to 4. The parameters of the potential are the height of the step (<math>A_{\mu}</math> set by {{TAG|FBIAS_A}}), the slope around the point <math>\xi_{0\mu}</math> (<math>D_{\mu}</math> set by {{TAG|FBIAS_D}}), and the position of the step (<math>\xi_{0\mu}</math> set by {{TAG|FBIAS_R0}}). These must be set in the {{FILE|INCAR}} file.  
{{NB|mind|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.|:}}  
{{NB|mind|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 restrictions on the upper (or lower) limit of the value of <math>\xi</math>.  
:This form of potential is suitable especially for imposing restrictions on the upper (or lower) limit of the value of <math>\xi</math>.


=== Gaussian potential ===
==== Gaussian potential ====


:A sum of Gauss functions (see curve (b) on the plot above)
:A sum of Gauss functions (curve (b) in Fig.1)
::<math>
::<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 ], \;
\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>
</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 coordinates are defined in the {{FILE|ICONST}} file by setting the <code>status</code> to 5. The parameters of the potentials, <math>h_{\nu}</math>, <math>w_{\nu}</math>, and <math>\xi_{0\nu,\mu}</math> are defined in the {{FILE|PENALTYPOT}} file.  
: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 potential is chosen in the {{FILE|ICONST}} file by setting the <code>status</code> to 5. The parameters of the potentials, <math>h_{\nu}</math>, <math>w_{\nu}</math>, and <math>\xi_{0\nu,\mu}</math> are defined in the {{FILE|PENALTYPOT}} file.  


:This type of bias potential is primarily intended for use in metadynamics, but since Gaussians can be used as basis functions for more general shapes, they can also be used to prepare various atypically-shaped bias potentials.
:This type of bias potential is primarily intended for use in metadynamics, but since Gaussians can be used as basis functions for more general shapes, they can also be used to prepare various atypically shaped bias potentials.
 
=== Output ===
The values of all collective variables defined in the {{FILE|ICONST}} file for each MD step are listed in the {{FILE|REPORT}} file. Check the lines after the string <tt>Metadynamics</tt>.


== Examples of usage ==
== Examples of usage ==
Line 102: Line 115:
  R 1 5 4   
  R 1 5 4   


Next, we need to specify the bias potential parameters {{TAG|FBIAS_A}}, {{TAG|FBIAS_D}}, and {{TAG|FBIAS_R0}} in the {{FILE|INCAR}} file. Since the bias potential acts only on one internal coordinate (<math>M_4=1</math>), we need to provide only one number for each of the tags. The setting
Next, we need to set the molecular dynamics parameters and specify the bias potential parameters {{TAG|FBIAS_A}}, {{TAG|FBIAS_D}}, and {{TAG|FBIAS_R0}} in the {{FILE|INCAR}} file:
# Molecular dynamics part
{{TAGBL|IBRION}} = 0
{{TAGBL|TEBEG}} = 300
{{TAGBL|TEEND}} = 300
{{TAGBL|MDALGO}} = 2
{{TAGBL|POTIM}} = 2.0
{{TAGBL|NSW}} = 10000
# Bias potential part
{{TAGBL|FBIAS_A}}  = 1
{{TAGBL|FBIAS_D}}  = 50
{{TAGBL|FBIAS_R0}} = 3.5


FBIAS_A  = 1
Since the bias potential acts only on one internal coordinate (<math>M_4=1</math>), we need to provide only one number for each of the tags. The chosen bias potential parameters  
  FBIAS_D  = 50
ensure that repulsive bias forces steeply increase when the C...Cl distance is increased beyond about <math>3.2 \AA</math>. This causes a shortening of the distance in the next MD step. Notice that the bias force is essentially negligible for distances below <math>3 \AA</math>.  
FBIAS_R0 = 3.5
 
ensures that repulsive bias forces steeply increase when the C...Cl distance is increased beyond about <math>3.2 \AA</math>. This causes a shortening of the distance in the next MD step. Notice that the bias force is essentially negligible for distances below <math>3 \AA</math>.  
A careful adjustment of {{TAG|FBIAS_A}} and {{TAG|FBIAS_D}} is needed to ensure that (i) the bias force is large enough to effectively limit the value of <math>\xi</math>, and (ii) the interval of <math>\xi</math> values for which the bias forces are significant is broad enough to avoid overcoming via random fluctuations.       
A careful adjustment of {{TAG|FBIAS_A}} and {{TAG|FBIAS_D}} is needed to ensure that (i) the bias force is large enough to effectively limit the value of <math>\xi</math>, and (ii) the interval of <math>\xi</math> values for which the bias forces are significant is broad enough to avoid overcoming via random fluctuations.       
A suitable setting can be found by noting that the maximal bias force of <math>\frac{D\,A}{4\xi_0}</math> is exerted on the system at the point <math>\xi = \xi_{0}</math>. This can be seen by inspecting the analytical expression for the potential.
A suitable setting can be found by noting that the maximal bias force of <math>\frac{D\,A}{4\xi_0}</math> is exerted on the system at the point <math>\xi = \xi_{0}</math>. This can be seen by inspecting the analytical expression for the potential.


== How to ==
== Related methods in VASP ==
 
*[[:Category:Metadynamics|Metadynamics]]: In contrast to the methods discussed on this page metadynamics continuously updates the bias potential of the system to push it into unvisited parts of phase space.
=== Biased MD run using the Andersen thermostat and a Gaussian bias potential ===
 
:For a biased MD 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}}
:#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
 
=== Biased MD run using the Nose-Hoover thermostat and a Gaussian bias potential ===
 
:For a biased MD 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}}
:#To avoid updating 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>.
*[[:Category:Interface pinning|Interface pinning]]: This employs a bias potential to pin the state of an interface between a solid and a liquid. This method uses entirely different {{TAG|INCAR}} tags than the bias potentials presented on this page.


== References ==
== References ==
<references/>
<references/>
----


[[Category:VASP|Biased molecular dynamics]][[Category:Molecular dynamics]][[Category:Theory]][[Category:Howto]]
[[Category:VASP|Biased molecular dynamics]][[Category:Molecular dynamics]]

Revision as of 07:57, 20 October 2023

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. In the following, a short theoretical overview, an introduction to the available bias potentials, and a description of how to use bias potentials in VASP is given. For a comprehensive description (especially about umbrella sampling), we refer the interested user to Ref. [3] written by D. Frenkel and B. Smit.

Theoretical background

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 Pi), e.g., transition regions of chemical reactions. The correct distributions are recovered afterward using the equation for above.

How to

Mind: Mind that for VASP 5.x the biased molecular dynamics runs have to be chosen by adding 10 to the chosen value of MDALGO. E.g. MDALGO=12 instead of MDALGO=2 has to be chosen for Nose-Hoover thermostat.
  • Then one needs to set the geometric parameters and bias potential types. The geometric parameters ξ, also called collective variables, that are controlled via the potentials are defined in the ICONST file. The type of bias potential is also set in this file. The format of this file follows this layout:
flag item(1) ... item(N) status

Here flag specifies the type of geometric parameters (bond lengths, angles, etc.), item(1) ... item(N) the actual geometric items (atom numbers, etc.) and the status sets the type of used bias potential. Here we advise the user to also look into the description of the full capabilities of the geometric parameters.

  • To avoid the update of the bias potential HILLS_BIN=NSW is set by default.
  • Next one needs to set the parameters of the potential in the INCAR (harmonic and step function) file or in the PENALTYPOT file (Gaussian). The following table summarizes the available potentials with their corresponding parameters:
Potential type status in ICONST INCAR parameters
Harmonic potential 8 SPRING_K, SPRING_R0, optional SPRING_V0
Step function 4 FBIAS_A, FBIAS_D, and FBIAS_R0
Gaussian potential 5 parameters set in PENALTYPOT file

Available potentials

Fig.1) Graphical representation of (a) harmonic, (b) Fermi function-shaped, and (c) and Gauss function-shaped bias potentials.

In the following details are given on how to control the available bias potentials in VASP that are plotted in Fig.1.

Harmonic potentials

A sum of Harmonic potentials (curve (a) in Fig.1)
where the sum runs over all () coordinates the potential acts upon. The potential is chosen in the ICONST file by setting the status to 8. The parameters of the potential are the force constant (SPRING_K) and the minimum or the potential SPRING_R0. These must be set in the INCAR file. Optionally, it is possible to change the value of every MD step at a constant rate defined via the INCAR tag SPRING_V0.
Mind: The number of items defined via SPRING_K, SPRING_R0, and SPRING_V0 must be equal to . Otherwise, the calculation terminates with an error message.
This form of bias potential is employed in several simulation protocols, such as the umbrella sampling[1], umbrella integration, or steered MD, and is useful also in cases where the values need to be restrained.

Step function

A sum of Fermi-like step functions (curve (b) in Fig.1)
where the sum runs over all () coordinates the potential acts upon. The potential is chosen in the ICONST file by setting the status to 4. The parameters of the potential are the height of the step ( set by FBIAS_A), the slope around the point ( set by FBIAS_D), and the position of the step ( set by FBIAS_R0). These must be set in the INCAR file.
Mind: The number of items defined via FBIAS_A, FBIAS_D, and FBIAS_R0 must be equal to . Otherwise, the calculation terminates with an error message.
This form of potential is suitable especially for imposing restrictions on the upper (or lower) limit of the value of .

Gaussian potential

A sum of Gauss functions (curve (b) in Fig.1)
where is the number of Gaussian functions and is the number of coordinates the potential acts upon. The potential is chosen in the ICONST file by setting the status to 5. The parameters of the potentials, , , and are defined in the PENALTYPOT file.
This type of bias potential is primarily intended for use in metadynamics, but since Gaussians can be used as basis functions for more general shapes, they can also be used to prepare various atypically shaped bias potentials.

Output

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

Examples of usage

Let us consider the nucleophile substitution reaction of CHCl with Cl. The reactant is a weak van-der-Waals complex. The corresponding POSCAR file reads

vdW complex CH3Cl...Cl 
1.00000000000000
12.0000000000000000    0.0000000000000000    0.0000000000000000
0.0000000000000000    12.0000000000000000    0.0000000000000000
0.0000000000000000    0.0000000000000000    12.0000000000000000
C H Cl
1 3 2
cart
5.91331371  7.11364924  5.78037960
5.81982231  8.15982106  5.46969017
4.92222130  6.65954232  5.88978969
6.47810398  7.03808479  6.71586385
4.32824726  8.75151396  7.80743202
6.84157897  6.18713289  4.46842049

Due to the weak interactions between CHCl and Cl, the complex can collapse at high temperatures. This can be avoided by setting an upper bound for the length of the non-bonding Cl...C interactions. This can be conveniently achieved by using a Fermi-like step-shaped bias potential. To this end, we need to define the Cl...C distance, i.e., the distance between the atoms 1 and 5, as a coordinate with status 4 in the ICONST file:

R 1 5 4   

Next, we need to set the molecular dynamics parameters and specify the bias potential parameters FBIAS_A, FBIAS_D, and FBIAS_R0 in the INCAR file:

# Molecular dynamics part
IBRION = 0
TEBEG = 300
TEEND = 300
MDALGO = 2
POTIM = 2.0
NSW = 10000
# Bias potential part
FBIAS_A  = 1
FBIAS_D  = 50
FBIAS_R0 = 3.5

Since the bias potential acts only on one internal coordinate (), we need to provide only one number for each of the tags. The chosen bias potential parameters ensure that repulsive bias forces steeply increase when the C...Cl distance is increased beyond about . This causes a shortening of the distance in the next MD step. Notice that the bias force is essentially negligible for distances below . A careful adjustment of FBIAS_A and FBIAS_D is needed to ensure that (i) the bias force is large enough to effectively limit the value of , and (ii) the interval of values for which the bias forces are significant is broad enough to avoid overcoming via random fluctuations. A suitable setting can be found by noting that the maximal bias force of is exerted on the system at the point . This can be seen by inspecting the analytical expression for the potential.

Related methods in VASP

  • Metadynamics: In contrast to the methods discussed on this page metadynamics continuously updates the bias potential of the system to push it into unvisited parts of phase space.
  • Interface pinning: This employs a bias potential to pin the state of an interface between a solid and a liquid. This method uses entirely different INCAR tags than the bias potentials presented on this page.

References

This category currently contains no pages or media.