[gmx-developers] free energy calculations with restraints

Tue Aug 20 16:52:49 CEST 2013

I wonder if the cleanest way to do this is to add this capability though the
pull code.  The pull code already cleanly handles multimolecule information,
and then the top file can stay with molecule specific indexing.

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: Tue, 20 Aug 2013 03:17:49 -0400
> To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
> Subject: Re: Fwd: [gmx-developers] free energy calculations with restraints
> 
> Hi,
> 
> I fully understand your sentiments.
> But I think this solution, although fine for personal use, it to ad-hoc to
> include in Gromacs. Currently the couple-moltype decouples a whole molecule,
> which is a very clear operation (although already somewhat complicated to
> understand due to all the different coupling mechanisms). I certainly don't
> want to replace this by a index group option. Adding an index group option is
> also messy.
> 
> I looked at the domain decomposition code and there is nothing molecule
> specific to how the bonded interactions are assigned. So on the mdrun level
> there is nothing that prevents the working of inter-molecular interactions.
> The only thing that needs to be done is adding something like an [
> intermolecular-interactions ] section to the topology file processing and the
> topology data structure. Having the current processing code go through this is
> a matter of adding a few lines. I estimate that implementing this would take
> me an hour or so, if we stick to simple, absolute system (not molecule) atom
> indices in this section. I can do this within a week, promise :) The only
> disadvantage is that the ligand will not automatically end up in the same
> periodic image as the protein.
> 
> I will discuss this with our local developers to get some feedback.
> 
> Cheers,
> 
> Berk
> 
> On 08/19/2013 08:31 PM, David Mobley wrote:
> 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<mailto: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<mailto:hess at kth.se>>
> To: Floris Buelens
> <floris_buelens at yahoo.com<mailto:floris_buelens at yahoo.com>>; Discussion list
> for GROMACS development
> <gmx-developers at gromacs.org<mailto: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<mailto:dmobley at gmail.com>>
> To: Floris Buelens <floris_buelens at yahoo.com<mailto:floris_buelens at yahoo.com>>
> Cc: Discussion list for GROMACS development
> <gmx-developers at gromacs.org<mailto: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<mailto: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<mailto:dmobley at gmail.com>>
> 
> To: Discussion list for GROMACS development
> <gmx-developers at gromacs.org<mailto: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<mailto:dmobley at gmail.com>
> 949-385-2436<tel:949-385-2436>
> 
> 
> --
> gmx-developers mailing list
> gmx-developers at gromacs.org<mailto: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<mailto:gmx-developers-request at gromacs.org>.
> 
> 
> --
> gmx-developers mailing list
> gmx-developers at gromacs.org<mailto: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<mailto:gmx-developers-request at gromacs.org>.
> 
> 
> 
> --
> David Mobley
> Assistant Professor
> Department of Pharmaceutical Sciences
> Department of Chemistry
> 3134B Natural Sciences I
> University of California, Irvine
> Irvine, CA 92697
> dmobley at uci.edu<mailto:dmobley at uci.edu>
> work (949) 824-6383<tel:%28949%29%20824-6383>
> cell (949) 385-2436<tel:%28949%29%20385-2436>
> 
> 
> 
> 
> --
> David Mobley
> dmobley at gmail.com<mailto:dmobley at gmail.com>
> 949-385-2436
> 
> 
> 