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