[gmx-developers] free energy calculations with restraints

Shirts, Michael (mrs5pt) mrs5pt at eservices.virginia.edu
Wed Aug 21 18:10:00 CEST 2013 

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