[gmx-users] keeping water (entirely) out of the bilayer core

[Christopher Neale]

Dear users:

I am interested in running simulations of lipid bilayers in which I keep all water molecules out of the bilayer core
(not just statistically, but absolutely). However, I have been unable to figure out how to do it. 
I'll list what I have tried in the hope that others have some ideas or even perhaps know how to do this 
with standard gromacs.

Everything that I have tried revolves around using the pull code, setting the entire lipid bilayer as the reference 
group, and having thousands of pulled groups -- one for each water molecule. Also, I modified the gromacs
source code in mdlib/pull.c (version 4.6.1) so that the force is only applied when the displacement is smaller than 
the desired distance and not when the displacement is larger than the specified distance (to keep water out but
then to otherwise let it go anywhere without bias):

static void do_pull_pot(int ePull,
                        t_pull *pull, t_pbc *pbc, double t, real lambda,
                        real *V, tensor vir, real *dVdl)
{
    int        g, j, m;
    dvec       dev;
    double     ndr, invdr;
    real       k, dkdl;
    t_pullgrp *pgrp;

    /* loop over the groups that are being pulled */
    *V    = 0;
    *dVdl = 0;
    for (g = 1; g < 1+pull->ngrp; g++)
    {
        pgrp = &pull->grp[g];
        get_pullgrp_distance(pull, pbc, g, t, pgrp->dr, dev);

/*  ########## HERE IS THE CODE CHANGE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        if(dev[0]>0.0){
                dev[0]=0.0;
        }
/* End code snippit */


The most relevant lines from the .mdp file were:

pull                     = umbrella
pull_geometry            = distance
pull_dim                 = N N Y
pull_start               = no
pull_ngroups             = 9000
pull_group0              = POPC

pull_group1              = r_1_&_OW
pull_init1               = 1.9
pull_k1                  = 500

pull_group2              = r_2_&_OW
pull_init2               = 1.9
pull_k2                  = 500

... etc...


The problem with this was that it crashed immediately with an error that my pulling distance was 
greater than 1/2 of the box vector.:

Fatal error:
Distance of pull group 165 (5.353939 nm) is larger than 0.49 times the box size (5.395960)

Pretty obvious in retrospect. If I could get this error to go away, then everything should be fine because I have modified the code so that forces are zero at large displacements. However, I am worried that if I simply modify grompp to bypass this error then I might get some type of strange and possibly silent error in mdrun. Any thoughts on this are really appreciated.

For my second try, I used pull_geometry = direction-periodic (see more details below), but that also didn't work 
because if I set pull_vec1 = 0 0 1, then everything gets pulled to the upper part of the box (and LINCS errors 
ensue), rather than simply pulling away from the bilayer.

Pcoupl                    = berendsen
pcoupltype           = semiisotropic
compressibility    = 4.5e-5 0
pull                     = umbrella
pull_geometry        = direction-periodic
pull_start               = no
pull_ngroups             = 9000
pull_group0              = POPC

pull_group1              = r_1_&_OW
pull_init1               = 1.9
pull_k1                  = 500
pull_vec1                = 0 0 1

pull_group2              = r_2_&_OW
pull_init2               = 1.9
pull_k2                  = 500
pull_vec2                = 0 0 1

... etc...


If anybody has an idea about how I could get this to work with standard gromacs or with a modified version, I would be very grateful to hear your thoughts. 

Thank you,
Chris.