src/gromacs/gmxana/gmx_densmap.c

   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, 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 <math.h>
  40 #include <string.h>
  41
  42 #include "gromacs/commandline/pargs.h"
  43 #include "gromacs/fileio/matio.h"
  44 #include "gromacs/fileio/tpxio.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/fatalerror.h"
  56 #include "gromacs/utility/futil.h"
  57 #include "gromacs/utility/smalloc.h"
  58
  59 int gmx_densmap(int argc, char *argv[])
  60 {
  61     const char        *desc[] = {
  62         "[THISMODULE] computes 2D number-density maps.",
  63         "It can make planar and axial-radial density maps.",
  64         "The output [TT].xpm[tt] file can be visualized with for instance xv",
  65         "and can be converted to postscript with [TT]xpm2ps[tt].",
  66         "Optionally, output can be in text form to a [TT].dat[tt] file with [TT]-od[tt], instead of the usual [TT].xpm[tt] file with [TT]-o[tt].",
  67         "[PAR]",
  68         "The default analysis is a 2-D number-density map for a selected",
  69         "group of atoms in the x-y plane.",
  70         "The averaging direction can be changed with the option [TT]-aver[tt].",
  71         "When [TT]-xmin[tt] and/or [TT]-xmax[tt] are set only atoms that are",
  72         "within the limit(s) in the averaging direction are taken into account.",
  73         "The grid spacing is set with the option [TT]-bin[tt].",
  74         "When [TT]-n1[tt] or [TT]-n2[tt] is non-zero, the grid",
  75         "size is set by this option.",
  76         "Box size fluctuations are properly taken into account.",
  77         "[PAR]",
  78         "When options [TT]-amax[tt] and [TT]-rmax[tt] are set, an axial-radial",
  79         "number-density map is made. Three groups should be supplied, the centers",
  80         "of mass of the first two groups define the axis, the third defines the",
  81         "analysis group. The axial direction goes from -amax to +amax, where",
  82         "the center is defined as the midpoint between the centers of mass and",
  83         "the positive direction goes from the first to the second center of mass.",
  84         "The radial direction goes from 0 to rmax or from -rmax to +rmax",
  85         "when the [TT]-mirror[tt] option has been set.",
  86         "[PAR]",
  87         "The normalization of the output is set with the [TT]-unit[tt] option.",
  88         "The default produces a true number density. Unit [TT]nm-2[tt] leaves out",
  89         "the normalization for the averaging or the angular direction.",
  90         "Option [TT]count[tt] produces the count for each grid cell.",
  91         "When you do not want the scale in the output to go",
  92         "from zero to the maximum density, you can set the maximum",
  93         "with the option [TT]-dmax[tt]."
  94     };
  95     static int         n1      = 0, n2 = 0;
  96     static real        xmin    = -1, xmax = -1, bin = 0.02, dmin = 0, dmax = 0, amax = 0, rmax = 0;
  97     static gmx_bool    bMirror = FALSE, bSums = FALSE;
  98     static const char *eaver[] = { NULL, "z", "y", "x", NULL };
  99     static const char *eunit[] = { NULL, "nm-3", "nm-2", "count", NULL };
 100
 101     t_pargs            pa[] = {
 102         { "-bin", FALSE, etREAL, {&bin},
 103           "Grid size (nm)" },
 104         { "-aver", FALSE, etENUM, {eaver},
 105           "The direction to average over" },
 106         { "-xmin", FALSE, etREAL, {&xmin},
 107           "Minimum coordinate for averaging" },
 108         { "-xmax", FALSE, etREAL, {&xmax},
 109           "Maximum coordinate for averaging" },
 110         { "-n1", FALSE, etINT, {&n1},
 111           "Number of grid cells in the first direction" },
 112         { "-n2", FALSE, etINT, {&n2},
 113           "Number of grid cells in the second direction" },
 114         { "-amax", FALSE, etREAL, {&amax},
 115           "Maximum axial distance from the center"},
 116         { "-rmax", FALSE, etREAL, {&rmax},
 117           "Maximum radial distance" },
 118         { "-mirror", FALSE, etBOOL, {&bMirror},
 119           "Add the mirror image below the axial axis" },
 120         { "-sums", FALSE, etBOOL, {&bSums},
 121           "Print density sums (1D map) to stdout" },
 122         { "-unit", FALSE, etENUM, {eunit},
 123           "Unit for the output" },
 124         { "-dmin", FALSE, etREAL, {&dmin},
 125           "Minimum density in output"},
 126         { "-dmax", FALSE, etREAL, {&dmax},
 127           "Maximum density in output (0 means calculate it)"},
 128     };
 129     gmx_bool           bXmin, bXmax, bRadial;
 130     FILE              *fp;
 131     t_trxstatus       *status;
 132     t_topology         top;
 133     int                ePBC = -1;
 134     rvec              *x, xcom[2], direction, center, dx;
 135     matrix             box;
 136     real               t, m, mtot;
 137     t_pbc              pbc;
 138     int                cav = 0, c1 = 0, c2 = 0, natoms;
 139     char             **grpname, title[256], buf[STRLEN];
 140     const char        *unit;
 141     int                i, j, k, l, ngrps, anagrp, *gnx = NULL, nindex, nradial = 0, nfr, nmpower;
 142     atom_id          **ind = NULL, *index;
 143     real             **grid, maxgrid, m1, m2, box1, box2, *tickx, *tickz, invcellvol;
 144     real               invspa = 0, invspz = 0, axial, r, vol_old, vol, rowsum;
 145     int                nlev   = 51;
 146     t_rgb              rlo    = {1, 1, 1}, rhi = {0, 0, 0};
 147     output_env_t       oenv;
 148     const char        *label[] = { "x (nm)", "y (nm)", "z (nm)" };
 149     t_filenm           fnm[]   = {
 150         { efTRX, "-f",   NULL,       ffREAD },
 151         { efTPS, NULL,   NULL,       ffOPTRD },
 152         { efNDX, NULL,   NULL,       ffOPTRD },
 153         { efDAT, "-od",  "densmap",   ffOPTWR },
 154         { efXPM, "-o",   "densmap",   ffWRITE }
 155     };
 156 #define NFILE asize(fnm)
 157     int                npargs;
 158
 159     npargs = asize(pa);
 160
 161     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
 162                            NFILE, fnm, npargs, pa, asize(desc), desc, 0, NULL, &oenv))
 163     {
 164         return 0;
 165     }
 166
 167     bXmin   = opt2parg_bSet("-xmin", npargs, pa);
 168     bXmax   = opt2parg_bSet("-xmax", npargs, pa);
 169     bRadial = (amax > 0 || rmax > 0);
 170     if (bRadial)
 171     {
 172         if (amax <= 0 || rmax <= 0)
 173         {
 174             gmx_fatal(FARGS, "Both amax and rmax should be larger than zero");
 175         }
 176     }
 177
 178     if (strcmp(eunit[0], "nm-3") == 0)
 179     {
 180         nmpower = -3;
 181         unit    = "(nm^-3)";
 182     }
 183     else if (strcmp(eunit[0], "nm-2") == 0)
 184     {
 185         nmpower = -2;
 186         unit    = "(nm^-2)";
 187     }
 188     else
 189     {
 190         nmpower = 0;
 191         unit    = "count";
 192     }
 193
 194     if (ftp2bSet(efTPS, NFILE, fnm) || !ftp2bSet(efNDX, NFILE, fnm))
 195     {
 196         read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &x, NULL, box,
 197                       bRadial);
 198     }
 199     if (!bRadial)
 200     {
 201         ngrps = 1;
 202         fprintf(stderr, "\nSelect an analysis group\n");
 203     }
 204     else
 205     {
 206         ngrps = 3;
 207         fprintf(stderr,
 208                 "\nSelect two groups to define the axis and an analysis group\n");
 209     }
 210     snew(gnx, ngrps);
 211     snew(grpname, ngrps);
 212     snew(ind, ngrps);
 213     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, gnx, ind, grpname);
 214     anagrp = ngrps - 1;
 215     nindex = gnx[anagrp];
 216     index  = ind[anagrp];
 217     if (bRadial)
 218     {
 219         if ((gnx[0] > 1 || gnx[1] > 1) && !ftp2bSet(efTPS, NFILE, fnm))
 220         {
 221             gmx_fatal(FARGS, "No run input file was supplied (option -s), this is required for the center of mass calculation");
 222         }
 223     }
 224
 225     switch (eaver[0][0])
 226     {
 227         case 'x': cav = XX; c1 = YY; c2 = ZZ; break;
 228         case 'y': cav = YY; c1 = XX; c2 = ZZ; break;
 229         case 'z': cav = ZZ; c1 = XX; c2 = YY; break;
 230     }
 231
 232     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 233
 234     if (!bRadial)
 235     {
 236         if (n1 == 0)
 237         {
 238             n1 = (int)(box[c1][c1]/bin + 0.5);
 239         }
 240         if (n2 == 0)
 241         {
 242             n2 = (int)(box[c2][c2]/bin + 0.5);
 243         }
 244     }
 245     else
 246     {
 247         n1      = (int)(2*amax/bin + 0.5);
 248         nradial = (int)(rmax/bin + 0.5);
 249         invspa  = n1/(2*amax);
 250         invspz  = nradial/rmax;
 251         if (bMirror)
 252         {
 253             n2 = 2*nradial;
 254         }
 255         else
 256         {
 257             n2 = nradial;
 258         }
 259     }
 260
 261     snew(grid, n1);
 262     for (i = 0; i < n1; i++)
 263     {
 264         snew(grid[i], n2);
 265     }
 266
 267     box1 = 0;
 268     box2 = 0;
 269     nfr  = 0;
 270     do
 271     {
 272         if (!bRadial)
 273         {
 274             box1      += box[c1][c1];
 275             box2      += box[c2][c2];
 276             invcellvol = n1*n2;
 277             if (nmpower == -3)
 278             {
 279                 invcellvol /= det(box);
 280             }
 281             else if (nmpower == -2)
 282             {
 283                 invcellvol /= box[c1][c1]*box[c2][c2];
 284             }
 285             for (i = 0; i < nindex; i++)
 286             {
 287                 j = index[i];
 288                 if ((!bXmin || x[j][cav] >= xmin) &&
 289                     (!bXmax || x[j][cav] <= xmax))
 290                 {
 291                     m1 = x[j][c1]/box[c1][c1];
 292                     if (m1 >= 1)
 293                     {
 294                         m1 -= 1;
 295                     }
 296                     if (m1 < 0)
 297                     {
 298                         m1 += 1;
 299                     }
 300                     m2 = x[j][c2]/box[c2][c2];
 301                     if (m2 >= 1)
 302                     {
 303                         m2 -= 1;
 304                     }
 305                     if (m2 < 0)
 306                     {
 307                         m2 += 1;
 308                     }
 309                     grid[(int)(m1*n1)][(int)(m2*n2)] += invcellvol;
 310                 }
 311             }
 312         }
 313         else
 314         {
 315             set_pbc(&pbc, ePBC, box);
 316             for (i = 0; i < 2; i++)
 317             {
 318                 if (gnx[i] == 1)
 319                 {
 320                     /* One atom, just copy the coordinates */
 321                     copy_rvec(x[ind[i][0]], xcom[i]);
 322                 }
 323                 else
 324                 {
 325                     /* Calculate the center of mass */
 326                     clear_rvec(xcom[i]);
 327                     mtot = 0;
 328                     for (j = 0; j < gnx[i]; j++)
 329                     {
 330                         k = ind[i][j];
 331                         m = top.atoms.atom[k].m;
 332                         for (l = 0; l < DIM; l++)
 333                         {
 334                             xcom[i][l] += m*x[k][l];
 335                         }
 336                         mtot += m;
 337                     }
 338                     svmul(1/mtot, xcom[i], xcom[i]);
 339                 }
 340             }
 341             pbc_dx(&pbc, xcom[1], xcom[0], direction);
 342             for (i = 0; i < DIM; i++)
 343             {
 344                 center[i] = xcom[0][i] + 0.5*direction[i];
 345             }
 346             unitv(direction, direction);
 347             for (i = 0; i < nindex; i++)
 348             {
 349                 j = index[i];
 350                 pbc_dx(&pbc, x[j], center, dx);
 351                 axial = iprod(dx, direction);
 352                 r     = sqrt(norm2(dx) - axial*axial);
 353                 if (axial >= -amax && axial < amax && r < rmax)
 354                 {
 355                     if (bMirror)
 356                     {
 357                         r += rmax;
 358                     }
 359                     grid[(int)((axial + amax)*invspa)][(int)(r*invspz)] += 1;
 360                 }
 361             }
 362         }
 363         nfr++;
 364     }
 365     while (read_next_x(oenv, status, &t, x, box));
 366     close_trj(status);
 367
 368     /* normalize gridpoints */
 369     maxgrid = 0;
 370     if (!bRadial)
 371     {
 372         for (i = 0; i < n1; i++)
 373         {
 374             for (j = 0; j < n2; j++)
 375             {
 376                 grid[i][j] /= nfr;
 377                 if (grid[i][j] > maxgrid)
 378                 {
 379                     maxgrid = grid[i][j];
 380                 }
 381             }
 382         }
 383     }
 384     else
 385     {
 386         for (i = 0; i < n1; i++)
 387         {
 388             vol_old = 0;
 389             for (j = 0; j < nradial; j++)
 390             {
 391                 switch (nmpower)
 392                 {
 393                     case -3: vol = M_PI*(j+1)*(j+1)/(invspz*invspz*invspa); break;
 394                     case -2: vol =            (j+1)/(invspz*invspa);        break;
 395                     default: vol =             j+1;                         break;
 396                 }
 397                 if (bMirror)
 398                 {
 399                     k = j + nradial;
 400                 }
 401                 else
 402                 {
 403                     k = j;
 404                 }
 405                 grid[i][k] /= nfr*(vol - vol_old);
 406                 if (bMirror)
 407                 {
 408                     grid[i][nradial-1-j] = grid[i][k];
 409                 }
 410                 vol_old = vol;
 411                 if (grid[i][k] > maxgrid)
 412                 {
 413                     maxgrid = grid[i][k];
 414                 }
 415             }
 416         }
 417     }
 418     fprintf(stdout, "\n  The maximum density is %f %s\n", maxgrid, unit);
 419     if (dmax > 0)
 420     {
 421         maxgrid = dmax;
 422     }
 423
 424     snew(tickx, n1+1);
 425     snew(tickz, n2+1);
 426     if (!bRadial)
 427     {
 428         /* normalize box-axes */
 429         box1 /= nfr;
 430         box2 /= nfr;
 431         for (i = 0; i <= n1; i++)
 432         {
 433             tickx[i] = i*box1/n1;
 434         }
 435         for (i = 0; i <= n2; i++)
 436         {
 437             tickz[i] = i*box2/n2;
 438         }
 439     }
 440     else
 441     {
 442         for (i = 0; i <= n1; i++)
 443         {
 444             tickx[i] = i/invspa - amax;
 445         }
 446         if (bMirror)
 447         {
 448             for (i = 0; i <= n2; i++)
 449             {
 450                 tickz[i] = i/invspz - rmax;
 451             }
 452         }
 453         else
 454         {
 455             for (i = 0; i <= n2; i++)
 456             {
 457                 tickz[i] = i/invspz;
 458             }
 459         }
 460     }
 461
 462     if (bSums)
 463     {
 464         for (i = 0; i < n1; ++i)
 465         {
 466             fprintf(stdout, "Density sums:\n");
 467             rowsum = 0;
 468             for (j = 0; j < n2; ++j)
 469             {
 470                 rowsum += grid[i][j];
 471             }
 472             fprintf(stdout, "%g\t", rowsum);
 473         }
 474         fprintf(stdout, "\n");
 475     }
 476
 477     sprintf(buf, "%s number density", grpname[anagrp]);
 478     if (!bRadial && (bXmin || bXmax))
 479     {
 480         if (!bXmax)
 481         {
 482             sprintf(buf+strlen(buf), ", %c > %g nm", eaver[0][0], xmin);
 483         }
 484         else if (!bXmin)
 485         {
 486             sprintf(buf+strlen(buf), ", %c < %g nm", eaver[0][0], xmax);
 487         }
 488         else
 489         {
 490             sprintf(buf+strlen(buf), ", %c: %g - %g nm", eaver[0][0], xmin, xmax);
 491         }
 492     }
 493     if (ftp2bSet(efDAT, NFILE, fnm))
 494     {
 495         fp = gmx_ffopen(ftp2fn(efDAT, NFILE, fnm), "w");
 496         /*optional text form output:  first row is tickz; first col is tickx */
 497         fprintf(fp, "0\t");
 498         for (j = 0; j < n2; ++j)
 499         {
 500             fprintf(fp, "%g\t", tickz[j]);
 501         }
 502         fprintf(fp, "\n");
 503
 504         for (i = 0; i < n1; ++i)
 505         {
 506             fprintf(fp, "%g\t", tickx[i]);
 507             for (j = 0; j < n2; ++j)
 508             {
 509                 fprintf(fp, "%g\t", grid[i][j]);
 510             }
 511             fprintf(fp, "\n");
 512         }
 513         gmx_ffclose(fp);
 514     }
 515     else
 516     {
 517         fp = gmx_ffopen(ftp2fn(efXPM, NFILE, fnm), "w");
 518         write_xpm(fp, MAT_SPATIAL_X | MAT_SPATIAL_Y, buf, unit,
 519                   bRadial ? "axial (nm)" : label[c1], bRadial ? "r (nm)" : label[c2],
 520                   n1, n2, tickx, tickz, grid, dmin, maxgrid, rlo, rhi, &nlev);
 521         gmx_ffclose(fp);
 522     }
 523
 524     do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);
 525
 526     return 0;
 527 }