src/tools/gmx_rdf.c

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