src/tools/gmx_wham.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
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40
  41 #include <stdio.h>
  42
  43 #include "statutil.h"
  44 #include "typedefs.h"
  45 #include "smalloc.h"
  46 #include "vec.h"
  47 #include "copyrite.h"
  48 #include "statutil.h"
  49 #include "tpxio.h"
  50 #include "names.h"
  51 #include "gmx_ana.h"
  52
  53
  54 #ifndef HAVE_STRDUP
  55 #define HAVE_STRDUP
  56 #endif
  57 #include "string2.h"
  58 #include "xvgr.h"
  59
  60
  61 /* enum for energy units */
  62 enum { enSel, en_kJ, en_kCal, en_kT, enNr };
  63 /* enum for type of input files (pdos, tpr, or pullf) */
  64 enum { whamin_unknown, whamin_tpr, whamin_pullxf, whamin_pdo };
  65 /* enum for methods to make profile cyclic/periodic */
  66 enum { enCycl, enCycl_no, enCycl_yes, enCycl_weighted, enCycl_nr};
  67
  68 typedef struct
  69 {
  70     /* umbrella with pull code of gromacs 4 */
  71     int npullgrps;      /* nr of pull groups in tpr file         */
  72     int pull_geometry;  /* such as distance, position            */
  73     ivec pull_dim;      /* pull dimension with geometry distance */
  74     int  pull_ndim;     /* nr of pull_dim != 0                   */
  75     real *k;            /* force constants in tpr file           */
  76     rvec *init_dist;    /* reference displacements               */
  77     real *umbInitDist;  /* referebce displacement in umbrella direction */
  78
  79     /* From here, old pdo stuff */
  80     int nSkip;
  81     char Reference[256];
  82     int nPull;
  83     int nDim;
  84     ivec Dims;
  85     char PullName[4][256];
  86     double UmbPos[4][3];
  87     double UmbCons[4][3];
  88     gmx_bool Flipped[4];
  89 } t_UmbrellaHeader;
  90
  91 typedef struct
  92 {
  93     int nPull;
  94     int nBin;
  95     double **Histo,**cum;
  96     double *k;
  97     double *pos;
  98     double *z;
  99     double * N, *Ntot;
 100     gmx_bool * Flipped;
 101     double dt;
 102     gmx_bool **bContrib;
 103 } t_UmbrellaWindow;
 104
 105 typedef struct
 106 {
 107     const char *fnTpr,*fnPullf,*fnPdo,*fnPullx;
 108     gmx_bool bTpr,bPullf,bPdo,bPullx;
 109     int bins,cycl;
 110     gmx_bool verbose,bShift,bAuto,bBoundsOnly;
 111     gmx_bool bFlipProf;
 112     real tmin, tmax, dt;
 113     real Temperature,Tolerance;
 114     int nBootStrap,histBootStrapBlockLength;
 115     real dtBootStrap,zProfZero,alpha;
 116     int bsSeed,stepchange;
 117     gmx_bool bHistBootStrap,bWeightedCycl,bHistOutOnly;
 118     gmx_bool bAutobounds,bNoprof;
 119     real min,max,dz;
 120     gmx_bool bLog;
 121     int unit;
 122     real zProf0;
 123     gmx_bool bProf0Set,bs_verbose;
 124     gmx_bool bHistEq, bTab;
 125     double *tabX,*tabY,tabMin,tabMax,tabDz;
 126     int tabNbins;
 127 } t_UmbrellaOptions;
 128
 129
 130 /* Return j such that xx[j] <= x < xx[j+1] */
 131 void searchOrderedTable(double xx[], int n, double x, int *j)
 132 {
 133     int ju,jm,jl;
 134     int ascending;
 135
 136
 137     jl=-1;
 138     ju=n;
 139     ascending=(xx[n-1] > xx[0]);
 140     while (ju-jl > 1)
 141     {
 142         jm=(ju+jl) >> 1;
 143         if ((x >= xx[jm]) == ascending)
 144             jl=jm;
 145         else
 146             ju=jm;
 147     }
 148     if (x==xx[0]) *j=0;
 149     else if (x==xx[n-1]) *j=n-2;
 150     else *j=jl;
 151 }
 152
 153
 154 /* Read and setup tabulated umbrella potential */
 155 void setup_tab(const char *fn,t_UmbrellaOptions *opt)
 156 {
 157     int i,ny,nl;
 158     double **y;
 159
 160
 161     printf("Setting up tabulated potential from file %s\n",fn);
 162     nl=read_xvg(fn,&y,&ny);
 163     opt->tabNbins=nl;
 164     if (ny!=2)
 165         gmx_fatal(FARGS,"Found %d columns in %s. Expected 2.\n",ny,fn);
 166     opt->tabMin=y[0][0];
 167     opt->tabMax=y[0][nl-1];
 168     opt->tabDz=(opt->tabMax-opt->tabMin)/(nl-1);
 169     if (opt->tabDz<=0)
 170         gmx_fatal(FARGS,"The tabulated potential in %s must be provided in \n"
 171                 "ascending z-direction",fn);
 172     for (i=0;i<nl-1;i++)
 173         if  (fabs(y[0][i+1]-y[0][i]-opt->tabDz) > opt->tabDz/1e6)
 174             gmx_fatal(FARGS,"z-values in %s are not equally spaced.\n",ny,fn);
 175     snew(opt->tabY,nl);
 176     snew(opt->tabX,nl);
 177     for (i=0;i<nl;i++){
 178         opt->tabX[i]=y[0][i];
 179         opt->tabY[i]=y[1][i];
 180     }
 181     printf("Found equally spaced tabulated potential from %g to %g, spacing %g\n",
 182             opt->tabMin,opt->tabMax,opt->tabDz);
 183 }
 184
 185
 186 void read_pdo_header(FILE * file,t_UmbrellaHeader * header,
 187                      t_UmbrellaOptions *opt)
 188 {
 189     char Buffer0[256];
 190     char Buffer1[256];
 191     char Buffer2[256];
 192     char Buffer3[256];
 193     char Buffer4[256];
 194     int i,j;
 195
 196
 197     /*  line 1 */
 198     if(3 != fscanf(file,"%s%s%s",Buffer0,Buffer1,Buffer2))
 199     {
 200         gmx_fatal(FARGS,"Error reading header from pdo file");
 201     }
 202     if(strcmp(Buffer1,"UMBRELLA"))
 203         gmx_fatal(FARGS,"This does not appear to be a valid pdo file. Found %s, expected %s",
 204                 Buffer1, "UMBRELLA");
 205     if(strcmp(Buffer2,"3.0"))
 206         gmx_fatal(FARGS,"This does not appear to be a version 3.0 pdo file");
 207
 208     /*  line 2 */
 209     if(6 != fscanf(file,"%s%s%s%d%d%d",Buffer0,Buffer1,Buffer2,
 210                    &(header->Dims[0]),&(header->Dims[1]),&(header->Dims[2])))
 211     {
 212         gmx_fatal(FARGS,"Error reading dimensions in header from pdo file");
 213     }
 214
 215     /* printf("%d %d %d\n", header->Dims[0],header->Dims[1],header->Dims[2]); */
 216
 217     header->nDim = header->Dims[0] + header->Dims[1] + header->Dims[2];
 218     if(header->nDim!=1)
 219         gmx_fatal(FARGS,"Currently only supports one dimension");
 220
 221     /* line3 */
 222     if(3 != fscanf(file,"%s%s%d",Buffer0,Buffer1,&(header->nSkip)))
 223     {
 224         gmx_fatal(FARGS,"Error reading header from pdo file");
 225     }
 226
 227     /* line 4 */
 228     if(4 != fscanf(file,"%s%s%s%s",Buffer0,Buffer1,Buffer2,header->Reference))
 229     {
 230         gmx_fatal(FARGS,"Error reading header from pdo file");
 231     }
 232
 233     /* line 5 */
 234     if(6 != fscanf(file,"%s%s%s%s%s%d",Buffer0,Buffer1,Buffer2,Buffer3,Buffer4,&(header->nPull)))
 235     {
 236         gmx_fatal(FARGS,"Error reading header from pdo file");
 237     }
 238
 239     if (opt->verbose)
 240         printf("Found nPull=%d , nSkip=%d, ref=%s\n",header->nPull,header->nSkip,
 241                 header->Reference);
 242
 243     for(i=0;i<header->nPull;++i)
 244     {
 245       if(4 != fscanf(file,"%s%s%s%s",Buffer0,Buffer1,Buffer2,header->PullName[i]))
 246       {
 247           gmx_fatal(FARGS,"Error reading header from pdo file");
 248       }
 249         if (opt->verbose)
 250             printf("pullgroup %d, pullname = %s\n",i,header->PullName[i]);
 251         for(j=0;j<header->nDim;++j)
 252         {
 253           if(6 != fscanf(file,"%s%s%lf%s%s%lf",Buffer0,Buffer1,&(header->UmbPos[i][j]),
 254                          Buffer2,Buffer3,&(header->UmbCons[i][j])))
 255           {
 256               gmx_fatal(FARGS,"Error reading header from pdo file");
 257           }
 258
 259             if (opt->bFlipProf)
 260             {
 261                 /* We want to combine both halves of a profile into one */
 262                 if(header->UmbPos[i][j]<0)
 263                 {
 264                     header->UmbPos[i][j]= -header->UmbPos[i][j];
 265                     header->Flipped[i]=TRUE;
 266                 }
 267             }
 268             else header->Flipped[i]=FALSE;
 269             /*printf("%f\t%f\n",header->UmbPos[i][j],header->UmbCons[i][j]);*/
 270         }
 271     }
 272
 273     if(1 != fscanf(file,"%s",Buffer3))
 274     {
 275         gmx_fatal(FARGS,"Error reading header from pdo file");
 276     }
 277
 278     if (strcmp(Buffer3,"#####") != 0)
 279         gmx_fatal(FARGS,"Expected '#####', found %s. Hick.\n",Buffer3);
 280 }
 281
 282
 283 static char *fgets3(FILE *fp,char ptr[],int *len)
 284 {
 285     char *p;
 286     int  slen;
 287
 288
 289     if (fgets(ptr,*len-1,fp) == NULL)
 290         return NULL;
 291     p = ptr;
 292     while ((strchr(ptr,'\n') == NULL) && (!feof(fp)))
 293     {
 294         /* This line is longer than len characters, let's increase len! */
 295         *len += STRLEN;
 296         p    += STRLEN;
 297         srenew(ptr,*len);
 298         if (fgets(p-1,STRLEN,fp) == NULL)
 299             break;
 300     }
 301     slen = strlen(ptr);
 302     if (ptr[slen-1] == '\n')
 303         ptr[slen-1] = '\0';
 304
 305     return ptr;
 306 }
 307
 308
 309 void read_pdo_data(FILE * file, t_UmbrellaHeader * header,
 310         int fileno, t_UmbrellaWindow * win,
 311         t_UmbrellaOptions *opt,
 312         gmx_bool bGetMinMax,real *mintmp,real *maxtmp)
 313 {
 314     int i,inttemp,bins,count;
 315     real min,max,minfound,maxfound;
 316     double temp,time,time0=0,dt;
 317     char *ptr;
 318     t_UmbrellaWindow * window=0;
 319     gmx_bool timeok,dt_ok=1;
 320     char  *tmpbuf,fmt[256],fmtign[256];
 321     int    len=STRLEN,dstep=1;
 322
 323         minfound=1e20;
 324         maxfound=-1e20;
 325
 326     if (!bGetMinMax)
 327     {
 328         bins=opt->bins;
 329         min=opt->min;
 330         max=opt->max;
 331
 332         window=win+fileno;
 333         /* Need to alocate memory and set up structure */
 334         window->nPull=header->nPull;
 335         window->nBin=bins;
 336
 337         snew(window->Histo,window->nPull);
 338         snew(window->z,window->nPull);
 339         snew(window->k,window->nPull);
 340         snew(window->pos,window->nPull);
 341         snew(window->Flipped,window->nPull);
 342         snew(window->N, window->nPull);
 343         snew(window->Ntot, window->nPull);
 344
 345         for(i=0;i<window->nPull;++i)
 346         {
 347             window->z[i]=1;
 348             snew(window->Histo[i],bins);
 349             window->k[i]=header->UmbCons[i][0];
 350             window->pos[i]=header->UmbPos[i][0];
 351             window->Flipped[i]=header->Flipped[i];
 352             window->N[i]=0;
 353             window->Ntot[i]=0;
 354         }
 355         /* Done with setup */
 356     }
 357     else
 358     {
 359         minfound=1e20;
 360         maxfound=-1e20;
 361         min=max=bins=0; /* Get rid of warnings */
 362     }
 363
 364     count=0;
 365     snew(tmpbuf,len);
 366     while ( (ptr=fgets3(file,tmpbuf,&len)) != NULL)
 367     {
 368         trim(ptr);
 369
 370         if (ptr[0] == '#' || strlen(ptr)<2)
 371             continue;
 372
 373         /* Initiate format string */
 374         fmtign[0] = '\0';
 375         strcat(fmtign,"%*s");
 376
 377         sscanf(ptr,"%lf",&time); /* printf("Time %f\n",time); */
 378         /* Round time to fs */
 379         time=1.0/1000*( (int) (time*1000+0.5) );
 380
 381         /* get time step of pdo file */
 382         if (count==0)
 383             time0=time;
 384         else if (count==1)
 385         {
 386             dt=time-time0;
 387             if (opt->dt>0.0)
 388             {
 389               dstep=(int)(opt->dt/dt+0.5);
 390                 if (dstep==0)
 391                     dstep=1;
 392             }
 393             if (!bGetMinMax)
 394                 window->dt=dt*dstep;
 395         }
 396         count++;
 397
 398         dt_ok=((count-1)%dstep == 0);
 399         timeok=(dt_ok && time >= opt->tmin && time <= opt->tmax);
 400         /* if (opt->verbose)
 401       printf(" time = %f, (tmin,tmax)=(%e,%e), dt_ok=%d timeok=%d\n",
 402       time,opt->tmin, opt->tmax, dt_ok,timeok); */
 403
 404         if (timeok)
 405         {
 406             for(i=0;i<header->nPull;++i)
 407             {
 408                 strcpy(fmt,fmtign);
 409                 strcat(fmt,"%lf");      /* Creating a format stings such as "%*s...%*s%lf" */
 410                 strcat(fmtign,"%*s");   /* ignoring one more entry in the next loop */
 411                 if(sscanf(ptr,fmt,&temp))
 412                 {
 413                     if(opt->bFlipProf)
 414                     {
 415                         if(header->Flipped[i]) temp=-temp;
 416                     }
 417
 418                     temp+=header->UmbPos[i][0];
 419                     if (bGetMinMax){
 420                         if (temp<minfound)
 421                             minfound=temp;
 422                         if (temp>maxfound)
 423                             maxfound=temp;
 424                     }
 425                     else
 426                     {
 427                         temp-=min;
 428                         temp/=(max-min);
 429                         temp*=bins;
 430                         temp=floor(temp);
 431
 432                         inttemp = (int)temp;
 433                         if (opt->cycl==enCycl_yes)
 434                         {
 435                             if (inttemp < 0)
 436                                 inttemp+=bins;
 437                             else if (inttemp >= bins)
 438                                 inttemp-=bins;
 439                         }
 440
 441                         if(inttemp >= 0 && inttemp < bins)
 442                         {
 443                             window->Histo[i][inttemp]+=1;
 444                             window->N[i]++;
 445                         }
 446                         window->Ntot[i]++;
 447                     }
 448                 }
 449             }
 450         }
 451         if (time>opt->tmax)
 452         {
 453             if (opt->verbose)
 454                 printf("time %f larger than tmax %f, stop reading pdo file\n",time,opt->tmax);
 455             break;
 456         }
 457     }
 458
 459     if (bGetMinMax)
 460     {
 461         *mintmp=minfound;
 462         *maxtmp=maxfound;
 463     }
 464 }
 465
 466
 467 void enforceEqualWeights(t_UmbrellaWindow * window,int nWindows)
 468 {
 469     int i,k,j,NEnforced;
 470     double ratio;
 471
 472
 473     NEnforced=window[0].Ntot[0];
 474     printf("\nFound -hist-eq. Enforcing equal weights for all histograms, \ni.e. doing a "
 475             "non-weighted histogram analysis method. Ndata = %d\n",NEnforced);
 476     /* enforce all histograms to have the same weight as the very first histogram */
 477
 478     for(j=0;j<nWindows;++j)
 479         for(k=0;k<window[j].nPull;++k)
 480         {
 481             ratio=1.0*NEnforced/window[j].Ntot[k];
 482             for(i=0;i<window[0].nBin;++i)
 483                 window[j].Histo[k][i]*=ratio;
 484             window[j].N[k]=(int)(ratio*window[j].N[k]+0.5);
 485         }
 486 }
 487
 488
 489 /* Simple linear interpolation between two given tabulated points */
 490 double tabulated_pot(double dist, t_UmbrellaOptions *opt)
 491 {
 492     int jl,ju;
 493     double pl,pu,dz,dp;
 494
 495
 496     jl=floor((dist-opt->tabMin)/opt->tabDz);
 497     ju=jl+1;
 498     if (jl<0 || ju>=opt->tabNbins)
 499         gmx_fatal(FARGS,"Distance %f out of bounds of tabulated potential (jl=%d, ju=%d).\n"
 500                 "Provide an extended table.",dist,jl,ju);
 501     pl=opt->tabY[jl];
 502     pu=opt->tabY[ju];
 503     dz=dist-opt->tabX[jl];
 504     dp=(pu-pl)*dz/opt->tabDz;
 505     return pl+dp;
 506 }
 507
 508
 509 /* Don't worry, that routine does not mean we compute the PMF in limited precision.
 510    After rapid convergence (using only substiantal contributions), we always switch to
 511    full precision. */
 512 #define WHAM_CONTRIB_LIM 1e-10
 513 void setup_acc_wham(t_UmbrellaWindow * window,int nWindows, t_UmbrellaOptions *opt)
 514 {
 515     int i,j,k;
 516     double U,min=opt->min,dz=opt->dz,temp,ztot_half,distance,ztot,contrib;
 517     gmx_bool bAnyContrib;
 518
 519
 520     ztot=opt->max-opt->min;
 521     ztot_half=ztot/2;
 522
 523     for(i=0;i<nWindows;++i)
 524     {
 525         snew(window[i].bContrib,window[i].nPull);
 526         for(j=0;j<window[i].nPull;++j)
 527         {
 528             snew(window[i].bContrib[j],opt->bins);
 529             bAnyContrib=FALSE;
 530             for(k=0;k<opt->bins;++k)
 531             {
 532                 temp=(1.0*k+0.5)*dz+min;
 533                 distance = temp - window[i].pos[j];   /* distance to umbrella center */
 534                 if (opt->cycl==enCycl_yes)
 535                 {                                     /* in cyclic wham:             */
 536                     if (distance > ztot_half)           /*    |distance| < ztot_half   */
 537                         distance-=ztot;
 538                     else if (distance < -ztot_half)
 539                         distance+=ztot;
 540                 }
 541                 if (!opt->bTab)
 542                     U=0.5*window[i].k[j]*sqr(distance);       /* harmonic potential assumed. */
 543                 else
 544                     U=tabulated_pot(distance,opt);            /* Use tabulated potential     */
 545
 546                 contrib=exp(- U/(8.314e-3*opt->Temperature));
 547                 window[i].bContrib[j][k] = (contrib > WHAM_CONTRIB_LIM);
 548                 bAnyContrib = (bAnyContrib | window[i].bContrib[j][k]);
 549             }
 550             /* If this histo is far outside min and max all bContrib may be FALSE,
 551              * causing a floating point exception later on. To avoid that, switch
 552              * them all to true.*/
 553             if (!bAnyContrib)
 554                 for(k=0;k<opt->bins;++k)
 555                     window[i].bContrib[j][k]=TRUE;
 556         }
 557     }
 558     printf("Initialized rapid wham stuff.\n");
 559 }
 560
 561
 562 void calc_profile(double *profile,t_UmbrellaWindow * window, int nWindows, t_UmbrellaOptions *opt,
 563         gmx_bool bExact)
 564 {
 565     int i,k,j;
 566     double num,ztot_half,ztot,distance,min=opt->min,dz=opt->dz;
 567     double denom,U=0,temp=0;
 568
 569
 570     ztot=opt->max-opt->min;
 571     ztot_half=ztot/2;
 572
 573
 574     for(i=0;i<opt->bins;++i)
 575     {
 576         num=denom=0;
 577         for(j=0;j<nWindows;++j)
 578         {
 579             for(k=0;k<window[j].nPull;++k)
 580             {
 581                 temp=(1.0*i+0.5)*dz+min;
 582                 num+=window[j].Histo[k][i];
 583
 584                 if (! (bExact || window[j].bContrib[k][i]))
 585                     continue;
 586                 distance = temp - window[j].pos[k];   /* distance to umbrella center */
 587                 if (opt->cycl==enCycl_yes){           /* in cyclic wham:             */
 588                     if (distance > ztot_half)           /*    |distance| < ztot_half   */
 589                         distance-=ztot;
 590                     else if (distance < -ztot_half)
 591                         distance+=ztot;
 592                 }
 593
 594                 if (!opt->bTab)
 595                     U=0.5*window[j].k[k]*sqr(distance);       /* harmonic potential assumed. */
 596                 else
 597                     U=tabulated_pot(distance,opt);            /* Use tabulated potential     */
 598                 denom+=window[j].N[k]*exp(- U/(8.314e-3*opt->Temperature) + window[j].z[k]);
 599             }
 600         }
 601         profile[i]=num/denom;
 602     }
 603 }
 604
 605
 606 double calc_z(double * profile,t_UmbrellaWindow * window, int nWindows, t_UmbrellaOptions *opt,
 607         gmx_bool bExact)
 608 {
 609     int i,j,k;
 610     double U=0,min=opt->min,dz=opt->dz,temp,ztot_half,distance,ztot;
 611     double MAX=-1e20, total=0;
 612
 613
 614     ztot=opt->max-opt->min;
 615     ztot_half=ztot/2;
 616
 617     for(i=0;i<nWindows;++i)
 618     {
 619         for(j=0;j<window[i].nPull;++j)
 620         {
 621             total=0;
 622             for(k=0;k<window[i].nBin;++k)
 623             {
 624                 if (! (bExact || window[i].bContrib[j][k]))
 625                     continue;
 626                 temp=(1.0*k+0.5)*dz+min;
 627                 distance = temp - window[i].pos[j];   /* distance to umbrella center */
 628                 if (opt->cycl==enCycl_yes)
 629                 {                                     /* in cyclic wham:             */
 630                     if (distance > ztot_half)           /*    |distance| < ztot_half   */
 631                         distance-=ztot;
 632                     else if (distance < -ztot_half)
 633                         distance+=ztot;
 634                 }
 635
 636                 if (!opt->bTab)
 637                     U=0.5*window[i].k[j]*sqr(distance);       /* harmonic potential assumed. */
 638                 else
 639                     U=tabulated_pot(distance,opt);            /* Use tabulated potential     */
 640
 641                 total+=profile[k]*exp(-U/(8.314e-3*opt->Temperature));
 642             }
 643             if (total > 0.0)
 644                 total = -log(total);
 645             else
 646                 total = 1000.0;
 647             temp = fabs(total - window[i].z[j]);
 648             if(temp > MAX) MAX=temp;
 649             window[i].z[j] = total;
 650         }
 651     }
 652     return MAX;
 653 }
 654
 655
 656 void cyclicProfByWeightedCorr(double *profile,t_UmbrellaWindow *window,
 657                               int nWindows, t_UmbrellaOptions * opt,
 658                               gmx_bool bAppendCorr2File, const char *fn,
 659                               const output_env_t oenv)
 660 {
 661     int i,j,k,bins=opt->bins;
 662     static int first=1;
 663     double *weight,sum=0.,diff,*histsum,*corr,sumCorr=0.,dCorr;
 664     FILE *fp;
 665     char buf[256];
 666
 667     if (first)
 668     {
 669         printf("\nEnforcing a periodic profile by a sampling wheighted correction.");
 670         please_cite(stdout,"Hub2008");
 671     }
 672
 673     snew(weight,bins-1);
 674     snew(histsum,bins);
 675     snew(corr,bins-1);
 676
 677     /* generate weights proportional to 1/(n(i)*n(i+1))^alpha
 678      where n(i) is the total nur of data points in bin i from all histograms */
 679     for(i=0;i<nWindows;++i)
 680         for(j=0;j<window[i].nPull;++j)
 681             for(k=0;k<bins;++k)
 682                 histsum[k]+=window[i].Histo[j][k];
 683
 684     for(k=0,sum=0.;k<bins-1;++k)
 685     {
 686         weight[k]=1./pow(histsum[k]*histsum[k+1],opt->alpha);
 687         sum+=weight[k];
 688     }
 689     for(k=0;k<bins-1;++k)
 690         weight[k]/=sum;
 691
 692     /* difference between last and first bin */
 693     diff=profile[bins-1]-profile[0];
 694     printf("Distributing %f between adjacent bins to enforce a cyclic profile\n",diff);
 695
 696     for (i=0;i<bins-1;i++)
 697     {
 698         dCorr=weight[i]*diff;
 699         sumCorr+=dCorr;
 700         corr[i]=sumCorr;
 701     }
 702
 703     for (i=0;i<bins-1;i++)
 704         profile[i+1]-=corr[i];
 705
 706     if (bAppendCorr2File)
 707     {
 708         fp=xvgropen(fn,"Corrections to enforce periodicity","z",
 709                     "\\f{12}D\\f{}G(z)",oenv);
 710         sprintf(buf,"corrections propotional to 1/(n\\si\\Nn\\si+1\\N)\\S%.2f",
 711                 opt->alpha);
 712         xvgr_subtitle(fp,buf,oenv);
 713         for (i=0;i<bins-1;i++)
 714             fprintf(fp,"%g %g\n",opt->min+opt->dz*(i+1),-corr[i]);
 715         ffclose(fp);
 716     }
 717
 718     sfree(histsum);
 719     sfree(corr);
 720     sfree(weight);
 721     first=0;
 722 }
 723
 724
 725 void prof_normalization_and_unit(double * profile, t_UmbrellaOptions *opt)
 726 {
 727     int i,bins,imin;
 728     double unit_factor=1., R_MolarGasConst, diff;
 729
 730
 731     R_MolarGasConst=8.314472e-3; /* in kJ/(mol*K) */
 732     bins=opt->bins;
 733
 734     /* No log? Nothing to do! */
 735     if (!opt->bLog)
 736         return;
 737
 738     /* Get profile in units of kT, kJ/mol, or kCal/mol */
 739     if (opt->unit == en_kT)
 740         unit_factor=1.0;
 741     else if (opt->unit == en_kJ)
 742         unit_factor=R_MolarGasConst*opt->Temperature;
 743     else if (opt->unit == en_kCal)
 744         unit_factor=R_MolarGasConst*opt->Temperature/4.1868;
 745     else
 746         gmx_fatal(FARGS,"Sorry I don't know this energy unit.");
 747
 748     for (i=0;i<bins;i++)
 749         if (profile[i]>0.0)
 750             profile[i]=-log(profile[i])*unit_factor;
 751
 752     /* shift to zero at z=opt->zProf0 */
 753     if (!opt->bProf0Set)
 754         diff=profile[0];
 755     else{
 756         /* Get bin with shortest distance to opt->zProf0 */
 757       imin=(int)((opt->zProf0-opt->min)/opt->dz);
 758         if (imin<0)
 759             imin=0;
 760         else if (imin>=bins)
 761             imin=bins-1;
 762         diff=profile[imin];
 763     }
 764
 765     /* Shift to zero */
 766     for (i=0;i<bins;i++)
 767         profile[i]-=diff;
 768 }
 769
 770
 771 void getRandomIntArray(int nPull,int blockLength,int* randomArray)
 772 {
 773     int ipull,blockBase,nr,ipullRandom;
 774
 775
 776     if (blockLength==0)
 777         blockLength=nPull;
 778
 779     for (ipull=0; ipull<nPull; ipull++)
 780     {
 781         blockBase=(ipull/blockLength)*blockLength;
 782         do{      /* make sure nothing bad happens in the last block */
 783             nr=(int)((1.0*rand()/RAND_MAX)*blockLength);
 784             ipullRandom = blockBase + nr;
 785         } while (ipullRandom >= nPull);
 786         if (ipullRandom<0 || ipullRandom>=nPull)
 787             gmx_fatal(FARGS,"Ups, random iWin = %d, nPull = %d, nr = %d, "
 788                     "blockLength = %d, blockBase = %d\n",
 789                     ipullRandom,nPull,nr,blockLength,blockBase);
 790         randomArray[ipull]=ipullRandom;
 791     }
 792     /*for (ipull=0; ipull<nPull; ipull++)
 793     printf("%d ",randomArray[ipull]); printf("\n"); */
 794 }
 795
 796
 797 void copy_pullgrp_to_synthwindow(t_UmbrellaWindow *synthWindow,
 798         t_UmbrellaWindow *thisWindow,int pullid)
 799 {
 800     synthWindow->N       [0]=thisWindow->N        [pullid];
 801     synthWindow->Histo   [0]=thisWindow->Histo    [pullid];
 802     synthWindow->pos     [0]=thisWindow->pos      [pullid];
 803     synthWindow->z       [0]=thisWindow->z        [pullid];
 804     synthWindow->k       [0]=thisWindow->k        [pullid];
 805     synthWindow->bContrib[0]=thisWindow->bContrib [pullid];
 806 }
 807
 808
 809 /* Calculate cummulative of all histograms. They allow to create random numbers
 810    which are distributed according to the histograms. Required to generate
 811    the "synthetic" histograms for the Bootstrap method */
 812 void calc_cummulants(t_UmbrellaWindow *window,int nWindows,
 813                      t_UmbrellaOptions *opt,const char *fnhist,
 814                      const output_env_t oenv)
 815 {
 816     int i,j,k,nbin;
 817     double last;
 818     char *fn=0,*buf=0;
 819     FILE *fp=0;
 820
 821
 822     if (opt->bs_verbose)
 823     {
 824         snew(fn,strlen(fnhist)+10);
 825         snew(buf,strlen(fnhist)+10);
 826         sprintf(fn,"%s_cumul.xvg",strncpy(buf,fnhist,strlen(fnhist)-4));
 827         fp=xvgropen(fn,"Cummulants of umbrella histograms","z","cummulant",
 828                     oenv);
 829     }
 830
 831     nbin=opt->bins;
 832     for (i=0; i<nWindows; i++)
 833     {
 834         snew(window[i].cum,window[i].nPull);
 835         for (j=0; j<window[i].nPull; j++)
 836         {
 837             snew(window[i].cum[j],nbin+1);
 838             window[i].cum[j][0]=0.;
 839             for (k=1; k<=nbin; k++)
 840                 window[i].cum[j][k] = window[i].cum[j][k-1]+window[i].Histo[j][k-1];
 841
 842             /* normalize cummulant. Ensure cum[nbin]==1 */
 843             last = window[i].cum[j][nbin];
 844             for (k=0; k<=nbin; k++)
 845                 window[i].cum[j][k] /= last;
 846         }
 847     }
 848
 849     printf("Cumulants of all histograms created.\n");
 850     if (opt->bs_verbose)
 851     {
 852         for (k=0; k<=nbin; k++)
 853         {
 854             fprintf(fp,"%g\t",opt->min+k*opt->dz);
 855             for (i=0; i<nWindows; i++)
 856                 for (j=0; j<window[i].nPull; j++)
 857                     fprintf(fp,"%g\t",window[i].cum[j][k]);
 858             fprintf(fp,"\n");
 859         }
 860         printf("Wrote cumulants to %s\n",fn);
 861         ffclose(fp);
 862         sfree(fn);
 863         sfree(buf);
 864     }
 865 }
 866
 867
 868 /* Return j such that xx[j] <= x < xx[j+1] */
 869 void searchCummulant(double xx[], int n, double x, int *j)
 870 {
 871     int ju,jm,jl;
 872
 873
 874     jl=-1;
 875     ju=n;
 876     while (ju-jl > 1)
 877     {
 878         jm=(ju+jl) >> 1;
 879         if (x >= xx[jm])
 880             jl=jm;
 881         else
 882             ju=jm;
 883     }
 884     if (x==xx[0])
 885         *j=0;
 886     else if (x==xx[n-1])
 887         *j=n-2;
 888     else
 889         *j=jl;
 890 }
 891
 892
 893 void create_synthetic_histo(t_UmbrellaWindow *synthWindow,
 894                             t_UmbrellaWindow *thisWindow,
 895                             int pullid,t_UmbrellaOptions *opt)
 896 {
 897     int nsynth,N,i,nbins,r_index;
 898     double r;
 899     static gmx_bool bWarnout=0;
 900
 901
 902     N=thisWindow->N[pullid];
 903     nbins=thisWindow->nBin;
 904
 905     /* nsynth = nr of data points in synthetic histo */
 906     if (opt->dtBootStrap==0.0)
 907         nsynth=N;
 908     else
 909     {
 910       nsynth=(int)(thisWindow->N[pullid]*thisWindow->dt/opt->dtBootStrap+0.5);
 911         if (nsynth>N)
 912             nsynth=N;
 913     }
 914
 915     if (!bWarnout && nsynth<10)
 916     {
 917         printf("\n++++ WARNING ++++\n\tOnly %d data points per synthetic histogram!\n"
 918                 "\tYou may want to consider option -bs-dt.\n\n",nsynth);
 919         bWarnout=1;
 920     }
 921
 922     synthWindow->N       [0]=nsynth;
 923     synthWindow->pos     [0]=thisWindow->pos[pullid];
 924     synthWindow->z       [0]=thisWindow->z[pullid];
 925     synthWindow->k       [0]=thisWindow->k[pullid];
 926     synthWindow->bContrib[0]=thisWindow->bContrib[pullid];
 927
 928     for (i=0;i<nbins;i++)
 929         synthWindow->Histo[0][i]=0.;
 930
 931     for (i=0;i<nsynth;i++)
 932     {
 933         r=1.0*rand()/RAND_MAX;
 934         searchCummulant(thisWindow->cum[pullid],nbins+1 ,r,&r_index);
 935         synthWindow->Histo[0][r_index]+=1.;
 936     }
 937 }
 938
 939
 940 void print_histograms(const char *fnhist, t_UmbrellaWindow * window,
 941                       int nWindows, int bs_index,t_UmbrellaOptions *opt,
 942                       const output_env_t oenv)
 943 {
 944     char *fn;
 945     char *buf=0,title[256];
 946     FILE *fp;
 947     int bins,l,i,j;
 948
 949
 950     if (bs_index<0)
 951     {
 952         fn=strdup(fnhist);
 953         strcpy(title,"Umbrella histograms");
 954     }
 955     else
 956     {
 957         snew(fn,strlen(fnhist)+6);
 958         snew(buf,strlen(fnhist)+1);
 959         sprintf(fn,"%s_bs%d.xvg",strncpy(buf,fnhist,strlen(fnhist)-4),bs_index);
 960         sprintf(title,"Umbrella histograms. Bootstrap #%d",bs_index);
 961     }
 962
 963     fp=xvgropen(fn,title,"z","count",oenv);
 964     bins=opt->bins;
 965
 966     /* Write histograms */
 967     for(l=0;l<bins;++l)
 968     {
 969         fprintf(fp,"%e\t",(double)(l+0.5)*opt->dz+opt->min);
 970         for(i=0;i<nWindows;++i)
 971         {
 972             for(j=0;j<window[i].nPull;++j)
 973             {
 974                 fprintf(fp,"%e\t",window[i].Histo[j][l]);
 975             }
 976         }
 977         fprintf(fp,"\n");
 978     }
 979
 980     ffclose(fp);
 981     printf("Wrote %s\n",fn);
 982     if (buf)
 983     {
 984         sfree(buf);
 985         sfree(fn);
 986     }
 987 }
 988
 989
 990 void do_bootstrapping(const char *fnres, const char* fnprof,
 991                       const char *fnhist, char* ylabel, double *profile,
 992                       t_UmbrellaWindow * window, int nWindows,
 993                       t_UmbrellaOptions *opt, const output_env_t oenv)
 994 {
 995     t_UmbrellaWindow * synthWindow;
 996     double *bsProfile,*bsProfiles_av, *bsProfiles_av2,maxchange=1e20,tmp,stddev;
 997     int i,j,*randomArray=0,winid,pullid,ib;
 998     int iAllPull,nAllPull,*allPull_winId,*allPull_pullId;
 999     FILE *fp;
1000     gmx_bool bExact=FALSE;
1001
1002
1003     /* init random */
1004     if (opt->bsSeed==-1)
1005         srand(time(NULL));
1006     else
1007         srand(opt->bsSeed);
1008
1009     snew(bsProfile,     opt->bins);
1010     snew(bsProfiles_av, opt->bins);
1011     snew(bsProfiles_av2,opt->bins);
1012
1013     /* Create array of all pull groups. Note that different windows
1014      may have different nr of pull groups
1015      First: Get total nr of pull groups */
1016     nAllPull=0;
1017     for (i=0;i<nWindows;i++)
1018         nAllPull+=window[i].nPull;
1019     snew(allPull_winId,nAllPull);
1020     snew(allPull_pullId,nAllPull);
1021     iAllPull=0;
1022     /* Setup one array of all pull groups */
1023     for (i=0;i<nWindows;i++)
1024         for (j=0;j<window[i].nPull;j++)
1025         {
1026             allPull_winId[iAllPull]=i;
1027             allPull_pullId[iAllPull]=j;
1028             iAllPull++;
1029         }
1030
1031     /* setup stuff for synthetic windows */
1032     snew(synthWindow,nAllPull);
1033     for (i=0;i<nAllPull;i++)
1034     {
1035         synthWindow[i].nPull=1;
1036         synthWindow[i].nBin=opt->bins;
1037         snew(synthWindow[i].Histo,1);
1038         if (!opt->bHistBootStrap)
1039             snew(synthWindow[i].Histo[0],opt->bins);
1040         snew(synthWindow[i].N,1);
1041         snew(synthWindow[i].pos,1);
1042         snew(synthWindow[i].z,1);
1043         snew(synthWindow[i].k,1);
1044         snew(synthWindow[i].bContrib,1);
1045     }
1046
1047     if (opt->bHistBootStrap)
1048     {
1049         snew(randomArray,nAllPull);
1050         printf("\n\nWhen computing statistical errors by bootstrapping entire histograms:\n");
1051         please_cite(stdout,"Hub2006");
1052     }
1053     else
1054     {
1055         calc_cummulants(window,nWindows,opt,fnhist,oenv);
1056     }
1057
1058     /* do bootstrapping */
1059     fp=xvgropen(fnprof,"Boot strap profiles","z",ylabel,oenv);
1060     for (ib=0;ib<opt->nBootStrap;ib++)
1061     {
1062         printf("  *******************************************\n"
1063                 "  ******** Start bootstrap nr %d ************\n"
1064                 "  *******************************************\n",ib+1);
1065
1066         if (opt->bHistBootStrap)
1067         {
1068             /* only mix given histos */
1069             getRandomIntArray(nAllPull,opt->histBootStrapBlockLength,randomArray);
1070             for (i=0;i<nAllPull;i++)
1071             {
1072                 winid =allPull_winId [randomArray[i]];
1073                 pullid=allPull_pullId[randomArray[i]];
1074                 copy_pullgrp_to_synthwindow(synthWindow+i,window+winid,pullid);
1075             }
1076         }
1077         else
1078         {
1079             /* create new histos from given histos */
1080             for (i=0;i<nAllPull;i++)
1081             {
1082                 winid=allPull_winId[i];
1083                 pullid=allPull_pullId[i];
1084                 create_synthetic_histo(synthWindow+i,window+winid,pullid,opt);
1085             }
1086         }
1087
1088         /* print histos in case of verbose output */
1089         if (opt->bs_verbose)
1090             print_histograms(fnhist,synthWindow,nAllPull,ib,opt,oenv);
1091
1092         /* do wham */
1093         i=0;
1094         bExact=FALSE;
1095         maxchange=1e20;
1096         memcpy(bsProfile,profile,opt->bins*sizeof(double)); /* use profile as guess */
1097         do {
1098             if (maxchange<opt->Tolerance)
1099                 bExact=TRUE;
1100             if (((i%opt->stepchange) == 0 || i==1) && !i==0)
1101                 printf("\t%4d) Maximum change %e\n",i,maxchange);
1102             calc_profile(bsProfile,synthWindow,nAllPull,opt,bExact);
1103             i++;
1104         } while( (maxchange=calc_z(bsProfile, synthWindow, nAllPull, opt,bExact)) > opt->Tolerance || !bExact);
1105         printf("\tConverged in %d iterations. Final maximum change %g\n",i,maxchange);
1106
1107         if (opt->bLog)
1108             prof_normalization_and_unit(bsProfile,opt);
1109         /* Force cyclic profile by wheighted correction */
1110         if (opt->cycl==enCycl_weighted)
1111             cyclicProfByWeightedCorr(bsProfile,synthWindow,nAllPull,opt,
1112                                      FALSE, 0,oenv);
1113
1114         /* save stuff to get average and stddev */
1115         for (i=0;i<opt->bins;i++)
1116         {
1117             tmp=bsProfile[i];
1118             bsProfiles_av[i]+=tmp;
1119             bsProfiles_av2[i]+=tmp*tmp;
1120             fprintf(fp,"%e\t%e\n",(i+0.5)*opt->dz+opt->min,tmp);
1121         }
1122         fprintf(fp,"&\n");
1123     }
1124     ffclose(fp);
1125
1126     /* write average and stddev */
1127     fp=ffopen(fnres,"w");
1128     fprintf(fp,"@    title \"%s\"\n","Average and stddev from bootstrapping");
1129     fprintf(fp,"@    xaxis  label \"%s\"\n","z");
1130     fprintf(fp,"@    yaxis  label \"%s\"\n",ylabel);
1131     fprintf(fp,"@TYPE xydy\n");
1132     for (i=0;i<opt->bins;i++)
1133     {
1134         bsProfiles_av [i]/=opt->nBootStrap;
1135         bsProfiles_av2[i]/=opt->nBootStrap;
1136         tmp=bsProfiles_av2[i]-sqr(bsProfiles_av[i]);
1137         stddev=(tmp>=0.) ? sqrt(tmp) : 0.; /* Catch rouding errors */
1138         fprintf(fp,"%e\t%e\t%e\n",(i+0.5)*opt->dz+opt->min,bsProfiles_av [i],stddev);
1139     }
1140     ffclose(fp);
1141     printf("Wrote boot strap result to %s\n",fnres);
1142 }
1143
1144
1145 int whaminFileType(const char *fn)
1146 {
1147     int len;
1148     len=strlen(fn);
1149     if (strcmp(fn+len-3,"tpr")==0)
1150         return whamin_tpr;
1151     else if (strcmp(fn+len-3,"xvg")==0 || strcmp(fn+len-6,"xvg.gz")==0)
1152         return whamin_pullxf;
1153     else if (strcmp(fn+len-3,"pdo")==0 || strcmp(fn+len-6,"pdo.gz")==0)
1154         return whamin_pdo;
1155     else
1156         gmx_fatal(FARGS,"Unknown file type of %s. Should be tpr, xvg, or pdo.\n",fn);
1157     return whamin_unknown;
1158 }
1159
1160
1161 void read_wham_in(const char *fn,char ***filenamesRet, int *nfilesRet,
1162                   t_UmbrellaOptions *opt)
1163 {
1164     char **filename,tmp[STRLEN];
1165     int nread,sizenow,i,block=10;
1166     FILE *fp;
1167 #define MAXFILELEN 512
1168
1169
1170     fp=ffopen(fn,"r");
1171     sizenow=block;
1172     snew(filename,sizenow);
1173     for (i=0;i<sizenow;i++)
1174         snew(filename[i],MAXFILELEN);
1175     nread=0;
1176     while (fscanf(fp,"%s",tmp) != EOF)
1177     {
1178         if (strlen(tmp)>=MAXFILELEN)
1179             gmx_fatal(FARGS,"Filename too long. Only %d characters allowed\n",MAXFILELEN);
1180         strcpy(filename[nread],tmp);
1181         if (opt->verbose)
1182             printf("Found file %s in %s\n",filename[nread],fn);
1183         nread++;
1184         if (nread>=sizenow)
1185         {
1186             sizenow+=block;
1187             srenew(filename,sizenow);
1188             for (i=sizenow-block;i<sizenow;i++)
1189                 snew(filename[i],MAXFILELEN);
1190         }
1191     }
1192     *filenamesRet=filename;
1193     *nfilesRet=nread;
1194 }
1195
1196
1197 FILE *pdo_open_file(const char *fn)
1198 {
1199     char Buffer[1024],gunzip[1024],*Path=0;
1200     FILE *fp;
1201
1202     if (!gmx_fexist(fn))
1203         {
1204         gmx_fatal(FARGS,"File %s does not exist.\n",fn);
1205         }
1206
1207     /* gzipped pdo file? */
1208     if (strcmp(fn+strlen(fn)-3,".gz")==0)
1209     {
1210 #ifdef HAVE_PIPES
1211         if(!(Path=getenv("GMX_PATH_GZIP")))
1212             sprintf(gunzip,"%s","/bin/gunzip");
1213         else
1214             sprintf(gunzip,"%s/gunzip",Path);
1215         if (!gmx_fexist(gunzip))
1216             gmx_fatal(FARGS,"Cannot find executable %s. You may want to define the path to gunzip "
1217                     "with the environment variable GMX_PATH_GZIP.",gunzip);
1218         sprintf(Buffer,"%s -c < %s",gunzip,fn);
1219                 if((fp=popen(Buffer,"r"))==NULL)
1220                 {
1221                         gmx_fatal(FARGS,"Unable to open pipe to `%s'\n",Buffer);
1222                 }
1223 #else
1224                 gmx_fatal(FARGS,"Cannot open a compressed file on platform without pipe support");
1225 #endif
1226     }
1227     else
1228         {
1229                 if((fp=ffopen(fn,"r"))==NULL)
1230                 {
1231                         gmx_fatal(FARGS,"Unable to open file %s\n",fn);
1232                 }
1233         }
1234         return fp;
1235 }
1236
1237 void
1238 pdo_close_file(FILE *fp)
1239 {
1240 #ifdef HAVE_PIPES
1241         pclose(fp);
1242 #else
1243         ffclose(fp);
1244 #endif
1245 }
1246
1247 /* Reading pdo files */
1248 void read_pdo_files(char **fn, int nfiles, t_UmbrellaHeader* header,
1249         t_UmbrellaWindow **window, t_UmbrellaOptions *opt)
1250 {
1251     FILE  * file;
1252     real mintmp,maxtmp;
1253     int i;
1254
1255
1256     if(nfiles<1)
1257         gmx_fatal(FARGS,"No files found. Hick.");
1258
1259     /* if min and max are not given, get min and max from the input files */
1260     if (opt->bAuto)
1261     {
1262         printf("Automatic determination of boundaries from %d pdo files...\n",nfiles);
1263         opt->min=1e20;
1264         opt->max=-1e20;
1265         for(i=0;i<nfiles;++i)
1266         {
1267             file=pdo_open_file(fn[i]);
1268             printf("\rOpening %s ...",fn[i]); fflush(stdout);
1269             if (opt->verbose)
1270                 printf("\n");
1271             read_pdo_header(file,header,opt);
1272             /* here only determine min and max of this window */
1273             read_pdo_data(file,header,i,NULL,opt,TRUE,&mintmp,&maxtmp);
1274             if (maxtmp>opt->max)
1275                 opt->max=maxtmp;
1276             if (mintmp<opt->min)
1277                 opt->min=mintmp;
1278             pdo_close_file(file);
1279         }
1280         printf("\n");
1281         printf("\nDetermined boundaries to %f and %f\n\n",opt->min,opt->max);
1282         if (opt->bBoundsOnly)
1283         {
1284             printf("Found option -boundsonly, now exiting.\n");
1285             exit (0);
1286         }
1287     }
1288     /* store stepsize in profile */
1289     opt->dz=(opt->max-opt->min)/opt->bins;
1290
1291     snew(*window,nfiles);
1292
1293     /* Having min and max, we read in all files */
1294     /* Loop over all files */
1295     for(i=0;i<nfiles;++i)
1296     {
1297         printf("\rOpening %s ...",fn[i]); fflush(stdout);
1298         if (opt->verbose)
1299             printf("\n");
1300         file=pdo_open_file(fn[i]);
1301         /* read in the headers */
1302         read_pdo_header(file,header,opt);
1303         /* load data into window */
1304         read_pdo_data(file,header,i,*window,opt,FALSE,NULL,NULL);
1305         pdo_close_file(file);
1306     }
1307     printf("\n");
1308 }
1309
1310
1311 #define int2YN(a) (((a)==0)?("N"):("Y"))
1312
1313 void read_tpr_header(const char *fn,t_UmbrellaHeader* header,
1314                      t_UmbrellaOptions *opt)
1315 {
1316     t_inputrec  ir;
1317     int         i,ngrp,d;
1318     t_state     state;
1319     static int first=1;
1320
1321
1322     /* printf("Reading %s \n",fn); */
1323     read_tpx_state(fn,&ir,&state,NULL,NULL);
1324
1325     if (ir.ePull != epullUMBRELLA)
1326         gmx_fatal(FARGS,"This is not a tpr of an umbrella simulation. Found ir.ePull = %s\n",
1327                 epull_names[ir.ePull]);
1328
1329     /* nr of pull groups */
1330     ngrp=ir.pull->ngrp;
1331     if (ngrp < 1)
1332         gmx_fatal(FARGS,"This is not a tpr of umbrella simulation. Found only %d pull groups\n",ngrp);
1333
1334     header->npullgrps=ir.pull->ngrp;
1335     header->pull_geometry=ir.pull->eGeom;
1336     copy_ivec(ir.pull->dim,header->pull_dim);
1337     header->pull_ndim=header->pull_dim[0]+header->pull_dim[1]+header->pull_dim[2];
1338     if (header->pull_geometry==epullgPOS && header->pull_ndim>1)
1339         gmx_fatal(FARGS,"Found pull geometry 'position' and more than 1 pull dimension (%d).\n"
1340                 "Hence, the pull potential does not correspond to a one-dimensional umbrella potential.\n"
1341                 "If you have some special umbrella setup you may want to write your own pdo files\n"
1342                 "and feed them into g_wham. Check g_wham -h !\n",header->pull_ndim);
1343     snew(header->k,ngrp);
1344     snew(header->init_dist,ngrp);
1345     snew(header->umbInitDist,ngrp);
1346
1347     for (i=0;i<ngrp;i++)
1348     {
1349         header->k[i]=ir.pull->grp[i+1].k;
1350         if (header->k[i]==0.0)
1351             gmx_fatal(FARGS,"Pull group %d has force constant of of 0.0 in %s.\n"
1352                     "That doesn't seem to be an Umbrella tpr.\n",
1353                     i,fn);
1354         copy_rvec(ir.pull->grp[i+1].init,header->init_dist[i]);
1355         header->Flipped[i]=opt->bFlipProf;
1356
1357         /* initial distance to reference */
1358         switch(header->pull_geometry)
1359         {
1360         case epullgPOS:
1361             for (d=0;d<DIM;d++)
1362                 if (header->pull_dim[d])
1363                     header->umbInitDist[i]=header->init_dist[i][d];
1364             break;
1365         case epullgDIST:
1366             header->umbInitDist[i]=header->init_dist[i][0];
1367             break;
1368         default:
1369             gmx_fatal(FARGS,"Pull geometry %s not supported\n",epullg_names[header->pull_geometry]);
1370         }
1371     }
1372
1373     if (opt->verbose || first)
1374     {
1375         printf("File %s, %d groups, geometry \"%s\", dimensions [%s %s %s], (%d dimensions)\n",
1376                 fn,header->npullgrps,epullg_names[header->pull_geometry],
1377                 int2YN(header->pull_dim[0]),int2YN(header->pull_dim[1]),int2YN(header->pull_dim[2]),
1378                 header->pull_ndim);
1379         for (i=0;i<ngrp;i++)
1380             printf("\tgrp %d) k = %.3f  inittial distance = %g\n",i,header->k[i],header->umbInitDist[i]);
1381     }
1382     first=0;
1383 }
1384
1385
1386 double dist_ndim(double **dx,int ndim,int line)
1387 {
1388     int i;
1389     double r2=0.;
1390
1391
1392     for (i=0;i<ndim;i++)
1393         r2+=sqr(dx[i][line]);
1394
1395     return sqrt(r2);
1396 }
1397
1398
1399 void read_pull_xf(const char *fn, const char *fntpr,
1400                   t_UmbrellaHeader * header, t_UmbrellaWindow * window,
1401                   t_UmbrellaOptions *opt, gmx_bool bGetMinMax,real *mintmp,
1402                   real *maxtmp)
1403 {
1404     double **y,pos=0.,t,force,time0=0.,dt;
1405     int ny,nt,bins,ibin,i,g,dstep=1,nColPerGrp,nColRef,nColExpect;
1406     real min,max,minfound,maxfound;
1407     gmx_bool dt_ok,timeok,bHaveForce;
1408     const char *quantity;
1409
1410         minfound=1e20;
1411         maxfound=-1e20;
1412
1413     /*
1414      in force    output pullf.xvg: No   reference, one  column  per pull group
1415      in position output pullx.xvg: ndim reference, ndim columns per pull group
1416      */
1417     nColPerGrp = opt->bPullx ? header->pull_ndim : 1;
1418     nColRef    = opt->bPullx ? header->pull_ndim : 0;
1419     quantity   = opt->bPullx ? "position" : "force";
1420     nColExpect = 1 + nColRef + header->npullgrps*nColPerGrp;
1421     bHaveForce = opt->bPullf;
1422
1423     nt=read_xvg(fn,&y,&ny);
1424
1425     /* Check consistency */
1426     if (nt<1)
1427         gmx_fatal(FARGS,"Empty pull %s file %s\n",quantity,fn);
1428     if (ny != nColExpect)
1429         gmx_fatal(FARGS,"Found %d pull groups in %s,\n but %d data columns in %s (expected %d)\n",
1430                 header->npullgrps,fntpr,ny-1,fn,nColExpect-1);
1431
1432     if (opt->verbose)
1433         printf("Found %d times and %d %s sets %s\n",nt,(ny-1)/nColPerGrp,quantity,fn);
1434
1435     if (!bGetMinMax)
1436     {
1437         bins=opt->bins;
1438         min=opt->min;
1439         max=opt->max;
1440         if (nt>1)
1441             window->dt=y[0][1]-y[0][0];
1442         else if (opt->nBootStrap && opt->dtBootStrap!=0.0)
1443             fprintf(stderr,"\n *** WARNING, Could not determine time step in %s\n",fn);
1444
1445         /* Need to alocate memory and set up structure */
1446         window->nPull=header->npullgrps;
1447         window->nBin=bins;
1448
1449         snew(window->Histo,window->nPull);
1450         snew(window->z,window->nPull);
1451         snew(window->k,window->nPull);
1452         snew(window->pos,window->nPull);
1453         snew(window->Flipped,window->nPull);
1454         snew(window->N, window->nPull);
1455         snew(window->Ntot, window->nPull);
1456
1457         for(g=0;g<window->nPull;++g)
1458         {
1459             window->z[g]=1;
1460             snew(window->Histo[g],bins);
1461             window->k[g]=header->k[g];
1462             window->Flipped[g]=header->Flipped[g];
1463             window->N[g]=0;
1464             window->Ntot[g]=0;
1465             window->pos[g]=header->umbInitDist[g];
1466         }
1467     }
1468     else
1469     {
1470         /* only determine min and max */
1471         minfound=1e20;
1472         maxfound=-1e20;
1473         min=max=bins=0; /* Get rid of warnings */
1474     }
1475
1476     if(header->Flipped[0])
1477         gmx_fatal(FARGS,"Sorry, flipping not supported for gmx4 output\n");
1478
1479     for (i=0;i<nt;i++)
1480     {
1481         /* Do you want that time frame? */
1482       t=1.0/1000*((int)(0.5+y[0][i]*1000)); /* round time to fs */
1483
1484         /* get time step of pdo file and get dstep from opt->dt */
1485         if (i==0)
1486             time0=t;
1487         else if (i==1)
1488         {
1489             dt=t-time0;
1490             if (opt->dt>0.0)
1491             {
1492                 dstep=(int)(opt->dt/dt+0.5);
1493                 if (dstep==0)
1494                     dstep=1;
1495             }
1496             if (!bGetMinMax)
1497                 window->dt=dt*dstep;
1498         }
1499
1500         dt_ok=(i%dstep == 0);
1501         timeok=(dt_ok && t >= opt->tmin && t <= opt->tmax);
1502         /*if (opt->verbose)
1503       printf(" time = %f, (tmin,tmax)=(%e,%e), dt_ok=%d timeok=%d\n",
1504       t,opt->tmin, opt->tmax, dt_ok,timeok); */
1505
1506         if (timeok)
1507         {
1508             for(g=0;g<header->npullgrps;++g)
1509             {
1510                 if (bHaveForce)
1511                 {
1512                     /* y has 1 time column y[0] and one column per force y[1],...,y[nGrps] */
1513                     force=y[g+1][i];
1514                     pos= -force/header->k[g] + header->umbInitDist[g];
1515                 }
1516                 else
1517                 {
1518                     switch (header->pull_geometry)
1519                     {
1520                     case epullgDIST:
1521                         /* y has 1 time column y[0] and nColPerGrps columns per pull group;
1522                                Distance to reference: */
1523                         pos=dist_ndim(y+1+nColRef+g*nColPerGrp,header->pull_ndim,i);
1524                         break;
1525                     case epullgPOS:
1526                         /* with geometry==position, we have always one column per group;
1527                                Distance to reference: */
1528                         pos=y[1+nColRef+g][i];
1529                         break;
1530                     default:
1531                         gmx_fatal(FARGS,"Bad error, this error should have been catched before. Ups.\n");
1532                     }
1533                 }
1534
1535                 /* printf("grp %d dpos %f poseq %f pos %f \n",g,dpos,poseq,pos); */
1536                 if (bGetMinMax)
1537                 {
1538                     if (pos<minfound)
1539                         minfound=pos;
1540                     if (pos>maxfound)
1541                         maxfound=pos;
1542                 }
1543                 else
1544                 {
1545                     ibin=(int) floor((pos-min)/(max-min)*bins);
1546                     if (opt->cycl==enCycl_yes)
1547                     {
1548                         if (ibin<0)
1549                             while ( (ibin+=bins) < 0);
1550                         else if (ibin>=bins)
1551                             while ( (ibin-=bins) >= bins);
1552                     }
1553                     if(ibin >= 0 && ibin < bins)
1554                     {
1555                         window->Histo[g][ibin]+=1.;
1556                         window->N[g]++;
1557                     }
1558                     window->Ntot[g]++;
1559                 }
1560             }
1561         }
1562         else if (t>opt->tmax)
1563         {
1564             if (opt->verbose)
1565                 printf("time %f larger than tmax %f, stop reading this pullx/pullf file\n",t,opt->tmax);
1566             break;
1567         }
1568     }
1569
1570     if (bGetMinMax)
1571     {
1572         *mintmp=minfound;
1573         *maxtmp=maxfound;
1574     }
1575 }
1576
1577
1578 void read_tpr_pullxf_files(char **fnTprs,char **fnPull,
1579                            int nfiles, t_UmbrellaHeader* header,
1580                            t_UmbrellaWindow **window, t_UmbrellaOptions *opt)
1581 {
1582     int i;
1583     real mintmp,maxtmp;
1584
1585
1586     printf("Reading %d tpr and pullf files\n",nfiles);
1587
1588     /* min and max not given? */
1589     if (opt->bAuto)
1590     {
1591         printf("Automatic determination of boundaries...\n");
1592         opt->min=1e20;
1593         opt->max=-1e20;
1594         for (i=0;i<nfiles; i++)
1595         {
1596             if (whaminFileType(fnTprs[i]) != whamin_tpr)
1597                 gmx_fatal(FARGS,"Expected the %d'th file in input file to be a tpr file\n",i);
1598             read_tpr_header(fnTprs[i],header,opt);
1599             if (whaminFileType(fnPull[i]) != whamin_pullxf)
1600                 gmx_fatal(FARGS,"Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n",i);
1601             read_pull_xf(fnPull[i],fnTprs[i],header,NULL,opt,TRUE,&mintmp,&maxtmp);
1602             if (maxtmp>opt->max)
1603                 opt->max=maxtmp;
1604             if (mintmp<opt->min)
1605                 opt->min=mintmp;
1606         }
1607         printf("\nDetermined boundaries to %f and %f\n\n",opt->min,opt->max);
1608         if (opt->bBoundsOnly){
1609             printf("Found option -boundsonly, now exiting.\n");
1610             exit (0);
1611         }
1612     }
1613     /* store stepsize in profile */
1614     opt->dz=(opt->max-opt->min)/opt->bins;
1615
1616     snew(*window,nfiles);
1617     for (i=0;i<nfiles; i++)
1618     {
1619         if (whaminFileType(fnTprs[i]) != whamin_tpr)
1620             gmx_fatal(FARGS,"Expected the %d'th file in input file to be a tpr file\n",i);
1621         read_tpr_header(fnTprs[i],header,opt);
1622         if (whaminFileType(fnPull[i]) != whamin_pullxf)
1623             gmx_fatal(FARGS,"Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n",i);
1624         read_pull_xf(fnPull[i],fnTprs[i],header,*window+i,opt,FALSE,NULL,NULL);
1625     }
1626 }
1627
1628
1629 int gmx_wham(int argc,char *argv[])
1630 {
1631     const char *desc[] = {
1632             "This is an analysis program that implements the Weighted",
1633             "Histogram Analysis Method (WHAM). It is intended to analyze",
1634             "output files generated by umbrella sampling simulations to ",
1635             "compute a potential of mean force (PMF). [PAR]",
1636             "At present, three input modes are supported:[BR]",
1637             "[TT]*[tt] With option [TT]-it[tt], the user provides a file which contains the",
1638             "  filenames of the umbrella simulation run-input files (tpr files),",
1639             "  AND, with option -ix, a file which contains filenames of",
1640             "  the pullx mdrun output files. The tpr and pullx files must",
1641             "  be in corresponding order, i.e. the first tpr created the",
1642             "  first pullx, etc.[BR]",
1643             "[TT]*[tt] Same as the previous input mode, except that the the user",
1644             "  provides the pull force ouput file names (pullf.xvg) with option -if.",
1645             "  From the pull force the position in the ubrella potential is",
1646             "  computed. This does not work with tabulated umbrella potentials.",
1647             "[TT]*[tt] With option [TT]-ip[tt], the user provides filenames of (gzipped) pdo files, i.e.",
1648             "  the gromacs 3.3 umbrella output files. If you have some unusual",
1649             "  reaction coordinate you may also generate your own pdo files and",
1650             "  feed them with the -ip option into to g_wham. The pdo file header",
1651             "  must be similar to the folowing:[BR]",
1652             "[TT]# UMBRELLA      3.0[BR]",
1653             "# Component selection: 0 0 1[BR]",
1654             "# nSkip 1[BR]",
1655             "# Ref. Group 'TestAtom'[BR]",
1656             "# Nr. of pull groups 2[BR]",
1657             "# Group 1 'GR1'  Umb. Pos. 5.0 Umb. Cons. 1000.0[BR]",
1658             "# Group 2 'GR2'  Umb. Pos. 2.0 Umb. Cons. 500.0[BR]",
1659             "#####[tt][BR]",
1660             "  Nr of pull groups, umbrella positions, force constants, and names",
1661             "  may (of course) differ. Following the header, a time column and",
1662             "  a data columns for each pull group follow (i.e. the displacement",
1663             "  with respect to the umbrella center). Up to four pull groups are possible",
1664             "  at present.[PAR]",
1665             "By default, the output files are[BR]",
1666             "  [TT]-o[tt]      PMF output file[BR]",
1667             "  [TT]-hist[tt]   histograms output file[PAR]",
1668             "The umbrella potential is assumed to be harmonic and the force constants are ",
1669             "read from the tpr or pdo files. If a non-harmonic umbrella force was applied ",
1670             "a tabulated potential can be provied with -tab.[PAR]",
1671             "WHAM OPTIONS[PAR]",
1672             "  [TT]-bins[tt]   Nr of bins used in analysis[BR]",
1673             "  [TT]-temp[tt]   Temperature in the simulations[BR]",
1674             "  [TT]-tol[tt]    Stop iteration if profile (probability) changed less than tolerance[BR]",
1675             "  [TT]-auto[tt]   Automatic determination of boudndaries[BR]",
1676             "  [TT]-min,-max[tt]   Boundaries of the profile [BR]",
1677             "The data points which are used ",
1678             "to compute the profile can be restricted with options -b, -e, and -dt. ",
1679             "Play particularly with -b to ensure sufficient equilibration in each ",
1680             "umbrella window![PAR]",
1681             "With -log (default) the profile is written in energy units, otherwise (-nolog) as ",
1682             "probability. The unit can be specified with -unit. With energy output, ",
1683             "the energy in the first bin is defined to be zero. If you want the free energy at a different ",
1684             "position to be zero, choose with -zprof0 (useful with bootstrapping, see below).[PAR]",
1685             "For cyclic (or periodic) reaction coordinates (dihedral angle, channel PMF",
1686             "without osmotic gradient), -cycl is useful.[BR]",
1687             "[TT]-cycl yes[tt]        min and max are assumed to",
1688             "be neighboring points and histogram points outside min and max are mapped into ",
1689             "the interval [min,max] (compare histogram output). [BR]",
1690             "[TT]-cycl weighted[tt]   First, a non-cyclic profile is computed. Subsequently, ",
1691             "periodicity is enforced by adding corrections dG(i) between neighboring bins",
1692             "i and i+1. The correction is chosen proportional to 1/[n(i)*n(i+1)]^alpha, where",
1693             "n(i) denotes the total nr of data points in bin i as collected from all histograms.",
1694             "alpha is defined with -alpha. The corrections are written to the file defined by -wcorr.",
1695             " (Compare Hub and de Groot, PNAS 105:1198 (2008))[PAR]",
1696             "ERROR ANALYSIS[BR]",
1697             "Statistical errors may be estimated with bootstrap analysis. Use it with care, ",
1698             "otherwise the statistical error may be substantially undererstimated !![BR]",
1699             "[TT]-nBootstrap[tt] defines the nr of bootstraps. Two bootstrapping modes are supported.[BR]",
1700             "[TT]-histbs[tt]    Complete histograms are considered as independent data points (default). For each",
1701             "bootstrap, N histograms are randomly chosen from the N given histograms (allowing duplication).",
1702             "To avoid gaps without data along the reaction coordinate blocks of histograms (-histbs-block)",
1703             "may be defined. In that case, the given histograms are divided into blocks and ",
1704             "only histograms within each block are mixed. Note that the histograms",
1705             "within each block must be representative for all possible histograms, otherwise the",
1706             "statistical error is undererstimated![BR]",
1707             "[TT]-nohistbs[tt]  The given histograms are used to generate new random histograms,",
1708             "such that the generated data points are distributed according the given histograms. The number",
1709             "of points generated for each bootstrap histogram can be controlled with -bs-dt.",
1710             "Note that one data point should be generated for each *independent* point in the given",
1711             "histograms. With the long autocorrelations in MD simulations, this procedure may ",
1712             "easily understimate the error![BR]",
1713             "Bootstrapping output:[BR]",
1714             "[TT]-bsres[tt]   Average profile and standard deviations[BR]",
1715             "[TT]-bsprof[tt]  All bootstrapping profiles[BR]",
1716             "With [TT]-vbs[tt] (verbose bootstrapping), the histograms of each bootstrap are written, and, ",
1717             "with [TT]-nohistBS[tt], the cummulants of the histogram.",
1718     };
1719
1720     static t_UmbrellaOptions opt;
1721     static gmx_bool bHistOnly=FALSE;
1722
1723     const char *en_unit[]={NULL,"kJ","kCal","kT",NULL};
1724     const char *en_unit_label[]={"","E (kJ mol\\S-1\\N)","E (kcal mol\\S-1\\N)","E (kT)",};
1725     const char *en_cycl[]={NULL,"no","yes","weighted",NULL};
1726
1727     t_pargs pa[] = {
1728             { "-min", FALSE, etREAL, {&opt.min},
1729               "Minimum coordinate in profile"},
1730             { "-max", FALSE, etREAL, {&opt.max},
1731               "Maximum coordinate in profile"},
1732             { "-auto", FALSE, etBOOL, {&opt.bAuto},
1733               "determine min and max automatically"},
1734             { "-bins",FALSE, etINT, {&opt.bins},
1735               "Number of bins in profile"},
1736             { "-temp", FALSE, etREAL, {&opt.Temperature},
1737               "Temperature"},
1738             { "-tol", FALSE, etREAL, {&opt.Tolerance},
1739               "Tolerance"},
1740             { "-v", FALSE, etBOOL, {&opt.verbose},
1741               "verbose mode"},
1742             { "-b", FALSE, etREAL, {&opt.tmin},
1743               "first time to analyse (ps)"},
1744             { "-e", FALSE, etREAL, {&opt.tmax},
1745               "last time to analyse (ps)"},
1746             { "-dt", FALSE, etREAL, {&opt.dt},
1747               "Analyse only every dt ps"},
1748             { "-histonly", FALSE, etBOOL, {&bHistOnly},
1749               "Write histograms and exit"},
1750             { "-boundsonly", FALSE, etBOOL, {&opt.bBoundsOnly},
1751               "Determine min and max and exit (with -auto)"},
1752             { "-log", FALSE, etBOOL, {&opt.bLog},
1753               "Calculate the log of the profile before printing"},
1754             { "-unit", FALSE,  etENUM, {en_unit},
1755               "energy unit in case of log output" },
1756             { "-zprof0", FALSE, etREAL, {&opt.zProf0},
1757               "Define profile to 0.0 at this position (with -log)"},
1758             { "-cycl", FALSE, etENUM, {en_cycl},
1759               "Create cyclic/periodic profile. Assumes min and max are the same point."},
1760             { "-alpha", FALSE, etREAL, {&opt.alpha},
1761               "for '-cycl weighted', set parameter alpha"},
1762             { "-flip", FALSE, etBOOL, {&opt.bFlipProf},
1763               "Combine halves of profile (not supported)"},
1764             { "-hist-eq", FALSE, etBOOL, {&opt.bHistEq},
1765               "Enforce equal weight for all histograms. (Non-Weighed-HAM)"},
1766             { "-nBootstrap", FALSE,  etINT, {&opt.nBootStrap},
1767               "nr of bootstraps to estimate statistical uncertainty" },
1768             { "-bs-dt", FALSE, etREAL, {&opt.dtBootStrap},
1769               "timestep for synthetic bootstrap histograms (ps). Ensure independent data points!"},
1770             { "-bs-seed", FALSE, etINT, {&opt.bsSeed},
1771               "seed for bootstrapping. (-1 = use time)"},
1772             { "-histbs", FALSE, etBOOL, {&opt.bHistBootStrap},
1773               "In bootstrapping, consider complete histograms as one data point. "
1774               "Accounts better for long autocorrelations."},
1775             { "-histbs-block", FALSE, etINT, {&opt.histBootStrapBlockLength},
1776               "when mixin histograms only mix within blocks of -histBS_block."},
1777             { "-vbs", FALSE, etBOOL, {&opt.bs_verbose},
1778                 "verbose bootstrapping. Print the cummulants and a histogram file for each bootstrap."},
1779     };
1780
1781     t_filenm fnm[] = {
1782             { efDAT, "-ix","pullx-files",ffOPTRD},     /* wham input: pullf.xvg's and tprs    */
1783             { efDAT, "-if","pullf-files",ffOPTRD},     /* wham input: pullf.xvg's and tprs    */
1784             { efDAT, "-it","tpr-files",ffOPTRD},       /* wham input: tprs                    */
1785             { efDAT, "-ip","pdo-files",ffOPTRD},       /* wham input: pdo files (gmx3 style)  */
1786             { efXVG, "-o", "profile", ffWRITE },       /* output file for profile */
1787             { efXVG, "-hist","histo", ffWRITE},        /* output file for histograms */
1788             { efXVG, "-bsres","bsResult", ffOPTWR},    /* average and errors of bootstrap analysis */
1789             { efXVG, "-bsprof","bsProfs", ffOPTWR},    /* output file for bootstrap profiles       */
1790             { efDAT, "-tab","umb-pot",ffOPTRD},        /* Tabulated umbrella potential (if not harmonic) */
1791             { efXVG, "-wcorr","cycl-corr",ffOPTRD},    /* Corrections to profile in case -cycl weighted */
1792     };
1793
1794     int i,j,l,nfiles,nwins,nfiles2;
1795     t_UmbrellaHeader header;
1796     t_UmbrellaWindow * window=NULL;
1797     double *profile,maxchange=1e20;
1798     gmx_bool bMinSet,bMaxSet,bAutoSet,bExact=FALSE;
1799     char **fninTpr,**fninPull,**fninPdo;
1800     const char *fnPull;
1801     FILE *histout,*profout;
1802     char ylabel[256],title[256];
1803     output_env_t oenv;
1804
1805     opt.bins=200;
1806     opt.verbose=FALSE;
1807     opt.cycl=enCycl_no;
1808     opt.tmin=50;
1809     opt.tmax=1e20;
1810     opt.dt=0.0;
1811     opt.bShift=TRUE;
1812     opt.nBootStrap=0;
1813     opt.dtBootStrap=0.0;
1814     opt.bsSeed=-1;
1815     opt.bHistBootStrap=TRUE;
1816     opt.histBootStrapBlockLength=12;
1817     opt.zProfZero=0.0;
1818     opt.bWeightedCycl=FALSE;
1819     opt.alpha=2;
1820     opt.bHistOutOnly=FALSE;
1821     opt.min=0;
1822     opt.max=0;
1823     opt.bLog=TRUE;
1824     opt.unit=en_kJ;
1825     opt.zProf0=0.0;
1826     opt.nBootStrap=0;
1827     opt.bsSeed=-1;
1828     opt.bHistBootStrap=TRUE;
1829     opt.histBootStrapBlockLength=8;
1830     opt.bs_verbose=FALSE;
1831     opt.Temperature=298;
1832     opt.bFlipProf=FALSE;
1833     opt.Tolerance=1e-6;
1834     opt.bAuto=TRUE;
1835     opt.bBoundsOnly=FALSE;
1836
1837
1838 #define NFILE asize(fnm)
1839
1840     CopyRight(stderr,argv[0]);
1841     parse_common_args(&argc,argv,PCA_BE_NICE,
1842             NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
1843
1844     opt.unit=nenum(en_unit);
1845     opt.cycl=nenum(en_cycl);
1846
1847     opt.bProf0Set=opt2parg_bSet("-zprof0",  asize(pa), pa);
1848
1849     opt.bTab=opt2bSet("-tab",NFILE,fnm);
1850     opt.bPdo=opt2bSet("-ip",NFILE,fnm);
1851     opt.bTpr=opt2bSet("-it",NFILE,fnm);
1852     opt.bPullx=opt2bSet("-ix",NFILE,fnm);
1853     opt.bPullf=opt2bSet("-if",NFILE,fnm);
1854     if  (opt.bTab && opt.bPullf)
1855         gmx_fatal(FARGS,"Force input does not work with tabulated potentials. "
1856                 "Provide pullx.xvg or pdo files!");
1857
1858 #define BOOLXOR(a,b) ( ((!(a))&&(b)) || ((a)&&(!(b))))
1859     if (!opt.bPdo && !BOOLXOR(opt.bPullx,opt.bPullf))
1860         gmx_fatal(FARGS,"Give either pullx (-ix) OR pullf (-if) data. Not both.");
1861     if ( !opt.bPdo && !(opt.bTpr || opt.bPullf || opt.bPullx))
1862         gmx_fatal(FARGS,"g_wham supports three input modes, pullx, pullf, or pdo file input."
1863                 "\n\n Check g_wham -h !");
1864
1865     opt.fnPdo=opt2fn("-ip",NFILE,fnm);
1866     opt.fnTpr=opt2fn("-it",NFILE,fnm);
1867     opt.fnPullf=opt2fn("-if",NFILE,fnm);
1868     opt.fnPullx=opt2fn("-ix",NFILE,fnm);
1869
1870     bMinSet = opt2parg_bSet("-min",  asize(pa), pa);
1871     bMaxSet = opt2parg_bSet("-max",  asize(pa), pa);
1872     bAutoSet = opt2parg_bSet("-auto",  asize(pa), pa);
1873     if ( (bMinSet || bMaxSet) && bAutoSet)
1874         gmx_fatal(FARGS,"With -auto, do not give -min or -max\n");
1875
1876     if ( (bMinSet && !bMaxSet) || (!bMinSet && bMaxSet))
1877         gmx_fatal(FARGS,"When giving -min, you must give -max (and vice versa), too\n");
1878
1879     if (bMinSet && opt.bAuto)
1880     {
1881         printf("Note: min and max given, switching off -auto.\n");
1882         opt.bAuto=FALSE;
1883     }
1884
1885     /* Reading gmx4 pull output and tpr files */
1886     if (opt.bTpr || opt.bPullf || opt.bPullx)
1887     {
1888         read_wham_in(opt.fnTpr,&fninTpr,&nfiles,&opt);
1889
1890         fnPull=opt.bPullf ? opt.fnPullf : opt.fnPullx;
1891         read_wham_in(fnPull,&fninPull,&nfiles2,&opt);
1892         printf("Found %d tpr and %d pull %s files in %s and %s, respectively\n",
1893                 nfiles,nfiles2,opt.bPullf ? "force" : "position",opt.fnTpr,fnPull);
1894         if (nfiles!=nfiles2)
1895             gmx_fatal(FARGS,"Found %d file names in %s, but %d in %s\n",nfiles,
1896                     opt.fnTpr,nfiles2,fnPull);
1897         read_tpr_pullxf_files(fninTpr,fninPull,nfiles, &header, &window, &opt);
1898     }
1899     else
1900     {
1901         /* reading pdo files */
1902         read_wham_in(opt.fnPdo,&fninPdo,&nfiles,&opt);
1903         printf("Found %d pdo files in %s\n",nfiles,opt.fnPdo);
1904         read_pdo_files(fninPdo,nfiles, &header, &window, &opt);
1905     }
1906     nwins=nfiles;
1907
1908     /* enforce equal weight for all histograms? */
1909     if (opt.bHistEq)
1910         enforceEqualWeights(window,nwins);
1911
1912     /* write histograms */
1913     histout=xvgropen(opt2fn("-hist",NFILE,fnm),"Umbrella histograms",
1914             "z","count",oenv);
1915     for(l=0;l<opt.bins;++l)
1916     {
1917         fprintf(histout,"%e\t",(double)(l+0.5)/opt.bins*(opt.max-opt.min)+opt.min);
1918         for(i=0;i<nwins;++i)
1919         {
1920             for(j=0;j<window[i].nPull;++j)
1921             {
1922                 fprintf(histout,"%e\t",window[i].Histo[j][l]);
1923             }
1924         }
1925         fprintf(histout,"\n");
1926     }
1927     ffclose(histout);
1928     if (bHistOnly)
1929     {
1930         printf("Wrote histograms to %s, now exiting.\n",opt2fn("-hist",NFILE,fnm));
1931         return 0;
1932     }
1933
1934     /* Using tabulated umbrella potential */
1935     if (opt.bTab)
1936         setup_tab(opt2fn("-tab",NFILE,fnm),&opt);
1937
1938     setup_acc_wham(window,nwins,&opt);
1939
1940     /* Calculate profile */
1941     snew(profile,opt.bins);
1942     opt.stepchange=(opt.verbose)? 1 : 100;
1943     i=0;
1944     do {
1945         if (maxchange<opt.Tolerance)
1946         {
1947             bExact=TRUE;
1948             /* if (opt.verbose) */
1949             printf("Switched to exact iteration in iteration %d\n",i);
1950         }
1951         calc_profile(profile,window,nwins,&opt,bExact);
1952         if (((i%opt.stepchange) == 0 || i==1) && !i==0)
1953             printf("\t%4d) Maximum change %e\n",i,maxchange);
1954         i++;
1955     } while( (maxchange=calc_z(profile, window, nwins, &opt,bExact)) > opt.Tolerance || !bExact);
1956     printf("Converged in %d iterations. Final maximum change %g\n",i,maxchange);
1957
1958     /* Write profile in energy units? */
1959     if (opt.bLog)
1960     {
1961         prof_normalization_and_unit(profile,&opt);
1962         strcpy(ylabel,en_unit_label[opt.unit]);
1963         strcpy(title,"Umbrella potential");
1964     }
1965     else
1966     {
1967         strcpy(ylabel,"Density of states");
1968         strcpy(title,"Density of states");
1969     }
1970     /* Force cyclic profile by wheighted correction? */
1971     if (opt.cycl==enCycl_weighted)
1972         cyclicProfByWeightedCorr(profile,window,nwins,&opt,TRUE,
1973                                  opt2fn("-wcorr",NFILE,fnm),oenv);
1974
1975     profout=xvgropen(opt2fn("-o",NFILE,fnm),title,"z",ylabel,oenv);
1976     for(i=0;i<opt.bins;++i)
1977         fprintf(profout,"%e\t%e\n",
1978                 (double)(i+0.5)/opt.bins*(opt.max-opt.min)+opt.min,profile[i]);
1979     ffclose(profout);
1980     printf("Wrote %s\n",opt2fn("-o",NFILE,fnm));
1981
1982     /* Bootstrap Method */
1983     if (opt.nBootStrap)
1984         do_bootstrapping(opt2fn("-bsres",NFILE,fnm),opt2fn("-bsprof",NFILE,fnm),
1985                 opt2fn("-hist",NFILE,fnm),
1986                 ylabel, profile, window, nwins, &opt,oenv);
1987
1988     thanx(stderr);
1989     return 0;
1990 }