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