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