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