src/tools/gmx_rdf.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <math.h>
  41 #include <ctype.h>
  42 #include "string2.h"
  43 #include "sysstuff.h"
  44 #include "typedefs.h"
  45 #include "macros.h"
  46 #include "vec.h"
  47 #include "pbc.h"
  48 #include "xvgr.h"
  49 #include "copyrite.h"
  50 #include "futil.h"
  51 #include "statutil.h"
  52 #include "tpxio.h"
  53 #include "physics.h"
  54 #include "index.h"
  55 #include "smalloc.h"
  56 #include "calcgrid.h"
  57 #include "nrnb.h"
  58 #include "coulomb.h"
  59 #include "gstat.h"
  60 #include "matio.h"
  61 #include "gmx_ana.h"
  62 #include "names.h"
  63 #include "sfactor.h"
  64
  65
  66 static void check_box_c(matrix box)
  67 {
  68   if (fabs(box[ZZ][XX]) > GMX_REAL_EPS*box[ZZ][ZZ] ||
  69       fabs(box[ZZ][YY]) > GMX_REAL_EPS*box[ZZ][ZZ])
  70     gmx_fatal(FARGS,
  71               "The last box vector is not parallel to the z-axis: %f %f %f",
  72               box[ZZ][XX],box[ZZ][YY],box[ZZ][ZZ]);
  73 }
  74
  75 static void calc_comg(int is,int *coi,int *index,gmx_bool bMass,t_atom *atom,
  76                       rvec *x,rvec *x_comg)
  77 {
  78   int  c,i,d;
  79   rvec xc;
  80   real mtot,m;
  81
  82   if (bMass && atom==NULL)
  83     gmx_fatal(FARGS,"No masses available while mass weighting was requested");
  84
  85   for(c=0; c<is; c++) {
  86     clear_rvec(xc);
  87     mtot = 0;
  88     for(i=coi[c]; i<coi[c+1]; i++) {
  89       if (bMass) {
  90         m = atom[index[i]].m;
  91         for(d=0; d<DIM; d++)
  92           xc[d] += m*x[index[i]][d];
  93         mtot += m;
  94       } else {
  95         rvec_inc(xc,x[index[i]]);
  96         mtot += 1.0;
  97       }
  98     }
  99     svmul(1/mtot,xc,x_comg[c]);
 100   }
 101 }
 102
 103 static void split_group(int isize,int *index,char *grpname,
 104                         t_topology *top,char type,
 105                         int *is_out,int **coi_out)
 106 {
 107   t_block *mols=NULL;
 108   t_atom  *atom=NULL;
 109   int     is,*coi;
 110   int     cur,mol,res,i,a,i1;
 111
 112   /* Split up the group in molecules or residues */
 113   switch (type) {
 114   case 'm':
 115     mols = &top->mols;
 116     break;
 117   case 'r':
 118     atom = top->atoms.atom;
 119     break;
 120   default:
 121     gmx_fatal(FARGS,"Unknown rdf option '%s'",type);
 122   }
 123   snew(coi,isize+1);
 124   is = 0;
 125   cur = -1;
 126   mol = 0;
 127   for(i=0; i<isize; i++) {
 128     a = index[i];
 129     if (type == 'm') {
 130       /* Check if the molecule number has changed */
 131       i1 = mols->index[mol+1];
 132       while(a >= i1) {
 133         mol++;
 134         i1 = mols->index[mol+1];
 135       }
 136       if (mol != cur) {
 137         coi[is++] = i;
 138         cur = mol;
 139       }
 140     } else if (type == 'r') {
 141       /* Check if the residue index has changed */
 142       res = atom[a].resind;
 143       if (res != cur) {
 144         coi[is++] = i;
 145         cur = res;
 146       }
 147     }
 148   }
 149   coi[is] = i;
 150   srenew(coi,is+1);
 151   printf("Group '%s' of %d atoms consists of %d %s\n",
 152          grpname,isize,is,
 153          (type=='m' ? "molecules" : "residues"));
 154
 155   *is_out  = is;
 156   *coi_out = coi;
 157 }
 158
 159 static void do_rdf(const char *fnNDX,const char *fnTPS,const char *fnTRX,
 160                    const char *fnRDF,const char *fnCNRDF, const char *fnHQ,
 161                    gmx_bool bCM,const char *close,
 162                    const char **rdft,gmx_bool bXY,gmx_bool bPBC,gmx_bool bNormalize,
 163                    real cutoff,real binwidth,real fade,int ng,
 164                    const output_env_t oenv)
 165 {
 166   FILE       *fp;
 167   t_trxstatus *status;
 168   char       outf1[STRLEN],outf2[STRLEN];
 169   char       title[STRLEN],gtitle[STRLEN],refgt[30];
 170   int        g,natoms,i,ii,j,k,nbin,j0,j1,n,nframes;
 171   int        **count;
 172   char       **grpname;
 173   int        *isize,isize_cm=0,nrdf=0,max_i,isize0,isize_g;
 174   atom_id    **index,*index_cm=NULL;
 175 #if (defined SIZEOF_LONG_LONG_INT) && (SIZEOF_LONG_LONG_INT >= 8)
 176   long long int *sum;
 177 #else
 178   double     *sum;
 179 #endif
 180   real       t,rmax2,cut2,r,r2,r2ii,invhbinw,normfac;
 181   real       segvol,spherevol,prev_spherevol,**rdf;
 182   rvec       *x,dx,*x0=NULL,*x_i1,xi;
 183   real       *inv_segvol,invvol,invvol_sum,rho;
 184   gmx_bool       bClose,*bExcl,bTop,bNonSelfExcl;
 185   matrix     box,box_pbc;
 186   int        **npairs;
 187   atom_id    ix,jx,***pairs;
 188   t_topology *top=NULL;
 189   int        ePBC=-1,ePBCrdf=-1;
 190   t_block    *mols=NULL;
 191   t_blocka   *excl;
 192   t_atom     *atom=NULL;
 193   t_pbc      pbc;
 194   gmx_rmpbc_t  gpbc=NULL;
 195   int        *is=NULL,**coi=NULL,cur,mol,i1,res,a;
 196
 197   excl=NULL;
 198
 199   bClose = (close[0] != 'n');
 200
 201   if (fnTPS) {
 202     snew(top,1);
 203     bTop=read_tps_conf(fnTPS,title,top,&ePBC,&x,NULL,box,TRUE);
 204     if (bTop && !(bCM || bClose))
 205       /* get exclusions from topology */
 206       excl = &(top->excls);
 207   }
 208   snew(grpname,ng+1);
 209   snew(isize,ng+1);
 210   snew(index,ng+1);
 211   fprintf(stderr,"\nSelect a reference group and %d group%s\n",
 212           ng,ng==1?"":"s");
 213   if (fnTPS) {
 214     get_index(&(top->atoms),fnNDX,ng+1,isize,index,grpname);
 215     atom = top->atoms.atom;
 216   } else {
 217     rd_index(fnNDX,ng+1,isize,index,grpname);
 218   }
 219
 220   if (bCM || bClose) {
 221     snew(is,1);
 222     snew(coi,1);
 223     if (bClose) {
 224       split_group(isize[0],index[0],grpname[0],top,close[0],&is[0],&coi[0]);
 225     }
 226   }
 227   if (rdft[0][0] != 'a') {
 228     /* Split up all the groups in molecules or residues */
 229     srenew(is,ng+1);
 230     srenew(coi,ng+1);
 231     for(g=((bCM || bClose) ? 1 : 0); g<ng+1; g++) {
 232       split_group(isize[g],index[g],grpname[g],top,rdft[0][0],&is[g],&coi[g]);
 233     }
 234   }
 235
 236   if (bCM) {
 237     is[0] = 1;
 238     snew(coi[0],is[0]+1);
 239     coi[0][0] = 0;
 240     coi[0][1] = isize[0];
 241     isize0 = is[0];
 242     snew(x0,isize0);
 243   } else if (bClose || rdft[0][0] != 'a') {
 244     isize0 = is[0];
 245     snew(x0,isize0);
 246   } else {
 247     isize0 = isize[0];
 248   }
 249
 250   natoms=read_first_x(oenv,&status,fnTRX,&t,&x,box);
 251   if ( !natoms )
 252     gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
 253   if (fnTPS)
 254     /* check with topology */
 255     if ( natoms > top->atoms.nr )
 256       gmx_fatal(FARGS,"Trajectory (%d atoms) does not match topology (%d atoms)",
 257                   natoms,top->atoms.nr);
 258   /* check with index groups */
 259   for (i=0; i<ng+1; i++)
 260     for (j=0; j<isize[i]; j++)
 261       if ( index[i][j] >= natoms )
 262         gmx_fatal(FARGS,"Atom index (%d) in index group %s (%d atoms) larger "
 263                   "than number of atoms in trajectory (%d atoms)",
 264                   index[i][j],grpname[i],isize[i],natoms);
 265
 266   /* initialize some handy things */
 267   if (ePBC == -1) {
 268     ePBC = guess_ePBC(box);
 269   }
 270   copy_mat(box,box_pbc);
 271   if (bXY) {
 272     check_box_c(box);
 273     switch (ePBC) {
 274     case epbcXYZ:
 275     case epbcXY:   ePBCrdf = epbcXY;   break;
 276     case epbcNONE: ePBCrdf = epbcNONE; break;
 277     default:
 278       gmx_fatal(FARGS,"xy-rdf's are not supported for pbc type'%s'",
 279                 EPBC(ePBC));
 280       break;
 281     }
 282     /* Make sure the z-height does not influence the cut-off */
 283     box_pbc[ZZ][ZZ] = 2*max(box[XX][XX],box[YY][YY]);
 284   } else {
 285     ePBCrdf = ePBC;
 286   }
 287   if (bPBC)
 288     rmax2   = 0.99*0.99*max_cutoff2(bXY ? epbcXY : epbcXYZ,box_pbc);
 289   else
 290     rmax2   = sqr(3*max(box[XX][XX],max(box[YY][YY],box[ZZ][ZZ])));
 291   if (debug)
 292     fprintf(debug,"rmax2 = %g\n",rmax2);
 293
 294   /* We use the double amount of bins, so we can correctly
 295    * write the rdf and rdf_cn output at i*binwidth values.
 296    */
 297   nbin     = (int)(sqrt(rmax2) * 2 / binwidth);
 298   invhbinw = 2.0 / binwidth;
 299   cut2   = sqr(cutoff);
 300
 301   snew(count,ng);
 302   snew(pairs,ng);
 303   snew(npairs,ng);
 304
 305   snew(bExcl,natoms);
 306   max_i = 0;
 307   for(g=0; g<ng; g++) {
 308     if (isize[g+1] > max_i)
 309       max_i = isize[g+1];
 310
 311     /* this is THE array */
 312     snew(count[g],nbin+1);
 313
 314     /* make pairlist array for groups and exclusions */
 315     snew(pairs[g],isize[0]);
 316     snew(npairs[g],isize[0]);
 317     for(i=0; i<isize[0]; i++) {
 318       /* We can only have exclusions with atomic rdfs */
 319       if (!(bCM || bClose || rdft[0][0] != 'a')) {
 320         ix = index[0][i];
 321         for(j=0; j < natoms; j++)
 322           bExcl[j] = FALSE;
 323         /* exclusions? */
 324         if (excl)
 325           for( j = excl->index[ix]; j < excl->index[ix+1]; j++)
 326             bExcl[excl->a[j]]=TRUE;
 327         k = 0;
 328         snew(pairs[g][i], isize[g+1]);
 329         bNonSelfExcl = FALSE;
 330         for(j=0; j<isize[g+1]; j++) {
 331           jx = index[g+1][j];
 332           if (!bExcl[jx])
 333             pairs[g][i][k++]=jx;
 334           else if (ix != jx)
 335             /* Check if we have exclusions other than self exclusions */
 336             bNonSelfExcl = TRUE;
 337         }
 338         if (bNonSelfExcl) {
 339           npairs[g][i]=k;
 340           srenew(pairs[g][i],npairs[g][i]);
 341         } else {
 342           /* Save a LOT of memory and some cpu cycles */
 343           npairs[g][i]=-1;
 344           sfree(pairs[g][i]);
 345         }
 346       } else {
 347         npairs[g][i]=-1;
 348       }
 349     }
 350   }
 351   sfree(bExcl);
 352
 353   snew(x_i1,max_i);
 354   nframes = 0;
 355   invvol_sum = 0;
 356   if (bPBC && (NULL != top))
 357     gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
 358
 359   do {
 360     /* Must init pbc every step because of pressure coupling */
 361     copy_mat(box,box_pbc);
 362     if (bPBC) {
 363       if (top != NULL)
 364         gmx_rmpbc(gpbc,natoms,box,x);
 365       if (bXY) {
 366         check_box_c(box);
 367         clear_rvec(box_pbc[ZZ]);
 368       }
 369       set_pbc(&pbc,ePBCrdf,box_pbc);
 370
 371       if (bXY)
 372         /* Set z-size to 1 so we get the surface iso the volume */
 373         box_pbc[ZZ][ZZ] = 1;
 374     }
 375     invvol = 1/det(box_pbc);
 376     invvol_sum += invvol;
 377
 378     if (bCM) {
 379       /* Calculate center of mass of the whole group */
 380       calc_comg(is[0],coi[0],index[0],TRUE           ,atom,x,x0);
 381     } else if (!bClose && rdft[0][0] != 'a') {
 382       calc_comg(is[0],coi[0],index[0],rdft[0][6]=='m',atom,x,x0);
 383     }
 384
 385     for(g=0; g<ng; g++) {
 386       if (rdft[0][0] == 'a') {
 387         /* Copy the indexed coordinates to a continuous array */
 388         for(i=0; i<isize[g+1]; i++)
 389           copy_rvec(x[index[g+1][i]],x_i1[i]);
 390       } else {
 391         /* Calculate the COMs/COGs and store in x_i1 */
 392         calc_comg(is[g+1],coi[g+1],index[g+1],rdft[0][6]=='m',atom,x,x_i1);
 393       }
 394
 395       for(i=0; i<isize0; i++) {
 396         if (bClose) {
 397           /* Special loop, since we need to determine the minimum distance
 398            * over all selected atoms in the reference molecule/residue.
 399            */
 400           if (rdft[0][0] == 'a')
 401             isize_g = isize[g+1];
 402           else
 403             isize_g = is[g+1];
 404           for(j=0; j<isize_g; j++) {
 405             r2 = 1e30;
 406             /* Loop over the selected atoms in the reference molecule */
 407             for(ii=coi[0][i]; ii<coi[0][i+1]; ii++) {
 408               if (bPBC)
 409                 pbc_dx(&pbc,x[index[0][ii]],x_i1[j],dx);
 410               else
 411                 rvec_sub(x[index[0][ii]],x_i1[j],dx);
 412               if (bXY)
 413                 r2ii = dx[XX]*dx[XX] + dx[YY]*dx[YY];
 414               else
 415                 r2ii = iprod(dx,dx);
 416               if (r2ii < r2)
 417                 r2 = r2ii;
 418             }
 419             if (r2>cut2 && r2<=rmax2)
 420               count[g][(int)(sqrt(r2)*invhbinw)]++;
 421           }
 422         } else {
 423           /* Real rdf between points in space */
 424           if (bCM || rdft[0][0] != 'a') {
 425             copy_rvec(x0[i],xi);
 426           } else {
 427             copy_rvec(x[index[0][i]],xi);
 428           }
 429           if (rdft[0][0] == 'a' && npairs[g][i] >= 0) {
 430             /* Expensive loop, because of indexing */
 431             for(j=0; j<npairs[g][i]; j++) {
 432               jx=pairs[g][i][j];
 433               if (bPBC)
 434                 pbc_dx(&pbc,xi,x[jx],dx);
 435               else
 436                 rvec_sub(xi,x[jx],dx);
 437
 438               if (bXY)
 439                 r2 = dx[XX]*dx[XX] + dx[YY]*dx[YY];
 440               else
 441                 r2=iprod(dx,dx);
 442               if (r2>cut2 && r2<=rmax2)
 443                 count[g][(int)(sqrt(r2)*invhbinw)]++;
 444             }
 445           } else {
 446             /* Cheaper loop, no exclusions */
 447             if (rdft[0][0] == 'a')
 448               isize_g = isize[g+1];
 449             else
 450               isize_g = is[g+1];
 451             for(j=0; j<isize_g; j++) {
 452               if (bPBC)
 453                 pbc_dx(&pbc,xi,x_i1[j],dx);
 454               else
 455                 rvec_sub(xi,x_i1[j],dx);
 456               if (bXY)
 457                 r2 = dx[XX]*dx[XX] + dx[YY]*dx[YY];
 458               else
 459                 r2=iprod(dx,dx);
 460               if (r2>cut2 && r2<=rmax2)
 461                 count[g][(int)(sqrt(r2)*invhbinw)]++;
 462             }
 463           }
 464         }
 465       }
 466     }
 467     nframes++;
 468   } while (read_next_x(oenv,status,&t,natoms,x,box));
 469   fprintf(stderr,"\n");
 470
 471   if (bPBC && (NULL != top))
 472     gmx_rmpbc_done(gpbc);
 473
 474   close_trj(status);
 475
 476   sfree(x);
 477
 478   /* Average volume */
 479   invvol = invvol_sum/nframes;
 480
 481   /* Calculate volume of sphere segments or length of circle segments */
 482   snew(inv_segvol,(nbin+1)/2);
 483   prev_spherevol=0;
 484   for(i=0; (i<(nbin+1)/2); i++) {
 485     r = (i + 0.5)*binwidth;
 486     if (bXY) {
 487       spherevol=M_PI*r*r;
 488     } else {
 489       spherevol=(4.0/3.0)*M_PI*r*r*r;
 490     }
 491     segvol=spherevol-prev_spherevol;
 492     inv_segvol[i]=1.0/segvol;
 493     prev_spherevol=spherevol;
 494   }
 495
 496   snew(rdf,ng);
 497   for(g=0; g<ng; g++) {
 498     /* We have to normalize by dividing by the number of frames */
 499     if (rdft[0][0] == 'a')
 500       normfac = 1.0/(nframes*invvol*isize0*isize[g+1]);
 501     else
 502       normfac = 1.0/(nframes*invvol*isize0*is[g+1]);
 503
 504     /* Do the normalization */
 505     nrdf = max((nbin+1)/2,1+2*fade/binwidth);
 506     snew(rdf[g],nrdf);
 507     for(i=0; i<(nbin+1)/2; i++) {
 508       r = i*binwidth;
 509       if (i == 0)
 510         j = count[g][0];
 511       else
 512         j = count[g][i*2-1] + count[g][i*2];
 513       if ((fade > 0) && (r >= fade))
 514         rdf[g][i] = 1 + (j*inv_segvol[i]*normfac-1)*exp(-16*sqr(r/fade-1));
 515       else {
 516         if (bNormalize)
 517           rdf[g][i] = j*inv_segvol[i]*normfac;
 518         else
 519           rdf[g][i] = j/(binwidth*isize0*nframes);
 520       }
 521     }
 522     for( ; (i<nrdf); i++)
 523       rdf[g][i] = 1.0;
 524   }
 525
 526   if (rdft[0][0] == 'a') {
 527     sprintf(gtitle,"Radial distribution");
 528   } else {
 529     sprintf(gtitle,"Radial distribution of %s %s",
 530             rdft[0][0]=='m' ? "molecule" : "residue",
 531             rdft[0][6]=='m' ? "COM" : "COG");
 532   }
 533   fp=xvgropen(fnRDF,gtitle,"r","",oenv);
 534   if (bCM) {
 535     sprintf(refgt," %s","COM");
 536   } else if (bClose) {
 537     sprintf(refgt," closest atom in %s.",close);
 538   } else {
 539     sprintf(refgt,"%s","");
 540   }
 541   if (ng==1) {
 542     if (output_env_get_print_xvgr_codes(oenv))
 543       fprintf(fp,"@ subtitle \"%s%s - %s\"\n",grpname[0],refgt,grpname[1]);
 544   }
 545   else {
 546     if (output_env_get_print_xvgr_codes(oenv))
 547         fprintf(fp,"@ subtitle \"reference %s%s\"\n",grpname[0],refgt);
 548     xvgr_legend(fp,ng,(const char**)(grpname+1),oenv);
 549   }
 550   for(i=0; (i<nrdf); i++) {
 551     fprintf(fp,"%10g",i*binwidth);
 552     for(g=0; g<ng; g++)
 553       fprintf(fp," %10g",rdf[g][i]);
 554     fprintf(fp,"\n");
 555   }
 556   ffclose(fp);
 557
 558   do_view(oenv,fnRDF,NULL);
 559
 560   /* h(Q) function: fourier transform of rdf */
 561   if (fnHQ) {
 562     int nhq = 401;
 563     real *hq,*integrand,Q;
 564
 565     /* Get a better number density later! */
 566     rho = isize[1]*invvol;
 567     snew(hq,nhq);
 568     snew(integrand,nrdf);
 569     for(i=0; (i<nhq); i++) {
 570       Q = i*0.5;
 571       integrand[0] = 0;
 572       for(j=1; (j<nrdf); j++) {
 573         r = j*binwidth;
 574         integrand[j]  = (Q == 0) ? 1.0 : sin(Q*r)/(Q*r);
 575         integrand[j] *= 4.0*M_PI*rho*r*r*(rdf[0][j]-1.0);
 576       }
 577       hq[i] = print_and_integrate(debug,nrdf,binwidth,integrand,NULL,0);
 578     }
 579     fp=xvgropen(fnHQ,"h(Q)","Q(/nm)","h(Q)",oenv);
 580     for(i=0; (i<nhq); i++)
 581       fprintf(fp,"%10g %10g\n",i*0.5,hq[i]);
 582     ffclose(fp);
 583     do_view(oenv,fnHQ,NULL);
 584     sfree(hq);
 585     sfree(integrand);
 586   }
 587
 588   if (fnCNRDF) {
 589     normfac = 1.0/(isize0*nframes);
 590     fp=xvgropen(fnCNRDF,"Cumulative Number RDF","r","number",oenv);
 591     if (ng==1) {
 592       if (output_env_get_print_xvgr_codes(oenv))
 593         fprintf(fp,"@ subtitle \"%s-%s\"\n",grpname[0],grpname[1]);
 594     }
 595     else {
 596       if (output_env_get_print_xvgr_codes(oenv))
 597         fprintf(fp,"@ subtitle \"reference %s\"\n",grpname[0]);
 598       xvgr_legend(fp,ng,(const char**)(grpname+1),oenv);
 599     }
 600     snew(sum,ng);
 601     for(i=0; (i<=nbin/2); i++) {
 602       fprintf(fp,"%10g",i*binwidth);
 603       for(g=0; g<ng; g++) {
 604         fprintf(fp," %10g",(real)((double)sum[g]*normfac));
 605         if (i*2+1 < nbin)
 606           sum[g] += count[g][i*2] + count[g][i*2+1];
 607       }
 608       fprintf(fp,"\n");
 609     }
 610     ffclose(fp);
 611     sfree(sum);
 612
 613     do_view(oenv,fnCNRDF,NULL);
 614   }
 615
 616   for(g=0; g<ng; g++)
 617     sfree(rdf[g]);
 618   sfree(rdf);
 619 }
 620
 621
 622 int gmx_rdf(int argc,char *argv[])
 623 {
 624   const char *desc[] = {
 625     "The structure of liquids can be studied by either neutron or X-ray",
 626     "scattering. The most common way to describe liquid structure is by a",
 627     "radial distribution function. However, this is not easy to obtain from",
 628     "a scattering experiment.[PAR]",
 629     "[TT]g_rdf[tt] calculates radial distribution functions in different ways.",
 630     "The normal method is around a (set of) particle(s), the other methods",
 631     "are around the center of mass of a set of particles ([TT]-com[tt])",
 632     "or to the closest particle in a set ([TT]-surf[tt]).",
 633     "With all methods rdf's can also be calculated around axes parallel",
 634     "to the z-axis with option [TT]-xy[tt].",
 635     "With option [TT]-surf[tt] normalization can not be used.[PAR]",
 636     "The option [TT]-rdf[tt] sets the type of rdf to be computed.",
 637     "Default is for atoms or particles, but one can also select center",
 638     "of mass or geometry of molecules or residues. In all cases only",
 639     "the atoms in the index groups are taken into account.",
 640     "For molecules and/or the center of mass option a run input file",
 641     "is required.",
 642     "Other weighting than COM or COG can currently only be achieved",
 643     "by providing a run input file with different masses.",
 644     "Options [TT]-com[tt] and [TT]-surf[tt] also work in conjunction",
 645     "with [TT]-rdf[tt].[PAR]",
 646     "If a run input file is supplied ([TT]-s[tt]) and [TT]-rdf[tt] is set",
 647     "to [TT]atom[tt], exclusions defined",
 648     "in that file are taken into account when calculating the rdf.",
 649     "The option [TT]-cut[tt] is meant as an alternative way to avoid",
 650     "intramolecular peaks in the rdf plot.",
 651     "It is however better to supply a run input file with a higher number of",
 652     "exclusions. For eg. benzene a topology with nrexcl set to 5",
 653     "would eliminate all intramolecular contributions to the rdf.",
 654     "Note that all atoms in the selected groups are used, also the ones",
 655     "that don't have Lennard-Jones interactions.[PAR]",
 656     "Option [TT]-cn[tt] produces the cumulative number rdf,",
 657     "i.e. the average number of particles within a distance r.[PAR]",
 658     "To bridge the gap between theory and experiment structure factors can",
 659     "be computed (option [TT]-sq[tt]). The algorithm uses FFT, the grid",
 660     "spacing of which is determined by option [TT]-grid[tt]."
 661   };
 662   static gmx_bool bCM=FALSE,bXY=FALSE,bPBC=TRUE,bNormalize=TRUE;
 663   static real cutoff=0,binwidth=0.002,grid=0.05,fade=0.0,lambda=0.1,distance=10;
 664   static int  npixel=256,nlevel=20,ngroups=1;
 665   static real start_q=0.0, end_q=60.0, energy=12.0;
 666
 667   static const char *closet[]= { NULL, "no", "mol", "res", NULL };
 668   static const char *rdft[]={ NULL, "atom", "mol_com", "mol_cog", "res_com", "res_cog", NULL };
 669
 670   t_pargs pa[] = {
 671     { "-bin",      FALSE, etREAL, {&binwidth},
 672       "Binwidth (nm)" },
 673     { "-com",      FALSE, etBOOL, {&bCM},
 674       "RDF with respect to the center of mass of first group" },
 675     { "-surf",     FALSE, etENUM, {closet},
 676       "RDF with respect to the surface of the first group" },
 677     { "-rdf",   FALSE, etENUM, {rdft},
 678       "RDF type" },
 679     { "-pbc",      FALSE, etBOOL, {&bPBC},
 680       "Use periodic boundary conditions for computing distances. Without PBC the maximum range will be three times the largest box edge." },
 681     { "-norm",     FALSE, etBOOL, {&bNormalize},
 682       "Normalize for volume and density" },
 683     { "-xy",       FALSE, etBOOL, {&bXY},
 684       "Use only the x and y components of the distance" },
 685     { "-cut",      FALSE, etREAL, {&cutoff},
 686       "Shortest distance (nm) to be considered"},
 687     { "-ng",       FALSE, etINT, {&ngroups},
 688       "Number of secondary groups to compute RDFs around a central group" },
 689     { "-fade",     FALSE, etREAL, {&fade},
 690       "From this distance onwards the RDF is tranformed by g'(r) = 1 + [g(r)-1] exp(-(r/fade-1)^2 to make it go to 1 smoothly. If fade is 0.0 nothing is done." },
 691     { "-grid",     FALSE, etREAL, {&grid},
 692       "[HIDDEN]Grid spacing (in nm) for FFTs when computing structure factors" },
 693     { "-npixel",   FALSE, etINT,  {&npixel},
 694       "[HIDDEN]# pixels per edge of the square detector plate" },
 695     { "-nlevel",   FALSE, etINT,  {&nlevel},
 696       "Number of different colors in the diffraction image" },
 697     { "-distance", FALSE, etREAL, {&distance},
 698       "[HIDDEN]Distance (in cm) from the sample to the detector" },
 699     { "-wave",     FALSE, etREAL, {&lambda},
 700       "[HIDDEN]Wavelength for X-rays/Neutrons for scattering. 0.1 nm corresponds to roughly 12 keV" },
 701
 702     {"-startq", FALSE, etREAL, {&start_q},
 703      "Starting q (1/nm) "},
 704     {"-endq", FALSE, etREAL, {&end_q},
 705      "Ending q (1/nm)"},
 706     {"-energy", FALSE, etREAL, {&energy},
 707      "Energy of the incoming X-ray (keV) "}
 708   };
 709 #define NPA asize(pa)
 710   const char *fnTPS,*fnNDX,*fnDAT=NULL;
 711   gmx_bool       bSQ,bRDF;
 712   output_env_t oenv;
 713
 714   t_filenm   fnm[] = {
 715     { efTRX, "-f",  NULL,     ffREAD },
 716     { efTPS, NULL,  NULL,     ffOPTRD },
 717     { efNDX, NULL,  NULL,     ffOPTRD },
 718     { efDAT, "-d",  "sfactor",     ffOPTRD },
 719     { efXVG, "-o",  "rdf",    ffOPTWR },
 720     { efXVG, "-sq", "sq",     ffOPTWR },
 721     { efXVG, "-cn", "rdf_cn", ffOPTWR },
 722     { efXVG, "-hq", "hq",     ffOPTWR },
 723 /*    { efXPM, "-image", "sq",  ffOPTWR }*/
 724   };
 725 #define NFILE asize(fnm)
 726
 727   CopyRight(stderr,argv[0]);
 728   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 729                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 730
 731   bSQ   = opt2bSet("-sq",NFILE,fnm);
 732   if (bSQ)
 733     please_cite(stdout,"Cromer1968a");
 734
 735   bRDF  = opt2bSet("-o",NFILE,fnm) || !bSQ;
 736   if (bSQ || bCM || closet[0][0]!='n' || rdft[0][0]=='m' || rdft[0][6]=='m') {
 737     fnTPS = ftp2fn(efTPS,NFILE,fnm);
 738     fnDAT = ftp2fn(efDAT,NFILE,fnm);
 739   } else {
 740     fnTPS = ftp2fn_null(efTPS,NFILE,fnm);
 741   }
 742   fnNDX = ftp2fn_null(efNDX,NFILE,fnm);
 743
 744   if (closet[0][0] != 'n') {
 745     if (bCM) {
 746       gmx_fatal(FARGS,"Can not have both -com and -surf");
 747     }
 748     if (bNormalize) {
 749       fprintf(stderr,"Turning of normalization because of option -surf\n");
 750       bNormalize = FALSE;
 751     }
 752   }
 753
 754   if (!bSQ && (!fnTPS && !fnNDX))
 755     gmx_fatal(FARGS,"Neither index file nor topology file specified\n"
 756               "Nothing to do!");
 757
 758   if  (bSQ)
 759    do_scattering_intensity(fnTPS,fnNDX,opt2fn("-sq",NFILE,fnm),
 760                            ftp2fn(efTRX,NFILE,fnm),fnDAT,
 761                            start_q, end_q, energy, ngroups,oenv);
 762
 763   if (bRDF)
 764     do_rdf(fnNDX,fnTPS,ftp2fn(efTRX,NFILE,fnm),
 765            opt2fn("-o",NFILE,fnm),opt2fn_null("-cn",NFILE,fnm),
 766            opt2fn_null("-hq",NFILE,fnm),
 767            bCM,closet[0],rdft,bXY,bPBC,bNormalize,cutoff,binwidth,fade,ngroups,
 768            oenv);
 769
 770   thanx(stderr);
 771
 772   return 0;
 773 }