using PLUMED to keep water molecules (entirely) out of the hydrophobic core of a lipid bilayer

[cand...@gmail.com]

Hello,

I am running gromacs simulations of a lipid bilayer in water. The bilayer extends along the Cartesian XY plane, with its normal along the Z-axis. I would like to add a potential that keeps water from penetrating within 2.0 nm of the bilayer center. The idea is to add a restraint on every water atom, acting only along the Z-axis, that forces water molecules to stay >=2.0 nm away from the bilayer center of mass alon the Z-axis.

I was keen to see that anything can be turned into a COM restraint in PLUMED 2.0, which is crucial for my usage (I want individual water molecules to be pushed away from the COM of the lipid bilayer).

If it makes things any clearer, below is what I have tried to do. I already know that it doesn't work, and I figure that DISTANCES simply don't have any label called ".z" (although I had hoped that they would because DISTANCE does appear to). Even in the absence of projecting the force along the Z-axis, I am still at a loss to figure out how one would use PLUMED to add a force to N water molecules based on their distance from a single COM).

COM ATOMS=22,23,24,25,26,27,28 LABEL=com1

DISTANCES GROUPA=com1 GROUPB=6678,6682,6686 LABEL=d1

UPPER_WALLS ARG=d1.z AT=2.0 KAPPA=200.0 EXP=2 EPS=1 OFFSET=0 LABEL=uwall

PRINT ARG=uwall.bias FILE=dat.wall STRIDE=10

Thank you for your help,

I'm an expert with gromacs, but very new to PLUMED.

PS: gromacs actually does have the facilities to do exactly what I want to do, but since water molecules can get very far away from the bilayer center, there seems to be some problem sorting out PBS in gromacs and it throws an error straight off, which is what started my trek to PLUMED.

Chris.

[cand...@gmail.com]

I was actually able to make some progress with the DISTANCE command (see below). However, regardless of whether I use UPPER_WALLS or LOWER_WALLS, the system compresses along the z-axis, making me think that what I am doing is not pulling along the negative of the vector of closest approach along z, but generally pulling along positive z. If I have to use this DISTANCE format, then is there a way that I can pull atoms in group 2 *away* from the COM of group 1, rather than in a specified direction?

INCLUDE FILE=popc2.group

COM ATOMS=popc LABEL=com1

DISTANCE ATOMS=com1,6678 COMPONENTS LABEL=d1

DISTANCE ATOMS=com1,6682 COMPONENTS LABEL=d2

DISTANCE ATOMS=com1,6686 COMPONENTS LABEL=d3

DISTANCE ATOMS=com1,6690 COMPONENTS LABEL=d4

DISTANCE ATOMS=com1,6694 COMPONENTS LABEL=d5

LOWER_WALLS ARG=d1.z,d2.z,d3.z,d4.z,d5.z AT=1.4,1.4,1.4,1.4,1.4 KAPPA=200,200,200,200,200 EXP=2,2,2,2,2 EPS=1,1,1,1,1 OFFSET=0,0,0,0,0 LABEL=uwall

PRINT ARG=uwall.bias FILE=wall STRIDE=10

Note that popc.group contains all atoms of the lipid bilayer and, in reality, I have ~10,000 water molecules and so ~10,000 DISTANCE calls, although I only put 5 here to describe what I am trying to do.

Thank you for your assistance,

Chris.

[Giovanni Bussi]

Hi Chris,

A possible mistake in the your input is that d1.z could be positive or negative (it is the component of the distance vector). So if you want its modulo (or squared modulo) you should take it explicitly (e.g. see directive COMBINE).

Anyway, I would do things in a different way, as using thousands of collective variables can be highly inefficient.

1. Choose a single collective variable that counts how many water molecules are within the region you want to avoid.

I would do it with COORDINATION between water oxygens and a selected set of atoms from the bilayer (e.g. all the lipid termini). E.g.

count: COORDINATION GROUPA=ow GROUPB=lipid R_0=2

Notice that this will not be a real counter, since if the lipid termini are closer than 2 nm you will count each water molecule more than once. But this should not be a problem for you. Also, instead of 2 nm you could use a smaller value for R_0.

You can easily make a simulation where water if free to enter just printing this variable and looking at the trajectory with vmd to see how well it counts the number of water molecules in the relevant region.

If this approach works, you can

* take advantage of the fact that COORDINATION is evaluated in parallel

* exploit advanced parameters for COORDINATION such as NLIST to accelerate the simulation

Also, there could be more efficient or more rigorous ways of counting waters in a region (see e.g. DENSITY and SUBCELL), see the manual. In any case, check everything visually to verify that you are really counting the water molecules inside the membrane.

2. Add a penalty on that counter.

The straightforward way is to restraint the counter to zero, but a direct linear penalty could be easier to tune:

restraint: RESTRAINT ARG=count AT=0 SLOPE=10 # this will add a penalty of 10 kj/mol for each water molecule entering the layer.

Ciao!

Giovanni

