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     if (!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         return 0;
 170     }
 171
 172     bXmin   = opt2parg_bSet("-xmin", npargs, pa);
 173     bXmax   = opt2parg_bSet("-xmax", npargs, pa);
 174     bRadial = (amax > 0 || rmax > 0);
 175     if (bRadial)
 176     {
 177         if (amax <= 0 || rmax <= 0)
 178         {
 179             gmx_fatal(FARGS, "Both amax and rmax should be larger than zero");
 180         }
 181     }
 182
 183     if (strcmp(eunit[0], "nm-3") == 0)
 184     {
 185         nmpower = -3;
 186         unit    = "(nm^-3)";
 187     }
 188     else if (strcmp(eunit[0], "nm-2") == 0)
 189     {
 190         nmpower = -2;
 191         unit    = "(nm^-2)";
 192     }
 193     else
 194     {
 195         nmpower = 0;
 196         unit    = "count";
 197     }
 198
 199     if (ftp2bSet(efTPS, NFILE, fnm) || !ftp2bSet(efNDX, NFILE, fnm))
 200     {
 201         read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &x, NULL, box,
 202                       bRadial);
 203     }
 204     if (!bRadial)
 205     {
 206         ngrps = 1;
 207         fprintf(stderr, "\nSelect an analysis group\n");
 208     }
 209     else
 210     {
 211         ngrps = 3;
 212         fprintf(stderr,
 213                 "\nSelect two groups to define the axis and an analysis group\n");
 214     }
 215     snew(gnx, ngrps);
 216     snew(grpname, ngrps);
 217     snew(ind, ngrps);
 218     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, gnx, ind, grpname);
 219     anagrp = ngrps - 1;
 220     nindex = gnx[anagrp];
 221     index  = ind[anagrp];
 222     if (bRadial)
 223     {
 224         if ((gnx[0] > 1 || gnx[1] > 1) && !ftp2bSet(efTPS, NFILE, fnm))
 225         {
 226             gmx_fatal(FARGS, "No run input file was supplied (option -s), this is required for the center of mass calculation");
 227         }
 228     }
 229
 230     switch (eaver[0][0])
 231     {
 232         case 'x': cav = XX; c1 = YY; c2 = ZZ; break;
 233         case 'y': cav = YY; c1 = XX; c2 = ZZ; break;
 234         case 'z': cav = ZZ; c1 = XX; c2 = YY; break;
 235     }
 236
 237     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 238
 239     if (!bRadial)
 240     {
 241         if (n1 == 0)
 242         {
 243             n1 = (int)(box[c1][c1]/bin + 0.5);
 244         }
 245         if (n2 == 0)
 246         {
 247             n2 = (int)(box[c2][c2]/bin + 0.5);
 248         }
 249     }
 250     else
 251     {
 252         n1      = (int)(2*amax/bin + 0.5);
 253         nradial = (int)(rmax/bin + 0.5);
 254         invspa  = n1/(2*amax);
 255         invspz  = nradial/rmax;
 256         if (bMirror)
 257         {
 258             n2 = 2*nradial;
 259         }
 260         else
 261         {
 262             n2 = nradial;
 263         }
 264     }
 265
 266     snew(grid, n1);
 267     for (i = 0; i < n1; i++)
 268     {
 269         snew(grid[i], n2);
 270     }
 271
 272     box1 = 0;
 273     box2 = 0;
 274     nfr  = 0;
 275     do
 276     {
 277         if (!bRadial)
 278         {
 279             box1      += box[c1][c1];
 280             box2      += box[c2][c2];
 281             invcellvol = n1*n2;
 282             if (nmpower == -3)
 283             {
 284                 invcellvol /= det(box);
 285             }
 286             else if (nmpower == -2)
 287             {
 288                 invcellvol /= box[c1][c1]*box[c2][c2];
 289             }
 290             for (i = 0; i < nindex; i++)
 291             {
 292                 j = index[i];
 293                 if ((!bXmin || x[j][cav] >= xmin) &&
 294                     (!bXmax || x[j][cav] <= xmax))
 295                 {
 296                     m1 = x[j][c1]/box[c1][c1];
 297                     if (m1 >= 1)
 298                     {
 299                         m1 -= 1;
 300                     }
 301                     if (m1 < 0)
 302                     {
 303                         m1 += 1;
 304                     }
 305                     m2 = x[j][c2]/box[c2][c2];
 306                     if (m2 >= 1)
 307                     {
 308                         m2 -= 1;
 309                     }
 310                     if (m2 < 0)
 311                     {
 312                         m2 += 1;
 313                     }
 314                     grid[(int)(m1*n1)][(int)(m2*n2)] += invcellvol;
 315                 }
 316             }
 317         }
 318         else
 319         {
 320             set_pbc(&pbc, ePBC, box);
 321             for (i = 0; i < 2; i++)
 322             {
 323                 if (gnx[i] == 1)
 324                 {
 325                     /* One atom, just copy the coordinates */
 326                     copy_rvec(x[ind[i][0]], xcom[i]);
 327                 }
 328                 else
 329                 {
 330                     /* Calculate the center of mass */
 331                     clear_rvec(xcom[i]);
 332                     mtot = 0;
 333                     for (j = 0; j < gnx[i]; j++)
 334                     {
 335                         k = ind[i][j];
 336                         m = top.atoms.atom[k].m;
 337                         for (l = 0; l < DIM; l++)
 338                         {
 339                             xcom[i][l] += m*x[k][l];
 340                         }
 341                         mtot += m;
 342                     }
 343                     svmul(1/mtot, xcom[i], xcom[i]);
 344                 }
 345             }
 346             pbc_dx(&pbc, xcom[1], xcom[0], direction);
 347             for (i = 0; i < DIM; i++)
 348             {
 349                 center[i] = xcom[0][i] + 0.5*direction[i];
 350             }
 351             unitv(direction, direction);
 352             for (i = 0; i < nindex; i++)
 353             {
 354                 j = index[i];
 355                 pbc_dx(&pbc, x[j], center, dx);
 356                 axial = iprod(dx, direction);
 357                 r     = sqrt(norm2(dx) - axial*axial);
 358                 if (axial >= -amax && axial < amax && r < rmax)
 359                 {
 360                     if (bMirror)
 361                     {
 362                         r += rmax;
 363                     }
 364                     grid[(int)((axial + amax)*invspa)][(int)(r*invspz)] += 1;
 365                 }
 366             }
 367         }
 368         nfr++;
 369     }
 370     while (read_next_x(oenv, status, &t, x, box));
 371     close_trj(status);
 372
 373     /* normalize gridpoints */
 374     maxgrid = 0;
 375     if (!bRadial)
 376     {
 377         for (i = 0; i < n1; i++)
 378         {
 379             for (j = 0; j < n2; j++)
 380             {
 381                 grid[i][j] /= nfr;
 382                 if (grid[i][j] > maxgrid)
 383                 {
 384                     maxgrid = grid[i][j];
 385                 }
 386             }
 387         }
 388     }
 389     else
 390     {
 391         for (i = 0; i < n1; i++)
 392         {
 393             vol_old = 0;
 394             for (j = 0; j < nradial; j++)
 395             {
 396                 switch (nmpower)
 397                 {
 398                     case -3: vol = M_PI*(j+1)*(j+1)/(invspz*invspz*invspa); break;
 399                     case -2: vol =            (j+1)/(invspz*invspa);        break;
 400                     default: vol =             j+1;                         break;
 401                 }
 402                 if (bMirror)
 403                 {
 404                     k = j + nradial;
 405                 }
 406                 else
 407                 {
 408                     k = j;
 409                 }
 410                 grid[i][k] /= nfr*(vol - vol_old);
 411                 if (bMirror)
 412                 {
 413                     grid[i][nradial-1-j] = grid[i][k];
 414                 }
 415                 vol_old = vol;
 416                 if (grid[i][k] > maxgrid)
 417                 {
 418                     maxgrid = grid[i][k];
 419                 }
 420             }
 421         }
 422     }
 423     fprintf(stdout, "\n  The maximum density is %f %s\n", maxgrid, unit);
 424     if (dmax > 0)
 425     {
 426         maxgrid = dmax;
 427     }
 428
 429     snew(tickx, n1+1);
 430     snew(tickz, n2+1);
 431     if (!bRadial)
 432     {
 433         /* normalize box-axes */
 434         box1 /= nfr;
 435         box2 /= nfr;
 436         for (i = 0; i <= n1; i++)
 437         {
 438             tickx[i] = i*box1/n1;
 439         }
 440         for (i = 0; i <= n2; i++)
 441         {
 442             tickz[i] = i*box2/n2;
 443         }
 444     }
 445     else
 446     {
 447         for (i = 0; i <= n1; i++)
 448         {
 449             tickx[i] = i/invspa - amax;
 450         }
 451         if (bMirror)
 452         {
 453             for (i = 0; i <= n2; i++)
 454             {
 455                 tickz[i] = i/invspz - rmax;
 456             }
 457         }
 458         else
 459         {
 460             for (i = 0; i <= n2; i++)
 461             {
 462                 tickz[i] = i/invspz;
 463             }
 464         }
 465     }
 466
 467     if (bSums)
 468     {
 469         for (i = 0; i < n1; ++i)
 470         {
 471             fprintf(stdout, "Density sums:\n");
 472             rowsum = 0;
 473             for (j = 0; j < n2; ++j)
 474             {
 475                 rowsum += grid[i][j];
 476             }
 477             fprintf(stdout, "%g\t", rowsum);
 478         }
 479         fprintf(stdout, "\n");
 480     }
 481
 482     sprintf(buf, "%s number density", grpname[anagrp]);
 483     if (!bRadial && (bXmin || bXmax))
 484     {
 485         if (!bXmax)
 486         {
 487             sprintf(buf+strlen(buf), ", %c > %g nm", eaver[0][0], xmin);
 488         }
 489         else if (!bXmin)
 490         {
 491             sprintf(buf+strlen(buf), ", %c < %g nm", eaver[0][0], xmax);
 492         }
 493         else
 494         {
 495             sprintf(buf+strlen(buf), ", %c: %g - %g nm", eaver[0][0], xmin, xmax);
 496         }
 497     }
 498     if (ftp2bSet(efDAT, NFILE, fnm))
 499     {
 500         fp = ffopen(ftp2fn(efDAT, NFILE, fnm), "w");
 501         /*optional text form output:  first row is tickz; first col is tickx */
 502         fprintf(fp, "0\t");
 503         for (j = 0; j < n2; ++j)
 504         {
 505             fprintf(fp, "%g\t", tickz[j]);
 506         }
 507         fprintf(fp, "\n");
 508
 509         for (i = 0; i < n1; ++i)
 510         {
 511             fprintf(fp, "%g\t", tickx[i]);
 512             for (j = 0; j < n2; ++j)
 513             {
 514                 fprintf(fp, "%g\t", grid[i][j]);
 515             }
 516             fprintf(fp, "\n");
 517         }
 518         ffclose(fp);
 519     }
 520     else
 521     {
 522         fp = ffopen(ftp2fn(efXPM, NFILE, fnm), "w");
 523         write_xpm(fp, MAT_SPATIAL_X | MAT_SPATIAL_Y, buf, unit,
 524                   bRadial ? "axial (nm)" : label[c1], bRadial ? "r (nm)" : label[c2],
 525                   n1, n2, tickx, tickz, grid, dmin, maxgrid, rlo, rhi, &nlev);
 526         ffclose(fp);
 527     }
 528
 529     do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);
 530
 531     return 0;
 532 }