Project

General

Profile

Feature #1340

LJ and Coul 1-4 are not turned off when coupl-intramol = yes

Added by Michael Shirts about 7 years ago. Updated over 6 years ago.

Status:
Closed
Priority:
Normal
Assignee:
Category:
mdrun
Target version:
Difficulty:
uncategorized
Close

Description

Not sure if it's a bug -- I would consider it such, but maybe it's just poorly documented.

Right now, if one sets coupl-intramol = yes, LJ and Coul 1-4's are left on at all lambda.

In most cases, people would expect that if the intermolecular nonbondeds are coupled, then
the LJ 1-4's are would be turned off in the decoupled end state. The LJ 1-5's, LJ 1-6's intramolecular
interactions are turned off, for example -- what's special about the 1-4's?

Another reason to change this effect is that if you were to set the LJ and Coul parameters
equal to zero manually in the topologies, then LJ and Coul 1-4's would be zero.
Coupl-intramol is supposed to replicate this behavior.

I've assigned this to Berk, since he wrote this code, but it can be switched to me as required.

prod.mdp (4.68 KB) prod.mdp Michael Shirts, 09/21/2013 02:46 AM
system_GMX.top (14.1 KB) system_GMX.top Michael Shirts, 09/21/2013 02:46 AM
equil_npt.gro (135 KB) equil_npt.gro Michael Shirts, 09/21/2013 02:46 AM

History

#1 Updated by Michael Shirts almost 7 years ago

It would be good to get feedback on this point. It's either a bug (my opinion) or need to be communicated to users much more strongly. It would be good to get this resolved for the next release.

#2 Updated by Michael Shirts almost 7 years ago

Thoughts, anyone? Is there a reason the behavior is like this?

#3 Updated by Michael Shirts almost 7 years ago

  • Status changed from New to Fix uploaded

Fix uploaded that should properly cause them to be lambda dependent.

#4 Updated by Mark Abraham almost 7 years ago

I don't understand the wording of the documentation of couple-intramol:

couple-intramol:
no
All intra-molecular non-bonded interactions for moleculetype couple-moltype
are replaced by exclusions and explicit pair interactions. In this manner the decou-
pled state of the molecule corresponds to the proper vacuum state without periodicity
effects.
yes
The intra-molecular Van der Waals and Coulomb interactions are also turned on/off.
This can be useful for partitioning free-energies of relatively large molecules, where
the intra-molecular non-bonded interactions might lead to kinetically trapped vacuum
conformations. The 1-4 pair interactions are not turned off.

The "no" case transforms a set of interactions into new forms, but doesn't say how those forms will vary with lambda. The "yes" case doesn't mention transformation, but implies variation with lambda. Berk, can you clarify this, please?

#5 Updated by Berk Hess almost 7 years ago

With "no", the intra-molecular intermolecular interactions are pure LJ and pure Coulomb. No cut-off is applied and there is no lambda dependence.
In the "yes" case, intra-molecular interactions couple in exactly the same way as the intermolecular ones.
The 1-4 pair interactions are not turned off, as those are considered bonded interactions, not non-bonded interactions (they are part of the dihedral potentials).
If you can phrase this better in the manual, pleas tell me!

#6 Updated by Michael Shirts almost 7 years ago

Berk Hess wrote:

The 1-4 pair interactions are not turned off, as those are considered bonded interactions, not non-bonded interactions (they are part of the dihedral potentials).

I would say a lot of people (the majority) consider nonbonded to refer to the the epsilon and sigma nonbonded parameters, and that when the intermolecular parameters are turned off, these should be as well. I mean, the 1-5's and 1-6's contribute to the dihedral potentials as well, though indirectly. What's that special about the 1-4's other than the magnitude of the effect?

If one wants to leave the 1-4's on, then the solution would be to use 2 pairs in the topology, which ignore the intramolecular coupling. But there is no current way to do what I would think is the 'standard' behavior without adding a bunch of atoms types.

I think there are a range of ways to resolve this, but it appears we are talking primarily about convention now, and that requires looking at more general definitions of what people. I may improperly be considering myself in the majority, but we at least do need to evaluate what the expected usage would be.

Note we could make the behavior more general by making Couple-Intramol more than a binary parameter, so that multiple behaviors can be supported.

I'll have David Mobley comment here so we get his perspective. Might be a little while since he's on West Coast.

#7 Updated by Erik Lindahl almost 7 years ago

