GOMC parser

[msoroush]

Hi,

I implemented GOMC parser for alchemlyb and you can find my implementation in my GitHub.

I have few questions about alchemlyb and I would appreciate your help and suggestion.

1. In alchemical-analysis, if PV value was available, was added to u_nk. However, in alchemlyb, Total Energy is also added to u_nk. Is there any reason you chose different approach?

2. In MBAR documentation it is mentioned that The MBAR method is only applicable to uncorrelated samples from probability distributions. Is there any appropriate series that can be applied to Monte Carlo results for equilibrium_detection, statistical_inefficiency, and more ?

[orbeckst]

Could you please create a pull request with your changes? That will make it easier to discuss it.

1. Can you point to the code that you're looking at.

2. How do you decide for MC data what counts as uncorrelated? The MC Field should have figured that out, or not?

[msoroush]

In line 70 and 72 of gmx.py:

# gromacs also gives us total/potential energy U directly; need this for reduced potential
    u_cols = [col for col in df.columns if any(single_u_col_match in col for single_u_col_match in u_col_match)]
    u = None
    if u_cols:
        u = df[u_cols[0]]

    u_k = dict()
    cols = list()
    for col in dH:
        u_col = eval(col.split('to')[1])
        # calculate reduced potential u_k = dH + pV + U
        u_k[u_col] = beta * dH[col].values
        if pv_cols:
            u_k[u_col] += beta * pv.values
        if u_cols:
            u_k[u_col] += beta * u.values
        cols.append(u_col)

Using alchemical-analysis tools, for some states, I get less than 50 uncorrelated samples (over 25000 samples). I am trying to find out if this is a limitation of MC sampling or not. Low uncorrelated samples is not limited to lambda near 1.0 (full interaction), but also it can happen at lower lambdas (~0.3), where the molecule can easily move around the simulation box.

[msoroush]

Here are the calculated FreeEnergy using TI and MBAR by skipping samples, using slicing function. I have 40,000 samples, each sample taken after 500 Monte Carlo steps.

Working on MBAR method.
Free Energy from TI by skipping    1 sample: -24.2326 (0.0448) kJ/mol 
Free Energy from TI by skipping    2 sample: -24.2066 (0.0632) kJ/mol 
Free Energy from TI by skipping    5 sample: -24.2195 (0.0999) kJ/mol 
Free Energy from TI by skipping   10 sample: -24.2230 (0.1404) kJ/mol 
Free Energy from TI by skipping   20 sample: -24.1591 (0.1990) kJ/mol 
Free Energy from TI by skipping   50 sample: -24.3787 (0.3156) kJ/mol 
Free Energy from TI by skipping  100 sample: -24.4605 (0.4439) kJ/mol 
Free Energy from TI by skipping  200 sample: -25.1494 (0.6071) kJ/mol 
Free Energy from TI by skipping  500 sample: -25.5409 (0.9411) kJ/mol 

Working on MBAR method.
Free Energy from MBAR by skipping    1 sample: -24.2724 (0.0439) kJ/mol 
Free Energy from MBAR by skipping    2 sample: -24.2323 (0.0621) kJ/mol 
Free Energy from MBAR by skipping    5 sample: -24.2198 (0.0982) kJ/mol 
Free Energy from MBAR by skipping   10 sample: -24.1747 (0.1389) kJ/mol 
Free Energy from MBAR by skipping   20 sample: -24.1215 (0.1966) kJ/mol 
Free Energy from MBAR by skipping   50 sample: -24.3783 (0.3097) kJ/mol 
Free Energy from MBAR by skipping  100 sample: -24.3614 (0.4387) kJ/mol 
Free Energy from MBAR by skipping  200 sample: -24.8116 (0.6188) kJ/mol 
Free Energy from MBAR by skipping  500 sample: -25.2521 (0.9738) kJ/mol 

All of the data for a given method agree within the uncertainties. These results suggest that even for the largest number of samples I have uncorrelated data.

[mrshirts]

All of the data for a given method agree within the uncertainties. These results suggest that even for the largest number of samples I have uncorrelated data.

