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, 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/fileio/confio.h"
  44 #include "gromacs/fileio/matio.h"
  45 #include "gromacs/fileio/trxio.h"
  46 #include "gromacs/gmxana/gmx_ana.h"
  47 #include "gromacs/gmxana/gstat.h"
  48 #include "gromacs/legacyheaders/macros.h"
  49 #include "gromacs/legacyheaders/txtdump.h"
  50 #include "gromacs/legacyheaders/typedefs.h"
  51 #include "gromacs/legacyheaders/viewit.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/pbcutil/pbc.h"
  54 #include "gromacs/topology/index.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[] = { NULL, "z", "y", "x", NULL };
 101     static const char *eunit[] = { NULL, "nm-3", "nm-2", "count", NULL };
 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, title[256], buf[STRLEN];
 142     const char        *unit;
 143     int                i, j, k, l, ngrps, anagrp, *gnx = NULL, nindex, nradial = 0, nfr, nmpower;
 144     atom_id          **ind = NULL, *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     output_env_t       oenv;
 150     const char        *label[] = { "x (nm)", "y (nm)", "z (nm)" };
 151     t_filenm           fnm[]   = {
 152         { efTRX, "-f",   NULL,       ffREAD },
 153         { efTPS, NULL,   NULL,       ffOPTRD },
 154         { efNDX, NULL,   NULL,       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, NULL, &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] != NULL, "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), title, &top, &ePBC, &x, NULL, 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] != NULL, "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 = static_cast<int>(box[c1][c1]/bin + 0.5);
 245         }
 246         if (n2 == 0)
 247         {
 248             n2 = static_cast<int>(box[c2][c2]/bin + 0.5);
 249         }
 250     }
 251     else
 252     {
 253         n1      = static_cast<int>(2*amax/bin + 0.5);
 254         nradial = static_cast<int>(rmax/bin + 0.5);
 255         /* cppcheck-suppress zerodiv fixed in 1.68-dev */
 256         invspa  = n1/(2*amax);
 257         /* cppcheck-suppress zerodiv fixed in 1.68-dev */
 258         invspz  = nradial/rmax;
 259         if (bMirror)
 260         {
 261             n2 = 2*nradial;
 262         }
 263         else
 264         {
 265             n2 = nradial;
 266         }
 267     }
 268
 269     snew(grid, n1);
 270     for (i = 0; i < n1; i++)
 271     {
 272         snew(grid[i], n2);
 273     }
 274
 275     box1 = 0;
 276     box2 = 0;
 277     nfr  = 0;
 278     do
 279     {
 280         if (!bRadial)
 281         {
 282             box1      += box[c1][c1];
 283             box2      += box[c2][c2];
 284             invcellvol = n1*n2;
 285             if (nmpower == -3)
 286             {
 287                 invcellvol /= det(box);
 288             }
 289             else if (nmpower == -2)
 290             {
 291                 invcellvol /= box[c1][c1]*box[c2][c2];
 292             }
 293             for (i = 0; i < nindex; i++)
 294             {
 295                 j = index[i];
 296                 if ((!bXmin || x[j][cav] >= xmin) &&
 297                     (!bXmax || x[j][cav] <= xmax))
 298                 {
 299                     m1 = x[j][c1]/box[c1][c1];
 300                     if (m1 >= 1)
 301                     {
 302                         m1 -= 1;
 303                     }
 304                     if (m1 < 0)
 305                     {
 306                         m1 += 1;
 307                     }
 308                     m2 = x[j][c2]/box[c2][c2];
 309                     if (m2 >= 1)
 310                     {
 311                         m2 -= 1;
 312                     }
 313                     if (m2 < 0)
 314                     {
 315                         m2 += 1;
 316                     }
 317                     grid[static_cast<int>(m1*n1)][static_cast<int>(m2*n2)] += invcellvol;
 318                 }
 319             }
 320         }
 321         else
 322         {
 323             set_pbc(&pbc, ePBC, box);
 324             for (i = 0; i < 2; i++)
 325             {
 326                 if (gnx[i] == 1)
 327                 {
 328                     /* One atom, just copy the coordinates */
 329                     copy_rvec(x[ind[i][0]], xcom[i]);
 330                 }
 331                 else
 332                 {
 333                     /* Calculate the center of mass */
 334                     clear_rvec(xcom[i]);
 335                     mtot = 0;
 336                     for (j = 0; j < gnx[i]; j++)
 337                     {
 338                         k = ind[i][j];
 339                         m = top.atoms.atom[k].m;
 340                         for (l = 0; l < DIM; l++)
 341                         {
 342                             xcom[i][l] += m*x[k][l];
 343                         }
 344                         mtot += m;
 345                     }
 346                     svmul(1/mtot, xcom[i], xcom[i]);
 347                 }
 348             }
 349             pbc_dx(&pbc, xcom[1], xcom[0], direction);
 350             for (i = 0; i < DIM; i++)
 351             {
 352                 center[i] = xcom[0][i] + 0.5*direction[i];
 353             }
 354             unitv(direction, direction);
 355             for (i = 0; i < nindex; i++)
 356             {
 357                 j = index[i];
 358                 pbc_dx(&pbc, x[j], center, dx);
 359                 axial = iprod(dx, direction);
 360                 r     = std::sqrt(norm2(dx) - axial*axial);
 361                 if (axial >= -amax && axial < amax && r < rmax)
 362                 {
 363                     if (bMirror)
 364                     {
 365                         r += rmax;
 366                     }
 367                     grid[static_cast<int>((axial + amax)*invspa)][static_cast<int>(r*invspz)] += 1;
 368                 }
 369             }
 370         }
 371         nfr++;
 372     }
 373     while (read_next_x(oenv, status, &t, x, box));
 374     close_trj(status);
 375
 376     /* normalize gridpoints */
 377     maxgrid = 0;
 378     if (!bRadial)
 379     {
 380         for (i = 0; i < n1; i++)
 381         {
 382             for (j = 0; j < n2; j++)
 383             {
 384                 grid[i][j] /= nfr;
 385                 if (grid[i][j] > maxgrid)
 386                 {
 387                     maxgrid = grid[i][j];
 388                 }
 389             }
 390         }
 391     }
 392     else
 393     {
 394         for (i = 0; i < n1; i++)
 395         {
 396             vol_old = 0;
 397             for (j = 0; j < nradial; j++)
 398             {
 399                 switch (nmpower)
 400                 {
 401                     case -3: vol = M_PI*(j+1)*(j+1)/(invspz*invspz*invspa); break;
 402                     case -2: vol =            (j+1)/(invspz*invspa);        break;
 403                     default: vol =             j+1;                         break;
 404                 }
 405                 if (bMirror)
 406                 {
 407                     k = j + nradial;
 408                 }
 409                 else
 410                 {
 411                     k = j;
 412                 }
 413                 grid[i][k] /= nfr*(vol - vol_old);
 414                 if (bMirror)
 415                 {
 416                     grid[i][nradial-1-j] = grid[i][k];
 417                 }
 418                 vol_old = vol;
 419                 if (grid[i][k] > maxgrid)
 420                 {
 421                     maxgrid = grid[i][k];
 422                 }
 423             }
 424         }
 425     }
 426     fprintf(stdout, "\n  The maximum density is %f %s\n", maxgrid, unit);
 427     if (dmax > 0)
 428     {
 429         maxgrid = dmax;
 430     }
 431
 432     snew(tickx, n1+1);
 433     snew(tickz, n2+1);
 434     if (!bRadial)
 435     {
 436         /* normalize box-axes */
 437         box1 /= nfr;
 438         box2 /= nfr;
 439         for (i = 0; i <= n1; i++)
 440         {
 441             tickx[i] = i*box1/n1;
 442         }
 443         for (i = 0; i <= n2; i++)
 444         {
 445             tickz[i] = i*box2/n2;
 446         }
 447     }
 448     else
 449     {
 450         for (i = 0; i <= n1; i++)
 451         {
 452             tickx[i] = i/invspa - amax;
 453         }
 454         if (bMirror)
 455         {
 456             for (i = 0; i <= n2; i++)
 457             {
 458                 tickz[i] = i/invspz - rmax;
 459             }
 460         }
 461         else
 462         {
 463             for (i = 0; i <= n2; i++)
 464             {
 465                 tickz[i] = i/invspz;
 466             }
 467         }
 468     }
 469
 470     if (bSums)
 471     {
 472         for (i = 0; i < n1; ++i)
 473         {
 474             fprintf(stdout, "Density sums:\n");
 475             rowsum = 0;
 476             for (j = 0; j < n2; ++j)
 477             {
 478                 rowsum += grid[i][j];
 479             }
 480             fprintf(stdout, "%g\t", rowsum);
 481         }
 482         fprintf(stdout, "\n");
 483     }
 484
 485     sprintf(buf, "%s number density", grpname[anagrp]);
 486     if (!bRadial && (bXmin || bXmax))
 487     {
 488         if (!bXmax)
 489         {
 490             sprintf(buf+std::strlen(buf), ", %c > %g nm", eaver[0][0], xmin);
 491         }
 492         else if (!bXmin)
 493         {
 494             sprintf(buf+std::strlen(buf), ", %c < %g nm", eaver[0][0], xmax);
 495         }
 496         else
 497         {
 498             sprintf(buf+std::strlen(buf), ", %c: %g - %g nm", eaver[0][0], xmin, xmax);
 499         }
 500     }
 501     if (ftp2bSet(efDAT, NFILE, fnm))
 502     {
 503         fp = gmx_ffopen(ftp2fn(efDAT, NFILE, fnm), "w");
 504         /*optional text form output:  first row is tickz; first col is tickx */
 505         fprintf(fp, "0\t");
 506         for (j = 0; j < n2; ++j)
 507         {
 508             fprintf(fp, "%g\t", tickz[j]);
 509         }
 510         fprintf(fp, "\n");
 511
 512         for (i = 0; i < n1; ++i)
 513         {
 514             fprintf(fp, "%g\t", tickx[i]);
 515             for (j = 0; j < n2; ++j)
 516             {
 517                 fprintf(fp, "%g\t", grid[i][j]);
 518             }
 519             fprintf(fp, "\n");
 520         }
 521         gmx_ffclose(fp);
 522     }
 523     else
 524     {
 525         fp = gmx_ffopen(ftp2fn(efXPM, NFILE, fnm), "w");
 526         write_xpm(fp, MAT_SPATIAL_X | MAT_SPATIAL_Y, buf, unit,
 527                   bRadial ? "axial (nm)" : label[c1], bRadial ? "r (nm)" : label[c2],
 528                   n1, n2, tickx, tickz, grid, dmin, maxgrid, rlo, rhi, &nlev);
 529         gmx_ffclose(fp);
 530     }
 531
 532     do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);
 533
 534     return 0;
 535 }