added Verlet scheme and NxN non-bonded functionality

[gromacs.git] / src / tools / gmx_sham.c

/*
 * 
 *                This source code is part of
 * 
 *                 G   R   O   M   A   C   S
 * 
 *          GROningen MAchine for Chemical Simulations
 * 
 *                        VERSION 3.2.0
 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
 * Copyright (c) 2001-2004, The GROMACS development team,
 * check out http://www.gromacs.org for more information.

 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 * 
 * If you want to redistribute modifications, please consider that
 * scientific software is very special. Version control is crucial -
 * bugs must be traceable. We will be happy to consider code for
 * inclusion in the official distribution, but derived work must not
 * be called official GROMACS. Details are found in the README & COPYING
 * files - if they are missing, get the official version at www.gromacs.org.
 * 
 * To help us fund GROMACS development, we humbly ask that you cite
 * the papers on the package - you can find them in the top README file.
 * 
 * For more info, check our website at http://www.gromacs.org
 * 
 * And Hey:
 * Green Red Orange Magenta Azure Cyan Skyblue
 */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <math.h>
#include <string.h>
#include "statutil.h"
#include "sysstuff.h"
#include "typedefs.h"
#include "smalloc.h"
#include "macros.h"
#include "gmx_fatal.h"
#include "vec.h"
#include "copyrite.h"
#include "futil.h"
#include "readinp.h"
#include "statutil.h"
#include "txtdump.h"
#include "gstat.h"
#include "xvgr.h"
#include "physics.h"
#include "pdbio.h"
#include "matio.h"
#include "gmx_ana.h"


static int index2(int *ibox,int x,int y) 
{
  return (ibox[1]*x+y);
}

static int index3(int *ibox,int x,int y,int z) 
{
  return (ibox[2]*(ibox[1]*x+y)+z);
}

static gmx_large_int_t indexn(int ndim, const int *ibox, const int *nxyz)
{
    gmx_large_int_t d, dd;
    int k, kk;
  
  /* Compute index in 1-D array */
  d = 0;
  for(k=0; (k<ndim); k++) {
    dd = nxyz[k];
    for(kk=k+1; (kk<ndim); kk++)
      dd = dd*ibox[kk];
    d += dd;
  }
  return d;
}

typedef struct{
  int Nx; /* x grid points in unit cell */
  int Ny; /* y grid points in unit cell */
  int Nz; /* z grid points in unit cell */
  int dmin[3]; /* starting point x,y,z*/
  int dmax[3]; /* ending point x,y,z */
  real cell[6]; /* usual cell parameters */
  real * ed; /* data */
} XplorMap;

static void lo_write_xplor(XplorMap * map,const char * file)
{
  FILE * fp;
  int z,i,j,n;
  
  fp = ffopen(file,"w");
  /* The REMARKS part is the worst part of the XPLOR format
   * and may cause problems with some programs 
   */
  fprintf(fp,"\n       2 !NTITLE\n") ;
  fprintf(fp," REMARKS Energy Landscape from GROMACS\n") ;
  fprintf(fp," REMARKS DATE: 2004-12-21 \n") ;
  fprintf(fp," %7d %7d %7d %7d %7d %7d %7d %7d %7d\n",
          map->Nx, map->dmin[0], map->dmax[0],
          map->Ny, map->dmin[1], map->dmax[1],
          map->Nz, map->dmin[2], map->dmax[2]);
  fprintf(fp,"%12.5E%12.5E%12.5E%12.5E%12.5E%12.5E\n",
          map->cell[0],map->cell[1],map->cell[2],
          map->cell[3],map->cell[4],map->cell[5]);
  fprintf(fp, "ZYX\n") ;
  
  z = map->dmin[2];
  for (n = 0; n < map->Nz; n++, z++) {
    fprintf(fp, "%8d\n", z) ;
    for (i = 0; i < map->Nx*map->Ny; i += 6) {
      for (j = 0; j < 6; j++)
        if (i+j < map->Nx*map->Ny)
          fprintf(fp, "%12.5E", map->ed[n*map->Nx*map->Ny+i+j]);
      fprintf(fp, "\n") ;
    }
  }
  fprintf(fp, "   -9999\n") ;
  ffclose(fp) ;
}

static void write_xplor(const char *file,real *data,int *ibox,real dmin[],real dmax[])
{
  XplorMap *xm;
  int i,j,k,n;
    
  snew(xm,1);
  xm->Nx = ibox[XX];
  xm->Ny = ibox[YY];
  xm->Nz = ibox[ZZ];
  snew(xm->ed,xm->Nx*xm->Ny*xm->Nz);
  n=0;
  for(k=0; (k<xm->Nz); k++)
    for(j=0; (j<xm->Ny); j++)
      for(i=0; (i<xm->Nx); i++)
        xm->ed[n++] = data[index3(ibox,i,j,k)];
  xm->cell[0] = dmax[XX]-dmin[XX];
  xm->cell[1] = dmax[YY]-dmin[YY];
  xm->cell[2] = dmax[ZZ]-dmin[ZZ];
  xm->cell[3] = xm->cell[4] = xm->cell[5] = 90;

  clear_ivec(xm->dmin);
  xm->dmax[XX] = ibox[XX]-1;
  xm->dmax[YY] = ibox[YY]-1;
  xm->dmax[ZZ] = ibox[ZZ]-1;
  
  lo_write_xplor(xm,file);
  
  sfree(xm->ed);
  sfree(xm);
}

