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
42 #include "gromacs/commandline/pargs.h"
49 #include "gromacs/fileio/futil.h"
56 #include "gromacs/fileio/tpxio.h"
57 #include "gromacs/fileio/trxio.h"
59 #include "gromacs/fileio/matio.h"
64 int gmx_densmap(int argc, char *argv[])
66 const char *desc[] = {
67 "[THISMODULE] computes 2D number-density maps.",
68 "It can make planar and axial-radial density maps.",
69 "The output [TT].xpm[tt] file can be visualized with for instance xv",
70 "and can be converted to postscript with [TT]xpm2ps[tt].",
71 "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].",
73 "The default analysis is a 2-D number-density map for a selected",
74 "group of atoms in the x-y plane.",
75 "The averaging direction can be changed with the option [TT]-aver[tt].",
76 "When [TT]-xmin[tt] and/or [TT]-xmax[tt] are set only atoms that are",
77 "within the limit(s) in the averaging direction are taken into account.",
78 "The grid spacing is set with the option [TT]-bin[tt].",
79 "When [TT]-n1[tt] or [TT]-n2[tt] is non-zero, the grid",
80 "size is set by this option.",
81 "Box size fluctuations are properly taken into account.",
83 "When options [TT]-amax[tt] and [TT]-rmax[tt] are set, an axial-radial",
84 "number-density map is made. Three groups should be supplied, the centers",
85 "of mass of the first two groups define the axis, the third defines the",
86 "analysis group. The axial direction goes from -amax to +amax, where",
87 "the center is defined as the midpoint between the centers of mass and",
88 "the positive direction goes from the first to the second center of mass.",
89 "The radial direction goes from 0 to rmax or from -rmax to +rmax",
90 "when the [TT]-mirror[tt] option has been set.",
92 "The normalization of the output is set with the [TT]-unit[tt] option.",
93 "The default produces a true number density. Unit [TT]nm-2[tt] leaves out",
94 "the normalization for the averaging or the angular direction.",
95 "Option [TT]count[tt] produces the count for each grid cell.",
96 "When you do not want the scale in the output to go",
97 "from zero to the maximum density, you can set the maximum",
98 "with the option [TT]-dmax[tt]."
100 static int n1 = 0, n2 = 0;
101 static real xmin = -1, xmax = -1, bin = 0.02, dmin = 0, dmax = 0, amax = 0, rmax = 0;
102 static gmx_bool bMirror = FALSE, bSums = FALSE;
103 static const char *eaver[] = { NULL, "z", "y", "x", NULL };
104 static const char *eunit[] = { NULL, "nm-3", "nm-2", "count", NULL };
107 { "-bin", FALSE, etREAL, {&bin},
109 { "-aver", FALSE, etENUM, {eaver},
110 "The direction to average over" },
111 { "-xmin", FALSE, etREAL, {&xmin},
112 "Minimum coordinate for averaging" },
113 { "-xmax", FALSE, etREAL, {&xmax},
114 "Maximum coordinate for averaging" },
115 { "-n1", FALSE, etINT, {&n1},
116 "Number of grid cells in the first direction" },
117 { "-n2", FALSE, etINT, {&n2},
118 "Number of grid cells in the second direction" },
119 { "-amax", FALSE, etREAL, {&amax},
120 "Maximum axial distance from the center"},
121 { "-rmax", FALSE, etREAL, {&rmax},
122 "Maximum radial distance" },
123 { "-mirror", FALSE, etBOOL, {&bMirror},
124 "Add the mirror image below the axial axis" },
125 { "-sums", FALSE, etBOOL, {&bSums},
126 "Print density sums (1D map) to stdout" },
127 { "-unit", FALSE, etENUM, {eunit},
128 "Unit for the output" },
129 { "-dmin", FALSE, etREAL, {&dmin},
130 "Minimum density in output"},
131 { "-dmax", FALSE, etREAL, {&dmax},
132 "Maximum density in output (0 means calculate it)"},
134 gmx_bool bXmin, bXmax, bRadial;
139 rvec *x, xcom[2], direction, center, dx;
143 int cav = 0, c1 = 0, c2 = 0, natoms;
144 char **grpname, title[256], buf[STRLEN];
146 int i, j, k, l, ngrps, anagrp, *gnx = NULL, nindex, nradial = 0, nfr, nmpower;
147 atom_id **ind = NULL, *index;
148 real **grid, maxgrid, m1, m2, box1, box2, *tickx, *tickz, invcellvol;
149 real invspa = 0, invspz = 0, axial, r, vol_old, vol, rowsum;
151 t_rgb rlo = {1, 1, 1}, rhi = {0, 0, 0};
153 const char *label[] = { "x (nm)", "y (nm)", "z (nm)" };
155 { efTRX, "-f", NULL, ffREAD },
156 { efTPS, NULL, NULL, ffOPTRD },
157 { efNDX, NULL, NULL, ffOPTRD },
158 { efDAT, "-od", "densmap", ffOPTWR },
159 { efXPM, "-o", "densmap", ffWRITE }
161 #define NFILE asize(fnm)
166 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
167 NFILE, fnm, npargs, pa, asize(desc), desc, 0, NULL, &oenv))
172 bXmin = opt2parg_bSet("-xmin", npargs, pa);
173 bXmax = opt2parg_bSet("-xmax", npargs, pa);
174 bRadial = (amax > 0 || rmax > 0);
177 if (amax <= 0 || rmax <= 0)
179 gmx_fatal(FARGS, "Both amax and rmax should be larger than zero");
183 if (strcmp(eunit[0], "nm-3") == 0)
188 else if (strcmp(eunit[0], "nm-2") == 0)
199 if (ftp2bSet(efTPS, NFILE, fnm) || !ftp2bSet(efNDX, NFILE, fnm))
201 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &x, NULL, box,
207 fprintf(stderr, "\nSelect an analysis group\n");
213 "\nSelect two groups to define the axis and an analysis group\n");
216 snew(grpname, ngrps);
218 get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, gnx, ind, grpname);
220 nindex = gnx[anagrp];
224 if ((gnx[0] > 1 || gnx[1] > 1) && !ftp2bSet(efTPS, NFILE, fnm))
226 gmx_fatal(FARGS, "No run input file was supplied (option -s), this is required for the center of mass calculation");
232 case 'x': cav = XX; c1 = YY; c2 = ZZ; break;
233 case 'y': cav = YY; c1 = XX; c2 = ZZ; break;
234 case 'z': cav = ZZ; c1 = XX; c2 = YY; break;
237 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
243 n1 = (int)(box[c1][c1]/bin + 0.5);
247 n2 = (int)(box[c2][c2]/bin + 0.5);
252 n1 = (int)(2*amax/bin + 0.5);
253 nradial = (int)(rmax/bin + 0.5);
254 invspa = n1/(2*amax);
255 invspz = nradial/rmax;
267 for (i = 0; i < n1; i++)
284 invcellvol /= det(box);
286 else if (nmpower == -2)
288 invcellvol /= box[c1][c1]*box[c2][c2];
290 for (i = 0; i < nindex; i++)
293 if ((!bXmin || x[j][cav] >= xmin) &&
294 (!bXmax || x[j][cav] <= xmax))
296 m1 = x[j][c1]/box[c1][c1];
305 m2 = x[j][c2]/box[c2][c2];
314 grid[(int)(m1*n1)][(int)(m2*n2)] += invcellvol;
320 set_pbc(&pbc, ePBC, box);
321 for (i = 0; i < 2; i++)
325 /* One atom, just copy the coordinates */
326 copy_rvec(x[ind[i][0]], xcom[i]);
330 /* Calculate the center of mass */
333 for (j = 0; j < gnx[i]; j++)
336 m = top.atoms.atom[k].m;
337 for (l = 0; l < DIM; l++)
339 xcom[i][l] += m*x[k][l];
343 svmul(1/mtot, xcom[i], xcom[i]);
346 pbc_dx(&pbc, xcom[1], xcom[0], direction);
347 for (i = 0; i < DIM; i++)
349 center[i] = xcom[0][i] + 0.5*direction[i];
351 unitv(direction, direction);
352 for (i = 0; i < nindex; i++)
355 pbc_dx(&pbc, x[j], center, dx);
356 axial = iprod(dx, direction);
357 r = sqrt(norm2(dx) - axial*axial);
358 if (axial >= -amax && axial < amax && r < rmax)
364 grid[(int)((axial + amax)*invspa)][(int)(r*invspz)] += 1;
370 while (read_next_x(oenv, status, &t, x, box));
373 /* normalize gridpoints */
377 for (i = 0; i < n1; i++)
379 for (j = 0; j < n2; j++)
382 if (grid[i][j] > maxgrid)
384 maxgrid = grid[i][j];
391 for (i = 0; i < n1; i++)
394 for (j = 0; j < nradial; j++)
398 case -3: vol = M_PI*(j+1)*(j+1)/(invspz*invspz*invspa); break;
399 case -2: vol = (j+1)/(invspz*invspa); break;
400 default: vol = j+1; break;
410 grid[i][k] /= nfr*(vol - vol_old);
413 grid[i][nradial-1-j] = grid[i][k];
416 if (grid[i][k] > maxgrid)
418 maxgrid = grid[i][k];
423 fprintf(stdout, "\n The maximum density is %f %s\n", maxgrid, unit);
433 /* normalize box-axes */
436 for (i = 0; i <= n1; i++)
438 tickx[i] = i*box1/n1;
440 for (i = 0; i <= n2; i++)
442 tickz[i] = i*box2/n2;
447 for (i = 0; i <= n1; i++)
449 tickx[i] = i/invspa - amax;
453 for (i = 0; i <= n2; i++)
455 tickz[i] = i/invspz - rmax;
460 for (i = 0; i <= n2; i++)
469 for (i = 0; i < n1; ++i)
471 fprintf(stdout, "Density sums:\n");
473 for (j = 0; j < n2; ++j)
475 rowsum += grid[i][j];
477 fprintf(stdout, "%g\t", rowsum);
479 fprintf(stdout, "\n");
482 sprintf(buf, "%s number density", grpname[anagrp]);
483 if (!bRadial && (bXmin || bXmax))
487 sprintf(buf+strlen(buf), ", %c > %g nm", eaver[0][0], xmin);
491 sprintf(buf+strlen(buf), ", %c < %g nm", eaver[0][0], xmax);
495 sprintf(buf+strlen(buf), ", %c: %g - %g nm", eaver[0][0], xmin, xmax);
498 if (ftp2bSet(efDAT, NFILE, fnm))
500 fp = ffopen(ftp2fn(efDAT, NFILE, fnm), "w");
501 /*optional text form output: first row is tickz; first col is tickx */
503 for (j = 0; j < n2; ++j)
505 fprintf(fp, "%g\t", tickz[j]);
509 for (i = 0; i < n1; ++i)
511 fprintf(fp, "%g\t", tickx[i]);
512 for (j = 0; j < n2; ++j)
514 fprintf(fp, "%g\t", grid[i][j]);
522 fp = ffopen(ftp2fn(efXPM, NFILE, fnm), "w");
523 write_xpm(fp, MAT_SPATIAL_X | MAT_SPATIAL_Y, buf, unit,
524 bRadial ? "axial (nm)" : label[c1], bRadial ? "r (nm)" : label[c2],
525 n1, n2, tickx, tickz, grid, dmin, maxgrid, rlo, rhi, &nlev);
529 do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);