src/gromacs/gmxana/gmx_rdf.c

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