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