src/gromacs/gmxana/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 "gromacs/fileio/futil.h"
  50 #include "statutil.h"
  51 #include "gromacs/fileio/tpxio.h"
  52 #include "gromacs/fileio/trxio.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 "gromacs/fileio/matio.h"
  61 #include "gmx_ana.h"
  62 #include "names.h"
  63
  64 static void check_box_c(matrix box)
  65 {
  66     if (fabs(box[ZZ][XX]) > GMX_REAL_EPS*box[ZZ][ZZ] ||
  67         fabs(box[ZZ][YY]) > GMX_REAL_EPS*box[ZZ][ZZ])
  68     {
  69         gmx_fatal(FARGS,
  70                   "The last box vector is not parallel to the z-axis: %f %f %f",
  71                   box[ZZ][XX], box[ZZ][YY], box[ZZ][ZZ]);
  72     }
  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     {
  84         gmx_fatal(FARGS, "No masses available while mass weighting was requested");
  85     }
  86
  87     for (c = 0; c < is; c++)
  88     {
  89         clear_rvec(xc);
  90         mtot = 0;
  91         for (i = coi[c]; i < coi[c+1]; i++)
  92         {
  93             if (bMass)
  94             {
  95                 m = atom[index[i]].m;
  96                 for (d = 0; d < DIM; d++)
  97                 {
  98                     xc[d] += m*x[index[i]][d];
  99                 }
 100                 mtot += m;
 101             }
 102             else
 103             {
 104                 rvec_inc(xc, x[index[i]]);
 105                 mtot += 1.0;
 106             }
 107         }
 108         svmul(1/mtot, xc, x_comg[c]);
 109     }
 110 }
 111
 112 static void split_group(int isize, int *index, char *grpname,
 113                         t_topology *top, char type,
 114                         int *is_out, int **coi_out)
 115 {
 116     t_block *mols = NULL;
 117     t_atom  *atom = NULL;
 118     int      is, *coi;
 119     int      cur, mol, res, i, a, i1;
 120
 121     /* Split up the group in molecules or residues */
 122     switch (type)
 123     {
 124         case 'm':
 125             mols = &top->mols;
 126             break;
 127         case 'r':
 128             atom = top->atoms.atom;
 129             break;
 130         default:
 131             gmx_fatal(FARGS, "Unknown rdf option '%s'", type);
 132     }
 133     snew(coi, isize+1);
 134     is  = 0;
 135     cur = -1;
 136     mol = 0;
 137     for (i = 0; i < isize; i++)
 138     {
 139         a = index[i];
 140         if (type == 'm')
 141         {
 142             /* Check if the molecule number has changed */
 143             i1 = mols->index[mol+1];
 144             while (a >= i1)
 145             {
 146                 mol++;
 147                 i1 = mols->index[mol+1];
 148             }
 149             if (mol != cur)
 150             {
 151                 coi[is++] = i;
 152                 cur       = mol;
 153             }
 154         }
 155         else if (type == 'r')
 156         {
 157             /* Check if the residue index has changed */
 158             res = atom[a].resind;
 159             if (res != cur)
 160             {
 161                 coi[is++] = i;
 162                 cur       = res;
 163             }
 164         }
 165     }
 166     coi[is] = i;
 167     srenew(coi, is+1);
 168     printf("Group '%s' of %d atoms consists of %d %s\n",
 169            grpname, isize, is,
 170            (type == 'm' ? "molecules" : "residues"));
 171
 172     *is_out  = is;
 173     *coi_out = coi;
 174 }
 175
 176 static void do_rdf(const char *fnNDX, const char *fnTPS, const char *fnTRX,
 177                    const char *fnRDF, const char *fnCNRDF, const char *fnHQ,
 178                    gmx_bool bCM, const char *close,
 179                    const char **rdft, gmx_bool bXY, gmx_bool bPBC, gmx_bool bNormalize,
 180                    real cutoff, real binwidth, real fade, int ng,
 181                    const output_env_t oenv)
 182 {
 183     FILE          *fp;
 184     t_trxstatus   *status;
 185     char           outf1[STRLEN], outf2[STRLEN];
 186     char           title[STRLEN], gtitle[STRLEN], refgt[30];
 187     int            g, natoms, i, ii, j, k, nbin, j0, j1, n, nframes;
 188     int          **count;
 189     char         **grpname;
 190     int           *isize, isize_cm = 0, nrdf = 0, max_i, isize0, isize_g;
 191     atom_id      **index, *index_cm = NULL;
 192 #if (defined SIZEOF_LONG_LONG_INT) && (SIZEOF_LONG_LONG_INT >= 8)
 193     long long int *sum;
 194 #else
 195     double        *sum;
 196 #endif
 197     real           t, rmax2, cut2, r, r2, r2ii, invhbinw, normfac;
 198     real           segvol, spherevol, prev_spherevol, **rdf;
 199     rvec          *x, dx, *x0 = NULL, *x_i1, xi;
 200     real          *inv_segvol, invvol, invvol_sum, rho;
 201     gmx_bool       bClose, *bExcl, bTop, bNonSelfExcl;
 202     matrix         box, box_pbc;
 203     int          **npairs;
 204     atom_id        ix, jx, ***pairs;
 205     t_topology    *top  = NULL;
 206     int            ePBC = -1, ePBCrdf = -1;
 207     t_block       *mols = NULL;
 208     t_blocka      *excl;
 209     t_atom        *atom = NULL;
 210     t_pbc          pbc;
 211     gmx_rmpbc_t    gpbc = NULL;
 212     int           *is   = NULL, **coi = NULL, cur, mol, i1, res, a;
 213
 214     excl = NULL;
 215
 216     bClose = (close[0] != 'n');
 217
 218     if (fnTPS)
 219     {
 220         snew(top, 1);
 221         bTop = read_tps_conf(fnTPS, title, top, &ePBC, &x, NULL, box, TRUE);
 222         if (bTop && !(bCM || bClose))
 223         {
 224             /* get exclusions from topology */
 225             excl = &(top->excls);
 226         }
 227     }
 228     snew(grpname, ng+1);
 229     snew(isize, ng+1);
 230     snew(index, ng+1);
 231     fprintf(stderr, "\nSelect a reference group and %d group%s\n",
 232             ng, ng == 1 ? "" : "s");
 233     if (fnTPS)
 234     {
 235         get_index(&(top->atoms), fnNDX, ng+1, isize, index, grpname);
 236         atom = top->atoms.atom;
 237     }
 238     else
 239     {
 240         rd_index(fnNDX, ng+1, isize, index, grpname);
 241     }
 242
 243     if (bCM || bClose)
 244     {
 245         snew(is, 1);
 246         snew(coi, 1);
 247         if (bClose)
 248         {
 249             split_group(isize[0], index[0], grpname[0], top, close[0], &is[0], &coi[0]);
 250         }
 251     }
 252     if (rdft[0][0] != 'a')
 253     {
 254         /* Split up all the groups in molecules or residues */
 255         srenew(is, ng+1);
 256         srenew(coi, ng+1);
 257         for (g = ((bCM || bClose) ? 1 : 0); g < ng+1; g++)
 258         {
 259             split_group(isize[g], index[g], grpname[g], top, rdft[0][0], &is[g], &coi[g]);
 260         }
 261     }
 262
 263     if (bCM)
 264     {
 265         is[0] = 1;
 266         snew(coi[0], is[0]+1);
 267         coi[0][0] = 0;
 268         coi[0][1] = isize[0];
 269         isize0    = is[0];
 270         snew(x0, isize0);
 271     }
 272     else if (bClose || rdft[0][0] != 'a')
 273     {
 274         isize0 = is[0];
 275         snew(x0, isize0);
 276     }
 277     else
 278     {
 279         isize0 = isize[0];
 280     }
 281
 282     natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box);
 283     if (!natoms)
 284     {
 285         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 286     }
 287     if (fnTPS)
 288     {
 289         /* check with topology */
 290         if (natoms > top->atoms.nr)
 291         {
 292             gmx_fatal(FARGS, "Trajectory (%d atoms) does not match topology (%d atoms)",
 293                       natoms, top->atoms.nr);
 294         }
 295     }
 296     /* check with index groups */
 297     for (i = 0; i < ng+1; i++)
 298     {
 299         for (j = 0; j < isize[i]; j++)
 300         {
 301             if (index[i][j] >= natoms)
 302             {
 303                 gmx_fatal(FARGS, "Atom index (%d) in index group %s (%d atoms) larger "
 304                           "than number of atoms in trajectory (%d atoms)",
 305                           index[i][j], grpname[i], isize[i], natoms);
 306             }
 307         }
 308     }
 309
 310     /* initialize some handy things */
 311     if (ePBC == -1)
 312     {
 313         ePBC = guess_ePBC(box);
 314     }
 315     copy_mat(box, box_pbc);
 316     if (bXY)
 317     {
 318         check_box_c(box);
 319         switch (ePBC)
 320         {
 321             case epbcXYZ:
 322             case epbcXY:   ePBCrdf = epbcXY;   break;
 323             case epbcNONE: ePBCrdf = epbcNONE; break;
 324             default:
 325                 gmx_fatal(FARGS, "xy-rdf's are not supported for pbc type'%s'",
 326                           EPBC(ePBC));
 327                 break;
 328         }
 329         /* Make sure the z-height does not influence the cut-off */
 330         box_pbc[ZZ][ZZ] = 2*max(box[XX][XX], box[YY][YY]);
 331     }
 332     else
 333     {
 334         ePBCrdf = ePBC;
 335     }
 336     if (bPBC)
 337     {
 338         rmax2   = 0.99*0.99*max_cutoff2(bXY ? epbcXY : epbcXYZ, box_pbc);
 339     }
 340     else
 341     {
 342         rmax2   = sqr(3*max(box[XX][XX], max(box[YY][YY], box[ZZ][ZZ])));
 343     }
 344     if (debug)
 345     {
 346         fprintf(debug, "rmax2 = %g\n", rmax2);
 347     }
 348
 349     /* We use the double amount of bins, so we can correctly
 350      * write the rdf and rdf_cn output at i*binwidth values.
 351      */
 352     nbin     = (int)(sqrt(rmax2) * 2 / binwidth);
 353     invhbinw = 2.0 / binwidth;
 354     cut2     = sqr(cutoff);
 355
 356     snew(count, ng);
 357     snew(pairs, ng);
 358     snew(npairs, ng);
 359
 360     snew(bExcl, natoms);
 361     max_i = 0;
 362     for (g = 0; g < ng; g++)
 363     {
 364         if (isize[g+1] > max_i)
 365         {
 366             max_i = isize[g+1];
 367         }
 368
 369         /* this is THE array */
 370         snew(count[g], nbin+1);
 371
 372         /* make pairlist array for groups and exclusions */
 373         snew(pairs[g], isize[0]);
 374         snew(npairs[g], isize[0]);
 375         for (i = 0; i < isize[0]; i++)
 376         {
 377             /* We can only have exclusions with atomic rdfs */
 378             if (!(bCM || bClose || rdft[0][0] != 'a'))
 379             {
 380                 ix = index[0][i];
 381                 for (j = 0; j < natoms; j++)
 382                 {
 383                     bExcl[j] = FALSE;
 384                 }
 385                 /* exclusions? */
 386                 if (excl)
 387                 {
 388                     for (j = excl->index[ix]; j < excl->index[ix+1]; j++)
 389                     {
 390                         bExcl[excl->a[j]] = TRUE;
 391                     }
 392                 }
 393                 k = 0;
 394                 snew(pairs[g][i], isize[g+1]);
 395                 bNonSelfExcl = FALSE;
 396                 for (j = 0; j < isize[g+1]; j++)
 397                 {
 398                     jx = index[g+1][j];
 399                     if (!bExcl[jx])
 400                     {
 401                         pairs[g][i][k++] = jx;
 402                     }
 403                     else if (ix != jx)
 404                     {
 405                         /* Check if we have exclusions other than self exclusions */
 406                         bNonSelfExcl = TRUE;
 407                     }
 408                 }
 409                 if (bNonSelfExcl)
 410                 {
 411                     npairs[g][i] = k;
 412                     srenew(pairs[g][i], npairs[g][i]);
 413                 }
 414                 else
 415                 {
 416                     /* Save a LOT of memory and some cpu cycles */
 417                     npairs[g][i] = -1;
 418                     sfree(pairs[g][i]);
 419                 }
 420             }
 421             else
 422             {
 423                 npairs[g][i] = -1;
 424             }
 425         }
 426     }
 427     sfree(bExcl);
 428
 429     snew(x_i1, max_i);
 430     nframes    = 0;
 431     invvol_sum = 0;
 432     if (bPBC && (NULL != top))
 433     {
 434         gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
 435     }
 436
 437     do
 438     {
 439         /* Must init pbc every step because of pressure coupling */
 440         copy_mat(box, box_pbc);
 441         if (bPBC)
 442         {
 443             if (top != NULL)
 444             {
 445                 gmx_rmpbc(gpbc, natoms, box, x);
 446             }
 447             if (bXY)
 448             {
 449                 check_box_c(box);
 450                 clear_rvec(box_pbc[ZZ]);
 451             }
 452             set_pbc(&pbc, ePBCrdf, box_pbc);
 453
 454             if (bXY)
 455             {
 456                 /* Set z-size to 1 so we get the surface iso the volume */
 457                 box_pbc[ZZ][ZZ] = 1;
 458             }
 459         }
 460         invvol      = 1/det(box_pbc);
 461         invvol_sum += invvol;
 462
 463         if (bCM)
 464         {
 465             /* Calculate center of mass of the whole group */
 466             calc_comg(is[0], coi[0], index[0], TRUE, atom, x, x0);
 467         }
 468         else if (!bClose && rdft[0][0] != 'a')
 469         {
 470             calc_comg(is[0], coi[0], index[0], rdft[0][6] == 'm', atom, x, x0);
 471         }
 472
 473         for (g = 0; g < ng; g++)
 474         {
 475             if (rdft[0][0] == 'a')
 476             {
 477                 /* Copy the indexed coordinates to a continuous array */
 478                 for (i = 0; i < isize[g+1]; i++)
 479                 {
 480                     copy_rvec(x[index[g+1][i]], x_i1[i]);
 481                 }
 482             }
 483             else
 484             {
 485                 /* Calculate the COMs/COGs and store in x_i1 */
 486                 calc_comg(is[g+1], coi[g+1], index[g+1], rdft[0][6] == 'm', atom, x, x_i1);
 487             }
 488
 489             for (i = 0; i < isize0; i++)
 490             {
 491                 if (bClose)
 492                 {
 493                     /* Special loop, since we need to determine the minimum distance
 494                      * over all selected atoms in the reference molecule/residue.
 495                      */
 496                     if (rdft[0][0] == 'a')
 497                     {
 498                         isize_g = isize[g+1];
 499                     }
 500                     else
 501                     {
 502                         isize_g = is[g+1];
 503                     }
 504                     for (j = 0; j < isize_g; j++)
 505                     {
 506                         r2 = 1e30;
 507                         /* Loop over the selected atoms in the reference molecule */
 508                         for (ii = coi[0][i]; ii < coi[0][i+1]; ii++)
 509                         {
 510                             if (bPBC)
 511                             {
 512                                 pbc_dx(&pbc, x[index[0][ii]], x_i1[j], dx);
 513                             }
 514                             else
 515                             {
 516                                 rvec_sub(x[index[0][ii]], x_i1[j], dx);
 517                             }
 518                             if (bXY)
 519                             {
 520                                 r2ii = dx[XX]*dx[XX] + dx[YY]*dx[YY];
 521                             }
 522                             else
 523                             {
 524                                 r2ii = iprod(dx, dx);
 525                             }
 526                             if (r2ii < r2)
 527                             {
 528                                 r2 = r2ii;
 529                             }
 530                         }
 531                         if (r2 > cut2 && r2 <= rmax2)
 532                         {
 533                             count[g][(int)(sqrt(r2)*invhbinw)]++;
 534                         }
 535                     }
 536                 }
 537                 else
 538                 {
 539                     /* Real rdf between points in space */
 540                     if (bCM || rdft[0][0] != 'a')
 541                     {
 542                         copy_rvec(x0[i], xi);
 543                     }
 544                     else
 545                     {
 546                         copy_rvec(x[index[0][i]], xi);
 547                     }
 548                     if (rdft[0][0] == 'a' && npairs[g][i] >= 0)
 549                     {
 550                         /* Expensive loop, because of indexing */
 551                         for (j = 0; j < npairs[g][i]; j++)
 552                         {
 553                             jx = pairs[g][i][j];
 554                             if (bPBC)
 555                             {
 556                                 pbc_dx(&pbc, xi, x[jx], dx);
 557                             }
 558                             else
 559                             {
 560                                 rvec_sub(xi, x[jx], dx);
 561                             }
 562
 563                             if (bXY)
 564                             {
 565                                 r2 = dx[XX]*dx[XX] + dx[YY]*dx[YY];
 566                             }
 567                             else
 568                             {
 569                                 r2 = iprod(dx, dx);
 570                             }
 571                             if (r2 > cut2 && r2 <= rmax2)
 572                             {
 573                                 count[g][(int)(sqrt(r2)*invhbinw)]++;
 574                             }
 575                         }
 576                     }
 577                     else
 578                     {
 579                         /* Cheaper loop, no exclusions */
 580                         if (rdft[0][0] == 'a')
 581                         {
 582                             isize_g = isize[g+1];
 583                         }
 584                         else
 585                         {
 586                             isize_g = is[g+1];
 587                         }
 588                         for (j = 0; j < isize_g; j++)
 589                         {
 590                             if (bPBC)
 591                             {
 592                                 pbc_dx(&pbc, xi, x_i1[j], dx);
 593                             }
 594                             else
 595                             {
 596                                 rvec_sub(xi, x_i1[j], dx);
 597                             }
 598                             if (bXY)
 599                             {
 600                                 r2 = dx[XX]*dx[XX] + dx[YY]*dx[YY];
 601                             }
 602                             else
 603                             {
 604                                 r2 = iprod(dx, dx);
 605                             }
 606                             if (r2 > cut2 && r2 <= rmax2)
 607                             {
 608                                 count[g][(int)(sqrt(r2)*invhbinw)]++;
 609                             }
 610                         }
 611                     }
 612                 }
 613             }
 614         }
 615         nframes++;
 616     }
 617     while (read_next_x(oenv, status, &t, x, box));
 618     fprintf(stderr, "\n");
 619
 620     if (bPBC && (NULL != top))
 621     {
 622         gmx_rmpbc_done(gpbc);
 623     }
 624
 625     close_trj(status);
 626
 627     sfree(x);
 628
 629     /* Average volume */
 630     invvol = invvol_sum/nframes;
 631
 632     /* Calculate volume of sphere segments or length of circle segments */
 633     snew(inv_segvol, (nbin+1)/2);
 634     prev_spherevol = 0;
 635     for (i = 0; (i < (nbin+1)/2); i++)
 636     {
 637         r = (i + 0.5)*binwidth;
 638         if (bXY)
 639         {
 640             spherevol = M_PI*r*r;
 641         }
 642         else
 643         {
 644             spherevol = (4.0/3.0)*M_PI*r*r*r;
 645         }
 646         segvol         = spherevol-prev_spherevol;
 647         inv_segvol[i]  = 1.0/segvol;
 648         prev_spherevol = spherevol;
 649     }
 650
 651     snew(rdf, ng);
 652     for (g = 0; g < ng; g++)
 653     {
 654         /* We have to normalize by dividing by the number of frames */
 655         if (rdft[0][0] == 'a')
 656         {
 657             normfac = 1.0/(nframes*invvol*isize0*isize[g+1]);
 658         }
 659         else
 660         {
 661             normfac = 1.0/(nframes*invvol*isize0*is[g+1]);
 662         }
 663
 664         /* Do the normalization */
 665         nrdf = max((nbin+1)/2, 1+2*fade/binwidth);
 666         snew(rdf[g], nrdf);
 667         for (i = 0; i < (nbin+1)/2; i++)
 668         {
 669             r = i*binwidth;
 670             if (i == 0)
 671             {
 672                 j = count[g][0];
 673             }
 674             else
 675             {
 676                 j = count[g][i*2-1] + count[g][i*2];
 677             }
 678             if ((fade > 0) && (r >= fade))
 679             {
 680                 rdf[g][i] = 1 + (j*inv_segvol[i]*normfac-1)*exp(-16*sqr(r/fade-1));
 681             }
 682             else
 683             {
 684                 if (bNormalize)
 685                 {
 686                     rdf[g][i] = j*inv_segvol[i]*normfac;
 687                 }
 688                 else
 689                 {
 690                     rdf[g][i] = j/(binwidth*isize0*nframes);
 691                 }
 692             }
 693         }
 694         for (; (i < nrdf); i++)
 695         {
 696             rdf[g][i] = 1.0;
 697         }
 698     }
 699
 700     if (rdft[0][0] == 'a')
 701     {
 702         sprintf(gtitle, "Radial distribution");
 703     }
 704     else
 705     {
 706         sprintf(gtitle, "Radial distribution of %s %s",
 707                 rdft[0][0] == 'm' ? "molecule" : "residue",
 708                 rdft[0][6] == 'm' ? "COM" : "COG");
 709     }
 710     fp = xvgropen(fnRDF, gtitle, "r", "", oenv);
 711     if (bCM)
 712     {
 713         sprintf(refgt, " %s", "COM");
 714     }
 715     else if (bClose)
 716     {
 717         sprintf(refgt, " closest atom in %s.", close);
 718     }
 719     else
 720     {
 721         sprintf(refgt, "%s", "");
 722     }
 723     if (ng == 1)
 724     {
 725         if (output_env_get_print_xvgr_codes(oenv))
 726         {
 727             fprintf(fp, "@ subtitle \"%s%s - %s\"\n", grpname[0], refgt, grpname[1]);
 728         }
 729     }
 730     else
 731     {
 732         if (output_env_get_print_xvgr_codes(oenv))
 733         {
 734             fprintf(fp, "@ subtitle \"reference %s%s\"\n", grpname[0], refgt);
 735         }
 736         xvgr_legend(fp, ng, (const char**)(grpname+1), oenv);
 737     }
 738     for (i = 0; (i < nrdf); i++)
 739     {
 740         fprintf(fp, "%10g", i*binwidth);
 741         for (g = 0; g < ng; g++)
 742         {
 743             fprintf(fp, " %10g", rdf[g][i]);
 744         }
 745         fprintf(fp, "\n");
 746     }
 747     ffclose(fp);
 748
 749     do_view(oenv, fnRDF, NULL);
 750
 751     /* h(Q) function: fourier transform of rdf */
 752     if (fnHQ)
 753     {
 754         int   nhq = 401;
 755         real *hq, *integrand, Q;
 756
 757         /* Get a better number density later! */
 758         rho = isize[1]*invvol;
 759         snew(hq, nhq);
 760         snew(integrand, nrdf);
 761         for (i = 0; (i < nhq); i++)
 762         {
 763             Q            = i*0.5;
 764             integrand[0] = 0;
 765             for (j = 1; (j < nrdf); j++)
 766             {
 767                 r             = j*binwidth;
 768                 integrand[j]  = (Q == 0) ? 1.0 : sin(Q*r)/(Q*r);
 769                 integrand[j] *= 4.0*M_PI*rho*r*r*(rdf[0][j]-1.0);
 770             }
 771             hq[i] = print_and_integrate(debug, nrdf, binwidth, integrand, NULL, 0);
 772         }
 773         fp = xvgropen(fnHQ, "h(Q)", "Q(/nm)", "h(Q)", oenv);
 774         for (i = 0; (i < nhq); i++)
 775         {
 776             fprintf(fp, "%10g %10g\n", i*0.5, hq[i]);
 777         }
 778         ffclose(fp);
 779         do_view(oenv, fnHQ, NULL);
 780         sfree(hq);
 781         sfree(integrand);
 782     }
 783
 784     if (fnCNRDF)
 785     {
 786         normfac = 1.0/(isize0*nframes);
 787         fp      = xvgropen(fnCNRDF, "Cumulative Number RDF", "r", "number", oenv);
 788         if (ng == 1)
 789         {
 790             if (output_env_get_print_xvgr_codes(oenv))
 791             {
 792                 fprintf(fp, "@ subtitle \"%s-%s\"\n", grpname[0], grpname[1]);
 793             }
 794         }
 795         else
 796         {
 797             if (output_env_get_print_xvgr_codes(oenv))
 798             {
 799                 fprintf(fp, "@ subtitle \"reference %s\"\n", grpname[0]);
 800             }
 801             xvgr_legend(fp, ng, (const char**)(grpname+1), oenv);
 802         }
 803         snew(sum, ng);
 804         for (i = 0; (i <= nbin/2); i++)
 805         {
 806             fprintf(fp, "%10g", i*binwidth);
 807             for (g = 0; g < ng; g++)
 808             {
 809                 fprintf(fp, " %10g", (real)((double)sum[g]*normfac));
 810                 if (i*2+1 < nbin)
 811                 {
 812                     sum[g] += count[g][i*2] + count[g][i*2+1];
 813                 }
 814             }
 815             fprintf(fp, "\n");
 816         }
 817         ffclose(fp);
 818         sfree(sum);
 819
 820         do_view(oenv, fnCNRDF, NULL);
 821     }
 822
 823     for (g = 0; g < ng; g++)
 824     {
 825         sfree(rdf[g]);
 826     }
 827     sfree(rdf);
 828 }
 829
 830
 831 int gmx_rdf(int argc, char *argv[])
 832 {
 833     const char        *desc[] = {
 834         "The structure of liquids can be studied by either neutron or X-ray",
 835         "scattering. The most common way to describe liquid structure is by a",
 836         "radial distribution function. However, this is not easy to obtain from",
 837         "a scattering experiment.[PAR]",
 838         "[THISMODULE] calculates radial distribution functions in different ways.",
 839         "The normal method is around a (set of) particle(s), the other methods",
 840         "are around the center of mass of a set of particles ([TT]-com[tt])",
 841         "or to the closest particle in a set ([TT]-surf[tt]).",
 842         "With all methods, the RDF can also be calculated around axes parallel",
 843         "to the [IT]z[it]-axis with option [TT]-xy[tt].",
 844         "With option [TT]-surf[tt] normalization can not be used.[PAR]",
 845         "The option [TT]-rdf[tt] sets the type of RDF to be computed.",
 846         "Default is for atoms or particles, but one can also select center",
 847         "of mass or geometry of molecules or residues. In all cases, only",
 848         "the atoms in the index groups are taken into account.",
 849         "For molecules and/or the center of mass option, a run input file",
 850         "is required.",
 851         "Weighting other than COM or COG can currently only be achieved",
 852         "by providing a run input file with different masses.",
 853         "Options [TT]-com[tt] and [TT]-surf[tt] also work in conjunction",
 854         "with [TT]-rdf[tt].[PAR]",
 855         "If a run input file is supplied ([TT]-s[tt]) and [TT]-rdf[tt] is set",
 856         "to [TT]atom[tt], exclusions defined",
 857         "in that file are taken into account when calculating the RDF.",
 858         "The option [TT]-cut[tt] is meant as an alternative way to avoid",
 859         "intramolecular peaks in the RDF plot.",
 860         "It is however better to supply a run input file with a higher number of",
 861         "exclusions. For e.g. benzene a topology, setting nrexcl to 5",
 862         "would eliminate all intramolecular contributions to the RDF.",
 863         "Note that all atoms in the selected groups are used, also the ones",
 864         "that don't have Lennard-Jones interactions.[PAR]",
 865         "Option [TT]-cn[tt] produces the cumulative number RDF,",
 866         "i.e. the average number of particles within a distance r.[PAR]"
 867     };
 868     static gmx_bool    bCM     = FALSE, bXY = FALSE, bPBC = TRUE, bNormalize = TRUE;
 869     static real        cutoff  = 0, binwidth = 0.002, fade = 0.0;
 870     static int         ngroups = 1;
 871
 872     static const char *closet[] = { NULL, "no", "mol", "res", NULL };
 873     static const char *rdft[]   = { NULL, "atom", "mol_com", "mol_cog", "res_com", "res_cog", NULL };
 874
 875     t_pargs            pa[] = {
 876         { "-bin",      FALSE, etREAL, {&binwidth},
 877           "Binwidth (nm)" },
 878         { "-com",      FALSE, etBOOL, {&bCM},
 879           "RDF with respect to the center of mass of first group" },
 880         { "-surf",     FALSE, etENUM, {closet},
 881           "RDF with respect to the surface of the first group" },
 882         { "-rdf",   FALSE, etENUM, {rdft},
 883           "RDF type" },
 884         { "-pbc",      FALSE, etBOOL, {&bPBC},
 885           "Use periodic boundary conditions for computing distances. Without PBC the maximum range will be three times the largest box edge." },
 886         { "-norm",     FALSE, etBOOL, {&bNormalize},
 887           "Normalize for volume and density" },
 888         { "-xy",       FALSE, etBOOL, {&bXY},
 889           "Use only the x and y components of the distance" },
 890         { "-cut",      FALSE, etREAL, {&cutoff},
 891           "Shortest distance (nm) to be considered"},
 892         { "-ng",       FALSE, etINT, {&ngroups},
 893           "Number of secondary groups to compute RDFs around a central group" },
 894         { "-fade",     FALSE, etREAL, {&fade},
 895           "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." }
 896     };
 897 #define NPA asize(pa)
 898     const char        *fnTPS, *fnNDX;
 899     output_env_t       oenv;
 900
 901     t_filenm           fnm[] = {
 902         { efTRX, "-f",  NULL,     ffREAD },
 903         { efTPS, NULL,  NULL,     ffOPTRD },
 904         { efNDX, NULL,  NULL,     ffOPTRD },
 905         { efXVG, "-o",  "rdf",    ffWRITE },
 906         { efXVG, "-cn", "rdf_cn", ffOPTWR },
 907         { efXVG, "-hq", "hq",     ffOPTWR },
 908     };
 909 #define NFILE asize(fnm)
 910     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 911                            NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
 912     {
 913         return 0;
 914     }
 915
 916     if (bCM || closet[0][0] != 'n' || rdft[0][0] == 'm' || rdft[0][6] == 'm')
 917     {
 918         fnTPS = ftp2fn(efTPS, NFILE, fnm);
 919     }
 920     else
 921     {
 922         fnTPS = ftp2fn_null(efTPS, NFILE, fnm);
 923     }
 924     fnNDX = ftp2fn_null(efNDX, NFILE, fnm);
 925
 926     if (!fnTPS && !fnNDX)
 927     {
 928         gmx_fatal(FARGS, "Neither index file nor topology file specified\n"
 929                   "Nothing to do!");
 930     }
 931
 932     if (closet[0][0] != 'n')
 933     {
 934         if (bCM)
 935         {
 936             gmx_fatal(FARGS, "Can not have both -com and -surf");
 937         }
 938         if (bNormalize)
 939         {
 940             fprintf(stderr, "Turning of normalization because of option -surf\n");
 941             bNormalize = FALSE;
 942         }
 943     }
 944
 945     do_rdf(fnNDX, fnTPS, ftp2fn(efTRX, NFILE, fnm),
 946            opt2fn("-o", NFILE, fnm), opt2fn_null("-cn", NFILE, fnm),
 947            opt2fn_null("-hq", NFILE, fnm),
 948            bCM, closet[0], rdft, bXY, bPBC, bNormalize, cutoff, binwidth, fade, ngroups,
 949            oenv);
 950
 951     return 0;
 952 }