Project

General

Profile

Feature #1237

Turning on free energy with GB in 4.6.2 beta yields incorrect (zero) short range Coulomb energies

Added by David Mobley over 4 years ago. Updated about 3 years ago.

Status:
Accepted
Priority:
Low
Assignee:
-
Category:
mdrun
Target version:
-
Difficulty:
uncategorized
Close

Description

Recently, I noticed that in 4.6 (and 4.6.2 beta), free_energy = yes and GB solvent are not incompatible, in that simulations at least run. However, computed free energies are wildly incorrect. On looking into what's going on, it seems that when the free energy code is turned on with GB, the short range coulomb interactions get turned off. Here's an average for my system with free energy on:

Statistics over 114791 steps using 1149 frames
Energies (kJ/mol)
Bond Angle Proper Dih. Improper Dih.GB Polarization
2.83949e+01 6.88145e+01 1.10921e+02 5.45136e+00 -7.09138e+02
Nonpolar Sol. LJ-14 Coulomb-14 LJ (SR) LJ (LR)
2.35202e+01 2.36118e+01 1.27878e+02 -1.91885e+01 0.00000e+00
Coulomb (SR) Coulomb (LR) Potential Kinetic En. Total Energy
0.00000e+00 0.00000e+00 -3.39736e+02 1.08668e+02 -2.31067e+02
Temperature Pressure (bar) dVremain/dl dVvdw/dl Constr. rmsd
3.00453e+02 0.00000e+00 0.00000e+00 6.43777e+00 0.00000e+00
Constr.2 rmsd
0.00000e+00

and with free energy off:

Bond          Angle    Proper Dih.  Improper Dih.GB Polarization
3.99457e+01 7.13394e+01 1.14065e+02 2.81320e+00 -1.02303e+02
Nonpolar Sol. LJ-14 Coulomb-14 LJ (SR) LJ (LR)
1.57390e+01 2.49700e+01 1.64072e+02 -2.16785e+01 0.00000e+00
Coulomb (SR) Coulomb (LR) Potential Kinetic En. Total Energy
-6.54423e+02 0.00000e+00 -3.45459e+02 8.69806e+01 -2.58479e+02
Temperature Pressure (bar) Constr. rmsd Constr.2 rmsd
2.40490e+02 0.00000e+00 0.00000e+00 0.00000e+00

Note the Coulomb (SR) energies, which are identically zero when free energy is on, and nonzero when off.

When computing a hydration free energy with annihilation (which is all I've tested so far) for a molecule with more than 1-4 interactions, the normal internal electrostatic interactions have a large contribution to the free energy for annihilation in vacuum. But in GB, the SR Coulomb interactions get turned off, and so this contribution is missing. This makes hydration free energies be incorrect by something on the order of the strength of the internal coulomb SR interactions, which in this case are a couple of hundred kJ/mol.

I do not know whether the other energies are correct or not -- this is the leading source of error, though.

Also, for solvation of a small molecule alone, it is possible that this bug would not affect decoupling, since all of the Coulomb contributions (except 1-4s) may be short range Coulomb.

Attached please find a test system which should reproduce this problem. But indeed, you should basically be able to take any system in GB and turn the free energy code on and see that the short range Coulomb energies end up being reported as zero.

It is not clear to me how easy it will be to make the free energy code work PROPERLY with GB solvent (it's entirely possible there are other, larger issues here) but I'd suggest that for 4.6.2 release, there be at least the following changes:
1) free_energy = yes should not make GB short range coulomb energies be incorrectly reported as zero
2) If there are larger reasons GB shouldn't be used with free_energy = yes, grompp should stop the user from doing this

implicit.tar.gz (7.04 KB) David Mobley, 05/01/2013 02:04 AM

Associated revisions

Revision 073fa1ed (diff)
Added by Mark Abraham over 4 years ago

Introduce fatal error for GB with FE

Old code silently computed zeroes for Coulomb(SR) when using
Generalized Born and the free energy kernel, because of the
inappropriate use of a default branch of a switch statement.