static void normalize_p_e(int len,double *P,int *nbin,real *E,real pmin)
{
  int i;
  double Ptot=0;
  
  for(i=0; (i<len); i++) {
    Ptot += P[i];
    if (nbin[i] > 0)
      E[i] = E[i]/nbin[i];
  }
  printf("Ptot = %g\n",Ptot);
  for(i=0; (i<len); i++) {
    P[i] = P[i]/Ptot;
    /* Have to check for pmin after normalizing to prevent "stretching"
     * the energies.
     */
    if (P[i] < pmin)
      P[i] = 0;
  }
}

typedef struct {
  gmx_large_int_t index;
  real ener;
} t_minimum;

static int comp_minima(const void *a,const void *b)
{
  t_minimum *ma = (t_minimum *) a;
  t_minimum *mb = (t_minimum *) b;
  
  if (ma->ener < mb->ener)
    return -1;
  else if (ma->ener > mb->ener)
    return 1;
  else
    return 0;
}

static inline
void print_minimum(FILE *fp, int num, const t_minimum *min)
{
    fprintf(fp,
            "Minimum %d at index " gmx_large_int_pfmt " energy %10.3f\n",
            num, min->index, min->ener);
}

static inline
void add_minimum(FILE *fp, int num, const t_minimum *min, t_minimum *mm)
{
    print_minimum(fp, num, min);
    mm[num].index = min->index;
    mm[num].ener  = min->ener;
}

static inline
gmx_bool is_local_minimum_from_below(const t_minimum *this_min,
                                     int dimension_index,
                                     int dimension_min,
                                     int neighbour_index,
                                     real *W)
{
    return ((dimension_index == dimension_min) ||
            ((dimension_index > dimension_min) &&
             (this_min->ener < W[neighbour_index])));
    /* Note over/underflow within W cannot occur. */
}

static inline
gmx_bool is_local_minimum_from_above(const t_minimum *this_min,
                                     int dimension_index,
                                     int dimension_max,
                                     int neighbour_index,
                                     real *W)
{
    return ((dimension_index == dimension_max) ||
            ((dimension_index < dimension_max) &&
             (this_min->ener < W[neighbour_index])));
    /* Note over/underflow within W cannot occur. */
}

static void pick_minima(const char *logfile,int *ibox,int ndim,int len,real W[])
{
    FILE *fp;
    int  i,j,k,nmin;
    t_minimum *mm, this_min;
    int *this_point;
    int loopmax, loopcounter;
  
    snew(mm,len);
    nmin = 0;
    fp = ffopen(logfile,"w");
    /* Loop over each element in the array of dimenion ndim seeking
     * minima with respect to every dimension. Specialized loops for
     * speed with ndim == 2 and ndim == 3. */
    switch(ndim)
    {
    case 0:
        /* This is probably impossible to reach anyway. */
        break;
    case 2:
        for(i=0; (i<ibox[0]); i++) {
            for(j=0; (j<ibox[1]); j++) {
                /* Get the index of this point in the flat array */
                this_min.index = index2(ibox,i,j);
                this_min.ener = W[this_min.index];
                if (is_local_minimum_from_below(&this_min, i, 0,         index2(ibox,i-1,j  ), W) &&
                    is_local_minimum_from_above(&this_min, i, ibox[0]-1, index2(ibox,i+1,j  ), W) &&
                    is_local_minimum_from_below(&this_min, j, 0,         index2(ibox,i  ,j-1), W) &&
                    is_local_minimum_from_above(&this_min, j, ibox[1]-1, index2(ibox,i  ,j+1), W))
                {
                    add_minimum(fp, nmin, &this_min, mm);
                    nmin++;
                }
            }
        }
        break;
    case 3:
        for(i=0; (i<ibox[0]); i++) {
            for(j=0; (j<ibox[1]); j++) {
                for(k=0; (k<ibox[2]); k++) {
                    /* Get the index of this point in the flat array */
                    this_min.index = index3(ibox,i,j,k);
                    this_min.ener = W[this_min.index];
                    if (is_local_minimum_from_below(&this_min, i, 0,         index3(ibox,i-1,j  ,k  ), W) &&
                        is_local_minimum_from_above(&this_min, i, ibox[0]-1, index3(ibox,i+1,j  ,k  ), W) &&
                        is_local_minimum_from_below(&this_min, j, 0,         index3(ibox,i  ,j-1,k  ), W) &&
                        is_local_minimum_from_above(&this_min, j, ibox[1]-1, index3(ibox,i  ,j+1,k  ), W) &&
                        is_local_minimum_from_below(&this_min, k, 0,         index3(ibox,i  ,j  ,k-1), W) &&
                        is_local_minimum_from_above(&this_min, k, ibox[2]-1, index3(ibox,i  ,j  ,k+1), W))
                    {
                        add_minimum(fp, nmin, &this_min, mm);
                        nmin++;
                    }
                }
            }
        }
        break;
    default:
        /* Note this treats ndim == 1 and ndim > 3 */

        /* Set up an ndim-dimensional vector to loop over the points
         * on the grid. (0,0,0, ... 0) is an acceptable place to
         * start. */
        snew(this_point, ndim);

        /* Determine the number of points of the ndim-dimensional
         * grid. */
        loopmax = ibox[0];
        for (i = 1; i < ndim; i++)
        {
            loopmax *= ibox[i];
        }

        loopcounter = 0;
        while (loopmax > loopcounter)
        {
            gmx_bool bMin = TRUE;

            /* Get the index of this_point in the flat array */
            this_min.index = indexn(ndim, ibox, this_point);
            this_min.ener = W[this_min.index];

            /* Is this_point a minimum from above and below in each
             * dimension? */
            for (i = 0; bMin && (i < ndim); i++)
            {
                /* Save the index of this_point within the curent
                 * dimension so we can change that index in the
                 * this_point array for use with indexn(). */
                int index = this_point[i];
                this_point[i]--;
                bMin = bMin &&
                    is_local_minimum_from_below(&this_min, index, 0,         indexn(ndim, ibox, this_point), W);
                this_point[i] += 2;
                bMin = bMin &&
                    is_local_minimum_from_above(&this_min, index, ibox[i]-1, indexn(ndim, ibox, this_point), W);
                this_point[i]--;
            }
            if (bMin)
            {
                add_minimum(fp, nmin, &this_min, mm);
                nmin++;
            }

            /* update global loop counter */
            loopcounter++;

            /* Avoid underflow of this_point[i] */
            if (loopmax > loopcounter)
            {
                /* update this_point non-recursively */
                i = ndim-1;
                this_point[i]++;
                while (ibox[i] == this_point[i])
                {
                    this_point[i] = 0;
                    i--;
                    /* this_point[i] cannot underflow because
                     * loopmax > loopcounter. */
                    this_point[i]++;
                }
            }
        }

        sfree(this_point);
        break;
    }
    qsort(mm,nmin,sizeof(mm[0]),comp_minima);
    fprintf(fp,"Minima sorted after energy\n");
    for(i=0; (i<nmin); i++)
    {
        print_minimum(fp, i, &mm[i]);
    }
    ffclose(fp);
    sfree(mm);
}

