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 "config.h"
  40
  41 #include <math.h>
  42 #include <string.h>
  43
  44 #include "gromacs/commandline/pargs.h"
  45 #include "gromacs/legacyheaders/typedefs.h"
  46 #include "gromacs/utility/smalloc.h"
  47 #include "gromacs/legacyheaders/macros.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/pbcutil/pbc.h"
  50 #include "gromacs/utility/futil.h"
  51 #include "gromacs/topology/index.h"
  52 #include "gromacs/legacyheaders/txtdump.h"
  53 #include "gromacs/fileio/tpxio.h"
  54 #include "gromacs/fileio/trxio.h"
  55 #include "gstat.h"
  56 #include "gromacs/fileio/matio.h"
  57 #include "gromacs/legacyheaders/viewit.h"
  58 #include "gmx_ana.h"
  59
  60 #include "gromacs/utility/fatalerror.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 [TT].xpm[tt] 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 [TT].dat[tt] file with [TT]-od[tt], instead of the usual [TT].xpm[tt] file with [TT]-o[tt].",
  70         "[PAR]",
  71         "The default analysis is a 2-D number-density map for a selected",
  72         "group of atoms in the x-y plane.",
  73         "The averaging direction can be changed with the option [TT]-aver[tt].",
  74         "When [TT]-xmin[tt] and/or [TT]-xmax[tt] are set only atoms that are",
  75         "within the limit(s) in the averaging direction are taken into account.",
  76         "The grid spacing is set with the option [TT]-bin[tt].",
  77         "When [TT]-n1[tt] or [TT]-n2[tt] is non-zero, the grid",
  78         "size is set by this option.",
  79         "Box size fluctuations are properly taken into account.",
  80         "[PAR]",
  81         "When options [TT]-amax[tt] and [TT]-rmax[tt] are set, an axial-radial",
  82         "number-density map is made. Three groups should be supplied, the centers",
  83         "of mass of the first two groups define the axis, the third defines the",
  84         "analysis group. The axial direction goes from -amax to +amax, where",
  85         "the center is defined as the midpoint between the centers of mass and",
  86         "the positive direction goes from the first to the second center of mass.",
  87         "The radial direction goes from 0 to rmax or from -rmax to +rmax",
  88         "when the [TT]-mirror[tt] option has been set.",
  89         "[PAR]",
  90         "The normalization of the output is set with the [TT]-unit[tt] option.",
  91         "The default produces a true number density. Unit [TT]nm-2[tt] leaves out",
  92         "the normalization for the averaging or the angular direction.",
  93         "Option [TT]count[tt] produces the count for each grid cell.",
  94         "When you do not want the scale in the output to go",
  95         "from zero to the maximum density, you can set the maximum",
  96         "with the option [TT]-dmax[tt]."
  97     };
  98     static int         n1      = 0, n2 = 0;
  99     static real        xmin    = -1, xmax = -1, bin = 0.02, dmin = 0, dmax = 0, amax = 0, rmax = 0;
 100     static gmx_bool    bMirror = FALSE, bSums = FALSE;
 101     static const char *eaver[] = { NULL, "z", "y", "x", NULL };
 102     static const char *eunit[] = { NULL, "nm-3", "nm-2", "count", NULL };
 103
 104     t_pargs            pa[] = {
 105         { "-bin", FALSE, etREAL, {&bin},
 106           "Grid size (nm)" },
 107         { "-aver", FALSE, etENUM, {eaver},
 108           "The direction to average over" },
 109         { "-xmin", FALSE, etREAL, {&xmin},
 110           "Minimum coordinate for averaging" },
 111         { "-xmax", FALSE, etREAL, {&xmax},
 112           "Maximum coordinate for averaging" },
 113         { "-n1", FALSE, etINT, {&n1},
 114           "Number of grid cells in the first direction" },
 115         { "-n2", FALSE, etINT, {&n2},
 116           "Number of grid cells in the second direction" },
 117         { "-amax", FALSE, etREAL, {&amax},
 118           "Maximum axial distance from the center"},
 119         { "-rmax", FALSE, etREAL, {&rmax},
 120           "Maximum radial distance" },
 121         { "-mirror", FALSE, etBOOL, {&bMirror},
 122           "Add the mirror image below the axial axis" },
 123         { "-sums", FALSE, etBOOL, {&bSums},
 124           "Print density sums (1D map) to stdout" },
 125         { "-unit", FALSE, etENUM, {eunit},
 126           "Unit for the output" },
 127         { "-dmin", FALSE, etREAL, {&dmin},
 128           "Minimum density in output"},
 129         { "-dmax", FALSE, etREAL, {&dmax},
 130           "Maximum density in output (0 means calculate it)"},
 131     };
 132     gmx_bool           bXmin, bXmax, bRadial;
 133     FILE              *fp;
 134     t_trxstatus       *status;
 135     t_topology         top;
 136     int                ePBC = -1;
 137     rvec              *x, xcom[2], direction, center, dx;
 138     matrix             box;
 139     real               t, m, mtot;
 140     t_pbc              pbc;
 141     int                cav = 0, c1 = 0, c2 = 0, natoms;
 142     char             **grpname, title[256], buf[STRLEN];
 143     const char        *unit;
 144     int                i, j, k, l, ngrps, anagrp, *gnx = NULL, nindex, nradial = 0, nfr, nmpower;
 145     atom_id          **ind = NULL, *index;
 146     real             **grid, maxgrid, m1, m2, box1, box2, *tickx, *tickz, invcellvol;
 147     real               invspa = 0, invspz = 0, axial, r, vol_old, vol, rowsum;
 148     int                nlev   = 51;
 149     t_rgb              rlo    = {1, 1, 1}, rhi = {0, 0, 0};
 150     output_env_t       oenv;
 151     const char        *label[] = { "x (nm)", "y (nm)", "z (nm)" };
 152     t_filenm           fnm[]   = {
 153         { efTRX, "-f",   NULL,       ffREAD },
 154         { efTPS, NULL,   NULL,       ffOPTRD },
 155         { efNDX, NULL,   NULL,       ffOPTRD },
 156         { efDAT, "-od",  "densmap",   ffOPTWR },
 157         { efXPM, "-o",   "densmap",   ffWRITE }
 158     };
 159 #define NFILE asize(fnm)
 160     int                npargs;
 161
 162     npargs = asize(pa);
 163
 164     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
 165                            NFILE, fnm, npargs, pa, asize(desc), desc, 0, NULL, &oenv))
 166     {
 167         return 0;
 168     }
 169
 170     bXmin   = opt2parg_bSet("-xmin", npargs, pa);
 171     bXmax   = opt2parg_bSet("-xmax", npargs, pa);
 172     bRadial = (amax > 0 || rmax > 0);
 173     if (bRadial)
 174     {
 175         if (amax <= 0 || rmax <= 0)
 176         {
 177             gmx_fatal(FARGS, "Both amax and rmax should be larger than zero");
 178         }
 179     }
 180
 181     if (strcmp(eunit[0], "nm-3") == 0)
 182     {
 183         nmpower = -3;
 184         unit    = "(nm^-3)";
 185     }
 186     else if (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), title, &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     switch (eaver[0][0])
 229     {
 230         case 'x': cav = XX; c1 = YY; c2 = ZZ; break;
 231         case 'y': cav = YY; c1 = XX; c2 = ZZ; break;
 232         case 'z': cav = ZZ; c1 = XX; c2 = YY; break;
 233     }
 234
 235     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 236
 237     if (!bRadial)
 238     {
 239         if (n1 == 0)
 240         {
 241             n1 = (int)(box[c1][c1]/bin + 0.5);
 242         }
 243         if (n2 == 0)
 244         {
 245             n2 = (int)(box[c2][c2]/bin + 0.5);
 246         }
 247     }
 248     else
 249     {
 250         n1      = (int)(2*amax/bin + 0.5);
 251         nradial = (int)(rmax/bin + 0.5);
 252         invspa  = n1/(2*amax);
 253         invspz  = nradial/rmax;
 254         if (bMirror)
 255         {
 256             n2 = 2*nradial;
 257         }
 258         else
 259         {
 260             n2 = nradial;
 261         }
 262     }
 263
 264     snew(grid, n1);
 265     for (i = 0; i < n1; i++)
 266     {
 267         snew(grid[i], n2);
 268     }
 269
 270     box1 = 0;
 271     box2 = 0;
 272     nfr  = 0;
 273     do
 274     {
 275         if (!bRadial)
 276         {
 277             box1      += box[c1][c1];
 278             box2      += box[c2][c2];
 279             invcellvol = n1*n2;
 280             if (nmpower == -3)
 281             {
 282                 invcellvol /= det(box);
 283             }
 284             else if (nmpower == -2)
 285             {
 286                 invcellvol /= box[c1][c1]*box[c2][c2];
 287             }
 288             for (i = 0; i < nindex; i++)
 289             {
 290                 j = index[i];
 291                 if ((!bXmin || x[j][cav] >= xmin) &&
 292                     (!bXmax || x[j][cav] <= xmax))
 293                 {
 294                     m1 = x[j][c1]/box[c1][c1];
 295                     if (m1 >= 1)
 296                     {
 297                         m1 -= 1;
 298                     }
 299                     if (m1 < 0)
 300                     {
 301                         m1 += 1;
 302                     }
 303                     m2 = x[j][c2]/box[c2][c2];
 304                     if (m2 >= 1)
 305                     {
 306                         m2 -= 1;
 307                     }
 308                     if (m2 < 0)
 309                     {
 310                         m2 += 1;
 311                     }
 312                     grid[(int)(m1*n1)][(int)(m2*n2)] += invcellvol;
 313                 }
 314             }
 315         }
 316         else
 317         {
 318             set_pbc(&pbc, ePBC, box);
 319             for (i = 0; i < 2; i++)
 320             {
 321                 if (gnx[i] == 1)
 322                 {
 323                     /* One atom, just copy the coordinates */
 324                     copy_rvec(x[ind[i][0]], xcom[i]);
 325                 }
 326                 else
 327                 {
 328                     /* Calculate the center of mass */
 329                     clear_rvec(xcom[i]);
 330                     mtot = 0;
 331                     for (j = 0; j < gnx[i]; j++)
 332                     {
 333                         k = ind[i][j];
 334                         m = top.atoms.atom[k].m;
 335                         for (l = 0; l < DIM; l++)
 336                         {
 337                             xcom[i][l] += m*x[k][l];
 338                         }
 339                         mtot += m;
 340                     }
 341                     svmul(1/mtot, xcom[i], xcom[i]);
 342                 }
 343             }
 344             pbc_dx(&pbc, xcom[1], xcom[0], direction);
 345             for (i = 0; i < DIM; i++)
 346             {
 347                 center[i] = xcom[0][i] + 0.5*direction[i];
 348             }
 349             unitv(direction, direction);
 350             for (i = 0; i < nindex; i++)
 351             {
 352                 j = index[i];
 353                 pbc_dx(&pbc, x[j], center, dx);
 354                 axial = iprod(dx, direction);
 355                 r     = sqrt(norm2(dx) - axial*axial);
 356                 if (axial >= -amax && axial < amax && r < rmax)
 357                 {
 358                     if (bMirror)
 359                     {
 360                         r += rmax;
 361                     }
 362                     grid[(int)((axial + amax)*invspa)][(int)(r*invspz)] += 1;
 363                 }
 364             }
 365         }
 366         nfr++;
 367     }
 368     while (read_next_x(oenv, status, &t, x, box));
 369     close_trj(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+strlen(buf), ", %c > %g nm", eaver[0][0], xmin);
 486         }
 487         else if (!bXmin)
 488         {
 489             sprintf(buf+strlen(buf), ", %c < %g nm", eaver[0][0], xmax);
 490         }
 491         else
 492         {
 493             sprintf(buf+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,
 522                   bRadial ? "axial (nm)" : label[c1], bRadial ? "r (nm)" : label[c2],
 523                   n1, n2, tickx, tickz, grid, dmin, maxgrid, rlo, rhi, &nlev);
 524         gmx_ffclose(fp);
 525     }
 526
 527     do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);
 528
 529     return 0;
 530 }