[gmx-developers] decoupling in free energy calculations for binding
chris.neale at utoronto.ca
chris.neale at utoronto.ca
Sat Jun 25 18:31:22 CEST 2011
Dear David:
The flat-bottom potential is not available in the pull code. I have
modified the source code to hard-code a flat-bottomed potential. You
can patch src/mdlib/pull.c in gromacs 4.5.4 like this:
--- pull.c 2011-03-15 08:44:34.000000000 -0400
+++ pull_flat.c 2011-06-25 12:18:40.439242000 -0400
@@ -61,6 +61,8 @@
#include "gmx_ga2la.h"
#include "copyrite.h"
+#define PFB 0.5
+
static void pull_print_x_grp(FILE *out,gmx_bool bRef,ivec
dim,t_pullgrp *pgrp)
{
int m;
@@ -365,6 +367,7 @@
int m;
dvec ref;
double drs,inpr;
+ double mydist;
pgrp = &pull->grp[g];
@@ -381,6 +384,16 @@
{
ref[0] = pgrp->init[0] + pgrp->rate*t;
}
+
+ mydist=sqrt(dr[0]*dr[0]+dr[1]*dr[1]+dr[2]*dr[2]);
+ if(mydist<=PFB){
+ dr[0]=0.0;dr[1]=0.0;dr[2]=0.0;
+ }else{
+ //reduce it's length by PFB
+ dr[0]*=(mydist-PFB)/mydist;
+ dr[1]*=(mydist-PFB)/mydist;
+ dr[2]*=(mydist-PFB)/mydist;
+ }
switch (pull->eGeom)
{
#################################################################
That modification was tested to work for pull-code options:
pull = umbrella
pull_geometry = position
But I think it should probably be general. Please do test it as it is
my extrapolation from a modification that I made to 4.0.5 a while ago.
That region of the source has not changed, so it seems like it should
be ok in 4.5.4, but I did not actually run any tests. The patch is for
4.5.4 though. Also note that the PFB flat bottom distance variable is
hard-coded as a defined value.
#########
As for the angles and the dihedrals, the pull code can not do this for
you. Further, one can only have a single reference group with the
pull-code so having done the flat-bottomed distance restraint with the
pull-code we have in effect used that up.
Dihedral restraints you could do like this:
http://www.gromacs.org/Documentation/How-tos/Dihedral_Restraints
although I haven't used those since gromacs 3.3.3.
Angles you might modify the topology.
But the dihedral and angle would run into the problem that they must
be in the same molecule and thus can not be decoupled as you desire.
I wonder, is it possible to decouple only a part of a "molecule" by
defining the B-state yourself with some scripting?
Chris.
-- original message --
Quoting David Mobley <dmobley at gmail.com>:
> Hi,
>
> The specific thing we're trying to achieve here (in this case) is add a flat
> bottom distance restraint between an atom in the ligand and a virtual site
> in the center of a protein binding site. Not sure if that can be achieved
> with the pull code?
>
> We are also (separately) interested in having distance, angle, and dihedral
> restraints between reference atoms in the ligand and reference atoms in the
> protein. These restraints need to be lambda dependent, and ideally we would
> like to be able to use them with decoupling (again, currently not possible).
> Is the pull code flexible enough to be able to do angle and dihedral
> restraints somehow? (There is a fairly specific set of restraints we need
> for technical reasons -- one distance, two angles, three torsions. Other
> forms won't get us what we need...).
>
> Thanks.
>
>
> On Thu, Jun 23, 2011 at 4:23 AM, Berk Hess <hess at cbr.su.se> wrote:
>
>> **
>> Hi,
>>
>> You can manually set up the decoupling, but that's very tedious and error
>> prone.
>> I think the best procedure would be to restrain using the pull code,
>> but I don't know if that's currently flexible enough to cover most cases.
>> If not, we should think about extend the functionality of the pull code.
>>
>> I though the decoupled state is described somewhere in the manual.
>> It is a non-periodic state with pure Coulomb and LJ interactions without
>> cut-offs
>> (unless you use couple-intramol).
>>
>> Berk
>>
>>
>> On 06/23/2011 04:46 AM, David Mobley wrote:
>>
>> Hi,
>>
>> We're trying to do absolute binding free energy calculations using
>> Michael Shirts' latest free energy code (to be included in 4.6 if I
>> understand correctly; right now it's a branch of 4.5). These require using
>> a restraint between the protein and the ligand, which currently we're doing
>> using virtual sites. I am interested in also doing these calculations using
>> "decoupling", wherein internal interactions of the perturbed molecule (in
>> this case the ligand) are retained. This involves something like the
>> following in the mdp file:
>>
>> couple-moltype = MOL
>> couple-lambda0 = vdw-coul
>> couple-lambda1 = none
>> couple-intramol = no
>>
>> if, for example, the ligand is a molecule named 'MOL'.
>>
>> My question is this: Is there any way to get this to work when the ligand
>> and the protein are part of the same "molecule"? Specifically, to have
>> restraints betweeen the protein and the ligand (such as using a virtual
>> site) which are NECESSARY for absolute binding free energy calculations, I
>> must have the protein and ligand as part of the same molecule (unless
>> there's a workaround I'm not aware of). But to get decoupling to work using
>> the above I seem to need to have the protein and ligand as separate
>> molecules, suggesting they are incompatible. Is there a workaround I'm not
>> aware of?
>>
>> Also, on a related note -- when decoupling is done, what is the end state
>> for the decoupled object? Is it the decoupled object in gas phase in a
>> periodic system (interacting with copies of itself), or in a nonperiodic
>> system?
>>
>> Thanks!
>>
>> --
>> David Mobley
>> dmobley at gmail.com
>> 504-383-3662
>>
>>
>>
>>
>> --
>> 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
> Assistant Professor
> Department of Chemistry
> University of New Orleans
> dlmobley at uno.edu
> ph. 504-383-3662
> fax 504-280-6860
>
>
>
>
>
>
> --
> David Mobley
> dmobley at gmail.com
> 504-383-3662
>
More information about the gromacs.org_gmx-developers
mailing list