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,2018, 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 = static_cast<t_electron *>(a); tmp2 = static_cast<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            reinterpret_cast<int(*)(const void*, const void*)>(compare));
 128
 129     return nr;
 130 }
 131
 132 static void center_coords(t_atoms *atoms, const 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, const 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 = static_cast<t_electron *>(bsearch((const void *)&sought,
 280                                                           (const void *)eltab, nr, sizeof(t_electron),
 281                                                           reinterpret_cast<int(*)(const void*, const void*)>(compare)));
 282
 283                 if (found == nullptr)
 284                 {
 285                     fprintf(stderr, "Couldn't find %s. Add it to the .dat file\n",
 286                             *(top->atoms.atomname[index[n][i]]));
 287                 }
 288                 else
 289                 {
 290                     (*slDensity)[n][slice] += (found->nr_el -
 291                                                top->atoms.atom[index[n][i]].q)*invvol;
 292                 }
 293                 free(sought.atomname);
 294             }
 295         }
 296         nr_frames++;
 297     }
 298     while (read_next_x(oenv, status, &t, x0, box));
 299     gmx_rmpbc_done(gpbc);
 300
 301     /*********** done with status file **********/
 302     close_trx(status);
 303
 304 /* slDensity now contains the total number of electrons per slice, summed
 305    over all frames. Now divide by nr_frames and volume of slice
 306  */
 307
 308     fprintf(stderr, "\nRead %d frames from trajectory. Counting electrons\n",
 309             nr_frames);
 310
 311     if (bRelative)
 312     {
 313         aveBox  /= nr_frames;
 314         *slWidth = aveBox/(*nslices);
 315     }
 316
 317     for (n = 0; n < nr_grps; n++)
 318     {
 319         for (i = 0; i < *nslices; i++)
 320         {
 321             (*slDensity)[n][i] /= nr_frames;
 322         }
 323     }
 324
 325     sfree(x0); /* free memory used by coordinate array */
 326 }
 327
 328 static void calc_density(const char *fn, int **index, const int gnx[],
 329                          double ***slDensity, int *nslices, t_topology *top, int ePBC,
 330                          int axis, int nr_grps, real *slWidth, gmx_bool bCenter,
 331                          int *index_center, int ncenter,
 332                          gmx_bool bRelative, const gmx_output_env_t *oenv, const char **dens_opt)
 333 {
 334     rvec        *x0;            /* coordinates without pbc */
 335     matrix       box;           /* box (3x3) */
 336     double       invvol;
 337     int          natoms;        /* nr. atoms in trj */
 338     t_trxstatus *status;
 339     int          i, n,          /* loop indices */
 340                  nr_frames = 0, /* number of frames */
 341                  slice;         /* current slice */
 342     real         t,
 343                  z;
 344     real         boxSz, aveBox;
 345     real        *den_val; /* values from which the density is calculated */
 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
 375     snew(den_val, top->atoms.nr);
 376     if (dens_opt[0][0] == 'n')
 377     {
 378         for (i = 0; (i < top->atoms.nr); i++)
 379         {
 380             den_val[i] = 1;
 381         }
 382     }
 383     else if (dens_opt[0][0] == 'c')
 384     {
 385         for (i = 0; (i < top->atoms.nr); i++)
 386         {
 387             den_val[i] = top->atoms.atom[i].q;
 388         }
 389     }
 390     else
 391     {
 392         for (i = 0; (i < top->atoms.nr); i++)
 393         {
 394             den_val[i] = top->atoms.atom[i].m;
 395         }
 396     }
 397
 398     do
 399     {
 400         gmx_rmpbc(gpbc, natoms, box, x0);
 401
 402         /* Translate atoms so the com of the center-group is in the
 403          * box geometrical center.
 404          */
 405         if (bCenter)
 406         {
 407             center_coords(&top->atoms, index_center, ncenter, box, x0);
 408         }
 409
 410         invvol   = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
 411
 412         if (bRelative)
 413         {
 414             *slWidth = 1.0/(*nslices);
 415             boxSz    = 1.0;
 416         }
 417         else
 418         {
 419             *slWidth = box[axis][axis]/(*nslices);
 420             boxSz    = box[axis][axis];
 421         }
 422
 423         aveBox += box[axis][axis];
 424
 425         for (n = 0; n < nr_grps; n++)
 426         {
 427             for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
 428             {
 429                 z = x0[index[n][i]][axis];
 430                 while (z < 0)
 431                 {
 432                     z += box[axis][axis];
 433                 }
 434                 while (z > box[axis][axis])
 435                 {
 436                     z -= box[axis][axis];
 437                 }
 438
 439                 if (bRelative)
 440                 {
 441                     z = z/box[axis][axis];
 442                 }
 443
 444                 /* determine which slice atom is in */
 445                 if (bCenter)
 446                 {
 447                     slice = static_cast<int>(std::floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2.);
 448                 }
 449                 else
 450                 {
 451                     slice = static_cast<int>(std::floor(z / (*slWidth)));
 452                 }
 453
 454                 /* Slice should already be 0<=slice<nslices, but we just make
 455                  * sure we are not hit by IEEE rounding errors since we do
 456                  * math operations after applying PBC above.
 457                  */
 458                 if (slice < 0)
 459                 {
 460                     slice += *nslices;
 461                 }
 462                 else if (slice >= *nslices)
 463                 {
 464                     slice -= *nslices;
 465                 }
 466
 467                 (*slDensity)[n][slice] += den_val[index[n][i]]*invvol;
 468             }
 469         }
 470         nr_frames++;
 471     }
 472     while (read_next_x(oenv, status, &t, x0, box));
 473     gmx_rmpbc_done(gpbc);
 474
 475     /*********** done with status file **********/
 476     close_trx(status);
 477
 478     /* slDensity now contains the total mass per slice, summed over all
 479        frames. Now divide by nr_frames and volume of slice
 480      */
 481
 482     fprintf(stderr, "\nRead %d frames from trajectory. Calculating density\n",
 483             nr_frames);
 484
 485     if (bRelative)
 486     {
 487         aveBox  /= nr_frames;
 488         *slWidth = aveBox/(*nslices);
 489     }
 490
 491     for (n = 0; n < nr_grps; n++)
 492     {
 493         for (i = 0; i < *nslices; i++)
 494         {
 495             (*slDensity)[n][i] /= nr_frames;
 496         }
 497     }
 498
 499     sfree(x0); /* free memory used by coordinate array */
 500     sfree(den_val);
 501 }
 502
 503 static void plot_density(double *slDensity[], const char *afile, int nslices,
 504                          int nr_grps, char *grpname[], real slWidth,
 505                          const char **dens_opt,
 506                          gmx_bool bCenter, gmx_bool bRelative, gmx_bool bSymmetrize,
 507                          const gmx_output_env_t *oenv)
 508 {
 509     FILE       *den;
 510     const char *title  = nullptr;
 511     const char *xlabel = nullptr;
 512     const char *ylabel = nullptr;
 513     int         slice, n;
 514     real        ddd;
 515     real        axispos;
 516
 517     title = bSymmetrize ? "Symmetrized partial density" : "Partial density";
 518
 519     if (bCenter)
 520     {
 521         xlabel = bRelative ?
 522             "Average relative position from center (nm)" :
 523             "Relative position from center (nm)";
 524     }
 525     else
 526     {
 527         xlabel = bRelative ? "Average coordinate (nm)" : "Coordinate (nm)";
 528     }
 529
 530     switch (dens_opt[0][0])
 531     {
 532         case 'm': ylabel = "Density (kg m\\S-3\\N)"; break;
 533         case 'n': ylabel = "Number density (nm\\S-3\\N)"; break;
 534         case 'c': ylabel = "Charge density (e nm\\S-3\\N)"; break;
 535         case 'e': ylabel = "Electron density (e nm\\S-3\\N)"; break;
 536     }
 537
 538     den = xvgropen(afile,
 539                    title, xlabel, ylabel, oenv);
 540
 541     xvgr_legend(den, nr_grps, (const char**)grpname, oenv);
 542
 543     for (slice = 0; (slice < nslices); slice++)
 544     {
 545         if (bCenter)
 546         {
 547             axispos = (slice - nslices/2.0 + 0.5) * slWidth;
 548         }
 549         else
 550         {
 551             axispos = (slice + 0.5) * slWidth;
 552         }
 553         fprintf(den, "%12g  ", axispos);
 554         for (n = 0; (n < nr_grps); n++)
 555         {
 556             if (bSymmetrize)
 557             {
 558                 ddd = (slDensity[n][slice]+slDensity[n][nslices-slice-1])*0.5;
 559             }
 560             else
 561             {
 562                 ddd = slDensity[n][slice];
 563             }
 564             if (dens_opt[0][0] == 'm')
 565             {
 566                 fprintf(den, "   %12g", ddd*AMU/(NANO*NANO*NANO));
 567             }
 568             else
 569             {
 570                 fprintf(den, "   %12g", ddd);
 571             }
 572         }
 573         fprintf(den, "\n");
 574     }
 575
 576     xvgrclose(den);
 577 }
 578
 579 int gmx_density(int argc, char *argv[])
 580 {
 581     const char        *desc[] = {
 582         "[THISMODULE] computes partial densities across the box, using an index file.[PAR]",
 583         "For the total density of NPT simulations, use [gmx-energy] instead.",
 584         "[PAR]",
 585
 586         "Option [TT]-center[tt] performs the histogram binning relative to the center",
 587         "of an arbitrary group, in absolute box coordinates. If you are calculating",
 588         "profiles along the Z axis box dimension bZ, output would be from -bZ/2 to",
 589         "bZ/2 if you center based on the entire system.",
 590         "Note that this behaviour has changed in GROMACS 5.0; earlier versions",
 591         "merely performed a static binning in (0,bZ) and shifted the output. Now",
 592         "we compute the center for each frame and bin in (-bZ/2,bZ/2).[PAR]",
 593
 594         "Option [TT]-symm[tt] symmetrizes the output around the center. This will",
 595         "automatically turn on [TT]-center[tt] too.",
 596
 597         "Option [TT]-relative[tt] performs the binning in relative instead of absolute",
 598         "box coordinates, and scales the final output with the average box dimension",
 599         "along the output axis. This can be used in combination with [TT]-center[tt].[PAR]",
 600
 601         "Densities are in kg/m^3, and number densities or electron densities can also be",
 602         "calculated. For electron densities, a file describing the number of",
 603         "electrons for each type of atom should be provided using [TT]-ei[tt].",
 604         "It should look like::",
 605         "",
 606         "   2",
 607         "   atomname = nrelectrons",
 608         "   atomname = nrelectrons",
 609         "",
 610         "The first line contains the number of lines to read from the file.",
 611         "There should be one line for each unique atom name in your system.",
 612         "The number of electrons for each atom is modified by its atomic",
 613         "partial charge.[PAR]",
 614
 615         "IMPORTANT CONSIDERATIONS FOR BILAYERS[PAR]",
 616         "One of the most common usage scenarios is to calculate the density of various",
 617         "groups across a lipid bilayer, typically with the z axis being the normal",
 618         "direction. For short simulations, small systems, and fixed box sizes this",
 619         "will work fine, but for the more general case lipid bilayers can be complicated.",
 620         "The first problem that while both proteins and lipids have low volume",
 621         "compressibility, lipids have quite high area compressiblity. This means the",
 622         "shape of the box (thickness and area/lipid) will fluctuate substantially even",
 623         "for a fully relaxed system. Since GROMACS places the box between the origin",
 624         "and positive coordinates, this in turn means that a bilayer centered in the",
 625         "box will move a bit up/down due to these fluctuations, and smear out your",
 626         "profile. The easiest way to fix this (if you want pressure coupling) is",
 627         "to use the [TT]-center[tt] option that calculates the density profile with",
 628         "respect to the center of the box. Note that you can still center on the",
 629         "bilayer part even if you have a complex non-symmetric system with a bilayer",
 630         "and, say, membrane proteins - then our output will simply have more values",
 631         "on one side of the (center) origin reference.[PAR]",
 632
 633         "Even the centered calculation will lead to some smearing out the output",
 634         "profiles, as lipids themselves are compressed and expanded. In most cases",
 635         "you probably want this (since it corresponds to macroscopic experiments),",
 636         "but if you want to look at molecular details you can use the [TT]-relative[tt]",
 637         "option to attempt to remove even more of the effects of volume fluctuations.[PAR]",
 638
 639         "Finally, large bilayers that are not subject to a surface tension will exhibit",
 640         "undulatory fluctuations, where there are 'waves' forming in the system.",
 641         "This is a fundamental property of the biological system, and if you are",
 642         "comparing against experiments you likely want to include the undulation",
 643         "smearing effect.",
 644         "",
 645     };
 646
 647     gmx_output_env_t  *oenv;
 648     static const char *dens_opt[] =
 649     { nullptr, "mass", "number", "charge", "electron", nullptr };
 650     static int         axis        = 2;  /* normal to memb. default z  */
 651     static const char *axtitle     = "Z";
 652     static int         nslices     = 50; /* nr of slices defined       */
 653     static int         ngrps       = 1;  /* nr. of groups              */
 654     static gmx_bool    bSymmetrize = FALSE;
 655     static gmx_bool    bCenter     = FALSE;
 656     static gmx_bool    bRelative   = FALSE;
 657
 658     t_pargs            pa[]        = {
 659         { "-d", FALSE, etSTR, {&axtitle},
 660           "Take the normal on the membrane in direction X, Y or Z." },
 661         { "-sl",  FALSE, etINT, {&nslices},
 662           "Divide the box in this number of slices." },
 663         { "-dens",    FALSE, etENUM, {dens_opt},
 664           "Density"},
 665         { "-ng",       FALSE, etINT, {&ngrps},
 666           "Number of groups of which to compute densities." },
 667         { "-center",   FALSE, etBOOL, {&bCenter},
 668           "Perform the binning relative to the center of the (changing) box. Useful for bilayers." },
 669         { "-symm",     FALSE, etBOOL, {&bSymmetrize},
 670           "Symmetrize the density along the axis, with respect to the center. Useful for bilayers." },
 671         { "-relative", FALSE, etBOOL, {&bRelative},
 672           "Use relative coordinates for changing boxes and scale output by average dimensions." }
 673     };
 674
 675     const char        *bugs[] = {
 676         "When calculating electron densities, atomnames are used instead of types. This is bad.",
 677     };
 678
 679     double           **density;        /* density per slice          */
 680     real               slWidth;        /* width of one slice         */
 681     char              *grpname_center; /* centering group name     */
 682     char             **grpname;        /* groupnames                 */
 683     int                nr_electrons;   /* nr. electrons              */
 684     int                ncenter;        /* size of centering group    */
 685     int               *ngx;            /* sizes of groups            */
 686     t_electron        *el_tab;         /* tabel with nr. of electrons*/
 687     t_topology        *top;            /* topology               */
 688     int                ePBC;
 689     int               *index_center;   /* index for centering group  */
 690     int              **index;          /* indices for all groups     */
 691
 692     t_filenm           fnm[] = {       /* files for g_density       */
 693         { efTRX, "-f", nullptr,  ffREAD },
 694         { efNDX, nullptr, nullptr,  ffOPTRD },
 695         { efTPR, nullptr, nullptr,  ffREAD },
 696         { efDAT, "-ei", "electrons", ffOPTRD }, /* file with nr. of electrons */
 697         { efXVG, "-o", "density", ffWRITE },
 698     };
 699
 700 #define NFILE asize(fnm)
 701
 702     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
 703                            NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs,
 704                            &oenv))
 705     {
 706         return 0;
 707     }
 708
 709     GMX_RELEASE_ASSERT(dens_opt[0] != nullptr, "Option setting inconsistency; dens_opt[0] is NULL");
 710
 711     if (bSymmetrize && !bCenter)
 712     {
 713         fprintf(stderr, "Can not symmetrize without centering. Turning on -center\n");
 714         bCenter = TRUE;
 715     }
 716     /* Calculate axis */
 717     axis = toupper(axtitle[0]) - 'X';
 718
 719     top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC); /* read topology file */
 720
 721     snew(grpname, ngrps);
 722     snew(index, ngrps);
 723     snew(ngx, ngrps);
 724
 725     if (bCenter)
 726     {
 727         fprintf(stderr,
 728                 "\nNote: that the center of mass is calculated inside the box without applying\n"
 729                 "any special periodicity. If necessary, it is your responsibility to first use\n"
 730                 "trjconv to make sure atoms in this group are placed in the right periodicity.\n\n"
 731                 "Select the group to center density profiles around:\n");
 732         get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &ncenter,
 733                   &index_center, &grpname_center);
 734     }
 735     else
 736     {
 737         ncenter      = 0;
 738         index_center = nullptr;
 739     }
 740
 741     fprintf(stderr, "\nSelect %d group%s to calculate density for:\n", ngrps, (ngrps > 1) ? "s" : "");
 742     get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, ngx, index, grpname);
 743
 744     if (dens_opt[0][0] == 'e')
 745     {
 746         nr_electrons =  get_electrons(&el_tab, ftp2fn(efDAT, NFILE, fnm));
 747         fprintf(stderr, "Read %d atomtypes from datafile\n", nr_electrons);
 748
 749         calc_electron_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density,
 750                               &nslices, top, ePBC, axis, ngrps, &slWidth, el_tab,
 751                               nr_electrons, bCenter, index_center, ncenter,
 752                               bRelative, oenv);
 753     }
 754     else
 755     {
 756         calc_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density, &nslices, top,
 757                      ePBC, axis, ngrps, &slWidth, bCenter, index_center, ncenter,
 758                      bRelative, oenv, dens_opt);
 759     }
 760
 761     plot_density(density, opt2fn("-o", NFILE, fnm),
 762                  nslices, ngrps, grpname, slWidth, dens_opt,
 763                  bCenter, bRelative, bSymmetrize, oenv);
 764
 765     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");  /* view xvgr file */
 766     return 0;
 767 }