--
You received this message because you are subscribed to the Google Groups "PLUMED users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to plumed-users...@googlegroups.com.
To post to this group, send email to plumed...@googlegroups.com.
Visit this group at http://groups.google.com/group/plumed-users?hl=en.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[Gareth Tribello]

Hello Chris

Further to what Giovanni was suggesting here is how you would do this using DENSITY and SUBCELL.

DENSITY SPECIES=6678,6682,6686,6690,6694 LABEL=d1
SUBCELL ARG=d1 ZLOWER=0 ZUPPER=0.5 LABEL=s1
RESTRAINT ARG=s1 AT=0 SLOPE=10

This is assuming that your membrane is in the lower half of the simulation box. 

If you are willing to do a little bit of programming you can do this even more cleverly by writing some code to create a volume of interest.  For instance, lets say the z position of the center of mass is z_m.  What you want to do is essentially calculate all the atoms that have | z - z_m |<2 nm.   This is the same as integrating the density of the species, p(z), between z_m-2 and z_m+2.

I can help you do this if you would like.

Gareth

[Giovanni Bussi]

Hi Chris,

variables which looks like "counters" in plumed are typically computed using smooth switching function. The default is usually (1-(r/r0)**6)/(1-(r/r0)**12), but you can change it. This is done exactly to obtain differentiable variables for MD. Notice that this means that counting is not exact (you will find "4.5 water molecules"), but it's very useful for biased MD.

Ciao!

Giovanni

On Wed, May 8, 2013 at 9:37 PM, <cand...@gmail.com> wrote:

Dear Giovanni and Gareth:

Thank you for your assistance. I started by getting my initial (inefficient) approach working. The following plumed.dat file does exactly what I want, although you are right that it is inefficient when there are many thousands of independently restrained water molecules (especially with MPI).

INCLUDE FILE=popc2.group
COM ATOMS=popc LABEL=com1
DISTANCE ATOMS=com1,6678 COMPONENTS LABEL=d1

COMBINE ARG=d1.z POWERS=2 PERIODIC=NO LABEL=e1
DISTANCE ATOMS=com1,6682 COMPONENTS LABEL=d2
COMBINE ARG=d2.z POWERS=2 PERIODIC=NO LABEL=e2
DISTANCE ATOMS=com1,6686 COMPONENTS LABEL=d3
COMBINE ARG=d3.z POWERS=2 PERIODIC=NO LABEL=e3
DISTANCE ATOMS=com1,6690 COMPONENTS LABEL=d4
COMBINE ARG=d4.z POWERS=2 PERIODIC=NO LABEL=e4
DISTANCE ATOMS=com1,6694 COMPONENTS LABEL=d5
COMBINE ARG=d5.z POWERS=2 PERIODIC=NO LABEL=e5
LOWER_WALLS ARG=e1,e2,e3,e4,e5 AT=2.25,2.25,2.25,2.25,2.25 KAPPA=200,200,200,200,200 EXP=2,2,2,2,2 EPS=1,1,1,1,1 OFFSET=0,0,0,0,0 LABEL=lwall
PRINT ARG=lwall.bias FILE=wall STRIDE=10


Just to provide numerical confirmation of your suggestion that this approach is inefficient, on a single core I get: 3 ns/day without plumed, 3 ns/day with the 5-water plumed.dat file above, and 2 ns/day with the whole 9000-water plumed.dat. However, on 8 cores the inefficiency is magnified: I get 17 ns/day without plumed, 12 ns/day with the 5-water plumed.dat, and 1.2 ns/day with the whole 9000-water plumed.dat.

I am now looking into the approaches that you suggested: (i) COORDINATION or (ii) DENSITY/SUBCELL. In this respect, I have a question about how these methods will add forces to bias the dynamics, especially the COORDINATION approach. It seems to me that COORDINATION is not differentiable since there are either 0 or 1 or 2, etc. water molecules that are coordinated. I have essentially the same question about the DENSITY/SUBCELL approach. I can see how either of these energetic terms could be used to drive a bias in a MC simulation, but would appreciate it if you can help me to make sure that these methods can drive a bias in an MD simulation before I work on them.

Thank you again for all of your help.

[Gareth Tribello]

Hi Chris

Subsequently to emailing you we realized that there are some issues with the density and subcell option that I mentioned yesterday.  At the moment I wouldn't recommend starting to use them.  I am going to work on them today and will keep you posted on how fixing them is going. 

Gareth

[cand...@gmail.com]

Thank you Giovanni and Gareth. I'll look into COORDINATION for now.

