diff --git a/tools/lammps-gui/update-wham-2.0.11.patch b/tools/lammps-gui/update-wham-2.0.11.patch index 3b6e0e372a..b884c25336 100644 --- a/tools/lammps-gui/update-wham-2.0.11.patch +++ b/tools/lammps-gui/update-wham-2.0.11.patch @@ -55,8 +55,79 @@ index 0000000..b4f0fe6 + TYPE DOC + PERMISSIONS OWNER_READ GROUP_READ WORLD_READ +) +diff --git a/nr/locate.c b/nr/locate.c +index 9f92dc0..f3bf294 100644 +--- a/nr/locate.c ++++ b/nr/locate.c +@@ -11,7 +11,7 @@ void locate(double xx[], int n, double x, int *j) + ascnd=(xx[n] > xx[0]); // I think this makes it zero based + while (ju-jl > 1) { + jm=(ju+jl) >> 1; +- if (x > xx[jm] == ascnd) ++ if ((x > xx[jm]) == ascnd) + jl=jm; + else + ju=jm; +diff --git a/wham-2d/histogram.c b/wham-2d/histogram.c +index 1bd1329..b5d1c01 100644 +--- a/wham-2d/histogram.c ++++ b/wham-2d/histogram.c +@@ -78,14 +78,14 @@ return hp; + } + + /* Get a value from a histogram structure +- * Given i and j, the indices into the global histogram, return ++ * Given i and j, the indices into the global histogram, return + * the correct histogram value. Since we only store range of histogram +- * indices containing the nonzero values, we have to check the index value ++ * indices containing the nonzero values, we have to check the index value + * against the range in the structure. + */ + double get_histval(struct histogram *hist, int i, int j) + { +-if ( (i < hist->first_x) || (i > hist->last_x) ++if ( (i < hist->first_x) || (i > hist->last_x) + ||(j < hist->first_y) || (j > hist->last_y)) + { + return 0.0; +@@ -239,7 +239,7 @@ return error; + + + // Calculate the free energy, setting the minimum value to 0 +-void calc_free(double **free, double **prob, double kT, ++void calc_free(double **free, double **prob, double kT, + int use_mask, int **mask) + { + int i,j; +@@ -257,7 +257,14 @@ for (i=0; i 0.0) ++ { ++ free[i][j] = -kT * log(prob[i][j]); ++ } ++ else ++ { ++ free[i][j] = 0.0; ++ } + if (free[i][j] < min) + { + min = free[i][j]; +@@ -321,8 +328,6 @@ return 0.5*(springx*(dx*dx) + springy*(dy*dy)); + + void calc_coor(int i, int j, double *coor) + { +-coor[0] = HIST_MINx + BIN_WIDTHx*((double)i+0.5); +-coor[1] = HIST_MINy + BIN_WIDTHy*((double)j+0.5); ++coor[0] = HIST_MINx + BIN_WIDTHx*((double)i+0.5); ++coor[1] = HIST_MINy + BIN_WIDTHy*((double)j+0.5); + } +- +- diff --git a/wham-2d/wham-2d.c b/wham-2d/wham-2d.c -index fb6e059..2c5594f 100644 +index fb6e059..a6b8483 100644 --- a/wham-2d/wham-2d.c +++ b/wham-2d/wham-2d.c @@ -25,7 +25,7 @@ @@ -77,6 +148,15 @@ index fb6e059..2c5594f 100644 int main(int argc, char *argv[]) { +@@ -57,7 +57,7 @@ double sum; + int iteration; + int max_iteration = 100000; + int numpad; +-int **mask; ++int **mask = NULL; + int use_mask; + + cpu1 = ((double) clock())/CLOCKS_PER_SEC; @@ -76,6 +76,61 @@ for (i=0; i 0.0) ++ { ++ free[i] = -kT * log(prob[i]); ++ } ++ else ++ { ++ free[i] = 0.0; ++ } ++ ++ if (free[i] < min) + { + min = free[i]; + bin_min = i; +@@ -252,5 +260,5 @@ return 0.5*dx*dx*spring; + + double calc_coor(int i) + { +-return HIST_MIN + BIN_WIDTH*((double)i+0.5); ++return HIST_MIN + BIN_WIDTH*((double)i+0.5); + } diff --git a/wham/wham.c b/wham/wham.c -index 487871b..1496eed 100644 +index 487871b..edb8125 100644 --- a/wham/wham.c +++ b/wham/wham.c +@@ -1,4 +1,4 @@ +-/* ++/* + * WHAM -- Weighted Histogram Analysis Method + * + * Reference 1: Computer Physics Communications, 91(1995) 275-282, Benoit Roux @@ -21,7 +21,7 @@ #include "wham.h" @@ -184,6 +458,15 @@ index 487871b..1496eed 100644 int main(int argc, char *argv[]) { +@@ -41,7 +42,7 @@ int first; + int bin_min; + int have_energy; + char *freefile; +-FILE *METAFILE, *FREEFILE; ++FILE *METAFILE, *FREEFILE; + struct hist_group *hist_group; + struct histogram *hp; + double coor; @@ -82,6 +83,61 @@ for (i=0; inum_windows); + // Figure out if we have trajectories at different temperatures. + // Missing temperatures are set to -1 in read_metadata, and + // since we require that either all trajectories specify a temperature +-// or all trajectories are assumed to be at the wham temperature, we only ++// or all trajectories are assumed to be at the wham temperature, we only + // have to check one of them + if (hist_group->kT[0] > 0) + { +@@ -257,7 +313,7 @@ if (hist_group->kT[0] > 0) + else + { + have_energy = 0; +- for (i=0; i< hist_group->num_windows; i++) ++ for (i=0; i< hist_group->num_windows; i++) + { + hist_group->kT[i] = kT; + } +@@ -269,7 +325,7 @@ if (!final_f) + { + printf("couldn't allocate space for final_f: %s\n", strerror(errno)); + exit(errno); +- } ++ } + + free(HISTOGRAM); + +@@ -305,7 +361,8 @@ while (! is_converged(hist_group) || first ) + for (i=0; i< NUM_BINS; i++) + { + coor = calc_coor(i); +- printf("%f\t%f\t%f\n", coor, free_ene[i], prob[i]); ++ if (prob[i] != 0.0) ++ printf("%f\t%f\t%f\n", coor, free_ene[i], prob[i]); + } + printf("\n"); + +@@ -319,7 +376,7 @@ while (! is_converged(hist_group) || first ) + } + } + // Cheesy bailout if we're going on too long +- if (iteration >= max_iteration) ++ if (iteration >= max_iteration) + { + printf("Too many iterations: %d\n", iteration); + break; +@@ -383,11 +440,11 @@ for (i=0; i< num_mc_runs; i++) + //printf("Faking %d: %d %d\n", i,j,hp->num_points); + num_used = hp->last - hp->first + 1; + mk_new_hist(hp->cum, hp->data, num_used, hp->num_mc_samples, &idum); +- ++ + hist_group->F_old[j] = 0.0; + hist_group->F[j] = 0.0; + } +- ++ + // perform WHAM iterations on the fake data sets + iteration = 0; + first = 1; +@@ -403,7 +460,7 @@ for (i=0; i< num_mc_runs; i++) + printf("Too many iterations: %d\n", iteration); + break; + } +- } ++ } + printf("#MC trial %d: %d iterations\n", i, iteration); + printf("#PMF values\n"); + // accumulate the average and stdev of the resulting probabilities +@@ -419,18 +476,19 @@ for (i=0; i< num_mc_runs; i++) + for (j=0; j < NUM_BINS; j++) + { + pdf = -kT*log(prob[j]); +- ++ + ave_p[j] += prob[j]; + ave_pdf[j] += pdf; + ave_p2[j] += prob[j] * prob[j]; + ave_pdf2[j] += pdf*pdf; + } +- for (j=0; jnum_windows;j++) ++ for (j=0; jnum_windows;j++) + { + ave_F[j] += hist_group->F[j] - hist_group->F[0]; +- ave_F2[j] += hist_group->F[j]*hist_group->F[j] ; ++ ave_F2[j] += hist_group->F[j]*hist_group->F[j] ; + } +- } ++ } ++ if (num_mc_runs == 0) num_mc_runs = 1; + for (i=0; i < NUM_BINS; i++) + { + ave_p[i] /= (double)num_mc_runs; +@@ -457,12 +515,12 @@ if (!FREEFILE) + for (i=0; i< NUM_BINS; i++) + { + coor = calc_coor(i); +- printf("%f\t%f\t%f\t%f\t%f\n", coor, free_ene[i], ave_pdf2[i], +- prob[i], ave_p2[i]); ++ if (prob[i] != 0.0) ++ printf("%f\t%f\t%f\t%f\t%f\n", coor, free_ene[i], ave_pdf2[i], prob[i], ave_p2[i]); + } + for (i=0; inum_windows; i++) + { +- fprintf(FREEFILE,"%d\t%f\t%f\n", i, final_f[i],ave_F2[i]); ++ printf("%d\t%f\t%f\n", i, final_f[i],ave_F2[i]); + } + + exit(errno); +@@ -470,38 +528,37 @@ if (!FREEFILE) + else + { + // write out header +- fprintf(FREEFILE, "#Coor\t\tFree\t+/-\t\tProb\t\t+/-\n"); ++ fprintf(FREEFILE, "#Coor\t\tFree\t\t+/-\t\tProb\t\t+/-\n"); + // write out the leading padded values + for (i=-numpad; i<0; i++) + { + coor = calc_coor(i); +- fprintf(FREEFILE,"%f\t%f\t%f\t%f\t%f\n", coor, free_ene[NUM_BINS+i], +- ave_pdf2[NUM_BINS+i], +- final_prob[NUM_BINS+i], +- ave_p2[NUM_BINS+i]); ++ if (prob[NUM_BINS+1] != 0.0) ++ fprintf(FREEFILE,"%f\t%f\t%f\t%f\t%f\n", coor, free_ene[NUM_BINS+i], ++ ave_pdf2[NUM_BINS+i], final_prob[NUM_BINS+i], ave_p2[NUM_BINS+i]); + } + // write out the center values + for (i=0; inum_windows; i++) + { +- fprintf(FREEFILE,"#%d\t%f\t%f\n", i, final_f[i],ave_F2[i]); ++ fprintf(FREEFILE,"#%d\t\t%f\t%f\n", i, final_f[i],ave_F2[i]); + } + } + +@@ -515,7 +572,7 @@ exit(0); + /* + * Perform a single WHAM iteration + */ +-void wham_iteration(struct hist_group* hist_group, double *prob, ++void wham_iteration(struct hist_group* hist_group, double *prob, + int have_energy) + { + int i,j; +@@ -535,9 +592,9 @@ for (i=0; ihists[j]),i); + bias = calc_bias(hist_group,j,coor); + bf = exp((hist_group->F_old[j] - bias) / hist_group->kT[j]); +- /* If we have energies, then the histograms contain boltzmann +- * factors. If not, they contain integer counts. Accordingly, +- * we either use a partition function to normalize, or simply the ++ /* If we have energies, then the histograms contain boltzmann ++ * factors. If not, they contain integer counts. Accordingly, ++ * we either use a partition function to normalize, or simply the + * number of points. + */ + if (have_energy) diff --git a/wham/wham.h b/wham/wham.h index aacc1e8..7d509f2 100644 --- a/wham/wham.h