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