Fixed the default branch for icoul and ivdw switch statements to issue
gmx_incons (similar to usage in the Verlet kernel code). Now we invoke
explicit code paths for all legal values of icoul and ivdw. Added
gmx_fatal for un-implemented GB case (similar to the gmx_fatal for
Buckingham with this free energy kernel).

Partial fix for #1237

Change-Id: I30684aed9fba9beac9b1b87eb99e2cbc7e7e374d

Revision 677062b5 (diff)
Added by Mark Abraham about 4 years ago

Introduce fatal error for GB with FE

Old code silently computed zeroes for Coulomb(SR) when using
Generalized Born and the free energy kernel, because of the
inappropriate use of a default branch of a switch statement.

Fixed the default branch for icoul and ivdw switch statements to issue
gmx_incons (similar to usage in the Verlet kernel code). Now we invoke
explicit code paths for all legal values of icoul and ivdw. Added
gmx_fatal for un-implemented GB case (similar to the gmx_fatal for
Buckingham with this free energy kernel).

Partial fix for #1237

Change-Id: I30684aed9fba9beac9b1b87eb99e2cbc7e7e374d

History

#1 Updated by Mark Abraham over 4 years ago

  • Status changed from New to Accepted

Reproduced with pre-4.6.2 git version.

#2 Updated by Mark Abraham over 4 years ago

  • Status changed from Accepted to In Progress
  • Assignee changed from Berk Hess to Mark Abraham
  • % Done changed from 0 to 50

It turns out that a switch statement in nb_free_energy.c had a default case to treat icoul==GMX_NBKERNEL_ELEC_NONE and this was also run for GMX_NBKERNEL_ELEC_GENERALIZED_BORN, generating the zeroes. Clear reason to not use default cases if ever we saw one!

A minimal acceptable solution is
  • to change the default branches to gmx_incons("something") to avoid this kind of bug in future,
  • have explicit case branches for GMX_NBKERNEL_*_NONE, and
  • gmx_fatal for GB (as the current code does for FE+Buckingham)
    and this is probably all we can do in time for 4.6.2.

I have started work on a proper GB case branch, but as I am a noob at both FE and GB, we will need such a patch to get serious review and testing (and a new test case).

#3 Updated by Mark Abraham over 4 years ago

  • Status changed from In Progress to Fix uploaded
  • Target version changed from 4.6.2 to 4.6.3

#4 Updated by David Mobley over 4 years ago

When you guys get to the point that you're ready (i.e. when you think the issue is resolved), I'm happy to test further. This specific test case only showed me that things were wildly wrong. However, there's a convenient property of GB (relative to explicit solvent) which makes it pretty easy to test that the free energy code is working properly, and that is that, for most small molecules, it's possible to run a "vanilla" (non free energy) simulation in GB and re-evaluate the energies with GB off, and compute the hydration free energy by applying BAR/MBAR to these energies. Normally this approach alone (for small relatively rigid molecules) is enough to get the hydration free energy, without using intermediates. So then I can repeat the same procedure using the free energy code and intermediate lambda values and compare...

Anyway, let me know when you think it's ready for that level of testing.

#5 Updated by Mark Abraham over 4 years ago

It was fairly easy to fix the short-range Coulomb potential and (presumably) forces, but I think someone will need to produce the right equations for the three components of G_solv before we can solve the other bits.

With mdrun -rerun doing a single-point calculation on the equil_npt2.0.gro file in the above tarball:

With no-FE .tpr

           Bond          Angle    Proper Dih.  Improper Dih.GB Polarization
    1.88790e+01    6.39214e+01    1.08525e+02    1.11358e+01   -1.37088e+02
  Nonpolar Sol.          LJ-14     Coulomb-14        LJ (SR)        LJ (LR)
    1.57793e+01    2.94219e+01    1.41759e+02   -2.22076e+01    0.00000e+00
   Coulomb (SR)   Coulomb (LR)      Potential    Kinetic En.   Total Energy
   -5.79862e+02    0.00000e+00   -3.49736e+02    9.45046e+01   -2.55232e+02
    Temperature Pressure (bar)   Constr. rmsd  Constr.2 rmsd
    2.61293e+02    0.00000e+00    3.51925e-06    2.87866e-06

