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