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