2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
45 #include "gromacs/commandline/pargs.h"
46 #include "gromacs/legacyheaders/typedefs.h"
47 #include "gromacs/utility/smalloc.h"
48 #include "gromacs/legacyheaders/macros.h"
49 #include "gromacs/utility/fatalerror.h"
50 #include "gromacs/math/vec.h"
51 #include "gromacs/utility/futil.h"
52 #include "gromacs/legacyheaders/readinp.h"
53 #include "gromacs/legacyheaders/txtdump.h"
55 #include "gromacs/fileio/xvgr.h"
56 #include "gromacs/math/units.h"
57 #include "gromacs/fileio/pdbio.h"
58 #include "gromacs/fileio/matio.h"
62 static int index2(int *ibox, int x, int y)
67 static int index3(int *ibox, int x, int y, int z)
69 return (ibox[2]*(ibox[1]*x+y)+z);
72 static gmx_int64_t indexn(int ndim, const int *ibox, const int *nxyz)
77 /* Compute index in 1-D array */
79 for (k = 0; (k < ndim); k++)
82 for (kk = k+1; (kk < ndim); kk++)
92 int Nx; /* x grid points in unit cell */
93 int Ny; /* y grid points in unit cell */
94 int Nz; /* z grid points in unit cell */
95 int dmin[3]; /* starting point x,y,z*/
96 int dmax[3]; /* ending point x,y,z */
97 real cell[6]; /* usual cell parameters */
101 static void lo_write_xplor(XplorMap * map, const char * file)
106 fp = gmx_ffopen(file, "w");
107 /* The REMARKS part is the worst part of the XPLOR format
108 * and may cause problems with some programs
110 fprintf(fp, "\n 2 !NTITLE\n");
111 fprintf(fp, " REMARKS Energy Landscape from GROMACS\n");
112 fprintf(fp, " REMARKS DATE: 2004-12-21 \n");
113 fprintf(fp, " %7d %7d %7d %7d %7d %7d %7d %7d %7d\n",
114 map->Nx, map->dmin[0], map->dmax[0],
115 map->Ny, map->dmin[1], map->dmax[1],
116 map->Nz, map->dmin[2], map->dmax[2]);
117 fprintf(fp, "%12.5E%12.5E%12.5E%12.5E%12.5E%12.5E\n",
118 map->cell[0], map->cell[1], map->cell[2],
119 map->cell[3], map->cell[4], map->cell[5]);
120 fprintf(fp, "ZYX\n");
123 for (n = 0; n < map->Nz; n++, z++)
125 fprintf(fp, "%8d\n", z);
126 for (i = 0; i < map->Nx*map->Ny; i += 6)
128 for (j = 0; j < 6; j++)
130 if (i+j < map->Nx*map->Ny)
132 fprintf(fp, "%12.5E", map->ed[n*map->Nx*map->Ny+i+j]);
138 fprintf(fp, " -9999\n");
142 static void write_xplor(const char *file, real *data, int *ibox, real dmin[], real dmax[])
151 snew(xm->ed, xm->Nx*xm->Ny*xm->Nz);
153 for (k = 0; (k < xm->Nz); k++)
155 for (j = 0; (j < xm->Ny); j++)
157 for (i = 0; (i < xm->Nx); i++)
159 xm->ed[n++] = data[index3(ibox, i, j, k)];
163 xm->cell[0] = dmax[XX]-dmin[XX];
164 xm->cell[1] = dmax[YY]-dmin[YY];
165 xm->cell[2] = dmax[ZZ]-dmin[ZZ];
166 xm->cell[3] = xm->cell[4] = xm->cell[5] = 90;
168 clear_ivec(xm->dmin);
169 xm->dmax[XX] = ibox[XX]-1;
170 xm->dmax[YY] = ibox[YY]-1;
171 xm->dmax[ZZ] = ibox[ZZ]-1;
173 lo_write_xplor(xm, file);
179 static void normalize_p_e(int len, double *P, int *nbin, real *E, real pmin)
184 for (i = 0; (i < len); i++)
192 printf("Ptot = %g\n", Ptot);
193 for (i = 0; (i < len); i++)
196 /* Have to check for pmin after normalizing to prevent "stretching"
211 static int comp_minima(const void *a, const void *b)
213 t_minimum *ma = (t_minimum *) a;
214 t_minimum *mb = (t_minimum *) b;
216 if (ma->ener < mb->ener)
220 else if (ma->ener > mb->ener)
231 void print_minimum(FILE *fp, int num, const t_minimum *min)
234 "Minimum %d at index " "%"GMX_PRId64 " energy %10.3f\n",
235 num, min->index, min->ener);
239 void add_minimum(FILE *fp, int num, const t_minimum *min, t_minimum *mm)
241 print_minimum(fp, num, min);
242 mm[num].index = min->index;
243 mm[num].ener = min->ener;
247 gmx_bool is_local_minimum_from_below(const t_minimum *this_min,
253 return ((dimension_index == dimension_min) ||
254 ((dimension_index > dimension_min) &&
255 (this_min->ener < W[neighbour_index])));
256 /* Note over/underflow within W cannot occur. */
260 gmx_bool is_local_minimum_from_above(const t_minimum *this_min,
266 return ((dimension_index == dimension_max) ||
267 ((dimension_index < dimension_max) &&
268 (this_min->ener < W[neighbour_index])));
269 /* Note over/underflow within W cannot occur. */
272 static void pick_minima(const char *logfile, int *ibox, int ndim, int len, real W[])
276 t_minimum *mm, this_min;
278 int loopmax, loopcounter;
282 fp = gmx_ffopen(logfile, "w");
283 /* Loop over each element in the array of dimenion ndim seeking
284 * minima with respect to every dimension. Specialized loops for
285 * speed with ndim == 2 and ndim == 3. */
289 /* This is probably impossible to reach anyway. */
292 for (i = 0; (i < ibox[0]); i++)
294 for (j = 0; (j < ibox[1]); j++)
296 /* Get the index of this point in the flat array */
297 this_min.index = index2(ibox, i, j);
298 this_min.ener = W[this_min.index];
299 if (is_local_minimum_from_below(&this_min, i, 0, index2(ibox, i-1, j ), W) &&
300 is_local_minimum_from_above(&this_min, i, ibox[0]-1, index2(ibox, i+1, j ), W) &&
301 is_local_minimum_from_below(&this_min, j, 0, index2(ibox, i, j-1), W) &&
302 is_local_minimum_from_above(&this_min, j, ibox[1]-1, index2(ibox, i, j+1), W))
304 add_minimum(fp, nmin, &this_min, mm);
311 for (i = 0; (i < ibox[0]); i++)
313 for (j = 0; (j < ibox[1]); j++)
315 for (k = 0; (k < ibox[2]); k++)
317 /* Get the index of this point in the flat array */
318 this_min.index = index3(ibox, i, j, k);
319 this_min.ener = W[this_min.index];
320 if (is_local_minimum_from_below(&this_min, i, 0, index3(ibox, i-1, j, k ), W) &&
321 is_local_minimum_from_above(&this_min, i, ibox[0]-1, index3(ibox, i+1, j, k ), W) &&
322 is_local_minimum_from_below(&this_min, j, 0, index3(ibox, i, j-1, k ), W) &&
323 is_local_minimum_from_above(&this_min, j, ibox[1]-1, index3(ibox, i, j+1, k ), W) &&
324 is_local_minimum_from_below(&this_min, k, 0, index3(ibox, i, j, k-1), W) &&
325 is_local_minimum_from_above(&this_min, k, ibox[2]-1, index3(ibox, i, j, k+1), W))
327 add_minimum(fp, nmin, &this_min, mm);
335 /* Note this treats ndim == 1 and ndim > 3 */
337 /* Set up an ndim-dimensional vector to loop over the points
338 * on the grid. (0,0,0, ... 0) is an acceptable place to
340 snew(this_point, ndim);
342 /* Determine the number of points of the ndim-dimensional
345 for (i = 1; i < ndim; i++)
351 while (loopmax > loopcounter)
353 gmx_bool bMin = TRUE;
355 /* Get the index of this_point in the flat array */
356 this_min.index = indexn(ndim, ibox, this_point);
357 this_min.ener = W[this_min.index];
359 /* Is this_point a minimum from above and below in each
361 for (i = 0; bMin && (i < ndim); i++)
363 /* Save the index of this_point within the curent
364 * dimension so we can change that index in the
365 * this_point array for use with indexn(). */
366 int index = this_point[i];
369 is_local_minimum_from_below(&this_min, index, 0, indexn(ndim, ibox, this_point), W);
372 is_local_minimum_from_above(&this_min, index, ibox[i]-1, indexn(ndim, ibox, this_point), W);
377 add_minimum(fp, nmin, &this_min, mm);
381 /* update global loop counter */
384 /* Avoid underflow of this_point[i] */
385 if (loopmax > loopcounter)
387 /* update this_point non-recursively */
390 while (ibox[i] == this_point[i])
394 /* this_point[i] cannot underflow because
395 * loopmax > loopcounter. */
404 qsort(mm, nmin, sizeof(mm[0]), comp_minima);
405 fprintf(fp, "Minima sorted after energy\n");
406 for (i = 0; (i < nmin); i++)
408 print_minimum(fp, i, &mm[i]);
414 static void do_sham(const char *fn, const char *ndx,
415 const char *xpmP, const char *xpm, const char *xpm2,
416 const char *xpm3, const char *pdb,
418 int n, int neig, real **eig,
419 gmx_bool bGE, int nenerT, real **enerT,
421 real pmax, real gmax,
422 real *emin, real *emax, int nlevels, real pmin,
423 int *idim, int *ibox,
424 gmx_bool bXmin, real *xmin, gmx_bool bXmax, real *xmax)
427 real *min_eig, *max_eig;
428 real *axis_x, *axis_y, *axis_z, *axis = NULL;
430 real **PP, *W, *E, **WW, **EE, *S, **SS, *M, *bE;
433 double *bfac, efac, bref, Pmax, Wmin, Wmax, Winf, Emin, Emax, Einf, Smin, Smax, Sinf, Mmin, Mmax;
435 int i, j, k, imin, len, index, d, *nbin, *bindex, bi;
440 t_rgb rlo = { 0, 0, 0 };
441 t_rgb rhi = { 1, 1, 1 };
443 /* Determine extremes for the eigenvectors */
450 for (i = 0; (i < neig); i++)
452 /* Check for input constraints */
453 min_eig[i] = max_eig[i] = eig[i][0];
454 for (j = 0; (j < n); j++)
456 min_eig[i] = min(min_eig[i], eig[i][j]);
457 max_eig[i] = max(max_eig[i], eig[i][j]);
458 delta[i] = (max_eig[i]-min_eig[i])/(2.0*ibox[i]);
460 /* Add some extra space, half a bin on each side, unless the
461 * user has set the limits.
465 if (max_eig[i] > xmax[i])
467 gmx_warning("Your xmax[%d] value %f is smaller than the largest data point %f", i, xmax[i], max_eig[i]);
469 max_eig[i] = xmax[i];
473 max_eig[i] += delta[i];
478 if (min_eig[i] < xmin[i])
480 gmx_warning("Your xmin[%d] value %f is larger than the smallest data point %f", i, xmin[i], min_eig[i]);
482 min_eig[i] = xmin[i];
486 min_eig[i] -= delta[i];
488 bfac[i] = ibox[i]/(max_eig[i]-min_eig[i]);
491 bref = 1/(BOLTZ*Tref);
493 if (bGE || nenerT == 2)
496 for (j = 0; (j < n); j++)
500 bE[j] = bref*enerT[0][j];
504 bE[j] = (bref - 1/(BOLTZ*enerT[1][j]))*enerT[0][j];
506 Emin = min(Emin, bE[j]);
514 for (i = 0; (i < neig); i++)
518 printf("There are %d bins in the %d-dimensional histogram. Beta-Emin = %g\n",
529 /* Loop over projections */
530 for (j = 0; (j < n); j++)
532 /* Loop over dimensions */
534 for (i = 0; (i < neig); i++)
536 nxyz[i] = bfac[i]*(eig[i][j]-min_eig[i]);
537 if (nxyz[i] < 0 || nxyz[i] >= ibox[i])
544 index = indexn(neig, ibox, nxyz);
545 range_check(index, 0, len);
546 /* Compute the exponential factor */
549 efac = exp(-bE[j]+Emin);
555 /* Apply the bin volume correction for a multi-dimensional distance */
556 for (i = 0; i < neig; i++)
562 else if (idim[i] == 3)
564 efac /= sqr(eig[i][j]);
566 else if (idim[i] == -1)
568 efac /= sin(DEG2RAD*eig[i][j]);
571 /* Update the probability */
573 /* Update the energy */
576 E[index] += enerT[0][j];
578 /* Statistics: which "structure" in which bin */
583 /* Normalize probability */
584 normalize_p_e(len, P, nbin, E, pmin);
586 /* Compute boundaries for the Free energy */
590 /* Recompute Emin: it may have changed due to averaging */
593 for (i = 0; (i < len); i++)
597 Pmax = max(P[i], Pmax);
598 W[i] = -BOLTZ*Tref*log(P[i]);
604 Emin = min(E[i], Emin);
605 Emax = max(E[i], Emax);
606 Wmax = max(W[i], Wmax);
626 /* Write out the free energy as a function of bin index */
627 fp = gmx_ffopen(fn, "w");
628 for (i = 0; (i < len); i++)
633 S[i] = E[i]-W[i]-Smin;
634 fprintf(fp, "%5d %10.5e %10.5e %10.5e\n", i, W[i], E[i], S[i]);
644 /* Organize the structures in the bins */
646 snew(b->index, len+1);
649 for (i = 0; (i < len); i++)
651 b->index[i+1] = b->index[i]+nbin[i];
654 for (i = 0; (i < n); i++)
657 b->a[b->index[bi]+nbin[bi]] = i;
660 /* Consistency check */
661 /* This no longer applies when we allow the plot to be smaller
662 than the sampled space.
663 for(i=0; (i<len); i++) {
664 if (nbin[i] != (b->index[i+1] - b->index[i]))
665 gmx_fatal(FARGS,"nbin[%d] = %d, should be %d",i,nbin[i],
666 b->index[i+1] - b->index[i]);
669 /* Write the index file */
670 fp = gmx_ffopen(ndx, "w");
671 for (i = 0; (i < len); i++)
675 fprintf(fp, "[ %d ]\n", i);
676 for (j = b->index[i]; (j < b->index[i+1]); j++)
678 fprintf(fp, "%d\n", b->a[j]+1);
683 snew(axis_x, ibox[0]+1);
684 snew(axis_y, ibox[1]+1);
685 snew(axis_z, ibox[2]+1);
686 maxbox = max(ibox[0], max(ibox[1], ibox[2]));
687 snew(PP, maxbox*maxbox);
688 snew(WW, maxbox*maxbox);
689 snew(EE, maxbox*maxbox);
690 snew(SS, maxbox*maxbox);
691 for (i = 0; (i < min(neig, 3)); i++)
695 case 0: axis = axis_x; break;
696 case 1: axis = axis_y; break;
697 case 2: axis = axis_z; break;
700 for (j = 0; j <= ibox[i]; j++)
702 axis[j] = min_eig[i] + j/bfac[i];
706 pick_minima(logf, ibox, neig, len, W);
711 flags = MAT_SPATIAL_X | MAT_SPATIAL_Y;
714 /* Dump to XPM file */
716 for (i = 0; (i < ibox[0]); i++)
718 snew(PP[i], ibox[1]);
719 for (j = 0; j < ibox[1]; j++)
721 PP[i][j] = P[i*ibox[1]+j];
723 WW[i] = &(W[i*ibox[1]]);
724 EE[i] = &(E[i*ibox[1]]);
725 SS[i] = &(S[i*ibox[1]]);
727 fp = gmx_ffopen(xpmP, "w");
728 write_xpm(fp, flags, "Probability Distribution", "", "PC1", "PC2",
729 ibox[0], ibox[1], axis_x, axis_y, PP, 0, Pmax, rlo, rhi, &nlevels);
731 fp = gmx_ffopen(xpm, "w");
732 write_xpm(fp, flags, "Gibbs Energy Landscape", "G (kJ/mol)", "PC1", "PC2",
733 ibox[0], ibox[1], axis_x, axis_y, WW, 0, gmax, rlo, rhi, &nlevels);
735 fp = gmx_ffopen(xpm2, "w");
736 write_xpm(fp, flags, "Enthalpy Landscape", "H (kJ/mol)", "PC1", "PC2",
737 ibox[0], ibox[1], axis_x, axis_y, EE,
738 emin ? *emin : Emin, emax ? *emax : Einf, rlo, rhi, &nlevels);
740 fp = gmx_ffopen(xpm3, "w");
741 write_xpm(fp, flags, "Entropy Landscape", "TDS (kJ/mol)", "PC1", "PC2",
742 ibox[0], ibox[1], axis_x, axis_y, SS, 0, Sinf, rlo, rhi, &nlevels);
747 /* Dump to PDB file */
748 fp = gmx_ffopen(pdb, "w");
749 for (i = 0; (i < ibox[0]); i++)
751 xxx[XX] = 3*(i+0.5-ibox[0]/2);
752 for (j = 0; (j < ibox[1]); j++)
754 xxx[YY] = 3*(j+0.5-ibox[1]/2);
755 for (k = 0; (k < ibox[2]); k++)
757 xxx[ZZ] = 3*(k+0.5-ibox[2]/2);
758 index = index3(ibox, i, j, k);
761 fprintf(fp, "%-6s%5u %-4.4s%3.3s %4d %8.3f%8.3f%8.3f%6.2f%6.2f\n",
762 "ATOM", (index+1) %10000, "H", "H", (index+1)%10000,
763 xxx[XX], xxx[YY], xxx[ZZ], 1.0, W[index]);
769 write_xplor("out.xplor", W, ibox, min_eig, max_eig);
770 nxyz[XX] = imin/(ibox[1]*ibox[2]);
771 nxyz[YY] = (imin-nxyz[XX]*ibox[1]*ibox[2])/ibox[2];
772 nxyz[ZZ] = imin % ibox[2];
773 for (i = 0; (i < ibox[0]); i++)
776 for (j = 0; (j < ibox[1]); j++)
778 WW[i][j] = W[index3(ibox, i, j, nxyz[ZZ])];
781 snew(buf, strlen(xpm)+4);
782 sprintf(buf, "%s", xpm);
783 sprintf(&buf[strlen(xpm)-4], "12.xpm");
784 fp = gmx_ffopen(buf, "w");
785 write_xpm(fp, flags, "Gibbs Energy Landscape", "W (kJ/mol)", "PC1", "PC2",
786 ibox[0], ibox[1], axis_x, axis_y, WW, 0, gmax, rlo, rhi, &nlevels);
788 for (i = 0; (i < ibox[0]); i++)
790 for (j = 0; (j < ibox[2]); j++)
792 WW[i][j] = W[index3(ibox, i, nxyz[YY], j)];
795 sprintf(&buf[strlen(xpm)-4], "13.xpm");
796 fp = gmx_ffopen(buf, "w");
797 write_xpm(fp, flags, "SHAM Energy Landscape", "kJ/mol", "PC1", "PC3",
798 ibox[0], ibox[2], axis_x, axis_z, WW, 0, gmax, rlo, rhi, &nlevels);
800 for (i = 0; (i < ibox[1]); i++)
802 for (j = 0; (j < ibox[2]); j++)
804 WW[i][j] = W[index3(ibox, nxyz[XX], i, j)];
807 sprintf(&buf[strlen(xpm)-4], "23.xpm");
808 fp = gmx_ffopen(buf, "w");
809 write_xpm(fp, flags, "SHAM Energy Landscape", "kJ/mol", "PC2", "PC3",
810 ibox[1], ibox[2], axis_y, axis_z, WW, 0, gmax, rlo, rhi, &nlevels);
816 static void ehisto(const char *fh, int n, real **enerT, const output_env_t oenv)
819 int i, j, k, nbin, blength;
821 real *T, bmin, bmax, bwidth;
829 for (j = 1; (j < n); j++)
831 for (k = 0; (k < nbin); k++)
833 if (T[k] == enerT[1][j])
842 T[nbin] = enerT[1][j];
845 bmin = min(enerT[0][j], bmin);
846 bmax = max(enerT[0][j], bmax);
849 blength = (bmax - bmin)/bwidth + 2;
851 for (i = 0; (i < nbin); i++)
853 snew(histo[i], blength);
855 for (j = 0; (j < n); j++)
857 k = (enerT[0][j]-bmin)/bwidth;
858 histo[bindex[j]][k]++;
860 fp = xvgropen(fh, "Energy distribution", "E (kJ/mol)", "", oenv);
861 for (j = 0; (j < blength); j++)
863 fprintf(fp, "%8.3f", bmin+j*bwidth);
864 for (k = 0; (k < nbin); k++)
866 fprintf(fp, " %6d", histo[k][j]);
873 int gmx_sham(int argc, char *argv[])
875 const char *desc[] = {
876 "[THISMODULE] makes multi-dimensional free-energy, enthalpy and entropy plots.",
877 "[THISMODULE] reads one or more [TT].xvg[tt] files and analyzes data sets.",
878 "The basic purpose of [THISMODULE] is to plot Gibbs free energy landscapes",
879 "(option [TT]-ls[tt])",
880 "by Bolzmann inverting multi-dimensional histograms (option [TT]-lp[tt]),",
882 "make enthalpy (option [TT]-lsh[tt]) and entropy (option [TT]-lss[tt])",
883 "plots. The histograms can be made for any quantities the user supplies.",
884 "A line in the input file may start with a time",
885 "(see option [TT]-time[tt]) and any number of [IT]y[it]-values may follow.",
886 "Multiple sets can also be",
887 "read when they are separated by & (option [TT]-n[tt]),",
888 "in this case only one [IT]y[it]-value is read from each line.",
889 "All lines starting with # and @ are skipped.",
891 "Option [TT]-ge[tt] can be used to supply a file with free energies",
892 "when the ensemble is not a Boltzmann ensemble, but needs to be biased",
893 "by this free energy. One free energy value is required for each",
894 "(multi-dimensional) data point in the [TT]-f[tt] input.",
896 "Option [TT]-ene[tt] can be used to supply a file with energies.",
897 "These energies are used as a weighting function in the single",
898 "histogram analysis method by Kumar et al. When temperatures",
899 "are supplied (as a second column in the file), an experimental",
900 "weighting scheme is applied. In addition the vales",
901 "are used for making enthalpy and entropy plots.",
903 "With option [TT]-dim[tt], dimensions can be gives for distances.",
904 "When a distance is 2- or 3-dimensional, the circumference or surface",
905 "sampled by two particles increases with increasing distance.",
906 "Depending on what one would like to show, one can choose to correct",
907 "the histogram and free-energy for this volume effect.",
908 "The probability is normalized by r and r^2 for dimensions of 2 and 3, ",
910 "A value of -1 is used to indicate an angle in degrees between two",
911 "vectors: a sin(angle) normalization will be applied.",
912 "[BB]Note[bb] that for angles between vectors the inner-product or cosine",
913 "is the natural quantity to use, as it will produce bins of the same",
916 static real tb = -1, te = -1, frac = 0.5, filtlen = 0;
917 static gmx_bool bHaveT = TRUE, bDer = FALSE, bSubAv = TRUE, bAverCorr = FALSE, bXYdy = FALSE;
918 static gmx_bool bEESEF = FALSE, bEENLC = FALSE, bEeFitAc = FALSE, bPower = FALSE;
919 static gmx_bool bShamEner = TRUE, bSham = TRUE;
920 static real Tref = 298.15, pmin = 0, ttol = 0, pmax = 0, gmax = 0, emin = 0, emax = 0;
921 static rvec nrdim = {1, 1, 1};
922 static rvec nrbox = {32, 32, 32};
923 static rvec xmin = {0, 0, 0}, xmax = {1, 1, 1};
924 static int nsets_in = 1, nb_min = 4, resol = 10, nlevels = 25;
926 { "-time", FALSE, etBOOL, {&bHaveT},
927 "Expect a time in the input" },
928 { "-b", FALSE, etREAL, {&tb},
929 "First time to read from set" },
930 { "-e", FALSE, etREAL, {&te},
931 "Last time to read from set" },
932 { "-ttol", FALSE, etREAL, {&ttol},
933 "Tolerance on time in appropriate units (usually ps)" },
934 { "-n", FALSE, etINT, {&nsets_in},
935 "Read this number of sets separated by lines containing only an ampersand" },
936 { "-d", FALSE, etBOOL, {&bDer},
937 "Use the derivative" },
938 { "-sham", FALSE, etBOOL, {&bSham},
939 "Turn off energy weighting even if energies are given" },
940 { "-tsham", FALSE, etREAL, {&Tref},
941 "Temperature for single histogram analysis" },
942 { "-pmin", FALSE, etREAL, {&pmin},
943 "Minimum probability. Anything lower than this will be set to zero" },
944 { "-dim", FALSE, etRVEC, {nrdim},
945 "Dimensions for distances, used for volume correction (max 3 values, dimensions > 3 will get the same value as the last)" },
946 { "-ngrid", FALSE, etRVEC, {nrbox},
947 "Number of bins for energy landscapes (max 3 values, dimensions > 3 will get the same value as the last)" },
948 { "-xmin", FALSE, etRVEC, {xmin},
949 "Minimum for the axes in energy landscape (see above for > 3 dimensions)" },
950 { "-xmax", FALSE, etRVEC, {xmax},
951 "Maximum for the axes in energy landscape (see above for > 3 dimensions)" },
952 { "-pmax", FALSE, etREAL, {&pmax},
953 "Maximum probability in output, default is calculate" },
954 { "-gmax", FALSE, etREAL, {&gmax},
955 "Maximum free energy in output, default is calculate" },
956 { "-emin", FALSE, etREAL, {&emin},
957 "Minimum enthalpy in output, default is calculate" },
958 { "-emax", FALSE, etREAL, {&emax},
959 "Maximum enthalpy in output, default is calculate" },
960 { "-nlevels", FALSE, etINT, {&nlevels},
961 "Number of levels for energy landscape" },
963 #define NPA asize(pa)
966 int n, e_n, nlast, s, nset, e_nset, d_nset, i, j = 0, *idim, *ibox;
967 real **val, **et_val, *t, *e_t, e_dt, d_dt, dt, tot, error;
969 double *av, *sig, cum1, cum2, cum3, cum4, db;
970 const char *fn_ge, *fn_ene;
972 gmx_int64_t num_grid_points;
975 { efXVG, "-f", "graph", ffREAD },
976 { efXVG, "-ge", "gibbs", ffOPTRD },
977 { efXVG, "-ene", "esham", ffOPTRD },
978 { efXVG, "-dist", "ener", ffOPTWR },
979 { efXVG, "-histo", "edist", ffOPTWR },
980 { efNDX, "-bin", "bindex", ffOPTWR },
981 { efXPM, "-lp", "prob", ffOPTWR },
982 { efXPM, "-ls", "gibbs", ffOPTWR },
983 { efXPM, "-lsh", "enthalpy", ffOPTWR },
984 { efXPM, "-lss", "entropy", ffOPTWR },
985 { efPDB, "-ls3", "gibbs3", ffOPTWR },
986 { efLOG, "-g", "shamlog", ffOPTWR }
988 #define NFILE asize(fnm)
993 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW,
994 NFILE, fnm, npargs, pa, asize(desc), desc, 0, NULL, &oenv))
999 val = read_xvg_time(opt2fn("-f", NFILE, fnm), bHaveT,
1000 opt2parg_bSet("-b", npargs, pa), tb-ttol,
1001 opt2parg_bSet("-e", npargs, pa), te+ttol,
1002 nsets_in, &nset, &n, &dt, &t);
1003 printf("Read %d sets of %d points, dt = %g\n\n", nset, n, dt);
1005 fn_ge = opt2fn_null("-ge", NFILE, fnm);
1006 fn_ene = opt2fn_null("-ene", NFILE, fnm);
1008 if (fn_ge && fn_ene)
1010 gmx_fatal(FARGS, "Can not do free energy and energy corrections at the same time");
1013 if (fn_ge || fn_ene)
1015 et_val = read_xvg_time(fn_ge ? fn_ge : fn_ene, bHaveT,
1016 opt2parg_bSet("-b", npargs, pa), tb-ttol,
1017 opt2parg_bSet("-e", npargs, pa), te+ttol,
1018 1, &e_nset, &e_n, &e_dt, &e_t);
1023 gmx_fatal(FARGS, "Can only handle one free energy component in %s",
1029 if (e_nset != 1 && e_nset != 2)
1031 gmx_fatal(FARGS, "Can only handle one energy component or one energy and one T in %s",
1037 gmx_fatal(FARGS, "Number of energies (%d) does not match number of entries (%d) in %s", e_n, n, opt2fn("-f", NFILE, fnm));
1045 if (fn_ene && et_val)
1047 ehisto(opt2fn("-histo", NFILE, fnm), e_n, et_val, oenv);
1050 snew(idim, max(3, nset));
1051 snew(ibox, max(3, nset));
1052 snew(rmin, max(3, nset));
1053 snew(rmax, max(3, nset));
1054 for (i = 0; (i < min(3, nset)); i++)
1061 for (; (i < nset); i++)
1069 /* Check that the grid size is manageable. */
1070 num_grid_points = ibox[0];
1071 for (i = 1; i < nset; i++)
1074 if (!check_int_multiply_for_overflow(num_grid_points, ibox[i], &result))
1077 "The number of dimensions and grid points is too large for this tool.\n");
1079 num_grid_points = result;
1081 /* The number of grid points fits in a gmx_int64_t. */
1083 do_sham(opt2fn("-dist", NFILE, fnm), opt2fn("-bin", NFILE, fnm),
1084 opt2fn("-lp", NFILE, fnm),
1085 opt2fn("-ls", NFILE, fnm), opt2fn("-lsh", NFILE, fnm),
1086 opt2fn("-lss", NFILE, fnm),
1087 opt2fn("-ls3", NFILE, fnm), opt2fn("-g", NFILE, fnm),
1088 n, nset, val, fn_ge != NULL, e_nset, et_val, Tref,
1090 opt2parg_bSet("-emin", NPA, pa) ? &emin : NULL,
1091 opt2parg_bSet("-emax", NPA, pa) ? &emax : NULL,
1094 opt2parg_bSet("-xmin", NPA, pa), rmin,
1095 opt2parg_bSet("-xmax", NPA, pa), rmax);