Blue-moon ensemble: Difference between revisions

Revision as of 16:06, 18 April 2022

In general, constrained molecular dynamics generates biased statistical averages. It can be shown that the correct average for a quantity $a(\xi )$ can be obtained using the formula (blue moon ensemble average):

a(\xi )={\frac {\langle |\mathbf {Z} |^{-1/2}a(\xi ^{*})\rangle _{\xi ^{*}}}{\langle |\mathbf {Z} |^{-1/2}\rangle _{\xi ^{*}}}},

where $\langle ...\rangle _{\xi ^{*}}$ stands for the statistical average of the quantity enclosed in angular parentheses computed for a constrained ensemble and $Z$ is a mass metric tensor defined as:

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]

{\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 $\lambda _{\xi _{k}}$ is the Lagrange multiplier associated with the parameter ${\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:

{\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 $\lambda _{\xi _{k}}$ , $|Z|^{-1/2}$ , $\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 $\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. 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, which can be performed 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. Here we extract $|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

[Carter89-1] E. A. Carter, G. Ciccotti, J. T. Hynes, and R. Kapral, Chem. Phys. Lett. 156, 472 (1989).

[Otter00-2] W. K. Den Otter and W. J. Briels, Mol. Phys. 98, 773 (2000).

[Darve02-3] E. Darve, M. A. Wilson, and A. Pohorille, Mol. Simul. 28, 113 (2002).

[Fleurat05-4] P. Fleurat-Lessard and T. Ziegler, J. Chem. Phys. 123, 084101 (2005).

[Ryckaert77-5] J. P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen, J. Comp. Phys. 23, 327 (1977).

[1]

[2]

[3]

[4]

[5]

@@ 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>
 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,
 </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>
 \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>
-where <math>\lambda_{\xi_k}</math> is the Lagrange multiplier associated with the parameter <math>{\xi_k}</math> used in the [[#SHAKE|SHAKE algorithm]].<ref name="Ryckaert77"/>
+where <math>\lambda_{\xi_k}</math> is the Lagrange multiplier associated with the parameter <math>{\xi_k}</math> used in the [[#SHAKE|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:
@@ Line 20: / Line 19: @@
 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}}:
   >Blue_moon
@@ Line 28: / Line 27: @@
-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 gradients can conveniently be determined by equation given above as ratio between averages of the last and second numerical terms, which can be performed as follows:
+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, which can be performed 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. 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==