Blue-moon ensemble: Difference between revisions
No edit summary |
m (Karsai moved page Blue moon ensemble to Blue-moon ensemble) |
||
(2 intermediate revisions by 2 users not shown) | |||
Line 11: | Line 11: | ||
\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 50: | 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).