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;
 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     { "-unit", FALSE, etENUM, {eunit},
 125       "Unit for the output" },
 126     { "-dmin", FALSE, etREAL, {&dmin},
 127       "Minimum density in output"},
 128     { "-dmax", FALSE, etREAL, {&dmax},
 129       "Maximum density in output (0 means calculate it)"},
 130   };
 131   bool       bXmin,bXmax,bRadial;
 132   FILE       *fp;
 133   int        status;
 134   t_topology top;
 135   int        ePBC=-1;
 136   rvec       *x,xcom[2],direction,center,dx;
 137   matrix     box;
 138   real       t,m,mtot;
 139   t_pbc      pbc;
 140   int        cav=0,c1=0,c2=0,natoms;
 141   char       **grpname,title[256],buf[STRLEN];
 142   const char *unit;
 143   int        i,j,k,l,ngrps,anagrp,*gnx=NULL,nindex,nradial=0,nfr,nmpower;
 144   atom_id    **ind=NULL,*index;
 145   real       **grid,maxgrid,m1,m2,box1,box2,*tickx,*tickz,invcellvol;
 146   real       invspa=0,invspz=0,axial,r,vol_old,vol;
 147   int        nlev=51;
 148   t_rgb rlo={1,1,1}, rhi={0,0,0};
 149   const char *label[]={ "x (nm)", "y (nm)", "z (nm)" };
 150   t_filenm fnm[] = {
 151     { efTRX, "-f",   NULL,       ffREAD },
 152     { efTPS, NULL,   NULL,       ffOPTRD },
 153     { efNDX, NULL,   NULL,       ffOPTRD },
 154     { efDAT, "-od",  "densmap",   ffOPTWR },
 155     { efXPM, "-o",   "densmap",   ffWRITE }
 156   };
 157 #define NFILE asize(fnm)
 158   int     npargs;
 159
 160   CopyRight(stderr,argv[0]);
 161   npargs = asize(pa);
 162
 163   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
 164                     NFILE,fnm,npargs,pa,asize(desc),desc,0,NULL);
 165
 166   bXmin = opt2parg_bSet("-xmin",npargs,pa);
 167   bXmax = opt2parg_bSet("-xmax",npargs,pa);
 168   bRadial = (amax>0 || rmax>0);
 169   if (bRadial) {
 170     if (amax<=0 || rmax<=0)
 171       gmx_fatal(FARGS,"Both amax and rmax should be larger than zero");
 172   }
 173
 174   if (strcmp(eunit[0],"nm-3") == 0) {
 175     nmpower = -3;
 176     unit = "(nm^-3)";
 177   } else if (strcmp(eunit[0],"nm-2") == 0) {
 178     nmpower = -2;
 179     unit = "(nm^-2)";
 180   } else {
 181     nmpower = 0;
 182     unit = "count";
 183   }
 184
 185   if (ftp2bSet(efTPS,NFILE,fnm) || !ftp2bSet(efNDX,NFILE,fnm))
 186     read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&x,NULL,box,
 187                   bRadial);
 188   if (!bRadial) {
 189     ngrps = 1;
 190     fprintf(stderr,"\nSelect an analysis group\n");
 191   } else {
 192     ngrps = 3;
 193     fprintf(stderr,
 194             "\nSelect two groups to define the axis and an analysis group\n");
 195   }
 196   snew(gnx,ngrps);
 197   snew(grpname,ngrps);
 198   snew(ind,ngrps);
 199   get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),ngrps,gnx,ind,grpname);
 200   anagrp = ngrps - 1;
 201   nindex = gnx[anagrp];
 202   index = ind[anagrp];
 203   if (bRadial) {
 204     if ((gnx[0]>1 || gnx[1]>1) && !ftp2bSet(efTPS,NFILE,fnm))
 205       gmx_fatal(FARGS,"No run input file was supplied (option -s), this is required for the center of mass calculation");
 206   }
 207
 208   switch (eaver[0][0]) {
 209   case 'x': cav = XX; c1 = YY; c2 = ZZ; break;
 210   case 'y': cav = YY; c1 = XX; c2 = ZZ; break;
 211   case 'z': cav = ZZ; c1 = XX; c2 = YY; break;
 212   }
 213
 214   natoms=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 215
 216   if (!bRadial) {
 217     if (n1 == 0)
 218       n1 = (int)(box[c1][c1]/bin + 0.5);
 219     if (n2 == 0)
 220       n2 = (int)(box[c2][c2]/bin + 0.5);
 221   } else {
 222     n1 = (int)(2*amax/bin + 0.5);
 223     nradial = (int)(rmax/bin + 0.5);
 224     invspa = n1/(2*amax);
 225     invspz = nradial/rmax;
 226     if (bMirror)
 227       n2 = 2*nradial;
 228     else
 229       n2 = nradial;
 230   }
 231
 232   snew(grid,n1);
 233   for(i=0; i<n1; i++)
 234     snew(grid[i],n2);
 235
 236   box1 = 0;
 237   box2 = 0;
 238   nfr = 0;
 239   do {
 240     if (!bRadial) {
 241       box1 += box[c1][c1];
 242       box2 += box[c2][c2];
 243       invcellvol = n1*n2;
 244       if (nmpower == -3)
 245         invcellvol /= det(box);
 246       else if (nmpower == -2)
 247         invcellvol /= box[c1][c1]*box[c2][c2];
 248       for(i=0; i<nindex; i++) {
 249         j = index[i];
 250         if ((!bXmin || x[j][cav] >= xmin) &&
 251             (!bXmax || x[j][cav] <= xmax)) {
 252           m1 = x[j][c1]/box[c1][c1];
 253           if (m1 >= 1)
 254             m1 -= 1;
 255           if (m1 < 0)
 256             m1 += 1;
 257           m2 = x[j][c2]/box[c2][c2];
 258           if (m2 >= 1)
 259             m2 -= 1;
 260           if (m2 < 0)
 261             m2 += 1;
 262           grid[(int)(m1*n1)][(int)(m2*n2)] += invcellvol;
 263         }
 264       }
 265     } else {
 266       set_pbc(&pbc,ePBC,box);
 267       for(i=0; i<2; i++) {
 268         if (gnx[i] == 1) {
 269           /* One atom, just copy the coordinates */
 270           copy_rvec(x[ind[i][0]],xcom[i]);
 271         } else {
 272           /* Calculate the center of mass */
 273           clear_rvec(xcom[i]);
 274           mtot = 0;
 275           for(j=0; j<gnx[i]; j++) {
 276             k = ind[i][j];
 277             m = top.atoms.atom[k].m;
 278             for(l=0; l<DIM; l++)
 279               xcom[i][l] += m*x[k][l];
 280             mtot += m;
 281           }
 282           svmul(1/mtot,xcom[i],xcom[i]);
 283         }
 284       }
 285       pbc_dx(&pbc,xcom[1],xcom[0],direction);
 286       for(i=0; i<DIM; i++)
 287         center[i] = xcom[0][i] + 0.5*direction[i];
 288       unitv(direction,direction);
 289       for(i=0; i<nindex; i++) {
 290         j = index[i];
 291         pbc_dx(&pbc,x[j],center,dx);
 292         axial = iprod(dx,direction);
 293         r = sqrt(norm2(dx) - axial*axial);
 294         if (axial>=-amax && axial<amax && r<rmax) {
 295           if (bMirror)
 296             r += rmax;
 297           grid[(int)((axial + amax)*invspa)][(int)(r*invspz)] += 1;
 298         }
 299       }
 300     }
 301     nfr++;
 302   } while(read_next_x(status,&t,natoms,x,box));
 303   close_trj(status);
 304
 305   /* normalize gridpoints */
 306   maxgrid = 0;
 307   if (!bRadial) {
 308     for (i=0; i<n1; i++) {
 309       for (j=0; j<n2; j++) {
 310         grid[i][j] /= nfr;
 311         if (grid[i][j] > maxgrid)
 312           maxgrid = grid[i][j];
 313       }
 314     }
 315   } else {
 316     for (i=0; i<n1; i++) {
 317       vol_old = 0;
 318       for (j=0; j<nradial; j++) {
 319         switch (nmpower) {
 320         case -3: vol = M_PI*(j+1)*(j+1)/(invspz*invspz*invspa); break;
 321         case -2: vol =            (j+1)/(invspz*invspa);        break;
 322         default: vol =             j+1;                         break;
 323         }
 324         if (bMirror)
 325           k = j + nradial;
 326         else
 327           k = j;
 328         grid[i][k] /= nfr*(vol - vol_old);
 329         if (bMirror)
 330           grid[i][nradial-1-j] = grid[i][k];
 331         vol_old = vol;
 332         if (grid[i][k] > maxgrid)
 333           maxgrid = grid[i][k];
 334       }
 335     }
 336   }
 337   fprintf(stdout,"\n  The maximum density is %f %s\n",maxgrid,unit);
 338   if (dmax > 0)
 339     maxgrid = dmax;
 340
 341   snew(tickx,n1+1);
 342   snew(tickz,n2+1);
 343   if (!bRadial) {
 344     /* normalize box-axes */
 345     box1 /= nfr;
 346     box2 /= nfr;
 347     for (i=0; i<=n1; i++)
 348       tickx[i] = i*box1/n1;
 349     for (i=0; i<=n2; i++)
 350       tickz[i] = i*box2/n2;
 351   } else {
 352     for (i=0; i<=n1; i++)
 353       tickx[i] = i/invspa - amax;
 354     if (bMirror) {
 355       for (i=0; i<=n2; i++)
 356         tickz[i] = i/invspz - rmax;
 357     } else {
 358       for (i=0; i<=n2; i++)
 359         tickz[i] = i/invspz;
 360     }
 361   }
 362
 363   sprintf(buf,"%s number density",grpname[anagrp]);
 364   if (!bRadial && (bXmin || bXmax)) {
 365     if (!bXmax)
 366       sprintf(buf+strlen(buf),", %c > %g nm",eaver[0][0],xmin);
 367     else if (!bXmin)
 368       sprintf(buf+strlen(buf),", %c < %g nm",eaver[0][0],xmax);
 369     else
 370       sprintf(buf+strlen(buf),", %c: %g - %g nm",eaver[0][0],xmin,xmax);
 371   }
 372   if (ftp2bSet(efDAT,NFILE,fnm))
 373   {
 374   fp = ffopen(ftp2fn(efDAT,NFILE,fnm),"w");
 375   /*optional text form output:  first row is tickz; first col is tickx */
 376   fprintf(fp,"0\t");
 377   for(j=0;j<n2;++j)
 378         fprintf(fp,"%g\t",tickz[j]);
 379   fprintf(fp,"\n");
 380
 381   for (i=0;i<n1;++i)
 382   {
 383     fprintf(fp,"%g\t",tickx[i]);
 384     for (j=0;j<n2;++j)
 385           fprintf(fp,"%g\t",grid[i][j]);
 386         fprintf(fp,"\n");
 387   }
 388   ffclose(fp);
 389   }
 390   else
 391   {
 392   fp = ffopen(ftp2fn(efXPM,NFILE,fnm),"w");
 393   write_xpm(fp,MAT_SPATIAL_X | MAT_SPATIAL_Y,buf,unit,
 394             bRadial ? "axial (nm)" : label[c1],bRadial ? "r (nm)" : label[c2],
 395             n1,n2,tickx,tickz,grid,dmin,maxgrid,rlo,rhi,&nlev);
 396   ffclose(fp);
 397   }
 398
 399   thanx(stderr);
 400
 401   do_view(opt2fn("-o",NFILE,fnm),NULL);
 402
 403   return 0;
 404 }