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