[gmx-developers] Understanding the g_rdf

[Christoph Freudenberger]

Hi there,

As I posted in the users list a couple of days before, I wanted
to write an extension to the g_rdf code to calculate the spatial
distrubtion function.
I've been working myself through the src during the last days.
I understand the main part of the code, but I hope
somebody might find the time to explain a few details:

   fprintf(stderr,"\nHow many groups do you want to calculate the RDF of?\n");
   do {
     scanf("%d",&ng);
   } while (ng < 1);
   snew(grpname,ng+1);
   snew(isize,ng+1);
   snew(index,ng+1);
   fprintf(stderr,"\nSelect a reference group and %d group%s\n",
           ng,ng==1?"":"s");
   if (fnTPS)
     get_index(&top.atoms,fnNDX,ng+1,isize,index,grpname);
   else
     rd_index(fnNDX,ng+1,isize,index,grpname);

I suppose this part reads in the index groups.
fnNDX is the filename, ng+1 is the number of groups to be read...

from rdgroups.h i gathered that
isize is of dim [ng+1] and stores the number of atoms in each group.
grpname is of equal size and stores the group-names (What are they
used for, besides the title string in the xvg-file?)
index should then be a 2D array of size [ng+1][isize[i]].
right so far?

   if (bCM) {
     /* move index[0] to index_cm and make 'dummy' index[0] */
     isize_cm=isize[0];
     snew(index_cm,isize_cm);
     for(i=0; i<isize[0]; i++)
       index_cm[i]=index[0][i];
     isize[0]=1;
     index[0][0]=natoms;
     srenew(index[0],isize[0]);
     /* make space for center of mass */
     srenew(x,natoms+1);
   }
This is the setup, for the com calculation. The data for the first
group is copied, isize[0] is set to 1 because the first group will
be represented by one single "atom"...
but why is index[0][0]=natoms?

I don't quite understand how the com is accually calculated:
     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]);
       }
This first loop seems to sum up the position vectors of all
atoms in the com group into the vector xcom.
       /* store it in the first 'group' */
       for(j=0; (j<DIM); j++)
         x[index[0][0]][j] = xcom[j] / isize_cm;
     }
Then xcom/isize_cm is stored... but why isn't there any weighting
by the masses of the atoms?

best regards
-- 
Christoph Freudenberger
Abt. Organische Chemie I AK Siehl - Uni Ulm -Tel: ++49-731-502-2785