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