2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
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"
59 int gmx_densmap(int argc, char *argv[])
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].",
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.",
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.",
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]."
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 };
102 { "-bin", FALSE, etREAL, {&bin},
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)"},
129 gmx_bool bXmin, bXmax, bRadial;
134 rvec *x, xcom[2], direction, center, dx;
138 int cav = 0, c1 = 0, c2 = 0, natoms;
139 char **grpname, title[256], buf[STRLEN];
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;
146 t_rgb rlo = {1, 1, 1}, rhi = {0, 0, 0};
148 const char *label[] = { "x (nm)", "y (nm)", "z (nm)" };
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 }
156 #define NFILE asize(fnm)
161 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
162 NFILE, fnm, npargs, pa, asize(desc), desc, 0, NULL, &oenv))
167 bXmin = opt2parg_bSet("-xmin", npargs, pa);
168 bXmax = opt2parg_bSet("-xmax", npargs, pa);
169 bRadial = (amax > 0 || rmax > 0);
172 if (amax <= 0 || rmax <= 0)
174 gmx_fatal(FARGS, "Both amax and rmax should be larger than zero");
178 if (strcmp(eunit[0], "nm-3") == 0)
183 else if (strcmp(eunit[0], "nm-2") == 0)
194 if (ftp2bSet(efTPS, NFILE, fnm) || !ftp2bSet(efNDX, NFILE, fnm))
196 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &x, NULL, box,
202 fprintf(stderr, "\nSelect an analysis group\n");
208 "\nSelect two groups to define the axis and an analysis group\n");
211 snew(grpname, ngrps);
213 get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, gnx, ind, grpname);
215 nindex = gnx[anagrp];
219 if ((gnx[0] > 1 || gnx[1] > 1) && !ftp2bSet(efTPS, NFILE, fnm))
221 gmx_fatal(FARGS, "No run input file was supplied (option -s), this is required for the center of mass calculation");
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;
232 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
238 n1 = (int)(box[c1][c1]/bin + 0.5);
242 n2 = (int)(box[c2][c2]/bin + 0.5);
247 n1 = (int)(2*amax/bin + 0.5);
248 nradial = (int)(rmax/bin + 0.5);
249 invspa = n1/(2*amax);
250 invspz = nradial/rmax;
262 for (i = 0; i < n1; i++)
279 invcellvol /= det(box);
281 else if (nmpower == -2)
283 invcellvol /= box[c1][c1]*box[c2][c2];
285 for (i = 0; i < nindex; i++)
288 if ((!bXmin || x[j][cav] >= xmin) &&
289 (!bXmax || x[j][cav] <= xmax))
291 m1 = x[j][c1]/box[c1][c1];
300 m2 = x[j][c2]/box[c2][c2];
309 grid[(int)(m1*n1)][(int)(m2*n2)] += invcellvol;
315 set_pbc(&pbc, ePBC, box);
316 for (i = 0; i < 2; i++)
320 /* One atom, just copy the coordinates */
321 copy_rvec(x[ind[i][0]], xcom[i]);
325 /* Calculate the center of mass */
328 for (j = 0; j < gnx[i]; j++)
331 m = top.atoms.atom[k].m;
332 for (l = 0; l < DIM; l++)
334 xcom[i][l] += m*x[k][l];
338 svmul(1/mtot, xcom[i], xcom[i]);
341 pbc_dx(&pbc, xcom[1], xcom[0], direction);
342 for (i = 0; i < DIM; i++)
344 center[i] = xcom[0][i] + 0.5*direction[i];
346 unitv(direction, direction);
347 for (i = 0; i < nindex; 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)
359 grid[(int)((axial + amax)*invspa)][(int)(r*invspz)] += 1;
365 while (read_next_x(oenv, status, &t, x, box));
368 /* normalize gridpoints */
372 for (i = 0; i < n1; i++)
374 for (j = 0; j < n2; j++)
377 if (grid[i][j] > maxgrid)
379 maxgrid = grid[i][j];
386 for (i = 0; i < n1; i++)
389 for (j = 0; j < nradial; j++)
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;
405 grid[i][k] /= nfr*(vol - vol_old);
408 grid[i][nradial-1-j] = grid[i][k];
411 if (grid[i][k] > maxgrid)
413 maxgrid = grid[i][k];
418 fprintf(stdout, "\n The maximum density is %f %s\n", maxgrid, unit);
428 /* normalize box-axes */
431 for (i = 0; i <= n1; i++)
433 tickx[i] = i*box1/n1;
435 for (i = 0; i <= n2; i++)
437 tickz[i] = i*box2/n2;
442 for (i = 0; i <= n1; i++)
444 tickx[i] = i/invspa - amax;
448 for (i = 0; i <= n2; i++)
450 tickz[i] = i/invspz - rmax;
455 for (i = 0; i <= n2; i++)
464 for (i = 0; i < n1; ++i)
466 fprintf(stdout, "Density sums:\n");
468 for (j = 0; j < n2; ++j)
470 rowsum += grid[i][j];
472 fprintf(stdout, "%g\t", rowsum);
474 fprintf(stdout, "\n");
477 sprintf(buf, "%s number density", grpname[anagrp]);
478 if (!bRadial && (bXmin || bXmax))
482 sprintf(buf+strlen(buf), ", %c > %g nm", eaver[0][0], xmin);
486 sprintf(buf+strlen(buf), ", %c < %g nm", eaver[0][0], xmax);
490 sprintf(buf+strlen(buf), ", %c: %g - %g nm", eaver[0][0], xmin, xmax);
493 if (ftp2bSet(efDAT, NFILE, fnm))
495 fp = gmx_ffopen(ftp2fn(efDAT, NFILE, fnm), "w");
496 /*optional text form output: first row is tickz; first col is tickx */
498 for (j = 0; j < n2; ++j)
500 fprintf(fp, "%g\t", tickz[j]);
504 for (i = 0; i < n1; ++i)
506 fprintf(fp, "%g\t", tickx[i]);
507 for (j = 0; j < n2; ++j)
509 fprintf(fp, "%g\t", grid[i][j]);
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);
524 do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);