static void do_sham(const char *fn,const char *ndx,
                    const char *xpmP,const char *xpm,const char *xpm2,
                    const char *xpm3,const char *xpm4,const char *pdb,
                    const char *logf,
                    int n,int neig,real **eig,
                    gmx_bool bGE,int nenerT,real **enerT,
                    int nmap,real *mapindex,real **map,
                    real Tref,
                    real pmax,real gmax,
                    real *emin,real *emax,int nlevels,real pmin,
                    const char *mname,gmx_bool bSham,int *idim,int *ibox,
                    gmx_bool bXmin,real *xmin,gmx_bool bXmax,real *xmax)
{
  FILE    *fp;
  real    *min_eig,*max_eig;
  real    *axis_x,*axis_y,*axis_z,*axis=NULL;
  double  *P;
  real    **PP,*W,*E,**WW,**EE,*S,**SS,*M,**MM,*bE;
  rvec    xxx;
  char    *buf;
  double  *bfac,efac,bref,Pmax,Wmin,Wmax,Winf,Emin,Emax,Einf,Smin,Smax,Sinf,Mmin,Mmax,Minf;
  real    *delta;
  int     i,j,k,imin,len,index,d,*nbin,*bindex,bi;
  int     *nxyz,maxbox;
  t_blocka *b;
  gmx_bool    bOutside;
  unsigned int flags;
  t_rgb   rlo  = { 0, 0, 0 };
  t_rgb   rhi  = { 1, 1, 1 };
  
  /* Determine extremes for the eigenvectors */
  snew(min_eig,neig);
  snew(max_eig,neig);
  snew(nxyz,neig);
  snew(bfac,neig);
  snew(delta,neig);
  
  for(i=0; (i<neig); i++) {
    /* Check for input constraints */
    min_eig[i] = max_eig[i] = eig[i][0];
    for(j=0; (j<n); j++) {
      min_eig[i] = min(min_eig[i],eig[i][j]);
      max_eig[i] = max(max_eig[i],eig[i][j]);
      delta[i]   = (max_eig[i]-min_eig[i])/(2.0*ibox[i]);
    }
    /* Add some extra space, half a bin on each side, unless the
     * user has set the limits.
     */
    if (bXmax) {
      if (max_eig[i] > xmax[i]) {
        gmx_warning("Your xmax[%d] value %f is smaller than the largest data point %f",i,xmax[i],max_eig[i]);
      }
      max_eig[i] = xmax[i];
    }
    else
      max_eig[i] += delta[i];
    
    if (bXmin) {
      if (min_eig[i] < xmin[i]) {
        gmx_warning("Your xmin[%d] value %f is larger than the smallest data point %f",i,xmin[i],min_eig[i]);
      }
      min_eig[i] = xmin[i];
    }
    else
      min_eig[i] -= delta[i];
    bfac[i]     = ibox[i]/(max_eig[i]-min_eig[i]);
  }
  /* Do the binning */ 
  bref = 1/(BOLTZ*Tref);
  snew(bE,n);
  if (bGE || nenerT==2) {
    Emin = 1e8;
    for(j=0; (j<n); j++) {
      if (bGE)
        bE[j] = bref*enerT[0][j];
      else
        bE[j] = (bref - 1/(BOLTZ*enerT[1][j]))*enerT[0][j];
      Emin  = min(Emin,bE[j]);
    }
  }
  else
    Emin = 0;
  len=1;
  for(i=0; (i<neig); i++) 
    len=len*ibox[i];
  printf("There are %d bins in the %d-dimensional histogram. Beta-Emin = %g\n",
         len,neig,Emin);
  snew(P,len);
  snew(W,len);
  snew(E,len);
  snew(S,len);
  snew(M,len);
  snew(nbin,len);
  snew(bindex,n);

  
  /* Loop over projections */
  for(j=0; (j<n); j++) {
    /* Loop over dimensions */
    bOutside = FALSE;
    for(i=0; (i<neig); i++) {
      nxyz[i] = bfac[i]*(eig[i][j]-min_eig[i]);
      if (nxyz[i] < 0 || nxyz[i] >= ibox[i])
        bOutside = TRUE;
    }
    if (!bOutside) {
      index = indexn(neig,ibox,nxyz); 
      range_check(index,0,len);
      /* Compute the exponential factor */
      if (enerT)
        efac = exp(-bE[j]+Emin);
      else
        efac = 1;
      /* Apply the bin volume correction for a multi-dimensional distance */
      for(i=0; i<neig; i++) {
        if (idim[i] == 2)
          efac /= eig[i][j];
        else if (idim[i] == 3)
          efac /= sqr(eig[i][j]);
        else if (idim[i] == -1)
          efac /= sin(DEG2RAD*eig[i][j]);
      }
      /* Update the probability */
      P[index] += efac;
      /* Update the energy */
      if (enerT)
        E[index] += enerT[0][j];
      /* Statistics: which "structure" in which bin */
      nbin[index]++;
      bindex[j]=index;
    }
  }
  /* Normalize probability */
  normalize_p_e(len,P,nbin,E,pmin);
  Pmax = 0;
  /* Compute boundaries for the Free energy */
  Wmin = 1e8;
  imin = -1;
  Wmax = -1e8;
  /* Recompute Emin: it may have changed due to averaging */
  Emin = 1e8;
  Emax = -1e8;
  for(i=0; (i<len); i++) {
    if (P[i] != 0) {
      Pmax = max(P[i],Pmax);
      W[i] = -BOLTZ*Tref*log(P[i]);
      if (W[i] < Wmin) {
        Wmin = W[i];
        imin = i;
      }
      Emin = min(E[i],Emin);
      Emax = max(E[i],Emax);
      Wmax = max(W[i],Wmax);
    }
  }
  if (pmax > 0) {
    Pmax = pmax;
  }
  if (gmax > 0) {
    Wmax = gmax;
  } else {
    Wmax -= Wmin;
  }
  Winf = Wmax+1;
  Einf = Emax+1;
  Smin = Emin-Wmax;
  Smax = Emax-Smin;
  Sinf = Smax+1;
  /* Write out the free energy as a function of bin index */
  fp = ffopen(fn,"w");
  for(i=0; (i<len); i++)
    if (P[i] != 0) {
      W[i] -= Wmin;
      S[i] = E[i]-W[i]-Smin;
      fprintf(fp,"%5d  %10.5e  %10.5e  %10.5e\n",i,W[i],E[i],S[i]);
    }
    else {
      W[i] = Winf;
      E[i] = Einf;
      S[i] = Sinf;
    }
  ffclose(fp);
  /* Organize the structures in the bins */
  snew(b,1);
  snew(b->index,len+1);
  snew(b->a,n);
  b->index[0] = 0;
  for(i=0; (i<len); i++) {
    b->index[i+1] = b->index[i]+nbin[i]; 
    nbin[i] = 0;
  }
  for(i=0; (i<n); i++) {
    bi = bindex[i];
    b->a[b->index[bi]+nbin[bi]] = i;
    nbin[bi]++;
  }
  /* Consistency check */
  /* This no longer applies when we allow the plot to be smaller
     than the sampled space.
  for(i=0; (i<len); i++) {
    if (nbin[i] != (b->index[i+1] - b->index[i]))
      gmx_fatal(FARGS,"nbin[%d] = %d, should be %d",i,nbin[i],
                b->index[i+1] - b->index[i]);
  }
  */
  /* Write the index file */
  fp = ffopen(ndx,"w");
  for(i=0; (i<len); i++) {
    if (nbin[i] > 0) {
      fprintf(fp,"[ %d ]\n",i);
      for(j=b->index[i]; (j<b->index[i+1]); j++)
        fprintf(fp,"%d\n",b->a[j]+1);
    }
  }  
  ffclose(fp);
  snew(axis_x,ibox[0]+1);
  snew(axis_y,ibox[1]+1);
  snew(axis_z,ibox[2]+1);
  maxbox = max(ibox[0],max(ibox[1],ibox[2]));
  snew(PP,maxbox*maxbox);
  snew(WW,maxbox*maxbox);
  snew(EE,maxbox*maxbox);
  snew(SS,maxbox*maxbox);
  for(i=0; (i<min(neig,3)); i++) {
    switch (i) {
    case 0: axis = axis_x; break;
    case 1: axis = axis_y; break;
    case 2: axis = axis_z; break;
    default: break;
    }
    for(j=0; j<=ibox[i]; j++)
      axis[j] = min_eig[i] + j/bfac[i];
  }
  if (map) {
    snew(M,len);
    snew(MM,maxbox*maxbox);
    for(i=0; (i<ibox[0]); i++) 
      MM[i] = &(M[i*ibox[1]]);
    Mmin = 1e8;
    Mmax = -1e8;
    for(i=0; (i<nmap); i++) {
      Mmin = min(Mmin,map[0][i]);
      Mmax = max(Mmax,map[0][i]);
    }
    Minf = Mmax*1.05;
    for(i=0; (i<len); i++) 
      M[i] = Minf;
    for(i=0; (i<nmap); i++) {
      index = gmx_nint(mapindex[i]);
      if (index >= len)
        gmx_fatal(FARGS,"Number of bins in file from -mdata option does not correspond to current analysis");
      
      if (P[index] != 0)
        M[index] = map[0][i];
    }
  }
  else {
    MM = NULL;
    Minf = NOTSET;
  }
  pick_minima(logf,ibox,neig,len,W);
  if (gmax <= 0)
    gmax = Winf;
  flags = MAT_SPATIAL_X | MAT_SPATIAL_Y;
  if (neig == 2) {
    /* Dump to XPM file */
    snew(PP,ibox[0]);
    for(i=0; (i<ibox[0]); i++) {
      snew(PP[i],ibox[1]);
      for(j=0; j<ibox[1]; j++) {
        PP[i][j] = P[i*ibox[1]+j];
      }
      WW[i] = &(W[i*ibox[1]]);
      EE[i] = &(E[i*ibox[1]]);
      SS[i] = &(S[i*ibox[1]]);
    }
    fp = ffopen(xpmP,"w");
    write_xpm(fp,flags,"Probability Distribution","","PC1","PC2",
              ibox[0],ibox[1],axis_x,axis_y,PP,0,Pmax,rlo,rhi,&nlevels);
    ffclose(fp);
    fp = ffopen(xpm,"w");
    write_xpm(fp,flags,"Gibbs Energy Landscape","G (kJ/mol)","PC1","PC2",
              ibox[0],ibox[1],axis_x,axis_y,WW,0,gmax,rlo,rhi,&nlevels);
    ffclose(fp);
    fp = ffopen(xpm2,"w");
    write_xpm(fp,flags,"Enthalpy Landscape","H (kJ/mol)","PC1","PC2",
              ibox[0],ibox[1],axis_x,axis_y,EE,
              emin ? *emin : Emin,emax ? *emax : Einf,rlo,rhi,&nlevels);
    ffclose(fp);
    fp = ffopen(xpm3,"w");
    write_xpm(fp,flags,"Entropy Landscape","TDS (kJ/mol)","PC1","PC2",
              ibox[0],ibox[1],axis_x,axis_y,SS,0,Sinf,rlo,rhi,&nlevels);
    ffclose(fp);
    if (map) {
      fp = ffopen(xpm4,"w");
      write_xpm(fp,flags,"Custom Landscape",mname,"PC1","PC2",
                ibox[0],ibox[1],axis_x,axis_y,MM,0,Minf,rlo,rhi,&nlevels);
      ffclose(fp);
    }
  }
  else if (neig == 3) {
    /* Dump to PDB file */
    fp = ffopen(pdb,"w");
    for(i=0; (i<ibox[0]); i++) {
      xxx[XX] = 3*(i+0.5-ibox[0]/2);
      for(j=0; (j<ibox[1]); j++) {
        xxx[YY] = 3*(j+0.5-ibox[1]/2);
        for(k=0; (k<ibox[2]); k++) {
          xxx[ZZ] = 3*(k+0.5-ibox[2]/2);
          index = index3(ibox,i,j,k);
          if (P[index] > 0)
            fprintf(fp,"%-6s%5u  %-4.4s%3.3s  %4d    %8.3f%8.3f%8.3f%6.2f%6.2f\n",
                    "ATOM",(index+1) %10000,"H","H",(index+1)%10000,
                    xxx[XX],xxx[YY],xxx[ZZ],1.0,W[index]);
        }
      }
    }
    ffclose(fp);
    write_xplor("out.xplor",W,ibox,min_eig,max_eig);
    if (map)
      write_xplor("user.xplor",M,ibox,min_eig,max_eig);
    nxyz[XX] = imin/(ibox[1]*ibox[2]);
    nxyz[YY] = (imin-nxyz[XX]*ibox[1]*ibox[2])/ibox[2];
    nxyz[ZZ] = imin % ibox[2];
    for(i=0; (i<ibox[0]); i++) {
      snew(WW[i],maxbox);
      for(j=0; (j<ibox[1]); j++)
        WW[i][j] = W[index3(ibox,i,j,nxyz[ZZ])];
    }
    snew(buf,strlen(xpm)+4);
    sprintf(buf,"%s",xpm);
    sprintf(&buf[strlen(xpm)-4],"12.xpm");
    fp = ffopen(buf,"w");
    write_xpm(fp,flags,"Gibbs Energy Landscape","W (kJ/mol)","PC1","PC2",
              ibox[0],ibox[1],axis_x,axis_y,WW,0,gmax,rlo,rhi,&nlevels);
    ffclose(fp);
    for(i=0; (i<ibox[0]); i++) {
      for(j=0; (j<ibox[2]); j++)
        WW[i][j] = W[index3(ibox,i,nxyz[YY],j)];
    }
    sprintf(&buf[strlen(xpm)-4],"13.xpm");
    fp = ffopen(buf,"w");
    write_xpm(fp,flags,"SHAM Energy Landscape","kJ/mol","PC1","PC3",
              ibox[0],ibox[2],axis_x,axis_z,WW,0,gmax,rlo,rhi,&nlevels);
    ffclose(fp);
    for(i=0; (i<ibox[1]); i++) {
      for(j=0; (j<ibox[2]); j++)
        WW[i][j] = W[index3(ibox,nxyz[XX],i,j)];
    }
    sprintf(&buf[strlen(xpm)-4],"23.xpm");
    fp = ffopen(buf,"w");
    write_xpm(fp,flags,"SHAM Energy Landscape","kJ/mol","PC2","PC3",
              ibox[1],ibox[2],axis_y,axis_z,WW,0,gmax,rlo,rhi,&nlevels);
    ffclose(fp);
    sfree(buf);
  }
  if (map) {
    sfree(MM);
    sfree(M);
  }
}