The averages should be consistent no matter what, assuming they are in equilibrium. It's the same time series the whole way though, so it has to be self-consistent. So that does not prove that the samples are uncorrelated. Why not do a time correlation analysis on an underlying observable like the potential energy?

[msoroush]

Why not do a time correlation analysis on an underlying observable like the potential energy?

@mrshirts How noisy is the correlation function for MD? In MBAR timeseries module,g is calculated when the correlation function crossed zero.

In paper by Gowers, they used different approach to calculate g because MC data are noisy. They calculated the normalized autocorrelation function (ACF), then fitted portion of it (AFC [0.1-1.0]) to an exponential decay to calculate tau and then g=(1+2*tau).

Here are the time correlation analysis on my system, using dH/dL value for VDW and Coulomb. In my system, lambda=0.0 correspond to no interaction and lambda=1.0 correspond to full interaction.

Lambda-VDW: [0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00]
dH/dL (VDW) for Lambda index:   0, Eq-index:      0, g:     1.000, Neff_max:  25001
dH/dL (VDW) for Lambda index:   1, Eq-index:      0, g:     1.000, Neff_max:  19618
dH/dL (VDW) for Lambda index:   2, Eq-index:    813, g:     1.028, Neff_max:  23517
dH/dL (VDW) for Lambda index:   3, Eq-index:     69, g:     7.205, Neff_max:   3460
dH/dL (VDW) for Lambda index:   4, Eq-index:    713, g:    41.300, Neff_max:    588
dH/dL (VDW) for Lambda index:   5, Eq-index:    258, g:    74.336, Neff_max:    332
dH/dL (VDW) for Lambda index:   6, Eq-index:   4162, g:    41.059, Neff_max:    507
dH/dL (VDW) for Lambda index:   7, Eq-index:    176, g:    31.137, Neff_max:    797
dH/dL (VDW) for Lambda index:   8, Eq-index:     63, g:    33.476, Neff_max:    744
dH/dL (VDW) for Lambda index:   9, Eq-index:   7067, g:    32.934, Neff_max:    544
dH/dL (VDW) for Lambda index:  10, Eq-index:    286, g:    22.882, Neff_max:   1080
dH/dL (VDW) for Lambda index:  11, Eq-index:  22621, g:     3.976, Neff_max:    598
dH/dL (VDW) for Lambda index:  12, Eq-index:    223, g:    22.200, Neff_max:   1116
dH/dL (VDW) for Lambda index:  13, Eq-index:     30, g:    27.768, Neff_max:    899
dH/dL (VDW) for Lambda index:  14, Eq-index:    847, g:    33.916, Neff_max:    712
dH/dL (VDW) for Lambda index:  15, Eq-index:     68, g:    55.386, Neff_max:    450
dH/dL (VDW) for Lambda index:  16, Eq-index:   2062, g:    25.103, Neff_max:    913
dH/dL (VDW) for Lambda index:  17, Eq-index:  16287, g:    15.269, Neff_max:    570
dH/dL (VDW) for Lambda index:  18, Eq-index:  17634, g:    20.325, Neff_max:    362
dH/dL (VDW) for Lambda index:  19, Eq-index:  22007, g:    14.227, Neff_max:    210
dH/dL (VDW) for Lambda index:  20, Eq-index:      0, g:    81.543, Neff_max:    306
dH/dL (VDW) for Lambda index:  21, Eq-index:  21324, g:    33.874, Neff_max:    108
dH/dL (VDW) for Lambda index:  22, Eq-index:   4753, g:   102.613, Neff_max:    197
dH/dL (VDW) for Lambda index:  23, Eq-index:  23934, g:    11.167, Neff_max:     95
dH/dL (VDW) for Lambda index:  24, Eq-index:   7455, g:   135.804, Neff_max:    129
dH/dL (VDW) for Lambda index:  25, Eq-index:  17690, g:    29.128, Neff_max:    250
dH/dL (VDW) for Lambda index:  26, Eq-index:   1070, g:   147.509, Neff_max:    162
dH/dL (VDW) for Lambda index:  27, Eq-index:  17852, g:    45.759, Neff_max:    156

