src/tools/gmx_spatial.c

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