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
64 /****************************************************************************/
65 /* This program calculates the partial density across the box. */
66 /* Peter Tieleman, Mei 1995 */
67 /****************************************************************************/
69 /* used for sorting the list */
70 int compare(void *a, void *b)
72 t_electron *tmp1, *tmp2;
73 tmp1 = (t_electron *)a; tmp2 = (t_electron *)b;
75 return strcmp(tmp1->atomname, tmp2->atomname);
78 int get_electrons(t_electron **eltab, const char *fn)
80 char buffer[256]; /* to read in a line */
81 char tempname[80]; /* buffer to hold name */
85 int nr; /* number of atomstypes to read */
88 if (!(in = ffopen(fn, "r")))
90 gmx_fatal(FARGS, "Couldn't open %s. Exiting.\n", fn);
93 if (NULL == fgets(buffer, 255, in))
95 gmx_fatal(FARGS, "Error reading from file %s", fn);
98 if (sscanf(buffer, "%d", &nr) != 1)
100 gmx_fatal(FARGS, "Invalid number of atomtypes in datafile\n");
105 for (i = 0; i < nr; i++)
107 if (fgets(buffer, 255, in) == NULL)
109 gmx_fatal(FARGS, "reading datafile. Check your datafile.\n");
111 if (sscanf(buffer, "%s = %d", tempname, &tempnr) != 2)
113 gmx_fatal(FARGS, "Invalid line in datafile at line %d\n", i+1);
115 (*eltab)[i].nr_el = tempnr;
116 (*eltab)[i].atomname = strdup(tempname);
121 fprintf(stderr, "Sorting list..\n");
122 qsort ((void*)*eltab, nr, sizeof(t_electron),
123 (int(*)(const void*, const void*))compare);
128 void center_coords(t_atoms *atoms, matrix box, rvec x0[], int axis)
132 rvec com, shift, box_center;
136 for (i = 0; (i < atoms->nr); i++)
138 mm = atoms->atom[i].m;
140 for (m = 0; (m < DIM); m++)
142 com[m] += mm*x0[i][m];
145 for (m = 0; (m < DIM); m++)
149 calc_box_center(ecenterDEF, box, box_center);
150 rvec_sub(box_center, com, shift);
151 shift[axis] -= box_center[axis];
153 for (i = 0; (i < atoms->nr); i++)
155 rvec_dec(x0[i], shift);
159 void calc_electron_density(const char *fn, atom_id **index, int gnx[],
160 double ***slDensity, int *nslices, t_topology *top,
162 int axis, int nr_grps, real *slWidth,
163 t_electron eltab[], int nr, gmx_bool bCenter,
164 const output_env_t oenv)
166 rvec *x0; /* coordinates without pbc */
167 matrix box; /* box (3x3) */
169 int natoms; /* nr. atoms in trj */
171 int i, n, /* loop indices */
172 nr_frames = 0, /* number of frames */
173 slice; /* current slice */
174 t_electron *found; /* found by bsearch */
175 t_electron sought; /* thingie thought by bsearch */
176 gmx_rmpbc_t gpbc = NULL;
181 if (axis < 0 || axis >= DIM)
183 gmx_fatal(FARGS, "Invalid axes. Terminating\n");
186 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
188 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
193 *nslices = (int)(box[axis][axis] * 10); /* default value */
195 fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
197 snew(*slDensity, nr_grps);
198 for (i = 0; i < nr_grps; i++)
200 snew((*slDensity)[i], *nslices);
203 gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
204 /*********** Start processing trajectory ***********/
207 gmx_rmpbc(gpbc, natoms, box, x0);
211 center_coords(&top->atoms, box, x0, axis);
214 *slWidth = box[axis][axis]/(*nslices);
215 invvol = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
217 for (n = 0; n < nr_grps; n++)
219 for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
221 z = x0[index[n][i]][axis];
224 z += box[axis][axis];
226 while (z > box[axis][axis])
228 z -= box[axis][axis];
231 /* determine which slice atom is in */
232 slice = (z / (*slWidth));
234 sought.atomname = strdup(*(top->atoms.atomname[index[n][i]]));
236 /* now find the number of electrons. This is not efficient. */
237 found = (t_electron *)
238 bsearch((const void *)&sought,
239 (const void *)eltab, nr, sizeof(t_electron),
240 (int(*)(const void*, const void*))compare);
244 fprintf(stderr, "Couldn't find %s. Add it to the .dat file\n",
245 *(top->atoms.atomname[index[n][i]]));
249 (*slDensity)[n][slice] += (found->nr_el -
250 top->atoms.atom[index[n][i]].q)*invvol;
252 free(sought.atomname);
257 while (read_next_x(oenv, status, &t, x0, box));
258 gmx_rmpbc_done(gpbc);
260 /*********** done with status file **********/
263 /* slDensity now contains the total number of electrons per slice, summed
264 over all frames. Now divide by nr_frames and volume of slice
267 fprintf(stderr, "\nRead %d frames from trajectory. Counting electrons\n",
270 for (n = 0; n < nr_grps; n++)
272 for (i = 0; i < *nslices; i++)
274 (*slDensity)[n][i] /= nr_frames;
278 sfree(x0); /* free memory used by coordinate array */
281 void calc_density(const char *fn, atom_id **index, int gnx[],
282 double ***slDensity, int *nslices, t_topology *top, int ePBC,
283 int axis, int nr_grps, real *slWidth, gmx_bool bCenter,
284 const output_env_t oenv)
286 rvec *x0; /* coordinates without pbc */
287 matrix box; /* box (3x3) */
289 int natoms; /* nr. atoms in trj */
291 int **slCount, /* nr. of atoms in one slice for a group */
292 i, j, n, /* loop indices */
295 nr_frames = 0, /* number of frames */
296 slice; /* current slice */
299 char *buf; /* for tmp. keeping atomname */
300 gmx_rmpbc_t gpbc = NULL;
302 if (axis < 0 || axis >= DIM)
304 gmx_fatal(FARGS, "Invalid axes. Terminating\n");
307 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
309 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
314 *nslices = (int)(box[axis][axis] * 10); /* default value */
315 fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
318 snew(*slDensity, nr_grps);
319 for (i = 0; i < nr_grps; i++)
321 snew((*slDensity)[i], *nslices);
324 gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
325 /*********** Start processing trajectory ***********/
328 gmx_rmpbc(gpbc, natoms, box, x0);
332 center_coords(&top->atoms, box, x0, axis);
335 *slWidth = box[axis][axis]/(*nslices);
336 invvol = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
339 for (n = 0; n < nr_grps; n++)
341 for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
343 z = x0[index[n][i]][axis];
346 z += box[axis][axis];
348 while (z > box[axis][axis])
350 z -= box[axis][axis];
353 /* determine which slice atom is in */
354 slice = (int)(z / (*slWidth));
355 (*slDensity)[n][slice] += top->atoms.atom[index[n][i]].m*invvol;
360 while (read_next_x(oenv, status, &t, x0, box));
361 gmx_rmpbc_done(gpbc);
363 /*********** done with status file **********/
366 /* slDensity now contains the total mass per slice, summed over all
367 frames. Now divide by nr_frames and volume of slice
370 fprintf(stderr, "\nRead %d frames from trajectory. Calculating density\n",
373 for (n = 0; n < nr_grps; n++)
375 for (i = 0; i < *nslices; i++)
377 (*slDensity)[n][i] /= nr_frames;
381 sfree(x0); /* free memory used by coordinate array */
384 void plot_density(double *slDensity[], const char *afile, int nslices,
385 int nr_grps, char *grpname[], real slWidth,
386 const char **dens_opt,
387 gmx_bool bSymmetrize, const output_env_t oenv)
390 const char *ylabel = NULL;
394 switch (dens_opt[0][0])
396 case 'm': ylabel = "Density (kg m\\S-3\\N)"; break;
397 case 'n': ylabel = "Number density (nm\\S-3\\N)"; break;
398 case 'c': ylabel = "Charge density (e nm\\S-3\\N)"; break;
399 case 'e': ylabel = "Electron density (e nm\\S-3\\N)"; break;
402 den = xvgropen(afile, "Partial densities", "Box (nm)", ylabel, oenv);
404 xvgr_legend(den, nr_grps, (const char**)grpname, oenv);
406 for (slice = 0; (slice < nslices); slice++)
408 fprintf(den, "%12g ", slice * slWidth);
409 for (n = 0; (n < nr_grps); n++)
413 ddd = (slDensity[n][slice]+slDensity[n][nslices-slice-1])*0.5;
417 ddd = slDensity[n][slice];
419 if (dens_opt[0][0] == 'm')
421 fprintf(den, " %12g", ddd*AMU/(NANO*NANO*NANO));
425 fprintf(den, " %12g", ddd);
434 int gmx_density(int argc, char *argv[])
436 const char *desc[] = {
437 "Compute partial densities across the box, using an index file.[PAR]",
438 "For the total density of NPT simulations, use [TT]g_energy[tt] instead.",
440 "Densities are in kg/m^3, and number densities or electron densities can also be",
441 "calculated. For electron densities, a file describing the number of",
442 "electrons for each type of atom should be provided using [TT]-ei[tt].",
443 "It should look like:[BR]",
445 " [TT]atomname = nrelectrons[tt][BR]",
446 " [TT]atomname = nrelectrons[tt][BR]",
447 "The first line contains the number of lines to read from the file.",
448 "There should be one line for each unique atom name in your system.",
449 "The number of electrons for each atom is modified by its atomic",
454 static const char *dens_opt[] =
455 { NULL, "mass", "number", "charge", "electron", NULL };
456 static int axis = 2; /* normal to memb. default z */
457 static const char *axtitle = "Z";
458 static int nslices = 50; /* nr of slices defined */
459 static int ngrps = 1; /* nr. of groups */
460 static gmx_bool bSymmetrize = FALSE;
461 static gmx_bool bCenter = FALSE;
463 { "-d", FALSE, etSTR, {&axtitle},
464 "Take the normal on the membrane in direction X, Y or Z." },
465 { "-sl", FALSE, etINT, {&nslices},
466 "Divide the box in this number of slices." },
467 { "-dens", FALSE, etENUM, {dens_opt},
469 { "-ng", FALSE, etINT, {&ngrps},
470 "Number of groups of which to compute densities." },
471 { "-symm", FALSE, etBOOL, {&bSymmetrize},
472 "Symmetrize the density along the axis, with respect to the center. Useful for bilayers." },
473 { "-center", FALSE, etBOOL, {&bCenter},
474 "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."}
477 const char *bugs[] = {
478 "When calculating electron densities, atomnames are used instead of types. This is bad.",
481 double **density; /* density per slice */
482 real slWidth; /* width of one slice */
483 char **grpname; /* groupnames */
484 int nr_electrons; /* nr. electrons */
485 int *ngx; /* sizes of groups */
486 t_electron *el_tab; /* tabel with nr. of electrons*/
487 t_topology *top; /* topology */
489 atom_id **index; /* indices for all groups */
492 t_filenm fnm[] = { /* files for g_density */
493 { efTRX, "-f", NULL, ffREAD },
494 { efNDX, NULL, NULL, ffOPTRD },
495 { efTPX, NULL, NULL, ffREAD },
496 { efDAT, "-ei", "electrons", ffOPTRD }, /* file with nr. of electrons */
497 { efXVG, "-o", "density", ffWRITE },
500 #define NFILE asize(fnm)
502 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
503 NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs,
509 if (bSymmetrize && !bCenter)
511 fprintf(stderr, "Can not symmetrize without centering. Turning on -center\n");
515 axis = toupper(axtitle[0]) - 'X';
517 top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC); /* read topology file */
518 if (dens_opt[0][0] == 'n')
520 for (i = 0; (i < top->atoms.nr); i++)
522 top->atoms.atom[i].m = 1;
525 else if (dens_opt[0][0] == 'c')
527 for (i = 0; (i < top->atoms.nr); i++)
529 top->atoms.atom[i].m = top->atoms.atom[i].q;
533 snew(grpname, ngrps);
537 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, ngx, index, grpname);
539 if (dens_opt[0][0] == 'e')
541 nr_electrons = get_electrons(&el_tab, ftp2fn(efDAT, NFILE, fnm));
542 fprintf(stderr, "Read %d atomtypes from datafile\n", nr_electrons);
544 calc_electron_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density,
545 &nslices, top, ePBC, axis, ngrps, &slWidth, el_tab,
546 nr_electrons, bCenter, oenv);
550 calc_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density, &nslices, top,
551 ePBC, axis, ngrps, &slWidth, bCenter, oenv);
554 plot_density(density, opt2fn("-o", NFILE, fnm),
555 nslices, ngrps, grpname, slWidth, dens_opt,
558 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy"); /* view xvgr file */