search

choderalab / alchemy

Alchemical tools for OpenMM
9 stars 3 forks source link

Add more testing to catch subtle bugs #38

Closed jchodera closed 8 years ago

andrrizzi commented 8 years ago

I've noticed that in several points in test_alchemy.py we use lambda_coulomb instead of lambda_electrostatics. When I switched the two in test_annihilated_states(), I got differences that can explain what happens with solvation free energies! I'll have some time to verify this on Tuesday/Wednesday.

In general, I think it is very easy to introduce this kind of bugs when passing arguments as kwargs without any type of checking (as we do in AlchemicalState). I usually prefer sticking to normal named_argument=default unless absolutely necessary.

After correcting the lambda_coulomb to lambda_electrostatics thing, this is ready to be merged for me!

jchodera commented 8 years ago

Great catch! The parameter had previously been named lambda_coulomb, but I migrated later to be more consistent and must have missed this.

In general, I think it is very easy to introduce this kind of bugs when passing arguments as kwargs without any type of checking (as we do in AlchemicalState). I usually prefer sticking to normal named_argument=default unless absolutely necessary.

In this case, we definitely want these arguments to be optional, since we want to be able to (1) only specify arguments that are not pinned to 1.0, and (2) add arguments without breaking existing code. But I should have inserted a check that made sure that no unknown arguments were provided to AlchemicalState.

After correcting the lambda_coulomb to lambda_electrostatics thing, this is ready to be merged for me!

Do you want to tackle that, or should I?

jchodera commented 8 years ago

Do you want to tackle that, or should I?

Oh, nevermind! This is my PR, so I'll fix it!

Thanks again!

jchodera commented 8 years ago

Hm, I am still only seeing the water box tests failing.

If you have a chance, could you give a bit more info about the new failures you saw? Did you make additional changes beyond changing lamda_coulomb? Or did you test on CUDA/OpenCL?

jchodera commented 8 years ago

Oh! Nevermind---I didn't actually add the check that throws an exception in test_annihilated_states. Thanks!

jchodera commented 8 years ago

OK, I can confirm I see this discrepancy. This is with CutoffPeriodic:

2016-08-13 09:42:49,785: Creating alchemically modified system...
2016-08-13 09:42:49,797: Elapsed time 0.012 s.
2016-08-13 09:42:49,797: Creating alchemically modified intermediate...
2016-08-13 09:42:49,800: Elapsed time 0.003 s.
2016-08-13 09:42:49,803: Creating alchemically modified system...
2016-08-13 09:42:49,812: Elapsed time 0.009 s.
2016-08-13 09:42:49,812: Creating alchemically modified intermediate...
2016-08-13 09:42:49,815: Elapsed time 0.003 s.
2016-08-13 09:42:50,018: vacuum   lambda = 1 :    2.999 kcal/mol
2016-08-13 09:42:50,018: vacuum   lambda = 0 :    0.168 kcal/mol
2016-08-13 09:42:50,018: difference          :    2.831 kcal/mol
2016-08-13 09:42:50,018: periodic lambda = 1 :    1.337 kcal/mol
2016-08-13 09:42:50,018: periodic lambda = 0 :    0.168 kcal/mol
2016-08-13 09:42:50,018: difference          :    1.169 kcal/mol

The lambda = 0 state has the same energy between NoCutoff and CutoffPeriodic, but the lambda = 1 state does not.

For PME, we recover the same energies at lambda = 0 and lambda = 1, however:

2016-08-13 09:46:29,330: Creating alchemically modified system...
2016-08-13 09:46:29,343: Elapsed time 0.013 s.
2016-08-13 09:46:29,344: Creating alchemically modified intermediate...
2016-08-13 09:46:29,347: Elapsed time 0.003 s.
2016-08-13 09:46:29,349: Creating alchemically modified system...
2016-08-13 09:46:29,358: Elapsed time 0.008 s.
2016-08-13 09:46:29,358: Creating alchemically modified intermediate...
2016-08-13 09:46:29,363: Elapsed time 0.005 s.
2016-08-13 09:46:29,734: vacuum   lambda = 1 :    2.999 kcal/mol
2016-08-13 09:46:29,734: vacuum   lambda = 0 :    0.168 kcal/mol
2016-08-13 09:46:29,735: difference          :    2.831 kcal/mol
2016-08-13 09:46:29,735: periodic lambda = 1 :    2.999 kcal/mol
2016-08-13 09:46:29,735: periodic lambda = 0 :    0.168 kcal/mol
2016-08-13 09:46:29,735: difference          :    2.831 kcal/mol
jchodera commented 8 years ago

It looks like the energy discrepancy in the fully-interacting state is due to a shift that is applied for reaction-field electrostatics. Reaction-field uses this equation, which adds a term like

E_shift = - ONE_4PI_EPS0 * q_i * q_j * c_rf

where for epsilon_solvent = 1, c_rf just goes to r_cutoff^(-3), meaning we have a shift of

E_shift = - ONE_4PI_EPS0 * q_i * q_j / r_cutoff^3

Now, the problem is that we scale the charges down to 0, meaning this shift goes away, causing a discrepancy in the energy difference from lambda = 1 -> 0 between vacuum and reaction field. The question is whether this discrepancy should remain or not.

In principle, reaction field models the effects of solvent outside the cutoff, so unless we take r_cutoff -> infinity (in which case the shift vanishes), we will see this effect, which may be something we should really be including.

So I think this is actually a correct effect to include. Making the cutoff larger or eliminating c_rf eliminates the discrepancy.

jchodera commented 8 years ago

The OpenMM bug was fixed in https://github.com/pandegroup/openmm/pull/1561