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

[chris.neale at utoronto.ca]

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.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.
>
>
>