src/tools/gmx_densmap.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
  39 #include <math.h>
  40 #include <string.h>
  41
  42 #include "statutil.h"
  43 #include "sysstuff.h"
  44 #include "typedefs.h"
  45 #include "smalloc.h"
  46 #include "macros.h"
  47 #include "vec.h"
  48 #include "pbc.h"
  49 #include "copyrite.h"
  50 #include "futil.h"
  51 #include "statutil.h"
  52 #include "index.h"
  53 #include "mshift.h"
  54 #include "xvgr.h"
  55 #include "princ.h"
  56 #include "rmpbc.h"
  57 #include "txtdump.h"
  58 #include "tpxio.h"
  59 #include "gstat.h"
  60 #include "matio.h"
  61 #include "pbc.h"
  62
  63 int gmx_densmap(int argc,char *argv[])
  64 {
  65   const char *desc[] = {
  66     "g_densmap computes 2D number-density maps.",
  67     "It can make planar and axial-radial density maps.",
  68     "The output [TT].xpm[tt] file can be visualized with for instance xv",
  69     "and can be converted to postscript with xpm2ps.",
  70         "Optionally, output can be in text form to a .dat file.",
  71     "[PAR]",
  72     "The default analysis is a 2-D number-density map for a selected",
  73     "group of atoms in the x-y plane.",
  74     "The averaging direction can be changed with the option [TT]-aver[tt].",
  75     "When [TT]-xmin[tt] and/or [TT]-xmax[tt] are set only atoms that are",
  76     "within the limit(s) in the averaging direction are taken into account.",
  77     "The grid spacing is set with the option [TT]-bin[tt].",
  78     "When [TT]-n1[tt] or [TT]-n2[tt] is non-zero, the grid",
  79     "size is set by this option.",
  80     "Box size fluctuations are properly taken into account.",
  81     "[PAR]",
  82     "When options [TT]-amax[tt] and [TT]-rmax[tt] are set, an axial-radial",
  83     "number-density map is made. Three groups should be supplied, the centers",
  84     "of mass of the first two groups define the axis, the third defines the",
  85     "analysis group. The axial direction goes from -amax to +amax, where",
  86     "the center is defined as the midpoint between the centers of mass and",
  87     "the positive direction goes from the first to the second center of mass.",
  88     "The radial direction goes from 0 to rmax or from -rmax to +rmax",
  89     "when the [TT]-mirror[tt] option has been set.",
  90     "[PAR]",
  91     "The normalization of the output is set with the [TT]-unit[tt] option.",
  92     "The default produces a true number density. Unit [TT]nm-2[tt] leaves out",
  93     "the normalization for the averaging or the angular direction.",
  94     "Option [TT]count[tt] produces the count for each grid cell.",
  95     "When you do not want the scale in the output to go",
  96     "from zero to the maximum density, you can set the maximum",
  97     "with the option [TT]-dmax[tt]."
  98   };
  99   static int n1=0,n2=0;
 100   static real xmin=-1,xmax=-1,bin=0.02,dmin=0,dmax=0,amax=0,rmax=0;
 101   static bool bMirror=FALSE, bSums=FALSE;
 102   static const char *eaver[]={ NULL, "z", "y", "x", NULL };
 103   static const char *eunit[]={ NULL, "nm-3", "nm-2", "count", NULL };
 104
 105   t_pargs pa[] = {
 106     { "-bin", FALSE, etREAL, {&bin},
 107       "Grid size (nm)" },
 108     { "-aver", FALSE, etENUM, {eaver},
 109       "The direction to average over" },
 110     { "-xmin", FALSE, etREAL, {&xmin},
 111       "Minimum coordinate for averaging" },
 112     { "-xmax", FALSE, etREAL, {&xmax},
 113       "Maximum coordinate for averaging" },
 114     { "-n1", FALSE, etINT, {&n1},
 115       "Number of grid cells in the first direction" },
 116     { "-n2", FALSE, etINT, {&n2},
 117       "Number of grid cells in the second direction" },
 118     { "-amax", FALSE, etREAL, {&amax},
 119       "Maximum axial distance from the center"},
 120     { "-rmax", FALSE, etREAL, {&rmax},
 121       "Maximum radial distance" },
 122     { "-mirror", FALSE, etBOOL, {&bMirror},
 123       "Add the mirror image below the axial axis" },
 124     { "-sums", FALSE, etBOOL, {&bSums},
 125       "Print density sums (1D map) to stdout" },
 126     { "-unit", FALSE, etENUM, {eunit},
 127       "Unit for the output" },
 128     { "-dmin", FALSE, etREAL, {&dmin},
 129       "Minimum density in output"},
 130     { "-dmax", FALSE, etREAL, {&dmax},
 131       "Maximum density in output (0 means calculate it)"},
 132   };
 133   bool       bXmin,bXmax,bRadial;
 134   FILE       *fp;
 135   int        status;
 136   t_topology top;
 137   int        ePBC=-1;
 138   rvec       *x,xcom[2],direction,center,dx;
 139   matrix     box;
 140   real       t,m,mtot;
 141   t_pbc      pbc;
 142   int        cav=0,c1=0,c2=0,natoms;
 143   char       **grpname,title[256],buf[STRLEN];
 144   const char *unit;
 145   int        i,j,k,l,ngrps,anagrp,*gnx=NULL,nindex,nradial=0,nfr,nmpower;
 146   atom_id    **ind=NULL,*index;
 147   real       **grid,maxgrid,m1,m2,box1,box2,*tickx,*tickz,invcellvol;
 148   real       invspa=0,invspz=0,axial,r,vol_old,vol,rowsum;
 149   int        nlev=51;
 150   t_rgb rlo={1,1,1}, rhi={0,0,0};
 151   output_env_t oenv;
 152   const char *label[]={ "x (nm)", "y (nm)", "z (nm)" };
 153   t_filenm fnm[] = {
 154     { efTRX, "-f",   NULL,       ffREAD },
 155     { efTPS, NULL,   NULL,       ffOPTRD },
 156     { efNDX, NULL,   NULL,       ffOPTRD },
 157     { efDAT, "-od",  "densmap",   ffOPTWR },
 158     { efXPM, "-o",   "densmap",   ffWRITE }
 159   };
 160 #define NFILE asize(fnm)
 161   int     npargs;
 162
 163   CopyRight(stderr,argv[0]);
 164   npargs = asize(pa);
 165
 166   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
 167                     NFILE,fnm,npargs,pa,asize(desc),desc,0,NULL,&oenv);
 168
 169   bXmin = opt2parg_bSet("-xmin",npargs,pa);
 170   bXmax = opt2parg_bSet("-xmax",npargs,pa);
 171   bRadial = (amax>0 || rmax>0);
 172   if (bRadial) {
 173     if (amax<=0 || rmax<=0)
 174       gmx_fatal(FARGS,"Both amax and rmax should be larger than zero");
 175   }
 176
 177   if (strcmp(eunit[0],"nm-3") == 0) {
 178     nmpower = -3;
 179     unit = "(nm^-3)";
 180   } else if (strcmp(eunit[0],"nm-2") == 0) {
 181     nmpower = -2;
 182     unit = "(nm^-2)";
 183   } else {
 184     nmpower = 0;
 185     unit = "count";
 186   }
 187
 188   if (ftp2bSet(efTPS,NFILE,fnm) || !ftp2bSet(efNDX,NFILE,fnm))
 189     read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&x,NULL,box,
 190                   bRadial);
 191   if (!bRadial) {
 192     ngrps = 1;
 193     fprintf(stderr,"\nSelect an analysis group\n");
 194   } else {
 195     ngrps = 3;
 196     fprintf(stderr,
 197             "\nSelect two groups to define the axis and an analysis group\n");
 198   }
 199   snew(gnx,ngrps);
 200   snew(grpname,ngrps);
 201   snew(ind,ngrps);
 202   get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),ngrps,gnx,ind,grpname);
 203   anagrp = ngrps - 1;
 204   nindex = gnx[anagrp];
 205   index = ind[anagrp];
 206   if (bRadial) {
 207     if ((gnx[0]>1 || gnx[1]>1) && !ftp2bSet(efTPS,NFILE,fnm))
 208       gmx_fatal(FARGS,"No run input file was supplied (option -s), this is required for the center of mass calculation");
 209   }
 210
 211   switch (eaver[0][0]) {
 212   case 'x': cav = XX; c1 = YY; c2 = ZZ; break;
 213   case 'y': cav = YY; c1 = XX; c2 = ZZ; break;
 214   case 'z': cav = ZZ; c1 = XX; c2 = YY; break;
 215   }
 216
 217   natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 218
 219   if (!bRadial) {
 220     if (n1 == 0)
 221       n1 = (int)(box[c1][c1]/bin + 0.5);
 222     if (n2 == 0)
 223       n2 = (int)(box[c2][c2]/bin + 0.5);
 224   } else {
 225     n1 = (int)(2*amax/bin + 0.5);
 226     nradial = (int)(rmax/bin + 0.5);
 227     invspa = n1/(2*amax);
 228     invspz = nradial/rmax;
 229     if (bMirror)
 230       n2 = 2*nradial;
 231     else
 232       n2 = nradial;
 233   }
 234
 235   snew(grid,n1);
 236   for(i=0; i<n1; i++)
 237     snew(grid[i],n2);
 238
 239   box1 = 0;
 240   box2 = 0;
 241   nfr = 0;
 242   do {
 243     if (!bRadial) {
 244       box1 += box[c1][c1];
 245       box2 += box[c2][c2];
 246       invcellvol = n1*n2;
 247       if (nmpower == -3)
 248         invcellvol /= det(box);
 249       else if (nmpower == -2)
 250         invcellvol /= box[c1][c1]*box[c2][c2];
 251       for(i=0; i<nindex; i++) {
 252         j = index[i];
 253         if ((!bXmin || x[j][cav] >= xmin) &&
 254             (!bXmax || x[j][cav] <= xmax)) {
 255           m1 = x[j][c1]/box[c1][c1];
 256           if (m1 >= 1)
 257             m1 -= 1;
 258           if (m1 < 0)
 259             m1 += 1;
 260           m2 = x[j][c2]/box[c2][c2];
 261           if (m2 >= 1)
 262             m2 -= 1;
 263           if (m2 < 0)
 264             m2 += 1;
 265           grid[(int)(m1*n1)][(int)(m2*n2)] += invcellvol;
 266         }
 267       }
 268     } else {
 269       set_pbc(&pbc,ePBC,box);
 270       for(i=0; i<2; i++) {
 271         if (gnx[i] == 1) {
 272           /* One atom, just copy the coordinates */
 273           copy_rvec(x[ind[i][0]],xcom[i]);
 274         } else {
 275           /* Calculate the center of mass */
 276           clear_rvec(xcom[i]);
 277           mtot = 0;
 278           for(j=0; j<gnx[i]; j++) {
 279             k = ind[i][j];
 280             m = top.atoms.atom[k].m;
 281             for(l=0; l<DIM; l++)
 282               xcom[i][l] += m*x[k][l];
 283             mtot += m;
 284           }
 285           svmul(1/mtot,xcom[i],xcom[i]);
 286         }
 287       }
 288       pbc_dx(&pbc,xcom[1],xcom[0],direction);
 289       for(i=0; i<DIM; i++)
 290         center[i] = xcom[0][i] + 0.5*direction[i];
 291       unitv(direction,direction);
 292       for(i=0; i<nindex; i++) {
 293         j = index[i];
 294         pbc_dx(&pbc,x[j],center,dx);
 295         axial = iprod(dx,direction);
 296         r = sqrt(norm2(dx) - axial*axial);
 297         if (axial>=-amax && axial<amax && r<rmax) {
 298           if (bMirror)
 299             r += rmax;
 300           grid[(int)((axial + amax)*invspa)][(int)(r*invspz)] += 1;
 301         }
 302       }
 303     }
 304     nfr++;
 305   } while(read_next_x(oenv,status,&t,natoms,x,box));
 306   close_trj(status);
 307
 308   /* normalize gridpoints */
 309   maxgrid = 0;
 310   if (!bRadial) {
 311     for (i=0; i<n1; i++) {
 312       for (j=0; j<n2; j++) {
 313         grid[i][j] /= nfr;
 314         if (grid[i][j] > maxgrid)
 315           maxgrid = grid[i][j];
 316       }
 317     }
 318   } else {
 319     for (i=0; i<n1; i++) {
 320       vol_old = 0;
 321       for (j=0; j<nradial; j++) {
 322         switch (nmpower) {
 323         case -3: vol = M_PI*(j+1)*(j+1)/(invspz*invspz*invspa); break;
 324         case -2: vol =            (j+1)/(invspz*invspa);        break;
 325         default: vol =             j+1;                         break;
 326         }
 327         if (bMirror)
 328           k = j + nradial;
 329         else
 330           k = j;
 331         grid[i][k] /= nfr*(vol - vol_old);
 332         if (bMirror)
 333           grid[i][nradial-1-j] = grid[i][k];
 334         vol_old = vol;
 335         if (grid[i][k] > maxgrid)
 336           maxgrid = grid[i][k];
 337       }
 338     }
 339   }
 340   fprintf(stdout,"\n  The maximum density is %f %s\n",maxgrid,unit);
 341   if (dmax > 0)
 342     maxgrid = dmax;
 343
 344   snew(tickx,n1+1);
 345   snew(tickz,n2+1);
 346   if (!bRadial) {
 347     /* normalize box-axes */
 348     box1 /= nfr;
 349     box2 /= nfr;
 350     for (i=0; i<=n1; i++)
 351       tickx[i] = i*box1/n1;
 352     for (i=0; i<=n2; i++)
 353       tickz[i] = i*box2/n2;
 354   } else {
 355     for (i=0; i<=n1; i++)
 356       tickx[i] = i/invspa - amax;
 357     if (bMirror) {
 358       for (i=0; i<=n2; i++)
 359         tickz[i] = i/invspz - rmax;
 360     } else {
 361       for (i=0; i<=n2; i++)
 362         tickz[i] = i/invspz;
 363     }
 364   }
 365
 366   if (bSums)
 367   {
 368     for (i=0;i<n1;++i)
 369     {
 370         fprintf(stdout,"Density sums:\n");
 371     rowsum=0;
 372     for (j=0;j<n2;++j)
 373           rowsum+=grid[i][j];
 374         fprintf(stdout,"%g\t",rowsum);
 375     }
 376         fprintf(stdout,"\n");
 377   }
 378
 379   sprintf(buf,"%s number density",grpname[anagrp]);
 380   if (!bRadial && (bXmin || bXmax)) {
 381     if (!bXmax)
 382       sprintf(buf+strlen(buf),", %c > %g nm",eaver[0][0],xmin);
 383     else if (!bXmin)
 384       sprintf(buf+strlen(buf),", %c < %g nm",eaver[0][0],xmax);
 385     else
 386       sprintf(buf+strlen(buf),", %c: %g - %g nm",eaver[0][0],xmin,xmax);
 387   }
 388   if (ftp2bSet(efDAT,NFILE,fnm))
 389   {
 390   fp = ffopen(ftp2fn(efDAT,NFILE,fnm),"w");
 391   /*optional text form output:  first row is tickz; first col is tickx */
 392   fprintf(fp,"0\t");
 393   for(j=0;j<n2;++j)
 394         fprintf(fp,"%g\t",tickz[j]);
 395   fprintf(fp,"\n");
 396
 397   for (i=0;i<n1;++i)
 398   {
 399     fprintf(fp,"%g\t",tickx[i]);
 400     for (j=0;j<n2;++j)
 401           fprintf(fp,"%g\t",grid[i][j]);
 402         fprintf(fp,"\n");
 403   }
 404   ffclose(fp);
 405   }
 406   else
 407   {
 408   fp = ffopen(ftp2fn(efXPM,NFILE,fnm),"w");
 409   write_xpm(fp,MAT_SPATIAL_X | MAT_SPATIAL_Y,buf,unit,
 410             bRadial ? "axial (nm)" : label[c1],bRadial ? "r (nm)" : label[c2],
 411             n1,n2,tickx,tickz,grid,dmin,maxgrid,rlo,rhi,&nlev);
 412   ffclose(fp);
 413   }
 414
 415   thanx(stderr);
 416
 417   do_view(oenv,opt2fn("-o",NFILE,fnm),NULL);
 418
 419   return 0;
 420 }