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,2019, 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], ",
  69         "instead of the usual [REF].xpm[ref] file with [TT]-o[tt].",
  70         "[PAR]",
  71         "The default analysis is a 2-D number-density map for a selected",
  72         "group of atoms in the x-y plane.",
  73         "The averaging direction can be changed with the option [TT]-aver[tt].",
  74         "When [TT]-xmin[tt] and/or [TT]-xmax[tt] are set only atoms that are",
  75         "within the limit(s) in the averaging direction are taken into account.",
  76         "The grid spacing is set with the option [TT]-bin[tt].",
  77         "When [TT]-n1[tt] or [TT]-n2[tt] is non-zero, the grid",
  78         "size is set by this option.",
  79         "Box size fluctuations are properly taken into account.",
  80         "[PAR]",
  81         "When options [TT]-amax[tt] and [TT]-rmax[tt] are set, an axial-radial",
  82         "number-density map is made. Three groups should be supplied, the centers",
  83         "of mass of the first two groups define the axis, the third defines the",
  84         "analysis group. The axial direction goes from -amax to +amax, where",
  85         "the center is defined as the midpoint between the centers of mass and",
  86         "the positive direction goes from the first to the second center of mass.",
  87         "The radial direction goes from 0 to rmax or from -rmax to +rmax",
  88         "when the [TT]-mirror[tt] option has been set.",
  89         "[PAR]",
  90         "The normalization of the output is set with the [TT]-unit[tt] option.",
  91         "The default produces a true number density. Unit [TT]nm-2[tt] leaves out",
  92         "the normalization for the averaging or the angular direction.",
  93         "Option [TT]count[tt] produces the count for each grid cell.",
  94         "When you do not want the scale in the output to go",
  95         "from zero to the maximum density, you can set the maximum",
  96         "with the option [TT]-dmax[tt]."
  97     };
  98     static int         n1 = 0, n2 = 0;
  99     static real        xmin = -1, xmax = -1, bin = 0.02, dmin = 0, dmax = 0, amax = 0, rmax = 0;
 100     static gmx_bool    bMirror = FALSE, bSums = FALSE;
 101     static const char* eaver[] = { nullptr, "z", "y", "x", nullptr };
 102     static const char* eunit[] = { nullptr, "nm-3", "nm-2", "count", nullptr };
 103
 104     t_pargs pa[] = {
 105         { "-bin", FALSE, etREAL, { &bin }, "Grid size (nm)" },
 106         { "-aver", FALSE, etENUM, { eaver }, "The direction to average over" },
 107         { "-xmin", FALSE, etREAL, { &xmin }, "Minimum coordinate for averaging" },
 108         { "-xmax", FALSE, etREAL, { &xmax }, "Maximum coordinate for averaging" },
 109         { "-n1", FALSE, etINT, { &n1 }, "Number of grid cells in the first direction" },
 110         { "-n2", FALSE, etINT, { &n2 }, "Number of grid cells in the second direction" },
 111         { "-amax", FALSE, etREAL, { &amax }, "Maximum axial distance from the center" },
 112         { "-rmax", FALSE, etREAL, { &rmax }, "Maximum radial distance" },
 113         { "-mirror", FALSE, etBOOL, { &bMirror }, "Add the mirror image below the axial axis" },
 114         { "-sums", FALSE, etBOOL, { &bSums }, "Print density sums (1D map) to stdout" },
 115         { "-unit", FALSE, etENUM, { eunit }, "Unit for the output" },
 116         { "-dmin", FALSE, etREAL, { &dmin }, "Minimum density in output" },
 117         { "-dmax", FALSE, etREAL, { &dmax }, "Maximum density in output (0 means calculate it)" },
 118     };
 119     gmx_bool          bXmin, bXmax, bRadial;
 120     FILE*             fp;
 121     t_trxstatus*      status;
 122     t_topology        top;
 123     int               ePBC = -1;
 124     rvec *            x, xcom[2], direction, center, dx;
 125     matrix            box;
 126     real              t, m, mtot;
 127     t_pbc             pbc;
 128     int               cav = 0, c1 = 0, c2 = 0;
 129     char **           grpname, buf[STRLEN];
 130     const char*       unit;
 131     int               i, j, k, l, ngrps, anagrp, *gnx = nullptr, nindex, nradial = 0, nfr, nmpower;
 132     int **            ind = nullptr, *index;
 133     real **           grid, maxgrid, m1, m2, box1, box2, *tickx, *tickz, invcellvol;
 134     real              invspa = 0, invspz = 0, axial, r, vol_old, vol, rowsum;
 135     int               nlev = 51;
 136     t_rgb             rlo = { 1, 1, 1 }, rhi = { 0, 0, 0 };
 137     gmx_output_env_t* oenv;
 138     const char*       label[] = { "x (nm)", "y (nm)", "z (nm)" };
 139     t_filenm          fnm[]   = { { efTRX, "-f", nullptr, ffREAD },
 140                        { efTPS, nullptr, nullptr, ffOPTRD },
 141                        { efNDX, nullptr, nullptr, ffOPTRD },
 142                        { efDAT, "-od", "densmap", ffOPTWR },
 143                        { efXPM, "-o", "densmap", ffWRITE } };
 144 #define NFILE asize(fnm)
 145     int npargs;
 146
 147     npargs = asize(pa);
 148
 149     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, npargs, pa,
 150                            asize(desc), desc, 0, nullptr, &oenv))
 151     {
 152         return 0;
 153     }
 154
 155     bXmin   = opt2parg_bSet("-xmin", npargs, pa);
 156     bXmax   = opt2parg_bSet("-xmax", npargs, pa);
 157     bRadial = (amax > 0 || rmax > 0);
 158     if (bRadial)
 159     {
 160         if (amax <= 0 || rmax <= 0)
 161         {
 162             gmx_fatal(FARGS, "Both amax and rmax should be larger than zero");
 163         }
 164     }
 165
 166     GMX_RELEASE_ASSERT(eunit[0] != nullptr, "Option setting inconsistency; eunit[0] is NULL");
 167
 168     if (std::strcmp(eunit[0], "nm-3") == 0)
 169     {
 170         nmpower = -3;
 171         unit    = "(nm^-3)";
 172     }
 173     else if (std::strcmp(eunit[0], "nm-2") == 0)
 174     {
 175         nmpower = -2;
 176         unit    = "(nm^-2)";
 177     }
 178     else
 179     {
 180         nmpower = 0;
 181         unit    = "count";
 182     }
 183
 184     if (ftp2bSet(efTPS, NFILE, fnm) || !ftp2bSet(efNDX, NFILE, fnm))
 185     {
 186         read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &x, nullptr, box, bRadial);
 187     }
 188     if (!bRadial)
 189     {
 190         ngrps = 1;
 191         fprintf(stderr, "\nSelect an analysis group\n");
 192     }
 193     else
 194     {
 195         ngrps = 3;
 196         fprintf(stderr, "\nSelect two groups to define the axis and an analysis group\n");
 197     }
 198     snew(gnx, ngrps);
 199     snew(grpname, ngrps);
 200     snew(ind, ngrps);
 201     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, gnx, ind, grpname);
 202     anagrp = ngrps - 1;
 203     nindex = gnx[anagrp];
 204     index  = ind[anagrp];
 205     if (bRadial)
 206     {
 207         if ((gnx[0] > 1 || gnx[1] > 1) && !ftp2bSet(efTPS, NFILE, fnm))
 208         {
 209             gmx_fatal(FARGS,
 210                       "No run input file was supplied (option -s), this is required for the center "
 211                       "of mass calculation");
 212         }
 213     }
 214
 215     GMX_RELEASE_ASSERT(eaver[0] != nullptr, "Option setting inconsistency; eaver[0] is NULL");
 216
 217     switch (eaver[0][0])
 218     {
 219         case 'x':
 220             cav = XX;
 221             c1  = YY;
 222             c2  = ZZ;
 223             break;
 224         case 'y':
 225             cav = YY;
 226             c1  = XX;
 227             c2  = ZZ;
 228             break;
 229         case 'z':
 230             cav = ZZ;
 231             c1  = XX;
 232             c2  = YY;
 233             break;
 234     }
 235
 236     read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 237
 238     if (!bRadial)
 239     {
 240         if (n1 == 0)
 241         {
 242             n1 = gmx::roundToInt(box[c1][c1] / bin);
 243         }
 244         if (n2 == 0)
 245         {
 246             n2 = gmx::roundToInt(box[c2][c2] / bin);
 247         }
 248     }
 249     else
 250     {
 251         n1      = gmx::roundToInt(2 * amax / bin);
 252         nradial = gmx::roundToInt(rmax / bin);
 253         invspa  = n1 / (2 * amax);
 254         invspz  = nradial / rmax;
 255         if (bMirror)
 256         {
 257             n2 = 2 * nradial;
 258         }
 259         else
 260         {
 261             n2 = nradial;
 262         }
 263     }
 264
 265     snew(grid, n1);
 266     for (i = 0; i < n1; i++)
 267     {
 268         snew(grid[i], n2);
 269     }
 270
 271     box1 = 0;
 272     box2 = 0;
 273     nfr  = 0;
 274     do
 275     {
 276         if (!bRadial)
 277         {
 278             box1 += box[c1][c1];
 279             box2 += box[c2][c2];
 280             invcellvol = n1 * n2;
 281             if (nmpower == -3)
 282             {
 283                 invcellvol /= det(box);
 284             }
 285             else if (nmpower == -2)
 286             {
 287                 invcellvol /= box[c1][c1] * box[c2][c2];
 288             }
 289             for (i = 0; i < nindex; i++)
 290             {
 291                 j = index[i];
 292                 if ((!bXmin || x[j][cav] >= xmin) && (!bXmax || x[j][cav] <= xmax))
 293                 {
 294                     m1 = x[j][c1] / box[c1][c1];
 295                     if (m1 >= 1)
 296                     {
 297                         m1 -= 1;
 298                     }
 299                     if (m1 < 0)
 300                     {
 301                         m1 += 1;
 302                     }
 303                     m2 = x[j][c2] / box[c2][c2];
 304                     if (m2 >= 1)
 305                     {
 306                         m2 -= 1;
 307                     }
 308                     if (m2 < 0)
 309                     {
 310                         m2 += 1;
 311                     }
 312                     grid[static_cast<int>(m1 * n1)][static_cast<int>(m2 * n2)] += invcellvol;
 313                 }
 314             }
 315         }
 316         else
 317         {
 318             set_pbc(&pbc, ePBC, box);
 319             for (i = 0; i < 2; i++)
 320             {
 321                 if (gnx[i] == 1)
 322                 {
 323                     /* One atom, just copy the coordinates */
 324                     copy_rvec(x[ind[i][0]], xcom[i]);
 325                 }
 326                 else
 327                 {
 328                     /* Calculate the center of mass */
 329                     clear_rvec(xcom[i]);
 330                     mtot = 0;
 331                     for (j = 0; j < gnx[i]; j++)
 332                     {
 333                         k = ind[i][j];
 334                         m = top.atoms.atom[k].m;
 335                         for (l = 0; l < DIM; l++)
 336                         {
 337                             xcom[i][l] += m * x[k][l];
 338                         }
 339                         mtot += m;
 340                     }
 341                     svmul(1 / mtot, xcom[i], xcom[i]);
 342                 }
 343             }
 344             pbc_dx(&pbc, xcom[1], xcom[0], direction);
 345             for (i = 0; i < DIM; i++)
 346             {
 347                 center[i] = xcom[0][i] + 0.5 * direction[i];
 348             }
 349             unitv(direction, direction);
 350             for (i = 0; i < nindex; i++)
 351             {
 352                 j = index[i];
 353                 pbc_dx(&pbc, x[j], center, dx);
 354                 axial = iprod(dx, direction);
 355                 r     = std::sqrt(norm2(dx) - axial * axial);
 356                 if (axial >= -amax && axial < amax && r < rmax)
 357                 {
 358                     if (bMirror)
 359                     {
 360                         r += rmax;
 361                     }
 362                     grid[static_cast<int>((axial + amax) * invspa)][static_cast<int>(r * invspz)] += 1;
 363                 }
 364             }
 365         }
 366         nfr++;
 367     } while (read_next_x(oenv, status, &t, x, box));
 368     close_trx(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 + std::strlen(buf), ", %c > %g nm", eaver[0][0], xmin);
 485         }
 486         else if (!bXmin)
 487         {
 488             sprintf(buf + std::strlen(buf), ", %c < %g nm", eaver[0][0], xmax);
 489         }
 490         else
 491         {
 492             sprintf(buf + std::strlen(buf), ", %c: %g - %g nm", eaver[0][0], xmin, xmax);
 493         }
 494     }
 495     if (ftp2bSet(efDAT, NFILE, fnm))
 496     {
 497         fp = gmx_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         gmx_ffclose(fp);
 516     }
 517     else
 518     {
 519         fp = gmx_ffopen(ftp2fn(efXPM, NFILE, fnm), "w");
 520         write_xpm(fp, MAT_SPATIAL_X | MAT_SPATIAL_Y, buf, unit, bRadial ? "axial (nm)" : label[c1],
 521                   bRadial ? "r (nm)" : label[c2], n1, n2, tickx, tickz, grid, dmin, maxgrid, rlo,
 522                   rhi, &nlev);
 523         gmx_ffclose(fp);
 524     }
 525
 526     do_view(oenv, opt2fn("-o", NFILE, fnm), nullptr);
 527
 528     return 0;
 529 }