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