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

[Dallas Warren]

Chris, not entirely sure if would work, but what about set of particles located at the center of the bilayer (might require anchoring/freezing) which only interacts (via repulsion only) with water, and setting the interaction such that it manages to exclude it from the bilayer region you require.

Catch ya,

Dr. Dallas Warren
Drug Discovery Biology
Monash Institute of Pharmaceutical Sciences, Monash University
381 Royal Parade, Parkville VIC 3052
dallas.warren at monash.edu
+61 3 9903 9304
---------------------------------
When the only tool you own is a hammer, every problem begins to resemble a nail. 

> -----Original Message-----
> From: gmx-users-bounces at gromacs.org [mailto:gmx-users-
> bounces at gromacs.org] On Behalf Of Christopher Neale
> Sent: Saturday, 4 May 2013 10:46 AM
> To: gmx-users at gromacs.org
> Subject: [gmx-users] keeping water (entirely) out of the bilayer core
> 
> 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