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/utility/cstringutil.h"
  45 #include "gromacs/legacyheaders/typedefs.h"
  46 #include "gromacs/legacyheaders/macros.h"
  47 #include "gstat.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"
  54 #include "gmx_ana.h"
  55 #include "gromacs/legacyheaders/macros.h"
  56
  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"
  64
  65 typedef struct {
  66     char *atomname;
  67     int   nr_el;
  68 } t_electron;
  69
  70 /****************************************************************************/
  71 /* This program calculates the partial density across the box.              */
  72 /* Peter Tieleman, Mei 1995                                                 */
  73 /****************************************************************************/
  74
  75 /* used for sorting the list */
  76 int compare(void *a, void *b)
  77 {
  78     t_electron *tmp1, *tmp2;
  79     tmp1 = (t_electron *)a; tmp2 = (t_electron *)b;
  80
  81     return strcmp(tmp1->atomname, tmp2->atomname);
  82 }
  83
  84 int get_electrons(t_electron **eltab, const char *fn)
  85 {
  86     char  buffer[256];  /* to read in a line   */
  87     char  tempname[80]; /* buffer to hold name */
  88     int   tempnr;
  89
  90     FILE *in;
  91     int   nr;        /* number of atomstypes to read */
  92     int   i;
  93
  94     if (!(in = gmx_ffopen(fn, "r")))
  95     {
  96         gmx_fatal(FARGS, "Couldn't open %s. Exiting.\n", fn);
  97     }
  98
  99     if (NULL == fgets(buffer, 255, in))
 100     {
 101         gmx_fatal(FARGS, "Error reading from file %s", fn);
 102     }
 103
 104     if (sscanf(buffer, "%d", &nr) != 1)
 105     {
 106         gmx_fatal(FARGS, "Invalid number of atomtypes in datafile\n");
 107     }
 108
 109     snew(*eltab, nr);
 110
 111     for (i = 0; i < nr; i++)
 112     {
 113         if (fgets(buffer, 255, in) == NULL)
 114         {
 115             gmx_fatal(FARGS, "reading datafile. Check your datafile.\n");
 116         }
 117         if (sscanf(buffer, "%s = %d", tempname, &tempnr) != 2)
 118         {
 119             gmx_fatal(FARGS, "Invalid line in datafile at line %d\n", i+1);
 120         }
 121         (*eltab)[i].nr_el    = tempnr;
 122         (*eltab)[i].atomname = gmx_strdup(tempname);
 123     }
 124     gmx_ffclose(in);
 125
 126     /* sort the list */
 127     fprintf(stderr, "Sorting list..\n");
 128     qsort ((void*)*eltab, nr, sizeof(t_electron),
 129            (int(*)(const void*, const void*))compare);
 130
 131     return nr;
 132 }
 133
 134 void center_coords(t_atoms *atoms, atom_id *index_center, int ncenter,
 135                    matrix box, rvec x0[])
 136 {
 137     int  i, k, m;
 138     real tmass, mm;
 139     rvec com, shift, box_center;
 140
 141     tmass = 0;
 142     clear_rvec(com);
 143     for (k = 0; (k < ncenter); k++)
 144     {
 145         i = index_center[k];
 146         if (i >= atoms->nr)
 147         {
 148             gmx_fatal(FARGS, "Index %d refers to atom %d, which is larger than natoms (%d).",
 149                       k+1, i+1, atoms->nr);
 150         }
 151         mm     = atoms->atom[i].m;
 152         tmass += mm;
 153         for (m = 0; (m < DIM); m++)
 154         {
 155             com[m] += mm*x0[i][m];
 156         }
 157     }
 158     for (m = 0; (m < DIM); m++)
 159     {
 160         com[m] /= tmass;
 161     }
 162     calc_box_center(ecenterDEF, box, box_center);
 163     rvec_sub(box_center, com, shift);
 164
 165     /* Important - while the center was calculated based on a group, we should move all atoms */
 166     for (i = 0; (i < atoms->nr); i++)
 167     {
 168         rvec_dec(x0[i], shift);
 169     }
 170 }
 171
 172 void calc_electron_density(const char *fn, atom_id **index, int gnx[],
 173                            double ***slDensity, int *nslices, t_topology *top,
 174                            int ePBC,
 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)
 179 {
 180     rvec        *x0;            /* coordinates without pbc */
 181     matrix       box;           /* box (3x3) */
 182     double       invvol;
 183     int          natoms;        /* nr. atoms in trj */
 184     t_trxstatus *status;
 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 */
 190     real         boxSz, aveBox;
 191     gmx_rmpbc_t  gpbc = NULL;
 192
 193     real         t,
 194                  z;
 195
 196     if (axis < 0 || axis >= DIM)
 197     {
 198         gmx_fatal(FARGS, "Invalid axes. Terminating\n");
 199     }
 200
 201     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
 202     {
 203         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 204     }
 205
 206     aveBox = 0;
 207
 208     if (!*nslices)
 209     {
 210         *nslices = (int)(box[axis][axis] * 10); /* default value */
 211         fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
 212     }
 213
 214     snew(*slDensity, nr_grps);
 215     for (i = 0; i < nr_grps; i++)
 216     {
 217         snew((*slDensity)[i], *nslices);
 218     }
 219
 220     gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
 221     /*********** Start processing trajectory ***********/
 222     do
 223     {
 224         gmx_rmpbc(gpbc, natoms, box, x0);
 225
 226         /* Translate atoms so the com of the center-group is in the
 227          * box geometrical center.
 228          */
 229         if (bCenter)
 230         {
 231             center_coords(&top->atoms, index_center, ncenter, box, x0);
 232         }
 233
 234         invvol   = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
 235
 236         if (bRelative)
 237         {
 238             *slWidth = 1.0/(*nslices);
 239             boxSz    = 1.0;
 240         }
 241         else
 242         {
 243             *slWidth = box[axis][axis]/(*nslices);
 244             boxSz    = box[axis][axis];
 245         }
 246
 247         aveBox += box[axis][axis];
 248
 249         for (n = 0; n < nr_grps; n++)
 250         {
 251             for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
 252             {
 253                 z = x0[index[n][i]][axis];
 254                 while (z < 0)
 255                 {
 256                     z += box[axis][axis];
 257                 }
 258                 while (z > box[axis][axis])
 259                 {
 260                     z -= box[axis][axis];
 261                 }
 262
 263                 if (bRelative)
 264                 {
 265                     z = z/box[axis][axis];
 266                 }
 267
 268                 /* determine which slice atom is in */
 269                 if (bCenter)
 270                 {
 271                     slice = floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2;
 272                 }
 273                 else
 274                 {
 275                     slice = (z / (*slWidth));
 276                 }
 277                 sought.nr_el    = 0;
 278                 sought.atomname = gmx_strdup(*(top->atoms.atomname[index[n][i]]));
 279
 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);
 285
 286                 if (found == NULL)
 287                 {
 288                     fprintf(stderr, "Couldn't find %s. Add it to the .dat file\n",
 289                             *(top->atoms.atomname[index[n][i]]));
 290                 }
 291                 else
 292                 {
 293                     (*slDensity)[n][slice] += (found->nr_el -
 294                                                top->atoms.atom[index[n][i]].q)*invvol;
 295                 }
 296                 free(sought.atomname);
 297             }
 298         }
 299         nr_frames++;
 300     }
 301     while (read_next_x(oenv, status, &t, x0, box));
 302     gmx_rmpbc_done(gpbc);
 303
 304     /*********** done with status file **********/
 305     close_trj(status);
 306
 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
 309  */
 310
 311     fprintf(stderr, "\nRead %d frames from trajectory. Counting electrons\n",
 312             nr_frames);
 313
 314     if (bRelative)
 315     {
 316         aveBox  /= nr_frames;
 317         *slWidth = aveBox/(*nslices);
 318     }
 319
 320     for (n = 0; n < nr_grps; n++)
 321     {
 322         for (i = 0; i < *nslices; i++)
 323         {
 324             (*slDensity)[n][i] /= nr_frames;
 325         }
 326     }
 327
 328     sfree(x0); /* free memory used by coordinate array */
 329 }
 330
 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)
 336 {
 337     rvec        *x0;            /* coordinates without pbc */
 338     matrix       box;           /* box (3x3) */
 339     double       invvol;
 340     int          natoms;        /* nr. atoms in trj */
 341     t_trxstatus *status;
 342     int        **slCount,       /* nr. of atoms in one slice for a group */
 343                  i, j, n,       /* loop indices */
 344                  ax1       = 0, ax2 = 0,
 345                  nr_frames = 0, /* number of frames */
 346                  slice;         /* current slice */
 347     real         t,
 348                  z;
 349     real         boxSz, aveBox;
 350     char        *buf;    /* for tmp. keeping atomname */
 351     gmx_rmpbc_t  gpbc = NULL;
 352
 353     if (axis < 0 || axis >= DIM)
 354     {
 355         gmx_fatal(FARGS, "Invalid axes. Terminating\n");
 356     }
 357
 358     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
 359     {
 360         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 361     }
 362
 363     aveBox = 0;
 364
 365     if (!*nslices)
 366     {
 367         *nslices = (int)(box[axis][axis] * 10); /* default value */
 368         fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
 369     }
 370
 371     snew(*slDensity, nr_grps);
 372     for (i = 0; i < nr_grps; i++)
 373     {
 374         snew((*slDensity)[i], *nslices);
 375     }
 376
 377     gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
 378     /*********** Start processing trajectory ***********/
 379     do
 380     {
 381         gmx_rmpbc(gpbc, natoms, box, x0);
 382
 383         /* Translate atoms so the com of the center-group is in the
 384          * box geometrical center.
 385          */
 386         if (bCenter)
 387         {
 388             center_coords(&top->atoms, index_center, ncenter, box, x0);
 389         }
 390
 391         invvol   = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
 392
 393         if (bRelative)
 394         {
 395             *slWidth = 1.0/(*nslices);
 396             boxSz    = 1.0;
 397         }
 398         else
 399         {
 400             *slWidth = box[axis][axis]/(*nslices);
 401             boxSz    = box[axis][axis];
 402         }
 403
 404         aveBox += box[axis][axis];
 405
 406         for (n = 0; n < nr_grps; n++)
 407         {
 408             for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
 409             {
 410                 z = x0[index[n][i]][axis];
 411                 while (z < 0)
 412                 {
 413                     z += box[axis][axis];
 414                 }
 415                 while (z > box[axis][axis])
 416                 {
 417                     z -= box[axis][axis];
 418                 }
 419
 420                 if (bRelative)
 421                 {
 422                     z = z/box[axis][axis];
 423                 }
 424
 425                 /* determine which slice atom is in */
 426                 if (bCenter)
 427                 {
 428                     slice = floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2;
 429                 }
 430                 else
 431                 {
 432                     slice = floor(z / (*slWidth));
 433                 }
 434
 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.
 438                  */
 439                 if (slice < 0)
 440                 {
 441                     slice += *nslices;
 442                 }
 443                 else if (slice >= *nslices)
 444                 {
 445                     slice -= *nslices;
 446                 }
 447
 448                 (*slDensity)[n][slice] += top->atoms.atom[index[n][i]].m*invvol;
 449             }
 450         }
 451         nr_frames++;
 452     }
 453     while (read_next_x(oenv, status, &t, x0, box));
 454     gmx_rmpbc_done(gpbc);
 455
 456     /*********** done with status file **********/
 457     close_trj(status);
 458
 459     /* slDensity now contains the total mass per slice, summed over all
 460        frames. Now divide by nr_frames and volume of slice
 461      */
 462
 463     fprintf(stderr, "\nRead %d frames from trajectory. Calculating density\n",
 464             nr_frames);
 465
 466     if (bRelative)
 467     {
 468         aveBox  /= nr_frames;
 469         *slWidth = aveBox/(*nslices);
 470     }
 471
 472     for (n = 0; n < nr_grps; n++)
 473     {
 474         for (i = 0; i < *nslices; i++)
 475         {
 476             (*slDensity)[n][i] /= nr_frames;
 477         }
 478     }
 479
 480     sfree(x0); /* free memory used by coordinate array */
 481 }
 482
 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)
 488 {
 489     FILE       *den;
 490     const char *title  = NULL;
 491     const char *xlabel = NULL;
 492     const char *ylabel = NULL;
 493     int         slice, n;
 494     real        ddd;
 495     real        axispos;
 496
 497     title = bSymmetrize ? "Symmetrized partial density" : "Partial density";
 498
 499     if (bCenter)
 500     {
 501         xlabel = bRelative ?
 502             "Average relative position from center (nm)" :
 503             "Relative position from center (nm)";
 504     }
 505     else
 506     {
 507         xlabel = bRelative ? "Average coordinate (nm)" : "Coordinate (nm)";
 508     }
 509
 510     switch (dens_opt[0][0])
 511     {
 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;
 516     }
 517
 518     den = xvgropen(afile,
 519                    title, xlabel, ylabel, oenv);
 520
 521     xvgr_legend(den, nr_grps, (const char**)grpname, oenv);
 522
 523     for (slice = 0; (slice < nslices); slice++)
 524     {
 525         if (bCenter)
 526         {
 527             axispos = (slice - nslices/2.0 + 0.5) * slWidth;
 528         }
 529         else
 530         {
 531             axispos = (slice + 0.5) * slWidth;
 532         }
 533         fprintf(den, "%12g  ", axispos);
 534         for (n = 0; (n < nr_grps); n++)
 535         {
 536             if (bSymmetrize)
 537             {
 538                 ddd = (slDensity[n][slice]+slDensity[n][nslices-slice-1])*0.5;
 539             }
 540             else
 541             {
 542                 ddd = slDensity[n][slice];
 543             }
 544             if (dens_opt[0][0] == 'm')
 545             {
 546                 fprintf(den, "   %12g", ddd*AMU/(NANO*NANO*NANO));
 547             }
 548             else
 549             {
 550                 fprintf(den, "   %12g", ddd);
 551             }
 552         }
 553         fprintf(den, "\n");
 554     }
 555
 556     gmx_ffclose(den);
 557 }
 558
 559 int gmx_density(int argc, char *argv[])
 560 {
 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.",
 564         "[PAR]",
 565
 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]",
 573
 574         "Option [TT]-symm[tt] symmetrizes the output around the center. This will",
 575         "automatically turn on [TT]-center[tt] too.",
 576
 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]",
 580
 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]",
 585         "   [TT]2[tt][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]",
 592
 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]",
 610
 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]",
 616
 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",
 621         "smearing effect.",
 622         "",
 623     };
 624
 625     output_env_t       oenv;
 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;
 635
 636     t_pargs            pa[]        = {
 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},
 642           "Density"},
 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." }
 651     };
 652
 653     const char        *bugs[] = {
 654         "When calculating electron densities, atomnames are used instead of types. This is bad.",
 655     };
 656
 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               */
 666     int                ePBC;
 667     atom_id           *index_center;   /* index for centering group  */
 668     atom_id          **index;          /* indices for all groups     */
 669     int                i;
 670
 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 },
 677     };
 678
 679 #define NFILE asize(fnm)
 680
 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,
 683                            &oenv))
 684     {
 685         return 0;
 686     }
 687
 688     if (bSymmetrize && !bCenter)
 689     {
 690         fprintf(stderr, "Can not symmetrize without centering. Turning on -center\n");
 691         bCenter = TRUE;
 692     }
 693     /* Calculate axis */
 694     axis = toupper(axtitle[0]) - 'X';
 695
 696     top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC); /* read topology file */
 697     if (dens_opt[0][0] == 'n')
 698     {
 699         for (i = 0; (i < top->atoms.nr); i++)
 700         {
 701             top->atoms.atom[i].m = 1;
 702         }
 703     }
 704     else if (dens_opt[0][0] == 'c')
 705     {
 706         for (i = 0; (i < top->atoms.nr); i++)
 707         {
 708             top->atoms.atom[i].m = top->atoms.atom[i].q;
 709         }
 710     }
 711
 712     snew(grpname, ngrps);
 713     snew(index, ngrps);
 714     snew(ngx, ngrps);
 715
 716     if (bCenter)
 717     {
 718         fprintf(stderr,
 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);
 725     }
 726     else
 727     {
 728         ncenter = 0;
 729     }
 730
 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);
 733
 734     if (dens_opt[0][0] == 'e')
 735     {
 736         nr_electrons =  get_electrons(&el_tab, ftp2fn(efDAT, NFILE, fnm));
 737         fprintf(stderr, "Read %d atomtypes from datafile\n", nr_electrons);
 738
 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,
 742                               bRelative, oenv);
 743     }
 744     else
 745     {
 746         calc_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density, &nslices, top,
 747                      ePBC, axis, ngrps, &slWidth, bCenter, index_center, ncenter,
 748                      bRelative, oenv);
 749     }
 750
 751     plot_density(density, opt2fn("-o", NFILE, fnm),
 752                  nslices, ngrps, grpname, slWidth, dens_opt,
 753                  bCenter, bRelative, bSymmetrize, oenv);
 754
 755     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");  /* view xvgr file */
 756     return 0;
 757 }