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