static void ehisto(const char *fh,int n,real **enerT, const output_env_t oenv)
{
  FILE *fp;
  int  i,j,k,nbin,blength;
  int  *bindex;
  real *T,bmin,bmax,bwidth;
  int  **histo;
  
  bmin =  1e8;
  bmax = -1e8;
  snew(bindex,n);
  snew(T,n);
  nbin = 0;
  for(j=1; (j<n); j++) {
    for(k=0; (k<nbin); k++) {
      if (T[k] == enerT[1][j]) {
        bindex[j] = k;
        break;
      }
    }
    if (k == nbin) {
      bindex[j] = nbin;
      T[nbin]   = enerT[1][j];
      nbin++;
    }
    bmin = min(enerT[0][j],bmin);
    bmax = max(enerT[0][j],bmax);
}
  bwidth  = 1.0;
  blength = (bmax - bmin)/bwidth + 2;
  snew(histo,nbin);
  for(i=0; (i<nbin); i++) {
    snew(histo[i],blength);
  }
  for(j=0; (j<n); j++) {
    k = (enerT[0][j]-bmin)/bwidth;
    histo[bindex[j]][k]++;
  }
  fp = xvgropen(fh,"Energy distribution","E (kJ/mol)","",oenv);
  for(j=0; (j<blength); j++) {
    fprintf(fp,"%8.3f",bmin+j*bwidth);
    for(k=0; (k<nbin); k++) {
      fprintf(fp,"  %6d",histo[k][j]);
    }
    fprintf(fp,"\n");
  }  
  ffclose(fp);
}

