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@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-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
