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 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40 #include <math.h>
  41 #include <ctype.h>
  42
  43 #include "sysstuff.h"
  44 #include <string.h>
  45 #include "gromacs/utility/cstringutil.h"
  46 #include "typedefs.h"
  47 #include "gromacs/utility/smalloc.h"
  48 #include "macros.h"
  49 #include "gstat.h"
  50 #include "vec.h"
  51 #include "xvgr.h"
  52 #include "pbc.h"
  53 #include "gromacs/fileio/futil.h"
  54 #include "gromacs/commandline/pargs.h"
  55 #include "index.h"
  56 #include "gromacs/fileio/tpxio.h"
  57 #include "gromacs/fileio/trxio.h"
  58 #include "physics.h"
  59 #include "gmx_ana.h"
  60 #include "macros.h"
  61
  62 #include "gromacs/legacyheaders/gmx_fatal.h"
  63
  64 typedef struct {
  65     char *atomname;
  66     int   nr_el;
  67 } t_electron;
  68
  69 /****************************************************************************/
  70 /* This program calculates the partial density across the box.              */
  71 /* Peter Tieleman, Mei 1995                                                 */
  72 /****************************************************************************/
  73
  74 /* used for sorting the list */
  75 int compare(void *a, void *b)
  76 {
  77     t_electron *tmp1, *tmp2;
  78     tmp1 = (t_electron *)a; tmp2 = (t_electron *)b;
  79
  80     return strcmp(tmp1->atomname, tmp2->atomname);
  81 }
  82
  83 int get_electrons(t_electron **eltab, const char *fn)
  84 {
  85     char  buffer[256];  /* to read in a line   */
  86     char  tempname[80]; /* buffer to hold name */
  87     int   tempnr;
  88
  89     FILE *in;
  90     int   nr;        /* number of atomstypes to read */
  91     int   i;
  92
  93     if (!(in = gmx_ffopen(fn, "r")))
  94     {
  95         gmx_fatal(FARGS, "Couldn't open %s. Exiting.\n", fn);
  96     }
  97
  98     if (NULL == fgets(buffer, 255, in))
  99     {
 100         gmx_fatal(FARGS, "Error reading from file %s", fn);
 101     }
 102
 103     if (sscanf(buffer, "%d", &nr) != 1)
 104     {
 105         gmx_fatal(FARGS, "Invalid number of atomtypes in datafile\n");
 106     }
 107
 108     snew(*eltab, nr);
 109
 110     for (i = 0; i < nr; i++)
 111     {
 112         if (fgets(buffer, 255, in) == NULL)
 113         {
 114             gmx_fatal(FARGS, "reading datafile. Check your datafile.\n");
 115         }
 116         if (sscanf(buffer, "%s = %d", tempname, &tempnr) != 2)
 117         {
 118             gmx_fatal(FARGS, "Invalid line in datafile at line %d\n", i+1);
 119         }
 120         (*eltab)[i].nr_el    = tempnr;
 121         (*eltab)[i].atomname = strdup(tempname);
 122     }
 123     gmx_ffclose(in);
 124
 125     /* sort the list */
 126     fprintf(stderr, "Sorting list..\n");
 127     qsort ((void*)*eltab, nr, sizeof(t_electron),
 128            (int(*)(const void*, const void*))compare);
 129
 130     return nr;
 131 }
 132
 133 void center_coords(t_atoms *atoms, atom_id *index_center, int ncenter,
 134                    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).",
 148                       k+1, 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(box_center, com, 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 void calc_electron_density(const char *fn, atom_id **index, int gnx[],
 172                            double ***slDensity, int *nslices, t_topology *top,
 173                            int ePBC,
 174                            int axis, int nr_grps, real *slWidth,
 175                            t_electron eltab[], int nr, gmx_bool bCenter,
 176                            atom_id *index_center, int ncenter,
 177                            gmx_bool bRelative, const output_env_t oenv)
 178 {
 179     rvec        *x0;            /* coordinates without pbc */
 180     matrix       box;           /* box (3x3) */
 181     double       invvol;
 182     int          natoms;        /* nr. atoms in trj */
 183     t_trxstatus *status;
 184     int          i, n,          /* loop indices */
 185                  nr_frames = 0, /* number of frames */
 186                  slice;         /* current slice */
 187     t_electron  *found;         /* found by bsearch */
 188     t_electron   sought;        /* thingie thought by bsearch */
 189     real         boxSz, aveBox;
 190     gmx_rmpbc_t  gpbc = NULL;
 191
 192     real         t,
 193                  z;
 194
 195     if (axis < 0 || axis >= DIM)
 196     {
 197         gmx_fatal(FARGS, "Invalid axes. Terminating\n");
 198     }
 199
 200     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
 201     {
 202         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 203     }
 204
 205     aveBox = 0;
 206
 207     if (!*nslices)
 208     {
 209         *nslices = (int)(box[axis][axis] * 10); /* default value */
 210         fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
 211     }
 212
 213     snew(*slDensity, nr_grps);
 214     for (i = 0; i < nr_grps; i++)
 215     {
 216         snew((*slDensity)[i], *nslices);
 217     }
 218
 219     gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
 220     /*********** Start processing trajectory ***********/
 221     do
 222     {
 223         gmx_rmpbc(gpbc, natoms, box, x0);
 224
 225         /* Translate atoms so the com of the center-group is in the
 226          * box geometrical center.
 227          */
 228         if (bCenter)
 229         {
 230             center_coords(&top->atoms, index_center, ncenter, box, x0);
 231         }
 232
 233         invvol   = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
 234
 235         if (bRelative)
 236         {
 237             *slWidth = 1.0/(*nslices);
 238             boxSz    = 1.0;
 239         }
 240         else
 241         {
 242             *slWidth = box[axis][axis]/(*nslices);
 243             boxSz    = box[axis][axis];
 244         }
 245
 246         aveBox += box[axis][axis];
 247
 248         for (n = 0; n < nr_grps; n++)
 249         {
 250             for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
 251             {
 252                 z = x0[index[n][i]][axis];
 253                 while (z < 0)
 254                 {
 255                     z += box[axis][axis];
 256                 }
 257                 while (z > box[axis][axis])
 258                 {
 259                     z -= box[axis][axis];
 260                 }
 261
 262                 if (bRelative)
 263                 {
 264                     z = z/box[axis][axis];
 265                 }
 266
 267                 /* determine which slice atom is in */
 268                 if (bCenter)
 269                 {
 270                     slice = floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2;
 271                 }
 272                 else
 273                 {
 274                     slice = (z / (*slWidth));
 275                 }
 276                 sought.nr_el    = 0;
 277                 sought.atomname = strdup(*(top->atoms.atomname[index[n][i]]));
 278
 279                 /* now find the number of electrons. This is not efficient. */
 280                 found = (t_electron *)
 281                     bsearch((const void *)&sought,
 282                             (const void *)eltab, nr, sizeof(t_electron),
 283                             (int(*)(const void*, const void*))compare);
 284
 285                 if (found == NULL)
 286                 {
 287                     fprintf(stderr, "Couldn't find %s. Add it to the .dat file\n",
 288                             *(top->atoms.atomname[index[n][i]]));
 289                 }
 290                 else
 291                 {
 292                     (*slDensity)[n][slice] += (found->nr_el -
 293                                                top->atoms.atom[index[n][i]].q)*invvol;
 294                 }
 295                 free(sought.atomname);
 296             }
 297         }
 298         nr_frames++;
 299     }
 300     while (read_next_x(oenv, status, &t, x0, box));
 301     gmx_rmpbc_done(gpbc);
 302
 303     /*********** done with status file **********/
 304     close_trj(status);
 305
 306 /* slDensity now contains the total number of electrons per slice, summed
 307    over all frames. Now divide by nr_frames and volume of slice
 308  */
 309
 310     fprintf(stderr, "\nRead %d frames from trajectory. Counting electrons\n",
 311             nr_frames);
 312
 313     if (bRelative)
 314     {
 315         aveBox  /= nr_frames;
 316         *slWidth = aveBox/(*nslices);
 317     }
 318
 319     for (n = 0; n < nr_grps; n++)
 320     {
 321         for (i = 0; i < *nslices; i++)
 322         {
 323             (*slDensity)[n][i] /= nr_frames;
 324         }
 325     }
 326
 327     sfree(x0); /* free memory used by coordinate array */
 328 }
 329
 330 void calc_density(const char *fn, atom_id **index, int gnx[],
 331                   double ***slDensity, int *nslices, t_topology *top, int ePBC,
 332                   int axis, int nr_grps, real *slWidth, gmx_bool bCenter,
 333                   atom_id *index_center, int ncenter,
 334                   gmx_bool bRelative, const output_env_t oenv)
 335 {
 336     rvec        *x0;            /* coordinates without pbc */
 337     matrix       box;           /* box (3x3) */
 338     double       invvol;
 339     int          natoms;        /* nr. atoms in trj */
 340     t_trxstatus *status;
 341     int        **slCount,       /* nr. of atoms in one slice for a group */
 342                  i, j, n,       /* loop indices */
 343                  ax1       = 0, ax2 = 0,
 344                  nr_frames = 0, /* number of frames */
 345                  slice;         /* current slice */
 346     real         t,
 347                  z;
 348     real         boxSz, aveBox;
 349     char        *buf;    /* for tmp. keeping atomname */
 350     gmx_rmpbc_t  gpbc = NULL;
 351
 352     if (axis < 0 || axis >= DIM)
 353     {
 354         gmx_fatal(FARGS, "Invalid axes. Terminating\n");
 355     }
 356
 357     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
 358     {
 359         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 360     }
 361
 362     aveBox = 0;
 363
 364     if (!*nslices)
 365     {
 366         *nslices = (int)(box[axis][axis] * 10); /* default value */
 367         fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
 368     }
 369
 370     snew(*slDensity, nr_grps);
 371     for (i = 0; i < nr_grps; i++)
 372     {
 373         snew((*slDensity)[i], *nslices);
 374     }
 375
 376     gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
 377     /*********** Start processing trajectory ***********/
 378     do
 379     {
 380         gmx_rmpbc(gpbc, natoms, box, x0);
 381
 382         /* Translate atoms so the com of the center-group is in the
 383          * box geometrical center.
 384          */
 385         if (bCenter)
 386         {
 387             center_coords(&top->atoms, index_center, ncenter, box, x0);
 388         }
 389
 390         invvol   = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
 391
 392         if (bRelative)
 393         {
 394             *slWidth = 1.0/(*nslices);
 395             boxSz    = 1.0;
 396         }
 397         else
 398         {
 399             *slWidth = box[axis][axis]/(*nslices);
 400             boxSz    = box[axis][axis];
 401         }
 402
 403         aveBox += box[axis][axis];
 404
 405         for (n = 0; n < nr_grps; n++)
 406         {
 407             for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
 408             {
 409                 z = x0[index[n][i]][axis];
 410                 while (z < 0)
 411                 {
 412                     z += box[axis][axis];
 413                 }
 414                 while (z > box[axis][axis])
 415                 {
 416                     z -= box[axis][axis];
 417                 }
 418
 419                 if (bRelative)
 420                 {
 421                     z = z/box[axis][axis];
 422                 }
 423
 424                 /* determine which slice atom is in */
 425                 if (bCenter)
 426                 {
 427                     slice = floor( (z-(boxSz/2.0)) / (*slWidth) ) + *nslices/2;
 428                 }
 429                 else
 430                 {
 431                     slice = floor(z / (*slWidth));
 432                 }
 433
 434                 /* Slice should already be 0<=slice<nslices, but we just make
 435                  * sure we are not hit by IEEE rounding errors since we do
 436                  * math operations after applying PBC above.
 437                  */
 438                 if (slice < 0)
 439                 {
 440                     slice += *nslices;
 441                 }
 442                 else if (slice >= *nslices)
 443                 {
 444                     slice -= *nslices;
 445                 }
 446
 447                 (*slDensity)[n][slice] += top->atoms.atom[index[n][i]].m*invvol;
 448             }
 449         }
 450         nr_frames++;
 451     }
 452     while (read_next_x(oenv, status, &t, x0, box));
 453     gmx_rmpbc_done(gpbc);
 454
 455     /*********** done with status file **********/
 456     close_trj(status);
 457
 458     /* slDensity now contains the total mass per slice, summed over all
 459        frames. Now divide by nr_frames and volume of slice
 460      */
 461
 462     fprintf(stderr, "\nRead %d frames from trajectory. Calculating density\n",
 463             nr_frames);
 464
 465     if (bRelative)
 466     {
 467         aveBox  /= nr_frames;
 468         *slWidth = aveBox/(*nslices);
 469     }
 470
 471     for (n = 0; n < nr_grps; n++)
 472     {
 473         for (i = 0; i < *nslices; i++)
 474         {
 475             (*slDensity)[n][i] /= nr_frames;
 476         }
 477     }
 478
 479     sfree(x0); /* free memory used by coordinate array */
 480 }
 481
 482 void plot_density(double *slDensity[], const char *afile, int nslices,
 483                   int nr_grps, char *grpname[], real slWidth,
 484                   const char **dens_opt,
 485                   gmx_bool bCenter, gmx_bool bRelative, gmx_bool bSymmetrize,
 486                   const output_env_t oenv)
 487 {
 488     FILE       *den;
 489     const char *title  = NULL;
 490     const char *xlabel = NULL;
 491     const char *ylabel = NULL;
 492     int         slice, n;
 493     real        ddd;
 494     real        axispos;
 495
 496     title = bSymmetrize ? "Symmetrized partial density" : "Partial density";
 497
 498     if (bCenter)
 499     {
 500         xlabel = bRelative ?
 501             "Average relative position from center (nm)" :
 502             "Relative position from center (nm)";
 503     }
 504     else
 505     {
 506         xlabel = bRelative ? "Average coordinate (nm)" : "Coordinate (nm)";
 507     }
 508
 509     switch (dens_opt[0][0])
 510     {
 511         case 'm': ylabel = "Density (kg m\\S-3\\N)"; break;
 512         case 'n': ylabel = "Number density (nm\\S-3\\N)"; break;
 513         case 'c': ylabel = "Charge density (e nm\\S-3\\N)"; break;
 514         case 'e': ylabel = "Electron density (e nm\\S-3\\N)"; break;
 515     }
 516
 517     den = xvgropen(afile,
 518                    title, xlabel, ylabel, oenv);
 519
 520     xvgr_legend(den, nr_grps, (const char**)grpname, oenv);
 521
 522     for (slice = 0; (slice < nslices); slice++)
 523     {
 524         if (bCenter)
 525         {
 526             axispos = (slice - nslices/2.0 + 0.5) * slWidth;
 527         }
 528         else
 529         {
 530             axispos = (slice + 0.5) * slWidth;
 531         }
 532         fprintf(den, "%12g  ", axispos);
 533         for (n = 0; (n < nr_grps); n++)
 534         {
 535             if (bSymmetrize)
 536             {
 537                 ddd = (slDensity[n][slice]+slDensity[n][nslices-slice-1])*0.5;
 538             }
 539             else
 540             {
 541                 ddd = slDensity[n][slice];
 542             }
 543             if (dens_opt[0][0] == 'm')
 544             {
 545                 fprintf(den, "   %12g", ddd*AMU/(NANO*NANO*NANO));
 546             }
 547             else
 548             {
 549                 fprintf(den, "   %12g", ddd);
 550             }
 551         }
 552         fprintf(den, "\n");
 553     }
 554
 555     gmx_ffclose(den);
 556 }
 557
 558 int gmx_density(int argc, char *argv[])
 559 {
 560     const char        *desc[] = {
 561         "[THISMODULE] computes partial densities across the box, using an index file.[PAR]",
 562         "For the total density of NPT simulations, use [gmx-energy] instead.",
 563         "[PAR]",
 564
 565         "Option [TT]-center[tt] performs the histogram binning relative to the center",
 566         "of an arbitrary group, in absolute box coordinates. If you are calculating",
 567         "profiles along the Z axis box dimension bZ, output would be from -bZ/2 to",
 568         "bZ/2 if you center based on the entire system.",
 569         "Note that this behaviour has changed in Gromacs 5.0; earlier versions",
 570         "merely performed a static binning in (0,bZ) and shifted the output. Now",
 571         "we compute the center for each frame and bin in (-bZ/2,bZ/2).[PAR]",
 572
 573         "Option [TT]-symm[tt] symmetrizes the output around the center. This will",
 574         "automatically turn on [TT]-center[tt] too.",
 575
 576         "Option [TT]-relative[tt] performs the binning in relative instead of absolute",
 577         "box coordinates, and scales the final output with the average box dimension",
 578         "along the output axis. This can be used in combination with [TT]-center[tt].[PAR]",
 579
 580         "Densities are in kg/m^3, and number densities or electron densities can also be",
 581         "calculated. For electron densities, a file describing the number of",
 582         "electrons for each type of atom should be provided using [TT]-ei[tt].",
 583         "It should look like:[BR]",
 584         "   [TT]2[tt][BR]",
 585         "   [TT]atomname = nrelectrons[tt][BR]",
 586         "   [TT]atomname = nrelectrons[tt][BR]",
 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         { efTPX, 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 | PCA_BE_NICE,
 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(efTPX, 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 }