Lambda-Coul: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.2, 0.4, 0.6, 0.7, 0.8, 0.9, 1.0]
dH/dL (Coul) for Lambda index:   0, Eq-index:      0, g:     1.046, Neff_max:  23898
dH/dL (Coul) for Lambda index:   1, Eq-index:    810, g:     1.037, Neff_max:  18134
dH/dL (Coul) for Lambda index:   2, Eq-index:    143, g:     1.060, Neff_max:  23441
dH/dL (Coul) for Lambda index:   3, Eq-index:      0, g:     1.003, Neff_max:  24916
dH/dL (Coul) for Lambda index:   4, Eq-index:      0, g:     1.008, Neff_max:  24802
dH/dL (Coul) for Lambda index:   5, Eq-index:   2550, g:     1.020, Neff_max:  22016
dH/dL (Coul) for Lambda index:   6, Eq-index:     66, g:     1.099, Neff_max:  22686
dH/dL (Coul) for Lambda index:   7, Eq-index:      1, g:     1.103, Neff_max:  22663
dH/dL (Coul) for Lambda index:   8, Eq-index:     67, g:     1.107, Neff_max:  22525
dH/dL (Coul) for Lambda index:   9, Eq-index:   1108, g:     1.638, Neff_max:  14583
dH/dL (Coul) for Lambda index:  10, Eq-index:   7056, g:     1.052, Neff_max:  17058
dH/dL (Coul) for Lambda index:  11, Eq-index:      3, g:     1.574, Neff_max:  15878
dH/dL (Coul) for Lambda index:  12, Eq-index:     20, g:     1.897, Neff_max:  13167
dH/dL (Coul) for Lambda index:  13, Eq-index:    761, g:     1.552, Neff_max:  15621
dH/dL (Coul) for Lambda index:  14, Eq-index:      7, g:     1.675, Neff_max:  14922
dH/dL (Coul) for Lambda index:  15, Eq-index:   1551, g:     2.436, Neff_max:   9625
dH/dL (Coul) for Lambda index:  16, Eq-index:    562, g:    12.512, Neff_max:   1953
dH/dL (Coul) for Lambda index:  17, Eq-index:   8472, g:     1.339, Neff_max:  12345
dH/dL (Coul) for Lambda index:  18, Eq-index:   9366, g:     1.384, Neff_max:  11297
dH/dL (Coul) for Lambda index:  19, Eq-index:    807, g:     1.397, Neff_max:  17312
dH/dL (Coul) for Lambda index:  20, Eq-index:   4135, g:     1.238, Neff_max:  16850
dH/dL (Coul) for Lambda index:  21, Eq-index:  19826, g:     8.163, Neff_max:    633
dH/dL (Coul) for Lambda index:  22, Eq-index:    665, g:    37.451, Neff_max:    649
dH/dL (Coul) for Lambda index:  23, Eq-index:  14579, g:    25.474, Neff_max:    409
dH/dL (Coul) for Lambda index:  24, Eq-index:  22660, g:     6.064, Neff_max:    386
dH/dL (Coul) for Lambda index:  25, Eq-index:   8148, g:   180.121, Neff_max:     93
dH/dL (Coul) for Lambda index:  26, Eq-index:  24933, g:     1.000, Neff_max:     68
dH/dL (Coul) for Lambda index:  27, Eq-index:  17770, g:   195.276, Neff_max:     37

In VDW case, for lambda index 20-27, lambda value is 1.0 and they start from same configuration, but equilibration steps and statistical inefficiencies are different.

[mrshirts]

@mrshirts https://github.com/mrshirts How noisy is the correlation function for MD? In MBAR timeseries module,g is calculated when the correlation function crossed zero.

It can be a fairly noisy. We are looking into improving it. I think the key is that SOME reasonable correlation analysis is done.

In paper by Gowers https://www.tandfonline.com/doi/full/10.1080/08927022.2017.1375492, they used different approach to calculate g because MC data are noisy. They calculated the normalized autocorrelation function (ACF), then fitted portion of it (AFC [0.1-1.0]) to an exponential decay to calculate tau and then g=(1+2*tau).

