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