src/tools/gmx_spatial.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.0
  10  *
  11  * Copyright (c) 1991-2001
  12  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  13  * University of Groningen, The Netherlands
  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  * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
  31  *
  32  * And Hey:
  33  * Gyas ROwers Mature At Cryogenic Speed
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39
  40
  41 #include "statutil.h"
  42 #include "typedefs.h"
  43 #include "smalloc.h"
  44 #include "vec.h"
  45 #include "copyrite.h"
  46 #include "statutil.h"
  47 #include "tpxio.h"
  48 #include "math.h"
  49 #include "index.h"
  50 #include "pbc.h"
  51 #include "rmpbc.h"
  52 #include "gmx_ana.h"
  53
  54
  55 static const double bohr=0.529177249;  /* conversion factor to compensate for VMD plugin conversion... */
  56
  57 static void mequit(void){
  58   printf("Memory allocation error\n");
  59   exit(1);
  60 }
  61
  62 int gmx_spatial(int argc,char *argv[])
  63 {
  64   const char *desc[] = {
  65     "g_spatial calculates the spatial distribution function and ",
  66     "outputs it in a form that can be read by VMD as Gaussian98 cube format. ",
  67     "This was developed from template.c (gromacs-3.3). ",
  68     "For a system of 32K atoms and a 50ns trajectory, the SDF can be generated ",
  69     "in about 30 minutes, with most of the time dedicated to the two runs through ",
  70     "trjconv that are required to center everything properly. ",
  71     "This also takes a whole bunch of space (3 copies of the xtc file). ",
  72     "Still, the pictures are pretty and very informative when the fitted selection is properly made. ",
  73     "3-4 atoms in a widely mobile group like a free amino acid in solution works ",
  74     "well, or select the protein backbone in a stable folded structure to get the SDF ",
  75     "of solvent and look at the time-averaged solvation shell. ",
  76     "It is also possible using this program to generate the SDF based on some arbitrarty ",
  77     "Cartesian coordinate. To do that, simply omit the preliminary trjconv steps. \n",
  78     "USAGE: \n",
  79     "1. Use make_ndx to create a group containing the atoms around which you want the SDF \n",
  80     "2. trjconv -s a.tpr -f a.xtc -o b.xtc -center tric -ur compact -pbc none \n",
  81     "3. trjconv -s a.tpr -f b.xtc -o c.xtc -fit rot+trans \n",
  82     "4. run g_spatial on the xtc output of step #3. \n",
  83     "5. Load grid.cube into VMD and view as an isosurface. \n",
  84     "*** Systems such as micelles will require trjconv -pbc cluster between steps 1 and 2\n",
  85     "WARNINGS: \n",
  86     "The SDF will be generated for a cube that contains all bins that have some non-zero occupancy. ",
  87     "However, the preparatory -fit rot+trans option to trjconv implies that your system will be rotating ",
  88     "and translating in space (in order that the selected group does not). Therefore the values that are ",
  89     "returned will only be valid for some region around your central group/coordinate that has full overlap ",
  90     "with system volume throughout the entire translated/rotated system over the course of the trajectory. ",
  91     "It is up to the user to ensure that this is the case. \n",
  92     "BUGS: \n",
  93     "When the allocated memory is not large enough, a segmentation fault may occur. This is usually detected ",
  94     "and the program is halted prior to the fault while displaying a warning message suggesting the use of the -nab ",
  95     "option. However, the program does not detect all such events. If you encounter a segmentation fault, run it again ",
  96     "with an increased -nab value. \n",
  97     "RISKY OPTIONS: \n",
  98     "To reduce the amount of space and time required, you can output only the coords ",
  99     "that are going to be used in the first and subsequent run through trjconv. ",
 100     "However, be sure to set the -nab option to a sufficiently high value since ",
 101     "memory is allocated for cube bins based on the initial coords and the -nab ",
 102     "(Number of Additional Bins) option value. \n"
 103   };
 104
 105   static gmx_bool bPBC=FALSE;
 106   static gmx_bool bSHIFT=FALSE;
 107   static int iIGNOREOUTER=-1; /*Positive values may help if the surface is spikey */
 108   static gmx_bool bCUTDOWN=TRUE;
 109   static real rBINWIDTH=0.05; /* nm */
 110   static gmx_bool bCALCDIV=TRUE;
 111   static int iNAB=4;
 112
 113   t_pargs pa[] = {
 114     { "-pbc",      FALSE, etBOOL, {&bPBC},
 115       "Use periodic boundary conditions for computing distances" },
 116     { "-div",      FALSE, etBOOL, {&bCALCDIV},
 117       "Calculate and apply the divisor for bin occupancies based on atoms/minimal cube size. Set as TRUE for visualization and as FALSE (-nodiv) to get accurate counts per frame" },
 118     { "-ign",      FALSE, etINT, {&iIGNOREOUTER},
 119       "Do not display this number of outer cubes (positive values may reduce boundary speckles; -1 ensures outer surface is visible)" },
 120     /*    { "-cut",      bCUTDOWN, etBOOL, {&bCUTDOWN},*/
 121     /*      "Display a total cube that is of minimal size" }, */
 122     { "-bin",      FALSE, etREAL, {&rBINWIDTH},
 123       "Width of the bins in nm" },
 124     { "-nab",      FALSE, etINT, {&iNAB},
 125       "Number of additional bins to ensure proper memory allocation" }
 126   };
 127
 128   double MINBIN[3];
 129   double MAXBIN[3];
 130   t_topology top;
 131   int        ePBC;
 132   char       title[STRLEN];
 133   t_trxframe fr;
 134   rvec       *xtop,*shx[26];
 135   matrix     box,box_pbc;
 136   t_trxstatus *status;
 137   int        flags = TRX_READ_X;
 138   t_pbc      pbc;
 139   t_atoms    *atoms;
 140   int        natoms;
 141   char        *grpnm,*grpnmp;
 142   atom_id     *index,*indexp;
 143   int         i,nidx,nidxp;
 144   int v;
 145   int j,k;
 146   long ***bin=(long ***)NULL;
 147   long nbin[3];
 148   FILE *flp;
 149   long x,y,z,minx,miny,minz,maxx,maxy,maxz;
 150   long numfr, numcu;
 151   long tot,max,min;
 152   double norm;
 153   output_env_t oenv;
 154   gmx_rmpbc_t  gpbc=NULL;
 155
 156   t_filenm fnm[] = {
 157     { efTPS,  NULL,  NULL, ffREAD },   /* this is for the topology */
 158     { efTRX, "-f", NULL, ffREAD },      /* and this for the trajectory */
 159     { efNDX, NULL, NULL, ffOPTRD }
 160   };
 161
 162 #define NFILE asize(fnm)
 163
 164   CopyRight(stderr,argv[0]);
 165
 166   /* This is the routine responsible for adding default options,
 167    * calling the X/motif interface, etc. */
 168   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW,
 169                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 170
 171   read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&xtop,NULL,box,TRUE);
 172   sfree(xtop);
 173
 174   atoms=&(top.atoms);
 175   printf("Select group to generate SDF:\n");
 176   get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&nidx,&index,&grpnm);
 177   printf("Select group to output coords (e.g. solute):\n");
 178   get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&nidxp,&indexp,&grpnmp);
 179
 180   /* The first time we read data is a little special */
 181   natoms=read_first_frame(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&fr,flags);
 182
 183   /* Memory Allocation */
 184   MINBIN[XX]=MAXBIN[XX]=fr.x[0][XX];
 185   MINBIN[YY]=MAXBIN[YY]=fr.x[0][YY];
 186   MINBIN[ZZ]=MAXBIN[ZZ]=fr.x[0][ZZ];
 187   for(i=1; i<top.atoms.nr; ++i) {
 188     if(fr.x[i][XX]<MINBIN[XX])MINBIN[XX]=fr.x[i][XX];
 189     if(fr.x[i][XX]>MAXBIN[XX])MAXBIN[XX]=fr.x[i][XX];
 190     if(fr.x[i][YY]<MINBIN[YY])MINBIN[YY]=fr.x[i][YY];
 191     if(fr.x[i][YY]>MAXBIN[YY])MAXBIN[YY]=fr.x[i][YY];
 192     if(fr.x[i][ZZ]<MINBIN[ZZ])MINBIN[ZZ]=fr.x[i][ZZ];
 193     if(fr.x[i][ZZ]>MAXBIN[ZZ])MAXBIN[ZZ]=fr.x[i][ZZ];
 194   }
 195   for (i=ZZ; i>=XX; --i){
 196     MAXBIN[i]=(ceil((MAXBIN[i]-MINBIN[i])/rBINWIDTH)+(double)iNAB)*rBINWIDTH+MINBIN[i];
 197     MINBIN[i]-=(double)iNAB*rBINWIDTH;
 198     nbin[i]=(long)ceil((MAXBIN[i]-MINBIN[i])/rBINWIDTH);
 199   }
 200   bin=(long ***)malloc(nbin[XX]*sizeof(long **));
 201   if(!bin)mequit();
 202   for(i=0; i<nbin[XX]; ++i){
 203     bin[i]=(long **)malloc(nbin[YY]*sizeof(long *));
 204     if(!bin[i])mequit();
 205     for(j=0; j<nbin[YY]; ++j){
 206       bin[i][j]=(long *)calloc(nbin[ZZ],sizeof(long));
 207       if(!bin[i][j])mequit();
 208     }
 209   }
 210   copy_mat(box,box_pbc);
 211   numfr=0;
 212   minx=miny=minz=999;
 213   maxx=maxy=maxz=0;
 214
 215   if (bPBC)
 216     gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
 217   /* This is the main loop over frames */
 218   do {
 219     /* Must init pbc every step because of pressure coupling */
 220
 221     copy_mat(box,box_pbc);
 222     if (bPBC) {
 223       gmx_rmpbc_trxfr(gpbc,&fr);
 224       set_pbc(&pbc,ePBC,box_pbc);
 225     }
 226
 227     for(i=0; i<nidx; i++) {
 228       if(fr.x[index[i]][XX]<MINBIN[XX]||fr.x[index[i]][XX]>MAXBIN[XX]||
 229          fr.x[index[i]][YY]<MINBIN[YY]||fr.x[index[i]][YY]>MAXBIN[YY]||
 230          fr.x[index[i]][ZZ]<MINBIN[ZZ]||fr.x[index[i]][ZZ]>MAXBIN[ZZ])
 231         {
 232           printf("There was an item outside of the allocated memory. Increase the value given with the -nab option.\n");
 233           printf("Memory was allocated for [%f,%f,%f]\tto\t[%f,%f,%f]\n",MINBIN[XX],MINBIN[YY],MINBIN[ZZ],MAXBIN[XX],MAXBIN[YY],MAXBIN[ZZ]);
 234           printf("Memory was required for [%f,%f,%f]\n",fr.x[index[i]][XX],fr.x[index[i]][YY],fr.x[index[i]][ZZ]);
 235           exit(1);
 236         }
 237       x=(long)ceil((fr.x[index[i]][XX]-MINBIN[XX])/rBINWIDTH);
 238       y=(long)ceil((fr.x[index[i]][YY]-MINBIN[YY])/rBINWIDTH);
 239       z=(long)ceil((fr.x[index[i]][ZZ]-MINBIN[ZZ])/rBINWIDTH);
 240       ++bin[x][y][z];
 241       if(x<minx)minx=x;
 242       if(x>maxx)maxx=x;
 243       if(y<miny)miny=y;
 244       if(y>maxy)maxy=y;
 245       if(z<minz)minz=z;
 246       if(z>maxz)maxz=z;
 247     }
 248     numfr++;
 249     /* printf("%f\t%f\t%f\n",box[XX][XX],box[YY][YY],box[ZZ][ZZ]); */
 250
 251   } while(read_next_frame(oenv,status,&fr));
 252
 253   if (bPBC)
 254     gmx_rmpbc_done(gpbc);
 255
 256   if(!bCUTDOWN){
 257     minx=miny=minz=0;
 258     maxx=nbin[XX];
 259     maxy=nbin[YY];
 260     maxz=nbin[ZZ];
 261   }
 262
 263   /* OUTPUT */
 264   flp=ffopen("grid.cube","w");
 265   fprintf(flp,"Spatial Distribution Function\n");
 266   fprintf(flp,"test\n");
 267   fprintf(flp,"%5d%12.6f%12.6f%12.6f\n",nidxp,(MINBIN[XX]+(minx+iIGNOREOUTER)*rBINWIDTH)*10./bohr,(MINBIN[YY]+(miny+iIGNOREOUTER)*rBINWIDTH)*10./bohr,(MINBIN[ZZ]+(minz+iIGNOREOUTER)*rBINWIDTH)*10./bohr);
 268   fprintf(flp,"%5ld%12.6f%12.6f%12.6f\n",maxx-minx+1-(2*iIGNOREOUTER),rBINWIDTH*10./bohr,0.,0.);
 269   fprintf(flp,"%5ld%12.6f%12.6f%12.6f\n",maxy-miny+1-(2*iIGNOREOUTER),0.,rBINWIDTH*10./bohr,0.);
 270   fprintf(flp,"%5ld%12.6f%12.6f%12.6f\n",maxz-minz+1-(2*iIGNOREOUTER),0.,0.,rBINWIDTH*10./bohr);
 271   for(i=0; i<nidxp; i++){
 272     v=2;
 273     if(*(top.atoms.atomname[indexp[i]][0])=='C')v=6;
 274     if(*(top.atoms.atomname[indexp[i]][0])=='N')v=7;
 275     if(*(top.atoms.atomname[indexp[i]][0])=='O')v=8;
 276     if(*(top.atoms.atomname[indexp[i]][0])=='H')v=1;
 277     if(*(top.atoms.atomname[indexp[i]][0])=='S')v=16;
 278     fprintf(flp,"%5d%12.6f%12.6f%12.6f%12.6f\n",v,0.,(double)fr.x[indexp[i]][XX]*10./bohr,(double)fr.x[indexp[i]][YY]*10./bohr,(double)fr.x[indexp[i]][ZZ]*10./bohr);
 279   }
 280
 281   tot=0;
 282   for(k=0;k<nbin[XX];k++) {
 283     if(!(k<minx||k>maxx))continue;
 284     for(j=0;j<nbin[YY];j++) {
 285       if(!(j<miny||j>maxy))continue;
 286       for(i=0;i<nbin[ZZ];i++) {
 287         if(!(i<minz||i>maxz))continue;
 288         if(bin[k][j][i]!=0){
 289           printf("A bin was not empty when it should have been empty. Programming error.\n");
 290           printf("bin[%d][%d][%d] was = %ld\n",k,j,i,bin[k][j][i]);
 291           exit(1);
 292         }
 293       }
 294     }
 295   }
 296
 297   min=999;
 298   max=0;
 299   for(k=0;k<nbin[XX];k++) {
 300     if(k<minx+iIGNOREOUTER||k>maxx-iIGNOREOUTER)continue;
 301     for(j=0;j<nbin[YY];j++) {
 302       if(j<miny+iIGNOREOUTER||j>maxy-iIGNOREOUTER)continue;
 303       for(i=0;i<nbin[ZZ];i++) {
 304         if(i<minz+iIGNOREOUTER||i>maxz-iIGNOREOUTER)continue;
 305         tot+=bin[k][j][i];
 306         if(bin[k][j][i]>max)max=bin[k][j][i];
 307         if(bin[k][j][i]<min)min=bin[k][j][i];
 308       }
 309     }
 310   }
 311
 312   numcu=(maxx-minx+1-(2*iIGNOREOUTER))*(maxy-miny+1-(2*iIGNOREOUTER))*(maxz-minz+1-(2*iIGNOREOUTER));
 313   if(bCALCDIV){
 314     norm=((double)numcu*(double)numfr) / (double)tot;
 315   }else{
 316     norm=1.0;
 317   }
 318
 319   for(k=0;k<nbin[XX];k++) {
 320     if(k<minx+iIGNOREOUTER||k>maxx-iIGNOREOUTER)continue;
 321     for(j=0;j<nbin[YY];j++) {
 322       if(j<miny+iIGNOREOUTER||j>maxy-iIGNOREOUTER)continue;
 323       for(i=0;i<nbin[ZZ];i++) {
 324         if(i<minz+iIGNOREOUTER||i>maxz-iIGNOREOUTER)continue;
 325         fprintf(flp,"%12.6f ",norm*(double)bin[k][j][i]/(double)numfr);
 326       }
 327       fprintf(flp,"\n");
 328     }
 329     fprintf(flp,"\n");
 330   }
 331   ffclose(flp);
 332
 333   /* printf("x=%d to %d\n",minx,maxx); */
 334   /* printf("y=%d to %d\n",miny,maxy); */
 335   /* printf("z=%d to %d\n",minz,maxz); */
 336
 337   if(bCALCDIV){
 338     printf("Counts per frame in all %ld cubes divided by %le\n",numcu,1.0/norm);
 339     printf("Normalized data: average %le, min %le, max %le\n",1.0,norm*(double)min/(double)numfr,norm*(double)max/(double)numfr);
 340   }else{
 341     printf("grid.cube contains counts per frame in all %ld cubes\n",numcu);
 342     printf("Raw data: average %le, min %le, max %le\n",1.0/norm,(double)min/(double)numfr,(double)max/(double)numfr);
 343   }
 344
 345   thanx(stderr);
 346
 347   return 0;
 348 }
 349