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