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