[gmx-developers] free energy calculations with restraints

[Shirts, Michael (mrs5pt)]

The other thing I would say is that the couple-mol code is a little
difficult to work with----there's known quirks, like handling the
intramolecular dispersion in a way that's slightly off from free energy off,
as well as some new quirks that I'll try to post about later this
week---that it would actually be cleaner to add clearly defined
intermolecular restraints than add complexity to the coupling feature.

Best,
~~~~~~~~~~~~
Michael Shirts
Assistant Professor
Department of Chemical Engineering
University of Virginia
michael.shirts at virginia.edu
(434)-243-1821


> From: Floris Buelens <floris_buelens at yahoo.com>
> Reply-To: Floris Buelens <floris_buelens at yahoo.com>, Discussion list for
> GROMACS development <gmx-developers at gromacs.org>
> Date: Tue, 20 Aug 2013 08:35:15 -0400
> To: Berk Hess <hess at kth.se>, Discussion list for GROMACS development
> <gmx-developers at gromacs.org>
> Subject: Re: [gmx-developers] free energy calculations with restraints
> 
> Hi Berk,
> 
> I think we might have gone out of sync there - I sent the message below before
> reading yours from this morning. From the misplaced context it maybe sounded
> like I was disagreeing with your '[ intermolecular-interactions ]' suggestion,
> which I'm definitely not - this would clearly be preferable and probably
> useful to a wider audience.
> 
> 
> 
> ----- Original Message -----
> From: Berk Hess <hess at kth.se>
> To: Floris Buelens <floris_buelens at yahoo.com>; Discussion list for GROMACS
> development <gmx-developers at gromacs.org>
> Cc: 
> Sent: Tuesday, 20 August 2013, 12:04
> Subject: Re: [gmx-developers] free energy calculations with restraints
> 
> Hi,
> 
> We could consider this.
> I don't know how you did this now, but this would need an extra mdp
> option, something like couple-group, to avoid one option with a moltype
> or index group as an input.
> Then we would also need to decide what to do with a group which is part
> of one molecule of a certain moleculetype with multiple copies. Disallow
> this? Or duplicate the moltype and only modify one copy (as is done now
> for qmmm)?
> 
> Cheers,
> 
> Berk
> 
> On 08/20/2013 09:24 AM, Floris Buelens wrote:
>> Hi,
>> 
>> I would also still argue for augmenting the current couple-moltype
>> functionality with this modification. I can see this is a niche requirement,
>> which of course does speak against inclusion. But in my view the replacement
>> code is no less elegant than the current implementation (modification of five
>> self-contained functions at the pre-processing stage, maybe 50-odd extra
>> lines), we can maintain full backward compatibility, and a single simple
>> check will remove any potential for accidental misuse.
>> 
>> 
>> 
>> ________________________________
>> From: David Mobley <dmobley at uci.edu>
>> To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
>> Cc: Floris Buelens <floris_buelens at yahoo.com>; Michael R. Shirts
>> <michael.shirts at virginia.edu>
>> Sent: Monday, 19 August 2013, 20:02
>> Subject: Re: [gmx-developers] free energy calculations with restraints
>> 
>> 
>> 
>> Hi, Berk,
>> 
>> I agree that this solution would be the "proper" one and would be nice.
>> However, it's been a long time coming (this is an issue I've been raising for
>> approximately 5 years, and there's no progress; normally no one even
>> answers). Currently, for the restraints I need to use, I *have* to merge the
>> ligand into my protein topology, inconvenient or not. Doing this also means
>> there are some features I simply can't use (couple-moltype for example).
>> 
>> It strikes me that Floris's solution is a good interim one, in that it at
>> least solves most of the typical use cases we're facing, and has the
>> advantage of being *already implemented*.
>> 
>> Is there any way we could get this into the main GROMACS, at least until
>> someone gets time to implementing the "proper" solution (which could be years
>> from now)?
>> 
>> Thanks!
>> David
>> 
>> 
>> 
>> 
>> 
>> On Mon, Aug 19, 2013 at 1:10 AM, Berk Hess <hess at kth.se> wrote:
>> 
>> Hi,
>>> I fully understood this point, maybe my reply wasn't well structured.
>>> I was trying to argue that the proper solution would be to implement
>>> inter-molecular restraints, instead of introducing an option which could me
>>> misused. Also merging a ligand topology into your protein topology is very
>>> inconvenient. This problem with the proper solution is, of course, that
>>> someone will need to implement inter-molecular restraints/potentials.
>>> 
>>> One issue used to be that we corrected molecules for PBC before calculating
>>> bonded interactions. But in 4.6 this is usually not faster and not used any
>>> more. That makes it easier to treat intra- and inter-molecular interactions
>>> the same way at the mdrun level.
>>> 
>>> Cheers,
>>> 
>>> Berk
>>> 
>>> 
>>> On 08/19/2013 09:55 AM, Floris Buelens wrote:
>>> 
>>> Hi Berk,
>>>> I think some clarification is needed. What we're talking about is only
>>>> relevant in the context of non-bonded interactions, typically between a
>>>> small molecule and a protein. To access the binding affinity of a ligand
>>>> through alchemical methods, it's useful to switch off only the interactions
>>>> of the ligand with its environment while maintaining intra-ligand
>>>> potentials - this is what's provided by the 'couple-moltype' code path.
>>>> 
>>>> The limitation that we're trying to circumvent comes from the fact that as
>>>> you scale down interactions with the environment, the ligand is no longer
>>>> held in the binding site. To counteract this, it's necessary to perturb in
>>>> restraining potentials as the intramolecular nonbonded interactions are
>>>> perturbed out. The most practical method makes use of a single distance
>>>> restraint, two angle and three dihedral restraints, whose effect can later
>>>> be accounted for analytically.
>>>> 
>>>> The current code only allows decoupling of a whole logical 'molecule' as
>>>> specified in the topology file. This precludes the use of regular (fully
>>>> perturbation-aware) bonded interactions. As far as I'm aware, the pull code
>>>> doesn't provide what's required. Inter-molecular bonded interactions would
>>>> be a great general solution but presumably won't show up any time soon.
>>>> 
>>>> The workaround here is instead to represent two physical molecules (e.g.
>>>> protein and ligand) as a single logical molecule (in the topology file), so
>>>> regular bonded potentials can be applied. Current Gromacs doesn't allow the
>>>> decoupling ('couple-moltype') code to be used in this scenario. My
>>>> modification allows you to target the ligand by its residue name and get
>>>> this functionality back.
>>>> 
>>>> Decoupling a single residue of a multi-residue chain is indeed probably not
>>>> correct in this framework (I did call it 'weird' in my last message :-) ).
>>>> An extra check would block users from trying this.
>>>> 
>>>> I hope that clarifies the problem we're trying to solve. I agree this will
>>>> be useful a relatively small number of users, but on the other hand it's a
>>>> very unobtrusive, self-contained modification which can maintain full
>>>> backwards compatibility.
>>>> 
>>>> thanks,
>>>> 
>>>> Floris
>>>> 
>>>> 
>>>> ----- Original Message -----
>>>> From: Berk Hess <hess at kth.se>
>>>> To: Floris Buelens <floris_buelens at yahoo.com>; Discussion list for GROMACS
>>>> development <gmx-developers at gromacs.org>
>>>> Cc:
>>>> Sent: Monday, 19 August 2013, 9:03
>>>> Subject: Re: [gmx-developers] free energy calculations with restraints
>>>> 
>>>> Hi,
>>>> 
>>>> There is the rotational pull code, but that might not provide the exact
>>>> functionality you want.
>>>> Already for some time we have been discussing inter-molecular
>>>> interactions, by specifying molecule type, molecule index and atom
>>>> index, but there are no concrete plans for implementing this yet.
>>>> 
>>>> An option for decoupling part of a molecule indeed sound useful. But in
>>>> practice you always need to replace that part by something else, at
>>>> least a hydrogen, and modify some potentials of connecting atoms, so I
>>>> don't know how generally useful such an option is.
>>>> 
>>>> Cheers,
>>>> 
>>>> Berk
>>>> 
>>>> On 08/19/2013 07:56 AM, Floris Buelens wrote:
>>>> 
>>>> I suppose that would make sense. The changes are fairly minor (more or less
>>>> just the function convert_moltype_couple and the four functions it calls)
>>>> and mainly consist of a bit of extra juggling with nonbonded exclusions
>>>> during preprocessing. The current functionality (decouple a whole molecule)
>>>> could be handled as a special case of the new code. A check for bonds to
>>>> the rest of the structure could stop people from trying weird stuff like
>>>> decoupling residues from a chain.
>>>>> 
>>>>> 
>>>>> 
>>>>> ________________________________
>>>>> From: David Mobley <dmobley at gmail.com>
>>>>> To: Floris Buelens <floris_buelens at yahoo.com>
>>>>> Cc: Discussion list for GROMACS development <gmx-developers at gromacs.org>
>>>>> Sent: Friday, 16 August 2013, 18:30
>>>>> Subject: Re: [gmx-developers] free energy calculations with restraints
>>>>> 
>>>>> 
>>>>> 
>>>>> This does sound useful, though it would be more useful if this could be
>>>>> implemented into the main GROMACS rather than a separate code (since
>>>>> otherwise it will go away as GROMACS is further developed).
>>>>> 
>>>>> Would this be a possibility going forward, devels?
>>>>> 
>>>>> Thanks!
>>>>> 
>>>>> 
>>>>> 
>>>>> 
>>>>> On Fri, Aug 16, 2013 at 3:07 AM, Floris Buelens <floris_buelens at yahoo.com>
>>>>> wrote:
>>>>> 
>>>>> Hi David,
>>>>> 
>>>>> I have the same requirement as you, but I've gone about it slightly
>>>>> differently. I've hacked the couple-moltype code path to allow decoupling
>>>>> of a specific residue (identified by name) instead of a molecule. This
>>>>> allows the residue of interest to be part of another molecule block so you
>>>>> can set up distance, angle and dihedral restraints using regular bonded
>>>>> interactions with full perturbation support.
>>>>>> If this is useful to you or to anyone else, let me know and I'll be happy
>>>>>> to share.
>>>>>> 
>>>>>> best,
>>>>>> 
>>>>>> Floris
>>>>>> 
>>>>>> ________________________________
>>>>>> From: David Mobley <dmobley at gmail.com>
>>>>>> 
>>>>>> To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
>>>>>> Sent: Friday, 14 June 2013, 21:47
>>>>>> Subject: [gmx-developers] free energy calculations with restraints
>>>>>> 
>>>>>> 
>>>>>> 
>>>>>> 
>>>>>> Hi, devs,
>>>>>> 
>>>>>> 
>>>>>> I'm writing with an issue relating to the interplay of new free energy
>>>>>> features with restraints.
>>>>>> 
>>>>>> I'm very much appreciating some of the new free energy features in
>>>>>> gromacs, such as the 'couple-moltype' option which provides a way to set
>>>>>> up decoupling or annihilation of a specific molecule via free energy
>>>>>> calculations without having to edit the topology file directly. This is
>>>>>> especially great when it comes to decoupling -- charge decoupling was not
>>>>>> previously possible via topology file editing, and vdW decoupling took
>>>>>> substantial manipulation of the topology file.
>>>>>> 
>>>>>> However, for binding free energies, my work employs orientational
>>>>>> restraints between the protein and ligand. I need to be able to impose
>>>>>> both dihedral and angle restraints on angles between the protein and
>>>>>> ligand. Currently, I do this using angle-restraints and
>>>>>> dihedral-restraints. This requires that both the protein and ligand be
>>>>>> within the same logical 'molecule', which (unfortunately) means that I
>>>>>> can't make use of the new free energy features above, since
>>>>>> couple-moltype has to apply to a whole molecule, not just part of a
>>>>>> molecule.
>>>>>> 
>>>>>> So, my I see two possible solutions to the problem, and hence have these
>>>>>> questions:
>>>>>> 1) Can dihedral and angle restraints be applied via the pull code? If
>>>>>> not, are there any plans to add that?
>>>>>> 2) Alternatively, what about modifying the restraints code so it uses (or
>>>>>> at least optionally allows) absolute atom numbering, rather than
>>>>>> numbering within a specific molecule, thus allowing restraints
>>>>>> (angle-restraints and dihedral-restraints) to be applied between
>>>>>> 'molecules'?
>>>>>> 
>>>>>> Thanks!
>>>>>> David
>>>>>> 
>>>>>> 
>>>>>> 
>>>>>> --
>>>>>> David Mobley
>>>>>> dmobley at gmail.com
>>>>>> 949-385-2436
>>>>>> 
>>>>>> 
>>>>>> --
>>>>>> gmx-developers mailing list
>>>>>> gmx-developers at gromacs.org
>>>>>> http://lists.gromacs.org/mailman/listinfo/gmx-developers
>>>>>> Please don't post (un)subscribe requests to the list. Use the
>>>>>> www interface or send it to gmx-developers-request at gromacs.org.
>>>>>> 
>>>>>> 
>>> -- 
>>> gmx-developers mailing list
>>> gmx-developers at gromacs.org
>>> http://lists.gromacs.org/mailman/listinfo/gmx-developers
>>> Please don't post (un)subscribe requests to the list. Use the www interface
>>> or send it to gmx-developers-request at gromacs.org.
>>> 
>> 
> -- 
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-developers-request at gromacs.org.