I think this is very reasonable to do if there is sufficient data to fit, which usually is the case for liquids. We are still investigating best practices.

Here are the time correlation analysis on my system, using dH/dL value for VDW and Coulomb. In my system, lambda=0.0 correspond to no interaction and lambda=1.0 correspond to full interaction.

Lambda-VDW: [0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00] dH/dL (VDW) for Lambda index: 0, Eq-index: 0, g: 1.000, Neff_max: 25001 dH/dL (VDW) for Lambda index: 1, Eq-index: 0, g: 1.000, Neff_max: 19618 dH/dL (VDW) for Lambda index: 2, Eq-index: 813, g: 1.028, Neff_max: 23517 dH/dL (VDW) for Lambda index: 3, Eq-index: 69, g: 7.205, Neff_max: 3460 dH/dL (VDW) for Lambda index: 4, Eq-index: 713, g: 41.300, Neff_max: 588 dH/dL (VDW) for Lambda index: 5, Eq-index: 258, g: 74.336, Neff_max: 332 dH/dL (VDW) for Lambda index: 6, Eq-index: 4162, g: 41.059, Neff_max: 507 dH/dL (VDW) for Lambda index: 7, Eq-index: 176, g: 31.137, Neff_max: 797 dH/dL (VDW) for Lambda index: 8, Eq-index: 63, g: 33.476, Neff_max: 744 dH/dL (VDW) for Lambda index: 9, Eq-index: 7067, g: 32.934, Neff_max: 544 dH/dL (VDW) for Lambda index: 10, Eq-index: 286, g: 22.882, Neff_max: 1080 dH/dL (VDW) for Lambda index: 11, Eq-index: 22621, g: 3.976, Neff_max: 598 dH/dL (VDW) for Lambda index: 12, Eq-index: 223, g: 22.200, Neff_max: 1116 dH/dL (VDW) for Lambda index: 13, Eq-index: 30, g: 27.768, Neff_max: 899 dH/dL (VDW) for Lambda index: 14, Eq-index: 847, g: 33.916, Neff_max: 712 dH/dL (VDW) for Lambda index: 15, Eq-index: 68, g: 55.386, Neff_max: 450 dH/dL (VDW) for Lambda index: 16, Eq-index: 2062, g: 25.103, Neff_max: 913 dH/dL (VDW) for Lambda index: 17, Eq-index: 16287, g: 15.269, Neff_max: 570 dH/dL (VDW) for Lambda index: 18, Eq-index: 17634, g: 20.325, Neff_max: 362 dH/dL (VDW) for Lambda index: 19, Eq-index: 22007, g: 14.227, Neff_max: 210 dH/dL (VDW) for Lambda index: 20, Eq-index: 0, g: 81.543, Neff_max: 306 dH/dL (VDW) for Lambda index: 21, Eq-index: 21324, g: 33.874, Neff_max: 108 dH/dL (VDW) for Lambda index: 22, Eq-index: 4753, g: 102.613, Neff_max: 197 dH/dL (VDW) for Lambda index: 23, Eq-index: 23934, g: 11.167, Neff_max: 95 dH/dL (VDW) for Lambda index: 24, Eq-index: 7455, g: 135.804, Neff_max: 129 dH/dL (VDW) for Lambda index: 25, Eq-index: 17690, g: 29.128, Neff_max: 250 dH/dL (VDW) for Lambda index: 26, Eq-index: 1070, g: 147.509, Neff_max: 162 dH/dL (VDW) for Lambda index: 27, Eq-index: 17852, g: 45.759, Neff_max: 156

