[gmx-developers] comment on g_rdf -com

[Jochen Hub]

Dear developers,

I just stepped over the piece of code in g_rdf that calculates the 
center of mass in case of -com usage:

     if (bCM) {
         /* calculate centre of mass */
         clear_rvec(xcom);
         for(i=0; (i < isize_cm); i++) {
             ix = index_cm[i];
             rvec_inc(xcom,x[ix]);
         }
         /* store it in the first 'group' */
         for(j=0; (j<DIM); j++)
             x[index[0][0]][j] = xcom[j] / isize_cm;
     }

This does, however, only give the center of mass if the masses of all 
atoms are the same. There is no mass information in this peace of code. 
Something like this would help:

     if (bCM) {
         /* calculate centre of mass */
         clear_rvec(xcom);
         M=0;
         for(i=0; (i < isize_cm); i++) {
             ix  = index_cm[i];
             m = top.atoms.atom[ix].m;
             M+=m;
             svmul(m, x[ix], vec_tmp);
             rvec_inc(xcom,vec_tmp);
         }
         /* store it in the first 'group' */
         for(j=0; (j<DIM); j++)
             x[index[0][0]][j] = xcom[j] / M;
     }

Probably this hardly makes any difference in the results (the COM is 
probably close the the center of atoms). So not a serious "bug"...

Cheers,
Jochen


-- 
************************************************
Jochen Hub
Max Planck Institute for Biophysical Chemistry
Computational biomolecular dynamics group
Am Fassberg 11
D-37077 Goettingen, Germany
Email: jhub[at]gwdg.de
************************************************