The historical reason for introducing 1,4 parameters was originally that the torsion potential was an approximation of the effect of the rotation around the two central atoms (2,3), while the 1,4 part described the influence of outer atoms. This way one did not have to provide as many different torsions, and even today I would argue it plays this role to a smaller extent since we do not reparameterize all torsions when partial charges change.

I have thought about whether we should change the nomenclature bonded/nonbonded interactions to listed vs. neigborlist interactions to make this clearer. This would be a better reflection of the setup that 1,4 interactions are first excluded and then specified explicitly (although pdb2gmx does that for you) in contrast to 1,5 or 1,6 that are only included through neighborlists.

#8 Updated by David Mobley almost 7 years ago

(I tried to e-mail this, but I think it doesn't take updates by e-mail, or at least not from me. Sorry).

Hi, all,

Just to chime in on this thread: I think the historical REASON for introducing 1,4 parameters is entirely a separate issue from whether these are considered bonded or nonbonded interactions. In general they are considered nonbonded interactions -- i.e., in the teaching context, LJ and Coulomb interactions are always listed as nonbonded (even 1,4's) and bonded interactions are just bonds, angles, and torsions. This suggests that if we're saying couple-intramol = yes corresponds to turning off intramolecular nonbonded interactions, the logical way to handle this is by turning off LJ 1,4s and Coulomb 1,4s.

Secondarily, the "couple-intramol" I believe ought to correspond roughly to the two main current approaches to free energy calculations seen in the literature:
1) "Decoupling" (couple-intramol = no) - maintain all internal LJ and electrostatic interactions, modulating only those in the environment. (Prior to adding couple-intramol, this couldn't be achieved for charges, and for LJ could be achieved via a hack in the topology file if one was careful)
2) "Annihilation" (couple-intramol = yes) - maintain NO internal LJ and electrostatic interactions. Can be achieved without couple-intramol by making the B state charges and LJ be zero.

I have seen ZERO publications where option #2 involves leaving on internal 1,4 interactions within the molecule, and many where it does not. I believe this argues for making sure "couple-intramol = yes" corrresponds to turning off all internal LJ and Coulomb interactions. This also makes sense, as otherwise, it's possible to follow one approach by modifying the B state explicitly to turn off charges and LJ (annihilation), and an entirely different approach by using couple-intramol = yes.

If some people feel strongly that there ought to be a way to turn off intramolecular interactions except 1,4s, I'm open to adding a third couple-intramol option. Perhaps their names then ought to be changed to:
a) couple-intramol = decoupling (option 1 above)
b) couple-intramol = annihilation (option 2 above)
c) couple-intramol = partial (turn off internal LJ and Coulomb interactions for atoms more distant than 1,4)

One further (perhaps minor) reason to include an 'annihilation' option is that this makes it a lot easier to test couple-intramol. Specifically, one should be able to reprocess a stored trajectory with couple-intramol = yes and lambda set to put us at the B state, and reproduce the same energies as with a topology set to put us at the B state with no charges and LJ. If I understand right, if we maintain internal LJ and Coulomb 1,4s in annihilation, it will be a lot harder to do that test.

Hope that helps.

#9 Updated by Berk Hess almost 7 years ago

I think there are two separable issues here:
1) a phrasing/nomenclature issue
2) useful setups for free-energy calculations.

My thoughts:
1) At the mdrun level, Gromacs doesn't actually know about 1-4 interactions. We have only explicitly listed pairs in the topology. We do have an explicit exclusion setup which usually runs up to 1-4. Thus the question would be if we want the explicitly listed pairs to be coupled or not.

2) 1-4 interaction contribute very significantly to the dihedral potential and I think these should not be turned off, as the dihedral distributions would then change significantly and you want to avoid changing distributions in free-energy calculations. If you think the 1-5 and 1-6 are beneficial, then you can use couple-intramol=no.

#10 Updated by David Mobley almost 7 years ago

Berk,

In terms of #2, there are different schools of thought on this and it is fairly controversial. I've used the "full annihilation" approach I describe on more than 700 molecules for hydration free energy calculations, and it works fairly reliably for me. Other people I know (including Bill Swope and Julia Rice) have revisited many of my calculations (a few hundred) with the "full decoupling" approach, and in general reproduce my values. However, the two approaches seem to have different challenges. As you note, removing these interactions can affect the free energy landscape you sample in the decoupled state, which could be bad (though as long as you sample it adequately, it in fact is fine). But, retaining these interactions also retains larger energy barriers between minima in the torsional landscape, for example. So you could imagine that in some cases, you might get better sampling (especially when using lambda exchange) without the 1,4 interactions, and that is exactly what we've seen.