Lambda-Coul: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.2, 0.4, 0.6, 0.7, 0.8, 0.9, 1.0] dH/dL (Coul) for Lambda index: 0, Eq-index: 0, g: 1.046, Neff_max: 23898 dH/dL (Coul) for Lambda index: 1, Eq-index: 810, g: 1.037, Neff_max: 18134 dH/dL (Coul) for Lambda index: 2, Eq-index: 143, g: 1.060, Neff_max: 23441 dH/dL (Coul) for Lambda index: 3, Eq-index: 0, g: 1.003, Neff_max: 24916 dH/dL (Coul) for Lambda index: 4, Eq-index: 0, g: 1.008, Neff_max: 24802 dH/dL (Coul) for Lambda index: 5, Eq-index: 2550, g: 1.020, Neff_max: 22016 dH/dL (Coul) for Lambda index: 6, Eq-index: 66, g: 1.099, Neff_max: 22686 dH/dL (Coul) for Lambda index: 7, Eq-index: 1, g: 1.103, Neff_max: 22663 dH/dL (Coul) for Lambda index: 8, Eq-index: 67, g: 1.107, Neff_max: 22525 dH/dL (Coul) for Lambda index: 9, Eq-index: 1108, g: 1.638, Neff_max: 14583 dH/dL (Coul) for Lambda index: 10, Eq-index: 7056, g: 1.052, Neff_max: 17058 dH/dL (Coul) for Lambda index: 11, Eq-index: 3, g: 1.574, Neff_max: 15878 dH/dL (Coul) for Lambda index: 12, Eq-index: 20, g: 1.897, Neff_max: 13167 dH/dL (Coul) for Lambda index: 13, Eq-index: 761, g: 1.552, Neff_max: 15621 dH/dL (Coul) for Lambda index: 14, Eq-index: 7, g: 1.675, Neff_max: 14922 dH/dL (Coul) for Lambda index: 15, Eq-index: 1551, g: 2.436, Neff_max: 9625 dH/dL (Coul) for Lambda index: 16, Eq-index: 562, g: 12.512, Neff_max: 1953 dH/dL (Coul) for Lambda index: 17, Eq-index: 8472, g: 1.339, Neff_max: 12345 dH/dL (Coul) for Lambda index: 18, Eq-index: 9366, g: 1.384, Neff_max: 11297 dH/dL (Coul) for Lambda index: 19, Eq-index: 807, g: 1.397, Neff_max: 17312 dH/dL (Coul) for Lambda index: 20, Eq-index: 4135, g: 1.238, Neff_max: 16850 dH/dL (Coul) for Lambda index: 21, Eq-index: 19826, g: 8.163, Neff_max: 633 dH/dL (Coul) for Lambda index: 22, Eq-index: 665, g: 37.451, Neff_max: 649 dH/dL (Coul) for Lambda index: 23, Eq-index: 14579, g: 25.474, Neff_max: 409 dH/dL (Coul) for Lambda index: 24, Eq-index: 22660, g: 6.064, Neff_max: 386 dH/dL (Coul) for Lambda index: 25, Eq-index: 8148, g: 180.121, Neff_max: 93 dH/dL (Coul) for Lambda index: 26, Eq-index: 24933, g: 1.000, Neff_max: 68 dH/dL (Coul) for Lambda index: 27, Eq-index: 17770, g: 195.276, Neff_max: 37

In VDW case, for lambda index 20-27, lambda value is 1.0 and they start from same configuration, but equilibration steps and statistical inefficiencies are different.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/alchemistry/alchemlyb/issues/77#issuecomment-488777933, or mute the thread https://github.com/notifications/unsubscribe-auth/ABATPVHZSUZL3NEA5W66FQ3PTMWENANCNFSM4HJGU5TQ .

[msoroush]

@mrshirts I calculated the correlation function (C) after equilibration step, using MBAR function and plot it versus number of samples. The last data point correspond to uncorrelated samples. Here are the figures for selected lambdas.

The data bellow C=0.1 are very noisy. I dont think I can use the MBAR time correlation analysis for my data.

[mrshirts]

I agree that things could be improved. I think that including an exponential fit option and an "area under the curve" option could be useful, and understanding the differences between the methods would be good. I would like to try to put them in this summer, though I doubt it will happen in the next few weeks.

[mrshirts]

For the near term, I would suggest decorrelating your data youself beforehand (using your preferred method), and then using alchemlyb to do the rest. Clearly, the data do have correlation times of 20-100 steps, and you definitely should subsample at about that level.

[msoroush]

@mrshirts Thank you for your help. I have few more questions, regarding the Free Energy calculation. Is that okay if I email you?

