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/legacyheaders/typedefs.h"
  44 #include "gromacs/utility/smalloc.h"
  45 #include "gromacs/legacyheaders/macros.h"
  46 #include "gromacs/math/vec.h"
  47 #include "gromacs/pbcutil/pbc.h"
  48 #include "gromacs/utility/futil.h"
  49 #include "gromacs/topology/index.h"
  50 #include "gromacs/legacyheaders/txtdump.h"
  51 #include "gromacs/fileio/tpxio.h"
  52 #include "gromacs/fileio/trxio.h"
  53 #include "gstat.h"
  54 #include "gromacs/fileio/matio.h"
  55 #include "gromacs/legacyheaders/viewit.h"
  56 #include "gmx_ana.h"
  57
  58 #include "gromacs/utility/fatalerror.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 [TT].xpm[tt] 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 [TT].dat[tt] file with [TT]-od[tt], instead of the usual [TT].xpm[tt] 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, natoms;
 140     char             **grpname, title[256], buf[STRLEN];
 141     const char        *unit;
 142     int                i, j, k, l, ngrps, anagrp, *gnx = NULL, nindex, nradial = 0, nfr, nmpower;
 143     atom_id          **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     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     if (strcmp(eunit[0], "nm-3") == 0)
 180     {
 181         nmpower = -3;
 182         unit    = "(nm^-3)";
 183     }
 184     else if (strcmp(eunit[0], "nm-2") == 0)
 185     {
 186         nmpower = -2;
 187         unit    = "(nm^-2)";
 188     }
 189     else
 190     {
 191         nmpower = 0;
 192         unit    = "count";
 193     }
 194
 195     if (ftp2bSet(efTPS, NFILE, fnm) || !ftp2bSet(efNDX, NFILE, fnm))
 196     {
 197         read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &x, NULL, box,
 198                       bRadial);
 199     }
 200     if (!bRadial)
 201     {
 202         ngrps = 1;
 203         fprintf(stderr, "\nSelect an analysis group\n");
 204     }
 205     else
 206     {
 207         ngrps = 3;
 208         fprintf(stderr,
 209                 "\nSelect two groups to define the axis and an analysis group\n");
 210     }
 211     snew(gnx, ngrps);
 212     snew(grpname, ngrps);
 213     snew(ind, ngrps);
 214     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, gnx, ind, grpname);
 215     anagrp = ngrps - 1;
 216     nindex = gnx[anagrp];
 217     index  = ind[anagrp];
 218     if (bRadial)
 219     {
 220         if ((gnx[0] > 1 || gnx[1] > 1) && !ftp2bSet(efTPS, NFILE, fnm))
 221         {
 222             gmx_fatal(FARGS, "No run input file was supplied (option -s), this is required for the center of mass calculation");
 223         }
 224     }
 225
 226     switch (eaver[0][0])
 227     {
 228         case 'x': cav = XX; c1 = YY; c2 = ZZ; break;
 229         case 'y': cav = YY; c1 = XX; c2 = ZZ; break;
 230         case 'z': cav = ZZ; c1 = XX; c2 = YY; break;
 231     }
 232
 233     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 234
 235     if (!bRadial)
 236     {
 237         if (n1 == 0)
 238         {
 239             n1 = (int)(box[c1][c1]/bin + 0.5);
 240         }
 241         if (n2 == 0)
 242         {
 243             n2 = (int)(box[c2][c2]/bin + 0.5);
 244         }
 245     }
 246     else
 247     {
 248         n1      = (int)(2*amax/bin + 0.5);
 249         nradial = (int)(rmax/bin + 0.5);
 250         invspa  = n1/(2*amax);
 251         invspz  = nradial/rmax;
 252         if (bMirror)
 253         {
 254             n2 = 2*nradial;
 255         }
 256         else
 257         {
 258             n2 = nradial;
 259         }
 260     }
 261
 262     snew(grid, n1);
 263     for (i = 0; i < n1; i++)
 264     {
 265         snew(grid[i], n2);
 266     }
 267
 268     box1 = 0;
 269     box2 = 0;
 270     nfr  = 0;
 271     do
 272     {
 273         if (!bRadial)
 274         {
 275             box1      += box[c1][c1];
 276             box2      += box[c2][c2];
 277             invcellvol = n1*n2;
 278             if (nmpower == -3)
 279             {
 280                 invcellvol /= det(box);
 281             }
 282             else if (nmpower == -2)
 283             {
 284                 invcellvol /= box[c1][c1]*box[c2][c2];
 285             }
 286             for (i = 0; i < nindex; i++)
 287             {
 288                 j = index[i];
 289                 if ((!bXmin || x[j][cav] >= xmin) &&
 290                     (!bXmax || x[j][cav] <= xmax))
 291                 {
 292                     m1 = x[j][c1]/box[c1][c1];
 293                     if (m1 >= 1)
 294                     {
 295                         m1 -= 1;
 296                     }
 297                     if (m1 < 0)
 298                     {
 299                         m1 += 1;
 300                     }
 301                     m2 = x[j][c2]/box[c2][c2];
 302                     if (m2 >= 1)
 303                     {
 304                         m2 -= 1;
 305                     }
 306                     if (m2 < 0)
 307                     {
 308                         m2 += 1;
 309                     }
 310                     grid[(int)(m1*n1)][(int)(m2*n2)] += invcellvol;
 311                 }
 312             }
 313         }
 314         else
 315         {
 316             set_pbc(&pbc, ePBC, box);
 317             for (i = 0; i < 2; i++)
 318             {
 319                 if (gnx[i] == 1)
 320                 {
 321                     /* One atom, just copy the coordinates */
 322                     copy_rvec(x[ind[i][0]], xcom[i]);
 323                 }
 324                 else
 325                 {
 326                     /* Calculate the center of mass */
 327                     clear_rvec(xcom[i]);
 328                     mtot = 0;
 329                     for (j = 0; j < gnx[i]; j++)
 330                     {
 331                         k = ind[i][j];
 332                         m = top.atoms.atom[k].m;
 333                         for (l = 0; l < DIM; l++)
 334                         {
 335                             xcom[i][l] += m*x[k][l];
 336                         }
 337                         mtot += m;
 338                     }
 339                     svmul(1/mtot, xcom[i], xcom[i]);
 340                 }
 341             }
 342             pbc_dx(&pbc, xcom[1], xcom[0], direction);
 343             for (i = 0; i < DIM; i++)
 344             {
 345                 center[i] = xcom[0][i] + 0.5*direction[i];
 346             }
 347             unitv(direction, direction);
 348             for (i = 0; i < nindex; i++)
 349             {
 350                 j = index[i];
 351                 pbc_dx(&pbc, x[j], center, dx);
 352                 axial = iprod(dx, direction);
 353                 r     = sqrt(norm2(dx) - axial*axial);
 354                 if (axial >= -amax && axial < amax && r < rmax)
 355                 {
 356                     if (bMirror)
 357                     {
 358                         r += rmax;
 359                     }
 360                     grid[(int)((axial + amax)*invspa)][(int)(r*invspz)] += 1;
 361                 }
 362             }
 363         }
 364         nfr++;
 365     }
 366     while (read_next_x(oenv, status, &t, x, box));
 367     close_trj(status);
 368
 369     /* normalize gridpoints */
 370     maxgrid = 0;
 371     if (!bRadial)
 372     {
 373         for (i = 0; i < n1; i++)
 374         {
 375             for (j = 0; j < n2; j++)
 376             {
 377                 grid[i][j] /= nfr;
 378                 if (grid[i][j] > maxgrid)
 379                 {
 380                     maxgrid = grid[i][j];
 381                 }
 382             }
 383         }
 384     }
 385     else
 386     {
 387         for (i = 0; i < n1; i++)
 388         {
 389             vol_old = 0;
 390             for (j = 0; j < nradial; j++)
 391             {
 392                 switch (nmpower)
 393                 {
 394                     case -3: vol = M_PI*(j+1)*(j+1)/(invspz*invspz*invspa); break;
 395                     case -2: vol =            (j+1)/(invspz*invspa);        break;
 396                     default: vol =             j+1;                         break;
 397                 }
 398                 if (bMirror)
 399                 {
 400                     k = j + nradial;
 401                 }
 402                 else
 403                 {
 404                     k = j;
 405                 }
 406                 grid[i][k] /= nfr*(vol - vol_old);
 407                 if (bMirror)
 408                 {
 409                     grid[i][nradial-1-j] = grid[i][k];
 410                 }
 411                 vol_old = vol;
 412                 if (grid[i][k] > maxgrid)
 413                 {
 414                     maxgrid = grid[i][k];
 415                 }
 416             }
 417         }
 418     }
 419     fprintf(stdout, "\n  The maximum density is %f %s\n", maxgrid, unit);
 420     if (dmax > 0)
 421     {
 422         maxgrid = dmax;
 423     }
 424
 425     snew(tickx, n1+1);
 426     snew(tickz, n2+1);
 427     if (!bRadial)
 428     {
 429         /* normalize box-axes */
 430         box1 /= nfr;
 431         box2 /= nfr;
 432         for (i = 0; i <= n1; i++)
 433         {
 434             tickx[i] = i*box1/n1;
 435         }
 436         for (i = 0; i <= n2; i++)
 437         {
 438             tickz[i] = i*box2/n2;
 439         }
 440     }
 441     else
 442     {
 443         for (i = 0; i <= n1; i++)
 444         {
 445             tickx[i] = i/invspa - amax;
 446         }
 447         if (bMirror)
 448         {
 449             for (i = 0; i <= n2; i++)
 450             {
 451                 tickz[i] = i/invspz - rmax;
 452             }
 453         }
 454         else
 455         {
 456             for (i = 0; i <= n2; i++)
 457             {
 458                 tickz[i] = i/invspz;
 459             }
 460         }
 461     }
 462
 463     if (bSums)
 464     {
 465         for (i = 0; i < n1; ++i)
 466         {
 467             fprintf(stdout, "Density sums:\n");
 468             rowsum = 0;
 469             for (j = 0; j < n2; ++j)
 470             {
 471                 rowsum += grid[i][j];
 472             }
 473             fprintf(stdout, "%g\t", rowsum);
 474         }
 475         fprintf(stdout, "\n");
 476     }
 477
 478     sprintf(buf, "%s number density", grpname[anagrp]);
 479     if (!bRadial && (bXmin || bXmax))
 480     {
 481         if (!bXmax)
 482         {
 483             sprintf(buf+strlen(buf), ", %c > %g nm", eaver[0][0], xmin);
 484         }
 485         else if (!bXmin)
 486         {
 487             sprintf(buf+strlen(buf), ", %c < %g nm", eaver[0][0], xmax);
 488         }
 489         else
 490         {
 491             sprintf(buf+strlen(buf), ", %c: %g - %g nm", eaver[0][0], xmin, xmax);
 492         }
 493     }
 494     if (ftp2bSet(efDAT, NFILE, fnm))
 495     {
 496         fp = gmx_ffopen(ftp2fn(efDAT, NFILE, fnm), "w");
 497         /*optional text form output:  first row is tickz; first col is tickx */
 498         fprintf(fp, "0\t");
 499         for (j = 0; j < n2; ++j)
 500         {
 501             fprintf(fp, "%g\t", tickz[j]);
 502         }
 503         fprintf(fp, "\n");
 504
 505         for (i = 0; i < n1; ++i)
 506         {
 507             fprintf(fp, "%g\t", tickx[i]);
 508             for (j = 0; j < n2; ++j)
 509             {
 510                 fprintf(fp, "%g\t", grid[i][j]);
 511             }
 512             fprintf(fp, "\n");
 513         }
 514         gmx_ffclose(fp);
 515     }
 516     else
 517     {
 518         fp = gmx_ffopen(ftp2fn(efXPM, NFILE, fnm), "w");
 519         write_xpm(fp, MAT_SPATIAL_X | MAT_SPATIAL_Y, buf, unit,
 520                   bRadial ? "axial (nm)" : label[c1], bRadial ? "r (nm)" : label[c2],
 521                   n1, n2, tickx, tickz, grid, dmin, maxgrid, rlo, rhi, &nlev);
 522         gmx_ffclose(fp);
 523     }
 524
 525     do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);
 526
 527     return 0;
 528 }