So, I would say at present, there is no clear conclusion about whether it is better/more efficient to retain 1,4s (or other internal interactions).

Since this in fact is an active topic of research, I think this is a further argument in favor of making it easy to use all three approaches, at least until someone is able to demonstrate whether one is more efficient than the others in general.

To reiterate -- while I've certainly talked to other people who agree with you in general about #2, this is in fact not clear to me, and we've looked at this to a limited extent and so far the data doesn't support that argument in general.

#11 Updated by Michael Shirts almost 7 years ago

As mentioned before, I side with David here.

I'm perfectly fine with adding something like 'couple-intramol=not14' which would preserve Berk's suggested behavior.

#12 Updated by Berk Hess almost 7 years ago

I think we shouldn't change code for nomenclature issues, we should change the nomenclature (which in this case is controversial anyhow).
For the actual usefulness, I agree with David's approach that different things could be different in different situations (or at least at present it is not clear if there is a single optimal approach). Thus having all three options is probably best. But I would suggest to leave this change for 5.0. We should think about a good naming scheme; I would like to avoid negative names such as "not14".
We can clarify the manual for 4.6.

#13 Updated by Michael Shirts almost 7 years ago

I think one reason for changing it now is that it's not currently possible to do it the most common way, and it is possible to do the less common way (use function type 2). David and I can't do it they way we'd like to, and nobody (as far as I can tell) is doing the way that is implemented.

I agree avoiding negatives. Perhaps couple-intermol='Unmodified14' would be better than 'Not14'?

Additionally, the the code needs fixing anyway, because it overwrites function type 2.

Obviously the manual needs to be much clearer here as well.

#14 Updated by Michael Shirts almost 7 years ago

  • Assignee changed from Berk Hess to Michael Shirts

#15 Updated by Berk Hess almost 7 years ago

I think we should never change mdp options in minor revisions, but if other don't object, I can reconsider.

I would suggest (since strictly speaking we don't have 1-4 interactions in the topology):
nonbonded: (de)couples nonbonded, but not pairs (usually these are the 1-4 interactions).
nonbonded-pairs: (de)couples nonbonded and pairs (usually these are the 1-4 interactions).
(we need to explain more extensively, of course)

#16 Updated by Berk Hess almost 7 years ago

  • Assignee changed from Michael Shirts to Berk Hess

I discussed this with Mark and we don't see why functionality that has been in Gromacs for 5 years needs to now be changed right now in a minor revision. I would suggest to add the requested functionality in 5.0.
Also I just actually read the mdp documentation and it clearly states that 1-4 interactions are not coupled. Since the option name "yes" doesn't even mention non-bondeds, I think the option is not even confusing as it is (which is not to say that it can be improved, but as said before, we should not change option names in minor revisions).

#17 Updated by Erik Lindahl almost 7 years ago

  • Tracker changed from Bug to Feature
  • Priority changed from High to Normal
  • Target version changed from 4.6.4 to 5.0

Agree with Berk. I fully understand (and agree) that it would be good to change it in 5.0, but it is not a bug since the code does exactly what it says in the documentation. I've changed the status to a feature request, with target 5.0. We will have a first beta of that in less than two months, so it will be possible to get the code out there almost as fast as a potential 4.6.5.

#18 Updated by Michael Shirts almost 7 years ago

Erik Lindahl wrote:

Agree with Berk. I fully understand (and agree) that it would be good to change it in 5.0, but it is not a bug since the code does exactly what it says in the documentation. I've changed the status to a feature request, with target 5.0. We will have a first beta of that in less than two months, so it will be possible to get the code out there almost as fast as a potential 4.6.5.

Sounds fine. I'll get the actual bug (#1315) involving this code fixed and checked in today.

#20 Updated by Rossen Apostolov over 6 years ago

Bump - can we close this one?

#21 Updated by Rossen Apostolov over 6 years ago

  • Status changed from Fix uploaded to Feedback wanted

#22 Updated by Erik Lindahl over 6 years ago

  • Status changed from Feedback wanted to Closed

No feedback, so assume Michael is happy.

Also available in: Atom PDF