[gmx-users] Re, density

[Christopher Neale]

Dear Ali:

I am not sure what ecenter is, it looks like it might be the position to which the user wants to set the center?
In any event, it doesn't mater for this modification.

I did not write the center_x() function, it is a standard part of gmx_trjconv.c What I did was to copy center_x() to 
center_x_com() and then modify center_x_com() such that it computes the center of mass and not the position of
(max-min)/2.

I then modified the call to center_x() to actually be a call to center_x_com(), where it requires the additional 
argument of "atoms->atom" at the end of the list.

What I posted here: http://lists.gromacs.org/pipermail/gmx-users/2011-June/062239.html is actually the non-
modified function call from the original gmx_trjconv.c in gromacs 4.5.3 (posted to show that there is a problem).

Here is a patch that you can apply to gromacs version 4.5.3 (but possibly not to 4.5.4 or 4.5.5 because they made 
slight other modifications to gmx_trjconv.c so I am not sure if this patch will work). Still, it is just 2 blocks of changes
so once you understand them you can make the changes by hand.


--- ../../../../gromacs-4.5.3b/source/src/tools/gmx_trjconv.c	2010-09-29 07:35:06.000000000 -0400
+++ ./gmx_trjconv.c	2011-08-02 13:13:00.000000000 -0400
@@ -385,6 +385,33 @@
     }
 }
 
+static void center_x_com(int ecenter,rvec x[],matrix box,int n,int nc,atom_id ci[],t_atom atom[])
+{
+    int i,m,ai;
+    rvec com,box_center,dx;
+    real mtot;
+
+    mtot=0.0;
+    com[0]=com[1]=com[2]=0.0;
+    if (nc > 0) {
+        for(i=0; i<nc; i++) {
+            ai=ci[i];
+            for(m=0; m<DIM; m++) {
+                com[m]+=atom[ai].m*x[ai][m];
+            }
+            mtot+=atom[ai].m;
+        }
+        calc_box_center(ecenter,box,box_center);
+        for(m=0; m<DIM; m++)
+            dx[m] = box_center[m]-com[m]/mtot;
+
+        for(i=0; i<n; i++)
+            rvec_inc(x[i],dx);
+
+//fprintf(stderr,"want %f %f %f but have %f %f %f\n",box_center[0],box_center[1],box_center[2],
+    }
+}
+
 static void mk_filenm(char *base,const char *ext,int ndigit,int file_nr,
                       char out_file[])
 {
@@ -1291,8 +1318,10 @@
                                         rvec_inc(fr.x[i],x_shift);
                             }
 
-                            if (bCenter)
-                                center_x(ecenter,fr.x,fr.box,natoms,ncent,cindex);
+                            if (bCenter){
+                                //center_x(ecenter,fr.x,fr.box,natoms,ncent,cindex);
+                                center_x_com(ecenter,fr.x,fr.box,natoms,ncent,cindex,atoms->atom);
+                            }
                         }
 
                         if (bPBCcomAtom) {





######## PATCH IS OVER

In case the patch is mangled by email formatting, you can try to just add the lines yourself. Just add the 
center_x_com() function and then use the call to center_x_com() and not center_x() in the "if(bCenter)" sataement.
Note that there is the additional argument "atoms->atom" in the call to center_x_com().


Chris.

-- original message --

I find out where you change the trjconv, I have a question

"ecenter" in this code, Where do you use it? and What is the ecenter?


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);
    }
}


-- 
Sincerely

Ali Alizadeh