With FE .tpr and code from commit d13fc483199 from release-4-6

           Bond          Angle    Proper Dih.  Improper Dih.GB Polarization
    1.88790e+01    6.39214e+01    1.08525e+02    1.11358e+01   -7.14803e+02
  Nonpolar Sol.          LJ-14     Coulomb-14        LJ (SR)        LJ (LR)
    2.33285e+01    2.94219e+01    1.41759e+02   -2.22076e+01    0.00000e+00
   Coulomb (SR)   Coulomb (LR)      Potential    Kinetic En.   Total Energy
    0.00000e+00    0.00000e+00   -3.40040e+02    9.44831e+01   -2.45557e+02
    Temperature Pressure (bar)    dVremain/dl       dVvdw/dl   Constr. rmsd
    2.61233e+02    0.00000e+00    0.00000e+00    1.07118e+01    2.99180e-06
  Constr.2 rmsd
    3.23584e-06

With FE .tpr and code from https://gerrit.gromacs.org/2355 patch 1

           Bond          Angle    Proper Dih.  Improper Dih.GB Polarization
    1.88790e+01    6.39214e+01    1.08525e+02    1.11358e+01   -2.83251e+02
  Nonpolar Sol.          LJ-14     Coulomb-14        LJ (SR)        LJ (LR)
    2.33285e+01    2.94219e+01    1.41759e+02   -2.22076e+01    0.00000e+00
   Coulomb (SR)   Coulomb (LR)      Potential    Kinetic En.   Total Energy
   -5.79862e+02    0.00000e+00   -4.88349e+02    9.63036e+01   -3.92046e+02
    Temperature Pressure (bar)    dVremain/dl       dVvdw/dl   Constr. rmsd
    2.66266e+02    0.00000e+00    5.79862e+02    1.07118e+01    3.19118e-06
  Constr.2 rmsd
    3.22377e-06

Coulomb (SR) is now right. It's possible I have not gotten the short-range Coulomb scalar force right yet.

Pol has changed, so I guess the computed Born radii have changed, but there's further stuff that is not working. No progress on non-pol yet.

I don't think there's any more I can do until someone puts pen to paper to show how G_solv and derivatives vary with lambda. (And even then I might not prioritise coding it.) The relevant code looks like it is in src/mdlib/genborn.c.

#6 Updated by Mark Abraham over 4 years ago

  • Assignee deleted (Mark Abraham)
  • Target version deleted (4.6.3)

#7 Updated by Mark Abraham over 4 years ago

  • Status changed from Fix uploaded to In Progress

#8 Updated by Mark Abraham over 4 years ago

  • Status changed from In Progress to Blocked, need info

#9 Updated by Michael Shirts over 4 years ago

I can volunteer to help work out the theory when we're at that point.

#10 Updated by Mark Abraham over 4 years ago

I think we're at the point we need a theory/derivation contribution. The GROMACS manual has the GB equations AFAICS. I don't know whether it has the relevant FE equations.

#11 Updated by Michael Shirts over 4 years ago

Mark Abraham wrote:

I think we're at the point we need a theory/derivation contribution. The GROMACS manual has the GB equations AFAICS. I don't know whether it has the relevant FE equations.

I'll probably have the time in a week or so, but I don't this week. Let's make sure this gets done for 4.6.3, though.

#12 Updated by Rossen Apostolov over 3 years ago

  • Assignee set to Michael Shirts

do we want this fixed eventually?

#13 Updated by Erik Lindahl about 3 years ago

  • Tracker changed from Bug to Feature

Gromacs has never supported FE with GB, so I'm changing this to a future feature.

#14 Updated by Mark Abraham about 3 years ago

  • Status changed from Blocked, need info to Accepted
  • Assignee deleted (Michael Shirts)
  • Priority changed from Normal to Low

If anyone ever wants to work on this, there's some code at https://gerrit.gromacs.org/#/c/2355/, based on release-4-6 branch.

Also available in: Atom PDF