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