[orbeckst]

@msoroush thanks for opening PR #78.

U was added in PR #62 – see there for discussion.

[mrshirts]

Having it here on alchemlyb might make it more useful to the community, but might belong in another issue, depending on the question.

[msoroush]

My question is about calculation of dH/dL or deltaH.

1- In MC, we are using hard cutoff to avoid atoms of the molecule get too close (when we have polar molecule, atoms with opposite charge stick together at close distance). For instance, to simulate water, we use usually use hard cutoff = 1 A. Do you think using hard cutoff would change the free energy value? For lambda > 0.5, naturally, atoms dont get closer than this distance. If it does change the free energy value, what is the safest hard cutoff?

2- Using soft core distance for coulomb interaction, would not cause error in long range electrostatic energy (Ewald summation method)?

3- Is there any contribution from intra-nonbonded energy of the solute molecule in free energy calculation? If yes, is it coming from total energy value or not?

Sorry about taking your time and I really appreciate your help.

[mrshirts]

2- Using soft core distance for coulomb interaction, would not cause error in long range electrostatic energy (Ewald summation method)?

Usually one would use soft core for the short range coulomb, and linear coupling for the long range Ewald summation.

[mrshirts]

1- In MC, we are using hard cutoff to avoid atoms of the molecule get too close (when we have polar molecule, atoms with opposite charge stick together at close distance). For instance, to simulate water, we use usually use hard cutoff = 1 A. Do you think using hard cutoff would change the free energy value? For lambda > 0.5, naturally, atoms dont get closer than this distance. If it does change the free energy value, what is the safest hard cutoff? 3- Is there any contribution from intra-nonbonded energy of the solute molecule in free energy calculation? If yes, is it coming from total energy value or not?

Any free energy calculation, if implemented correctly, will give the free energy from changing the Hamiltonian from the initial Hamiltonian to the final Hamiltonian. What you describe above are questions of whether certain choices of the initial and final Hamiltonian are consistent with each other - that's not really a free energy question, but a molecular model question. Whether one needs an intramolecular nonbonded energy is a more general question.

[msoroush]

Whether one needs an intramolecular nonbonded energy is a more general question

To be more specific, if someone is calculating the free energy of n-octane (with no electrostatic interaction) in vacuum using NVT ensemble, what is the expected value? I only calculate the derivative of energy (nonbonded and long range corrections) with respect to lambda. So, in this example, the only contribution must come from long range correction of VDW. In paper by Garrido they mentioned:

The separate calculation in vacuum is necessary to compensate for changes in solute-solute intramolecular nonbonded interactions that take place when the intermolecular interactions are switched off

and long range correction of VDW is independent from molecular conformation. I just want to make sure I include all necessary energies.

[mrshirts]

To be more specific, if someone is calculating the free energy of n-octane (with no electrostatic interaction) in vacuum using NVT ensemble, what is the expected value?

I really don't know. It depends on the choice of parameters.

Neglecting the intermolecular van der Waals is an approximation, and I really don't know how big it is. But that is question that should be independent of the analysis. If your energy differences that you stick into MBAR are not actually the differences between the energies that are used in the Boltzmann weighting of the configurations (which could happen by neglecting intramolecular interactions), the free energies coming out of MBAR won't be correct.

[msoroush]

Neglecting the intermolecular van der Waals is an approximation, and I really don't know how big it is. But that is question that should be independent of the analysis.

Just to be clear, I am not neglecting the intra nonbonded energy in my simulation. I only couple the nonbonded energy with lambda. When I evaluate the Boltzmann weight (-Beta * change in energy acting on the same configuration space), bonded and intra nonbonded energy will cancel out, since the solute's configuration is same.

This is a snapshot of Bennett paper:

If your energy differences that you stick into MBAR are not actually the differences between the energies that are used in the Boltzmann weighting of the configurations (which could happen by neglecting intramolecular interactions), the free energies coming out of MBAR won't be correct.

I feed (deltaE + PV + totalEnergy) to MBAR, where deltaE is change in energy (nonbonded and long range correction), PV (if we are doing NPT simulation), and total energy (as discussed in issue #62)