src/gromacs/gmxana/gmx_density.c

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