[gmx-developers] free energy calculations with restraints

Shirts, Michael (mrs5pt) mrs5pt at eservices.virginia.edu
Fri Aug 30 23:51:48 CEST 2013 

Berk, just following up on the requests below to make it easier to test the
code:

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


> From: "Michael R. Shirts" <mrs5pt at eservices.virginia.edu>
> Reply-To: <michael.shirts at virginia.edu>
> Date: Wed, 21 Aug 2013 10:09:55 -0400
> To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
> Subject: Re: [gmx-developers] free energy calculations with restraints
> 
> Hi, Berk-
> 
> * Would it be possible to post (perhaps to redmine?) some example files of
> calling usage?
> 
> * Perhaps you can provide a bit more of the logic for how this is
> implemented? I was expecting a new interaction_function type, but this seems
> to be done somewhat differently, and I'm not quite sure how it works -- I
> suspect other people will be confused as well.
> 
> Best,
> ~~~~~~~~~~~~
> Michael Shirts
> Assistant Professor
> Department of Chemical Engineering
> University of Virginia
> michael.shirts at virginia.edu
> (434)-243-1821
> 
> 
>> From: Berk Hess <hess at kth.se>
>> Reply-To: Discussion list for GROMACS development
>> <gmx-developers at gromacs.org>
>> Date: Wed, 21 Aug 2013 17:47:42 +0200
>> To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
>> Subject: Re: [gmx-developers] free energy calculations with restraints
>> 
>> 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.
>> 
>> -- 
>> 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