src/gromacs/gmxana/gmx_densmap.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 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include <string.h>
  41
  42 #include "statutil.h"
  43 #include "sysstuff.h"
  44 #include "typedefs.h"
  45 #include "smalloc.h"
  46 #include "macros.h"
  47 #include "vec.h"
  48 #include "pbc.h"
  49 #include "gromacs/fileio/futil.h"
  50 #include "statutil.h"
  51 #include "index.h"
  52 #include "mshift.h"
  53 #include "xvgr.h"
  54 #include "princ.h"
  55 #include "rmpbc.h"
  56 #include "txtdump.h"
  57 #include "gromacs/fileio/tpxio.h"
  58 #include "gromacs/fileio/trxio.h"
  59 #include "gstat.h"
  60 #include "gromacs/fileio/matio.h"
  61 #include "pbc.h"
  62 #include "gmx_ana.h"
  63
  64
  65 int gmx_densmap(int argc, char *argv[])
  66 {
  67     const char        *desc[] = {
  68         "[TT]g_densmap[tt] computes 2D number-density maps.",
  69         "It can make planar and axial-radial density maps.",
  70         "The output [TT].xpm[tt] file can be visualized with for instance xv",
  71         "and can be converted to postscript with [TT]xpm2ps[tt].",
  72         "Optionally, output can be in text form to a [TT].dat[tt] file with [TT]-od[tt], instead of the usual [TT].xpm[tt] file with [TT]-o[tt].",
  73         "[PAR]",
  74         "The default analysis is a 2-D number-density map for a selected",
  75         "group of atoms in the x-y plane.",
  76         "The averaging direction can be changed with the option [TT]-aver[tt].",
  77         "When [TT]-xmin[tt] and/or [TT]-xmax[tt] are set only atoms that are",
  78         "within the limit(s) in the averaging direction are taken into account.",
  79         "The grid spacing is set with the option [TT]-bin[tt].",
  80         "When [TT]-n1[tt] or [TT]-n2[tt] is non-zero, the grid",
  81         "size is set by this option.",
  82         "Box size fluctuations are properly taken into account.",
  83         "[PAR]",
  84         "When options [TT]-amax[tt] and [TT]-rmax[tt] are set, an axial-radial",
  85         "number-density map is made. Three groups should be supplied, the centers",
  86         "of mass of the first two groups define the axis, the third defines the",
  87         "analysis group. The axial direction goes from -amax to +amax, where",
  88         "the center is defined as the midpoint between the centers of mass and",
  89         "the positive direction goes from the first to the second center of mass.",
  90         "The radial direction goes from 0 to rmax or from -rmax to +rmax",
  91         "when the [TT]-mirror[tt] option has been set.",
  92         "[PAR]",
  93         "The normalization of the output is set with the [TT]-unit[tt] option.",
  94         "The default produces a true number density. Unit [TT]nm-2[tt] leaves out",
  95         "the normalization for the averaging or the angular direction.",
  96         "Option [TT]count[tt] produces the count for each grid cell.",
  97         "When you do not want the scale in the output to go",
  98         "from zero to the maximum density, you can set the maximum",
  99         "with the option [TT]-dmax[tt]."
 100     };
 101     static int         n1      = 0, n2 = 0;
 102     static real        xmin    = -1, xmax = -1, bin = 0.02, dmin = 0, dmax = 0, amax = 0, rmax = 0;
 103     static gmx_bool    bMirror = FALSE, bSums = FALSE;
 104     static const char *eaver[] = { NULL, "z", "y", "x", NULL };
 105     static const char *eunit[] = { NULL, "nm-3", "nm-2", "count", NULL };
 106
 107     t_pargs            pa[] = {
 108         { "-bin", FALSE, etREAL, {&bin},
 109           "Grid size (nm)" },
 110         { "-aver", FALSE, etENUM, {eaver},
 111           "The direction to average over" },
 112         { "-xmin", FALSE, etREAL, {&xmin},
 113           "Minimum coordinate for averaging" },
 114         { "-xmax", FALSE, etREAL, {&xmax},
 115           "Maximum coordinate for averaging" },
 116         { "-n1", FALSE, etINT, {&n1},
 117           "Number of grid cells in the first direction" },
 118         { "-n2", FALSE, etINT, {&n2},
 119           "Number of grid cells in the second direction" },
 120         { "-amax", FALSE, etREAL, {&amax},
 121           "Maximum axial distance from the center"},
 122         { "-rmax", FALSE, etREAL, {&rmax},
 123           "Maximum radial distance" },
 124         { "-mirror", FALSE, etBOOL, {&bMirror},
 125           "Add the mirror image below the axial axis" },
 126         { "-sums", FALSE, etBOOL, {&bSums},
 127           "Print density sums (1D map) to stdout" },
 128         { "-unit", FALSE, etENUM, {eunit},
 129           "Unit for the output" },
 130         { "-dmin", FALSE, etREAL, {&dmin},
 131           "Minimum density in output"},
 132         { "-dmax", FALSE, etREAL, {&dmax},
 133           "Maximum density in output (0 means calculate it)"},
 134     };
 135     gmx_bool           bXmin, bXmax, bRadial;
 136     FILE              *fp;
 137     t_trxstatus       *status;
 138     t_topology         top;
 139     int                ePBC = -1;
 140     rvec              *x, xcom[2], direction, center, dx;
 141     matrix             box;
 142     real               t, m, mtot;
 143     t_pbc              pbc;
 144     int                cav = 0, c1 = 0, c2 = 0, natoms;
 145     char             **grpname, title[256], buf[STRLEN];
 146     const char        *unit;
 147     int                i, j, k, l, ngrps, anagrp, *gnx = NULL, nindex, nradial = 0, nfr, nmpower;
 148     atom_id          **ind = NULL, *index;
 149     real             **grid, maxgrid, m1, m2, box1, box2, *tickx, *tickz, invcellvol;
 150     real               invspa = 0, invspz = 0, axial, r, vol_old, vol, rowsum;
 151     int                nlev   = 51;
 152     t_rgb              rlo    = {1, 1, 1}, rhi = {0, 0, 0};
 153     output_env_t       oenv;
 154     const char        *label[] = { "x (nm)", "y (nm)", "z (nm)" };
 155     t_filenm           fnm[]   = {
 156         { efTRX, "-f",   NULL,       ffREAD },
 157         { efTPS, NULL,   NULL,       ffOPTRD },
 158         { efNDX, NULL,   NULL,       ffOPTRD },
 159         { efDAT, "-od",  "densmap",   ffOPTWR },
 160         { efXPM, "-o",   "densmap",   ffWRITE }
 161     };
 162 #define NFILE asize(fnm)
 163     int                npargs;
 164
 165     npargs = asize(pa);
 166
 167     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
 168                            NFILE, fnm, npargs, pa, asize(desc), desc, 0, NULL, &oenv))
 169     {
 170         return 0;
 171     }
 172
 173     bXmin   = opt2parg_bSet("-xmin", npargs, pa);
 174     bXmax   = opt2parg_bSet("-xmax", npargs, pa);
 175     bRadial = (amax > 0 || rmax > 0);
 176     if (bRadial)
 177     {
 178         if (amax <= 0 || rmax <= 0)
 179         {
 180             gmx_fatal(FARGS, "Both amax and rmax should be larger than zero");
 181         }
 182     }
 183
 184     if (strcmp(eunit[0], "nm-3") == 0)
 185     {
 186         nmpower = -3;
 187         unit    = "(nm^-3)";
 188     }
 189     else if (strcmp(eunit[0], "nm-2") == 0)
 190     {
 191         nmpower = -2;
 192         unit    = "(nm^-2)";
 193     }
 194     else
 195     {
 196         nmpower = 0;
 197         unit    = "count";
 198     }
 199
 200     if (ftp2bSet(efTPS, NFILE, fnm) || !ftp2bSet(efNDX, NFILE, fnm))
 201     {
 202         read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &x, NULL, box,
 203                       bRadial);
 204     }
 205     if (!bRadial)
 206     {
 207         ngrps = 1;
 208         fprintf(stderr, "\nSelect an analysis group\n");
 209     }
 210     else
 211     {
 212         ngrps = 3;
 213         fprintf(stderr,
 214                 "\nSelect two groups to define the axis and an analysis group\n");
 215     }
 216     snew(gnx, ngrps);
 217     snew(grpname, ngrps);
 218     snew(ind, ngrps);
 219     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, gnx, ind, grpname);
 220     anagrp = ngrps - 1;
 221     nindex = gnx[anagrp];
 222     index  = ind[anagrp];
 223     if (bRadial)
 224     {
 225         if ((gnx[0] > 1 || gnx[1] > 1) && !ftp2bSet(efTPS, NFILE, fnm))
 226         {
 227             gmx_fatal(FARGS, "No run input file was supplied (option -s), this is required for the center of mass calculation");
 228         }
 229     }
 230
 231     switch (eaver[0][0])
 232     {
 233         case 'x': cav = XX; c1 = YY; c2 = ZZ; break;
 234         case 'y': cav = YY; c1 = XX; c2 = ZZ; break;
 235         case 'z': cav = ZZ; c1 = XX; c2 = YY; break;
 236     }
 237
 238     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 239
 240     if (!bRadial)
 241     {
 242         if (n1 == 0)
 243         {
 244             n1 = (int)(box[c1][c1]/bin + 0.5);
 245         }
 246         if (n2 == 0)
 247         {
 248             n2 = (int)(box[c2][c2]/bin + 0.5);
 249         }
 250     }
 251     else
 252     {
 253         n1      = (int)(2*amax/bin + 0.5);
 254         nradial = (int)(rmax/bin + 0.5);
 255         invspa  = n1/(2*amax);
 256         invspz  = nradial/rmax;
 257         if (bMirror)
 258         {
 259             n2 = 2*nradial;
 260         }
 261         else
 262         {
 263             n2 = nradial;
 264         }
 265     }
 266
 267     snew(grid, n1);
 268     for (i = 0; i < n1; i++)
 269     {
 270         snew(grid[i], n2);
 271     }
 272
 273     box1 = 0;
 274     box2 = 0;
 275     nfr  = 0;
 276     do
 277     {
 278         if (!bRadial)
 279         {
 280             box1      += box[c1][c1];
 281             box2      += box[c2][c2];
 282             invcellvol = n1*n2;
 283             if (nmpower == -3)
 284             {
 285                 invcellvol /= det(box);
 286             }
 287             else if (nmpower == -2)
 288             {
 289                 invcellvol /= box[c1][c1]*box[c2][c2];
 290             }
 291             for (i = 0; i < nindex; i++)
 292             {
 293                 j = index[i];
 294                 if ((!bXmin || x[j][cav] >= xmin) &&
 295                     (!bXmax || x[j][cav] <= xmax))
 296                 {
 297                     m1 = x[j][c1]/box[c1][c1];
 298                     if (m1 >= 1)
 299                     {
 300                         m1 -= 1;
 301                     }
 302                     if (m1 < 0)
 303                     {
 304                         m1 += 1;
 305                     }
 306                     m2 = x[j][c2]/box[c2][c2];
 307                     if (m2 >= 1)
 308                     {
 309                         m2 -= 1;
 310                     }
 311                     if (m2 < 0)
 312                     {
 313                         m2 += 1;
 314                     }
 315                     grid[(int)(m1*n1)][(int)(m2*n2)] += invcellvol;
 316                 }
 317             }
 318         }
 319         else
 320         {
 321             set_pbc(&pbc, ePBC, box);
 322             for (i = 0; i < 2; i++)
 323             {
 324                 if (gnx[i] == 1)
 325                 {
 326                     /* One atom, just copy the coordinates */
 327                     copy_rvec(x[ind[i][0]], xcom[i]);
 328                 }
 329                 else
 330                 {
 331                     /* Calculate the center of mass */
 332                     clear_rvec(xcom[i]);
 333                     mtot = 0;
 334                     for (j = 0; j < gnx[i]; j++)
 335                     {
 336                         k = ind[i][j];
 337                         m = top.atoms.atom[k].m;
 338                         for (l = 0; l < DIM; l++)
 339                         {
 340                             xcom[i][l] += m*x[k][l];
 341                         }
 342                         mtot += m;
 343                     }
 344                     svmul(1/mtot, xcom[i], xcom[i]);
 345                 }
 346             }
 347             pbc_dx(&pbc, xcom[1], xcom[0], direction);
 348             for (i = 0; i < DIM; i++)
 349             {
 350                 center[i] = xcom[0][i] + 0.5*direction[i];
 351             }
 352             unitv(direction, direction);
 353             for (i = 0; i < nindex; i++)
 354             {
 355                 j = index[i];
 356                 pbc_dx(&pbc, x[j], center, dx);
 357                 axial = iprod(dx, direction);
 358                 r     = sqrt(norm2(dx) - axial*axial);
 359                 if (axial >= -amax && axial < amax && r < rmax)
 360                 {
 361                     if (bMirror)
 362                     {
 363                         r += rmax;
 364                     }
 365                     grid[(int)((axial + amax)*invspa)][(int)(r*invspz)] += 1;
 366                 }
 367             }
 368         }
 369         nfr++;
 370     }
 371     while (read_next_x(oenv, status, &t, x, box));
 372     close_trj(status);
 373
 374     /* normalize gridpoints */
 375     maxgrid = 0;
 376     if (!bRadial)
 377     {
 378         for (i = 0; i < n1; i++)
 379         {
 380             for (j = 0; j < n2; j++)
 381             {
 382                 grid[i][j] /= nfr;
 383                 if (grid[i][j] > maxgrid)
 384                 {
 385                     maxgrid = grid[i][j];
 386                 }
 387             }
 388         }
 389     }
 390     else
 391     {
 392         for (i = 0; i < n1; i++)
 393         {
 394             vol_old = 0;
 395             for (j = 0; j < nradial; j++)
 396             {
 397                 switch (nmpower)
 398                 {
 399                     case -3: vol = M_PI*(j+1)*(j+1)/(invspz*invspz*invspa); break;
 400                     case -2: vol =            (j+1)/(invspz*invspa);        break;
 401                     default: vol =             j+1;                         break;
 402                 }
 403                 if (bMirror)
 404                 {
 405                     k = j + nradial;
 406                 }
 407                 else
 408                 {
 409                     k = j;
 410                 }
 411                 grid[i][k] /= nfr*(vol - vol_old);
 412                 if (bMirror)
 413                 {
 414                     grid[i][nradial-1-j] = grid[i][k];
 415                 }
 416                 vol_old = vol;
 417                 if (grid[i][k] > maxgrid)
 418                 {
 419                     maxgrid = grid[i][k];
 420                 }
 421             }
 422         }
 423     }
 424     fprintf(stdout, "\n  The maximum density is %f %s\n", maxgrid, unit);
 425     if (dmax > 0)
 426     {
 427         maxgrid = dmax;
 428     }
 429
 430     snew(tickx, n1+1);
 431     snew(tickz, n2+1);
 432     if (!bRadial)
 433     {
 434         /* normalize box-axes */
 435         box1 /= nfr;
 436         box2 /= nfr;
 437         for (i = 0; i <= n1; i++)
 438         {
 439             tickx[i] = i*box1/n1;
 440         }
 441         for (i = 0; i <= n2; i++)
 442         {
 443             tickz[i] = i*box2/n2;
 444         }
 445     }
 446     else
 447     {
 448         for (i = 0; i <= n1; i++)
 449         {
 450             tickx[i] = i/invspa - amax;
 451         }
 452         if (bMirror)
 453         {
 454             for (i = 0; i <= n2; i++)
 455             {
 456                 tickz[i] = i/invspz - rmax;
 457             }
 458         }
 459         else
 460         {
 461             for (i = 0; i <= n2; i++)
 462             {
 463                 tickz[i] = i/invspz;
 464             }
 465         }
 466     }
 467
 468     if (bSums)
 469     {
 470         for (i = 0; i < n1; ++i)
 471         {
 472             fprintf(stdout, "Density sums:\n");
 473             rowsum = 0;
 474             for (j = 0; j < n2; ++j)
 475             {
 476                 rowsum += grid[i][j];
 477             }
 478             fprintf(stdout, "%g\t", rowsum);
 479         }
 480         fprintf(stdout, "\n");
 481     }
 482
 483     sprintf(buf, "%s number density", grpname[anagrp]);
 484     if (!bRadial && (bXmin || bXmax))
 485     {
 486         if (!bXmax)
 487         {
 488             sprintf(buf+strlen(buf), ", %c > %g nm", eaver[0][0], xmin);
 489         }
 490         else if (!bXmin)
 491         {
 492             sprintf(buf+strlen(buf), ", %c < %g nm", eaver[0][0], xmax);
 493         }
 494         else
 495         {
 496             sprintf(buf+strlen(buf), ", %c: %g - %g nm", eaver[0][0], xmin, xmax);
 497         }
 498     }
 499     if (ftp2bSet(efDAT, NFILE, fnm))
 500     {
 501         fp = ffopen(ftp2fn(efDAT, NFILE, fnm), "w");
 502         /*optional text form output:  first row is tickz; first col is tickx */
 503         fprintf(fp, "0\t");
 504         for (j = 0; j < n2; ++j)
 505         {
 506             fprintf(fp, "%g\t", tickz[j]);
 507         }
 508         fprintf(fp, "\n");
 509
 510         for (i = 0; i < n1; ++i)
 511         {
 512             fprintf(fp, "%g\t", tickx[i]);
 513             for (j = 0; j < n2; ++j)
 514             {
 515                 fprintf(fp, "%g\t", grid[i][j]);
 516             }
 517             fprintf(fp, "\n");
 518         }
 519         ffclose(fp);
 520     }
 521     else
 522     {
 523         fp = ffopen(ftp2fn(efXPM, NFILE, fnm), "w");
 524         write_xpm(fp, MAT_SPATIAL_X | MAT_SPATIAL_Y, buf, unit,
 525                   bRadial ? "axial (nm)" : label[c1], bRadial ? "r (nm)" : label[c2],
 526                   n1, n2, tickx, tickz, grid, dmin, maxgrid, rlo, rhi, &nlev);
 527         ffclose(fp);
 528     }
 529
 530     do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);
 531
 532     return 0;
 533 }