src/gromacs/gmxana/gmx_density.cpp

   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,2015,2017, 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 <cctype>
  40 #include <cmath>
  41 #include <cstdlib>
  42 #include <cstring>
  43
  44 #include "gromacs/commandline/pargs.h"
  45 #include "gromacs/commandline/viewit.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/math/units.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/pbcutil/pbc.h"
  53 #include "gromacs/pbcutil/rmpbc.h"
  54 #include "gromacs/topology/index.h"
  55 #include "gromacs/topology/topology.h"
  56 #include "gromacs/utility/arraysize.h"
  57 #include "gromacs/utility/cstringutil.h"
  58 #include "gromacs/utility/fatalerror.h"
  59 #include "gromacs/utility/futil.h"
  60 #include "gromacs/utility/gmxassert.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 static int compare(void *a, void *b)
  75 {
  76     t_electron *tmp1, *tmp2;
  77     tmp1 = (t_electron *)a; tmp2 = (t_electron *)b;
  78
  79     return std::strcmp(tmp1->atomname, tmp2->atomname);
  80 }
  81
  82 static 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 (nullptr == 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) == nullptr)
 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 static void center_coords(t_atoms *atoms, int *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(com, box_center, 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 static void calc_electron_density(const char *fn, int **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                                   int *index_center, int ncenter,
 176                                   gmx_bool bRelative, const gmx_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 = nullptr;
 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 = static_cast<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 = static_cast<int>(std::floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2);
 270                 }
 271                 else
 272                 {
 273                     slice = static_cast<int>(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 == nullptr)
 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_trx(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 static void calc_density(const char *fn, int **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                          int *index_center, int ncenter,
 333                          gmx_bool bRelative, const gmx_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          i, n,          /* loop indices */
 341                  nr_frames = 0, /* number of frames */
 342                  slice;         /* current slice */
 343     real         t,
 344                  z;
 345     real         boxSz, aveBox;
 346     gmx_rmpbc_t  gpbc = nullptr;
 347
 348     if (axis < 0 || axis >= DIM)
 349     {
 350         gmx_fatal(FARGS, "Invalid axes. Terminating\n");
 351     }
 352
 353     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
 354     {
 355         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 356     }
 357
 358     aveBox = 0;
 359
 360     if (!*nslices)
 361     {
 362         *nslices = static_cast<int>(box[axis][axis] * 10); /* default value */
 363         fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
 364     }
 365
 366     snew(*slDensity, nr_grps);
 367     for (i = 0; i < nr_grps; i++)
 368     {
 369         snew((*slDensity)[i], *nslices);
 370     }
 371
 372     gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
 373     /*********** Start processing trajectory ***********/
 374     do
 375     {
 376         gmx_rmpbc(gpbc, natoms, box, x0);
 377
 378         /* Translate atoms so the com of the center-group is in the
 379          * box geometrical center.
 380          */
 381         if (bCenter)
 382         {
 383             center_coords(&top->atoms, index_center, ncenter, box, x0);
 384         }
 385
 386         invvol   = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
 387
 388         if (bRelative)
 389         {
 390             *slWidth = 1.0/(*nslices);
 391             boxSz    = 1.0;
 392         }
 393         else
 394         {
 395             *slWidth = box[axis][axis]/(*nslices);
 396             boxSz    = box[axis][axis];
 397         }
 398
 399         aveBox += box[axis][axis];
 400
 401         for (n = 0; n < nr_grps; n++)
 402         {
 403             for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
 404             {
 405                 z = x0[index[n][i]][axis];
 406                 while (z < 0)
 407                 {
 408                     z += box[axis][axis];
 409                 }
 410                 while (z > box[axis][axis])
 411                 {
 412                     z -= box[axis][axis];
 413                 }
 414
 415                 if (bRelative)
 416                 {
 417                     z = z/box[axis][axis];
 418                 }
 419
 420                 /* determine which slice atom is in */
 421                 if (bCenter)
 422                 {
 423                     slice = static_cast<int>(std::floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2);
 424                 }
 425                 else
 426                 {
 427                     slice = static_cast<int>(std::floor(z / (*slWidth)));
 428                 }
 429
 430                 /* Slice should already be 0<=slice<nslices, but we just make
 431                  * sure we are not hit by IEEE rounding errors since we do
 432                  * math operations after applying PBC above.
 433                  */
 434                 if (slice < 0)
 435                 {
 436                     slice += *nslices;
 437                 }
 438                 else if (slice >= *nslices)
 439                 {
 440                     slice -= *nslices;
 441                 }
 442
 443                 (*slDensity)[n][slice] += top->atoms.atom[index[n][i]].m*invvol;
 444             }
 445         }
 446         nr_frames++;
 447     }
 448     while (read_next_x(oenv, status, &t, x0, box));
 449     gmx_rmpbc_done(gpbc);
 450
 451     /*********** done with status file **********/
 452     close_trx(status);
 453
 454     /* slDensity now contains the total mass per slice, summed over all
 455        frames. Now divide by nr_frames and volume of slice
 456      */
 457
 458     fprintf(stderr, "\nRead %d frames from trajectory. Calculating density\n",
 459             nr_frames);
 460
 461     if (bRelative)
 462     {
 463         aveBox  /= nr_frames;
 464         *slWidth = aveBox/(*nslices);
 465     }
 466
 467     for (n = 0; n < nr_grps; n++)
 468     {
 469         for (i = 0; i < *nslices; i++)
 470         {
 471             (*slDensity)[n][i] /= nr_frames;
 472         }
 473     }
 474
 475     sfree(x0); /* free memory used by coordinate array */
 476 }
 477
 478 static void plot_density(double *slDensity[], const char *afile, int nslices,
 479                          int nr_grps, char *grpname[], real slWidth,
 480                          const char **dens_opt,
 481                          gmx_bool bCenter, gmx_bool bRelative, gmx_bool bSymmetrize,
 482                          const gmx_output_env_t *oenv)
 483 {
 484     FILE       *den;
 485     const char *title  = nullptr;
 486     const char *xlabel = nullptr;
 487     const char *ylabel = nullptr;
 488     int         slice, n;
 489     real        ddd;
 490     real        axispos;
 491
 492     title = bSymmetrize ? "Symmetrized partial density" : "Partial density";
 493
 494     if (bCenter)
 495     {
 496         xlabel = bRelative ?
 497             "Average relative position from center (nm)" :
 498             "Relative position from center (nm)";
 499     }
 500     else
 501     {
 502         xlabel = bRelative ? "Average coordinate (nm)" : "Coordinate (nm)";
 503     }
 504
 505     switch (dens_opt[0][0])
 506     {
 507         case 'm': ylabel = "Density (kg m\\S-3\\N)"; break;
 508         case 'n': ylabel = "Number density (nm\\S-3\\N)"; break;
 509         case 'c': ylabel = "Charge density (e nm\\S-3\\N)"; break;
 510         case 'e': ylabel = "Electron density (e nm\\S-3\\N)"; break;
 511     }
 512
 513     den = xvgropen(afile,
 514                    title, xlabel, ylabel, oenv);
 515
 516     xvgr_legend(den, nr_grps, (const char**)grpname, oenv);
 517
 518     for (slice = 0; (slice < nslices); slice++)
 519     {
 520         if (bCenter)
 521         {
 522             axispos = (slice - nslices/2.0 + 0.5) * slWidth;
 523         }
 524         else
 525         {
 526             axispos = (slice + 0.5) * slWidth;
 527         }
 528         fprintf(den, "%12g  ", axispos);
 529         for (n = 0; (n < nr_grps); n++)
 530         {
 531             if (bSymmetrize)
 532             {
 533                 ddd = (slDensity[n][slice]+slDensity[n][nslices-slice-1])*0.5;
 534             }
 535             else
 536             {
 537                 ddd = slDensity[n][slice];
 538             }
 539             if (dens_opt[0][0] == 'm')
 540             {
 541                 fprintf(den, "   %12g", ddd*AMU/(NANO*NANO*NANO));
 542             }
 543             else
 544             {
 545                 fprintf(den, "   %12g", ddd);
 546             }
 547         }
 548         fprintf(den, "\n");
 549     }
 550
 551     xvgrclose(den);
 552 }
 553
 554 int gmx_density(int argc, char *argv[])
 555 {
 556     const char        *desc[] = {
 557         "[THISMODULE] computes partial densities across the box, using an index file.[PAR]",
 558         "For the total density of NPT simulations, use [gmx-energy] instead.",
 559         "[PAR]",
 560
 561         "Option [TT]-center[tt] performs the histogram binning relative to the center",
 562         "of an arbitrary group, in absolute box coordinates. If you are calculating",
 563         "profiles along the Z axis box dimension bZ, output would be from -bZ/2 to",
 564         "bZ/2 if you center based on the entire system.",
 565         "Note that this behaviour has changed in GROMACS 5.0; earlier versions",
 566         "merely performed a static binning in (0,bZ) and shifted the output. Now",
 567         "we compute the center for each frame and bin in (-bZ/2,bZ/2).[PAR]",
 568
 569         "Option [TT]-symm[tt] symmetrizes the output around the center. This will",
 570         "automatically turn on [TT]-center[tt] too.",
 571
 572         "Option [TT]-relative[tt] performs the binning in relative instead of absolute",
 573         "box coordinates, and scales the final output with the average box dimension",
 574         "along the output axis. This can be used in combination with [TT]-center[tt].[PAR]",
 575
 576         "Densities are in kg/m^3, and number densities or electron densities can also be",
 577         "calculated. For electron densities, a file describing the number of",
 578         "electrons for each type of atom should be provided using [TT]-ei[tt].",
 579         "It should look like::",
 580         "",
 581         "   2",
 582         "   atomname = nrelectrons",
 583         "   atomname = nrelectrons",
 584         "",
 585         "The first line contains the number of lines to read from the file.",
 586         "There should be one line for each unique atom name in your system.",
 587         "The number of electrons for each atom is modified by its atomic",
 588         "partial charge.[PAR]",
 589
 590         "IMPORTANT CONSIDERATIONS FOR BILAYERS[PAR]",
 591         "One of the most common usage scenarios is to calculate the density of various",
 592         "groups across a lipid bilayer, typically with the z axis being the normal",
 593         "direction. For short simulations, small systems, and fixed box sizes this",
 594         "will work fine, but for the more general case lipid bilayers can be complicated.",
 595         "The first problem that while both proteins and lipids have low volume",
 596         "compressibility, lipids have quite high area compressiblity. This means the",
 597         "shape of the box (thickness and area/lipid) will fluctuate substantially even",
 598         "for a fully relaxed system. Since GROMACS places the box between the origin",
 599         "and positive coordinates, this in turn means that a bilayer centered in the",
 600         "box will move a bit up/down due to these fluctuations, and smear out your",
 601         "profile. The easiest way to fix this (if you want pressure coupling) is",
 602         "to use the [TT]-center[tt] option that calculates the density profile with",
 603         "respect to the center of the box. Note that you can still center on the",
 604         "bilayer part even if you have a complex non-symmetric system with a bilayer",
 605         "and, say, membrane proteins - then our output will simply have more values",
 606         "on one side of the (center) origin reference.[PAR]",
 607
 608         "Even the centered calculation will lead to some smearing out the output",
 609         "profiles, as lipids themselves are compressed and expanded. In most cases",
 610         "you probably want this (since it corresponds to macroscopic experiments),",
 611         "but if you want to look at molecular details you can use the [TT]-relative[tt]",
 612         "option to attempt to remove even more of the effects of volume fluctuations.[PAR]",
 613
 614         "Finally, large bilayers that are not subject to a surface tension will exhibit",
 615         "undulatory fluctuations, where there are 'waves' forming in the system.",
 616         "This is a fundamental property of the biological system, and if you are",
 617         "comparing against experiments you likely want to include the undulation",
 618         "smearing effect.",
 619         "",
 620     };
 621
 622     gmx_output_env_t  *oenv;
 623     static const char *dens_opt[] =
 624     { nullptr, "mass", "number", "charge", "electron", nullptr };
 625     static int         axis        = 2;  /* normal to memb. default z  */
 626     static const char *axtitle     = "Z";
 627     static int         nslices     = 50; /* nr of slices defined       */
 628     static int         ngrps       = 1;  /* nr. of groups              */
 629     static gmx_bool    bSymmetrize = FALSE;
 630     static gmx_bool    bCenter     = FALSE;
 631     static gmx_bool    bRelative   = FALSE;
 632
 633     t_pargs            pa[]        = {
 634         { "-d", FALSE, etSTR, {&axtitle},
 635           "Take the normal on the membrane in direction X, Y or Z." },
 636         { "-sl",  FALSE, etINT, {&nslices},
 637           "Divide the box in this number of slices." },
 638         { "-dens",    FALSE, etENUM, {dens_opt},
 639           "Density"},
 640         { "-ng",       FALSE, etINT, {&ngrps},
 641           "Number of groups of which to compute densities." },
 642         { "-center",   FALSE, etBOOL, {&bCenter},
 643           "Perform the binning relative to the center of the (changing) box. Useful for bilayers." },
 644         { "-symm",     FALSE, etBOOL, {&bSymmetrize},
 645           "Symmetrize the density along the axis, with respect to the center. Useful for bilayers." },
 646         { "-relative", FALSE, etBOOL, {&bRelative},
 647           "Use relative coordinates for changing boxes and scale output by average dimensions." }
 648     };
 649
 650     const char        *bugs[] = {
 651         "When calculating electron densities, atomnames are used instead of types. This is bad.",
 652     };
 653
 654     double           **density;        /* density per slice          */
 655     real               slWidth;        /* width of one slice         */
 656     char              *grpname_center; /* centering group name     */
 657     char             **grpname;        /* groupnames                 */
 658     int                nr_electrons;   /* nr. electrons              */
 659     int                ncenter;        /* size of centering group    */
 660     int               *ngx;            /* sizes of groups            */
 661     t_electron        *el_tab;         /* tabel with nr. of electrons*/
 662     t_topology        *top;            /* topology               */
 663     int                ePBC;
 664     int               *index_center;   /* index for centering group  */
 665     int              **index;          /* indices for all groups     */
 666     int                i;
 667
 668     t_filenm           fnm[] = { /* files for g_density       */
 669         { efTRX, "-f", nullptr,  ffREAD },
 670         { efNDX, nullptr, nullptr,  ffOPTRD },
 671         { efTPR, nullptr, nullptr,  ffREAD },
 672         { efDAT, "-ei", "electrons", ffOPTRD }, /* file with nr. of electrons */
 673         { efXVG, "-o", "density", ffWRITE },
 674     };
 675
 676 #define NFILE asize(fnm)
 677
 678     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
 679                            NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs,
 680                            &oenv))
 681     {
 682         return 0;
 683     }
 684
 685     GMX_RELEASE_ASSERT(dens_opt[0] != nullptr, "Option setting inconsistency; dens_opt[0] is NULL");
 686
 687     if (bSymmetrize && !bCenter)
 688     {
 689         fprintf(stderr, "Can not symmetrize without centering. Turning on -center\n");
 690         bCenter = TRUE;
 691     }
 692     /* Calculate axis */
 693     axis = toupper(axtitle[0]) - 'X';
 694
 695     top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC); /* read topology file */
 696     if (dens_opt[0][0] == 'n')
 697     {
 698         for (i = 0; (i < top->atoms.nr); i++)
 699         {
 700             top->atoms.atom[i].m = 1;
 701         }
 702     }
 703     else if (dens_opt[0][0] == 'c')
 704     {
 705         for (i = 0; (i < top->atoms.nr); i++)
 706         {
 707             top->atoms.atom[i].m = top->atoms.atom[i].q;
 708         }
 709     }
 710
 711     snew(grpname, ngrps);
 712     snew(index, ngrps);
 713     snew(ngx, ngrps);
 714
 715     if (bCenter)
 716     {
 717         fprintf(stderr,
 718                 "\nNote: that the center of mass is calculated inside the box without applying\n"
 719                 "any special periodicity. If necessary, it is your responsibility to first use\n"
 720                 "trjconv to make sure atoms in this group are placed in the right periodicity.\n\n"
 721                 "Select the group to center density profiles around:\n");
 722         get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &ncenter,
 723                   &index_center, &grpname_center);
 724     }
 725     else
 726     {
 727         ncenter      = 0;
 728         index_center = nullptr;
 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 }