MobleyLab / alchemical-analysis

An open tool implementing some recommended practices for analyzing alchemical free energy calculations
MIT License
114 stars 58 forks source link

Unable to reproduce results from Gromacs data sample #112

Closed chtat12 closed 5 years ago

chtat12 commented 5 years ago

I recently obtained the code via git clone https://github.com/MobleyLab/alchemical-analysis.git The file /samples/gromacs/output_11steps/screen_printout.txt contains a block of text below:

 Number of correlated and uncorrelated samples:

 State            N          N_k        N/N_k

     0          501          402         1.25
     1          501          501         1.00
     2          501          414         1.21
     3          501          248         2.02
     4          501          501         1.00
     5          501          501         1.00
     6          501          501         1.00
     7          501          484         1.03
     8          501          501         1.00
     9          501          445         1.13
    10          501          501         1.00

The total change in free energy estimated using MBAR method located at the bottom of this file is 20.017320 +- 0.608113 kJ/mol.

However, I got (slightly) different results when I applied the script /alchemical-analysis/alchemical_analysis/alchemical_analysis.py to the data set located at /alchemical-analysis/samples/gromacs/data/3-methylindole-11steps

Below is the corresponding text from locally generated results:

Number of correlated and uncorrelated samples:

 State            N          N_k        N/N_k

     0          341          341         1.47
     1          407          407         1.23
     2          388          388         1.29
     3          248          248         2.02
     4          501          501         1.00
     5          480          480         1.04
     6          329          329         1.53
     7          183          183         2.74
     8          275          275         1.82
     9          331          331         1.51
    10          501          501         1.00

The total change in free energy estimated using MBAR method is 20.021332 +- 0.727413 kJ/mol.

All the numbers in results.txt also seem to differ from the pre-made output file by a small amount.

Is this behavior expected? It appears to me that all the numbers in column N above should be 501, since each xvg file contains 501 data points. Or is there any change in the way the code determines data correlation?

harlor commented 5 years ago

Or is there any change in the way the code determines data correlation?

Yes, the default uncorrelation method changed in c7f227a7ed7fc43f2deeca235c7972c499d75b3f

To reproduce the previous behavior try: alchemical_analysis.py -n dhdl_all

davidlmobley commented 5 years ago

Thanks, @harlor !