src/gromacs/gmxana/gmx_sham.c

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