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