Biased molecular dynamics: Difference between revisions

Revision as of 09:37, 6 April 2023

Simulation methods such as umbrella sampling^[1] use a bias potentials ( $\tilde{V}(\xi)$ ) acting on one or more selected internal coordinates of the system ξ=ξ(q) to enhance sampling of parts of configuration space with low Boltzmann weights. Furthermore, 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.

Presently, the following types of $\tilde{V}(\xi)$ are supported:

Gauss function

{\displaystyle \tilde{V}(\xi)= h\,\text{exp}\left [-\frac{(\xi(q)-\xi_0)^2}{2w^2} \right ], \; }

harmonic potential

{\displaystyle \tilde{V}(\xi) = \frac{1}{2}\kappa (\xi(q)-\xi_0)^2 \; }

Fermi function

{\displaystyle \tilde{V}(\xi)= \frac{A}{1+\text{exp}\left [-D\frac{\xi(q)}{\xi_0} -1 \right ]} \; }

The Hamiltonian of the physical system:

{\displaystyle H(q,p) = T(p) + V(q), \; }

with T(p), and V(q) being kinetic, and potential energies, respectively, can be extended by adding a bias potential $\tilde{V}(\xi)$ acting only on one or more selected internal coordinates of the system ξ=ξ(q):

{\displaystyle \tilde{H}(q,p) = H(q,p) + \tilde{V}(\xi). }

Presently, the following types of $\tilde{V}(\xi)$ are supported:

Gauss function

{\displaystyle \tilde{V}(\xi)= h\,\text{exp}\left [-\frac{(\xi(q)-\xi_0)^2}{2w^2} \right ], \; }

harmonic potential

{\displaystyle \tilde{V}(\xi) = \frac{1}{2}\kappa (\xi(q)-\xi_0)^2 \; }

Fermi function

{\displaystyle \tilde{V}(\xi)= \frac{A}{1+\text{exp}\left [-D\frac{\xi(q)}{\xi_0} -1 \right ]} \; }

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

{\displaystyle H(q,p) = T(p) + V(q), \; }

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

{\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}. }

The term $\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 $\tilde{V}(\xi)$ acting only on a selected internal parameter of the system ξ=ξ(q), the Hamiltonian takes a form:

{\displaystyle \tilde{H}(q,p) = H(q,p) + \tilde{V}(\xi), }

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

{\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}} }

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

{\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}}}. }

More generally, an observable $\langle A \rangle_{H}$ :

{\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} }

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

{\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}}}. }

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 $\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:

Gauss function

{\displaystyle V_{bias}(\xi) = h\,\text{exp}\left [-\frac{(\xi(q)-\xi_0)^2}{2w^2} \right ], \; }

harmonic potential

{\displaystyle V_{bias}(\xi) = \frac{1}{2}\kappa (\xi(q)-\xi_0)^2 \; }

Fermi function

{\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

↑ ^a ^b 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]