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 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include <cmath>
  41 #include <cstring>
  42
  43 #include "gromacs/commandline/pargs.h"
  44 #include "gromacs/commandline/viewit.h"
  45 #include "gromacs/fileio/confio.h"
  46 #include "gromacs/fileio/matio.h"
  47 #include "gromacs/fileio/trxio.h"
  48 #include "gromacs/gmxana/gmx_ana.h"
  49 #include "gromacs/gmxana/gstat.h"
  50 #include "gromacs/math/utilities.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/pbcutil/pbc.h"
  53 #include "gromacs/topology/index.h"
  54 #include "gromacs/topology/topology.h"
  55 #include "gromacs/utility/arraysize.h"
  56 #include "gromacs/utility/cstringutil.h"
  57 #include "gromacs/utility/fatalerror.h"
  58 #include "gromacs/utility/futil.h"
  59 #include "gromacs/utility/gmxassert.h"
  60 #include "gromacs/utility/smalloc.h"
  61
  62 int gmx_densmap(int argc, char* argv[])
  63 {
  64     const char* desc[] = {
  65         "[THISMODULE] computes 2D number-density maps.",
  66         "It can make planar and axial-radial density maps.",
  67         "The output [REF].xpm[ref] file can be visualized with for instance xv",
  68         "and can be converted to postscript with [TT]xpm2ps[tt].",
  69         "Optionally, output can be in text form to a [REF].dat[ref] file with [TT]-od[tt], ",
  70         "instead of the usual [REF].xpm[ref] file with [TT]-o[tt].",
  71         "[PAR]",
  72         "The default analysis is a 2-D number-density map for a selected",
  73         "group of atoms in the x-y plane.",
  74         "The averaging direction can be changed with the option [TT]-aver[tt].",
  75         "When [TT]-xmin[tt] and/or [TT]-xmax[tt] are set only atoms that are",
  76         "within the limit(s) in the averaging direction are taken into account.",
  77         "The grid spacing is set with the option [TT]-bin[tt].",
  78         "When [TT]-n1[tt] or [TT]-n2[tt] is non-zero, the grid",
  79         "size is set by this option.",
  80         "Box size fluctuations are properly taken into account.",
  81         "[PAR]",
  82         "When options [TT]-amax[tt] and [TT]-rmax[tt] are set, an axial-radial",
  83         "number-density map is made. Three groups should be supplied, the centers",
  84         "of mass of the first two groups define the axis, the third defines the",
  85         "analysis group. The axial direction goes from -amax to +amax, where",
  86         "the center is defined as the midpoint between the centers of mass and",
  87         "the positive direction goes from the first to the second center of mass.",
  88         "The radial direction goes from 0 to rmax or from -rmax to +rmax",
  89         "when the [TT]-mirror[tt] option has been set.",
  90         "[PAR]",
  91         "The normalization of the output is set with the [TT]-unit[tt] option.",
  92         "The default produces a true number density. Unit [TT]nm-2[tt] leaves out",
  93         "the normalization for the averaging or the angular direction.",
  94         "Option [TT]count[tt] produces the count for each grid cell.",
  95         "When you do not want the scale in the output to go",
  96         "from zero to the maximum density, you can set the maximum",
  97         "with the option [TT]-dmax[tt]."
  98     };
  99     static int         n1 = 0, n2 = 0;
 100     static real        xmin = -1, xmax = -1, bin = 0.02, dmin = 0, dmax = 0, amax = 0, rmax = 0;
 101     static gmx_bool    bMirror = FALSE, bSums = FALSE;
 102     static const char* eaver[] = { nullptr, "z", "y", "x", nullptr };
 103     static const char* eunit[] = { nullptr, "nm-3", "nm-2", "count", nullptr };
 104
 105     t_pargs pa[] = {
 106         { "-bin", FALSE, etREAL, { &bin }, "Grid size (nm)" },
 107         { "-aver", FALSE, etENUM, { eaver }, "The direction to average over" },
 108         { "-xmin", FALSE, etREAL, { &xmin }, "Minimum coordinate for averaging" },
 109         { "-xmax", FALSE, etREAL, { &xmax }, "Maximum coordinate for averaging" },
 110         { "-n1", FALSE, etINT, { &n1 }, "Number of grid cells in the first direction" },
 111         { "-n2", FALSE, etINT, { &n2 }, "Number of grid cells in the second direction" },
 112         { "-amax", FALSE, etREAL, { &amax }, "Maximum axial distance from the center" },
 113         { "-rmax", FALSE, etREAL, { &rmax }, "Maximum radial distance" },
 114         { "-mirror", FALSE, etBOOL, { &bMirror }, "Add the mirror image below the axial axis" },
 115         { "-sums", FALSE, etBOOL, { &bSums }, "Print density sums (1D map) to stdout" },
 116         { "-unit", FALSE, etENUM, { eunit }, "Unit for the output" },
 117         { "-dmin", FALSE, etREAL, { &dmin }, "Minimum density in output" },
 118         { "-dmax", FALSE, etREAL, { &dmax }, "Maximum density in output (0 means calculate it)" },
 119     };
 120     gmx_bool          bXmin, bXmax, bRadial;
 121     FILE*             fp;
 122     t_trxstatus*      status;
 123     t_topology        top;
 124     PbcType           pbcType = PbcType::Unset;
 125     rvec *            x, xcom[2], direction, center, dx;
 126     matrix            box;
 127     real              t, m, mtot;
 128     t_pbc             pbc;
 129     int               cav = 0, c1 = 0, c2 = 0;
 130     char **           grpname, buf[STRLEN];
 131     const char*       unit;
 132     int               i, j, k, l, ngrps, anagrp, *gnx = nullptr, nindex, nradial = 0, nfr, nmpower;
 133     int **            ind = nullptr, *index;
 134     real **           grid, maxgrid, m1, m2, box1, box2, *tickx, *tickz, invcellvol;
 135     real              invspa = 0, invspz = 0, axial, r, vol_old, vol, rowsum;
 136     int               nlev = 51;
 137     t_rgb             rlo = { 1, 1, 1 }, rhi = { 0, 0, 0 };
 138     gmx_output_env_t* oenv;
 139     const char*       label[] = { "x (nm)", "y (nm)", "z (nm)" };
 140     t_filenm          fnm[]   = { { efTRX, "-f", nullptr, ffREAD },
 141                        { efTPS, nullptr, nullptr, ffOPTRD },
 142                        { efNDX, nullptr, nullptr, ffOPTRD },
 143                        { efDAT, "-od", "densmap", ffOPTWR },
 144                        { efXPM, "-o", "densmap", ffWRITE } };
 145 #define NFILE asize(fnm)
 146     int npargs;
 147
 148     npargs = asize(pa);
 149
 150     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, npargs, pa,
 151                            asize(desc), desc, 0, nullptr, &oenv))
 152     {
 153         return 0;
 154     }
 155
 156     bXmin   = opt2parg_bSet("-xmin", npargs, pa);
 157     bXmax   = opt2parg_bSet("-xmax", npargs, pa);
 158     bRadial = (amax > 0 || rmax > 0);
 159     if (bRadial)
 160     {
 161         if (amax <= 0 || rmax <= 0)
 162         {
 163             gmx_fatal(FARGS, "Both amax and rmax should be larger than zero");
 164         }
 165     }
 166
 167     GMX_RELEASE_ASSERT(eunit[0] != nullptr, "Option setting inconsistency; eunit[0] is NULL");
 168
 169     if (std::strcmp(eunit[0], "nm-3") == 0)
 170     {
 171         nmpower = -3;
 172         unit    = "(nm^-3)";
 173     }
 174     else if (std::strcmp(eunit[0], "nm-2") == 0)
 175     {
 176         nmpower = -2;
 177         unit    = "(nm^-2)";
 178     }
 179     else
 180     {
 181         nmpower = 0;
 182         unit    = "count";
 183     }
 184
 185     if (ftp2bSet(efTPS, NFILE, fnm) || !ftp2bSet(efNDX, NFILE, fnm))
 186     {
 187         read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, &x, nullptr, box, bRadial);
 188     }
 189     if (!bRadial)
 190     {
 191         ngrps = 1;
 192         fprintf(stderr, "\nSelect an analysis group\n");
 193     }
 194     else
 195     {
 196         ngrps = 3;
 197         fprintf(stderr, "\nSelect two groups to define the axis and an analysis group\n");
 198     }
 199     snew(gnx, ngrps);
 200     snew(grpname, ngrps);
 201     snew(ind, ngrps);
 202     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, gnx, ind, grpname);
 203     anagrp = ngrps - 1;
 204     nindex = gnx[anagrp];
 205     index  = ind[anagrp];
 206     if (bRadial)
 207     {
 208         if ((gnx[0] > 1 || gnx[1] > 1) && !ftp2bSet(efTPS, NFILE, fnm))
 209         {
 210             gmx_fatal(FARGS,
 211                       "No run input file was supplied (option -s), this is required for the center "
 212                       "of mass calculation");
 213         }
 214     }
 215
 216     GMX_RELEASE_ASSERT(eaver[0] != nullptr, "Option setting inconsistency; eaver[0] is NULL");
 217
 218     switch (eaver[0][0])
 219     {
 220         case 'x':
 221             cav = XX;
 222             c1  = YY;
 223             c2  = ZZ;
 224             break;
 225         case 'y':
 226             cav = YY;
 227             c1  = XX;
 228             c2  = ZZ;
 229             break;
 230         case 'z':
 231             cav = ZZ;
 232             c1  = XX;
 233             c2  = YY;
 234             break;
 235     }
 236
 237     read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 238
 239     if (!bRadial)
 240     {
 241         if (n1 == 0)
 242         {
 243             n1 = gmx::roundToInt(box[c1][c1] / bin);
 244         }
 245         if (n2 == 0)
 246         {
 247             n2 = gmx::roundToInt(box[c2][c2] / bin);
 248         }
 249     }
 250     else
 251     {
 252         n1      = gmx::roundToInt(2 * amax / bin);
 253         nradial = gmx::roundToInt(rmax / bin);
 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) && (!bXmax || x[j][cav] <= xmax))
 294                 {
 295                     m1 = x[j][c1] / box[c1][c1];
 296                     if (m1 >= 1)
 297                     {
 298                         m1 -= 1;
 299                     }
 300                     if (m1 < 0)
 301                     {
 302                         m1 += 1;
 303                     }
 304                     m2 = x[j][c2] / box[c2][c2];
 305                     if (m2 >= 1)
 306                     {
 307                         m2 -= 1;
 308                     }
 309                     if (m2 < 0)
 310                     {
 311                         m2 += 1;
 312                     }
 313                     grid[static_cast<int>(m1 * n1)][static_cast<int>(m2 * n2)] += invcellvol;
 314                 }
 315             }
 316         }
 317         else
 318         {
 319             set_pbc(&pbc, pbcType, box);
 320             for (i = 0; i < 2; i++)
 321             {
 322                 if (gnx[i] == 1)
 323                 {
 324                     /* One atom, just copy the coordinates */
 325                     copy_rvec(x[ind[i][0]], xcom[i]);
 326                 }
 327                 else
 328                 {
 329                     /* Calculate the center of mass */
 330                     clear_rvec(xcom[i]);
 331                     mtot = 0;
 332                     for (j = 0; j < gnx[i]; j++)
 333                     {
 334                         k = ind[i][j];
 335                         m = top.atoms.atom[k].m;
 336                         for (l = 0; l < DIM; l++)
 337                         {
 338                             xcom[i][l] += m * x[k][l];
 339                         }
 340                         mtot += m;
 341                     }
 342                     svmul(1 / mtot, xcom[i], xcom[i]);
 343                 }
 344             }
 345             pbc_dx(&pbc, xcom[1], xcom[0], direction);
 346             for (i = 0; i < DIM; i++)
 347             {
 348                 center[i] = xcom[0][i] + 0.5 * direction[i];
 349             }
 350             unitv(direction, direction);
 351             for (i = 0; i < nindex; i++)
 352             {
 353                 j = index[i];
 354                 pbc_dx(&pbc, x[j], center, dx);
 355                 axial = iprod(dx, direction);
 356                 r     = std::sqrt(norm2(dx) - axial * axial);
 357                 if (axial >= -amax && axial < amax && r < rmax)
 358                 {
 359                     if (bMirror)
 360                     {
 361                         r += rmax;
 362                     }
 363                     grid[static_cast<int>((axial + amax) * invspa)][static_cast<int>(r * invspz)] += 1;
 364                 }
 365             }
 366         }
 367         nfr++;
 368     } while (read_next_x(oenv, status, &t, x, box));
 369     close_trx(status);
 370
 371     /* normalize gridpoints */
 372     maxgrid = 0;
 373     if (!bRadial)
 374     {
 375         for (i = 0; i < n1; i++)
 376         {
 377             for (j = 0; j < n2; j++)
 378             {
 379                 grid[i][j] /= nfr;
 380                 if (grid[i][j] > maxgrid)
 381                 {
 382                     maxgrid = grid[i][j];
 383                 }
 384             }
 385         }
 386     }
 387     else
 388     {
 389         for (i = 0; i < n1; i++)
 390         {
 391             vol_old = 0;
 392             for (j = 0; j < nradial; j++)
 393             {
 394                 switch (nmpower)
 395                 {
 396                     case -3: vol = M_PI * (j + 1) * (j + 1) / (invspz * invspz * invspa); break;
 397                     case -2: vol = (j + 1) / (invspz * invspa); break;
 398                     default: vol = j + 1; break;
 399                 }
 400                 if (bMirror)
 401                 {
 402                     k = j + nradial;
 403                 }
 404                 else
 405                 {
 406                     k = j;
 407                 }
 408                 grid[i][k] /= nfr * (vol - vol_old);
 409                 if (bMirror)
 410                 {
 411                     grid[i][nradial - 1 - j] = grid[i][k];
 412                 }
 413                 vol_old = vol;
 414                 if (grid[i][k] > maxgrid)
 415                 {
 416                     maxgrid = grid[i][k];
 417                 }
 418             }
 419         }
 420     }
 421     fprintf(stdout, "\n  The maximum density is %f %s\n", maxgrid, unit);
 422     if (dmax > 0)
 423     {
 424         maxgrid = dmax;
 425     }
 426
 427     snew(tickx, n1 + 1);
 428     snew(tickz, n2 + 1);
 429     if (!bRadial)
 430     {
 431         /* normalize box-axes */
 432         box1 /= nfr;
 433         box2 /= nfr;
 434         for (i = 0; i <= n1; i++)
 435         {
 436             tickx[i] = i * box1 / n1;
 437         }
 438         for (i = 0; i <= n2; i++)
 439         {
 440             tickz[i] = i * box2 / n2;
 441         }
 442     }
 443     else
 444     {
 445         for (i = 0; i <= n1; i++)
 446         {
 447             tickx[i] = i / invspa - amax;
 448         }
 449         if (bMirror)
 450         {
 451             for (i = 0; i <= n2; i++)
 452             {
 453                 tickz[i] = i / invspz - rmax;
 454             }
 455         }
 456         else
 457         {
 458             for (i = 0; i <= n2; i++)
 459             {
 460                 tickz[i] = i / invspz;
 461             }
 462         }
 463     }
 464
 465     if (bSums)
 466     {
 467         for (i = 0; i < n1; ++i)
 468         {
 469             fprintf(stdout, "Density sums:\n");
 470             rowsum = 0;
 471             for (j = 0; j < n2; ++j)
 472             {
 473                 rowsum += grid[i][j];
 474             }
 475             fprintf(stdout, "%g\t", rowsum);
 476         }
 477         fprintf(stdout, "\n");
 478     }
 479
 480     sprintf(buf, "%s number density", grpname[anagrp]);
 481     if (!bRadial && (bXmin || bXmax))
 482     {
 483         if (!bXmax)
 484         {
 485             sprintf(buf + std::strlen(buf), ", %c > %g nm", eaver[0][0], xmin);
 486         }
 487         else if (!bXmin)
 488         {
 489             sprintf(buf + std::strlen(buf), ", %c < %g nm", eaver[0][0], xmax);
 490         }
 491         else
 492         {
 493             sprintf(buf + std::strlen(buf), ", %c: %g - %g nm", eaver[0][0], xmin, xmax);
 494         }
 495     }
 496     if (ftp2bSet(efDAT, NFILE, fnm))
 497     {
 498         fp = gmx_ffopen(ftp2fn(efDAT, NFILE, fnm), "w");
 499         /*optional text form output:  first row is tickz; first col is tickx */
 500         fprintf(fp, "0\t");
 501         for (j = 0; j < n2; ++j)
 502         {
 503             fprintf(fp, "%g\t", tickz[j]);
 504         }
 505         fprintf(fp, "\n");
 506
 507         for (i = 0; i < n1; ++i)
 508         {
 509             fprintf(fp, "%g\t", tickx[i]);
 510             for (j = 0; j < n2; ++j)
 511             {
 512                 fprintf(fp, "%g\t", grid[i][j]);
 513             }
 514             fprintf(fp, "\n");
 515         }
 516         gmx_ffclose(fp);
 517     }
 518     else
 519     {
 520         fp = gmx_ffopen(ftp2fn(efXPM, NFILE, fnm), "w");
 521         write_xpm(fp, MAT_SPATIAL_X | MAT_SPATIAL_Y, buf, unit, bRadial ? "axial (nm)" : label[c1],
 522                   bRadial ? "r (nm)" : label[c2], n1, n2, tickx, tickz, grid, dmin, maxgrid, rlo,
 523                   rhi, &nlev);
 524         gmx_ffclose(fp);
 525     }
 526
 527     do_view(oenv, opt2fn("-o", NFILE, fnm), nullptr);
 528
 529     return 0;
 530 }