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