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.
46 #include "gromacs/utility/cstringutil.h"
51 #include "gromacs/pbcutil/pbc.h"
52 #include "gromacs/topology/index.h"
53 #include "gromacs/fileio/tpxio.h"
54 #include "gromacs/fileio/trxio.h"
55 #include "gromacs/math/units.h"
59 #include "gromacs/commandline/pargs.h"
60 #include "gromacs/fileio/xvgr.h"
61 #include "gromacs/math/vec.h"
62 #include "gromacs/pbcutil/rmpbc.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/futil.h"
65 #include "gromacs/utility/smalloc.h"
72 /****************************************************************************/
73 /* This program calculates the partial density across the box. */
74 /* Peter Tieleman, Mei 1995 */
75 /****************************************************************************/
77 /* used for sorting the list */
78 int compare(void *a, void *b)
80 t_electron *tmp1, *tmp2;
81 tmp1 = (t_electron *)a; tmp2 = (t_electron *)b;
83 return strcmp(tmp1->atomname, tmp2->atomname);
86 int get_electrons(t_electron **eltab, const char *fn)
88 char buffer[256]; /* to read in a line */
89 char tempname[80]; /* buffer to hold name */
93 int nr; /* number of atomstypes to read */
96 if (!(in = gmx_ffopen(fn, "r")))
98 gmx_fatal(FARGS, "Couldn't open %s. Exiting.\n", fn);
101 if (NULL == fgets(buffer, 255, in))
103 gmx_fatal(FARGS, "Error reading from file %s", fn);
106 if (sscanf(buffer, "%d", &nr) != 1)
108 gmx_fatal(FARGS, "Invalid number of atomtypes in datafile\n");
113 for (i = 0; i < nr; i++)
115 if (fgets(buffer, 255, in) == NULL)
117 gmx_fatal(FARGS, "reading datafile. Check your datafile.\n");
119 if (sscanf(buffer, "%s = %d", tempname, &tempnr) != 2)
121 gmx_fatal(FARGS, "Invalid line in datafile at line %d\n", i+1);
123 (*eltab)[i].nr_el = tempnr;
124 (*eltab)[i].atomname = strdup(tempname);
129 fprintf(stderr, "Sorting list..\n");
130 qsort ((void*)*eltab, nr, sizeof(t_electron),
131 (int(*)(const void*, const void*))compare);
136 void center_coords(t_atoms *atoms, matrix box, rvec x0[], int axis)
140 rvec com, shift, box_center;
144 for (i = 0; (i < atoms->nr); i++)
146 mm = atoms->atom[i].m;
148 for (m = 0; (m < DIM); m++)
150 com[m] += mm*x0[i][m];
153 for (m = 0; (m < DIM); m++)
157 calc_box_center(ecenterDEF, box, box_center);
158 rvec_sub(box_center, com, shift);
159 shift[axis] -= box_center[axis];
161 for (i = 0; (i < atoms->nr); i++)
163 rvec_dec(x0[i], shift);
167 void calc_electron_density(const char *fn, atom_id **index, int gnx[],
168 double ***slDensity, int *nslices, t_topology *top,
170 int axis, int nr_grps, real *slWidth,
171 t_electron eltab[], int nr, gmx_bool bCenter,
172 const output_env_t oenv)
174 rvec *x0; /* coordinates without pbc */
175 matrix box; /* box (3x3) */
177 int natoms; /* nr. atoms in trj */
179 int i, n, /* loop indices */
180 nr_frames = 0, /* number of frames */
181 slice; /* current slice */
182 t_electron *found; /* found by bsearch */
183 t_electron sought; /* thingie thought by bsearch */
184 gmx_rmpbc_t gpbc = NULL;
189 if (axis < 0 || axis >= DIM)
191 gmx_fatal(FARGS, "Invalid axes. Terminating\n");
194 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
196 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
201 *nslices = (int)(box[axis][axis] * 10); /* default value */
203 fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
205 snew(*slDensity, nr_grps);
206 for (i = 0; i < nr_grps; i++)
208 snew((*slDensity)[i], *nslices);
211 gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
212 /*********** Start processing trajectory ***********/
215 gmx_rmpbc(gpbc, natoms, box, x0);
219 center_coords(&top->atoms, box, x0, axis);
222 *slWidth = box[axis][axis]/(*nslices);
223 invvol = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
225 for (n = 0; n < nr_grps; n++)
227 for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
229 z = x0[index[n][i]][axis];
232 z += box[axis][axis];
234 while (z > box[axis][axis])
236 z -= box[axis][axis];
239 /* determine which slice atom is in */
240 slice = (z / (*slWidth));
242 sought.atomname = strdup(*(top->atoms.atomname[index[n][i]]));
244 /* now find the number of electrons. This is not efficient. */
245 found = (t_electron *)
246 bsearch((const void *)&sought,
247 (const void *)eltab, nr, sizeof(t_electron),
248 (int(*)(const void*, const void*))compare);
252 fprintf(stderr, "Couldn't find %s. Add it to the .dat file\n",
253 *(top->atoms.atomname[index[n][i]]));
257 (*slDensity)[n][slice] += (found->nr_el -
258 top->atoms.atom[index[n][i]].q)*invvol;
260 free(sought.atomname);
265 while (read_next_x(oenv, status, &t, x0, box));
266 gmx_rmpbc_done(gpbc);
268 /*********** done with status file **********/
271 /* slDensity now contains the total number of electrons per slice, summed
272 over all frames. Now divide by nr_frames and volume of slice
275 fprintf(stderr, "\nRead %d frames from trajectory. Counting electrons\n",
278 for (n = 0; n < nr_grps; n++)
280 for (i = 0; i < *nslices; i++)
282 (*slDensity)[n][i] /= nr_frames;
286 sfree(x0); /* free memory used by coordinate array */
289 void calc_density(const char *fn, atom_id **index, int gnx[],
290 double ***slDensity, int *nslices, t_topology *top, int ePBC,
291 int axis, int nr_grps, real *slWidth, gmx_bool bCenter,
292 const output_env_t oenv)
294 rvec *x0; /* coordinates without pbc */
295 matrix box; /* box (3x3) */
297 int natoms; /* nr. atoms in trj */
299 int **slCount, /* nr. of atoms in one slice for a group */
300 i, j, n, /* loop indices */
303 nr_frames = 0, /* number of frames */
304 slice; /* current slice */
307 char *buf; /* for tmp. keeping atomname */
308 gmx_rmpbc_t gpbc = NULL;
310 if (axis < 0 || axis >= DIM)
312 gmx_fatal(FARGS, "Invalid axes. Terminating\n");
315 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
317 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
322 *nslices = (int)(box[axis][axis] * 10); /* default value */
323 fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
326 snew(*slDensity, nr_grps);
327 for (i = 0; i < nr_grps; i++)
329 snew((*slDensity)[i], *nslices);
332 gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
333 /*********** Start processing trajectory ***********/
336 gmx_rmpbc(gpbc, natoms, box, x0);
340 center_coords(&top->atoms, box, x0, axis);
343 *slWidth = box[axis][axis]/(*nslices);
344 invvol = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
347 for (n = 0; n < nr_grps; n++)
349 for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
351 z = x0[index[n][i]][axis];
354 z += box[axis][axis];
356 while (z > box[axis][axis])
358 z -= box[axis][axis];
361 /* determine which slice atom is in */
362 slice = (int)(z / (*slWidth));
363 (*slDensity)[n][slice] += top->atoms.atom[index[n][i]].m*invvol;
368 while (read_next_x(oenv, status, &t, x0, box));
369 gmx_rmpbc_done(gpbc);
371 /*********** done with status file **********/
374 /* slDensity now contains the total mass per slice, summed over all
375 frames. Now divide by nr_frames and volume of slice
378 fprintf(stderr, "\nRead %d frames from trajectory. Calculating density\n",
381 for (n = 0; n < nr_grps; n++)
383 for (i = 0; i < *nslices; i++)
385 (*slDensity)[n][i] /= nr_frames;
389 sfree(x0); /* free memory used by coordinate array */
392 void plot_density(double *slDensity[], const char *afile, int nslices,
393 int nr_grps, char *grpname[], real slWidth,
394 const char **dens_opt,
395 gmx_bool bSymmetrize, const output_env_t oenv)
398 const char *ylabel = NULL;
402 switch (dens_opt[0][0])
404 case 'm': ylabel = "Density (kg m\\S-3\\N)"; break;
405 case 'n': ylabel = "Number density (nm\\S-3\\N)"; break;
406 case 'c': ylabel = "Charge density (e nm\\S-3\\N)"; break;
407 case 'e': ylabel = "Electron density (e nm\\S-3\\N)"; break;
410 den = xvgropen(afile, "Partial densities", "Box (nm)", ylabel, oenv);
412 xvgr_legend(den, nr_grps, (const char**)grpname, oenv);
414 for (slice = 0; (slice < nslices); slice++)
416 fprintf(den, "%12g ", slice * slWidth);
417 for (n = 0; (n < nr_grps); n++)
421 ddd = (slDensity[n][slice]+slDensity[n][nslices-slice-1])*0.5;
425 ddd = slDensity[n][slice];
427 if (dens_opt[0][0] == 'm')
429 fprintf(den, " %12g", ddd*AMU/(NANO*NANO*NANO));
433 fprintf(den, " %12g", ddd);
442 int gmx_density(int argc, char *argv[])
444 const char *desc[] = {
445 "[THISMODULE] computes partial densities across the box, using an index file.[PAR]",
446 "For the total density of NPT simulations, use [gmx-energy] instead.",
448 "Densities are in kg/m^3, and number densities or electron densities can also be",
449 "calculated. For electron densities, a file describing the number of",
450 "electrons for each type of atom should be provided using [TT]-ei[tt].",
451 "It should look like:[BR]",
453 " [TT]atomname = nrelectrons[tt][BR]",
454 " [TT]atomname = nrelectrons[tt][BR]",
455 "The first line contains the number of lines to read from the file.",
456 "There should be one line for each unique atom name in your system.",
457 "The number of electrons for each atom is modified by its atomic",
462 static const char *dens_opt[] =
463 { NULL, "mass", "number", "charge", "electron", NULL };
464 static int axis = 2; /* normal to memb. default z */
465 static const char *axtitle = "Z";
466 static int nslices = 50; /* nr of slices defined */
467 static int ngrps = 1; /* nr. of groups */
468 static gmx_bool bSymmetrize = FALSE;
469 static gmx_bool bCenter = FALSE;
471 { "-d", FALSE, etSTR, {&axtitle},
472 "Take the normal on the membrane in direction X, Y or Z." },
473 { "-sl", FALSE, etINT, {&nslices},
474 "Divide the box in this number of slices." },
475 { "-dens", FALSE, etENUM, {dens_opt},
477 { "-ng", FALSE, etINT, {&ngrps},
478 "Number of groups of which to compute densities." },
479 { "-symm", FALSE, etBOOL, {&bSymmetrize},
480 "Symmetrize the density along the axis, with respect to the center. Useful for bilayers." },
481 { "-center", FALSE, etBOOL, {&bCenter},
482 "Shift the center of mass along the axis to zero. This means if your axis is Z and your box is bX, bY, bZ, the center of mass will be at bX/2, bY/2, 0."}
485 const char *bugs[] = {
486 "When calculating electron densities, atomnames are used instead of types. This is bad.",
489 double **density; /* density per slice */
490 real slWidth; /* width of one slice */
491 char **grpname; /* groupnames */
492 int nr_electrons; /* nr. electrons */
493 int *ngx; /* sizes of groups */
494 t_electron *el_tab; /* tabel with nr. of electrons*/
495 t_topology *top; /* topology */
497 atom_id **index; /* indices for all groups */
500 t_filenm fnm[] = { /* files for g_density */
501 { efTRX, "-f", NULL, ffREAD },
502 { efNDX, NULL, NULL, ffOPTRD },
503 { efTPX, NULL, NULL, ffREAD },
504 { efDAT, "-ei", "electrons", ffOPTRD }, /* file with nr. of electrons */
505 { efXVG, "-o", "density", ffWRITE },
508 #define NFILE asize(fnm)
510 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
511 NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs,
517 if (bSymmetrize && !bCenter)
519 fprintf(stderr, "Can not symmetrize without centering. Turning on -center\n");
523 axis = toupper(axtitle[0]) - 'X';
525 top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC); /* read topology file */
526 if (dens_opt[0][0] == 'n')
528 for (i = 0; (i < top->atoms.nr); i++)
530 top->atoms.atom[i].m = 1;
533 else if (dens_opt[0][0] == 'c')
535 for (i = 0; (i < top->atoms.nr); i++)
537 top->atoms.atom[i].m = top->atoms.atom[i].q;
541 snew(grpname, ngrps);
545 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, ngx, index, grpname);
547 if (dens_opt[0][0] == 'e')
549 nr_electrons = get_electrons(&el_tab, ftp2fn(efDAT, NFILE, fnm));
550 fprintf(stderr, "Read %d atomtypes from datafile\n", nr_electrons);
552 calc_electron_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density,
553 &nslices, top, ePBC, axis, ngrps, &slWidth, el_tab,
554 nr_electrons, bCenter, oenv);
558 calc_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density, &nslices, top,
559 ePBC, axis, ngrps, &slWidth, bCenter, oenv);
562 plot_density(density, opt2fn("-o", NFILE, fnm),
563 nslices, ngrps, grpname, slWidth, dens_opt,
566 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy"); /* view xvgr file */