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