Dear Jianguo, XAvier, and Dallas:
Thank you for your great suggestions. I am making progress, but still have not
found an efficient way
to do this. I don;t have any particular questions in this post, but wanted to
provide an update and
also to see if anybody has any additional ideas based on what I have tried.
Jianguo: I have successfully interfaced plumed with gromacs and was able to do
what I wanted.
My initial (naive) approach is inefficient (speed reduced from 26 ns/day to 1.7
ns.day) although I
am getting some help from the plumed developers to find a more efficient
approach.
XAvier: I was unable to find a repulsive potential that I could apply to water
molecules, excepting
simply increasing the LJ repulsion. I suspect that this approach is not going
to work for me since the
lipid-water interface is rugged and lipid tails sometimes approach the surface
and, therefore,
a LJ repulsion approach might end up creating large vaccuums at the interfacial
region.
Dallas and Jianguo: This might work. I would need to find a way to anchor
virtual atoms at the bilayer center.
Unfortunately, other parts of my approach make it undesirable to have an r^12
(LJ) or exponential (Buckingham)
repulsion. Specifically, I'd like to be able to not only keep water molecules
out of the bilayer core
but also force them out if they are already inside, and these rapidly rising
potentials will likely lead
to expolding systems if water molecules are already deeply burried in the
bilayer core. However, there might
be a way forward with the free energy code using some type of soft core
potential on the LJ via the
free energy code. I'm going to try this and will report back.
Also, I did get gromacs to do this with the following 2 code changes in
src/mdlib/pull.c:
(i) in static void do_pull_pot(), immediately after calling
get_pullgrp_distance(), I
added the following code
if(g>1&&dev[0]>0.0){
dev[0]=0.0;
}
When used with ePull == epullUMBRELLA and pull->eGeom == epullgDIST, this
ensures that the pull force only
contributes when the displacement is smaller than the reference distance, but
is not applied when
the displacement is larger than the reference distance.
(ii) I removed the following code from get_pullgrps_dr():
if (max_dist2 >= 0 && dr2 > 0.98*0.98*max_dist2)
{
gmx_fatal(FARGS, "Distance of pull group %d (%f nm) is larger than 0.49
times the box size (%f)", g, sqrt(dr2), sqrt(max_dist2));
}
Normally, I suppose that this code is there so that atoms are not being pulled
in oscillating directions when
they are 1/2 box dimension away from the reference position (although I am not
exactly sure why that would
be an error). Nevertheless, this error is not necessary in my case because I
have also removed any force when atoms are farther away than the reference
distance (in modificaion i, above).
After making these modification, I then added about 10,000 separate pull
groups, one for each water oxygen
atom, each of which keeps that water molecule out of the bilayer core along Z.
The only problem with this approach is that it is really inefficient: with 16
threads over 8 cores, I get
26 ns/day without the 10,000 pull groups and only 4 ns/day with the pull
groups, probably because so much
information is being passed around between threads for the water molecules that
are really far away from
the bilayer center.
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
