[gmx-developers] free energy calculations with restraints

Berk Hess hess at kth.se
Wed Aug 21 17:47:42 CEST 2013 

Hi,

I pushed up a patch to gerrit (master) which is complete afaik, but with 
a change as extensive as this, there is always a chance I overlooked 
something:
https://gerrit.gromacs.org/#/c/2566/

Please try and report back!

Cheers,

Berk

On 08/21/2013 10:08 AM, Mark Abraham wrote:
> Hi,
>
> I've been following with interest, and Berk and I discussed these
> issues yesterday. I'd been afraid there would be artifacts of
> "[moleculetype] is self-contained and exclusive" all through the code,
> but Berk seems pretty sure that's not the case. The kind of
> [intermolecule-interactions] suggested above sounds great to me, as a
> potential building block for all sorts of things.
>
> Mark
>
> On Tue, Aug 20, 2013 at 7:12 PM, David Mobley <dmobley at gmail.com> wrote:
>> Yes, I'm very enthusiastic about adding defined intermolecular restraints
>> using absolute atom indices. It seems like a straightforward, correct,
>> easy-to-understand way to deal with this simple problem.
>>
>> David
>>
>>
>>
>> On Tue, Aug 20, 2013 at 9:39 AM, Shirts, Michael (mrs5pt)
>> <mrs5pt at eservices.virginia.edu> wrote:
>>>
>>> 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.
>>> --
>>> 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.
>>
>>
>>
>> --
>> 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.

More information about the gromacs.org_gmx-developers mailing list