Blue-moon ensemble: Difference between revisions
No edit summary |
m (Karsai moved page Blue moon ensemble to Blue-moon ensemble) |
||
(6 intermediate revisions by 3 users not shown) | |||
Line 1: | Line 1: | ||
In general, constrained molecular dynamics generates biased statistical averages. | In general, constrained molecular dynamics generates biased statistical averages. It can be shown that the correct average for a quantity <math>a(\xi)</math> can be obtained using the formula (blue moon ensemble average): | ||
It can be shown that the correct average for a quantity <math>a(\xi)</math> can be obtained using the formula: | |||
:<math> | :<math> | ||
a(\xi)=\frac{\langle |\mathbf{Z}|^{-1/2} a(\xi^*) \rangle_{\xi^*}}{\langle |\mathbf{Z}|^{-1/2}\rangle_{\xi^*}}, | a(\xi)=\frac{\langle |\mathbf{Z}|^{-1/2} a(\xi^*) \rangle_{\xi^*}}{\langle |\mathbf{Z}|^{-1/2}\rangle_{\xi^*}}, | ||
Line 8: | Line 7: | ||
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, | 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, | ||
</math> | </math> | ||
It can be shown that the free energy gradient can be computed using the equation:<ref name="Carter89"/><ref name="Otter00"/><ref name="Darve02"/><ref name="Fleurat05"/> | It can be shown that the free energy gradient can be computed using the equation:<ref name="Carter89" /><ref name="Otter00" /><ref name="Darve02" /><ref name="Fleurat05" /> | ||
:<math> | :<math> | ||
\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^*}, | \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^*}, | ||
</math> | </math> | ||
where <math>\lambda_{\xi_k}</math> is the Lagrange multiplier associated with the parameter <math>{\xi_k}</math> used in the [[ | where <math>\lambda_{\xi_k}</math> is the Lagrange multiplier associated with the parameter <math>{\xi_k}</math> used in the [[Constrained molecular dynamics|SHAKE algorithm]].<ref name="Ryckaert77" /> | ||
The free-energy difference between states (1) and (2) can be computed by integrating the free-energy gradients over a connecting path: | The free-energy difference between states (1) and (2) can be computed by integrating the free-energy gradients over a connecting path: | ||
Line 20: | Line 19: | ||
Note that as the free-energy is a state quantity, the choice of path connecting (1) with (2) is irrelevant. | Note that as the free-energy is a state quantity, the choice of path connecting (1) with (2) is irrelevant. | ||
== 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, free energy | 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 the string '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. In the simplest case when only one 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 <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 | |||
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== | ==References== | ||
Line 43: | Line 50: | ||
[[Category:Molecular dynamics]] | [[Category:Molecular dynamics]] | ||
[[Category:Theory]] | [[Category:Theory]] | ||
[[Category:Howto]] | [[Category:Howto]] |
Latest revision as of 08:18, 19 July 2022
In general, constrained molecular dynamics generates biased statistical averages. It can be shown that the correct average for a quantity can be obtained using the formula (blue moon ensemble average):
where stands for the statistical average of the quantity enclosed in angular parentheses computed for a constrained ensemble and is a mass metric tensor defined as:
It can be shown that the free energy gradient can be computed using the equation:[1][2][3][4]
where is the Lagrange multiplier associated with the parameter 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:
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 , , , and , respectively. Note that one line introduced by the string 'b_m>' is written for each constrained coordinate. With this output, the free energy gradient with respected to the fixed coordinate can conveniently be determined (by equation given above) as a ratio between averages of the last and the second numerical terms. In the simplest case when only one 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 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).