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