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