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