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

[Christopher Neale]

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.