I solved it. It was my own error to assume that trjconv -center
centers the COM. Actually, it centers the value of (min+max)/2 in each
dimension.
I can modify trjconv locally to do what I want.
Thank you,
Chris.
static void center_x(int ecenter,rvec x[],matrix box,int n,int
nc,atom_id ci[])
{
int i,m,ai;
rvec cmin,cmax,box_center,dx;
if (nc > 0) {
copy_rvec(x[ci[0]],cmin);
copy_rvec(x[ci[0]],cmax);
for(i=0; i<nc; i++) {
ai=ci[i];
for(m=0; m<DIM; m++) {
if (x[ai][m] < cmin[m])
cmin[m]=x[ai][m];
else if (x[ai][m] > cmax[m])
cmax[m]=x[ai][m];
}
}
calc_box_center(ecenter,box,box_center);
for(m=0; m<DIM; m++)
dx[m] = box_center[m]-(cmin[m]+cmax[m])*0.5;
for(i=0; i<n; i++)
rvec_inc(x[i],dx);
}
}
Quoting chris.ne...@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.
--
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
