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.
45 #include "gromacs/utility/cstringutil.h"
47 #include "gromacs/utility/smalloc.h"
53 #include "gromacs/fileio/futil.h"
54 #include "gromacs/commandline/pargs.h"
56 #include "gromacs/fileio/tpxio.h"
57 #include "gromacs/fileio/trxio.h"
62 #include "gromacs/legacyheaders/gmx_fatal.h"
69 /****************************************************************************/
70 /* This program calculates the partial density across the box. */
71 /* Peter Tieleman, Mei 1995 */
72 /****************************************************************************/
74 /* used for sorting the list */
75 int compare(void *a, void *b)
77 t_electron *tmp1, *tmp2;
78 tmp1 = (t_electron *)a; tmp2 = (t_electron *)b;
80 return strcmp(tmp1->atomname, tmp2->atomname);
83 int get_electrons(t_electron **eltab, const char *fn)
85 char buffer[256]; /* to read in a line */
86 char tempname[80]; /* buffer to hold name */
90 int nr; /* number of atomstypes to read */
93 if (!(in = gmx_ffopen(fn, "r")))
95 gmx_fatal(FARGS, "Couldn't open %s. Exiting.\n", fn);
98 if (NULL == fgets(buffer, 255, in))
100 gmx_fatal(FARGS, "Error reading from file %s", fn);
103 if (sscanf(buffer, "%d", &nr) != 1)
105 gmx_fatal(FARGS, "Invalid number of atomtypes in datafile\n");
110 for (i = 0; i < nr; i++)
112 if (fgets(buffer, 255, in) == NULL)
114 gmx_fatal(FARGS, "reading datafile. Check your datafile.\n");
116 if (sscanf(buffer, "%s = %d", tempname, &tempnr) != 2)
118 gmx_fatal(FARGS, "Invalid line in datafile at line %d\n", i+1);
120 (*eltab)[i].nr_el = tempnr;
121 (*eltab)[i].atomname = strdup(tempname);
126 fprintf(stderr, "Sorting list..\n");
127 qsort ((void*)*eltab, nr, sizeof(t_electron),
128 (int(*)(const void*, const void*))compare);
133 void center_coords(t_atoms *atoms, matrix box, rvec x0[], int axis)
137 rvec com, shift, box_center;
141 for (i = 0; (i < atoms->nr); i++)
143 mm = atoms->atom[i].m;
145 for (m = 0; (m < DIM); m++)
147 com[m] += mm*x0[i][m];
150 for (m = 0; (m < DIM); m++)
154 calc_box_center(ecenterDEF, box, box_center);
155 rvec_sub(box_center, com, shift);
156 shift[axis] -= box_center[axis];
158 for (i = 0; (i < atoms->nr); i++)
160 rvec_dec(x0[i], shift);
164 void calc_electron_density(const char *fn, atom_id **index, int gnx[],
165 double ***slDensity, int *nslices, t_topology *top,
167 int axis, int nr_grps, real *slWidth,
168 t_electron eltab[], int nr, gmx_bool bCenter,
169 const output_env_t oenv)
171 rvec *x0; /* coordinates without pbc */
172 matrix box; /* box (3x3) */
174 int natoms; /* nr. atoms in trj */
176 int i, n, /* loop indices */
177 nr_frames = 0, /* number of frames */
178 slice; /* current slice */
179 t_electron *found; /* found by bsearch */
180 t_electron sought; /* thingie thought by bsearch */
181 gmx_rmpbc_t gpbc = NULL;
186 if (axis < 0 || axis >= DIM)
188 gmx_fatal(FARGS, "Invalid axes. Terminating\n");
191 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
193 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
198 *nslices = (int)(box[axis][axis] * 10); /* default value */
200 fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
202 snew(*slDensity, nr_grps);
203 for (i = 0; i < nr_grps; i++)
205 snew((*slDensity)[i], *nslices);
208 gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
209 /*********** Start processing trajectory ***********/
212 gmx_rmpbc(gpbc, natoms, box, x0);
216 center_coords(&top->atoms, box, x0, axis);
219 *slWidth = box[axis][axis]/(*nslices);
220 invvol = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
222 for (n = 0; n < nr_grps; n++)
224 for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
226 z = x0[index[n][i]][axis];
229 z += box[axis][axis];
231 while (z > box[axis][axis])
233 z -= box[axis][axis];
236 /* determine which slice atom is in */
237 slice = (z / (*slWidth));
239 sought.atomname = strdup(*(top->atoms.atomname[index[n][i]]));
241 /* now find the number of electrons. This is not efficient. */
242 found = (t_electron *)
243 bsearch((const void *)&sought,
244 (const void *)eltab, nr, sizeof(t_electron),
245 (int(*)(const void*, const void*))compare);
249 fprintf(stderr, "Couldn't find %s. Add it to the .dat file\n",
250 *(top->atoms.atomname[index[n][i]]));
254 (*slDensity)[n][slice] += (found->nr_el -
255 top->atoms.atom[index[n][i]].q)*invvol;
257 free(sought.atomname);
262 while (read_next_x(oenv, status, &t, x0, box));
263 gmx_rmpbc_done(gpbc);
265 /*********** done with status file **********/
268 /* slDensity now contains the total number of electrons per slice, summed
269 over all frames. Now divide by nr_frames and volume of slice
272 fprintf(stderr, "\nRead %d frames from trajectory. Counting electrons\n",
275 for (n = 0; n < nr_grps; n++)
277 for (i = 0; i < *nslices; i++)
279 (*slDensity)[n][i] /= nr_frames;
283 sfree(x0); /* free memory used by coordinate array */
286 void calc_density(const char *fn, atom_id **index, int gnx[],
287 double ***slDensity, int *nslices, t_topology *top, int ePBC,
288 int axis, int nr_grps, real *slWidth, gmx_bool bCenter,
289 const output_env_t oenv)
291 rvec *x0; /* coordinates without pbc */
292 matrix box; /* box (3x3) */
294 int natoms; /* nr. atoms in trj */
296 int **slCount, /* nr. of atoms in one slice for a group */
297 i, j, n, /* loop indices */
300 nr_frames = 0, /* number of frames */
301 slice; /* current slice */
304 char *buf; /* for tmp. keeping atomname */
305 gmx_rmpbc_t gpbc = NULL;
307 if (axis < 0 || axis >= DIM)
309 gmx_fatal(FARGS, "Invalid axes. Terminating\n");
312 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
314 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
319 *nslices = (int)(box[axis][axis] * 10); /* default value */
320 fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
323 snew(*slDensity, nr_grps);
324 for (i = 0; i < nr_grps; i++)
326 snew((*slDensity)[i], *nslices);
329 gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
330 /*********** Start processing trajectory ***********/
333 gmx_rmpbc(gpbc, natoms, box, x0);
337 center_coords(&top->atoms, box, x0, axis);
340 *slWidth = box[axis][axis]/(*nslices);
341 invvol = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
344 for (n = 0; n < nr_grps; n++)
346 for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
348 z = x0[index[n][i]][axis];
351 z += box[axis][axis];
353 while (z > box[axis][axis])
355 z -= box[axis][axis];
358 /* determine which slice atom is in */
359 slice = (int)(z / (*slWidth));
360 (*slDensity)[n][slice] += top->atoms.atom[index[n][i]].m*invvol;
365 while (read_next_x(oenv, status, &t, x0, box));
366 gmx_rmpbc_done(gpbc);
368 /*********** done with status file **********/
371 /* slDensity now contains the total mass per slice, summed over all
372 frames. Now divide by nr_frames and volume of slice
375 fprintf(stderr, "\nRead %d frames from trajectory. Calculating density\n",
378 for (n = 0; n < nr_grps; n++)
380 for (i = 0; i < *nslices; i++)
382 (*slDensity)[n][i] /= nr_frames;
386 sfree(x0); /* free memory used by coordinate array */
389 void plot_density(double *slDensity[], const char *afile, int nslices,
390 int nr_grps, char *grpname[], real slWidth,
391 const char **dens_opt,
392 gmx_bool bSymmetrize, const output_env_t oenv)
395 const char *ylabel = NULL;
399 switch (dens_opt[0][0])
401 case 'm': ylabel = "Density (kg m\\S-3\\N)"; break;
402 case 'n': ylabel = "Number density (nm\\S-3\\N)"; break;
403 case 'c': ylabel = "Charge density (e nm\\S-3\\N)"; break;
404 case 'e': ylabel = "Electron density (e nm\\S-3\\N)"; break;
407 den = xvgropen(afile, "Partial densities", "Box (nm)", ylabel, oenv);
409 xvgr_legend(den, nr_grps, (const char**)grpname, oenv);
411 for (slice = 0; (slice < nslices); slice++)
413 fprintf(den, "%12g ", slice * slWidth);
414 for (n = 0; (n < nr_grps); n++)
418 ddd = (slDensity[n][slice]+slDensity[n][nslices-slice-1])*0.5;
422 ddd = slDensity[n][slice];
424 if (dens_opt[0][0] == 'm')
426 fprintf(den, " %12g", ddd*AMU/(NANO*NANO*NANO));
430 fprintf(den, " %12g", ddd);
439 int gmx_density(int argc, char *argv[])
441 const char *desc[] = {
442 "[THISMODULE] computes partial densities across the box, using an index file.[PAR]",
443 "For the total density of NPT simulations, use [gmx-energy] instead.",
445 "Densities are in kg/m^3, and number densities or electron densities can also be",
446 "calculated. For electron densities, a file describing the number of",
447 "electrons for each type of atom should be provided using [TT]-ei[tt].",
448 "It should look like:[BR]",
450 " [TT]atomname = nrelectrons[tt][BR]",
451 " [TT]atomname = nrelectrons[tt][BR]",
452 "The first line contains the number of lines to read from the file.",
453 "There should be one line for each unique atom name in your system.",
454 "The number of electrons for each atom is modified by its atomic",
459 static const char *dens_opt[] =
460 { NULL, "mass", "number", "charge", "electron", NULL };
461 static int axis = 2; /* normal to memb. default z */
462 static const char *axtitle = "Z";
463 static int nslices = 50; /* nr of slices defined */
464 static int ngrps = 1; /* nr. of groups */
465 static gmx_bool bSymmetrize = FALSE;
466 static gmx_bool bCenter = FALSE;
468 { "-d", FALSE, etSTR, {&axtitle},
469 "Take the normal on the membrane in direction X, Y or Z." },
470 { "-sl", FALSE, etINT, {&nslices},
471 "Divide the box in this number of slices." },
472 { "-dens", FALSE, etENUM, {dens_opt},
474 { "-ng", FALSE, etINT, {&ngrps},
475 "Number of groups of which to compute densities." },
476 { "-symm", FALSE, etBOOL, {&bSymmetrize},
477 "Symmetrize the density along the axis, with respect to the center. Useful for bilayers." },
478 { "-center", FALSE, etBOOL, {&bCenter},
479 "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."}
482 const char *bugs[] = {
483 "When calculating electron densities, atomnames are used instead of types. This is bad.",
486 double **density; /* density per slice */
487 real slWidth; /* width of one slice */
488 char **grpname; /* groupnames */
489 int nr_electrons; /* nr. electrons */
490 int *ngx; /* sizes of groups */
491 t_electron *el_tab; /* tabel with nr. of electrons*/
492 t_topology *top; /* topology */
494 atom_id **index; /* indices for all groups */
497 t_filenm fnm[] = { /* files for g_density */
498 { efTRX, "-f", NULL, ffREAD },
499 { efNDX, NULL, NULL, ffOPTRD },
500 { efTPX, NULL, NULL, ffREAD },
501 { efDAT, "-ei", "electrons", ffOPTRD }, /* file with nr. of electrons */
502 { efXVG, "-o", "density", ffWRITE },
505 #define NFILE asize(fnm)
507 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
508 NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs,
514 if (bSymmetrize && !bCenter)
516 fprintf(stderr, "Can not symmetrize without centering. Turning on -center\n");
520 axis = toupper(axtitle[0]) - 'X';
522 top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC); /* read topology file */
523 if (dens_opt[0][0] == 'n')
525 for (i = 0; (i < top->atoms.nr); i++)
527 top->atoms.atom[i].m = 1;
530 else if (dens_opt[0][0] == 'c')
532 for (i = 0; (i < top->atoms.nr); i++)
534 top->atoms.atom[i].m = top->atoms.atom[i].q;
538 snew(grpname, ngrps);
542 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, ngx, index, grpname);
544 if (dens_opt[0][0] == 'e')
546 nr_electrons = get_electrons(&el_tab, ftp2fn(efDAT, NFILE, fnm));
547 fprintf(stderr, "Read %d atomtypes from datafile\n", nr_electrons);
549 calc_electron_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density,
550 &nslices, top, ePBC, axis, ngrps, &slWidth, el_tab,
551 nr_electrons, bCenter, oenv);
555 calc_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density, &nslices, top,
556 ePBC, axis, ngrps, &slWidth, bCenter, oenv);
559 plot_density(density, opt2fn("-o", NFILE, fnm),
560 nslices, ngrps, grpname, slWidth, dens_opt,
563 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy"); /* view xvgr file */