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

[XAvier Periole]

You could define a repulsive potential that would apply to water molecules. I am not sure how it is called but it is available and described in the manual. 

On May 4, 2013, at 2:46, Christopher Neale <chris.neale at mail.utoronto.ca> wrote:

> 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.
> 
> -- 
> gmx-users mailing list    gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-users-request at gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists