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