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