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, 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.
53 #include "gromacs/fileio/futil.h"
54 #include "gromacs/commandline/pargs.h"
56 #include "gromacs/fileio/tpxio.h"
57 #include "gromacs/fileio/trxio.h"
67 /****************************************************************************/
68 /* This program calculates the partial density across the box. */
69 /* Peter Tieleman, Mei 1995 */
70 /****************************************************************************/
72 /* used for sorting the list */
73 int compare(void *a, void *b)
75 t_electron *tmp1, *tmp2;
76 tmp1 = (t_electron *)a; tmp2 = (t_electron *)b;
78 return strcmp(tmp1->atomname, tmp2->atomname);
81 int get_electrons(t_electron **eltab, const char *fn)
83 char buffer[256]; /* to read in a line */
84 char tempname[80]; /* buffer to hold name */
88 int nr; /* number of atomstypes to read */
91 if (!(in = ffopen(fn, "r")))
93 gmx_fatal(FARGS, "Couldn't open %s. Exiting.\n", fn);
96 if (NULL == fgets(buffer, 255, in))
98 gmx_fatal(FARGS, "Error reading from file %s", fn);
101 if (sscanf(buffer, "%d", &nr) != 1)
103 gmx_fatal(FARGS, "Invalid number of atomtypes in datafile\n");
108 for (i = 0; i < nr; i++)
110 if (fgets(buffer, 255, in) == NULL)
112 gmx_fatal(FARGS, "reading datafile. Check your datafile.\n");
114 if (sscanf(buffer, "%s = %d", tempname, &tempnr) != 2)
116 gmx_fatal(FARGS, "Invalid line in datafile at line %d\n", i+1);
118 (*eltab)[i].nr_el = tempnr;
119 (*eltab)[i].atomname = strdup(tempname);
124 fprintf(stderr, "Sorting list..\n");
125 qsort ((void*)*eltab, nr, sizeof(t_electron),
126 (int(*)(const void*, const void*))compare);
131 void center_coords(t_atoms *atoms, matrix box, rvec x0[], int axis)
135 rvec com, shift, box_center;
139 for (i = 0; (i < atoms->nr); i++)
141 mm = atoms->atom[i].m;
143 for (m = 0; (m < DIM); m++)
145 com[m] += mm*x0[i][m];
148 for (m = 0; (m < DIM); m++)
152 calc_box_center(ecenterDEF, box, box_center);
153 rvec_sub(box_center, com, shift);
154 shift[axis] -= box_center[axis];
156 for (i = 0; (i < atoms->nr); i++)
158 rvec_dec(x0[i], shift);
162 void calc_electron_density(const char *fn, atom_id **index, int gnx[],
163 double ***slDensity, int *nslices, t_topology *top,
165 int axis, int nr_grps, real *slWidth,
166 t_electron eltab[], int nr, gmx_bool bCenter,
167 const output_env_t oenv)
169 rvec *x0; /* coordinates without pbc */
170 matrix box; /* box (3x3) */
172 int natoms; /* nr. atoms in trj */
174 int i, n, /* loop indices */
175 nr_frames = 0, /* number of frames */
176 slice; /* current slice */
177 t_electron *found; /* found by bsearch */
178 t_electron sought; /* thingie thought by bsearch */
179 gmx_rmpbc_t gpbc = NULL;
184 if (axis < 0 || axis >= DIM)
186 gmx_fatal(FARGS, "Invalid axes. Terminating\n");
189 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
191 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
196 *nslices = (int)(box[axis][axis] * 10); /* default value */
198 fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
200 snew(*slDensity, nr_grps);
201 for (i = 0; i < nr_grps; i++)
203 snew((*slDensity)[i], *nslices);
206 gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
207 /*********** Start processing trajectory ***********/
210 gmx_rmpbc(gpbc, natoms, box, x0);
214 center_coords(&top->atoms, box, x0, axis);
217 *slWidth = box[axis][axis]/(*nslices);
218 invvol = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
220 for (n = 0; n < nr_grps; n++)
222 for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
224 z = x0[index[n][i]][axis];
227 z += box[axis][axis];
229 while (z > box[axis][axis])
231 z -= box[axis][axis];
234 /* determine which slice atom is in */
235 slice = (z / (*slWidth));
237 sought.atomname = strdup(*(top->atoms.atomname[index[n][i]]));
239 /* now find the number of electrons. This is not efficient. */
240 found = (t_electron *)
241 bsearch((const void *)&sought,
242 (const void *)eltab, nr, sizeof(t_electron),
243 (int(*)(const void*, const void*))compare);
247 fprintf(stderr, "Couldn't find %s. Add it to the .dat file\n",
248 *(top->atoms.atomname[index[n][i]]));
252 (*slDensity)[n][slice] += (found->nr_el -
253 top->atoms.atom[index[n][i]].q)*invvol;
255 free(sought.atomname);
260 while (read_next_x(oenv, status, &t, x0, box));
261 gmx_rmpbc_done(gpbc);
263 /*********** done with status file **********/
266 /* slDensity now contains the total number of electrons per slice, summed
267 over all frames. Now divide by nr_frames and volume of slice
270 fprintf(stderr, "\nRead %d frames from trajectory. Counting electrons\n",
273 for (n = 0; n < nr_grps; n++)
275 for (i = 0; i < *nslices; i++)
277 (*slDensity)[n][i] /= nr_frames;
281 sfree(x0); /* free memory used by coordinate array */
284 void calc_density(const char *fn, atom_id **index, int gnx[],
285 double ***slDensity, int *nslices, t_topology *top, int ePBC,
286 int axis, int nr_grps, real *slWidth, gmx_bool bCenter,
287 const output_env_t oenv)
289 rvec *x0; /* coordinates without pbc */
290 matrix box; /* box (3x3) */
292 int natoms; /* nr. atoms in trj */
294 int **slCount, /* nr. of atoms in one slice for a group */
295 i, j, n, /* loop indices */
298 nr_frames = 0, /* number of frames */
299 slice; /* current slice */
302 char *buf; /* for tmp. keeping atomname */
303 gmx_rmpbc_t gpbc = NULL;
305 if (axis < 0 || axis >= DIM)
307 gmx_fatal(FARGS, "Invalid axes. Terminating\n");
310 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
312 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
317 *nslices = (int)(box[axis][axis] * 10); /* default value */
318 fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
321 snew(*slDensity, nr_grps);
322 for (i = 0; i < nr_grps; i++)
324 snew((*slDensity)[i], *nslices);
327 gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
328 /*********** Start processing trajectory ***********/
331 gmx_rmpbc(gpbc, natoms, box, x0);
335 center_coords(&top->atoms, box, x0, axis);
338 *slWidth = box[axis][axis]/(*nslices);
339 invvol = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
342 for (n = 0; n < nr_grps; n++)
344 for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
346 z = x0[index[n][i]][axis];
349 z += box[axis][axis];
351 while (z > box[axis][axis])
353 z -= box[axis][axis];
356 /* determine which slice atom is in */
357 slice = (int)(z / (*slWidth));
358 (*slDensity)[n][slice] += top->atoms.atom[index[n][i]].m*invvol;
363 while (read_next_x(oenv, status, &t, x0, box));
364 gmx_rmpbc_done(gpbc);
366 /*********** done with status file **********/
369 /* slDensity now contains the total mass per slice, summed over all
370 frames. Now divide by nr_frames and volume of slice
373 fprintf(stderr, "\nRead %d frames from trajectory. Calculating density\n",
376 for (n = 0; n < nr_grps; n++)
378 for (i = 0; i < *nslices; i++)
380 (*slDensity)[n][i] /= nr_frames;
384 sfree(x0); /* free memory used by coordinate array */
387 void plot_density(double *slDensity[], const char *afile, int nslices,
388 int nr_grps, char *grpname[], real slWidth,
389 const char **dens_opt,
390 gmx_bool bSymmetrize, const output_env_t oenv)
393 const char *ylabel = NULL;
397 switch (dens_opt[0][0])
399 case 'm': ylabel = "Density (kg m\\S-3\\N)"; break;
400 case 'n': ylabel = "Number density (nm\\S-3\\N)"; break;
401 case 'c': ylabel = "Charge density (e nm\\S-3\\N)"; break;
402 case 'e': ylabel = "Electron density (e nm\\S-3\\N)"; break;
405 den = xvgropen(afile, "Partial densities", "Box (nm)", ylabel, oenv);
407 xvgr_legend(den, nr_grps, (const char**)grpname, oenv);
409 for (slice = 0; (slice < nslices); slice++)
411 fprintf(den, "%12g ", slice * slWidth);
412 for (n = 0; (n < nr_grps); n++)
416 ddd = (slDensity[n][slice]+slDensity[n][nslices-slice-1])*0.5;
420 ddd = slDensity[n][slice];
422 if (dens_opt[0][0] == 'm')
424 fprintf(den, " %12g", ddd*AMU/(NANO*NANO*NANO));
428 fprintf(den, " %12g", ddd);
437 int gmx_density(int argc, char *argv[])
439 const char *desc[] = {
440 "[THISMODULE] computes partial densities across the box, using an index file.[PAR]",
441 "For the total density of NPT simulations, use [gmx-energy] instead.",
443 "Densities are in kg/m^3, and number densities or electron densities can also be",
444 "calculated. For electron densities, a file describing the number of",
445 "electrons for each type of atom should be provided using [TT]-ei[tt].",
446 "It should look like:[BR]",
448 " [TT]atomname = nrelectrons[tt][BR]",
449 " [TT]atomname = nrelectrons[tt][BR]",
450 "The first line contains the number of lines to read from the file.",
451 "There should be one line for each unique atom name in your system.",
452 "The number of electrons for each atom is modified by its atomic",
457 static const char *dens_opt[] =
458 { NULL, "mass", "number", "charge", "electron", NULL };
459 static int axis = 2; /* normal to memb. default z */
460 static const char *axtitle = "Z";
461 static int nslices = 50; /* nr of slices defined */
462 static int ngrps = 1; /* nr. of groups */
463 static gmx_bool bSymmetrize = FALSE;
464 static gmx_bool bCenter = FALSE;
466 { "-d", FALSE, etSTR, {&axtitle},
467 "Take the normal on the membrane in direction X, Y or Z." },
468 { "-sl", FALSE, etINT, {&nslices},
469 "Divide the box in this number of slices." },
470 { "-dens", FALSE, etENUM, {dens_opt},
472 { "-ng", FALSE, etINT, {&ngrps},
473 "Number of groups of which to compute densities." },
474 { "-symm", FALSE, etBOOL, {&bSymmetrize},
475 "Symmetrize the density along the axis, with respect to the center. Useful for bilayers." },
476 { "-center", FALSE, etBOOL, {&bCenter},
477 "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."}
480 const char *bugs[] = {
481 "When calculating electron densities, atomnames are used instead of types. This is bad.",
484 double **density; /* density per slice */
485 real slWidth; /* width of one slice */
486 char **grpname; /* groupnames */
487 int nr_electrons; /* nr. electrons */
488 int *ngx; /* sizes of groups */
489 t_electron *el_tab; /* tabel with nr. of electrons*/
490 t_topology *top; /* topology */
492 atom_id **index; /* indices for all groups */
495 t_filenm fnm[] = { /* files for g_density */
496 { efTRX, "-f", NULL, ffREAD },
497 { efNDX, NULL, NULL, ffOPTRD },
498 { efTPX, NULL, NULL, ffREAD },
499 { efDAT, "-ei", "electrons", ffOPTRD }, /* file with nr. of electrons */
500 { efXVG, "-o", "density", ffWRITE },
503 #define NFILE asize(fnm)
505 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
506 NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs,
512 if (bSymmetrize && !bCenter)
514 fprintf(stderr, "Can not symmetrize without centering. Turning on -center\n");
518 axis = toupper(axtitle[0]) - 'X';
520 top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC); /* read topology file */
521 if (dens_opt[0][0] == 'n')
523 for (i = 0; (i < top->atoms.nr); i++)
525 top->atoms.atom[i].m = 1;
528 else if (dens_opt[0][0] == 'c')
530 for (i = 0; (i < top->atoms.nr); i++)
532 top->atoms.atom[i].m = top->atoms.atom[i].q;
536 snew(grpname, ngrps);
540 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, ngx, index, grpname);
542 if (dens_opt[0][0] == 'e')
544 nr_electrons = get_electrons(&el_tab, ftp2fn(efDAT, NFILE, fnm));
545 fprintf(stderr, "Read %d atomtypes from datafile\n", nr_electrons);
547 calc_electron_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density,
548 &nslices, top, ePBC, axis, ngrps, &slWidth, el_tab,
549 nr_electrons, bCenter, oenv);
553 calc_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density, &nslices, top,
554 ePBC, axis, ngrps, &slWidth, bCenter, oenv);
557 plot_density(density, opt2fn("-o", NFILE, fnm),
558 nslices, ngrps, grpname, slWidth, dens_opt,
561 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy"); /* view xvgr file */