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
65 /****************************************************************************/
66 /* This program calculates the partial density across the box. */
67 /* Peter Tieleman, Mei 1995 */
68 /****************************************************************************/
70 /* used for sorting the list */
71 int compare(void *a, void *b)
73 t_electron *tmp1, *tmp2;
74 tmp1 = (t_electron *)a; tmp2 = (t_electron *)b;
76 return strcmp(tmp1->atomname, tmp2->atomname);
79 int get_electrons(t_electron **eltab, const char *fn)
81 char buffer[256]; /* to read in a line */
82 char tempname[80]; /* buffer to hold name */
86 int nr; /* number of atomstypes to read */
89 if (!(in = ffopen(fn, "r")))
91 gmx_fatal(FARGS, "Couldn't open %s. Exiting.\n", fn);
94 if (NULL == fgets(buffer, 255, in))
96 gmx_fatal(FARGS, "Error reading from file %s", fn);
99 if (sscanf(buffer, "%d", &nr) != 1)
101 gmx_fatal(FARGS, "Invalid number of atomtypes in datafile\n");
106 for (i = 0; i < nr; i++)
108 if (fgets(buffer, 255, in) == NULL)
110 gmx_fatal(FARGS, "reading datafile. Check your datafile.\n");
112 if (sscanf(buffer, "%s = %d", tempname, &tempnr) != 2)
114 gmx_fatal(FARGS, "Invalid line in datafile at line %d\n", i+1);
116 (*eltab)[i].nr_el = tempnr;
117 (*eltab)[i].atomname = strdup(tempname);
122 fprintf(stderr, "Sorting list..\n");
123 qsort ((void*)*eltab, nr, sizeof(t_electron),
124 (int(*)(const void*, const void*))compare);
129 void center_coords(t_atoms *atoms, matrix box, rvec x0[], int axis)
133 rvec com, shift, box_center;
137 for (i = 0; (i < atoms->nr); i++)
139 mm = atoms->atom[i].m;
141 for (m = 0; (m < DIM); m++)
143 com[m] += mm*x0[i][m];
146 for (m = 0; (m < DIM); m++)
150 calc_box_center(ecenterDEF, box, box_center);
151 rvec_sub(box_center, com, shift);
152 shift[axis] -= box_center[axis];
154 for (i = 0; (i < atoms->nr); i++)
156 rvec_dec(x0[i], shift);
160 void calc_electron_density(const char *fn, atom_id **index, int gnx[],
161 double ***slDensity, int *nslices, t_topology *top,
163 int axis, int nr_grps, real *slWidth,
164 t_electron eltab[], int nr, gmx_bool bCenter,
165 const output_env_t oenv)
167 rvec *x0; /* coordinates without pbc */
168 matrix box; /* box (3x3) */
170 int natoms; /* nr. atoms in trj */
172 int i, n, /* loop indices */
173 nr_frames = 0, /* number of frames */
174 slice; /* current slice */
175 t_electron *found; /* found by bsearch */
176 t_electron sought; /* thingie thought by bsearch */
177 gmx_rmpbc_t gpbc = NULL;
182 if (axis < 0 || axis >= DIM)
184 gmx_fatal(FARGS, "Invalid axes. Terminating\n");
187 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
189 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
194 *nslices = (int)(box[axis][axis] * 10); /* default value */
196 fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
198 snew(*slDensity, nr_grps);
199 for (i = 0; i < nr_grps; i++)
201 snew((*slDensity)[i], *nslices);
204 gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr, box);
205 /*********** Start processing trajectory ***********/
208 gmx_rmpbc(gpbc, natoms, box, x0);
212 center_coords(&top->atoms, box, x0, axis);
215 *slWidth = box[axis][axis]/(*nslices);
216 invvol = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
218 for (n = 0; n < nr_grps; n++)
220 for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
222 z = x0[index[n][i]][axis];
225 z += box[axis][axis];
227 while (z > box[axis][axis])
229 z -= box[axis][axis];
232 /* determine which slice atom is in */
233 slice = (z / (*slWidth));
235 sought.atomname = strdup(*(top->atoms.atomname[index[n][i]]));
237 /* now find the number of electrons. This is not efficient. */
238 found = (t_electron *)
239 bsearch((const void *)&sought,
240 (const void *)eltab, nr, sizeof(t_electron),
241 (int(*)(const void*, const void*))compare);
245 fprintf(stderr, "Couldn't find %s. Add it to the .dat file\n",
246 *(top->atoms.atomname[index[n][i]]));
250 (*slDensity)[n][slice] += (found->nr_el -
251 top->atoms.atom[index[n][i]].q)*invvol;
253 free(sought.atomname);
258 while (read_next_x(oenv, status, &t, natoms, x0, box));
259 gmx_rmpbc_done(gpbc);
261 /*********** done with status file **********/
264 /* slDensity now contains the total number of electrons per slice, summed
265 over all frames. Now divide by nr_frames and volume of slice
268 fprintf(stderr, "\nRead %d frames from trajectory. Counting electrons\n",
271 for (n = 0; n < nr_grps; n++)
273 for (i = 0; i < *nslices; i++)
275 (*slDensity)[n][i] /= nr_frames;
279 sfree(x0); /* free memory used by coordinate array */
282 void calc_density(const char *fn, atom_id **index, int gnx[],
283 double ***slDensity, int *nslices, t_topology *top, int ePBC,
284 int axis, int nr_grps, real *slWidth, gmx_bool bCenter,
285 const output_env_t oenv)
287 rvec *x0; /* coordinates without pbc */
288 matrix box; /* box (3x3) */
290 int natoms; /* nr. atoms in trj */
292 int **slCount, /* nr. of atoms in one slice for a group */
293 i, j, n, /* loop indices */
296 nr_frames = 0, /* number of frames */
297 slice; /* current slice */
300 char *buf; /* for tmp. keeping atomname */
301 gmx_rmpbc_t gpbc = NULL;
303 if (axis < 0 || axis >= DIM)
305 gmx_fatal(FARGS, "Invalid axes. Terminating\n");
308 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
310 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
315 *nslices = (int)(box[axis][axis] * 10); /* default value */
316 fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
319 snew(*slDensity, nr_grps);
320 for (i = 0; i < nr_grps; i++)
322 snew((*slDensity)[i], *nslices);
325 gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr, box);
326 /*********** Start processing trajectory ***********/
329 gmx_rmpbc(gpbc, natoms, box, x0);
333 center_coords(&top->atoms, box, x0, axis);
336 *slWidth = box[axis][axis]/(*nslices);
337 invvol = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
340 for (n = 0; n < nr_grps; n++)
342 for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
344 z = x0[index[n][i]][axis];
347 z += box[axis][axis];
349 while (z > box[axis][axis])
351 z -= box[axis][axis];
354 /* determine which slice atom is in */
355 slice = (int)(z / (*slWidth));
356 (*slDensity)[n][slice] += top->atoms.atom[index[n][i]].m*invvol;
361 while (read_next_x(oenv, status, &t, natoms, x0, box));
362 gmx_rmpbc_done(gpbc);
364 /*********** done with status file **********/
367 /* slDensity now contains the total mass per slice, summed over all
368 frames. Now divide by nr_frames and volume of slice
371 fprintf(stderr, "\nRead %d frames from trajectory. Calculating density\n",
374 for (n = 0; n < nr_grps; n++)
376 for (i = 0; i < *nslices; i++)
378 (*slDensity)[n][i] /= nr_frames;
382 sfree(x0); /* free memory used by coordinate array */
385 void plot_density(double *slDensity[], const char *afile, int nslices,
386 int nr_grps, char *grpname[], real slWidth,
387 const char **dens_opt,
388 gmx_bool bSymmetrize, const output_env_t oenv)
391 const char *ylabel = NULL;
395 switch (dens_opt[0][0])
397 case 'm': ylabel = "Density (kg m\\S-3\\N)"; break;
398 case 'n': ylabel = "Number density (nm\\S-3\\N)"; break;
399 case 'c': ylabel = "Charge density (e nm\\S-3\\N)"; break;
400 case 'e': ylabel = "Electron density (e nm\\S-3\\N)"; break;
403 den = xvgropen(afile, "Partial densities", "Box (nm)", ylabel, oenv);
405 xvgr_legend(den, nr_grps, (const char**)grpname, oenv);
407 for (slice = 0; (slice < nslices); slice++)
409 fprintf(den, "%12g ", slice * slWidth);
410 for (n = 0; (n < nr_grps); n++)
414 ddd = (slDensity[n][slice]+slDensity[n][nslices-slice-1])*0.5;
418 ddd = slDensity[n][slice];
420 if (dens_opt[0][0] == 'm')
422 fprintf(den, " %12g", ddd*AMU/(NANO*NANO*NANO));
426 fprintf(den, " %12g", ddd);
435 int gmx_density(int argc, char *argv[])
437 const char *desc[] = {
438 "Compute partial densities across the box, using an index file.[PAR]",
439 "For the total density of NPT simulations, use [TT]g_energy[tt] instead.",
441 "Densities are in kg/m^3, and number densities or electron densities can also be",
442 "calculated. For electron densities, a file describing the number of",
443 "electrons for each type of atom should be provided using [TT]-ei[tt].",
444 "It should look like:[BR]",
446 " [TT]atomname = nrelectrons[tt][BR]",
447 " [TT]atomname = nrelectrons[tt][BR]",
448 "The first line contains the number of lines to read from the file.",
449 "There should be one line for each unique atom name in your system.",
450 "The number of electrons for each atom is modified by its atomic",
455 static const char *dens_opt[] =
456 { NULL, "mass", "number", "charge", "electron", NULL };
457 static int axis = 2; /* normal to memb. default z */
458 static const char *axtitle = "Z";
459 static int nslices = 50; /* nr of slices defined */
460 static int ngrps = 1; /* nr. of groups */
461 static gmx_bool bSymmetrize = FALSE;
462 static gmx_bool bCenter = FALSE;
464 { "-d", FALSE, etSTR, {&axtitle},
465 "Take the normal on the membrane in direction X, Y or Z." },
466 { "-sl", FALSE, etINT, {&nslices},
467 "Divide the box in this number of slices." },
468 { "-dens", FALSE, etENUM, {dens_opt},
470 { "-ng", FALSE, etINT, {&ngrps},
471 "Number of groups of which to compute densities." },
472 { "-symm", FALSE, etBOOL, {&bSymmetrize},
473 "Symmetrize the density along the axis, with respect to the center. Useful for bilayers." },
474 { "-center", FALSE, etBOOL, {&bCenter},
475 "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."}
478 const char *bugs[] = {
479 "When calculating electron densities, atomnames are used instead of types. This is bad.",
482 double **density; /* density per slice */
483 real slWidth; /* width of one slice */
484 char **grpname; /* groupnames */
485 int nr_electrons; /* nr. electrons */
486 int *ngx; /* sizes of groups */
487 t_electron *el_tab; /* tabel with nr. of electrons*/
488 t_topology *top; /* topology */
490 atom_id **index; /* indices for all groups */
493 t_filenm fnm[] = { /* files for g_density */
494 { efTRX, "-f", NULL, ffREAD },
495 { efNDX, NULL, NULL, ffOPTRD },
496 { efTPX, NULL, NULL, ffREAD },
497 { efDAT, "-ei", "electrons", ffOPTRD }, /* file with nr. of electrons */
498 { efXVG, "-o", "density", ffWRITE },
501 #define NFILE asize(fnm)
503 parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
504 NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs,
507 if (bSymmetrize && !bCenter)
509 fprintf(stderr, "Can not symmetrize without centering. Turning on -center\n");
513 axis = toupper(axtitle[0]) - 'X';
515 top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC); /* read topology file */
516 if (dens_opt[0][0] == 'n')
518 for (i = 0; (i < top->atoms.nr); i++)
520 top->atoms.atom[i].m = 1;
523 else if (dens_opt[0][0] == 'c')
525 for (i = 0; (i < top->atoms.nr); i++)
527 top->atoms.atom[i].m = top->atoms.atom[i].q;
531 snew(grpname, ngrps);
535 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, ngx, index, grpname);
537 if (dens_opt[0][0] == 'e')
539 nr_electrons = get_electrons(&el_tab, ftp2fn(efDAT, NFILE, fnm));
540 fprintf(stderr, "Read %d atomtypes from datafile\n", nr_electrons);
542 calc_electron_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density,
543 &nslices, top, ePBC, axis, ngrps, &slWidth, el_tab,
544 nr_electrons, bCenter, oenv);
548 calc_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density, &nslices, top,
549 ePBC, axis, ngrps, &slWidth, bCenter, oenv);
552 plot_density(density, opt2fn("-o", NFILE, fnm),
553 nslices, ngrps, grpname, slWidth, dens_opt,
556 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy"); /* view xvgr file */