Biased molecular dynamics
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 only on a selected internal parameter of the system ξ=ξ(q), the Hamiltonian takes a 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 a simple formula:
- [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} }[/math]:
- [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:
- [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]
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 [math]\displaystyle{ \langle A \rangle_{H} }[/math] 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.
- 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=11, 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
The values of all collective variables for each MD step are listed in the REPORT-file, check the lines after the string Metadynamics.