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, matrix box, rvec x0[], int axis)
 134 {
 135     int  i, m;
 136     real tmass, mm;
 137     rvec com, shift, box_center;
 138
 139     tmass = 0;
 140     clear_rvec(com);
 141     for (i = 0; (i < atoms->nr); i++)
 142     {
 143         mm     = atoms->atom[i].m;
 144         tmass += mm;
 145         for (m = 0; (m < DIM); m++)
 146         {
 147             com[m] += mm*x0[i][m];
 148         }
 149     }
 150     for (m = 0; (m < DIM); m++)
 151     {
 152         com[m] /= tmass;
 153     }
 154     calc_box_center(ecenterDEF, box, box_center);
 155     rvec_sub(box_center, com, shift);
 156     shift[axis] -= box_center[axis];
 157
 158     for (i = 0; (i < atoms->nr); i++)
 159     {
 160         rvec_dec(x0[i], shift);
 161     }
 162 }
 163
 164 void calc_electron_density(const char *fn, atom_id **index, int gnx[],
 165                            double ***slDensity, int *nslices, t_topology *top,
 166                            int ePBC,
 167                            int axis, int nr_grps, real *slWidth,
 168                            t_electron eltab[], int nr, gmx_bool bCenter,
 169                            const output_env_t oenv)
 170 {
 171     rvec        *x0;            /* coordinates without pbc */
 172     matrix       box;           /* box (3x3) */
 173     double       invvol;
 174     int          natoms;        /* nr. atoms in trj */
 175     t_trxstatus *status;
 176     int          i, n,          /* loop indices */
 177                  nr_frames = 0, /* number of frames */
 178                  slice;         /* current slice */
 179     t_electron  *found;         /* found by bsearch */
 180     t_electron   sought;        /* thingie thought by bsearch */
 181     gmx_rmpbc_t  gpbc = NULL;
 182
 183     real         t,
 184                  z;
 185
 186     if (axis < 0 || axis >= DIM)
 187     {
 188         gmx_fatal(FARGS, "Invalid axes. Terminating\n");
 189     }
 190
 191     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
 192     {
 193         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 194     }
 195
 196     if (!*nslices)
 197     {
 198         *nslices = (int)(box[axis][axis] * 10); /* default value */
 199     }
 200     fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
 201
 202     snew(*slDensity, nr_grps);
 203     for (i = 0; i < nr_grps; i++)
 204     {
 205         snew((*slDensity)[i], *nslices);
 206     }
 207
 208     gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
 209     /*********** Start processing trajectory ***********/
 210     do
 211     {
 212         gmx_rmpbc(gpbc, natoms, box, x0);
 213
 214         if (bCenter)
 215         {
 216             center_coords(&top->atoms, box, x0, axis);
 217         }
 218
 219         *slWidth = box[axis][axis]/(*nslices);
 220         invvol   = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
 221
 222         for (n = 0; n < nr_grps; n++)
 223         {
 224             for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
 225             {
 226                 z = x0[index[n][i]][axis];
 227                 while (z < 0)
 228                 {
 229                     z += box[axis][axis];
 230                 }
 231                 while (z > box[axis][axis])
 232                 {
 233                     z -= box[axis][axis];
 234                 }
 235
 236                 /* determine which slice atom is in */
 237                 slice           = (z / (*slWidth));
 238                 sought.nr_el    = 0;
 239                 sought.atomname = strdup(*(top->atoms.atomname[index[n][i]]));
 240
 241                 /* now find the number of electrons. This is not efficient. */
 242                 found = (t_electron *)
 243                     bsearch((const void *)&sought,
 244                             (const void *)eltab, nr, sizeof(t_electron),
 245                             (int(*)(const void*, const void*))compare);
 246
 247                 if (found == NULL)
 248                 {
 249                     fprintf(stderr, "Couldn't find %s. Add it to the .dat file\n",
 250                             *(top->atoms.atomname[index[n][i]]));
 251                 }
 252                 else
 253                 {
 254                     (*slDensity)[n][slice] += (found->nr_el -
 255                                                top->atoms.atom[index[n][i]].q)*invvol;
 256                 }
 257                 free(sought.atomname);
 258             }
 259         }
 260         nr_frames++;
 261     }
 262     while (read_next_x(oenv, status, &t, x0, box));
 263     gmx_rmpbc_done(gpbc);
 264
 265     /*********** done with status file **********/
 266     close_trj(status);
 267
 268 /* slDensity now contains the total number of electrons per slice, summed
 269    over all frames. Now divide by nr_frames and volume of slice
 270  */
 271
 272     fprintf(stderr, "\nRead %d frames from trajectory. Counting electrons\n",
 273             nr_frames);
 274
 275     for (n = 0; n < nr_grps; n++)
 276     {
 277         for (i = 0; i < *nslices; i++)
 278         {
 279             (*slDensity)[n][i] /= nr_frames;
 280         }
 281     }
 282
 283     sfree(x0); /* free memory used by coordinate array */
 284 }
 285
 286 void calc_density(const char *fn, atom_id **index, int gnx[],
 287                   double ***slDensity, int *nslices, t_topology *top, int ePBC,
 288                   int axis, int nr_grps, real *slWidth, gmx_bool bCenter,
 289                   const output_env_t oenv)
 290 {
 291     rvec        *x0;      /* coordinates without pbc */
 292     matrix       box;     /* box (3x3) */
 293     double       invvol;
 294     int          natoms;  /* nr. atoms in trj */
 295     t_trxstatus *status;
 296     int        **slCount, /* nr. of atoms in one slice for a group */
 297                  i, j, n, /* loop indices */
 298                  teller    = 0,
 299                  ax1       = 0, ax2 = 0,
 300                  nr_frames = 0, /* number of frames */
 301                  slice;         /* current slice */
 302     real         t,
 303                  z;
 304     char        *buf;    /* for tmp. keeping atomname */
 305     gmx_rmpbc_t  gpbc = NULL;
 306
 307     if (axis < 0 || axis >= DIM)
 308     {
 309         gmx_fatal(FARGS, "Invalid axes. Terminating\n");
 310     }
 311
 312     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
 313     {
 314         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 315     }
 316
 317     if (!*nslices)
 318     {
 319         *nslices = (int)(box[axis][axis] * 10); /* default value */
 320         fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
 321     }
 322
 323     snew(*slDensity, nr_grps);
 324     for (i = 0; i < nr_grps; i++)
 325     {
 326         snew((*slDensity)[i], *nslices);
 327     }
 328
 329     gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
 330     /*********** Start processing trajectory ***********/
 331     do
 332     {
 333         gmx_rmpbc(gpbc, natoms, box, x0);
 334
 335         if (bCenter)
 336         {
 337             center_coords(&top->atoms, box, x0, axis);
 338         }
 339
 340         *slWidth = box[axis][axis]/(*nslices);
 341         invvol   = *nslices/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
 342         teller++;
 343
 344         for (n = 0; n < nr_grps; n++)
 345         {
 346             for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
 347             {
 348                 z = x0[index[n][i]][axis];
 349                 while (z < 0)
 350                 {
 351                     z += box[axis][axis];
 352                 }
 353                 while (z > box[axis][axis])
 354                 {
 355                     z -= box[axis][axis];
 356                 }
 357
 358                 /* determine which slice atom is in */
 359                 slice                   = (int)(z / (*slWidth));
 360                 (*slDensity)[n][slice] += top->atoms.atom[index[n][i]].m*invvol;
 361             }
 362         }
 363         nr_frames++;
 364     }
 365     while (read_next_x(oenv, status, &t, x0, box));
 366     gmx_rmpbc_done(gpbc);
 367
 368     /*********** done with status file **********/
 369     close_trj(status);
 370
 371     /* slDensity now contains the total mass per slice, summed over all
 372        frames. Now divide by nr_frames and volume of slice
 373      */
 374
 375     fprintf(stderr, "\nRead %d frames from trajectory. Calculating density\n",
 376             nr_frames);
 377
 378     for (n = 0; n < nr_grps; n++)
 379     {
 380         for (i = 0; i < *nslices; i++)
 381         {
 382             (*slDensity)[n][i] /= nr_frames;
 383         }
 384     }
 385
 386     sfree(x0); /* free memory used by coordinate array */
 387 }
 388
 389 void plot_density(double *slDensity[], const char *afile, int nslices,
 390                   int nr_grps, char *grpname[], real slWidth,
 391                   const char **dens_opt,
 392                   gmx_bool bSymmetrize, const output_env_t oenv)
 393 {
 394     FILE       *den;
 395     const char *ylabel = NULL;
 396     int         slice, n;
 397     real        ddd;
 398
 399     switch (dens_opt[0][0])
 400     {
 401         case 'm': ylabel = "Density (kg m\\S-3\\N)"; break;
 402         case 'n': ylabel = "Number density (nm\\S-3\\N)"; break;
 403         case 'c': ylabel = "Charge density (e nm\\S-3\\N)"; break;
 404         case 'e': ylabel = "Electron density (e nm\\S-3\\N)"; break;
 405     }
 406
 407     den = xvgropen(afile, "Partial densities", "Box (nm)", ylabel, oenv);
 408
 409     xvgr_legend(den, nr_grps, (const char**)grpname, oenv);
 410
 411     for (slice = 0; (slice < nslices); slice++)
 412     {
 413         fprintf(den, "%12g  ", slice * slWidth);
 414         for (n = 0; (n < nr_grps); n++)
 415         {
 416             if (bSymmetrize)
 417             {
 418                 ddd = (slDensity[n][slice]+slDensity[n][nslices-slice-1])*0.5;
 419             }
 420             else
 421             {
 422                 ddd = slDensity[n][slice];
 423             }
 424             if (dens_opt[0][0] == 'm')
 425             {
 426                 fprintf(den, "   %12g", ddd*AMU/(NANO*NANO*NANO));
 427             }
 428             else
 429             {
 430                 fprintf(den, "   %12g", ddd);
 431             }
 432         }
 433         fprintf(den, "\n");
 434     }
 435
 436     gmx_ffclose(den);
 437 }
 438
 439 int gmx_density(int argc, char *argv[])
 440 {
 441     const char        *desc[] = {
 442         "[THISMODULE] computes partial densities across the box, using an index file.[PAR]",
 443         "For the total density of NPT simulations, use [gmx-energy] instead.",
 444         "[PAR]",
 445         "Densities are in kg/m^3, and number densities or electron densities can also be",
 446         "calculated. For electron densities, a file describing the number of",
 447         "electrons for each type of atom should be provided using [TT]-ei[tt].",
 448         "It should look like:[BR]",
 449         "   [TT]2[tt][BR]",
 450         "   [TT]atomname = nrelectrons[tt][BR]",
 451         "   [TT]atomname = nrelectrons[tt][BR]",
 452         "The first line contains the number of lines to read from the file.",
 453         "There should be one line for each unique atom name in your system.",
 454         "The number of electrons for each atom is modified by its atomic",
 455         "partial charge."
 456     };
 457
 458     output_env_t       oenv;
 459     static const char *dens_opt[] =
 460     { NULL, "mass", "number", "charge", "electron", NULL };
 461     static int         axis        = 2;  /* normal to memb. default z  */
 462     static const char *axtitle     = "Z";
 463     static int         nslices     = 50; /* nr of slices defined       */
 464     static int         ngrps       = 1;  /* nr. of groups              */
 465     static gmx_bool    bSymmetrize = FALSE;
 466     static gmx_bool    bCenter     = FALSE;
 467     t_pargs            pa[]        = {
 468         { "-d", FALSE, etSTR, {&axtitle},
 469           "Take the normal on the membrane in direction X, Y or Z." },
 470         { "-sl",  FALSE, etINT, {&nslices},
 471           "Divide the box in this number of slices." },
 472         { "-dens",    FALSE, etENUM, {dens_opt},
 473           "Density"},
 474         { "-ng",       FALSE, etINT, {&ngrps},
 475           "Number of groups of which to compute densities." },
 476         { "-symm",    FALSE, etBOOL, {&bSymmetrize},
 477           "Symmetrize the density along the axis, with respect to the center. Useful for bilayers." },
 478         { "-center",  FALSE, etBOOL, {&bCenter},
 479           "Shift the center of mass along the axis to zero. This means if your axis is Z and your box is bX, bY, bZ, the center of mass will be at bX/2, bY/2, 0."}
 480     };
 481
 482     const char        *bugs[] = {
 483         "When calculating electron densities, atomnames are used instead of types. This is bad.",
 484     };
 485
 486     double           **density;      /* density per slice          */
 487     real               slWidth;      /* width of one slice         */
 488     char             **grpname;      /* groupnames                 */
 489     int                nr_electrons; /* nr. electrons              */
 490     int               *ngx;          /* sizes of groups            */
 491     t_electron        *el_tab;       /* tabel with nr. of electrons*/
 492     t_topology        *top;          /* topology               */
 493     int                ePBC;
 494     atom_id          **index;        /* indices for all groups     */
 495     int                i;
 496
 497     t_filenm           fnm[] = { /* files for g_density       */
 498         { efTRX, "-f", NULL,  ffREAD },
 499         { efNDX, NULL, NULL,  ffOPTRD },
 500         { efTPX, NULL, NULL,  ffREAD },
 501         { efDAT, "-ei", "electrons", ffOPTRD }, /* file with nr. of electrons */
 502         { efXVG, "-o", "density", ffWRITE },
 503     };
 504
 505 #define NFILE asize(fnm)
 506
 507     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 508                            NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs,
 509                            &oenv))
 510     {
 511         return 0;
 512     }
 513
 514     if (bSymmetrize && !bCenter)
 515     {
 516         fprintf(stderr, "Can not symmetrize without centering. Turning on -center\n");
 517         bCenter = TRUE;
 518     }
 519     /* Calculate axis */
 520     axis = toupper(axtitle[0]) - 'X';
 521
 522     top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC); /* read topology file */
 523     if (dens_opt[0][0] == 'n')
 524     {
 525         for (i = 0; (i < top->atoms.nr); i++)
 526         {
 527             top->atoms.atom[i].m = 1;
 528         }
 529     }
 530     else if (dens_opt[0][0] == 'c')
 531     {
 532         for (i = 0; (i < top->atoms.nr); i++)
 533         {
 534             top->atoms.atom[i].m = top->atoms.atom[i].q;
 535         }
 536     }
 537
 538     snew(grpname, ngrps);
 539     snew(index, ngrps);
 540     snew(ngx, ngrps);
 541
 542     get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, ngx, index, grpname);
 543
 544     if (dens_opt[0][0] == 'e')
 545     {
 546         nr_electrons =  get_electrons(&el_tab, ftp2fn(efDAT, NFILE, fnm));
 547         fprintf(stderr, "Read %d atomtypes from datafile\n", nr_electrons);
 548
 549         calc_electron_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density,
 550                               &nslices, top, ePBC, axis, ngrps, &slWidth, el_tab,
 551                               nr_electrons, bCenter, oenv);
 552     }
 553     else
 554     {
 555         calc_density(ftp2fn(efTRX, NFILE, fnm), index, ngx, &density, &nslices, top,
 556                      ePBC, axis, ngrps, &slWidth, bCenter, oenv);
 557     }
 558
 559     plot_density(density, opt2fn("-o", NFILE, fnm),
 560                  nslices, ngrps, grpname, slWidth, dens_opt,
 561                  bSymmetrize, oenv);
 562
 563     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");  /* view xvgr file */
 564     return 0;
 565 }