src/tools/gmx_sham.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 #include <math.h>
  39 #include <string.h>
  40 #include "statutil.h"
  41 #include "sysstuff.h"
  42 #include "typedefs.h"
  43 #include "smalloc.h"
  44 #include "macros.h"
  45 #include "gmx_fatal.h"
  46 #include "vec.h"
  47 #include "copyrite.h"
  48 #include "futil.h"
  49 #include "readinp.h"
  50 #include "statutil.h"
  51 #include "txtdump.h"
  52 #include "gstat.h"
  53 #include "xvgr.h"
  54 #include "physics.h"
  55 #include "pdbio.h"
  56 #include "matio.h"
  57 #include "gmx_ana.h"
  58
  59
  60 static int index2(int *ibox,int x,int y)
  61 {
  62   return (ibox[1]*x+y);
  63 }
  64
  65 static int index3(int *ibox,int x,int y,int z)
  66 {
  67   return (ibox[2]*(ibox[1]*x+y)+z);
  68 }
  69
  70 static int indexn(int ndim,int *ibox,int *nxyz)
  71 {
  72   int d,dd,k,kk;
  73
  74   /* Compute index in 1-D array */
  75   d = 0;
  76   for(k=0; (k<ndim); k++) {
  77     dd = nxyz[k];
  78     for(kk=k+1; (kk<ndim); kk++)
  79       dd = dd*ibox[kk];
  80     d += dd;
  81   }
  82   return d;
  83 }
  84
  85 typedef struct{
  86   int Nx; /* x grid points in unit cell */
  87   int Ny; /* y grid points in unit cell */
  88   int Nz; /* z grid points in unit cell */
  89   int dmin[3]; /* starting point x,y,z*/
  90   int dmax[3]; /* ending point x,y,z */
  91   real cell[6]; /* usual cell parameters */
  92   real * ed; /* data */
  93 } XplorMap;
  94
  95 static void lo_write_xplor(XplorMap * map,const char * file)
  96 {
  97   FILE * fp;
  98   int z,i,j,n;
  99
 100   fp = ffopen(file,"w");
 101   /* The REMARKS part is the worst part of the XPLOR format
 102    * and may cause problems with some programs
 103    */
 104   fprintf(fp,"\n       2 !NTITLE\n") ;
 105   fprintf(fp," REMARKS Energy Landscape from GROMACS\n") ;
 106   fprintf(fp," REMARKS DATE: 2004-12-21 \n") ;
 107   fprintf(fp," %7d %7d %7d %7d %7d %7d %7d %7d %7d\n",
 108           map->Nx, map->dmin[0], map->dmax[0],
 109           map->Ny, map->dmin[1], map->dmax[1],
 110           map->Nz, map->dmin[2], map->dmax[2]);
 111   fprintf(fp,"%12.5E%12.5E%12.5E%12.5E%12.5E%12.5E\n",
 112           map->cell[0],map->cell[1],map->cell[2],
 113           map->cell[3],map->cell[4],map->cell[5]);
 114   fprintf(fp, "ZYX\n") ;
 115
 116   z = map->dmin[2];
 117   for (n = 0; n < map->Nz; n++, z++) {
 118     fprintf(fp, "%8d\n", z) ;
 119     for (i = 0; i < map->Nx*map->Ny; i += 6) {
 120       for (j = 0; j < 6; j++)
 121         if (i+j < map->Nx*map->Ny)
 122           fprintf(fp, "%12.5E", map->ed[n*map->Nx*map->Ny+i+j]);
 123       fprintf(fp, "\n") ;
 124     }
 125   }
 126   fprintf(fp, "   -9999\n") ;
 127   ffclose(fp) ;
 128 }
 129
 130 static void write_xplor(const char *file,real *data,int *ibox,real dmin[],real dmax[])
 131 {
 132   XplorMap *xm;
 133   int i,j,k,n;
 134
 135   snew(xm,1);
 136   xm->Nx = ibox[XX];
 137   xm->Ny = ibox[YY];
 138   xm->Nz = ibox[ZZ];
 139   snew(xm->ed,xm->Nx*xm->Ny*xm->Nz);
 140   n=0;
 141   for(k=0; (k<xm->Nz); k++)
 142     for(j=0; (j<xm->Ny); j++)
 143       for(i=0; (i<xm->Nx); i++)
 144         xm->ed[n++] = data[index3(ibox,i,j,k)];
 145   xm->cell[0] = dmax[XX]-dmin[XX];
 146   xm->cell[1] = dmax[YY]-dmin[YY];
 147   xm->cell[2] = dmax[ZZ]-dmin[ZZ];
 148   xm->cell[3] = xm->cell[4] = xm->cell[5] = 90;
 149
 150   clear_ivec(xm->dmin);
 151   xm->dmax[XX] = ibox[XX]-1;
 152   xm->dmax[YY] = ibox[YY]-1;
 153   xm->dmax[ZZ] = ibox[ZZ]-1;
 154
 155   lo_write_xplor(xm,file);
 156
 157   sfree(xm->ed);
 158   sfree(xm);
 159 }
 160
 161 static void normalize_p_e(int len,double *P,int *nbin,real *E,real pmin)
 162 {
 163   int i;
 164   double Ptot=0;
 165
 166   for(i=0; (i<len); i++) {
 167     Ptot += P[i];
 168     if (nbin[i] > 0)
 169       E[i] = E[i]/nbin[i];
 170   }
 171   printf("Ptot = %g\n",Ptot);
 172   for(i=0; (i<len); i++) {
 173     P[i] = P[i]/Ptot;
 174     /* Have to check for pmin after normalizing to prevent "stretching"
 175      * the energies.
 176      */
 177     if (P[i] < pmin)
 178       P[i] = 0;
 179   }
 180 }
 181
 182 typedef struct {
 183   int index;
 184   real ener;
 185 } t_minimum;
 186
 187 static int comp_minima(const void *a,const void *b)
 188 {
 189   t_minimum *ma = (t_minimum *) a;
 190   t_minimum *mb = (t_minimum *) b;
 191
 192   if (ma->ener < mb->ener)
 193     return -1;
 194   else if (ma->ener > mb->ener)
 195     return 1;
 196   else
 197     return 0;
 198 }
 199
 200 static void pick_minima(const char *logfile,int *ibox,int ndim,int len,real W[])
 201 {
 202   FILE *fp;
 203   int  i,j,k,ijk,nmin;
 204   gmx_bool bMin;
 205   t_minimum *mm;
 206
 207   snew(mm,len);
 208   nmin = 0;
 209   fp = ffopen(logfile,"w");
 210   for(i=0; (i<ibox[0]); i++) {
 211     for(j=0; (j<ibox[1]); j++) {
 212       if (ndim == 3) {
 213         for(k=0; (k<ibox[2]); k++) {
 214           ijk    = index3(ibox,i,j,k);
 215           bMin   = (((i == 0)       || ((i > 0)       &&
 216                                         (W[ijk] < W[index3(ibox,i-1,j,k)]))) &&
 217                     ((i == ibox[0]-1) || ((i < ibox[0]-1) &&
 218                                         (W[ijk] < W[index3(ibox,i+1,j,k)]))) &&
 219                     ((j == 0)       || ((j > 0)       &&
 220                                         (W[ijk] < W[index3(ibox,i,j-1,k)]))) &&
 221                     ((j == ibox[1]-1) || ((j < ibox[1]-1) &&
 222                                         (W[ijk] < W[index3(ibox,i,j+1,k)]))) &&
 223                     ((k == 0)       || ((k > 0)       &&
 224                                         (W[ijk] < W[index3(ibox,i,j,k-1)]))) &&
 225                     ((k == ibox[2]-1) || ((k < ibox[2]-1) &&
 226                                         (W[ijk] < W[index3(ibox,i,j,k+1)]))));
 227           if (bMin) {
 228             fprintf(fp,"Minimum %d at index %6d energy %10.3f\n",
 229                     nmin,ijk,W[ijk]);
 230             mm[nmin].index = ijk;
 231             mm[nmin].ener  = W[ijk];
 232             nmin++;
 233           }
 234         }
 235       }
 236       else {
 237         ijk    = index2(ibox,i,j);
 238         bMin   = (((i == 0)       || ((i > 0)       &&
 239                                       (W[ijk] < W[index2(ibox,i-1,j)]))) &&
 240                   ((i == ibox[0]-1) || ((i < ibox[0]-1) &&
 241                                       (W[ijk] < W[index2(ibox,i+1,j)]))) &&
 242                   ((j == 0)       || ((j > 0)       &&
 243                                       (W[ijk] < W[index2(ibox,i,j-1)]))) &&
 244                   ((j == ibox[1]-1) || ((j < ibox[1]-1) &&
 245                                       (W[ijk] < W[index2(ibox,i,j+1)]))));
 246         if (bMin) {
 247           fprintf(fp,"Minimum %d at index %6d energy %10.3f\n",
 248                   nmin,ijk,W[ijk]);
 249           mm[nmin].index = ijk;
 250           mm[nmin].ener  = W[ijk];
 251           nmin++;
 252         }
 253       }
 254     }
 255   }
 256   qsort(mm,nmin,sizeof(mm[0]),comp_minima);
 257   fprintf(fp,"Minima sorted after energy\n");
 258   for(i=0; (i<nmin); i++) {
 259     fprintf(fp,"Minimum %d at index %6d energy %10.3f\n",
 260             i,mm[i].index,mm[i].ener);
 261   }
 262   ffclose(fp);
 263   sfree(mm);
 264 }
 265
 266 static void do_sham(const char *fn,const char *ndx,
 267                     const char *xpmP,const char *xpm,const char *xpm2,
 268                     const char *xpm3,const char *xpm4,const char *pdb,
 269                     const char *logf,
 270                     int n,int neig,real **eig,
 271                     gmx_bool bGE,int nenerT,real **enerT,
 272                     int nmap,real *mapindex,real **map,
 273                     real Tref,
 274                     real pmax,real gmax,
 275                     real *emin,real *emax,int nlevels,real pmin,
 276                     const char *mname,gmx_bool bSham,int *idim,int *ibox,
 277                     gmx_bool bXmin,real *xmin,gmx_bool bXmax,real *xmax)
 278 {
 279   FILE    *fp;
 280   real    *min_eig,*max_eig;
 281   real    *axis_x,*axis_y,*axis_z,*axis=NULL;
 282   double  *P;
 283   real    **PP,*W,*E,**WW,**EE,*S,**SS,*M,**MM,*bE;
 284   rvec    xxx;
 285   char    *buf;
 286   double  *bfac,efac,bref,Pmax,Wmin,Wmax,Winf,Emin,Emax,Einf,Smin,Smax,Sinf,Mmin,Mmax,Minf;
 287   real    *delta;
 288   int     i,j,k,imin,len,index,d,*nbin,*bindex,bi;
 289   int     *nxyz,maxbox;
 290   t_blocka *b;
 291   gmx_bool    bOutside;
 292   unsigned int flags;
 293   t_rgb   rlo  = { 0, 0, 0 };
 294   t_rgb   rhi  = { 1, 1, 1 };
 295
 296   /* Determine extremes for the eigenvectors */
 297   snew(min_eig,neig);
 298   snew(max_eig,neig);
 299   snew(nxyz,neig);
 300   snew(bfac,neig);
 301   snew(delta,neig);
 302
 303   for(i=0; (i<neig); i++) {
 304     /* Check for input constraints */
 305     min_eig[i] = max_eig[i] = eig[i][0];
 306     for(j=0; (j<n); j++) {
 307       min_eig[i] = min(min_eig[i],eig[i][j]);
 308       max_eig[i] = max(max_eig[i],eig[i][j]);
 309       delta[i]   = (max_eig[i]-min_eig[i])/(2.0*ibox[i]);
 310     }
 311     /* Add some extra space, half a bin on each side, unless the
 312      * user has set the limits.
 313      */
 314     if (bXmax) {
 315       if (max_eig[i] > xmax[i]) {
 316         gmx_warning("Your xmax[%d] value %f is smaller than the largest data point %f",i,xmax[i],max_eig[i]);
 317       }
 318       max_eig[i] = xmax[i];
 319     }
 320     else
 321       max_eig[i] += delta[i];
 322
 323     if (bXmin) {
 324       if (min_eig[i] < xmin[i]) {
 325         gmx_warning("Your xmin[%d] value %f is larger than the smallest data point %f",i,xmin[i],min_eig[i]);
 326       }
 327       min_eig[i] = xmin[i];
 328     }
 329     else
 330       min_eig[i] -= delta[i];
 331     bfac[i]     = ibox[i]/(max_eig[i]-min_eig[i]);
 332   }
 333   /* Do the binning */
 334   bref = 1/(BOLTZ*Tref);
 335   snew(bE,n);
 336   if (bGE || nenerT==2) {
 337     Emin = 1e8;
 338     for(j=0; (j<n); j++) {
 339       if (bGE)
 340         bE[j] = bref*enerT[0][j];
 341       else
 342         bE[j] = (bref - 1/(BOLTZ*enerT[1][j]))*enerT[0][j];
 343       Emin  = min(Emin,bE[j]);
 344     }
 345   }
 346   else
 347     Emin = 0;
 348   len=1;
 349   for(i=0; (i<neig); i++)
 350     len=len*ibox[i];
 351   printf("There are %d bins in the %d-dimensional histogram. Beta-Emin = %g\n",
 352          len,neig,Emin);
 353   snew(P,len);
 354   snew(W,len);
 355   snew(E,len);
 356   snew(S,len);
 357   snew(M,len);
 358   snew(nbin,len);
 359   snew(bindex,n);
 360
 361
 362   /* Loop over projections */
 363   for(j=0; (j<n); j++) {
 364     /* Loop over dimensions */
 365     bOutside = FALSE;
 366     for(i=0; (i<neig); i++) {
 367       nxyz[i] = bfac[i]*(eig[i][j]-min_eig[i]);
 368       if (nxyz[i] < 0 || nxyz[i] >= ibox[i])
 369         bOutside = TRUE;
 370     }
 371     if (!bOutside) {
 372       index = indexn(neig,ibox,nxyz);
 373       range_check(index,0,len);
 374       /* Compute the exponential factor */
 375       if (enerT)
 376         efac = exp(-bE[j]+Emin);
 377       else
 378         efac = 1;
 379       /* Apply the bin volume correction for a multi-dimensional distance */
 380       for(i=0; i<neig; i++) {
 381         if (idim[i] == 2)
 382           efac /= eig[i][j];
 383         else if (idim[i] == 3)
 384           efac /= sqr(eig[i][j]);
 385         else if (idim[i] == -1)
 386           efac /= sin(DEG2RAD*eig[i][j]);
 387       }
 388       /* Update the probability */
 389       P[index] += efac;
 390       /* Update the energy */
 391       if (enerT)
 392         E[index] += enerT[0][j];
 393       /* Statistics: which "structure" in which bin */
 394       nbin[index]++;
 395       bindex[j]=index;
 396     }
 397   }
 398   /* Normalize probability */
 399   normalize_p_e(len,P,nbin,E,pmin);
 400   Pmax = 0;
 401   /* Compute boundaries for the Free energy */
 402   Wmin = 1e8;
 403   imin = -1;
 404   Wmax = -1e8;
 405   /* Recompute Emin: it may have changed due to averaging */
 406   Emin = 1e8;
 407   Emax = -1e8;
 408   for(i=0; (i<len); i++) {
 409     if (P[i] != 0) {
 410       Pmax = max(P[i],Pmax);
 411       W[i] = -BOLTZ*Tref*log(P[i]);
 412       if (W[i] < Wmin) {
 413         Wmin = W[i];
 414         imin = i;
 415       }
 416       Emin = min(E[i],Emin);
 417       Emax = max(E[i],Emax);
 418       Wmax = max(W[i],Wmax);
 419     }
 420   }
 421   if (pmax > 0) {
 422     Pmax = pmax;
 423   }
 424   if (gmax > 0) {
 425     Wmax = gmax;
 426   } else {
 427     Wmax -= Wmin;
 428   }
 429   Winf = Wmax+1;
 430   Einf = Emax+1;
 431   Smin = Emin-Wmax;
 432   Smax = Emax-Smin;
 433   Sinf = Smax+1;
 434   /* Write out the free energy as a function of bin index */
 435   fp = ffopen(fn,"w");
 436   for(i=0; (i<len); i++)
 437     if (P[i] != 0) {
 438       W[i] -= Wmin;
 439       S[i] = E[i]-W[i]-Smin;
 440       fprintf(fp,"%5d  %10.5e  %10.5e  %10.5e\n",i,W[i],E[i],S[i]);
 441     }
 442     else {
 443       W[i] = Winf;
 444       E[i] = Einf;
 445       S[i] = Sinf;
 446     }
 447   ffclose(fp);
 448   /* Organize the structures in the bins */
 449   snew(b,1);
 450   snew(b->index,len+1);
 451   snew(b->a,n);
 452   b->index[0] = 0;
 453   for(i=0; (i<len); i++) {
 454     b->index[i+1] = b->index[i]+nbin[i];
 455     nbin[i] = 0;
 456   }
 457   for(i=0; (i<n); i++) {
 458     bi = bindex[i];
 459     b->a[b->index[bi]+nbin[bi]] = i;
 460     nbin[bi]++;
 461   }
 462   /* Consistency check */
 463   /* This no longer applies when we allow the plot to be smaller
 464      than the sampled space.
 465   for(i=0; (i<len); i++) {
 466     if (nbin[i] != (b->index[i+1] - b->index[i]))
 467       gmx_fatal(FARGS,"nbin[%d] = %d, should be %d",i,nbin[i],
 468                 b->index[i+1] - b->index[i]);
 469   }
 470   */
 471   /* Write the index file */
 472   fp = ffopen(ndx,"w");
 473   for(i=0; (i<len); i++) {
 474     if (nbin[i] > 0) {
 475       fprintf(fp,"[ %d ]\n",i);
 476       for(j=b->index[i]; (j<b->index[i+1]); j++)
 477         fprintf(fp,"%d\n",b->a[j]+1);
 478     }
 479   }
 480   ffclose(fp);
 481   snew(axis_x,ibox[0]+1);
 482   snew(axis_y,ibox[1]+1);
 483   snew(axis_z,ibox[2]+1);
 484   maxbox = max(ibox[0],max(ibox[1],ibox[2]));
 485   snew(PP,maxbox*maxbox);
 486   snew(WW,maxbox*maxbox);
 487   snew(EE,maxbox*maxbox);
 488   snew(SS,maxbox*maxbox);
 489   for(i=0; (i<min(neig,3)); i++) {
 490     switch (i) {
 491     case 0: axis = axis_x; break;
 492     case 1: axis = axis_y; break;
 493     case 2: axis = axis_z; break;
 494     default: break;
 495     }
 496     for(j=0; j<=ibox[i]; j++)
 497       axis[j] = min_eig[i] + j/bfac[i];
 498   }
 499   if (map) {
 500     snew(M,len);
 501     snew(MM,maxbox*maxbox);
 502     for(i=0; (i<ibox[0]); i++)
 503       MM[i] = &(M[i*ibox[1]]);
 504     Mmin = 1e8;
 505     Mmax = -1e8;
 506     for(i=0; (i<nmap); i++) {
 507       Mmin = min(Mmin,map[0][i]);
 508       Mmax = max(Mmax,map[0][i]);
 509     }
 510     Minf = Mmax*1.05;
 511     for(i=0; (i<len); i++)
 512       M[i] = Minf;
 513     for(i=0; (i<nmap); i++) {
 514       index = gmx_nint(mapindex[i]);
 515       if (index >= len)
 516         gmx_fatal(FARGS,"Number of bins in file from -mdata option does not correspond to current analysis");
 517
 518       if (P[index] != 0)
 519         M[index] = map[0][i];
 520     }
 521   }
 522   else {
 523     MM = NULL;
 524     Minf = NOTSET;
 525   }
 526   pick_minima(logf,ibox,neig,len,W);
 527   if (gmax <= 0)
 528     gmax = Winf;
 529   flags = MAT_SPATIAL_X | MAT_SPATIAL_Y;
 530   if (neig == 2) {
 531     /* Dump to XPM file */
 532     snew(PP,ibox[0]);
 533     for(i=0; (i<ibox[0]); i++) {
 534       snew(PP[i],ibox[1]);
 535       for(j=0; j<ibox[1]; j++) {
 536         PP[i][j] = P[i*ibox[1]+j];
 537       }
 538       WW[i] = &(W[i*ibox[1]]);
 539       EE[i] = &(E[i*ibox[1]]);
 540       SS[i] = &(S[i*ibox[1]]);
 541     }
 542     fp = ffopen(xpmP,"w");
 543     write_xpm(fp,flags,"Probability Distribution","","PC1","PC2",
 544               ibox[0],ibox[1],axis_x,axis_y,PP,0,Pmax,rlo,rhi,&nlevels);
 545     ffclose(fp);
 546     fp = ffopen(xpm,"w");
 547     write_xpm(fp,flags,"Gibbs Energy Landscape","G (kJ/mol)","PC1","PC2",
 548               ibox[0],ibox[1],axis_x,axis_y,WW,0,gmax,rlo,rhi,&nlevels);
 549     ffclose(fp);
 550     fp = ffopen(xpm2,"w");
 551     write_xpm(fp,flags,"Enthalpy Landscape","H (kJ/mol)","PC1","PC2",
 552               ibox[0],ibox[1],axis_x,axis_y,EE,
 553               emin ? *emin : Emin,emax ? *emax : Einf,rlo,rhi,&nlevels);
 554     ffclose(fp);
 555     fp = ffopen(xpm3,"w");
 556     write_xpm(fp,flags,"Entropy Landscape","TDS (kJ/mol)","PC1","PC2",
 557               ibox[0],ibox[1],axis_x,axis_y,SS,0,Sinf,rlo,rhi,&nlevels);
 558     ffclose(fp);
 559     if (map) {
 560       fp = ffopen(xpm4,"w");
 561       write_xpm(fp,flags,"Custom Landscape",mname,"PC1","PC2",
 562                 ibox[0],ibox[1],axis_x,axis_y,MM,0,Minf,rlo,rhi,&nlevels);
 563       ffclose(fp);
 564     }
 565   }
 566   else if (neig == 3) {
 567     /* Dump to PDB file */
 568     fp = ffopen(pdb,"w");
 569     for(i=0; (i<ibox[0]); i++) {
 570       xxx[XX] = 3*(i+0.5-ibox[0]/2);
 571       for(j=0; (j<ibox[1]); j++) {
 572         xxx[YY] = 3*(j+0.5-ibox[1]/2);
 573         for(k=0; (k<ibox[2]); k++) {
 574           xxx[ZZ] = 3*(k+0.5-ibox[2]/2);
 575           index = index3(ibox,i,j,k);
 576           if (P[index] > 0)
 577             fprintf(fp,"%-6s%5u  %-4.4s%3.3s  %4d    %8.3f%8.3f%8.3f%6.2f%6.2f\n",
 578                     "ATOM",(index+1) %10000,"H","H",(index+1)%10000,
 579                     xxx[XX],xxx[YY],xxx[ZZ],1.0,W[index]);
 580         }
 581       }
 582     }
 583     ffclose(fp);
 584     write_xplor("out.xplor",W,ibox,min_eig,max_eig);
 585     if (map)
 586       write_xplor("user.xplor",M,ibox,min_eig,max_eig);
 587     nxyz[XX] = imin/(ibox[1]*ibox[2]);
 588     nxyz[YY] = (imin-nxyz[XX]*ibox[1]*ibox[2])/ibox[2];
 589     nxyz[ZZ] = imin % ibox[2];
 590     for(i=0; (i<ibox[0]); i++) {
 591       snew(WW[i],maxbox);
 592       for(j=0; (j<ibox[1]); j++)
 593         WW[i][j] = W[index3(ibox,i,j,nxyz[ZZ])];
 594     }
 595     snew(buf,strlen(xpm)+4);
 596     sprintf(buf,"%s",xpm);
 597     sprintf(&buf[strlen(xpm)-4],"12.xpm");
 598     fp = ffopen(buf,"w");
 599     write_xpm(fp,flags,"Gibbs Energy Landscape","W (kJ/mol)","PC1","PC2",
 600               ibox[0],ibox[1],axis_x,axis_y,WW,0,gmax,rlo,rhi,&nlevels);
 601     ffclose(fp);
 602     for(i=0; (i<ibox[0]); i++) {
 603       for(j=0; (j<ibox[2]); j++)
 604         WW[i][j] = W[index3(ibox,i,nxyz[YY],j)];
 605     }
 606     sprintf(&buf[strlen(xpm)-4],"13.xpm");
 607     fp = ffopen(buf,"w");
 608     write_xpm(fp,flags,"SHAM Energy Landscape","kJ/mol","PC1","PC3",
 609               ibox[0],ibox[2],axis_x,axis_z,WW,0,gmax,rlo,rhi,&nlevels);
 610     ffclose(fp);
 611     for(i=0; (i<ibox[1]); i++) {
 612       for(j=0; (j<ibox[2]); j++)
 613         WW[i][j] = W[index3(ibox,nxyz[XX],i,j)];
 614     }
 615     sprintf(&buf[strlen(xpm)-4],"23.xpm");
 616     fp = ffopen(buf,"w");
 617     write_xpm(fp,flags,"SHAM Energy Landscape","kJ/mol","PC2","PC3",
 618               ibox[1],ibox[2],axis_y,axis_z,WW,0,gmax,rlo,rhi,&nlevels);
 619     ffclose(fp);
 620     sfree(buf);
 621   }
 622   if (map) {
 623     sfree(MM);
 624     sfree(M);
 625   }
 626 }
 627
 628 static void ehisto(const char *fh,int n,real **enerT, const output_env_t oenv)
 629 {
 630   FILE *fp;
 631   int  i,j,k,nbin,blength;
 632   int  *bindex;
 633   real *T,bmin,bmax,bwidth;
 634   int  **histo;
 635
 636   bmin =  1e8;
 637   bmax = -1e8;
 638   snew(bindex,n);
 639   snew(T,n);
 640   nbin = 0;
 641   for(j=1; (j<n); j++) {
 642     for(k=0; (k<nbin); k++) {
 643       if (T[k] == enerT[1][j]) {
 644         bindex[j] = k;
 645         break;
 646       }
 647     }
 648     if (k == nbin) {
 649       bindex[j] = nbin;
 650       T[nbin]   = enerT[1][j];
 651       nbin++;
 652     }
 653     bmin = min(enerT[0][j],bmin);
 654     bmax = max(enerT[0][j],bmax);
 655 }
 656   bwidth  = 1.0;
 657   blength = (bmax - bmin)/bwidth + 2;
 658   snew(histo,nbin);
 659   for(i=0; (i<nbin); i++) {
 660     snew(histo[i],blength);
 661   }
 662   for(j=0; (j<n); j++) {
 663     k = (enerT[0][j]-bmin)/bwidth;
 664     histo[bindex[j]][k]++;
 665   }
 666   fp = xvgropen(fh,"Energy distribution","E (kJ/mol)","",oenv);
 667   for(j=0; (j<blength); j++) {
 668     fprintf(fp,"%8.3f",bmin+j*bwidth);
 669     for(k=0; (k<nbin); k++) {
 670       fprintf(fp,"  %6d",histo[k][j]);
 671     }
 672     fprintf(fp,"\n");
 673   }
 674   ffclose(fp);
 675 }
 676
 677 int gmx_sham(int argc,char *argv[])
 678 {
 679   const char *desc[] = {
 680     "g_sham makes multi-dimensional free-energy, enthalpy and entropy plots.",
 681     "g_sham reads one or more xvg files and analyzes data sets.",
 682     "g_sham basic purpose is plotting Gibbs free energy landscapes",
 683     "(option [TT]-ls[tt])",
 684     "by Bolzmann inverting multi-dimensional histograms (option [TT]-lp[tt])",
 685     "but it can also",
 686     "make enthalpy (option [TT]-lsh[tt]) and entropy (option [TT]-lss[tt])",
 687     "plots. The histograms can be made for any quantities the user supplies.",
 688     "A line in the input file may start with a time",
 689     "(see option [TT]-time[tt]) and any number of y values may follow.",
 690     "Multiple sets can also be",
 691     "read when they are separated by & (option [TT]-n[tt]),",
 692     "in this case only one y value is read from each line.",
 693     "All lines starting with # and @ are skipped.",
 694     "[PAR]",
 695     "Option [TT]-ge[tt] can be used to supply a file with free energies",
 696     "when the ensemble is not a Boltzmann ensemble, but needs to be biased",
 697     "by this free energy. One free energy value is required for each",
 698     "(multi-dimensional) data point in the [TT]-f[tt] input.",
 699     "[PAR]",
 700     "Option [TT]-ene[tt] can be used to supply a file with energies.",
 701     "These energies are used as a weighting function in the single",
 702     "histogram analysis method due to Kumar et. al. When also temperatures",
 703     "are supplied (as a second column in the file) an experimental",
 704     "weighting scheme is applied. In addition the vales",
 705     "are used for making enthalpy and entropy plots.",
 706     "[PAR]",
 707     "With option [TT]-dim[tt] dimensions can be gives for distances.",
 708     "When a distance is 2- or 3-dimensional, the circumference or surface",
 709     "sampled by two particles increases with increasing distance.",
 710     "Depending on what one would like to show, one can choose to correct",
 711     "the histogram and free-energy for this volume effect.",
 712     "The probability is normalized by r and r^2 for a dimension of 2 and 3",
 713     "respectively.",
 714     "A value of -1 is used to indicate an angle in degrees between two",
 715     "vectors: a sin(angle) normalization will be applied.",
 716     "Note that for angles between vectors the inner-product or cosine",
 717     "is the natural quantity to use, as it will produce bins of the same",
 718     "volume."
 719   };
 720   static real tb=-1,te=-1,frac=0.5,filtlen=0,binwidth=0.1;
 721   static gmx_bool bHaveT=TRUE,bDer=FALSE,bSubAv=TRUE,bAverCorr=FALSE,bXYdy=FALSE;
 722   static gmx_bool bEESEF=FALSE,bEENLC=FALSE,bEeFitAc=FALSE,bPower=FALSE;
 723   static gmx_bool bShamEner=TRUE,bSham=TRUE;
 724   static real Tref=298.15,pmin=0,ttol=0,pmax=0,gmax=0,emin=0,emax=0;
 725   static rvec nrdim = {1,1,1};
 726   static rvec nrbox = {32,32,32};
 727   static rvec xmin  = {0,0,0}, xmax={1,1,1};
 728   static int  nsets_in=1,nb_min=4,resol=10,nlevels=25;
 729   static const char *mname="";
 730   t_pargs pa[] = {
 731     { "-time",    FALSE, etBOOL, {&bHaveT},
 732       "Expect a time in the input" },
 733     { "-b",       FALSE, etREAL, {&tb},
 734       "First time to read from set" },
 735     { "-e",       FALSE, etREAL, {&te},
 736       "Last time to read from set" },
 737     { "-ttol",     FALSE, etREAL, {&ttol},
 738       "Tolerance on time in appropriate units (usually ps)" },
 739     { "-n",       FALSE, etINT, {&nsets_in},
 740       "Read # sets separated by &" },
 741     { "-d",       FALSE, etBOOL, {&bDer},
 742         "Use the derivative" },
 743     { "-bw",      FALSE, etREAL, {&binwidth},
 744       "Binwidth for the distribution" },
 745     { "-sham",    FALSE, etBOOL, {&bSham},
 746       "Turn off energy weighting even if energies are given" },
 747     { "-tsham",   FALSE, etREAL, {&Tref},
 748       "Temperature for single histogram analysis" },
 749     { "-pmin",    FALSE, etREAL, {&pmin},
 750       "Minimum probability. Anything lower than this will be set to zero" },
 751     { "-dim",     FALSE, etRVEC, {nrdim},
 752       "Dimensions for distances, used for volume correction (max 3 values, dimensions > 3 will get the same value as the last)" },
 753     { "-ngrid",   FALSE, etRVEC, {nrbox},
 754       "Number of bins for energy landscapes (max 3 values, dimensions > 3 will get the same value as the last)" },
 755     { "-xmin",    FALSE, etRVEC, {xmin},
 756       "Minimum for the axes in energy landscape (see above for > 3 dimensions)" },
 757     { "-xmax",    FALSE, etRVEC, {xmax},
 758       "Maximum for the axes in energy landscape (see above for > 3 dimensions)" },
 759     { "-pmax",    FALSE, etREAL, {&pmax},
 760       "Maximum probability in output, default is calculate" },
 761     { "-gmax",    FALSE, etREAL, {&gmax},
 762       "Maximum free energy in output, default is calculate" },
 763     { "-emin",    FALSE, etREAL, {&emin},
 764       "Minimum enthalpy in output, default is calculate" },
 765     { "-emax",    FALSE, etREAL, {&emax},
 766       "Maximum enthalpy in output, default is calculate" },
 767     { "-nlevels", FALSE, etINT,  {&nlevels},
 768       "Number of levels for energy landscape" },
 769     { "-mname",   FALSE, etSTR,  {&mname},
 770       "Legend label for the custom landscape" },
 771   };
 772 #define NPA asize(pa)
 773
 774   FILE     *out;
 775   int      n,e_n,d_n,nlast,s,nset,e_nset,d_nset,i,j=0,*idim,*ibox;
 776   real     **val,**et_val,**dt_val,*t,*e_t,e_dt,d_dt,*d_t,dt,tot,error;
 777   real     *rmin,*rmax;
 778   double   *av,*sig,cum1,cum2,cum3,cum4,db;
 779   const char     *fn_ge,*fn_ene;
 780   output_env_t oenv;
 781
 782   t_filenm fnm[] = {
 783     { efXVG, "-f",    "graph",    ffREAD   },
 784     { efXVG, "-ge",   "gibbs",    ffOPTRD  },
 785     { efXVG, "-ene",  "esham",    ffOPTRD  },
 786     { efXVG, "-dist", "ener",     ffOPTWR  },
 787     { efXVG, "-histo","edist",    ffOPTWR  },
 788     { efNDX, "-bin",  "bindex",   ffOPTWR  },
 789     { efXPM, "-lp",   "prob",     ffOPTWR  },
 790     { efXPM, "-ls",   "gibbs",    ffOPTWR  },
 791     { efXPM, "-lsh",  "enthalpy", ffOPTWR  },
 792     { efXPM, "-lss",  "entropy",  ffOPTWR  },
 793     { efXPM, "-map",  "map",      ffOPTWR  },
 794     { efPDB, "-ls3",  "gibbs3",   ffOPTWR  },
 795     { efXVG, "-mdata","mapdata",  ffOPTWR  },
 796     { efLOG, "-g",    "shamlog",  ffOPTWR  }
 797   };
 798 #define NFILE asize(fnm)
 799
 800   int     npargs;
 801
 802   npargs = asize(pa);
 803
 804   CopyRight(stderr,argv[0]);
 805   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_BE_NICE ,
 806                     NFILE,fnm,npargs,pa,asize(desc),desc,0,NULL,&oenv);
 807
 808   val=read_xvg_time(opt2fn("-f",NFILE,fnm),bHaveT,
 809                     opt2parg_bSet("-b",npargs,pa),tb-ttol,
 810                     opt2parg_bSet("-e",npargs,pa),te+ttol,
 811                     nsets_in,&nset,&n,&dt,&t);
 812   printf("Read %d sets of %d points, dt = %g\n\n",nset,n,dt);
 813
 814   fn_ge  = opt2fn_null("-ge",NFILE,fnm);
 815   fn_ene = opt2fn_null("-ene",NFILE,fnm);
 816
 817   if (fn_ge && fn_ene)
 818     gmx_fatal(FARGS,"Can not do free energy and energy corrections at the same time");
 819
 820   if (fn_ge || fn_ene) {
 821     et_val=read_xvg_time(fn_ge ? fn_ge : fn_ene,bHaveT,
 822                          opt2parg_bSet("-b",npargs,pa),tb-ttol,
 823                          opt2parg_bSet("-e",npargs,pa),te+ttol,
 824                          1,&e_nset,&e_n,&e_dt,&e_t);
 825     if (fn_ge) {
 826       if (e_nset != 1)
 827         gmx_fatal(FARGS,"Can only handle one free energy component in %s",
 828                 fn_ge);
 829     } else {
 830       if (e_nset!=1 && e_nset!=2)
 831         gmx_fatal(FARGS,"Can only handle one energy component or one energy and one T in %s",
 832                   fn_ene);
 833     }
 834     if (e_n != n)
 835       gmx_fatal(FARGS,"Number of energies (%d) does not match number of entries (%d) in %s",e_n,n,opt2fn("-f",NFILE,fnm));
 836   }
 837   else
 838     et_val = NULL;
 839
 840   if (opt2fn_null("-mdata",NFILE,fnm) != NULL) {
 841     dt_val=read_xvg_time(opt2fn("-mdata",NFILE,fnm),bHaveT,
 842                          FALSE,tb,FALSE,te,
 843                          nsets_in,&d_nset,&d_n,&d_dt,&d_t);
 844     if (d_nset != 1)
 845       gmx_fatal(FARGS,"Can only handle one mapping data column in %s",
 846                 opt2fn("-mdata",NFILE,fnm));
 847   }
 848   else
 849     dt_val = NULL;
 850
 851   if (fn_ene && et_val)
 852     ehisto(opt2fn("-histo",NFILE,fnm),e_n,et_val,oenv);
 853
 854   snew(idim,nset);
 855   snew(ibox,nset);
 856   snew(rmin,nset);
 857   snew(rmax,nset);
 858   for(i=0; (i<min(3,nset)); i++) {
 859     idim[i] = nrdim[i];
 860     ibox[i] = nrbox[i];
 861     rmin[i] = xmin[i];
 862     rmax[i] = xmax[i];
 863   }
 864   for(; (i<nset); i++) {
 865     idim[i] = nrdim[2];
 866     ibox[i] = nrbox[2];
 867     rmin[i] = xmin[2];
 868     rmax[i] = xmax[2];
 869   }
 870
 871   do_sham(opt2fn("-dist",NFILE,fnm),opt2fn("-bin",NFILE,fnm),
 872           opt2fn("-lp",NFILE,fnm),
 873           opt2fn("-ls",NFILE,fnm),opt2fn("-lsh",NFILE,fnm),
 874           opt2fn("-lss",NFILE,fnm),opt2fn("-map",NFILE,fnm),
 875           opt2fn("-ls3",NFILE,fnm),opt2fn("-g",NFILE,fnm),
 876           n,nset,val,fn_ge!=NULL,e_nset,et_val,d_n,d_t,dt_val,Tref,
 877           pmax,gmax,
 878           opt2parg_bSet("-emin",NPA,pa) ? &emin : NULL,
 879           opt2parg_bSet("-emax",NPA,pa) ? &emax : NULL,
 880           nlevels,pmin,
 881           mname,bSham,idim,ibox,
 882           opt2parg_bSet("-xmin",NPA,pa),rmin,
 883           opt2parg_bSet("-xmax",NPA,pa),rmax);
 884
 885   thanx(stderr);
 886
 887   return 0;
 888 }
 889