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