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.
44 #include "gromacs/utility/cstringutil.h"
45 #include "gromacs/legacyheaders/typedefs.h"
46 #include "gromacs/legacyheaders/macros.h"
48 #include "gromacs/legacyheaders/viewit.h"
49 #include "gromacs/pbcutil/pbc.h"
50 #include "gromacs/topology/index.h"
51 #include "gromacs/fileio/tpxio.h"
52 #include "gromacs/fileio/trxio.h"
53 #include "gromacs/math/units.h"
55 #include "gromacs/legacyheaders/macros.h"
57 #include "gromacs/commandline/pargs.h"
58 #include "gromacs/fileio/xvgr.h"
59 #include "gromacs/math/vec.h"
60 #include "gromacs/pbcutil/rmpbc.h"
61 #include "gromacs/utility/fatalerror.h"
62 #include "gromacs/utility/futil.h"
63 #include "gromacs/utility/smalloc.h"
70 /****************************************************************************/
71 /* This program calculates the partial density across the box. */
72 /* Peter Tieleman, Mei 1995 */
73 /****************************************************************************/
75 /* used for sorting the list */
76 int compare(void *a, void *b)
78 t_electron *tmp1, *tmp2;
79 tmp1 = (t_electron *)a; tmp2 = (t_electron *)b;
81 return strcmp(tmp1->atomname, tmp2->atomname);
84 int get_electrons(t_electron **eltab, const char *fn)
86 char buffer[256]; /* to read in a line */
87 char tempname[80]; /* buffer to hold name */
91 int nr; /* number of atomstypes to read */
94 if (!(in = gmx_ffopen(fn, "r")))
96 gmx_fatal(FARGS, "Couldn't open %s. Exiting.\n", fn);
99 if (NULL == fgets(buffer, 255, in))
101 gmx_fatal(FARGS, "Error reading from file %s", fn);
104 if (sscanf(buffer, "%d", &nr) != 1)
106 gmx_fatal(FARGS, "Invalid number of atomtypes in datafile\n");
111 for (i = 0; i < nr; i++)
113 if (fgets(buffer, 255, in) == NULL)
115 gmx_fatal(FARGS, "reading datafile. Check your datafile.\n");
117 if (sscanf(buffer, "%s = %d", tempname, &tempnr) != 2)
119 gmx_fatal(FARGS, "Invalid line in datafile at line %d\n", i+1);
121 (*eltab)[i].nr_el = tempnr;
122 (*eltab)[i].atomname = gmx_strdup(tempname);
127 fprintf(stderr, "Sorting list..\n");
128 qsort ((void*)*eltab, nr, sizeof(t_electron),
129 (int(*)(const void*, const void*))compare);
134 void center_coords(t_atoms *atoms, atom_id *index_center, int ncenter,
135 matrix box, rvec x0[])
139 rvec com, shift, box_center;
143 for (k = 0; (k < ncenter); k++)
148 gmx_fatal(FARGS, "Index %d refers to atom %d, which is larger than natoms (%d).",
149 k+1, i+1, atoms->nr);
151 mm = atoms->atom[i].m;
153 for (m = 0; (m < DIM); m++)
155 com[m] += mm*x0[i][m];
158 for (m = 0; (m < DIM); m++)
162 calc_box_center(ecenterDEF, box, box_center);
163 rvec_sub(box_center, com, shift);
165 /* Important - while the center was calculated based on a group, we should move all atoms */
166 for (i = 0; (i < atoms->nr); i++)
168 rvec_dec(x0[i], shift);
172 void calc_electron_density(const char *fn, atom_id **index, int gnx[],
173 double ***slDensity, int *nslices, t_topology *top,
175 int axis, int nr_grps, real *slWidth,
176 t_electron eltab[], int nr, gmx_bool bCenter,
177 atom_id *index_center, int ncenter,
178 gmx_bool bRelative, const output_env_t oenv)
180 rvec *x0; /* coordinates without pbc */
181 matrix box; /* box (3x3) */
183 int natoms; /* nr. atoms in trj */
185 int i, n, /* loop indices */
186 nr_frames = 0, /* number of frames */
187 slice; /* current slice */
188 t_electron *found; /* found by bsearch */
189 t_electron sought; /* thingie thought by bsearch */
191 gmx_rmpbc_t gpbc = NULL;
196 if (axis < 0 || axis >= DIM)
198 gmx_fatal(FARGS, "Invalid axes. Terminating\n");
201 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
203 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
210 *nslices = (int)(box[axis][axis] * 10); /* default value */
211 fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
214 snew(*slDensity, nr_grps);
215 for (i = 0; i < nr_grps; i++)
217 snew((*slDensity)[i], *nslices);
220 gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
221 /*********** Start processing trajectory ***********/
224 gmx_rmpbc(gpbc, natoms, box, x0);
226 /* Translate atoms so the com of the center-group is in the
227 * box geometrical center.
231 center_coords(&top->atoms, index_center, ncenter, box, x0);
234 invvol = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
238 *slWidth = 1.0/(*nslices);
243 *slWidth = box[axis][axis]/(*nslices);
244 boxSz = box[axis][axis];
247 aveBox += box[axis][axis];
249 for (n = 0; n < nr_grps; n++)
251 for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
253 z = x0[index[n][i]][axis];
256 z += box[axis][axis];
258 while (z > box[axis][axis])
260 z -= box[axis][axis];
265 z = z/box[axis][axis];
268 /* determine which slice atom is in */
271 slice = floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2;
275 slice = (z / (*slWidth));
278 sought.atomname = gmx_strdup(*(top->atoms.atomname[index[n][i]]));
280 /* now find the number of electrons. This is not efficient. */
281 found = (t_electron *)
282 bsearch((const void *)&sought,
283 (const void *)eltab, nr, sizeof(t_electron),
284 (int(*)(const void*, const void*))compare);
288 fprintf(stderr, "Couldn't find %s. Add it to the .dat file\n",
289 *(top->atoms.atomname[index[n][i]]));
293 (*slDensity)[n][slice] += (found->nr_el -
294 top->atoms.atom[index[n][i]].q)*invvol;
296 free(sought.atomname);
301 while (read_next_x(oenv, status, &t, x0, box));
302 gmx_rmpbc_done(gpbc);
304 /*********** done with status file **********/
307 /* slDensity now contains the total number of electrons per slice, summed
308 over all frames. Now divide by nr_frames and volume of slice
311 fprintf(stderr, "\nRead %d frames from trajectory. Counting electrons\n",
317 *slWidth = aveBox/(*nslices);
320 for (n = 0; n < nr_grps; n++)
322 for (i = 0; i < *nslices; i++)
324 (*slDensity)[n][i] /= nr_frames;
328 sfree(x0); /* free memory used by coordinate array */
331 void calc_density(const char *fn, atom_id **index, int gnx[],
332 double ***slDensity, int *nslices, t_topology *top, int ePBC,
333 int axis, int nr_grps, real *slWidth, gmx_bool bCenter,
334 atom_id *index_center, int ncenter,
335 gmx_bool bRelative, const output_env_t oenv)
337 rvec *x0; /* coordinates without pbc */
338 matrix box; /* box (3x3) */
340 int natoms; /* nr. atoms in trj */
342 int **slCount, /* nr. of atoms in one slice for a group */
343 i, j, n, /* loop indices */
345 nr_frames = 0, /* number of frames */
346 slice; /* current slice */
350 char *buf; /* for tmp. keeping atomname */
351 gmx_rmpbc_t gpbc = NULL;
353 if (axis < 0 || axis >= DIM)
355 gmx_fatal(FARGS, "Invalid axes. Terminating\n");
358 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
360 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
367 *nslices = (int)(box[axis][axis] * 10); /* default value */
368 fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
371 snew(*slDensity, nr_grps);
372 for (i = 0; i < nr_grps; i++)
374 snew((*slDensity)[i], *nslices);
377 gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
378 /*********** Start processing trajectory ***********/
381 gmx_rmpbc(gpbc, natoms, box, x0);
383 /* Translate atoms so the com of the center-group is in the
384 * box geometrical center.
388 center_coords(&top->atoms, index_center, ncenter, box, x0);
391 invvol = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
395 *slWidth = 1.0/(*nslices);
400 *slWidth = box[axis][axis]/(*nslices);
401 boxSz = box[axis][axis];
404 aveBox += box[axis][axis];
406 for (n = 0; n < nr_grps; n++)
408 for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
410 z = x0[index[n][i]][axis];
413 z += box[axis][axis];
415 while (z > box[axis][axis])
417 z -= box[axis][axis];
422 z = z/box[axis][axis];
425 /* determine which slice atom is in */
428 slice = floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2;
432 slice = floor(z / (*slWidth));
435 /* Slice should already be 0<=slice<nslices, but we just make
436 * sure we are not hit by IEEE rounding errors since we do
437 * math operations after applying PBC above.
443 else if (slice >= *nslices)
448 (*slDensity)[n][slice] += top->atoms.atom[index[n][i]].m*invvol;
453 while (read_next_x(oenv, status, &t, x0, box));
454 gmx_rmpbc_done(gpbc);
456 /*********** done with status file **********/
459 /* slDensity now contains the total mass per slice, summed over all
460 frames. Now divide by nr_frames and volume of slice
463 fprintf(stderr, "\nRead %d frames from trajectory. Calculating density\n",
469 *slWidth = aveBox/(*nslices);
472 for (n = 0; n < nr_grps; n++)
474 for (i = 0; i < *nslices; i++)
476 (*slDensity)[n][i] /= nr_frames;
480 sfree(x0); /* free memory used by coordinate array */
483 void plot_density(double *slDensity[], const char *afile, int nslices,
484 int nr_grps, char *grpname[], real slWidth,
485 const char **dens_opt,
486 gmx_bool bCenter, gmx_bool bRelative, gmx_bool bSymmetrize,
487 const output_env_t oenv)
490 const char *title = NULL;
491 const char *xlabel = NULL;
492 const char *ylabel = NULL;
497 title = bSymmetrize ? "Symmetrized partial density" : "Partial density";
502 "Average relative position from center (nm)" :
503 "Relative position from center (nm)";
507 xlabel = bRelative ? "Average coordinate (nm)" : "Coordinate (nm)";
510 switch (dens_opt[0][0])
512 case 'm': ylabel = "Density (kg m\\S-3\\N)"; break;
513 case 'n': ylabel = "Number density (nm\\S-3\\N)"; break;
514 case 'c': ylabel = "Charge density (e nm\\S-3\\N)"; break;
515 case 'e': ylabel = "Electron density (e nm\\S-3\\N)"; break;
518 den = xvgropen(afile,
519 title, xlabel, ylabel, oenv);
521 xvgr_legend(den, nr_grps, (const char**)grpname, oenv);
523 for (slice = 0; (slice < nslices); slice++)
527 axispos = (slice - nslices/2.0 + 0.5) * slWidth;
531 axispos = (slice + 0.5) * slWidth;
533 fprintf(den, "%12g ", axispos);
534 for (n = 0; (n < nr_grps); n++)
538 ddd = (slDensity[n][slice]+slDensity[n][nslices-slice-1])*0.5;
542 ddd = slDensity[n][slice];
544 if (dens_opt[0][0] == 'm')
546 fprintf(den, " %12g", ddd*AMU/(NANO*NANO*NANO));
550 fprintf(den, " %12g", ddd);
559 int gmx_density(int argc, char *argv[])
561 const char *desc[] = {
562 "[THISMODULE] computes partial densities across the box, using an index file.[PAR]",
563 "For the total density of NPT simulations, use [gmx-energy] instead.",
566 "Option [TT]-center[tt] performs the histogram binning relative to the center",
567 "of an arbitrary group, in absolute box coordinates. If you are calculating",
568 "profiles along the Z axis box dimension bZ, output would be from -bZ/2 to",
569 "bZ/2 if you center based on the entire system.",
570 "Note that this behaviour has changed in Gromacs 5.0; earlier versions",
571 "merely performed a static binning in (0,bZ) and shifted the output. Now",
572 "we compute the center for each frame and bin in (-bZ/2,bZ/2).[PAR]",
574 "Option [TT]-symm[tt] symmetrizes the output around the center. This will",
575 "automatically turn on [TT]-center[tt] too.",
577 "Option [TT]-relative[tt] performs the binning in relative instead of absolute",
578 "box coordinates, and scales the final output with the average box dimension",
579 "along the output axis. This can be used in combination with [TT]-center[tt].[PAR]",
581 "Densities are in kg/m^3, and number densities or electron densities can also be",
582 "calculated. For electron densities, a file describing the number of",
583 "electrons for each type of atom should be provided using [TT]-ei[tt].",
584 "It should look like:[BR]",
586 " [TT]atomname = nrelectrons[tt][BR]",
587 " [TT]atomname = nrelectrons[tt][BR]",
588 "The first line contains the number of lines to read from the file.",
589 "There should be one line for each unique atom name in your system.",
590 "The number of electrons for each atom is modified by its atomic",
591 "partial charge.[PAR]",
593 "IMPORTANT CONSIDERATIONS FOR BILAYERS[PAR]",
594 "One of the most common usage scenarios is to calculate the density of various",
595 "groups across a lipid bilayer, typically with the z axis being the normal",
596 "direction. For short simulations, small systems, and fixed box sizes this",
597 "will work fine, but for the more general case lipid bilayers can be complicated.",
598 "The first problem that while both proteins and lipids have low volume",
599 "compressibility, lipids have quite high area compressiblity. This means the",
600 "shape of the box (thickness and area/lipid) will fluctuate substantially even",
601 "for a fully relaxed system. Since Gromacs places the box between the origin",
602 "and positive coordinates, this in turn means that a bilayer centered in the",
603 "box will move a bit up/down due to these fluctuations, and smear out your",
604 "profile. The easiest way to fix this (if you want pressure coupling) is",
605 "to use the [TT]-center[tt] option that calculates the density profile with",
606 "respect to the center of the box. Note that you can still center on the",
607 "bilayer part even if you have a complex non-symmetric system with a bilayer",
608 "and, say, membrane proteins - then our output will simply have more values",
609 "on one side of the (center) origin reference.[PAR]",
611 "Even the centered calculation will lead to some smearing out the output",
612 "profiles, as lipids themselves are compressed and expanded. In most cases",
613 "you probably want this (since it corresponds to macroscopic experiments),",
614 "but if you want to look at molecular details you can use the [TT]-relative[tt]",
615 "option to attempt to remove even more of the effects of volume fluctuations.[PAR]",
617 "Finally, large bilayers that are not subject to a surface tension will exhibit",
618 "undulatory fluctuations, where there are 'waves' forming in the system.",
619 "This is a fundamental property of the biological system, and if you are",
620 "comparing against experiments you likely want to include the undulation",
626 static const char *dens_opt[] =
627 { NULL, "mass", "number", "charge", "electron", NULL };
628 static int axis = 2; /* normal to memb. default z */
629 static const char *axtitle = "Z";
630 static int nslices = 50; /* nr of slices defined */
631 static int ngrps = 1; /* nr. of groups */
632 static gmx_bool bSymmetrize = FALSE;
633 static gmx_bool bCenter = FALSE;
634 static gmx_bool bRelative = FALSE;
637 { "-d", FALSE, etSTR, {&axtitle},
638 "Take the normal on the membrane in direction X, Y or Z." },
639 { "-sl", FALSE, etINT, {&nslices},
640 "Divide the box in this number of slices." },
641 { "-dens", FALSE, etENUM, {dens_opt},
643 { "-ng", FALSE, etINT, {&ngrps},
644 "Number of groups of which to compute densities." },
645 { "-center", FALSE, etBOOL, {&bCenter},
646 "Perform the binning relative to the center of the (changing) box. Useful for bilayers." },
647 { "-symm", FALSE, etBOOL, {&bSymmetrize},
648 "Symmetrize the density along the axis, with respect to the center. Useful for bilayers." },
649 { "-relative", FALSE, etBOOL, {&bRelative},
650 "Use relative coordinates for changing boxes and scale output by average dimensions." }
653 const char *bugs[] = {
654 "When calculating electron densities, atomnames are used instead of types. This is bad.",
657 double **density; /* density per slice */
658 real slWidth; /* width of one slice */
659 char *grpname_center; /* centering group name */
660 char **grpname; /* groupnames */
661 int nr_electrons; /* nr. electrons */
662 int ncenter; /* size of centering group */
663 int *ngx; /* sizes of groups */
664 t_electron *el_tab; /* tabel with nr. of electrons*/
665 t_topology *top; /* topology */
667 atom_id *index_center; /* index for centering group */
668 atom_id **index; /* indices for all groups */
671 t_filenm fnm[] = { /* files for g_density */
672 { efTRX, "-f", NULL, ffREAD },
673 { efNDX, NULL, NULL, ffOPTRD },
674 { efTPR, NULL, NULL, ffREAD },
675 { efDAT, "-ei", "electrons", ffOPTRD }, /* file with nr. of electrons */
676 { efXVG, "-o", "density", ffWRITE },
679 #define NFILE asize(fnm)
681 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
682 NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs,
688 if (bSymmetrize && !bCenter)
690 fprintf(stderr, "Can not symmetrize without centering. Turning on -center\n");
694 axis = toupper(axtitle[0]) - 'X';
696 top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC); /* read topology file */
697 if (dens_opt[0][0] == 'n')
699 for (i = 0; (i < top->atoms.nr); i++)
701 top->atoms.atom[i].m = 1;
704 else if (dens_opt[0][0] == 'c')
706 for (i = 0; (i < top->atoms.nr); i++)
708 top->atoms.atom[i].m = top->atoms.atom[i].q;
712 snew(grpname, ngrps);
719 "\nNote: that the center of mass is calculated inside the box without applying\n"
720 "any special periodicity. If necessary, it is your responsibility to first use\n"
721 "trjconv to make sure atoms in this group are placed in the right periodicity.\n\n"
722 "Select the group to center density profiles around:\n");
723 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &ncenter,
724 &index_center, &grpname_center);
731 fprintf(stderr, "\nSelect %d group%s to calculate density for:\n", ngrps, (ngrps > 1) ? "s" : "");
732 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, ngx, index, grpname);
734 if (dens_opt[0][0] == 'e')
736 nr_electrons = get_electrons(&el_tab, ftp2fn(efDAT, NFILE, fnm));
737 fprintf(stderr, "Read %d atomtypes from datafile\n", nr_electrons);
739 calc_electron_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density,
740 &nslices, top, ePBC, axis, ngrps, &slWidth, el_tab,
741 nr_electrons, bCenter, index_center, ncenter,
746 calc_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density, &nslices, top,
747 ePBC, axis, ngrps, &slWidth, bCenter, index_center, ncenter,
751 plot_density(density, opt2fn("-o", NFILE, fnm),
752 nslices, ngrps, grpname, slWidth, dens_opt,
753 bCenter, bRelative, bSymmetrize, oenv);
755 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy"); /* view xvgr file */