Bug #1315

Unphysical conformations in decoupled free energy simulation

Added by Joerg Sauter about 7 years ago. Updated almost 7 years ago.

Target version:
Affected version - extra info:
Affected version:


This follows the discussion

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 (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.

fe.tar.gz (2.38 MB) fe.tar.gz Joerg Sauter, 08/06/2013 10:03 AM

Associated revisions

Revision 8839a4f9 (diff)
Added by Michael Shirts about 7 years ago

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.

Fixes #1315

Change-Id: I240479a8dc083f7a355917ed9f74f4337fa3448f


#1 Updated by Michael Shirts about 7 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 7 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 . . .

#3 Updated by Michael Shirts about 7 years ago

  • Status changed from In Progress to Fix uploaded
  • Assignee set to Michael Shirts

#4 Updated by Michael Shirts about 7 years ago

Fix uploaded that makes pair function type 2 always on, regardless of coupling type.

#5 Updated by Michael Shirts about 7 years ago

  • Status changed from Fix uploaded to Resolved
  • % Done changed from 0 to 100

#6 Updated by Rossen Apostolov almost 7 years ago

  • Status changed from Resolved to Closed

Also available in: Atom PDF