src/tools/gmx_densmap.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <math.h>
  43 #include <string.h>
  44
  45 #include "statutil.h"
  46 #include "sysstuff.h"
  47 #include "typedefs.h"
  48 #include "smalloc.h"
  49 #include "macros.h"
  50 #include "vec.h"
  51 #include "pbc.h"
  52 #include "copyrite.h"
  53 #include "futil.h"
  54 #include "statutil.h"
  55 #include "index.h"
  56 #include "mshift.h"
  57 #include "xvgr.h"
  58 #include "princ.h"
  59 #include "rmpbc.h"
  60 #include "txtdump.h"
  61 #include "tpxio.h"
  62 #include "gstat.h"
  63 #include "matio.h"
  64 #include "pbc.h"
  65 #include "gmx_ana.h"
  66
  67
  68 int gmx_densmap(int argc,char *argv[])
  69 {
  70   const char *desc[] = {
  71     "[TT]g_densmap[tt] computes 2D number-density maps.",
  72     "It can make planar and axial-radial density maps.",
  73     "The output [TT].xpm[tt] file can be visualized with for instance xv",
  74     "and can be converted to postscript with [TT]xpm2ps[tt].",
  75     "Optionally, output can be in text form to a [TT].dat[tt] file with [TT]-od[tt], instead of the usual [TT].xpm[tt] file with [TT]-o[tt].",
  76     "[PAR]",
  77     "The default analysis is a 2-D number-density map for a selected",
  78     "group of atoms in the x-y plane.",
  79     "The averaging direction can be changed with the option [TT]-aver[tt].",
  80     "When [TT]-xmin[tt] and/or [TT]-xmax[tt] are set only atoms that are",
  81     "within the limit(s) in the averaging direction are taken into account.",
  82     "The grid spacing is set with the option [TT]-bin[tt].",
  83     "When [TT]-n1[tt] or [TT]-n2[tt] is non-zero, the grid",
  84     "size is set by this option.",
  85     "Box size fluctuations are properly taken into account.",
  86     "[PAR]",
  87     "When options [TT]-amax[tt] and [TT]-rmax[tt] are set, an axial-radial",
  88     "number-density map is made. Three groups should be supplied, the centers",
  89     "of mass of the first two groups define the axis, the third defines the",
  90     "analysis group. The axial direction goes from -amax to +amax, where",
  91     "the center is defined as the midpoint between the centers of mass and",
  92     "the positive direction goes from the first to the second center of mass.",
  93     "The radial direction goes from 0 to rmax or from -rmax to +rmax",
  94     "when the [TT]-mirror[tt] option has been set.",
  95     "[PAR]",
  96     "The normalization of the output is set with the [TT]-unit[tt] option.",
  97     "The default produces a true number density. Unit [TT]nm-2[tt] leaves out",
  98     "the normalization for the averaging or the angular direction.",
  99     "Option [TT]count[tt] produces the count for each grid cell.",
 100     "When you do not want the scale in the output to go",
 101     "from zero to the maximum density, you can set the maximum",
 102     "with the option [TT]-dmax[tt]."
 103   };
 104   static int n1=0,n2=0;
 105   static real xmin=-1,xmax=-1,bin=0.02,dmin=0,dmax=0,amax=0,rmax=0;
 106   static gmx_bool bMirror=FALSE, bSums=FALSE;
 107   static const char *eaver[]={ NULL, "z", "y", "x", NULL };
 108   static const char *eunit[]={ NULL, "nm-3", "nm-2", "count", NULL };
 109
 110   t_pargs pa[] = {
 111     { "-bin", FALSE, etREAL, {&bin},
 112       "Grid size (nm)" },
 113     { "-aver", FALSE, etENUM, {eaver},
 114       "The direction to average over" },
 115     { "-xmin", FALSE, etREAL, {&xmin},
 116       "Minimum coordinate for averaging" },
 117     { "-xmax", FALSE, etREAL, {&xmax},
 118       "Maximum coordinate for averaging" },
 119     { "-n1", FALSE, etINT, {&n1},
 120       "Number of grid cells in the first direction" },
 121     { "-n2", FALSE, etINT, {&n2},
 122       "Number of grid cells in the second direction" },
 123     { "-amax", FALSE, etREAL, {&amax},
 124       "Maximum axial distance from the center"},
 125     { "-rmax", FALSE, etREAL, {&rmax},
 126       "Maximum radial distance" },
 127     { "-mirror", FALSE, etBOOL, {&bMirror},
 128       "Add the mirror image below the axial axis" },
 129     { "-sums", FALSE, etBOOL, {&bSums},
 130       "Print density sums (1D map) to stdout" },
 131     { "-unit", FALSE, etENUM, {eunit},
 132       "Unit for the output" },
 133     { "-dmin", FALSE, etREAL, {&dmin},
 134       "Minimum density in output"},
 135     { "-dmax", FALSE, etREAL, {&dmax},
 136       "Maximum density in output (0 means calculate it)"},
 137   };
 138   gmx_bool       bXmin,bXmax,bRadial;
 139   FILE       *fp;
 140   t_trxstatus *status;
 141   t_topology top;
 142   int        ePBC=-1;
 143   rvec       *x,xcom[2],direction,center,dx;
 144   matrix     box;
 145   real       t,m,mtot;
 146   t_pbc      pbc;
 147   int        cav=0,c1=0,c2=0,natoms;
 148   char       **grpname,title[256],buf[STRLEN];
 149   const char *unit;
 150   int        i,j,k,l,ngrps,anagrp,*gnx=NULL,nindex,nradial=0,nfr,nmpower;
 151   atom_id    **ind=NULL,*index;
 152   real       **grid,maxgrid,m1,m2,box1,box2,*tickx,*tickz,invcellvol;
 153   real       invspa=0,invspz=0,axial,r,vol_old,vol,rowsum;
 154   int        nlev=51;
 155   t_rgb rlo={1,1,1}, rhi={0,0,0};
 156   output_env_t oenv;
 157   const char *label[]={ "x (nm)", "y (nm)", "z (nm)" };
 158   t_filenm fnm[] = {
 159     { efTRX, "-f",   NULL,       ffREAD },
 160     { efTPS, NULL,   NULL,       ffOPTRD },
 161     { efNDX, NULL,   NULL,       ffOPTRD },
 162     { efDAT, "-od",  "densmap",   ffOPTWR },
 163     { efXPM, "-o",   "densmap",   ffWRITE }
 164   };
 165 #define NFILE asize(fnm)
 166   int     npargs;
 167
 168   CopyRight(stderr,argv[0]);
 169   npargs = asize(pa);
 170
 171   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
 172                     NFILE,fnm,npargs,pa,asize(desc),desc,0,NULL,&oenv);
 173
 174   bXmin = opt2parg_bSet("-xmin",npargs,pa);
 175   bXmax = opt2parg_bSet("-xmax",npargs,pa);
 176   bRadial = (amax>0 || rmax>0);
 177   if (bRadial) {
 178     if (amax<=0 || rmax<=0)
 179       gmx_fatal(FARGS,"Both amax and rmax should be larger than zero");
 180   }
 181
 182   if (strcmp(eunit[0],"nm-3") == 0) {
 183     nmpower = -3;
 184     unit = "(nm^-3)";
 185   } else if (strcmp(eunit[0],"nm-2") == 0) {
 186     nmpower = -2;
 187     unit = "(nm^-2)";
 188   } else {
 189     nmpower = 0;
 190     unit = "count";
 191   }
 192
 193   if (ftp2bSet(efTPS,NFILE,fnm) || !ftp2bSet(efNDX,NFILE,fnm))
 194     read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&x,NULL,box,
 195                   bRadial);
 196   if (!bRadial) {
 197     ngrps = 1;
 198     fprintf(stderr,"\nSelect an analysis group\n");
 199   } else {
 200     ngrps = 3;
 201     fprintf(stderr,
 202             "\nSelect two groups to define the axis and an analysis group\n");
 203   }
 204   snew(gnx,ngrps);
 205   snew(grpname,ngrps);
 206   snew(ind,ngrps);
 207   get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),ngrps,gnx,ind,grpname);
 208   anagrp = ngrps - 1;
 209   nindex = gnx[anagrp];
 210   index = ind[anagrp];
 211   if (bRadial) {
 212     if ((gnx[0]>1 || gnx[1]>1) && !ftp2bSet(efTPS,NFILE,fnm))
 213       gmx_fatal(FARGS,"No run input file was supplied (option -s), this is required for the center of mass calculation");
 214   }
 215
 216   switch (eaver[0][0]) {
 217   case 'x': cav = XX; c1 = YY; c2 = ZZ; break;
 218   case 'y': cav = YY; c1 = XX; c2 = ZZ; break;
 219   case 'z': cav = ZZ; c1 = XX; c2 = YY; break;
 220   }
 221
 222   natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 223
 224   if (!bRadial) {
 225     if (n1 == 0)
 226       n1 = (int)(box[c1][c1]/bin + 0.5);
 227     if (n2 == 0)
 228       n2 = (int)(box[c2][c2]/bin + 0.5);
 229   } else {
 230     n1 = (int)(2*amax/bin + 0.5);
 231     nradial = (int)(rmax/bin + 0.5);
 232     invspa = n1/(2*amax);
 233     invspz = nradial/rmax;
 234     if (bMirror)
 235       n2 = 2*nradial;
 236     else
 237       n2 = nradial;
 238   }
 239
 240   snew(grid,n1);
 241   for(i=0; i<n1; i++)
 242     snew(grid[i],n2);
 243
 244   box1 = 0;
 245   box2 = 0;
 246   nfr = 0;
 247   do {
 248     if (!bRadial) {
 249       box1 += box[c1][c1];
 250       box2 += box[c2][c2];
 251       invcellvol = n1*n2;
 252       if (nmpower == -3)
 253         invcellvol /= det(box);
 254       else if (nmpower == -2)
 255         invcellvol /= box[c1][c1]*box[c2][c2];
 256       for(i=0; i<nindex; i++) {
 257         j = index[i];
 258         if ((!bXmin || x[j][cav] >= xmin) &&
 259             (!bXmax || x[j][cav] <= xmax)) {
 260           m1 = x[j][c1]/box[c1][c1];
 261           if (m1 >= 1)
 262             m1 -= 1;
 263           if (m1 < 0)
 264             m1 += 1;
 265           m2 = x[j][c2]/box[c2][c2];
 266           if (m2 >= 1)
 267             m2 -= 1;
 268           if (m2 < 0)
 269             m2 += 1;
 270           grid[(int)(m1*n1)][(int)(m2*n2)] += invcellvol;
 271         }
 272       }
 273     } else {
 274       set_pbc(&pbc,ePBC,box);
 275       for(i=0; i<2; i++) {
 276         if (gnx[i] == 1) {
 277           /* One atom, just copy the coordinates */
 278           copy_rvec(x[ind[i][0]],xcom[i]);
 279         } else {
 280           /* Calculate the center of mass */
 281           clear_rvec(xcom[i]);
 282           mtot = 0;
 283           for(j=0; j<gnx[i]; j++) {
 284             k = ind[i][j];
 285             m = top.atoms.atom[k].m;
 286             for(l=0; l<DIM; l++)
 287               xcom[i][l] += m*x[k][l];
 288             mtot += m;
 289           }
 290           svmul(1/mtot,xcom[i],xcom[i]);
 291         }
 292       }
 293       pbc_dx(&pbc,xcom[1],xcom[0],direction);
 294       for(i=0; i<DIM; i++)
 295         center[i] = xcom[0][i] + 0.5*direction[i];
 296       unitv(direction,direction);
 297       for(i=0; i<nindex; i++) {
 298         j = index[i];
 299         pbc_dx(&pbc,x[j],center,dx);
 300         axial = iprod(dx,direction);
 301         r = sqrt(norm2(dx) - axial*axial);
 302         if (axial>=-amax && axial<amax && r<rmax) {
 303           if (bMirror)
 304             r += rmax;
 305           grid[(int)((axial + amax)*invspa)][(int)(r*invspz)] += 1;
 306         }
 307       }
 308     }
 309     nfr++;
 310   } while(read_next_x(oenv,status,&t,natoms,x,box));
 311   close_trj(status);
 312
 313   /* normalize gridpoints */
 314   maxgrid = 0;
 315   if (!bRadial) {
 316     for (i=0; i<n1; i++) {
 317       for (j=0; j<n2; j++) {
 318         grid[i][j] /= nfr;
 319         if (grid[i][j] > maxgrid)
 320           maxgrid = grid[i][j];
 321       }
 322     }
 323   } else {
 324     for (i=0; i<n1; i++) {
 325       vol_old = 0;
 326       for (j=0; j<nradial; j++) {
 327         switch (nmpower) {
 328         case -3: vol = M_PI*(j+1)*(j+1)/(invspz*invspz*invspa); break;
 329         case -2: vol =            (j+1)/(invspz*invspa);        break;
 330         default: vol =             j+1;                         break;
 331         }
 332         if (bMirror)
 333           k = j + nradial;
 334         else
 335           k = j;
 336         grid[i][k] /= nfr*(vol - vol_old);
 337         if (bMirror)
 338           grid[i][nradial-1-j] = grid[i][k];
 339         vol_old = vol;
 340         if (grid[i][k] > maxgrid)
 341           maxgrid = grid[i][k];
 342       }
 343     }
 344   }
 345   fprintf(stdout,"\n  The maximum density is %f %s\n",maxgrid,unit);
 346   if (dmax > 0)
 347     maxgrid = dmax;
 348
 349   snew(tickx,n1+1);
 350   snew(tickz,n2+1);
 351   if (!bRadial) {
 352     /* normalize box-axes */
 353     box1 /= nfr;
 354     box2 /= nfr;
 355     for (i=0; i<=n1; i++)
 356       tickx[i] = i*box1/n1;
 357     for (i=0; i<=n2; i++)
 358       tickz[i] = i*box2/n2;
 359   } else {
 360     for (i=0; i<=n1; i++)
 361       tickx[i] = i/invspa - amax;
 362     if (bMirror) {
 363       for (i=0; i<=n2; i++)
 364         tickz[i] = i/invspz - rmax;
 365     } else {
 366       for (i=0; i<=n2; i++)
 367         tickz[i] = i/invspz;
 368     }
 369   }
 370
 371   if (bSums)
 372   {
 373     for (i=0;i<n1;++i)
 374     {
 375         fprintf(stdout,"Density sums:\n");
 376     rowsum=0;
 377     for (j=0;j<n2;++j)
 378           rowsum+=grid[i][j];
 379         fprintf(stdout,"%g\t",rowsum);
 380     }
 381         fprintf(stdout,"\n");
 382   }
 383
 384   sprintf(buf,"%s number density",grpname[anagrp]);
 385   if (!bRadial && (bXmin || bXmax)) {
 386     if (!bXmax)
 387       sprintf(buf+strlen(buf),", %c > %g nm",eaver[0][0],xmin);
 388     else if (!bXmin)
 389       sprintf(buf+strlen(buf),", %c < %g nm",eaver[0][0],xmax);
 390     else
 391       sprintf(buf+strlen(buf),", %c: %g - %g nm",eaver[0][0],xmin,xmax);
 392   }
 393   if (ftp2bSet(efDAT,NFILE,fnm))
 394   {
 395   fp = ffopen(ftp2fn(efDAT,NFILE,fnm),"w");
 396   /*optional text form output:  first row is tickz; first col is tickx */
 397   fprintf(fp,"0\t");
 398   for(j=0;j<n2;++j)
 399         fprintf(fp,"%g\t",tickz[j]);
 400   fprintf(fp,"\n");
 401
 402   for (i=0;i<n1;++i)
 403   {
 404     fprintf(fp,"%g\t",tickx[i]);
 405     for (j=0;j<n2;++j)
 406           fprintf(fp,"%g\t",grid[i][j]);
 407         fprintf(fp,"\n");
 408   }
 409   ffclose(fp);
 410   }
 411   else
 412   {
 413   fp = ffopen(ftp2fn(efXPM,NFILE,fnm),"w");
 414   write_xpm(fp,MAT_SPATIAL_X | MAT_SPATIAL_Y,buf,unit,
 415             bRadial ? "axial (nm)" : label[c1],bRadial ? "r (nm)" : label[c2],
 416             n1,n2,tickx,tickz,grid,dmin,maxgrid,rlo,rhi,&nlev);
 417   ffclose(fp);
 418   }
 419
 420   thanx(stderr);
 421
 422   do_view(oenv,opt2fn("-o",NFILE,fnm),NULL);
 423
 424   return 0;
 425 }