Blue-moon ensemble: Difference between revisions
No edit summary |
No edit summary |
||
| Line 20: | Line 20: | ||
==How to == | ==How to == | ||
The information needed to determine the blue moon ensemble averages within a [[:Category:Constrained molecular dynamics|Constrained molecular dynamics]] can be obtained by setting {{TAG|LBLUEOUT}}=.TRUE. The following output is written for each MD step in the file {{FILE|REPORT}}: | The information needed to determine the blue moon ensemble averages within a [[:Category:Constrained molecular dynamics|Constrained molecular dynamics]] can be obtained by setting {{TAG|LBLUEOUT}}=.TRUE. The following output is written for each MD step in the file {{FILE|REPORT}}: | ||
>Blue_moon | >Blue_moon | ||
lambda |z|^(-1/2) GkT |z|^(-1/2)*(lambda+GkT) | lambda |z|^(-1/2) GkT |z|^(-1/2)*(lambda+GkT) | ||
b_m> 0.585916E+01 0.215200E+02 -0.117679E+00 0.123556E+03 | b_m> 0.585916E+01 0.215200E+02 -0.117679E+00 0.123556E+03 | ||
with the four numerical terms indicating <math>\lambda_{\xi_k}</math>, <math>|Z|^{-1/2}</math>, <math>\left ( \frac{k_B T}{2 |Z|} \sum_{j=1}^{r}(Z^{-1})_{kj} \sum_{i=1}^{3N} m_i^{-1}\nabla_i \xi_j \cdot \nabla_i |Z| \right ) </math>, and <math>\left ( |Z|^{-1/2} [\lambda_k +\frac{k_B T}{2 |Z|} \sum_{j=1}^{r}(Z^{-1})_{kj} \sum_{i=1}^{3N} m_i^{-1}\nabla_i \xi_j \cdot \nabla_i |Z|] \right ) </math>, respectively. With this output, the free energy gradient can conveniently be determined (by equation given above) as a ratio between averages of the last and the second numerical terms, | with the four numerical terms indicating <math>\lambda_{\xi_k}</math>, <math>|Z|^{-1/2}</math>, <math>\left ( \frac{k_B T}{2 |Z|} \sum_{j=1}^{r}(Z^{-1})_{kj} \sum_{i=1}^{3N} m_i^{-1}\nabla_i \xi_j \cdot \nabla_i |Z| \right ) </math>, and <math>\left ( |Z|^{-1/2} [\lambda_k +\frac{k_B T}{2 |Z|} \sum_{j=1}^{r}(Z^{-1})_{kj} \sum_{i=1}^{3N} m_i^{-1}\nabla_i \xi_j \cdot \nabla_i |Z|] \right ) </math>, respectively. Note that one line introduced by b_m is written for each constrained coordinate. With this output, the free energy gradient with respected to the fixed coordinate <math>{\xi_k}</math> can conveniently be determined (by equation given above) as a ratio between averages of the last and the second numerical terms. For instance, if only a single constraint is used, the free energy gradient can be obtained as follows: | ||
grep b_m REPORT |awk 'BEGIN {a=0.;b=0.} {a+=$5;b+=$3} END {print a/b}' | grep b_m REPORT |awk 'BEGIN {a=0.;b=0.} {a+=$5;b+=$3} END {print a/b}' | ||
As example of a blue moon ensemble average, let us consider calculation of unbiased potential energy average from constrained MD. Here we extract <math>|Z|^{-1/2}</math> for each step and store the data in an auxiliary file zet.dat: | As example of a blue moon ensemble average, let us consider calculation of unbiased potential energy average from constrained MD. For simplicity, only a single constraint is assumed. Here we extract <math>|Z|^{-1/2}</math> for each step and store the data in an auxiliary file zet.dat: | ||
grep b_m REPORT |awk '{print $3}' > zet.dat | grep b_m REPORT |awk '{print $3}' > zet.dat | ||
Here we extract potential energy for each step and store the data in an auxiliary file energy.dat: | Here we extract potential energy for each step and store the data in an auxiliary file energy.dat: | ||
Revision as of 16:13, 18 April 2022
In general, constrained molecular dynamics generates biased statistical averages. It can be shown that the correct average for a quantity 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/":): a(\xi ) can be obtained using the formula (blue moon ensemble average):
- 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/":): a(\xi )={\frac {\langle |{\mathbf {Z}}|^{{-1/2}}a(\xi ^{*})\rangle _{{\xi ^{*}}}}{\langle |{\mathbf {Z}}|^{{-1/2}}\rangle _{{\xi ^{*}}}}},
where 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 ...\rangle _{{\xi ^{*}}} stands for the statistical average of the quantity enclosed in angular parentheses computed for a constrained ensemble 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/":): Z is a mass metric tensor defined 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/":): Z_{{\alpha ,\beta }}={\sum }_{{i=1}}^{{3N}}m_{i}^{{-1}}\nabla _{i}\xi _{\alpha }\cdot \nabla _{i}\xi _{\beta },\,\alpha =1,...,r,\,\beta =1,...,r,
It can be shown that the free energy gradient can be computed using the equation:[1][2][3][4]
- 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/":): {\Bigl (}{\frac {\partial A}{\partial \xi _{k}}}{\Bigr )}_{{\xi ^{*}}}={\frac {1}{\langle |Z|^{{-1/2}}\rangle _{{\xi ^{*}}}}}\langle |Z|^{{-1/2}}[\lambda _{k}+{\frac {k_{B}T}{2|Z|}}\sum _{{j=1}}^{{r}}(Z^{{-1}})_{{kj}}\sum _{{i=1}}^{{3N}}m_{i}^{{-1}}\nabla _{i}\xi _{j}\cdot \nabla _{i}|Z|]\rangle _{{\xi ^{*}}},
where 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/":): \lambda _{{\xi _{k}}} is the Lagrange multiplier associated with the parameter 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/":): {\xi _{k}} used in the SHAKE algorithm.[5]
The free-energy difference between states (1) and (2) can be computed by integrating the free-energy gradients over a connecting path:
- 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/":): {\Delta }A_{{1\rightarrow 2}}=\int _{{{\xi (1)}}}^{{{\xi (2)}}}{\Bigl (}{\frac {\partial {A}}{\partial \xi }}{\Bigr )}_{{\xi ^{*}}}\cdot d{\xi }.
Note that as the free-energy is a state quantity, the choice of path connecting (1) with (2) is irrelevant.
How to
The information needed to determine the blue moon ensemble averages within a Constrained molecular dynamics can be obtained by setting LBLUEOUT=.TRUE. The following output is written for each MD step in the file REPORT:
>Blue_moon
lambda |z|^(-1/2) GkT |z|^(-1/2)*(lambda+GkT)
b_m> 0.585916E+01 0.215200E+02 -0.117679E+00 0.123556E+03
with the four numerical terms indicating 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/":): \lambda _{{\xi _{k}}}
, 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 |Z|^{-1/2}}
, 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 \left ( \frac{k_B T}{2 |Z|} \sum_{j=1}^{r}(Z^{-1})_{kj} \sum_{i=1}^{3N} m_i^{-1}\nabla_i \xi_j \cdot \nabla_i |Z| \right ) }
, 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 \left ( |Z|^{-1/2} [\lambda_k +\frac{k_B T}{2 |Z|} \sum_{j=1}^{r}(Z^{-1})_{kj} \sum_{i=1}^{3N} m_i^{-1}\nabla_i \xi_j \cdot \nabla_i |Z|] \right ) }
, respectively. Note that one line introduced by b_m is written for each constrained coordinate. With this output, the free energy gradient with respected to the fixed coordinate 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/":): {\xi _{k}}
can conveniently be determined (by equation given above) as a ratio between averages of the last and the second numerical terms. For instance, if only a single constraint is used, the free energy gradient can be obtained as follows:
grep b_m REPORT |awk 'BEGIN {a=0.;b=0.} {a+=$5;b+=$3} END {print a/b}'
As example of a blue moon ensemble average, let us consider calculation of unbiased potential energy average from constrained MD. For simplicity, only a single constraint is assumed. Here we extract 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 |Z|^{-1/2}} for each step and store the data in an auxiliary file zet.dat:
grep b_m REPORT |awk '{print $3}' > zet.dat
Here we extract potential energy for each step and store the data in an auxiliary file energy.dat:
grep e_b REPORT |awk '{print $3}' > energy.dat
Finally, the weighted average is determined according to the first formula shown above:
paste energy.dat zet.dat |awk 'BEGIN {a=0.;b=0.} {a+=$1*$2;b+=$2} END {print a/b}'
References
- ↑ E. A. Carter, G. Ciccotti, J. T. Hynes, and R. Kapral, Chem. Phys. Lett. 156, 472 (1989).
- ↑ W. K. Den Otter and W. J. Briels, Mol. Phys. 98, 773 (2000).
- ↑ E. Darve, M. A. Wilson, and A. Pohorille, Mol. Simul. 28, 113 (2002).
- ↑ P. Fleurat-Lessard and T. Ziegler, J. Chem. Phys. 123, 084101 (2005).
- ↑ J. P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen, J. Comp. Phys. 23, 327 (1977).