src/gromacs/gmxana/gmx_density.c

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