src/gromacs/gmxana/gmx_sham.c

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