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