Biased molecular dynamics: Difference between revisions

Revision as of 10:01, 6 April 2023

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

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): H(q,p)=T(p)+V(q),\;

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

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): 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}}.

The term Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): \langle X\rangle _{H} 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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): {\tilde {V}}(\xi ) acting on one or multiple selected internal coordinates of the system ξ=ξ(q), the Hamiltonian takes a form:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): {\tilde {H}}(q,p)=H(q,p)+{\tilde {V}}(\xi ),

and the probability density of ξ in the biased ensemble is:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): {\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}}}}

It can be shown that the biased and unbiased averages are related via a simple formula:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): 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}}}}}}.

More generally, an observable Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): \langle A\rangle _{{H}} :

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): \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}}

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

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): \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}}}}}}.

Simulation methods such as umbrella sampling^[1] use a bias potential to enhance sampling of ξ in regions with low P(ξ_i) such as transition regions of chemical reactions. The correct distributions are recovered afterwards using the equation for Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): \langle A\rangle _{{H}} above.

A more detailed description of the method can be found in Ref.^[2]. 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.

Supported types of bias potentials

Presently, the following types of bias potential are supported:

harmonic potential (see curve (a) on the plot below)

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): {\displaystyle V_{bias}(\xi) = \sum_{\mu=1}^{M}\frac{1}{2}\kappa_{\mu} (\xi(q)-\xi_{0,\mu})^2 \; }

where the sum runs all (Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): M ) coordinates the potential acts upon, which are defined in the ICONST-file by setting the status to 8. The parameters of the potential, Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): \kappa (force constant) and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): {\displaystyle \xi_0} (minimum of potential), are defined, respectively, via the keywords SPRING_K and SPRING_R0 in INCAR. Optionally, it is also possible to change the value of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): {\displaystyle \xi_{0,\mu}} every MD step at a constant rate defined via parameter SPRING_V. Note that the number of items defined via SPRING_K, SPRING_R0, and SPRING_V must be equal to Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): M , otherwise the calculation terminates with an error message.

Gauss function (see curve (b) on the plot below)

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): {\displaystyle V_{bias}(\xi) = h\,\text{exp}\left [-\frac{(\xi(q)-\xi_0)^2}{2w^2} \right ], \; }

Fermi function (see curve (b) on the plot below)

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://www.vasp.at/wiki/restbase/vasp.at/v1/":): {\displaystyle V_{bias}(\xi) = \frac{A}{1+\text{exp}\left [-D\frac{\xi(q)}{\xi_0} -1 \right ]} \; }

Andersen thermostat

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=1 (MDALGO=11 in VASP 5.x), 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

Nose-Hoover thermostat

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=2 (MDALGO=21 in VASP 5.x), 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.

References

↑ G. M. Torrie and J. P. Valleau, J. Comp. Phys. 23, 187 (1977).
↑ D. Frenkel and B. Smit, Understanding molecular simulations: from algorithms to applications, Academic Press: San Diego, 2002.

[Torrie77-1] G. M. Torrie and J. P. Valleau, J. Comp. Phys. 23, 187 (1977).

[FrenkelSmit-2] D. Frenkel and B. Smit, Understanding molecular simulations: from algorithms to applications, Academic Press: San Diego, 2002.

[1]

[2]

@@ Line 45: / Line 45: @@
 *harmonic potential (see curve (a) on the plot below)
 :<math>
-V_{bias}(\xi) = \sum_{\mu}^{M}\frac{1}{2}\kappa_{\mu} (\xi(q)-\xi_{0,\mu})^2 \;
+V_{bias}(\xi) = \sum_{\mu=1}^{M}\frac{1}{2}\kappa_{\mu} (\xi(q)-\xi_{0,\mu})^2 \;
 </math>
-where the sum runs all coordinates the potential acts upon, which are defined in the {{FILE|ICONST}}-file by setting the status to 8.
+where the sum runs all (<math>M</math>) coordinates the potential acts upon, which are defined in the {{FILE|ICONST}}-file by setting the status to 8.
-The parameters <math>\kappa</math> (force constant) and <math>\xi_0</math> (minimum of potential)
+The parameters of the potential, <math>\kappa</math> (force constant) and <math>\xi_0</math> (minimum of potential), are defined, respectively, via the keywords {{TAG|SPRING_K}} and {{TAG|SPRING_R0}} in {{FILE|INCAR}}.
-The potential acts on coordinates with status 8 defined in the {{FILE|ICONST}}-file.
+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}}.
+Note that the number of items defined via {{TAG|SPRING_K}}, {{TAG|SPRING_R0}}, and {{TAG|SPRING_V}} must be equal to <math>M</math>, otherwise the calculation terminates with an error message.
 *Gauss function (see curve (b) on the plot below)