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