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