int gmx_sham(int argc,char *argv[])
{
  const char *desc[] = {
    "[TT]g_sham[tt] makes multi-dimensional free-energy, enthalpy and entropy plots.",
    "[TT]g_sham[tt] reads one or more [TT].xvg[tt] files and analyzes data sets.",
    "The basic purpose of [TT]g_sham[tt] is to plot Gibbs free energy landscapes",
    "(option [TT]-ls[tt])",
    "by Bolzmann inverting multi-dimensional histograms (option [TT]-lp[tt]),",
    "but it can also",
    "make enthalpy (option [TT]-lsh[tt]) and entropy (option [TT]-lss[tt])",
    "plots. The histograms can be made for any quantities the user supplies.",
    "A line in the input file may start with a time",
    "(see option [TT]-time[tt]) and any number of [IT]y[it]-values may follow.",
    "Multiple sets can also be",
    "read when they are separated by & (option [TT]-n[tt]),",
    "in this case only one [IT]y[it]-value is read from each line.",
    "All lines starting with # and @ are skipped.",
    "[PAR]",
    "Option [TT]-ge[tt] can be used to supply a file with free energies",
    "when the ensemble is not a Boltzmann ensemble, but needs to be biased",
    "by this free energy. One free energy value is required for each",
    "(multi-dimensional) data point in the [TT]-f[tt] input.",
    "[PAR]",
    "Option [TT]-ene[tt] can be used to supply a file with energies.",
    "These energies are used as a weighting function in the single",
    "histogram analysis method by Kumar et al. When temperatures",
    "are supplied (as a second column in the file), an experimental",
    "weighting scheme is applied. In addition the vales",
    "are used for making enthalpy and entropy plots.",
    "[PAR]",
    "With option [TT]-dim[tt], dimensions can be gives for distances.",
    "When a distance is 2- or 3-dimensional, the circumference or surface",
    "sampled by two particles increases with increasing distance.",
    "Depending on what one would like to show, one can choose to correct",
    "the histogram and free-energy for this volume effect.",
    "The probability is normalized by r and r^2 for dimensions of 2 and 3, ",
    "respectively.",
    "A value of -1 is used to indicate an angle in degrees between two",
    "vectors: a sin(angle) normalization will be applied.",
    "[BB]Note[bb] that for angles between vectors the inner-product or cosine",
    "is the natural quantity to use, as it will produce bins of the same",
    "volume."
  };
  static real tb=-1,te=-1,frac=0.5,filtlen=0,binwidth=0.1;
  static gmx_bool bHaveT=TRUE,bDer=FALSE,bSubAv=TRUE,bAverCorr=FALSE,bXYdy=FALSE;
  static gmx_bool bEESEF=FALSE,bEENLC=FALSE,bEeFitAc=FALSE,bPower=FALSE;
  static gmx_bool bShamEner=TRUE,bSham=TRUE; 
  static real Tref=298.15,pmin=0,ttol=0,pmax=0,gmax=0,emin=0,emax=0;
  static rvec nrdim = {1,1,1};
  static rvec nrbox = {32,32,32};
  static rvec xmin  = {0,0,0}, xmax={1,1,1};
  static int  nsets_in=1,nb_min=4,resol=10,nlevels=25;
  static const char *mname="";
  t_pargs pa[] = {
    { "-time",    FALSE, etBOOL, {&bHaveT},
      "Expect a time in the input" },
    { "-b",       FALSE, etREAL, {&tb},
      "First time to read from set" },
    { "-e",       FALSE, etREAL, {&te},
      "Last time to read from set" },
    { "-ttol",     FALSE, etREAL, {&ttol},
      "Tolerance on time in appropriate units (usually ps)" },
    { "-n",       FALSE, etINT, {&nsets_in},
      "Read this number of sets separated by lines containing only an ampersand" },
    { "-d",       FALSE, etBOOL, {&bDer},
        "Use the derivative" },
    { "-bw",      FALSE, etREAL, {&binwidth},
      "Binwidth for the distribution" },
    { "-sham",    FALSE, etBOOL, {&bSham},
      "Turn off energy weighting even if energies are given" },
    { "-tsham",   FALSE, etREAL, {&Tref},
      "Temperature for single histogram analysis" },
    { "-pmin",    FALSE, etREAL, {&pmin},
      "Minimum probability. Anything lower than this will be set to zero" },
    { "-dim",     FALSE, etRVEC, {nrdim},
      "Dimensions for distances, used for volume correction (max 3 values, dimensions > 3 will get the same value as the last)" },
    { "-ngrid",   FALSE, etRVEC, {nrbox},   
      "Number of bins for energy landscapes (max 3 values, dimensions > 3 will get the same value as the last)" },
    { "-xmin",    FALSE, etRVEC, {xmin},
      "Minimum for the axes in energy landscape (see above for > 3 dimensions)" },
    { "-xmax",    FALSE, etRVEC, {xmax},
      "Maximum for the axes in energy landscape (see above for > 3 dimensions)" },
    { "-pmax",    FALSE, etREAL, {&pmax},
      "Maximum probability in output, default is calculate" },
    { "-gmax",    FALSE, etREAL, {&gmax},
      "Maximum free energy in output, default is calculate" },
    { "-emin",    FALSE, etREAL, {&emin},
      "Minimum enthalpy in output, default is calculate" },
    { "-emax",    FALSE, etREAL, {&emax},
      "Maximum enthalpy in output, default is calculate" },
    { "-nlevels", FALSE, etINT,  {&nlevels},
      "Number of levels for energy landscape" },
    { "-mname",   FALSE, etSTR,  {&mname},
      "Legend label for the custom landscape" },
  };
#define NPA asize(pa)

  FILE     *out;
  int      n,e_n,d_n,nlast,s,nset,e_nset,d_nset,i,j=0,*idim,*ibox;
  real     **val,**et_val,**dt_val,*t,*e_t,e_dt,d_dt,*d_t,dt,tot,error;
  real     *rmin,*rmax;
  double   *av,*sig,cum1,cum2,cum3,cum4,db;
  const char     *fn_ge,*fn_ene;
  output_env_t oenv;
  gmx_large_int_t num_grid_points;
    
  t_filenm fnm[] = { 
    { efXVG, "-f",    "graph",    ffREAD   },
    { efXVG, "-ge",   "gibbs",    ffOPTRD  },
    { efXVG, "-ene",  "esham",    ffOPTRD  },
    { efXVG, "-dist", "ener",     ffOPTWR  },
    { efXVG, "-histo","edist",    ffOPTWR  },
    { efNDX, "-bin",  "bindex",   ffOPTWR  },
    { efXPM, "-lp",   "prob",     ffOPTWR  },
    { efXPM, "-ls",   "gibbs",    ffOPTWR  },
    { efXPM, "-lsh",  "enthalpy", ffOPTWR  },
    { efXPM, "-lss",  "entropy",  ffOPTWR  },
    { efXPM, "-map",  "map",      ffOPTWR  },
    { efPDB, "-ls3",  "gibbs3",   ffOPTWR  },
    { efXVG, "-mdata","mapdata",  ffOPTRD  },
    { efLOG, "-g",    "shamlog",  ffOPTWR  }
  }; 
#define NFILE asize(fnm) 

  int     npargs;

  npargs = asize(pa); 
  
  CopyRight(stderr,argv[0]); 
  parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_BE_NICE ,
                    NFILE,fnm,npargs,pa,asize(desc),desc,0,NULL,&oenv); 

  val=read_xvg_time(opt2fn("-f",NFILE,fnm),bHaveT,
                    opt2parg_bSet("-b",npargs,pa),tb-ttol,
                    opt2parg_bSet("-e",npargs,pa),te+ttol,
                    nsets_in,&nset,&n,&dt,&t);
  printf("Read %d sets of %d points, dt = %g\n\n",nset,n,dt);
  
  fn_ge  = opt2fn_null("-ge",NFILE,fnm);
  fn_ene = opt2fn_null("-ene",NFILE,fnm);

  if (fn_ge && fn_ene)
    gmx_fatal(FARGS,"Can not do free energy and energy corrections at the same time");

  if (fn_ge || fn_ene) {
    et_val=read_xvg_time(fn_ge ? fn_ge : fn_ene,bHaveT,
                         opt2parg_bSet("-b",npargs,pa),tb-ttol,
                         opt2parg_bSet("-e",npargs,pa),te+ttol,
                         1,&e_nset,&e_n,&e_dt,&e_t);
    if (fn_ge) {
      if (e_nset != 1)
        gmx_fatal(FARGS,"Can only handle one free energy component in %s",
                fn_ge);
    } else {
      if (e_nset!=1 && e_nset!=2)
        gmx_fatal(FARGS,"Can only handle one energy component or one energy and one T in %s",
                  fn_ene);
    }
    if (e_n != n)
      gmx_fatal(FARGS,"Number of energies (%d) does not match number of entries (%d) in %s",e_n,n,opt2fn("-f",NFILE,fnm));
  }
  else 
    et_val = NULL;
    
  if (opt2fn_null("-mdata",NFILE,fnm) != NULL) {
    dt_val=read_xvg_time(opt2fn("-mdata",NFILE,fnm),bHaveT,
                         FALSE,tb,FALSE,te,
                         nsets_in,&d_nset,&d_n,&d_dt,&d_t);
    if (d_nset != 1)
      gmx_fatal(FARGS,"Can only handle one mapping data column in %s",
                opt2fn("-mdata",NFILE,fnm));
  }
  else
    dt_val = NULL;

  if (fn_ene && et_val)
    ehisto(opt2fn("-histo",NFILE,fnm),e_n,et_val,oenv);

  snew(idim,max(3,nset));
  snew(ibox,max(3,nset));
  snew(rmin,max(3,nset));
  snew(rmax,max(3,nset));
  for(i=0; (i<min(3,nset)); i++) {
    idim[i] = nrdim[i];
    ibox[i] = nrbox[i];
    rmin[i] = xmin[i];
    rmax[i] = xmax[i];
  }
  for(; (i<nset); i++) {
    idim[i] = nrdim[2];
    ibox[i] = nrbox[2];
    rmin[i] = xmin[2];
    rmax[i] = xmax[2];
  }

  /* Check that the grid size is manageable. */
  num_grid_points = ibox[0];
  for(i = 1; i < nset; i++)
  {
      gmx_large_int_t result;
      if (!check_int_multiply_for_overflow(num_grid_points, ibox[i], &result))
      {
          gmx_fatal(FARGS,
                    "The number of dimensions and grid points is too large for this tool\n"
                    "to handle with what it knows about the architecture upon which it\n"
                    "is running. Use a different machine or consult the GROMACS mailing list.");
      }
      num_grid_points = result;
  }
  /* The number of grid points fits in a gmx_large_int_t. */

  do_sham(opt2fn("-dist",NFILE,fnm),opt2fn("-bin",NFILE,fnm),
          opt2fn("-lp",NFILE,fnm),
          opt2fn("-ls",NFILE,fnm),opt2fn("-lsh",NFILE,fnm),
          opt2fn("-lss",NFILE,fnm),opt2fn("-map",NFILE,fnm),
          opt2fn("-ls3",NFILE,fnm),opt2fn("-g",NFILE,fnm),
          n,nset,val,fn_ge!=NULL,e_nset,et_val,d_n,d_t,dt_val,Tref,
          pmax,gmax,
          opt2parg_bSet("-emin",NPA,pa) ? &emin : NULL,
          opt2parg_bSet("-emax",NPA,pa) ? &emax : NULL,
          nlevels,pmin,
          mname,bSham,idim,ibox,
          opt2parg_bSet("-xmin",NPA,pa),rmin,
          opt2parg_bSet("-xmax",NPA,pa),rmax);
  
  thanx(stderr);

  return 0;
}