[gmx-developers] free energy calculations with restraints

Mark Abraham mark.j.abraham at gmail.com
Wed Aug 21 10:08:06 CEST 2013 

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