src/gromacs/gmxana/gmx_density.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
  10  *
  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.
  15  *
  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.
  20  *
  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.
  25  *
  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.
  33  *
  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.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <ctype.h>
  40 #include <math.h>
  41 #include <stdlib.h>
  42 #include <string.h>
  43
  44 #include "gromacs/commandline/pargs.h"
  45 #include "gromacs/fileio/tpxio.h"
  46 #include "gromacs/fileio/trxio.h"
  47 #include "gromacs/fileio/xvgr.h"
  48 #include "gromacs/gmxana/gmx_ana.h"
  49 #include "gromacs/gmxana/gstat.h"
  50 #include "gromacs/legacyheaders/macros.h"
  51 #include "gromacs/legacyheaders/typedefs.h"
  52 #include "gromacs/legacyheaders/viewit.h"
  53 #include "gromacs/math/units.h"
  54 #include "gromacs/math/vec.h"
  55 #include "gromacs/pbcutil/pbc.h"
  56 #include "gromacs/pbcutil/rmpbc.h"
  57 #include "gromacs/topology/index.h"
  58 #include "gromacs/utility/cstringutil.h"
  59 #include "gromacs/utility/fatalerror.h"
  60 #include "gromacs/utility/futil.h"
  61 #include "gromacs/utility/smalloc.h"
  62
  63 typedef struct {
  64     char *atomname;
  65     int   nr_el;
  66 } t_electron;
  67
  68 /****************************************************************************/
  69 /* This program calculates the partial density across the box.              */
  70 /* Peter Tieleman, Mei 1995                                                 */
  71 /****************************************************************************/
  72
  73 /* used for sorting the list */
  74 int compare(void *a, void *b)
  75 {
  76     t_electron *tmp1, *tmp2;
  77     tmp1 = (t_electron *)a; tmp2 = (t_electron *)b;
  78
  79     return strcmp(tmp1->atomname, tmp2->atomname);
  80 }
  81
  82 int get_electrons(t_electron **eltab, const char *fn)
  83 {
  84     char  buffer[256];  /* to read in a line   */
  85     char  tempname[80]; /* buffer to hold name */
  86     int   tempnr;
  87
  88     FILE *in;
  89     int   nr;        /* number of atomstypes to read */
  90     int   i;
  91
  92     if (!(in = gmx_ffopen(fn, "r")))
  93     {
  94         gmx_fatal(FARGS, "Couldn't open %s. Exiting.\n", fn);
  95     }
  96
  97     if (NULL == fgets(buffer, 255, in))
  98     {
  99         gmx_fatal(FARGS, "Error reading from file %s", fn);
 100     }
 101
 102     if (sscanf(buffer, "%d", &nr) != 1)
 103     {
 104         gmx_fatal(FARGS, "Invalid number of atomtypes in datafile\n");
 105     }
 106
 107     snew(*eltab, nr);
 108
 109     for (i = 0; i < nr; i++)
 110     {
 111         if (fgets(buffer, 255, in) == NULL)
 112         {
 113             gmx_fatal(FARGS, "reading datafile. Check your datafile.\n");
 114         }
 115         if (sscanf(buffer, "%s = %d", tempname, &tempnr) != 2)
 116         {
 117             gmx_fatal(FARGS, "Invalid line in datafile at line %d\n", i+1);
 118         }
 119         (*eltab)[i].nr_el    = tempnr;
 120         (*eltab)[i].atomname = gmx_strdup(tempname);
 121     }
 122     gmx_ffclose(in);
 123
 124     /* sort the list */
 125     fprintf(stderr, "Sorting list..\n");
 126     qsort ((void*)*eltab, nr, sizeof(t_electron),
 127            (int(*)(const void*, const void*))compare);
 128
 129     return nr;
 130 }
 131
 132 void center_coords(t_atoms *atoms, atom_id *index_center, int ncenter,
 133                    matrix box, rvec x0[])
 134 {
 135     int  i, k, m;
 136     real tmass, mm;
 137     rvec com, shift, box_center;
 138
 139     tmass = 0;
 140     clear_rvec(com);
 141     for (k = 0; (k < ncenter); k++)
 142     {
 143         i = index_center[k];
 144         if (i >= atoms->nr)
 145         {
 146             gmx_fatal(FARGS, "Index %d refers to atom %d, which is larger than natoms (%d).",
 147                       k+1, i+1, atoms->nr);
 148         }
 149         mm     = atoms->atom[i].m;
 150         tmass += mm;
 151         for (m = 0; (m < DIM); m++)
 152         {
 153             com[m] += mm*x0[i][m];
 154         }
 155     }
 156     for (m = 0; (m < DIM); m++)
 157     {
 158         com[m] /= tmass;
 159     }
 160     calc_box_center(ecenterDEF, box, box_center);
 161     rvec_sub(box_center, com, shift);
 162
 163     /* Important - while the center was calculated based on a group, we should move all atoms */
 164     for (i = 0; (i < atoms->nr); i++)
 165     {
 166         rvec_dec(x0[i], shift);
 167     }
 168 }
 169
 170 void calc_electron_density(const char *fn, atom_id **index, int gnx[],
 171                            double ***slDensity, int *nslices, t_topology *top,
 172                            int ePBC,
 173                            int axis, int nr_grps, real *slWidth,
 174                            t_electron eltab[], int nr, gmx_bool bCenter,
 175                            atom_id *index_center, int ncenter,
 176                            gmx_bool bRelative, const output_env_t oenv)
 177 {
 178     rvec        *x0;            /* coordinates without pbc */
 179     matrix       box;           /* box (3x3) */
 180     double       invvol;
 181     int          natoms;        /* nr. atoms in trj */
 182     t_trxstatus *status;
 183     int          i, n,          /* loop indices */
 184                  nr_frames = 0, /* number of frames */
 185                  slice;         /* current slice */
 186     t_electron  *found;         /* found by bsearch */
 187     t_electron   sought;        /* thingie thought by bsearch */
 188     real         boxSz, aveBox;
 189     gmx_rmpbc_t  gpbc = NULL;
 190
 191     real         t,
 192                  z;
 193
 194     if (axis < 0 || axis >= DIM)
 195     {
 196         gmx_fatal(FARGS, "Invalid axes. Terminating\n");
 197     }
 198
 199     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
 200     {
 201         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 202     }
 203
 204     aveBox = 0;
 205
 206     if (!*nslices)
 207     {
 208         *nslices = (int)(box[axis][axis] * 10); /* default value */
 209         fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
 210     }
 211
 212     snew(*slDensity, nr_grps);
 213     for (i = 0; i < nr_grps; i++)
 214     {
 215         snew((*slDensity)[i], *nslices);
 216     }
 217
 218     gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
 219     /*********** Start processing trajectory ***********/
 220     do
 221     {
 222         gmx_rmpbc(gpbc, natoms, box, x0);
 223
 224         /* Translate atoms so the com of the center-group is in the
 225          * box geometrical center.
 226          */
 227         if (bCenter)
 228         {
 229             center_coords(&top->atoms, index_center, ncenter, box, x0);
 230         }
 231
 232         invvol   = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
 233
 234         if (bRelative)
 235         {
 236             *slWidth = 1.0/(*nslices);
 237             boxSz    = 1.0;
 238         }
 239         else
 240         {
 241             *slWidth = box[axis][axis]/(*nslices);
 242             boxSz    = box[axis][axis];
 243         }
 244
 245         aveBox += box[axis][axis];
 246
 247         for (n = 0; n < nr_grps; n++)
 248         {
 249             for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
 250             {
 251                 z = x0[index[n][i]][axis];
 252                 while (z < 0)
 253                 {
 254                     z += box[axis][axis];
 255                 }
 256                 while (z > box[axis][axis])
 257                 {
 258                     z -= box[axis][axis];
 259                 }
 260
 261                 if (bRelative)
 262                 {
 263                     z = z/box[axis][axis];
 264                 }
 265
 266                 /* determine which slice atom is in */
 267                 if (bCenter)
 268                 {
 269                     slice = floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2;
 270                 }
 271                 else
 272                 {
 273                     slice = (z / (*slWidth));
 274                 }
 275                 sought.nr_el    = 0;
 276                 sought.atomname = gmx_strdup(*(top->atoms.atomname[index[n][i]]));
 277
 278                 /* now find the number of electrons. This is not efficient. */
 279                 found = (t_electron *)
 280                     bsearch((const void *)&sought,
 281                             (const void *)eltab, nr, sizeof(t_electron),
 282                             (int(*)(const void*, const void*))compare);
 283
 284                 if (found == NULL)
 285                 {
 286                     fprintf(stderr, "Couldn't find %s. Add it to the .dat file\n",
 287                             *(top->atoms.atomname[index[n][i]]));
 288                 }
 289                 else
 290                 {
 291                     (*slDensity)[n][slice] += (found->nr_el -
 292                                                top->atoms.atom[index[n][i]].q)*invvol;
 293                 }
 294                 free(sought.atomname);
 295             }
 296         }
 297         nr_frames++;
 298     }
 299     while (read_next_x(oenv, status, &t, x0, box));
 300     gmx_rmpbc_done(gpbc);
 301
 302     /*********** done with status file **********/
 303     close_trj(status);
 304
 305 /* slDensity now contains the total number of electrons per slice, summed
 306    over all frames. Now divide by nr_frames and volume of slice
 307  */
 308
 309     fprintf(stderr, "\nRead %d frames from trajectory. Counting electrons\n",
 310             nr_frames);
 311
 312     if (bRelative)
 313     {
 314         aveBox  /= nr_frames;
 315         *slWidth = aveBox/(*nslices);
 316     }
 317
 318     for (n = 0; n < nr_grps; n++)
 319     {
 320         for (i = 0; i < *nslices; i++)
 321         {
 322             (*slDensity)[n][i] /= nr_frames;
 323         }
 324     }
 325
 326     sfree(x0); /* free memory used by coordinate array */
 327 }
 328
 329 void calc_density(const char *fn, atom_id **index, int gnx[],
 330                   double ***slDensity, int *nslices, t_topology *top, int ePBC,
 331                   int axis, int nr_grps, real *slWidth, gmx_bool bCenter,
 332                   atom_id *index_center, int ncenter,
 333                   gmx_bool bRelative, const output_env_t oenv)
 334 {
 335     rvec        *x0;            /* coordinates without pbc */
 336     matrix       box;           /* box (3x3) */
 337     double       invvol;
 338     int          natoms;        /* nr. atoms in trj */
 339     t_trxstatus *status;
 340     int        **slCount,       /* nr. of atoms in one slice for a group */
 341                  i, j, n,       /* loop indices */
 342                  ax1       = 0, ax2 = 0,
 343                  nr_frames = 0, /* number of frames */
 344                  slice;         /* current slice */
 345     real         t,
 346                  z;
 347     real         boxSz, aveBox;
 348     char        *buf;    /* for tmp. keeping atomname */
 349     gmx_rmpbc_t  gpbc = NULL;
 350
 351     if (axis < 0 || axis >= DIM)
 352     {
 353         gmx_fatal(FARGS, "Invalid axes. Terminating\n");
 354     }
 355
 356     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
 357     {
 358         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 359     }
 360
 361     aveBox = 0;
 362
 363     if (!*nslices)
 364     {
 365         *nslices = (int)(box[axis][axis] * 10); /* default value */
 366         fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
 367     }
 368
 369     snew(*slDensity, nr_grps);
 370     for (i = 0; i < nr_grps; i++)
 371     {
 372         snew((*slDensity)[i], *nslices);
 373     }
 374
 375     gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
 376     /*********** Start processing trajectory ***********/
 377     do
 378     {
 379         gmx_rmpbc(gpbc, natoms, box, x0);
 380
 381         /* Translate atoms so the com of the center-group is in the
 382          * box geometrical center.
 383          */
 384         if (bCenter)
 385         {
 386             center_coords(&top->atoms, index_center, ncenter, box, x0);
 387         }
 388
 389         invvol   = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
 390
 391         if (bRelative)
 392         {
 393             *slWidth = 1.0/(*nslices);
 394             boxSz    = 1.0;
 395         }
 396         else
 397         {
 398             *slWidth = box[axis][axis]/(*nslices);
 399             boxSz    = box[axis][axis];
 400         }
 401
 402         aveBox += box[axis][axis];
 403
 404         for (n = 0; n < nr_grps; n++)
 405         {
 406             for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
 407             {
 408                 z = x0[index[n][i]][axis];
 409                 while (z < 0)
 410                 {
 411                     z += box[axis][axis];
 412                 }
 413                 while (z > box[axis][axis])
 414                 {
 415                     z -= box[axis][axis];
 416                 }
 417
 418                 if (bRelative)
 419                 {
 420                     z = z/box[axis][axis];
 421                 }
 422
 423                 /* determine which slice atom is in */
 424                 if (bCenter)
 425                 {
 426                     slice = floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2;
 427                 }
 428                 else
 429                 {
 430                     slice = floor(z / (*slWidth));
 431                 }
 432
 433                 /* Slice should already be 0<=slice<nslices, but we just make
 434                  * sure we are not hit by IEEE rounding errors since we do
 435                  * math operations after applying PBC above.
 436                  */
 437                 if (slice < 0)
 438                 {
 439                     slice += *nslices;
 440                 }
 441                 else if (slice >= *nslices)
 442                 {
 443                     slice -= *nslices;
 444                 }
 445
 446                 (*slDensity)[n][slice] += top->atoms.atom[index[n][i]].m*invvol;
 447             }
 448         }
 449         nr_frames++;
 450     }
 451     while (read_next_x(oenv, status, &t, x0, box));
 452     gmx_rmpbc_done(gpbc);
 453
 454     /*********** done with status file **********/
 455     close_trj(status);
 456
 457     /* slDensity now contains the total mass per slice, summed over all
 458        frames. Now divide by nr_frames and volume of slice
 459      */
 460
 461     fprintf(stderr, "\nRead %d frames from trajectory. Calculating density\n",
 462             nr_frames);
 463
 464     if (bRelative)
 465     {
 466         aveBox  /= nr_frames;
 467         *slWidth = aveBox/(*nslices);
 468     }
 469
 470     for (n = 0; n < nr_grps; n++)
 471     {
 472         for (i = 0; i < *nslices; i++)
 473         {
 474             (*slDensity)[n][i] /= nr_frames;
 475         }
 476     }
 477
 478     sfree(x0); /* free memory used by coordinate array */
 479 }
 480
 481 void plot_density(double *slDensity[], const char *afile, int nslices,
 482                   int nr_grps, char *grpname[], real slWidth,
 483                   const char **dens_opt,
 484                   gmx_bool bCenter, gmx_bool bRelative, gmx_bool bSymmetrize,
 485                   const output_env_t oenv)
 486 {
 487     FILE       *den;
 488     const char *title  = NULL;
 489     const char *xlabel = NULL;
 490     const char *ylabel = NULL;
 491     int         slice, n;
 492     real        ddd;
 493     real        axispos;
 494
 495     title = bSymmetrize ? "Symmetrized partial density" : "Partial density";
 496
 497     if (bCenter)
 498     {
 499         xlabel = bRelative ?
 500             "Average relative position from center (nm)" :
 501             "Relative position from center (nm)";
 502     }
 503     else
 504     {
 505         xlabel = bRelative ? "Average coordinate (nm)" : "Coordinate (nm)";
 506     }
 507
 508     switch (dens_opt[0][0])
 509     {
 510         case 'm': ylabel = "Density (kg m\\S-3\\N)"; break;
 511         case 'n': ylabel = "Number density (nm\\S-3\\N)"; break;
 512         case 'c': ylabel = "Charge density (e nm\\S-3\\N)"; break;
 513         case 'e': ylabel = "Electron density (e nm\\S-3\\N)"; break;
 514     }
 515
 516     den = xvgropen(afile,
 517                    title, xlabel, ylabel, oenv);
 518
 519     xvgr_legend(den, nr_grps, (const char**)grpname, oenv);
 520
 521     for (slice = 0; (slice < nslices); slice++)
 522     {
 523         if (bCenter)
 524         {
 525             axispos = (slice - nslices/2.0 + 0.5) * slWidth;
 526         }
 527         else
 528         {
 529             axispos = (slice + 0.5) * slWidth;
 530         }
 531         fprintf(den, "%12g  ", axispos);
 532         for (n = 0; (n < nr_grps); n++)
 533         {
 534             if (bSymmetrize)
 535             {
 536                 ddd = (slDensity[n][slice]+slDensity[n][nslices-slice-1])*0.5;
 537             }
 538             else
 539             {
 540                 ddd = slDensity[n][slice];
 541             }
 542             if (dens_opt[0][0] == 'm')
 543             {
 544                 fprintf(den, "   %12g", ddd*AMU/(NANO*NANO*NANO));
 545             }
 546             else
 547             {
 548                 fprintf(den, "   %12g", ddd);
 549             }
 550         }
 551         fprintf(den, "\n");
 552     }
 553
 554     gmx_ffclose(den);
 555 }
 556
 557 int gmx_density(int argc, char *argv[])
 558 {
 559     const char        *desc[] = {
 560         "[THISMODULE] computes partial densities across the box, using an index file.[PAR]",
 561         "For the total density of NPT simulations, use [gmx-energy] instead.",
 562         "[PAR]",
 563
 564         "Option [TT]-center[tt] performs the histogram binning relative to the center",
 565         "of an arbitrary group, in absolute box coordinates. If you are calculating",
 566         "profiles along the Z axis box dimension bZ, output would be from -bZ/2 to",
 567         "bZ/2 if you center based on the entire system.",
 568         "Note that this behaviour has changed in Gromacs 5.0; earlier versions",
 569         "merely performed a static binning in (0,bZ) and shifted the output. Now",
 570         "we compute the center for each frame and bin in (-bZ/2,bZ/2).[PAR]",
 571
 572         "Option [TT]-symm[tt] symmetrizes the output around the center. This will",
 573         "automatically turn on [TT]-center[tt] too.",
 574
 575         "Option [TT]-relative[tt] performs the binning in relative instead of absolute",
 576         "box coordinates, and scales the final output with the average box dimension",
 577         "along the output axis. This can be used in combination with [TT]-center[tt].[PAR]",
 578
 579         "Densities are in kg/m^3, and number densities or electron densities can also be",
 580         "calculated. For electron densities, a file describing the number of",
 581         "electrons for each type of atom should be provided using [TT]-ei[tt].",
 582         "It should look like:[BR]",
 583         "   [TT]2[tt][BR]",
 584         "   [TT]atomname = nrelectrons[tt][BR]",
 585         "   [TT]atomname = nrelectrons[tt][BR]",
 586         "The first line contains the number of lines to read from the file.",
 587         "There should be one line for each unique atom name in your system.",
 588         "The number of electrons for each atom is modified by its atomic",
 589         "partial charge.[PAR]",
 590
 591         "IMPORTANT CONSIDERATIONS FOR BILAYERS[PAR]",
 592         "One of the most common usage scenarios is to calculate the density of various",
 593         "groups across a lipid bilayer, typically with the z axis being the normal",
 594         "direction. For short simulations, small systems, and fixed box sizes this",
 595         "will work fine, but for the more general case lipid bilayers can be complicated.",
 596         "The first problem that while both proteins and lipids have low volume",
 597         "compressibility, lipids have quite high area compressiblity. This means the",
 598         "shape of the box (thickness and area/lipid) will fluctuate substantially even",
 599         "for a fully relaxed system. Since Gromacs places the box between the origin",
 600         "and positive coordinates, this in turn means that a bilayer centered in the",
 601         "box will move a bit up/down due to these fluctuations, and smear out your",
 602         "profile. The easiest way to fix this (if you want pressure coupling) is",
 603         "to use the [TT]-center[tt] option that calculates the density profile with",
 604         "respect to the center of the box. Note that you can still center on the",
 605         "bilayer part even if you have a complex non-symmetric system with a bilayer",
 606         "and, say, membrane proteins - then our output will simply have more values",
 607         "on one side of the (center) origin reference.[PAR]",
 608
 609         "Even the centered calculation will lead to some smearing out the output",
 610         "profiles, as lipids themselves are compressed and expanded. In most cases",
 611         "you probably want this (since it corresponds to macroscopic experiments),",
 612         "but if you want to look at molecular details you can use the [TT]-relative[tt]",
 613         "option to attempt to remove even more of the effects of volume fluctuations.[PAR]",
 614
 615         "Finally, large bilayers that are not subject to a surface tension will exhibit",
 616         "undulatory fluctuations, where there are 'waves' forming in the system.",
 617         "This is a fundamental property of the biological system, and if you are",
 618         "comparing against experiments you likely want to include the undulation",
 619         "smearing effect.",
 620         "",
 621     };
 622
 623     output_env_t       oenv;
 624     static const char *dens_opt[] =
 625     { NULL, "mass", "number", "charge", "electron", NULL };
 626     static int         axis        = 2;  /* normal to memb. default z  */
 627     static const char *axtitle     = "Z";
 628     static int         nslices     = 50; /* nr of slices defined       */
 629     static int         ngrps       = 1;  /* nr. of groups              */
 630     static gmx_bool    bSymmetrize = FALSE;
 631     static gmx_bool    bCenter     = FALSE;
 632     static gmx_bool    bRelative   = FALSE;
 633
 634     t_pargs            pa[]        = {
 635         { "-d", FALSE, etSTR, {&axtitle},
 636           "Take the normal on the membrane in direction X, Y or Z." },
 637         { "-sl",  FALSE, etINT, {&nslices},
 638           "Divide the box in this number of slices." },
 639         { "-dens",    FALSE, etENUM, {dens_opt},
 640           "Density"},
 641         { "-ng",       FALSE, etINT, {&ngrps},
 642           "Number of groups of which to compute densities." },
 643         { "-center",   FALSE, etBOOL, {&bCenter},
 644           "Perform the binning relative to the center of the (changing) box. Useful for bilayers." },
 645         { "-symm",     FALSE, etBOOL, {&bSymmetrize},
 646           "Symmetrize the density along the axis, with respect to the center. Useful for bilayers." },
 647         { "-relative", FALSE, etBOOL, {&bRelative},
 648           "Use relative coordinates for changing boxes and scale output by average dimensions." }
 649     };
 650
 651     const char        *bugs[] = {
 652         "When calculating electron densities, atomnames are used instead of types. This is bad.",
 653     };
 654
 655     double           **density;        /* density per slice          */
 656     real               slWidth;        /* width of one slice         */
 657     char              *grpname_center; /* centering group name     */
 658     char             **grpname;        /* groupnames                 */
 659     int                nr_electrons;   /* nr. electrons              */
 660     int                ncenter;        /* size of centering group    */
 661     int               *ngx;            /* sizes of groups            */
 662     t_electron        *el_tab;         /* tabel with nr. of electrons*/
 663     t_topology        *top;            /* topology               */
 664     int                ePBC;
 665     atom_id           *index_center;   /* index for centering group  */
 666     atom_id          **index;          /* indices for all groups     */
 667     int                i;
 668
 669     t_filenm           fnm[] = { /* files for g_density       */
 670         { efTRX, "-f", NULL,  ffREAD },
 671         { efNDX, NULL, NULL,  ffOPTRD },
 672         { efTPR, NULL, NULL,  ffREAD },
 673         { efDAT, "-ei", "electrons", ffOPTRD }, /* file with nr. of electrons */
 674         { efXVG, "-o", "density", ffWRITE },
 675     };
 676
 677 #define NFILE asize(fnm)
 678
 679     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
 680                            NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs,
 681                            &oenv))
 682     {
 683         return 0;
 684     }
 685
 686     if (bSymmetrize && !bCenter)
 687     {
 688         fprintf(stderr, "Can not symmetrize without centering. Turning on -center\n");
 689         bCenter = TRUE;
 690     }
 691     /* Calculate axis */
 692     axis = toupper(axtitle[0]) - 'X';
 693
 694     top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC); /* read topology file */
 695     if (dens_opt[0][0] == 'n')
 696     {
 697         for (i = 0; (i < top->atoms.nr); i++)
 698         {
 699             top->atoms.atom[i].m = 1;
 700         }
 701     }
 702     else if (dens_opt[0][0] == 'c')
 703     {
 704         for (i = 0; (i < top->atoms.nr); i++)
 705         {
 706             top->atoms.atom[i].m = top->atoms.atom[i].q;
 707         }
 708     }
 709
 710     snew(grpname, ngrps);
 711     snew(index, ngrps);
 712     snew(ngx, ngrps);
 713
 714     if (bCenter)
 715     {
 716         fprintf(stderr,
 717                 "\nNote: that the center of mass is calculated inside the box without applying\n"
 718                 "any special periodicity. If necessary, it is your responsibility to first use\n"
 719                 "trjconv to make sure atoms in this group are placed in the right periodicity.\n\n"
 720                 "Select the group to center density profiles around:\n");
 721         get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &ncenter,
 722                   &index_center, &grpname_center);
 723     }
 724     else
 725     {
 726         ncenter = 0;
 727     }
 728
 729     fprintf(stderr, "\nSelect %d group%s to calculate density for:\n", ngrps, (ngrps > 1) ? "s" : "");
 730     get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, ngx, index, grpname);
 731
 732     if (dens_opt[0][0] == 'e')
 733     {
 734         nr_electrons =  get_electrons(&el_tab, ftp2fn(efDAT, NFILE, fnm));
 735         fprintf(stderr, "Read %d atomtypes from datafile\n", nr_electrons);
 736
 737         calc_electron_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density,
 738                               &nslices, top, ePBC, axis, ngrps, &slWidth, el_tab,
 739                               nr_electrons, bCenter, index_center, ncenter,
 740                               bRelative, oenv);
 741     }
 742     else
 743     {
 744         calc_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density, &nslices, top,
 745                      ePBC, axis, ngrps, &slWidth, bCenter, index_center, ncenter,
 746                      bRelative, oenv);
 747     }
 748
 749     plot_density(density, opt2fn("-o", NFILE, fnm),
 750                  nslices, ngrps, grpname, slWidth, dens_opt,
 751                  bCenter, bRelative, bSymmetrize, oenv);
 752
 753     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");  /* view xvgr file */
 754     return 0;
 755 }