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