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