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