[gmx-developers] free energy calculations with restraints

Wed Aug 21 17:58:35 CEST 2013

This changes the tpr format, we would like to avoid that in a patches 
branch.
You can use master branch now and hopefully during most of the 
development process.
It would actually be good if more people (or even some) use master for 
production, so we get some testing of the code. Note that all code still 
needs to pass our tests before it gets merged into master.

Cheers,

Berk

On 08/21/2013 05:49 PM, Shirts, Michael (mrs5pt) wrote:
> One question:
>
> If it is in master, it may be a 6-12 months before anyone can really use it.
> Would it make sense to port into the 4.6 repository?
>
> ~~~~~~~~~~~~
> 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.