[gmx-users] centering a bilayer at a cartesian coordinate

[chris.neale at utoronto.ca]

Dear users:

I am trying to process a trajectory so that a group of molecules has  
their center of mass at a constant position in Cartesian space. I have  
not been able to figure out how to do this.

The reason that I would like this is that I have conducted umbrella  
sampling of a solute along the normal to a lipid bilayer and I would  
like to create a movie of a false-trajectory of the immersion profile.  
Such processing is also important for generating SDFs and for use with  
g_density.

I expect that the problems that I am having are due to volume  
fluctuations in the NPT ensemble.

The best idea that I have tried was (gromacs 4.5.3):

trjconv -center -box 10 10 15

where the original dimensions are much less that 10 10 15

But still, when I run g_traj on the processed .gro file I do not get  
similar values for the COM.

In more detail:

echo 0 | trjconv -s md1.tpr -f md1_50ps_md1.part0001.xtc -o tmp.xtc -e  
100000 -pbc mol
# select POPC for centering and System for output and use a .tpr with pbc=no
echo -e "13\n0\n" | trjconv -s empty.tpr -f tmp.xtc -box 10 10 15  
-center -pbc none -o tmp2.xtc
# select POPC for coordinate output
echo 13 | g_traj  -s empty.tpr -f tmp2.xtc -ox -com -nopbc

  99550	4.8424	5.00025	7.49519
  99600	4.79691	4.94304	7.39868
  99650	4.75789	4.85189	7.41581
  99700	4.74974	4.80423	7.48257
  99750	5.07858	5.04788	7.51654
  99800	4.99091	4.78976	7.47778
  99850	5.11406	4.84656	7.48769
  99900	5.09739	5.12395	7.47889
  99950	5.0411	4.88659	7.52737
  100000	4.91271	5.07101	7.50397

For Z, the value is 7.48037+/-0.072562 (MIN=7.3 and MAX=7.8)

If I look at the frames from the min and max z values, they do clearly  
differ in the placement of the bilayer along z.

###################

To simplify things, I redid this while only using the phosphorous  
atoms (1 per lipid) in the centering and coordinate output.

Here, Z=7.49323+/-0.0658727 (MIN=7.3 and MAX=7.7)

To verify this and since all the atoms have the same mass, I checked  
it with awk:

echo 0 | trjconv -s empty.tpr -f tmp2.xtc -dump 47550 -b 40000 -o look_low.gro
echo 0 | trjconv -s empty.tpr -f tmp2.xtc -dump 80750 -b 80000 -o  
look_high.gro
cat look_low.gro|grep P8|awk '{s+=$NF;n++} END{print s/n}'
7.39197
cat look_high.gro|grep P8|awk '{s+=$NF;n++} END{print s/n}'
7.62207

###################

Finally, I rechecked with a single atom and here it worked exactly as  
expected. The atoms is placed at 5 5 7.5 in every single frame.

###################

Then I checked using 2 phosphorous atoms from the same (or,  
separately, different) leaflets, and the result was the same:

$ tail coord.xvg
  99550	5	5	7.5
  99600	5	5.0005	7.5
  99650	5.0005	5.0005	7.5005
  99700	5	5.0005	7.5
  99750	5.0005	5	7.5005
  99800	5	5.0005	7.5005
  99850	5	5	7.5005
  99900	5.0005	5	7.5005
  99950	5.0005	5	7.4995
  100000	5	5	7.5


Then using 3 phosphorous atoms from (containing atoms from both  
leaflets), and now there is more deviation:

$ tail coord.xvg
  99550	5.48833	5.30233	6.972
  99600	5.48733	5.29233	6.97167
  99650	5.50033	5.29567	6.93467
  99700	5.52133	5.27433	6.93967
  99750	5.489	5.253	6.99133
  99800	5.481	5.319	6.91533
  99850	5.47867	5.309	6.827
  99900	5.48733	5.292	6.83133
  99950	5.52767	5.27	6.86867
  100000	5.518	5.319	6.89167

(MIN=6.7 and MAX=7.1)

I originally had problems with pbc during trjconv in which the -center  
option would place the bilayer at the peripheries of the unit cell in  
some frames (COM should still be at the center). This is what lead me  
to the route that I have used above (a .tpr that was generated with  
pbc=no and trjconv -pbc none -center). Nevertheless, I think that  
perhaps trjconv is using pbc information even in this case where ti  
should not.

##############################################

If anybody can see what I am doing wrong or has a method to center a  
bilayer at a cartesian coordinate, then I am very interested.

Thank you for your time,
Chris.