src/gromacs/gmxana/gmx_sham.cpp

   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,2015,2016,2017,2018, 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 <cassert>
  40 #include <cmath>
  41 #include <cstdlib>
  42 #include <cstring>
  43
  44 #include <algorithm>
  45
  46 #include "gromacs/commandline/pargs.h"
  47 #include "gromacs/fileio/matio.h"
  48 #include "gromacs/fileio/pdbio.h"
  49 #include "gromacs/fileio/xvgr.h"
  50 #include "gromacs/gmxana/gmx_ana.h"
  51 #include "gromacs/gmxana/gstat.h"
  52 #include "gromacs/math/functions.h"
  53 #include "gromacs/math/units.h"
  54 #include "gromacs/math/vec.h"
  55 #include "gromacs/topology/topology.h"
  56 #include "gromacs/utility/arraysize.h"
  57 #include "gromacs/utility/fatalerror.h"
  58 #include "gromacs/utility/futil.h"
  59 #include "gromacs/utility/smalloc.h"
  60
  61
  62 static int index2(const int *ibox, int x, int y)
  63 {
  64     return (ibox[1]*x+y);
  65 }
  66
  67 static int index3(const 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, const real *data, int *ibox, const real dmin[], const 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, const 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 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 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 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                                      const 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 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                                      const 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                     const real *emin, const real *emax, int nlevels, real pmin,
 423                     const 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 = nullptr;
 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;
 434     real        *delta;
 435     int          i, j, k, imin, len, index, *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] = std::min(min_eig[i], eig[i][j]);
 457             max_eig[i] = std::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         assert(enerT);
 496         Emin = 1e8;
 497         for (j = 0; (j < n); j++)
 498         {
 499             if (bGE)
 500             {
 501                 bE[j] = bref*enerT[0][j];
 502             }
 503             else
 504             {
 505                 bE[j] = (bref - 1/(BOLTZ*enerT[1][j]))*enerT[0][j];
 506             }
 507             Emin  = std::min(Emin, static_cast<double>(bE[j]));
 508         }
 509     }
 510     else
 511     {
 512         Emin = 0;
 513     }
 514     len = 1;
 515     for (i = 0; (i < neig); i++)
 516     {
 517         len = len*ibox[i];
 518     }
 519     printf("There are %d bins in the %d-dimensional histogram. Beta-Emin = %g\n",
 520            len, neig, Emin);
 521     snew(P, len);
 522     snew(W, len);
 523     snew(E, len);
 524     snew(S, len);
 525     snew(M, len);
 526     snew(nbin, len);
 527     snew(bindex, n);
 528
 529
 530     /* Loop over projections */
 531     for (j = 0; (j < n); j++)
 532     {
 533         /* Loop over dimensions */
 534         bOutside = FALSE;
 535         for (i = 0; (i < neig); i++)
 536         {
 537             nxyz[i] = static_cast<int>(bfac[i]*(eig[i][j]-min_eig[i]));
 538             if (nxyz[i] < 0 || nxyz[i] >= ibox[i])
 539             {
 540                 bOutside = TRUE;
 541             }
 542         }
 543         if (!bOutside)
 544         {
 545             index = indexn(neig, ibox, nxyz);
 546             range_check(index, 0, len);
 547             /* Compute the exponential factor */
 548             if (enerT)
 549             {
 550                 efac = std::exp(-bE[j]+Emin);
 551             }
 552             else
 553             {
 554                 efac = 1;
 555             }
 556             /* Apply the bin volume correction for a multi-dimensional distance */
 557             for (i = 0; i < neig; i++)
 558             {
 559                 if (idim[i] == 2)
 560                 {
 561                     efac /= eig[i][j];
 562                 }
 563                 else if (idim[i] == 3)
 564                 {
 565                     efac /= gmx::square(eig[i][j]);
 566                 }
 567                 else if (idim[i] == -1)
 568                 {
 569                     efac /= std::sin(DEG2RAD*eig[i][j]);
 570                 }
 571             }
 572             /* Update the probability */
 573             P[index] += efac;
 574             /* Update the energy */
 575             if (enerT)
 576             {
 577                 E[index] += enerT[0][j];
 578             }
 579             /* Statistics: which "structure" in which bin */
 580             nbin[index]++;
 581             bindex[j] = index;
 582         }
 583     }
 584     /* Normalize probability */
 585     normalize_p_e(len, P, nbin, E, pmin);
 586     Pmax = 0;
 587     /* Compute boundaries for the Free energy */
 588     Wmin = 1e8;
 589     imin = -1;
 590     Wmax = -1e8;
 591     /* Recompute Emin: it may have changed due to averaging */
 592     Emin = 1e8;
 593     Emax = -1e8;
 594     for (i = 0; (i < len); i++)
 595     {
 596         if (P[i] != 0)
 597         {
 598             Pmax = std::max(P[i], Pmax);
 599             W[i] = -BOLTZ*Tref*std::log(P[i]);
 600             if (W[i] < Wmin)
 601             {
 602                 Wmin = W[i];
 603                 imin = i;
 604             }
 605             Emin = std::min(static_cast<double>(E[i]), Emin);
 606             Emax = std::max(static_cast<double>(E[i]), Emax);
 607             Wmax = std::max(static_cast<double>(W[i]), Wmax);
 608         }
 609     }
 610     if (pmax > 0)
 611     {
 612         Pmax = pmax;
 613     }
 614     if (gmax > 0)
 615     {
 616         Wmax = gmax;
 617     }
 618     else
 619     {
 620         Wmax -= Wmin;
 621     }
 622     Winf = Wmax+1;
 623     Einf = Emax+1;
 624     Smin = Emin-Wmax;
 625     Smax = Emax-Smin;
 626     Sinf = Smax+1;
 627     /* Write out the free energy as a function of bin index */
 628     fp = gmx_ffopen(fn, "w");
 629     for (i = 0; (i < len); i++)
 630     {
 631         if (P[i] != 0)
 632         {
 633             W[i] -= Wmin;
 634             S[i]  = E[i]-W[i]-Smin;
 635             fprintf(fp, "%5d  %10.5e  %10.5e  %10.5e\n", i, W[i], E[i], S[i]);
 636         }
 637         else
 638         {
 639             W[i] = Winf;
 640             E[i] = Einf;
 641             S[i] = Sinf;
 642         }
 643     }
 644     gmx_ffclose(fp);
 645     /* Organize the structures in the bins */
 646     snew(b, 1);
 647     snew(b->index, len+1);
 648     snew(b->a, n);
 649     b->index[0] = 0;
 650     for (i = 0; (i < len); i++)
 651     {
 652         b->index[i+1] = b->index[i]+nbin[i];
 653         nbin[i]       = 0;
 654     }
 655     for (i = 0; (i < n); i++)
 656     {
 657         bi = bindex[i];
 658         b->a[b->index[bi]+nbin[bi]] = i;
 659         nbin[bi]++;
 660     }
 661     /* Consistency check */
 662     /* This no longer applies when we allow the plot to be smaller
 663        than the sampled space.
 664        for(i=0; (i<len); i++) {
 665        if (nbin[i] != (b->index[i+1] - b->index[i]))
 666         gmx_fatal(FARGS,"nbin[%d] = %d, should be %d",i,nbin[i],
 667           b->index[i+1] - b->index[i]);
 668        }
 669      */
 670     /* Write the index file */
 671     fp = gmx_ffopen(ndx, "w");
 672     for (i = 0; (i < len); i++)
 673     {
 674         if (nbin[i] > 0)
 675         {
 676             fprintf(fp, "[ %d ]\n", i);
 677             for (j = b->index[i]; (j < b->index[i+1]); j++)
 678             {
 679                 fprintf(fp, "%d\n", b->a[j]+1);
 680             }
 681         }
 682     }
 683     gmx_ffclose(fp);
 684     snew(axis_x, ibox[0]+1);
 685     snew(axis_y, ibox[1]+1);
 686     snew(axis_z, ibox[2]+1);
 687     maxbox = std::max(ibox[0], std::max(ibox[1], ibox[2]));
 688     snew(PP, maxbox*maxbox);
 689     snew(WW, maxbox*maxbox);
 690     snew(EE, maxbox*maxbox);
 691     snew(SS, maxbox*maxbox);
 692     for (i = 0; (i < std::min(neig, 3)); i++)
 693     {
 694         switch (i)
 695         {
 696             case 0: axis = axis_x; break;
 697             case 1: axis = axis_y; break;
 698             case 2: axis = axis_z; break;
 699             default: break;
 700         }
 701         for (j = 0; j <= ibox[i]; j++)
 702         {
 703             axis[j] = min_eig[i] + j/bfac[i];
 704         }
 705     }
 706
 707     pick_minima(logf, ibox, neig, len, W);
 708     if (gmax <= 0)
 709     {
 710         gmax = Winf;
 711     }
 712     flags = MAT_SPATIAL_X | MAT_SPATIAL_Y;
 713     if (neig == 2)
 714     {
 715         /* Dump to XPM file */
 716         snew(PP, ibox[0]);
 717         for (i = 0; (i < ibox[0]); i++)
 718         {
 719             snew(PP[i], ibox[1]);
 720             for (j = 0; j < ibox[1]; j++)
 721             {
 722                 PP[i][j] = P[i*ibox[1]+j];
 723             }
 724             WW[i] = &(W[i*ibox[1]]);
 725             EE[i] = &(E[i*ibox[1]]);
 726             SS[i] = &(S[i*ibox[1]]);
 727         }
 728         fp = gmx_ffopen(xpmP, "w");
 729         write_xpm(fp, flags, "Probability Distribution", "", "PC1", "PC2",
 730                   ibox[0], ibox[1], axis_x, axis_y, PP, 0, Pmax, rlo, rhi, &nlevels);
 731         gmx_ffclose(fp);
 732         fp = gmx_ffopen(xpm, "w");
 733         write_xpm(fp, flags, "Gibbs Energy Landscape", "G (kJ/mol)", "PC1", "PC2",
 734                   ibox[0], ibox[1], axis_x, axis_y, WW, 0, gmax, rlo, rhi, &nlevels);
 735         gmx_ffclose(fp);
 736         fp = gmx_ffopen(xpm2, "w");
 737         write_xpm(fp, flags, "Enthalpy Landscape", "H (kJ/mol)", "PC1", "PC2",
 738                   ibox[0], ibox[1], axis_x, axis_y, EE,
 739                   emin ? *emin : Emin, emax ? *emax : Einf, rlo, rhi, &nlevels);
 740         gmx_ffclose(fp);
 741         fp = gmx_ffopen(xpm3, "w");
 742         write_xpm(fp, flags, "Entropy Landscape", "TDS (kJ/mol)", "PC1", "PC2",
 743                   ibox[0], ibox[1], axis_x, axis_y, SS, 0, Sinf, rlo, rhi, &nlevels);
 744         gmx_ffclose(fp);
 745     }
 746     else if (neig == 3)
 747     {
 748         /* Dump to PDB file */
 749         fp = gmx_ffopen(pdb, "w");
 750         for (i = 0; (i < ibox[0]); i++)
 751         {
 752             xxx[XX] = 3*i+1.5*(1-ibox[0]);
 753             for (j = 0; (j < ibox[1]); j++)
 754             {
 755                 xxx[YY] = 3*j+1.5*(1-ibox[1]);
 756                 for (k = 0; (k < ibox[2]); k++)
 757                 {
 758                     xxx[ZZ] = 3*k+1.5*(1-ibox[2]);
 759                     index   = index3(ibox, i, j, k);
 760                     if (P[index] > 0)
 761                     {
 762                         fprintf(fp, "%-6s%5d  %-4.4s%3.3s  %4d    %8.3f%8.3f%8.3f%6.2f%6.2f\n",
 763                                 "ATOM", (index+1) %10000, "H", "H", (index+1)%10000,
 764                                 xxx[XX], xxx[YY], xxx[ZZ], 1.0, W[index]);
 765                     }
 766                 }
 767             }
 768         }
 769         gmx_ffclose(fp);
 770         write_xplor("out.xplor", W, ibox, min_eig, max_eig);
 771         nxyz[XX] = imin/(ibox[1]*ibox[2]);
 772         nxyz[YY] = (imin-nxyz[XX]*ibox[1]*ibox[2])/ibox[2];
 773         nxyz[ZZ] = imin % ibox[2];
 774         for (i = 0; (i < ibox[0]); i++)
 775         {
 776             snew(WW[i], maxbox);
 777             for (j = 0; (j < ibox[1]); j++)
 778             {
 779                 WW[i][j] = W[index3(ibox, i, j, nxyz[ZZ])];
 780             }
 781         }
 782         snew(buf, std::strlen(xpm)+4);
 783         sprintf(buf, "%s", xpm);
 784         sprintf(&buf[std::strlen(xpm)-4], "12.xpm");
 785         fp = gmx_ffopen(buf, "w");
 786         write_xpm(fp, flags, "Gibbs Energy Landscape", "W (kJ/mol)", "PC1", "PC2",
 787                   ibox[0], ibox[1], axis_x, axis_y, WW, 0, gmax, rlo, rhi, &nlevels);
 788         gmx_ffclose(fp);
 789         for (i = 0; (i < ibox[0]); i++)
 790         {
 791             for (j = 0; (j < ibox[2]); j++)
 792             {
 793                 WW[i][j] = W[index3(ibox, i, nxyz[YY], j)];
 794             }
 795         }
 796         sprintf(&buf[std::strlen(xpm)-4], "13.xpm");
 797         fp = gmx_ffopen(buf, "w");
 798         write_xpm(fp, flags, "SHAM Energy Landscape", "kJ/mol", "PC1", "PC3",
 799                   ibox[0], ibox[2], axis_x, axis_z, WW, 0, gmax, rlo, rhi, &nlevels);
 800         gmx_ffclose(fp);
 801         for (i = 0; (i < ibox[1]); i++)
 802         {
 803             for (j = 0; (j < ibox[2]); j++)
 804             {
 805                 WW[i][j] = W[index3(ibox, nxyz[XX], i, j)];
 806             }
 807         }
 808         sprintf(&buf[std::strlen(xpm)-4], "23.xpm");
 809         fp = gmx_ffopen(buf, "w");
 810         write_xpm(fp, flags, "SHAM Energy Landscape", "kJ/mol", "PC2", "PC3",
 811                   ibox[1], ibox[2], axis_y, axis_z, WW, 0, gmax, rlo, rhi, &nlevels);
 812         gmx_ffclose(fp);
 813         sfree(buf);
 814     }
 815 }
 816
 817 static void ehisto(const char *fh, int n, real **enerT, const gmx_output_env_t *oenv)
 818 {
 819     FILE  *fp;
 820     int    i, j, k, nbin, blength;
 821     int   *bindex;
 822     real  *T, bmin, bmax, bwidth;
 823     int  **histo;
 824
 825     bmin =  1e8;
 826     bmax = -1e8;
 827     snew(bindex, n);
 828     snew(T, n);
 829     nbin = 0;
 830     for (j = 1; (j < n); j++)
 831     {
 832         for (k = 0; (k < nbin); k++)
 833         {
 834             if (T[k] == enerT[1][j])
 835             {
 836                 bindex[j] = k;
 837                 break;
 838             }
 839         }
 840         if (k == nbin)
 841         {
 842             bindex[j] = nbin;
 843             T[nbin]   = enerT[1][j];
 844             nbin++;
 845         }
 846         bmin = std::min(enerT[0][j], bmin);
 847         bmax = std::max(enerT[0][j], bmax);
 848     }
 849     bwidth  = 1.0;
 850     blength = static_cast<int>((bmax - bmin)/bwidth + 2);
 851     snew(histo, nbin);
 852     for (i = 0; (i < nbin); i++)
 853     {
 854         snew(histo[i], blength);
 855     }
 856     for (j = 0; (j < n); j++)
 857     {
 858         k = static_cast<int>((enerT[0][j]-bmin)/bwidth);
 859         histo[bindex[j]][k]++;
 860     }
 861     fp = xvgropen(fh, "Energy distribution", "E (kJ/mol)", "", oenv);
 862     for (j = 0; (j < blength); j++)
 863     {
 864         fprintf(fp, "%8.3f", bmin+j*bwidth);
 865         for (k = 0; (k < nbin); k++)
 866         {
 867             fprintf(fp, "  %6d", histo[k][j]);
 868         }
 869         fprintf(fp, "\n");
 870     }
 871     xvgrclose(fp);
 872 }
 873
 874 int gmx_sham(int argc, char *argv[])
 875 {
 876     const char        *desc[] = {
 877         "[THISMODULE] makes multi-dimensional free-energy, enthalpy and entropy plots.",
 878         "[THISMODULE] reads one or more [REF].xvg[ref] files and analyzes data sets.",
 879         "The basic purpose of [THISMODULE] is to plot Gibbs free energy landscapes",
 880         "(option [TT]-ls[tt])",
 881         "by Bolzmann inverting multi-dimensional histograms (option [TT]-lp[tt]),",
 882         "but it can also",
 883         "make enthalpy (option [TT]-lsh[tt]) and entropy (option [TT]-lss[tt])",
 884         "plots. The histograms can be made for any quantities the user supplies.",
 885         "A line in the input file may start with a time",
 886         "(see option [TT]-time[tt]) and any number of [IT]y[it]-values may follow.",
 887         "Multiple sets can also be",
 888         "read when they are separated by & (option [TT]-n[tt]),",
 889         "in this case only one [IT]y[it]-value is read from each line.",
 890         "All lines starting with # and @ are skipped.",
 891         "[PAR]",
 892         "Option [TT]-ge[tt] can be used to supply a file with free energies",
 893         "when the ensemble is not a Boltzmann ensemble, but needs to be biased",
 894         "by this free energy. One free energy value is required for each",
 895         "(multi-dimensional) data point in the [TT]-f[tt] input.",
 896         "[PAR]",
 897         "Option [TT]-ene[tt] can be used to supply a file with energies.",
 898         "These energies are used as a weighting function in the single",
 899         "histogram analysis method by Kumar et al. When temperatures",
 900         "are supplied (as a second column in the file), an experimental",
 901         "weighting scheme is applied. In addition the vales",
 902         "are used for making enthalpy and entropy plots.",
 903         "[PAR]",
 904         "With option [TT]-dim[tt], dimensions can be gives for distances.",
 905         "When a distance is 2- or 3-dimensional, the circumference or surface",
 906         "sampled by two particles increases with increasing distance.",
 907         "Depending on what one would like to show, one can choose to correct",
 908         "the histogram and free-energy for this volume effect.",
 909         "The probability is normalized by r and r^2 for dimensions of 2 and 3, ",
 910         "respectively.",
 911         "A value of -1 is used to indicate an angle in degrees between two",
 912         "vectors: a sin(angle) normalization will be applied.",
 913         "[BB]Note[bb] that for angles between vectors the inner-product or cosine",
 914         "is the natural quantity to use, as it will produce bins of the same",
 915         "volume."
 916     };
 917     static real        tb        = -1, te = -1;
 918     static gmx_bool    bHaveT    = TRUE, bDer = FALSE;
 919     static gmx_bool    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, 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     int               n, e_n, nset, e_nset = 0, i, *idim, *ibox;
 966     real            **val, **et_val, *t, *e_t, e_dt, dt;
 967     real             *rmin, *rmax;
 968     const char       *fn_ge, *fn_ene;
 969     gmx_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, nullptr, &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 = nullptr;
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, std::max(3, nset));
1049     snew(ibox, std::max(3, nset));
1050     snew(rmin, std::max(3, nset));
1051     snew(rmax, std::max(3, nset));
1052     for (i = 0; (i < std::min(3, nset)); i++)
1053     {
1054         idim[i] = static_cast<int>(nrdim[i]);
1055         ibox[i] = static_cast<int>(nrbox[i]);
1056         rmin[i] = xmin[i];
1057         rmax[i] = xmax[i];
1058     }
1059     for (; (i < nset); i++)
1060     {
1061         idim[i] = static_cast<int>(nrdim[2]);
1062         ibox[i] = static_cast<int>(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 != nullptr, e_nset, et_val, Tref,
1087             pmax, gmax,
1088             opt2parg_bSet("-emin", NPA, pa) ? &emin : nullptr,
1089             opt2parg_bSet("-emax", NPA, pa) ? &emax : nullptr,
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 }