src/gromacs/gmxana/gmx_densmap.cpp

   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,2014,2015,2016,2017,2018, 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 #include "gmxpre.h"
  38
  39 #include <cmath>
  40 #include <cstring>
  41
  42 #include "gromacs/commandline/pargs.h"
  43 #include "gromacs/commandline/viewit.h"
  44 #include "gromacs/fileio/confio.h"
  45 #include "gromacs/fileio/matio.h"
  46 #include "gromacs/fileio/trxio.h"
  47 #include "gromacs/gmxana/gmx_ana.h"
  48 #include "gromacs/gmxana/gstat.h"
  49 #include "gromacs/math/utilities.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/pbcutil/pbc.h"
  52 #include "gromacs/topology/index.h"
  53 #include "gromacs/topology/topology.h"
  54 #include "gromacs/utility/arraysize.h"
  55 #include "gromacs/utility/cstringutil.h"
  56 #include "gromacs/utility/fatalerror.h"
  57 #include "gromacs/utility/futil.h"
  58 #include "gromacs/utility/gmxassert.h"
  59 #include "gromacs/utility/smalloc.h"
  60
  61 int gmx_densmap(int argc, char *argv[])
  62 {
  63     const char        *desc[] = {
  64         "[THISMODULE] computes 2D number-density maps.",
  65         "It can make planar and axial-radial density maps.",
  66         "The output [REF].xpm[ref] file can be visualized with for instance xv",
  67         "and can be converted to postscript with [TT]xpm2ps[tt].",
  68         "Optionally, output can be in text form to a [REF].dat[ref] file with [TT]-od[tt], instead of the usual [REF].xpm[ref] file with [TT]-o[tt].",
  69         "[PAR]",
  70         "The default analysis is a 2-D number-density map for a selected",
  71         "group of atoms in the x-y plane.",
  72         "The averaging direction can be changed with the option [TT]-aver[tt].",
  73         "When [TT]-xmin[tt] and/or [TT]-xmax[tt] are set only atoms that are",
  74         "within the limit(s) in the averaging direction are taken into account.",
  75         "The grid spacing is set with the option [TT]-bin[tt].",
  76         "When [TT]-n1[tt] or [TT]-n2[tt] is non-zero, the grid",
  77         "size is set by this option.",
  78         "Box size fluctuations are properly taken into account.",
  79         "[PAR]",
  80         "When options [TT]-amax[tt] and [TT]-rmax[tt] are set, an axial-radial",
  81         "number-density map is made. Three groups should be supplied, the centers",
  82         "of mass of the first two groups define the axis, the third defines the",
  83         "analysis group. The axial direction goes from -amax to +amax, where",
  84         "the center is defined as the midpoint between the centers of mass and",
  85         "the positive direction goes from the first to the second center of mass.",
  86         "The radial direction goes from 0 to rmax or from -rmax to +rmax",
  87         "when the [TT]-mirror[tt] option has been set.",
  88         "[PAR]",
  89         "The normalization of the output is set with the [TT]-unit[tt] option.",
  90         "The default produces a true number density. Unit [TT]nm-2[tt] leaves out",
  91         "the normalization for the averaging or the angular direction.",
  92         "Option [TT]count[tt] produces the count for each grid cell.",
  93         "When you do not want the scale in the output to go",
  94         "from zero to the maximum density, you can set the maximum",
  95         "with the option [TT]-dmax[tt]."
  96     };
  97     static int         n1      = 0, n2 = 0;
  98     static real        xmin    = -1, xmax = -1, bin = 0.02, dmin = 0, dmax = 0, amax = 0, rmax = 0;
  99     static gmx_bool    bMirror = FALSE, bSums = FALSE;
 100     static const char *eaver[] = { nullptr, "z", "y", "x", nullptr };
 101     static const char *eunit[] = { nullptr, "nm-3", "nm-2", "count", nullptr };
 102
 103     t_pargs            pa[] = {
 104         { "-bin", FALSE, etREAL, {&bin},
 105           "Grid size (nm)" },
 106         { "-aver", FALSE, etENUM, {eaver},
 107           "The direction to average over" },
 108         { "-xmin", FALSE, etREAL, {&xmin},
 109           "Minimum coordinate for averaging" },
 110         { "-xmax", FALSE, etREAL, {&xmax},
 111           "Maximum coordinate for averaging" },
 112         { "-n1", FALSE, etINT, {&n1},
 113           "Number of grid cells in the first direction" },
 114         { "-n2", FALSE, etINT, {&n2},
 115           "Number of grid cells in the second direction" },
 116         { "-amax", FALSE, etREAL, {&amax},
 117           "Maximum axial distance from the center"},
 118         { "-rmax", FALSE, etREAL, {&rmax},
 119           "Maximum radial distance" },
 120         { "-mirror", FALSE, etBOOL, {&bMirror},
 121           "Add the mirror image below the axial axis" },
 122         { "-sums", FALSE, etBOOL, {&bSums},
 123           "Print density sums (1D map) to stdout" },
 124         { "-unit", FALSE, etENUM, {eunit},
 125           "Unit for the output" },
 126         { "-dmin", FALSE, etREAL, {&dmin},
 127           "Minimum density in output"},
 128         { "-dmax", FALSE, etREAL, {&dmax},
 129           "Maximum density in output (0 means calculate it)"},
 130     };
 131     gmx_bool           bXmin, bXmax, bRadial;
 132     FILE              *fp;
 133     t_trxstatus       *status;
 134     t_topology         top;
 135     int                ePBC = -1;
 136     rvec              *x, xcom[2], direction, center, dx;
 137     matrix             box;
 138     real               t, m, mtot;
 139     t_pbc              pbc;
 140     int                cav = 0, c1 = 0, c2 = 0;
 141     char             **grpname, buf[STRLEN];
 142     const char        *unit;
 143     int                i, j, k, l, ngrps, anagrp, *gnx = nullptr, nindex, nradial = 0, nfr, nmpower;
 144     int              **ind = nullptr, *index;
 145     real             **grid, maxgrid, m1, m2, box1, box2, *tickx, *tickz, invcellvol;
 146     real               invspa = 0, invspz = 0, axial, r, vol_old, vol, rowsum;
 147     int                nlev   = 51;
 148     t_rgb              rlo    = {1, 1, 1}, rhi = {0, 0, 0};
 149     gmx_output_env_t  *oenv;
 150     const char        *label[] = { "x (nm)", "y (nm)", "z (nm)" };
 151     t_filenm           fnm[]   = {
 152         { efTRX, "-f",   nullptr,       ffREAD },
 153         { efTPS, nullptr,   nullptr,       ffOPTRD },
 154         { efNDX, nullptr,   nullptr,       ffOPTRD },
 155         { efDAT, "-od",  "densmap",   ffOPTWR },
 156         { efXPM, "-o",   "densmap",   ffWRITE }
 157     };
 158 #define NFILE asize(fnm)
 159     int                npargs;
 160
 161     npargs = asize(pa);
 162
 163     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
 164                            NFILE, fnm, npargs, pa, asize(desc), desc, 0, nullptr, &oenv))
 165     {
 166         return 0;
 167     }
 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     GMX_RELEASE_ASSERT(eunit[0] != nullptr, "Option setting inconsistency; eunit[0] is NULL");
 181
 182     if (std::strcmp(eunit[0], "nm-3") == 0)
 183     {
 184         nmpower = -3;
 185         unit    = "(nm^-3)";
 186     }
 187     else if (std::strcmp(eunit[0], "nm-2") == 0)
 188     {
 189         nmpower = -2;
 190         unit    = "(nm^-2)";
 191     }
 192     else
 193     {
 194         nmpower = 0;
 195         unit    = "count";
 196     }
 197
 198     if (ftp2bSet(efTPS, NFILE, fnm) || !ftp2bSet(efNDX, NFILE, fnm))
 199     {
 200         read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &x, nullptr, box,
 201                       bRadial);
 202     }
 203     if (!bRadial)
 204     {
 205         ngrps = 1;
 206         fprintf(stderr, "\nSelect an analysis group\n");
 207     }
 208     else
 209     {
 210         ngrps = 3;
 211         fprintf(stderr,
 212                 "\nSelect two groups to define the axis and an analysis group\n");
 213     }
 214     snew(gnx, ngrps);
 215     snew(grpname, ngrps);
 216     snew(ind, ngrps);
 217     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, gnx, ind, grpname);
 218     anagrp = ngrps - 1;
 219     nindex = gnx[anagrp];
 220     index  = ind[anagrp];
 221     if (bRadial)
 222     {
 223         if ((gnx[0] > 1 || gnx[1] > 1) && !ftp2bSet(efTPS, NFILE, fnm))
 224         {
 225             gmx_fatal(FARGS, "No run input file was supplied (option -s), this is required for the center of mass calculation");
 226         }
 227     }
 228
 229     GMX_RELEASE_ASSERT(eaver[0] != nullptr, "Option setting inconsistency; eaver[0] is NULL");
 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     read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 239
 240     if (!bRadial)
 241     {
 242         if (n1 == 0)
 243         {
 244             n1 = gmx::roundToInt(box[c1][c1]/bin);
 245         }
 246         if (n2 == 0)
 247         {
 248             n2 = gmx::roundToInt(box[c2][c2]/bin);
 249         }
 250     }
 251     else
 252     {
 253         n1      = gmx::roundToInt(2*amax/bin);
 254         nradial = gmx::roundToInt(rmax/bin);
 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[static_cast<int>(m1*n1)][static_cast<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     = std::sqrt(norm2(dx) - axial*axial);
 359                 if (axial >= -amax && axial < amax && r < rmax)
 360                 {
 361                     if (bMirror)
 362                     {
 363                         r += rmax;
 364                     }
 365                     grid[static_cast<int>((axial + amax)*invspa)][static_cast<int>(r*invspz)] += 1;
 366                 }
 367             }
 368         }
 369         nfr++;
 370     }
 371     while (read_next_x(oenv, status, &t, x, box));
 372     close_trx(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+std::strlen(buf), ", %c > %g nm", eaver[0][0], xmin);
 489         }
 490         else if (!bXmin)
 491         {
 492             sprintf(buf+std::strlen(buf), ", %c < %g nm", eaver[0][0], xmax);
 493         }
 494         else
 495         {
 496             sprintf(buf+std::strlen(buf), ", %c: %g - %g nm", eaver[0][0], xmin, xmax);
 497         }
 498     }
 499     if (ftp2bSet(efDAT, NFILE, fnm))
 500     {
 501         fp = gmx_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         gmx_ffclose(fp);
 520     }
 521     else
 522     {
 523         fp = gmx_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         gmx_ffclose(fp);
 528     }
 529
 530     do_view(oenv, opt2fn("-o", NFILE, fnm), nullptr);
 531
 532     return 0;
 533 }