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.
43 #include "gromacs/commandline/pargs.h"
45 #include "gromacs/utility/smalloc.h"
47 #include "gromacs/utility/fatalerror.h"
48 #include "gromacs/math/vec.h"
49 #include "gromacs/utility/futil.h"
53 #include "gromacs/fileio/xvgr.h"
54 #include "gromacs/math/units.h"
55 #include "gromacs/fileio/pdbio.h"
56 #include "gromacs/fileio/matio.h"
60 static int index2(int *ibox, int x, int y)
65 static int index3(int *ibox, int x, int y, int z)
67 return (ibox[2]*(ibox[1]*x+y)+z);
70 static gmx_int64_t indexn(int ndim, const int *ibox, const int *nxyz)
75 /* Compute index in 1-D array */
77 for (k = 0; (k < ndim); k++)
80 for (kk = k+1; (kk < ndim); kk++)
90 int Nx; /* x grid points in unit cell */
91 int Ny; /* y grid points in unit cell */
92 int Nz; /* z grid points in unit cell */
93 int dmin[3]; /* starting point x,y,z*/
94 int dmax[3]; /* ending point x,y,z */
95 real cell[6]; /* usual cell parameters */
99 static void lo_write_xplor(XplorMap * map, const char * file)
104 fp = gmx_ffopen(file, "w");
105 /* The REMARKS part is the worst part of the XPLOR format
106 * and may cause problems with some programs
108 fprintf(fp, "\n 2 !NTITLE\n");
109 fprintf(fp, " REMARKS Energy Landscape from GROMACS\n");
110 fprintf(fp, " REMARKS DATE: 2004-12-21 \n");
111 fprintf(fp, " %7d %7d %7d %7d %7d %7d %7d %7d %7d\n",
112 map->Nx, map->dmin[0], map->dmax[0],
113 map->Ny, map->dmin[1], map->dmax[1],
114 map->Nz, map->dmin[2], map->dmax[2]);
115 fprintf(fp, "%12.5E%12.5E%12.5E%12.5E%12.5E%12.5E\n",
116 map->cell[0], map->cell[1], map->cell[2],
117 map->cell[3], map->cell[4], map->cell[5]);
118 fprintf(fp, "ZYX\n");
121 for (n = 0; n < map->Nz; n++, z++)
123 fprintf(fp, "%8d\n", z);
124 for (i = 0; i < map->Nx*map->Ny; i += 6)
126 for (j = 0; j < 6; j++)
128 if (i+j < map->Nx*map->Ny)
130 fprintf(fp, "%12.5E", map->ed[n*map->Nx*map->Ny+i+j]);
136 fprintf(fp, " -9999\n");
140 static void write_xplor(const char *file, real *data, int *ibox, real dmin[], real dmax[])
149 snew(xm->ed, xm->Nx*xm->Ny*xm->Nz);
151 for (k = 0; (k < xm->Nz); k++)
153 for (j = 0; (j < xm->Ny); j++)
155 for (i = 0; (i < xm->Nx); i++)
157 xm->ed[n++] = data[index3(ibox, i, j, k)];
161 xm->cell[0] = dmax[XX]-dmin[XX];
162 xm->cell[1] = dmax[YY]-dmin[YY];
163 xm->cell[2] = dmax[ZZ]-dmin[ZZ];
164 xm->cell[3] = xm->cell[4] = xm->cell[5] = 90;
166 clear_ivec(xm->dmin);
167 xm->dmax[XX] = ibox[XX]-1;
168 xm->dmax[YY] = ibox[YY]-1;
169 xm->dmax[ZZ] = ibox[ZZ]-1;
171 lo_write_xplor(xm, file);
177 static void normalize_p_e(int len, double *P, int *nbin, real *E, real pmin)
182 for (i = 0; (i < len); i++)
190 printf("Ptot = %g\n", Ptot);
191 for (i = 0; (i < len); i++)
194 /* Have to check for pmin after normalizing to prevent "stretching"
209 static int comp_minima(const void *a, const void *b)
211 t_minimum *ma = (t_minimum *) a;
212 t_minimum *mb = (t_minimum *) b;
214 if (ma->ener < mb->ener)
218 else if (ma->ener > mb->ener)
229 void print_minimum(FILE *fp, int num, const t_minimum *min)
232 "Minimum %d at index " "%"GMX_PRId64 " energy %10.3f\n",
233 num, min->index, min->ener);
237 void add_minimum(FILE *fp, int num, const t_minimum *min, t_minimum *mm)
239 print_minimum(fp, num, min);
240 mm[num].index = min->index;
241 mm[num].ener = min->ener;
245 gmx_bool is_local_minimum_from_below(const t_minimum *this_min,
251 return ((dimension_index == dimension_min) ||
252 ((dimension_index > dimension_min) &&
253 (this_min->ener < W[neighbour_index])));
254 /* Note over/underflow within W cannot occur. */
258 gmx_bool is_local_minimum_from_above(const t_minimum *this_min,
264 return ((dimension_index == dimension_max) ||
265 ((dimension_index < dimension_max) &&
266 (this_min->ener < W[neighbour_index])));
267 /* Note over/underflow within W cannot occur. */
270 static void pick_minima(const char *logfile, int *ibox, int ndim, int len, real W[])
274 t_minimum *mm, this_min;
276 int loopmax, loopcounter;
280 fp = gmx_ffopen(logfile, "w");
281 /* Loop over each element in the array of dimenion ndim seeking
282 * minima with respect to every dimension. Specialized loops for
283 * speed with ndim == 2 and ndim == 3. */
287 /* This is probably impossible to reach anyway. */
290 for (i = 0; (i < ibox[0]); i++)
292 for (j = 0; (j < ibox[1]); j++)
294 /* Get the index of this point in the flat array */
295 this_min.index = index2(ibox, i, j);
296 this_min.ener = W[this_min.index];
297 if (is_local_minimum_from_below(&this_min, i, 0, index2(ibox, i-1, j ), W) &&
298 is_local_minimum_from_above(&this_min, i, ibox[0]-1, index2(ibox, i+1, j ), W) &&
299 is_local_minimum_from_below(&this_min, j, 0, index2(ibox, i, j-1), W) &&
300 is_local_minimum_from_above(&this_min, j, ibox[1]-1, index2(ibox, i, j+1), W))
302 add_minimum(fp, nmin, &this_min, mm);
309 for (i = 0; (i < ibox[0]); i++)
311 for (j = 0; (j < ibox[1]); j++)
313 for (k = 0; (k < ibox[2]); k++)
315 /* Get the index of this point in the flat array */
316 this_min.index = index3(ibox, i, j, k);
317 this_min.ener = W[this_min.index];
318 if (is_local_minimum_from_below(&this_min, i, 0, index3(ibox, i-1, j, k ), W) &&
319 is_local_minimum_from_above(&this_min, i, ibox[0]-1, index3(ibox, i+1, j, k ), W) &&
320 is_local_minimum_from_below(&this_min, j, 0, index3(ibox, i, j-1, k ), W) &&
321 is_local_minimum_from_above(&this_min, j, ibox[1]-1, index3(ibox, i, j+1, k ), W) &&
322 is_local_minimum_from_below(&this_min, k, 0, index3(ibox, i, j, k-1), W) &&
323 is_local_minimum_from_above(&this_min, k, ibox[2]-1, index3(ibox, i, j, k+1), W))
325 add_minimum(fp, nmin, &this_min, mm);
333 /* Note this treats ndim == 1 and ndim > 3 */
335 /* Set up an ndim-dimensional vector to loop over the points
336 * on the grid. (0,0,0, ... 0) is an acceptable place to
338 snew(this_point, ndim);
340 /* Determine the number of points of the ndim-dimensional
343 for (i = 1; i < ndim; i++)
349 while (loopmax > loopcounter)
351 gmx_bool bMin = TRUE;
353 /* Get the index of this_point in the flat array */
354 this_min.index = indexn(ndim, ibox, this_point);
355 this_min.ener = W[this_min.index];
357 /* Is this_point a minimum from above and below in each
359 for (i = 0; bMin && (i < ndim); i++)
361 /* Save the index of this_point within the curent
362 * dimension so we can change that index in the
363 * this_point array for use with indexn(). */
364 int index = this_point[i];
367 is_local_minimum_from_below(&this_min, index, 0, indexn(ndim, ibox, this_point), W);
370 is_local_minimum_from_above(&this_min, index, ibox[i]-1, indexn(ndim, ibox, this_point), W);
375 add_minimum(fp, nmin, &this_min, mm);
379 /* update global loop counter */
382 /* Avoid underflow of this_point[i] */
383 if (loopmax > loopcounter)
385 /* update this_point non-recursively */
388 while (ibox[i] == this_point[i])
392 /* this_point[i] cannot underflow because
393 * loopmax > loopcounter. */
402 qsort(mm, nmin, sizeof(mm[0]), comp_minima);
403 fprintf(fp, "Minima sorted after energy\n");
404 for (i = 0; (i < nmin); i++)
406 print_minimum(fp, i, &mm[i]);
412 static void do_sham(const char *fn, const char *ndx,
413 const char *xpmP, const char *xpm, const char *xpm2,
414 const char *xpm3, const char *pdb,
416 int n, int neig, real **eig,
417 gmx_bool bGE, int nenerT, real **enerT,
419 real pmax, real gmax,
420 real *emin, real *emax, int nlevels, real pmin,
421 int *idim, int *ibox,
422 gmx_bool bXmin, real *xmin, gmx_bool bXmax, real *xmax)
425 real *min_eig, *max_eig;
426 real *axis_x, *axis_y, *axis_z, *axis = NULL;
428 real **PP, *W, *E, **WW, **EE, *S, **SS, *M, *bE;
431 double *bfac, efac, bref, Pmax, Wmin, Wmax, Winf, Emin, Emax, Einf, Smin, Smax, Sinf, Mmin, Mmax;
433 int i, j, k, imin, len, index, d, *nbin, *bindex, bi;
438 t_rgb rlo = { 0, 0, 0 };
439 t_rgb rhi = { 1, 1, 1 };
441 /* Determine extremes for the eigenvectors */
448 for (i = 0; (i < neig); i++)
450 /* Check for input constraints */
451 min_eig[i] = max_eig[i] = eig[i][0];
452 for (j = 0; (j < n); j++)
454 min_eig[i] = min(min_eig[i], eig[i][j]);
455 max_eig[i] = max(max_eig[i], eig[i][j]);
456 delta[i] = (max_eig[i]-min_eig[i])/(2.0*ibox[i]);
458 /* Add some extra space, half a bin on each side, unless the
459 * user has set the limits.
463 if (max_eig[i] > xmax[i])
465 gmx_warning("Your xmax[%d] value %f is smaller than the largest data point %f", i, xmax[i], max_eig[i]);
467 max_eig[i] = xmax[i];
471 max_eig[i] += delta[i];
476 if (min_eig[i] < xmin[i])
478 gmx_warning("Your xmin[%d] value %f is larger than the smallest data point %f", i, xmin[i], min_eig[i]);
480 min_eig[i] = xmin[i];
484 min_eig[i] -= delta[i];
486 bfac[i] = ibox[i]/(max_eig[i]-min_eig[i]);
489 bref = 1/(BOLTZ*Tref);
491 if (bGE || nenerT == 2)
494 for (j = 0; (j < n); j++)
498 bE[j] = bref*enerT[0][j];
502 bE[j] = (bref - 1/(BOLTZ*enerT[1][j]))*enerT[0][j];
504 Emin = min(Emin, bE[j]);
512 for (i = 0; (i < neig); i++)
516 printf("There are %d bins in the %d-dimensional histogram. Beta-Emin = %g\n",
527 /* Loop over projections */
528 for (j = 0; (j < n); j++)
530 /* Loop over dimensions */
532 for (i = 0; (i < neig); i++)
534 nxyz[i] = bfac[i]*(eig[i][j]-min_eig[i]);
535 if (nxyz[i] < 0 || nxyz[i] >= ibox[i])
542 index = indexn(neig, ibox, nxyz);
543 range_check(index, 0, len);
544 /* Compute the exponential factor */
547 efac = exp(-bE[j]+Emin);
553 /* Apply the bin volume correction for a multi-dimensional distance */
554 for (i = 0; i < neig; i++)
560 else if (idim[i] == 3)
562 efac /= sqr(eig[i][j]);
564 else if (idim[i] == -1)
566 efac /= sin(DEG2RAD*eig[i][j]);
569 /* Update the probability */
571 /* Update the energy */
574 E[index] += enerT[0][j];
576 /* Statistics: which "structure" in which bin */
581 /* Normalize probability */
582 normalize_p_e(len, P, nbin, E, pmin);
584 /* Compute boundaries for the Free energy */
588 /* Recompute Emin: it may have changed due to averaging */
591 for (i = 0; (i < len); i++)
595 Pmax = max(P[i], Pmax);
596 W[i] = -BOLTZ*Tref*log(P[i]);
602 Emin = min(E[i], Emin);
603 Emax = max(E[i], Emax);
604 Wmax = max(W[i], Wmax);
624 /* Write out the free energy as a function of bin index */
625 fp = gmx_ffopen(fn, "w");
626 for (i = 0; (i < len); i++)
631 S[i] = E[i]-W[i]-Smin;
632 fprintf(fp, "%5d %10.5e %10.5e %10.5e\n", i, W[i], E[i], S[i]);
642 /* Organize the structures in the bins */
644 snew(b->index, len+1);
647 for (i = 0; (i < len); i++)
649 b->index[i+1] = b->index[i]+nbin[i];
652 for (i = 0; (i < n); i++)
655 b->a[b->index[bi]+nbin[bi]] = i;
658 /* Consistency check */
659 /* This no longer applies when we allow the plot to be smaller
660 than the sampled space.
661 for(i=0; (i<len); i++) {
662 if (nbin[i] != (b->index[i+1] - b->index[i]))
663 gmx_fatal(FARGS,"nbin[%d] = %d, should be %d",i,nbin[i],
664 b->index[i+1] - b->index[i]);
667 /* Write the index file */
668 fp = gmx_ffopen(ndx, "w");
669 for (i = 0; (i < len); i++)
673 fprintf(fp, "[ %d ]\n", i);
674 for (j = b->index[i]; (j < b->index[i+1]); j++)
676 fprintf(fp, "%d\n", b->a[j]+1);
681 snew(axis_x, ibox[0]+1);
682 snew(axis_y, ibox[1]+1);
683 snew(axis_z, ibox[2]+1);
684 maxbox = max(ibox[0], max(ibox[1], ibox[2]));
685 snew(PP, maxbox*maxbox);
686 snew(WW, maxbox*maxbox);
687 snew(EE, maxbox*maxbox);
688 snew(SS, maxbox*maxbox);
689 for (i = 0; (i < min(neig, 3)); i++)
693 case 0: axis = axis_x; break;
694 case 1: axis = axis_y; break;
695 case 2: axis = axis_z; break;
698 for (j = 0; j <= ibox[i]; j++)
700 axis[j] = min_eig[i] + j/bfac[i];
704 pick_minima(logf, ibox, neig, len, W);
709 flags = MAT_SPATIAL_X | MAT_SPATIAL_Y;
712 /* Dump to XPM file */
714 for (i = 0; (i < ibox[0]); i++)
716 snew(PP[i], ibox[1]);
717 for (j = 0; j < ibox[1]; j++)
719 PP[i][j] = P[i*ibox[1]+j];
721 WW[i] = &(W[i*ibox[1]]);
722 EE[i] = &(E[i*ibox[1]]);
723 SS[i] = &(S[i*ibox[1]]);
725 fp = gmx_ffopen(xpmP, "w");
726 write_xpm(fp, flags, "Probability Distribution", "", "PC1", "PC2",
727 ibox[0], ibox[1], axis_x, axis_y, PP, 0, Pmax, rlo, rhi, &nlevels);
729 fp = gmx_ffopen(xpm, "w");
730 write_xpm(fp, flags, "Gibbs Energy Landscape", "G (kJ/mol)", "PC1", "PC2",
731 ibox[0], ibox[1], axis_x, axis_y, WW, 0, gmax, rlo, rhi, &nlevels);
733 fp = gmx_ffopen(xpm2, "w");
734 write_xpm(fp, flags, "Enthalpy Landscape", "H (kJ/mol)", "PC1", "PC2",
735 ibox[0], ibox[1], axis_x, axis_y, EE,
736 emin ? *emin : Emin, emax ? *emax : Einf, rlo, rhi, &nlevels);
738 fp = gmx_ffopen(xpm3, "w");
739 write_xpm(fp, flags, "Entropy Landscape", "TDS (kJ/mol)", "PC1", "PC2",
740 ibox[0], ibox[1], axis_x, axis_y, SS, 0, Sinf, rlo, rhi, &nlevels);
745 /* Dump to PDB file */
746 fp = gmx_ffopen(pdb, "w");
747 for (i = 0; (i < ibox[0]); i++)
749 xxx[XX] = 3*(i+0.5-ibox[0]/2);
750 for (j = 0; (j < ibox[1]); j++)
752 xxx[YY] = 3*(j+0.5-ibox[1]/2);
753 for (k = 0; (k < ibox[2]); k++)
755 xxx[ZZ] = 3*(k+0.5-ibox[2]/2);
756 index = index3(ibox, i, j, k);
759 fprintf(fp, "%-6s%5u %-4.4s%3.3s %4d %8.3f%8.3f%8.3f%6.2f%6.2f\n",
760 "ATOM", (index+1) %10000, "H", "H", (index+1)%10000,
761 xxx[XX], xxx[YY], xxx[ZZ], 1.0, W[index]);
767 write_xplor("out.xplor", W, ibox, min_eig, max_eig);
768 nxyz[XX] = imin/(ibox[1]*ibox[2]);
769 nxyz[YY] = (imin-nxyz[XX]*ibox[1]*ibox[2])/ibox[2];
770 nxyz[ZZ] = imin % ibox[2];
771 for (i = 0; (i < ibox[0]); i++)
774 for (j = 0; (j < ibox[1]); j++)
776 WW[i][j] = W[index3(ibox, i, j, nxyz[ZZ])];
779 snew(buf, strlen(xpm)+4);
780 sprintf(buf, "%s", xpm);
781 sprintf(&buf[strlen(xpm)-4], "12.xpm");
782 fp = gmx_ffopen(buf, "w");
783 write_xpm(fp, flags, "Gibbs Energy Landscape", "W (kJ/mol)", "PC1", "PC2",
784 ibox[0], ibox[1], axis_x, axis_y, WW, 0, gmax, rlo, rhi, &nlevels);
786 for (i = 0; (i < ibox[0]); i++)
788 for (j = 0; (j < ibox[2]); j++)
790 WW[i][j] = W[index3(ibox, i, nxyz[YY], j)];
793 sprintf(&buf[strlen(xpm)-4], "13.xpm");
794 fp = gmx_ffopen(buf, "w");
795 write_xpm(fp, flags, "SHAM Energy Landscape", "kJ/mol", "PC1", "PC3",
796 ibox[0], ibox[2], axis_x, axis_z, WW, 0, gmax, rlo, rhi, &nlevels);
798 for (i = 0; (i < ibox[1]); i++)
800 for (j = 0; (j < ibox[2]); j++)
802 WW[i][j] = W[index3(ibox, nxyz[XX], i, j)];
805 sprintf(&buf[strlen(xpm)-4], "23.xpm");
806 fp = gmx_ffopen(buf, "w");
807 write_xpm(fp, flags, "SHAM Energy Landscape", "kJ/mol", "PC2", "PC3",
808 ibox[1], ibox[2], axis_y, axis_z, WW, 0, gmax, rlo, rhi, &nlevels);
814 static void ehisto(const char *fh, int n, real **enerT, const output_env_t oenv)
817 int i, j, k, nbin, blength;
819 real *T, bmin, bmax, bwidth;
827 for (j = 1; (j < n); j++)
829 for (k = 0; (k < nbin); k++)
831 if (T[k] == enerT[1][j])
840 T[nbin] = enerT[1][j];
843 bmin = min(enerT[0][j], bmin);
844 bmax = max(enerT[0][j], bmax);
847 blength = (bmax - bmin)/bwidth + 2;
849 for (i = 0; (i < nbin); i++)
851 snew(histo[i], blength);
853 for (j = 0; (j < n); j++)
855 k = (enerT[0][j]-bmin)/bwidth;
856 histo[bindex[j]][k]++;
858 fp = xvgropen(fh, "Energy distribution", "E (kJ/mol)", "", oenv);
859 for (j = 0; (j < blength); j++)
861 fprintf(fp, "%8.3f", bmin+j*bwidth);
862 for (k = 0; (k < nbin); k++)
864 fprintf(fp, " %6d", histo[k][j]);
871 int gmx_sham(int argc, char *argv[])
873 const char *desc[] = {
874 "[THISMODULE] makes multi-dimensional free-energy, enthalpy and entropy plots.",
875 "[THISMODULE] reads one or more [TT].xvg[tt] files and analyzes data sets.",
876 "The basic purpose of [THISMODULE] is to plot Gibbs free energy landscapes",
877 "(option [TT]-ls[tt])",
878 "by Bolzmann inverting multi-dimensional histograms (option [TT]-lp[tt]),",
880 "make enthalpy (option [TT]-lsh[tt]) and entropy (option [TT]-lss[tt])",
881 "plots. The histograms can be made for any quantities the user supplies.",
882 "A line in the input file may start with a time",
883 "(see option [TT]-time[tt]) and any number of [IT]y[it]-values may follow.",
884 "Multiple sets can also be",
885 "read when they are separated by & (option [TT]-n[tt]),",
886 "in this case only one [IT]y[it]-value is read from each line.",
887 "All lines starting with # and @ are skipped.",
889 "Option [TT]-ge[tt] can be used to supply a file with free energies",
890 "when the ensemble is not a Boltzmann ensemble, but needs to be biased",
891 "by this free energy. One free energy value is required for each",
892 "(multi-dimensional) data point in the [TT]-f[tt] input.",
894 "Option [TT]-ene[tt] can be used to supply a file with energies.",
895 "These energies are used as a weighting function in the single",
896 "histogram analysis method by Kumar et al. When temperatures",
897 "are supplied (as a second column in the file), an experimental",
898 "weighting scheme is applied. In addition the vales",
899 "are used for making enthalpy and entropy plots.",
901 "With option [TT]-dim[tt], dimensions can be gives for distances.",
902 "When a distance is 2- or 3-dimensional, the circumference or surface",
903 "sampled by two particles increases with increasing distance.",
904 "Depending on what one would like to show, one can choose to correct",
905 "the histogram and free-energy for this volume effect.",
906 "The probability is normalized by r and r^2 for dimensions of 2 and 3, ",
908 "A value of -1 is used to indicate an angle in degrees between two",
909 "vectors: a sin(angle) normalization will be applied.",
910 "[BB]Note[bb] that for angles between vectors the inner-product or cosine",
911 "is the natural quantity to use, as it will produce bins of the same",
914 static real tb = -1, te = -1, frac = 0.5, filtlen = 0;
915 static gmx_bool bHaveT = TRUE, bDer = FALSE, bSubAv = TRUE, bAverCorr = FALSE, bXYdy = FALSE;
916 static gmx_bool bEESEF = FALSE, bEENLC = FALSE, bEeFitAc = FALSE, bPower = FALSE;
917 static gmx_bool bShamEner = TRUE, bSham = TRUE;
918 static real Tref = 298.15, pmin = 0, ttol = 0, pmax = 0, gmax = 0, emin = 0, emax = 0;
919 static rvec nrdim = {1, 1, 1};
920 static rvec nrbox = {32, 32, 32};
921 static rvec xmin = {0, 0, 0}, xmax = {1, 1, 1};
922 static int nsets_in = 1, nb_min = 4, resol = 10, nlevels = 25;
924 { "-time", FALSE, etBOOL, {&bHaveT},
925 "Expect a time in the input" },
926 { "-b", FALSE, etREAL, {&tb},
927 "First time to read from set" },
928 { "-e", FALSE, etREAL, {&te},
929 "Last time to read from set" },
930 { "-ttol", FALSE, etREAL, {&ttol},
931 "Tolerance on time in appropriate units (usually ps)" },
932 { "-n", FALSE, etINT, {&nsets_in},
933 "Read this number of sets separated by lines containing only an ampersand" },
934 { "-d", FALSE, etBOOL, {&bDer},
935 "Use the derivative" },
936 { "-sham", FALSE, etBOOL, {&bSham},
937 "Turn off energy weighting even if energies are given" },
938 { "-tsham", FALSE, etREAL, {&Tref},
939 "Temperature for single histogram analysis" },
940 { "-pmin", FALSE, etREAL, {&pmin},
941 "Minimum probability. Anything lower than this will be set to zero" },
942 { "-dim", FALSE, etRVEC, {nrdim},
943 "Dimensions for distances, used for volume correction (max 3 values, dimensions > 3 will get the same value as the last)" },
944 { "-ngrid", FALSE, etRVEC, {nrbox},
945 "Number of bins for energy landscapes (max 3 values, dimensions > 3 will get the same value as the last)" },
946 { "-xmin", FALSE, etRVEC, {xmin},
947 "Minimum for the axes in energy landscape (see above for > 3 dimensions)" },
948 { "-xmax", FALSE, etRVEC, {xmax},
949 "Maximum for the axes in energy landscape (see above for > 3 dimensions)" },
950 { "-pmax", FALSE, etREAL, {&pmax},
951 "Maximum probability in output, default is calculate" },
952 { "-gmax", FALSE, etREAL, {&gmax},
953 "Maximum free energy in output, default is calculate" },
954 { "-emin", FALSE, etREAL, {&emin},
955 "Minimum enthalpy in output, default is calculate" },
956 { "-emax", FALSE, etREAL, {&emax},
957 "Maximum enthalpy in output, default is calculate" },
958 { "-nlevels", FALSE, etINT, {&nlevels},
959 "Number of levels for energy landscape" },
961 #define NPA asize(pa)
964 int n, e_n, nlast, s, nset, e_nset, d_nset, i, j = 0, *idim, *ibox;
965 real **val, **et_val, *t, *e_t, e_dt, d_dt, dt, tot, error;
967 double *av, *sig, cum1, cum2, cum3, cum4, db;
968 const char *fn_ge, *fn_ene;
970 gmx_int64_t num_grid_points;
973 { efXVG, "-f", "graph", ffREAD },
974 { efXVG, "-ge", "gibbs", ffOPTRD },
975 { efXVG, "-ene", "esham", ffOPTRD },
976 { efXVG, "-dist", "ener", ffOPTWR },
977 { efXVG, "-histo", "edist", ffOPTWR },
978 { efNDX, "-bin", "bindex", ffOPTWR },
979 { efXPM, "-lp", "prob", ffOPTWR },
980 { efXPM, "-ls", "gibbs", ffOPTWR },
981 { efXPM, "-lsh", "enthalpy", ffOPTWR },
982 { efXPM, "-lss", "entropy", ffOPTWR },
983 { efPDB, "-ls3", "gibbs3", ffOPTWR },
984 { efLOG, "-g", "shamlog", ffOPTWR }
986 #define NFILE asize(fnm)
991 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_BE_NICE,
992 NFILE, fnm, npargs, pa, asize(desc), desc, 0, NULL, &oenv))
997 val = read_xvg_time(opt2fn("-f", NFILE, fnm), bHaveT,
998 opt2parg_bSet("-b", npargs, pa), tb-ttol,
999 opt2parg_bSet("-e", npargs, pa), te+ttol,
1000 nsets_in, &nset, &n, &dt, &t);
1001 printf("Read %d sets of %d points, dt = %g\n\n", nset, n, dt);
1003 fn_ge = opt2fn_null("-ge", NFILE, fnm);
1004 fn_ene = opt2fn_null("-ene", NFILE, fnm);
1006 if (fn_ge && fn_ene)
1008 gmx_fatal(FARGS, "Can not do free energy and energy corrections at the same time");
1011 if (fn_ge || fn_ene)
1013 et_val = read_xvg_time(fn_ge ? fn_ge : fn_ene, bHaveT,
1014 opt2parg_bSet("-b", npargs, pa), tb-ttol,
1015 opt2parg_bSet("-e", npargs, pa), te+ttol,
1016 1, &e_nset, &e_n, &e_dt, &e_t);
1021 gmx_fatal(FARGS, "Can only handle one free energy component in %s",
1027 if (e_nset != 1 && e_nset != 2)
1029 gmx_fatal(FARGS, "Can only handle one energy component or one energy and one T in %s",
1035 gmx_fatal(FARGS, "Number of energies (%d) does not match number of entries (%d) in %s", e_n, n, opt2fn("-f", NFILE, fnm));
1043 if (fn_ene && et_val)
1045 ehisto(opt2fn("-histo", NFILE, fnm), e_n, et_val, oenv);
1048 snew(idim, max(3, nset));
1049 snew(ibox, max(3, nset));
1050 snew(rmin, max(3, nset));
1051 snew(rmax, max(3, nset));
1052 for (i = 0; (i < min(3, nset)); i++)
1059 for (; (i < nset); i++)
1067 /* Check that the grid size is manageable. */
1068 num_grid_points = ibox[0];
1069 for (i = 1; i < nset; i++)
1072 if (!check_int_multiply_for_overflow(num_grid_points, ibox[i], &result))
1075 "The number of dimensions and grid points is too large for this tool.\n");
1077 num_grid_points = result;
1079 /* The number of grid points fits in a gmx_int64_t. */
1081 do_sham(opt2fn("-dist", NFILE, fnm), opt2fn("-bin", NFILE, fnm),
1082 opt2fn("-lp", NFILE, fnm),
1083 opt2fn("-ls", NFILE, fnm), opt2fn("-lsh", NFILE, fnm),
1084 opt2fn("-lss", NFILE, fnm),
1085 opt2fn("-ls3", NFILE, fnm), opt2fn("-g", NFILE, fnm),
1086 n, nset, val, fn_ge != NULL, e_nset, et_val, Tref,
1088 opt2parg_bSet("-emin", NPA, pa) ? &emin : NULL,
1089 opt2parg_bSet("-emax", NPA, pa) ? &emax : NULL,
1092 opt2parg_bSet("-xmin", NPA, pa), rmin,
1093 opt2parg_bSet("-xmax", NPA, pa), rmax);