3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
49 #include "gromacs/fileio/futil.h"
57 #include "gromacs/fileio/tpxio.h"
58 #include "gromacs/fileio/trxio.h"
60 #include "gromacs/fileio/matio.h"
65 int gmx_densmap(int argc, char *argv[])
67 const char *desc[] = {
68 "[THISMODULE] computes 2D number-density maps.",
69 "It can make planar and axial-radial density maps.",
70 "The output [TT].xpm[tt] file can be visualized with for instance xv",
71 "and can be converted to postscript with [TT]xpm2ps[tt].",
72 "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].",
74 "The default analysis is a 2-D number-density map for a selected",
75 "group of atoms in the x-y plane.",
76 "The averaging direction can be changed with the option [TT]-aver[tt].",
77 "When [TT]-xmin[tt] and/or [TT]-xmax[tt] are set only atoms that are",
78 "within the limit(s) in the averaging direction are taken into account.",
79 "The grid spacing is set with the option [TT]-bin[tt].",
80 "When [TT]-n1[tt] or [TT]-n2[tt] is non-zero, the grid",
81 "size is set by this option.",
82 "Box size fluctuations are properly taken into account.",
84 "When options [TT]-amax[tt] and [TT]-rmax[tt] are set, an axial-radial",
85 "number-density map is made. Three groups should be supplied, the centers",
86 "of mass of the first two groups define the axis, the third defines the",
87 "analysis group. The axial direction goes from -amax to +amax, where",
88 "the center is defined as the midpoint between the centers of mass and",
89 "the positive direction goes from the first to the second center of mass.",
90 "The radial direction goes from 0 to rmax or from -rmax to +rmax",
91 "when the [TT]-mirror[tt] option has been set.",
93 "The normalization of the output is set with the [TT]-unit[tt] option.",
94 "The default produces a true number density. Unit [TT]nm-2[tt] leaves out",
95 "the normalization for the averaging or the angular direction.",
96 "Option [TT]count[tt] produces the count for each grid cell.",
97 "When you do not want the scale in the output to go",
98 "from zero to the maximum density, you can set the maximum",
99 "with the option [TT]-dmax[tt]."
101 static int n1 = 0, n2 = 0;
102 static real xmin = -1, xmax = -1, bin = 0.02, dmin = 0, dmax = 0, amax = 0, rmax = 0;
103 static gmx_bool bMirror = FALSE, bSums = FALSE;
104 static const char *eaver[] = { NULL, "z", "y", "x", NULL };
105 static const char *eunit[] = { NULL, "nm-3", "nm-2", "count", NULL };
108 { "-bin", FALSE, etREAL, {&bin},
110 { "-aver", FALSE, etENUM, {eaver},
111 "The direction to average over" },
112 { "-xmin", FALSE, etREAL, {&xmin},
113 "Minimum coordinate for averaging" },
114 { "-xmax", FALSE, etREAL, {&xmax},
115 "Maximum coordinate for averaging" },
116 { "-n1", FALSE, etINT, {&n1},
117 "Number of grid cells in the first direction" },
118 { "-n2", FALSE, etINT, {&n2},
119 "Number of grid cells in the second direction" },
120 { "-amax", FALSE, etREAL, {&amax},
121 "Maximum axial distance from the center"},
122 { "-rmax", FALSE, etREAL, {&rmax},
123 "Maximum radial distance" },
124 { "-mirror", FALSE, etBOOL, {&bMirror},
125 "Add the mirror image below the axial axis" },
126 { "-sums", FALSE, etBOOL, {&bSums},
127 "Print density sums (1D map) to stdout" },
128 { "-unit", FALSE, etENUM, {eunit},
129 "Unit for the output" },
130 { "-dmin", FALSE, etREAL, {&dmin},
131 "Minimum density in output"},
132 { "-dmax", FALSE, etREAL, {&dmax},
133 "Maximum density in output (0 means calculate it)"},
135 gmx_bool bXmin, bXmax, bRadial;
140 rvec *x, xcom[2], direction, center, dx;
144 int cav = 0, c1 = 0, c2 = 0, natoms;
145 char **grpname, title[256], buf[STRLEN];
147 int i, j, k, l, ngrps, anagrp, *gnx = NULL, nindex, nradial = 0, nfr, nmpower;
148 atom_id **ind = NULL, *index;
149 real **grid, maxgrid, m1, m2, box1, box2, *tickx, *tickz, invcellvol;
150 real invspa = 0, invspz = 0, axial, r, vol_old, vol, rowsum;
152 t_rgb rlo = {1, 1, 1}, rhi = {0, 0, 0};
154 const char *label[] = { "x (nm)", "y (nm)", "z (nm)" };
156 { efTRX, "-f", NULL, ffREAD },
157 { efTPS, NULL, NULL, ffOPTRD },
158 { efNDX, NULL, NULL, ffOPTRD },
159 { efDAT, "-od", "densmap", ffOPTWR },
160 { efXPM, "-o", "densmap", ffWRITE }
162 #define NFILE asize(fnm)
167 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
168 NFILE, fnm, npargs, pa, asize(desc), desc, 0, NULL, &oenv))
173 bXmin = opt2parg_bSet("-xmin", npargs, pa);
174 bXmax = opt2parg_bSet("-xmax", npargs, pa);
175 bRadial = (amax > 0 || rmax > 0);
178 if (amax <= 0 || rmax <= 0)
180 gmx_fatal(FARGS, "Both amax and rmax should be larger than zero");
184 if (strcmp(eunit[0], "nm-3") == 0)
189 else if (strcmp(eunit[0], "nm-2") == 0)
200 if (ftp2bSet(efTPS, NFILE, fnm) || !ftp2bSet(efNDX, NFILE, fnm))
202 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &x, NULL, box,
208 fprintf(stderr, "\nSelect an analysis group\n");
214 "\nSelect two groups to define the axis and an analysis group\n");
217 snew(grpname, ngrps);
219 get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, gnx, ind, grpname);
221 nindex = gnx[anagrp];
225 if ((gnx[0] > 1 || gnx[1] > 1) && !ftp2bSet(efTPS, NFILE, fnm))
227 gmx_fatal(FARGS, "No run input file was supplied (option -s), this is required for the center of mass calculation");
233 case 'x': cav = XX; c1 = YY; c2 = ZZ; break;
234 case 'y': cav = YY; c1 = XX; c2 = ZZ; break;
235 case 'z': cav = ZZ; c1 = XX; c2 = YY; break;
238 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
244 n1 = (int)(box[c1][c1]/bin + 0.5);
248 n2 = (int)(box[c2][c2]/bin + 0.5);
253 n1 = (int)(2*amax/bin + 0.5);
254 nradial = (int)(rmax/bin + 0.5);
255 invspa = n1/(2*amax);
256 invspz = nradial/rmax;
268 for (i = 0; i < n1; i++)
285 invcellvol /= det(box);
287 else if (nmpower == -2)
289 invcellvol /= box[c1][c1]*box[c2][c2];
291 for (i = 0; i < nindex; i++)
294 if ((!bXmin || x[j][cav] >= xmin) &&
295 (!bXmax || x[j][cav] <= xmax))
297 m1 = x[j][c1]/box[c1][c1];
306 m2 = x[j][c2]/box[c2][c2];
315 grid[(int)(m1*n1)][(int)(m2*n2)] += invcellvol;
321 set_pbc(&pbc, ePBC, box);
322 for (i = 0; i < 2; i++)
326 /* One atom, just copy the coordinates */
327 copy_rvec(x[ind[i][0]], xcom[i]);
331 /* Calculate the center of mass */
334 for (j = 0; j < gnx[i]; j++)
337 m = top.atoms.atom[k].m;
338 for (l = 0; l < DIM; l++)
340 xcom[i][l] += m*x[k][l];
344 svmul(1/mtot, xcom[i], xcom[i]);
347 pbc_dx(&pbc, xcom[1], xcom[0], direction);
348 for (i = 0; i < DIM; i++)
350 center[i] = xcom[0][i] + 0.5*direction[i];
352 unitv(direction, direction);
353 for (i = 0; i < nindex; i++)
356 pbc_dx(&pbc, x[j], center, dx);
357 axial = iprod(dx, direction);
358 r = sqrt(norm2(dx) - axial*axial);
359 if (axial >= -amax && axial < amax && r < rmax)
365 grid[(int)((axial + amax)*invspa)][(int)(r*invspz)] += 1;
371 while (read_next_x(oenv, status, &t, x, box));
374 /* normalize gridpoints */
378 for (i = 0; i < n1; i++)
380 for (j = 0; j < n2; j++)
383 if (grid[i][j] > maxgrid)
385 maxgrid = grid[i][j];
392 for (i = 0; i < n1; i++)
395 for (j = 0; j < nradial; j++)
399 case -3: vol = M_PI*(j+1)*(j+1)/(invspz*invspz*invspa); break;
400 case -2: vol = (j+1)/(invspz*invspa); break;
401 default: vol = j+1; break;
411 grid[i][k] /= nfr*(vol - vol_old);
414 grid[i][nradial-1-j] = grid[i][k];
417 if (grid[i][k] > maxgrid)
419 maxgrid = grid[i][k];
424 fprintf(stdout, "\n The maximum density is %f %s\n", maxgrid, unit);
434 /* normalize box-axes */
437 for (i = 0; i <= n1; i++)
439 tickx[i] = i*box1/n1;
441 for (i = 0; i <= n2; i++)
443 tickz[i] = i*box2/n2;
448 for (i = 0; i <= n1; i++)
450 tickx[i] = i/invspa - amax;
454 for (i = 0; i <= n2; i++)
456 tickz[i] = i/invspz - rmax;
461 for (i = 0; i <= n2; i++)
470 for (i = 0; i < n1; ++i)
472 fprintf(stdout, "Density sums:\n");
474 for (j = 0; j < n2; ++j)
476 rowsum += grid[i][j];
478 fprintf(stdout, "%g\t", rowsum);
480 fprintf(stdout, "\n");
483 sprintf(buf, "%s number density", grpname[anagrp]);
484 if (!bRadial && (bXmin || bXmax))
488 sprintf(buf+strlen(buf), ", %c > %g nm", eaver[0][0], xmin);
492 sprintf(buf+strlen(buf), ", %c < %g nm", eaver[0][0], xmax);
496 sprintf(buf+strlen(buf), ", %c: %g - %g nm", eaver[0][0], xmin, xmax);
499 if (ftp2bSet(efDAT, NFILE, fnm))
501 fp = ffopen(ftp2fn(efDAT, NFILE, fnm), "w");
502 /*optional text form output: first row is tickz; first col is tickx */
504 for (j = 0; j < n2; ++j)
506 fprintf(fp, "%g\t", tickz[j]);
510 for (i = 0; i < n1; ++i)
512 fprintf(fp, "%g\t", tickx[i]);
513 for (j = 0; j < n2; ++j)
515 fprintf(fp, "%g\t", grid[i][j]);
523 fp = ffopen(ftp2fn(efXPM, NFILE, fnm), "w");
524 write_xpm(fp, MAT_SPATIAL_X | MAT_SPATIAL_Y, buf, unit,
525 bRadial ? "axial (nm)" : label[c1], bRadial ? "r (nm)" : label[c2],
526 n1, n2, tickx, tickz, grid, dmin, maxgrid, rlo, rhi, &nlev);
530 do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);