Unphysical conformations in decoupled free energy simulation
This follows the discussion http://gromacs.5086.x6.nabble.com/Unphysical-conformations-in-decoupled-free-energy-simulation-tt5010337.html
I am trying to compute the free energy of hydration for cellobiose (a beta (1-4) glucose dimer) using BAR in Gromacs 4.6.3. However, I encounter a problem.
I find that the vacuum conformations of the molecule in a regular vacuum simulation differ from the conformations in the decoupled simulation in the free energy case i.e., with the additional mdp entries:
free-energy = yes
init-lambda = 0
couple-lambda0 = none
couple-moltype = solute
couple-intramol = no
Here is a histogram of the dihedral angles of the glycosidic linkage in
the vacuum case
https://dl.dropboxusercontent.com/u/70358077/reg.pdf (it stays in the global free energy minimum)
and this is the decoupled free energy case (same starting conformation in the global free energy minimum)
I understand that Gromacs replaces the intramolecular interactions with explicit pair interactions. Therefore, I had to increase table-extension but that did not change much. I was thinking that maybe this could be a problem specific to this topology (a GLYCAM06h conversion from Amber), however, I do not see how this can occur. I hope someone has an idea what is going wrong.
This is sort of a minimal example. The same problem occurs in a regular simulation when decoupling from water using multiple lambdas.
All data is attached in the tarball.
Fixes a problem with pair type 2 interactions with free energy
Pair type 2 interactions, which should remain on regardless
of couple-intramol=yes, were being turned off. Currently, when free
energies were turned on, they were just ignored, because the (empty)
pair one 1 type list was copied over them. This fix adresses
this problem by adding onto the list instead of copying it over.
#1 Updated by Michael Shirts about 4 years ago
- Category set to mdrun
- Status changed from New to In Progress
I've think figured out the root case -- currently free energy calculations can't handle the pair = 2 functional type, so those are completely ignored, rather than just keeping them on, or issuing a warning/error.
Clearly this is not the desired behavior. I'll have to dig in and figure out if the 'solution' is to issue an error that this topology can't be handled with free energy, or to actually handle the case properly.
#2 Updated by Michael Shirts about 4 years ago
OK, A bit more information here:
With pairs of type 2, then the information is written into the F_LJC14_Q pairlist.
Then a problem occurs at toppush.c, lines 2524, in the function convert_pairs_to_pairsQ, called when we are using couple-moltype.
/* Copy the pair list to the pairQ list */
plist[F_LJC14_Q] = plist[F_LJ14];
With pairs function type 1, this moves the list in plist[F_LJ14] over to plist[F_LJC14_Q], and then stuff happens.
With pairs function type 2, this erases everything in the F_LJC14_Q pairlist, eliminating these interactions. Nothing is currently in plist[F_LJ14].
I'd appreciate any insights as to how this should be fixed . . .