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