BTW: I hacked the gromacs source code to get it to do what I want, but that is also very inefficient  (for more details, see http://lists.gromacs.org/pipermail/gmx-users/2013-May/081124.html ). Briefly, gromacs without this restraint is 26 ns/day, gromacs + DISTANCES in plumed gets 1.7 ns/day and gromacs hacked to do these restraints gets 4 ns/day. Hopefully COORDINATION will give the best performance.

Thanks again for all of your help,

Chris.

[cand...@gmail.com]

Thanks again for the help. I did get COORDINATION to work. The following plumed.dat file lead to a time of 19.6 ns/day (vs 26 ns/day without plumed). This efficiency is 75%. Not perfect, but much better than the DISTANCES plumed routine and also much better than thousands of restraints within gromacs (via my source code hack).

INCLUDE FILE=central_lipid_atoms.group

INCLUDE FILE=water.group

COORDINATION GROUPA=central_lipid_atoms GROUPB=water R_0=0.8 NLIST NL_CUTOFF=1.0 NL_STRIDE=10 LABEL=d1

UPPER_WALLS ARG=d1 AT=0 KAPPA=200 OFFSET=0 EXP=1 EPS=1 LABEL=uwall

PRINT ARG=uwall.bias FILE=wall STRIDE=10

I had to play with KAPPA and EXP to find a combination that didn;t crash when there was water in the bilayer at the start. (i.e., KAPPA=1000 EXP=1 crashed, as did KAPPA=10, EXP=2).

This is a workable solution for me and I am grateful for your assistance. If you do get the DENSITY/SUBCELL method working and think that it might be more efficient, then please let me know.

(Also, it would be nice if I could somehow force the UPPER_WALLS command to act only along the z dimension. Is that possible? I don;t see the option to use d1.z when d1 is a COORDINATION).

Thank you,

Chris.

[Giovanni Bussi]

Hi Chris,

a few comments on your input are below

On Fri, May 10, 2013 at 12:50 AM, <cand...@gmail.com> wrote:

Thanks again for the help. I did get COORDINATION to work. The following plumed.dat file lead to a time of 19.6 ns/day (vs 26 ns/day without plumed). This efficiency is 75%. Not perfect, but much better than the DISTANCES plumed routine and also much better than thousands of restraints within gromacs (via my source code hack).

INCLUDE FILE=central_lipid_atoms.group

INCLUDE FILE=water.group

COORDINATION GROUPA=central_lipid_atoms GROUPB=water R_0=0.8 NLIST NL_CUTOFF=1.0 NL_STRIDE=10 LABEL=d1

Be careful with NL_CUTOFF. It seems too small (should be at least twice R_0, since you want the switching function to be close to zero at that distance).

A simple way to check the neighbor list parameters is to perform a calculation with:

COORDINATION GROUPA=central_lipid_atoms GROUPB=water R_0=0.8 NLIST NL_CUTOFF=1.0 NL_STRIDE=10 LABEL=d1

COORDINATION GROUPA=central_lipid_atoms GROUPB=water R_0=0.8 LABEL=d1check

PRINT FILE=test ARG=d1,d1check

and verify in file "test" that with/without the NLIST you get the same result.

Also, there has been a few optimizations for neighbor lists in the v2.0 branch on our git repository, so this could run a bit faster if you download the latest update (or if you wait for the next beta snapshot).

 

UPPER_WALLS ARG=d1 AT=0 KAPPA=200 OFFSET=0 EXP=1 EPS=1 LABEL=uwall

Since coordination is always positive, this should be completely equivalent to

RESTRAINT ARG=d1 AT=0 SLOPE=200. I think the KAPPA (or the SLOPE) is very large, since it corresponds to a penalty of 200 kj/mol per water molecule in the bilayer. The fact that the simulation does not crash is a necessary condition, not a sufficient one.

 

PRINT ARG=uwall.bias FILE=wall STRIDE=10

I had to play with KAPPA and EXP to find a combination that didn;t crash when there was water in the bilayer at the start. (i.e., KAPPA=1000 EXP=1 crashed, as did KAPPA=10, EXP=2).

This is a workable solution for me and I am grateful for your assistance. If you do get the DENSITY/SUBCELL method working and think that it might be more efficient, then please let me know.

(Also, it would be nice if I could somehow force the UPPER_WALLS command to act only along the z dimension. Is that possible? I don;t see the option to use d1.z when d1 is a COORDINATION).

The point is not with the wall, but with the coordination itself. COORDINATION computes a pairwise sum of a switching function. What you want to do is change the switching function such that it only measures z displacement. This is currently not possible, and I am not sure it is a good idea, but it's trivial to code it: if you want I can give you some hint.

Giovanni

 

Thank you,

[cand...@gmail.com]

Thank you very much, Giovanni.

Thank you for the caution about NL_CUTOFF. I had actually assumed that the function would go to zero at R_0 and that the NL_CUTOFF was only greater than R_0 because NL_STRIDE was greater than 1 (This is what one does for neighbour lists in MD simulations, so I got confused). After receiving your advice, I have plotted the function of s and now understand it much better.

Thanks for the note about efficiency in neighbourlists in the git repo. I'll check that out once I feel comfortable with the system I am now working on.

I'll look into using RESTRAINT instead of UPPER_WALLS. I presume that it is not faster, just a clarity issue.

Indeed my KAPPA is very large and I appreciate your caution. I am not into production yet and am simply verifying that my approach can do what I want it to do. The next step will be to select R_0, KAPPA, etc. so that they are minimally perturbing.

Regarding COORDINATION.z, thank youf rot eh offer to help, but since I am new to plumed I think that I'll stick with what is available out of the box for now. I think that I can do everything that I want with what is available.