<div dir="ltr"><div>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.<br>
<br></div>David<br><br><div><div></div></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Tue, Aug 20, 2013 at 9:39 AM, Shirts, Michael (mrs5pt) <span dir="ltr"><<a href="mailto:mrs5pt@eservices.virginia.edu" target="_blank">mrs5pt@eservices.virginia.edu</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br>
The other thing I would say is that the couple-mol code is a little<br>
difficult to work with----there's known quirks, like handling the<br>
intramolecular dispersion in a way that's slightly off from free energy off,<br>
as well as some new quirks that I'll try to post about later this<br>
week---that it would actually be cleaner to add clearly defined<br>
intermolecular restraints than add complexity to the coupling feature.<br>
<div class="im"><br>
Best,<br>
~~~~~~~~~~~~<br>
Michael Shirts<br>
Assistant Professor<br>
Department of Chemical Engineering<br>
University of Virginia<br>
<a href="mailto:michael.shirts@virginia.edu">michael.shirts@virginia.edu</a><br>
<a href="tel:%28434%29-243-1821" value="+14342431821">(434)-243-1821</a><br>
<br>
<br>
</div>> From: Floris Buelens <<a href="mailto:floris_buelens@yahoo.com">floris_buelens@yahoo.com</a>><br>
> Reply-To: Floris Buelens <<a href="mailto:floris_buelens@yahoo.com">floris_buelens@yahoo.com</a>>, Discussion list for<br>
> GROMACS development <<a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a>><br>
> Date: Tue, 20 Aug 2013 08:35:15 -0400<br>
> To: Berk Hess <<a href="mailto:hess@kth.se">hess@kth.se</a>>, Discussion list for GROMACS development<br>
> <<a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a>><br>
<div class="HOEnZb"><div class="h5">> Subject: Re: [gmx-developers] free energy calculations with restraints<br>
><br>
> Hi Berk,<br>
><br>
> I think we might have gone out of sync there - I sent the message below before<br>
> reading yours from this morning. From the misplaced context it maybe sounded<br>
> like I was disagreeing with your '[ intermolecular-interactions ]' suggestion,<br>
> which I'm definitely not - this would clearly be preferable and probably<br>
> useful to a wider audience.<br>
><br>
><br>
><br>
> ----- Original Message -----<br>
> From: Berk Hess <<a href="mailto:hess@kth.se">hess@kth.se</a>><br>
> To: Floris Buelens <<a href="mailto:floris_buelens@yahoo.com">floris_buelens@yahoo.com</a>>; Discussion list for GROMACS<br>
> development <<a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a>><br>
> Cc:<br>
> Sent: Tuesday, 20 August 2013, 12:04<br>
> Subject: Re: [gmx-developers] free energy calculations with restraints<br>
><br>
> Hi,<br>
><br>
> We could consider this.<br>
> I don't know how you did this now, but this would need an extra mdp<br>
> option, something like couple-group, to avoid one option with a moltype<br>
> or index group as an input.<br>
> Then we would also need to decide what to do with a group which is part<br>
> of one molecule of a certain moleculetype with multiple copies. Disallow<br>
> this? Or duplicate the moltype and only modify one copy (as is done now<br>
> for qmmm)?<br>
><br>
> Cheers,<br>
><br>
> Berk<br>
><br>
> On 08/20/2013 09:24 AM, Floris Buelens wrote:<br>
>> Hi,<br>
>><br>
>> I would also still argue for augmenting the current couple-moltype<br>
>> functionality with this modification. I can see this is a niche requirement,<br>
>> which of course does speak against inclusion. But in my view the replacement<br>
>> code is no less elegant than the current implementation (modification of five<br>
>> self-contained functions at the pre-processing stage, maybe 50-odd extra<br>
>> lines), we can maintain full backward compatibility, and a single simple<br>
>> check will remove any potential for accidental misuse.<br>
>><br>
>><br>
>><br>
>> ________________________________<br>
>> From: David Mobley <<a href="mailto:dmobley@uci.edu">dmobley@uci.edu</a>><br>
>> To: Discussion list for GROMACS development <<a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a>><br>
>> Cc: Floris Buelens <<a href="mailto:floris_buelens@yahoo.com">floris_buelens@yahoo.com</a>>; Michael R. Shirts<br>
>> <<a href="mailto:michael.shirts@virginia.edu">michael.shirts@virginia.edu</a>><br>
>> Sent: Monday, 19 August 2013, 20:02<br>
>> Subject: Re: [gmx-developers] free energy calculations with restraints<br>
>><br>
>><br>
>><br>
>> Hi, Berk,<br>
>><br>
>> I agree that this solution would be the "proper" one and would be nice.<br>
>> However, it's been a long time coming (this is an issue I've been raising for<br>
>> approximately 5 years, and there's no progress; normally no one even<br>
>> answers). Currently, for the restraints I need to use, I *have* to merge the<br>
>> ligand into my protein topology, inconvenient or not. Doing this also means<br>
>> there are some features I simply can't use (couple-moltype for example).<br>
>><br>
>> It strikes me that Floris's solution is a good interim one, in that it at<br>
>> least solves most of the typical use cases we're facing, and has the<br>
>> advantage of being *already implemented*.<br>
>><br>
>> Is there any way we could get this into the main GROMACS, at least until<br>
>> someone gets time to implementing the "proper" solution (which could be years<br>
>> from now)?<br>
>><br>
>> Thanks!<br>
>> David<br>
>><br>
>><br>
>><br>
>><br>
>><br>
>> On Mon, Aug 19, 2013 at 1:10 AM, Berk Hess <<a href="mailto:hess@kth.se">hess@kth.se</a>> wrote:<br>
>><br>
>> Hi,<br>
>>> I fully understood this point, maybe my reply wasn't well structured.<br>
>>> I was trying to argue that the proper solution would be to implement<br>
>>> inter-molecular restraints, instead of introducing an option which could me<br>
>>> misused. Also merging a ligand topology into your protein topology is very<br>
>>> inconvenient. This problem with the proper solution is, of course, that<br>
>>> someone will need to implement inter-molecular restraints/potentials.<br>
>>><br>
>>> One issue used to be that we corrected molecules for PBC before calculating<br>
>>> bonded interactions. But in 4.6 this is usually not faster and not used any<br>
>>> more. That makes it easier to treat intra- and inter-molecular interactions<br>
>>> the same way at the mdrun level.<br>
>>><br>
>>> Cheers,<br>
>>><br>
>>> Berk<br>
>>><br>
>>><br>
>>> On 08/19/2013 09:55 AM, Floris Buelens wrote:<br>
>>><br>
>>> Hi Berk,<br>
>>>> I think some clarification is needed. What we're talking about is only<br>
>>>> relevant in the context of non-bonded interactions, typically between a<br>
>>>> small molecule and a protein. To access the binding affinity of a ligand<br>
>>>> through alchemical methods, it's useful to switch off only the interactions<br>
>>>> of the ligand with its environment while maintaining intra-ligand<br>
>>>> potentials - this is what's provided by the 'couple-moltype' code path.<br>
>>>><br>
>>>> The limitation that we're trying to circumvent comes from the fact that as<br>
>>>> you scale down interactions with the environment, the ligand is no longer<br>
>>>> held in the binding site. To counteract this, it's necessary to perturb in<br>
>>>> restraining potentials as the intramolecular nonbonded interactions are<br>
>>>> perturbed out. The most practical method makes use of a single distance<br>
>>>> restraint, two angle and three dihedral restraints, whose effect can later<br>
>>>> be accounted for analytically.<br>
>>>><br>
>>>> The current code only allows decoupling of a whole logical 'molecule' as<br>
>>>> specified in the topology file. This precludes the use of regular (fully<br>
>>>> perturbation-aware) bonded interactions. As far as I'm aware, the pull code<br>
>>>> doesn't provide what's required. Inter-molecular bonded interactions would<br>
>>>> be a great general solution but presumably won't show up any time soon.<br>
>>>><br>
>>>> The workaround here is instead to represent two physical molecules (e.g.<br>
>>>> protein and ligand) as a single logical molecule (in the topology file), so<br>
>>>> regular bonded potentials can be applied. Current Gromacs doesn't allow the<br>
>>>> decoupling ('couple-moltype') code to be used in this scenario. My<br>
>>>> modification allows you to target the ligand by its residue name and get<br>
>>>> this functionality back.<br>
>>>><br>
>>>> Decoupling a single residue of a multi-residue chain is indeed probably not<br>
>>>> correct in this framework (I did call it 'weird' in my last message :-) ).<br>
>>>> An extra check would block users from trying this.<br>
>>>><br>
>>>> I hope that clarifies the problem we're trying to solve. I agree this will<br>
>>>> be useful a relatively small number of users, but on the other hand it's a<br>
>>>> very unobtrusive, self-contained modification which can maintain full<br>
>>>> backwards compatibility.<br>
>>>><br>
>>>> thanks,<br>
>>>><br>
>>>> Floris<br>
>>>><br>
>>>><br>
>>>> ----- Original Message -----<br>
>>>> From: Berk Hess <<a href="mailto:hess@kth.se">hess@kth.se</a>><br>
>>>> To: Floris Buelens <<a href="mailto:floris_buelens@yahoo.com">floris_buelens@yahoo.com</a>>; Discussion list for GROMACS<br>
>>>> development <<a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a>><br>
>>>> Cc:<br>
>>>> Sent: Monday, 19 August 2013, 9:03<br>
>>>> Subject: Re: [gmx-developers] free energy calculations with restraints<br>
>>>><br>
>>>> Hi,<br>
>>>><br>
>>>> There is the rotational pull code, but that might not provide the exact<br>
>>>> functionality you want.<br>
>>>> Already for some time we have been discussing inter-molecular<br>
>>>> interactions, by specifying molecule type, molecule index and atom<br>
>>>> index, but there are no concrete plans for implementing this yet.<br>
>>>><br>
>>>> An option for decoupling part of a molecule indeed sound useful. But in<br>
>>>> practice you always need to replace that part by something else, at<br>
>>>> least a hydrogen, and modify some potentials of connecting atoms, so I<br>
>>>> don't know how generally useful such an option is.<br>
>>>><br>
>>>> Cheers,<br>
>>>><br>
>>>> Berk<br>
>>>><br>
>>>> On 08/19/2013 07:56 AM, Floris Buelens wrote:<br>
>>>><br>
>>>> I suppose that would make sense. The changes are fairly minor (more or less<br>
>>>> just the function convert_moltype_couple and the four functions it calls)<br>
>>>> and mainly consist of a bit of extra juggling with nonbonded exclusions<br>
>>>> during preprocessing. The current functionality (decouple a whole molecule)<br>
>>>> could be handled as a special case of the new code. A check for bonds to<br>
>>>> the rest of the structure could stop people from trying weird stuff like<br>
>>>> decoupling residues from a chain.<br>
>>>>><br>
>>>>><br>
>>>>><br>
>>>>> ________________________________<br>
>>>>> From: David Mobley <<a href="mailto:dmobley@gmail.com">dmobley@gmail.com</a>><br>
>>>>> To: Floris Buelens <<a href="mailto:floris_buelens@yahoo.com">floris_buelens@yahoo.com</a>><br>
>>>>> Cc: Discussion list for GROMACS development <<a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a>><br>
>>>>> Sent: Friday, 16 August 2013, 18:30<br>
>>>>> Subject: Re: [gmx-developers] free energy calculations with restraints<br>
>>>>><br>
>>>>><br>
>>>>><br>
>>>>> This does sound useful, though it would be more useful if this could be<br>
>>>>> implemented into the main GROMACS rather than a separate code (since<br>
>>>>> otherwise it will go away as GROMACS is further developed).<br>
>>>>><br>
>>>>> Would this be a possibility going forward, devels?<br>
>>>>><br>
>>>>> Thanks!<br>
>>>>><br>
>>>>><br>
>>>>><br>
>>>>><br>
>>>>> On Fri, Aug 16, 2013 at 3:07 AM, Floris Buelens <<a href="mailto:floris_buelens@yahoo.com">floris_buelens@yahoo.com</a>><br>
>>>>> wrote:<br>
>>>>><br>
>>>>> Hi David,<br>
>>>>><br>
>>>>> I have the same requirement as you, but I've gone about it slightly<br>
>>>>> differently. I've hacked the couple-moltype code path to allow decoupling<br>
>>>>> of a specific residue (identified by name) instead of a molecule. This<br>
>>>>> allows the residue of interest to be part of another molecule block so you<br>
>>>>> can set up distance, angle and dihedral restraints using regular bonded<br>
>>>>> interactions with full perturbation support.<br>
>>>>>> If this is useful to you or to anyone else, let me know and I'll be happy<br>
>>>>>> to share.<br>
>>>>>><br>
>>>>>> best,<br>
>>>>>><br>
>>>>>> Floris<br>
>>>>>><br>
>>>>>> ________________________________<br>
>>>>>> From: David Mobley <<a href="mailto:dmobley@gmail.com">dmobley@gmail.com</a>><br>
>>>>>><br>
>>>>>> To: Discussion list for GROMACS development <<a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a>><br>
>>>>>> Sent: Friday, 14 June 2013, 21:47<br>
>>>>>> Subject: [gmx-developers] free energy calculations with restraints<br>
>>>>>><br>
>>>>>><br>
>>>>>><br>
>>>>>><br>
>>>>>> Hi, devs,<br>
>>>>>><br>
>>>>>><br>
>>>>>> I'm writing with an issue relating to the interplay of new free energy<br>
>>>>>> features with restraints.<br>
>>>>>><br>
>>>>>> I'm very much appreciating some of the new free energy features in<br>
>>>>>> gromacs, such as the 'couple-moltype' option which provides a way to set<br>
>>>>>> up decoupling or annihilation of a specific molecule via free energy<br>
>>>>>> calculations without having to edit the topology file directly. This is<br>
>>>>>> especially great when it comes to decoupling -- charge decoupling was not<br>
>>>>>> previously possible via topology file editing, and vdW decoupling took<br>
>>>>>> substantial manipulation of the topology file.<br>
>>>>>><br>
>>>>>> However, for binding free energies, my work employs orientational<br>
>>>>>> restraints between the protein and ligand. I need to be able to impose<br>
>>>>>> both dihedral and angle restraints on angles between the protein and<br>
>>>>>> ligand. Currently, I do this using angle-restraints and<br>
>>>>>> dihedral-restraints. This requires that both the protein and ligand be<br>
>>>>>> within the same logical 'molecule', which (unfortunately) means that I<br>
>>>>>> can't make use of the new free energy features above, since<br>
>>>>>> couple-moltype has to apply to a whole molecule, not just part of a<br>
>>>>>> molecule.<br>
>>>>>><br>
>>>>>> So, my I see two possible solutions to the problem, and hence have these<br>
>>>>>> questions:<br>
>>>>>> 1) Can dihedral and angle restraints be applied via the pull code? If<br>
>>>>>> not, are there any plans to add that?<br>
>>>>>> 2) Alternatively, what about modifying the restraints code so it uses (or<br>
>>>>>> at least optionally allows) absolute atom numbering, rather than<br>
>>>>>> numbering within a specific molecule, thus allowing restraints<br>
>>>>>> (angle-restraints and dihedral-restraints) to be applied between<br>
>>>>>> 'molecules'?<br>
>>>>>><br>
>>>>>> Thanks!<br>
>>>>>> David<br>
>>>>>><br>
>>>>>><br>
>>>>>><br>
>>>>>> --<br>
>>>>>> David Mobley<br>
>>>>>> <a href="mailto:dmobley@gmail.com">dmobley@gmail.com</a><br>
>>>>>> <a href="tel:949-385-2436" value="+19493852436">949-385-2436</a><br>
>>>>>><br>
>>>>>><br>
>>>>>> --<br>
>>>>>> gmx-developers mailing list<br>
>>>>>> <a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a><br>
>>>>>> <a href="http://lists.gromacs.org/mailman/listinfo/gmx-developers" target="_blank">http://lists.gromacs.org/mailman/listinfo/gmx-developers</a><br>
>>>>>> Please don't post (un)subscribe requests to the list. Use the<br>
>>>>>> www interface or send it to <a href="mailto:gmx-developers-request@gromacs.org">gmx-developers-request@gromacs.org</a>.<br>
>>>>>><br>
>>>>>><br>
>>> --<br>
>>> gmx-developers mailing list<br>
>>> <a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a><br>
>>> <a href="http://lists.gromacs.org/mailman/listinfo/gmx-developers" target="_blank">http://lists.gromacs.org/mailman/listinfo/gmx-developers</a><br>
>>> Please don't post (un)subscribe requests to the list. Use the www interface<br>
>>> or send it to <a href="mailto:gmx-developers-request@gromacs.org">gmx-developers-request@gromacs.org</a>.<br>
>>><br>
>><br>
> --<br>
> gmx-developers mailing list<br>
> <a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a><br>
> <a href="http://lists.gromacs.org/mailman/listinfo/gmx-developers" target="_blank">http://lists.gromacs.org/mailman/listinfo/gmx-developers</a><br>
> Please don't post (un)subscribe requests to the list. Use the<br>
> www interface or send it to <a href="mailto:gmx-developers-request@gromacs.org">gmx-developers-request@gromacs.org</a>.<br>
<br>
--<br>
gmx-developers mailing list<br>
<a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-developers" target="_blank">http://lists.gromacs.org/mailman/listinfo/gmx-developers</a><br>
Please don't post (un)subscribe requests to the list. Use the<br>
www interface or send it to <a href="mailto:gmx-developers-request@gromacs.org">gmx-developers-request@gromacs.org</a>.<br>
</div></div></blockquote></div><br><br clear="all"><br>-- <br>David Mobley<br><a href="mailto:dmobley@gmail.com" target="_blank">dmobley@gmail.com</a><br>949-385-2436<br>
</div>