3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
45 #include "gmx_fatal.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_large_int_t indexn(int ndim, const int *ibox, const int *nxyz)
72 gmx_large_int_t d, dd;
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 = 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"
205 gmx_large_int_t index;
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_large_int_pfmt " 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 = 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 *xpm4, const char *pdb,
416 int n, int neig, real **eig,
417 gmx_bool bGE, int nenerT, real **enerT,
418 int nmap, real *mapindex, real **map,
420 real pmax, real gmax,
421 real *emin, real *emax, int nlevels, real pmin,
422 const char *mname, gmx_bool bSham, int *idim, int *ibox,
423 gmx_bool bXmin, real *xmin, gmx_bool bXmax, real *xmax)
426 real *min_eig, *max_eig;
427 real *axis_x, *axis_y, *axis_z, *axis = NULL;
429 real **PP, *W, *E, **WW, **EE, *S, **SS, *M, **MM, *bE;
432 double *bfac, efac, bref, Pmax, Wmin, Wmax, Winf, Emin, Emax, Einf, Smin, Smax, Sinf, Mmin, Mmax, Minf;
434 int i, j, k, imin, len, index, d, *nbin, *bindex, bi;
439 t_rgb rlo = { 0, 0, 0 };
440 t_rgb rhi = { 1, 1, 1 };
442 /* Determine extremes for the eigenvectors */
449 for (i = 0; (i < neig); i++)
451 /* Check for input constraints */
452 min_eig[i] = max_eig[i] = eig[i][0];
453 for (j = 0; (j < n); j++)
455 min_eig[i] = min(min_eig[i], eig[i][j]);
456 max_eig[i] = max(max_eig[i], eig[i][j]);
457 delta[i] = (max_eig[i]-min_eig[i])/(2.0*ibox[i]);
459 /* Add some extra space, half a bin on each side, unless the
460 * user has set the limits.
464 if (max_eig[i] > xmax[i])
466 gmx_warning("Your xmax[%d] value %f is smaller than the largest data point %f", i, xmax[i], max_eig[i]);
468 max_eig[i] = xmax[i];
472 max_eig[i] += delta[i];
477 if (min_eig[i] < xmin[i])
479 gmx_warning("Your xmin[%d] value %f is larger than the smallest data point %f", i, xmin[i], min_eig[i]);
481 min_eig[i] = xmin[i];
485 min_eig[i] -= delta[i];
487 bfac[i] = ibox[i]/(max_eig[i]-min_eig[i]);
490 bref = 1/(BOLTZ*Tref);
492 if (bGE || nenerT == 2)
495 for (j = 0; (j < n); j++)
499 bE[j] = bref*enerT[0][j];
503 bE[j] = (bref - 1/(BOLTZ*enerT[1][j]))*enerT[0][j];
505 Emin = min(Emin, bE[j]);
513 for (i = 0; (i < neig); i++)
517 printf("There are %d bins in the %d-dimensional histogram. Beta-Emin = %g\n",
528 /* Loop over projections */
529 for (j = 0; (j < n); j++)
531 /* Loop over dimensions */
533 for (i = 0; (i < neig); i++)
535 nxyz[i] = bfac[i]*(eig[i][j]-min_eig[i]);
536 if (nxyz[i] < 0 || nxyz[i] >= ibox[i])
543 index = indexn(neig, ibox, nxyz);
544 range_check(index, 0, len);
545 /* Compute the exponential factor */
548 efac = exp(-bE[j]+Emin);
554 /* Apply the bin volume correction for a multi-dimensional distance */
555 for (i = 0; i < neig; i++)
561 else if (idim[i] == 3)
563 efac /= sqr(eig[i][j]);
565 else if (idim[i] == -1)
567 efac /= sin(DEG2RAD*eig[i][j]);
570 /* Update the probability */
572 /* Update the energy */
575 E[index] += enerT[0][j];
577 /* Statistics: which "structure" in which bin */
582 /* Normalize probability */
583 normalize_p_e(len, P, nbin, E, pmin);
585 /* Compute boundaries for the Free energy */
589 /* Recompute Emin: it may have changed due to averaging */
592 for (i = 0; (i < len); i++)
596 Pmax = max(P[i], Pmax);
597 W[i] = -BOLTZ*Tref*log(P[i]);
603 Emin = min(E[i], Emin);
604 Emax = max(E[i], Emax);
605 Wmax = max(W[i], Wmax);
625 /* Write out the free energy as a function of bin index */
626 fp = ffopen(fn, "w");
627 for (i = 0; (i < len); i++)
632 S[i] = E[i]-W[i]-Smin;
633 fprintf(fp, "%5d %10.5e %10.5e %10.5e\n", i, W[i], E[i], S[i]);
643 /* Organize the structures in the bins */
645 snew(b->index, len+1);
648 for (i = 0; (i < len); i++)
650 b->index[i+1] = b->index[i]+nbin[i];
653 for (i = 0; (i < n); i++)
656 b->a[b->index[bi]+nbin[bi]] = i;
659 /* Consistency check */
660 /* This no longer applies when we allow the plot to be smaller
661 than the sampled space.
662 for(i=0; (i<len); i++) {
663 if (nbin[i] != (b->index[i+1] - b->index[i]))
664 gmx_fatal(FARGS,"nbin[%d] = %d, should be %d",i,nbin[i],
665 b->index[i+1] - b->index[i]);
668 /* Write the index file */
669 fp = ffopen(ndx, "w");
670 for (i = 0; (i < len); i++)
674 fprintf(fp, "[ %d ]\n", i);
675 for (j = b->index[i]; (j < b->index[i+1]); j++)
677 fprintf(fp, "%d\n", b->a[j]+1);
682 snew(axis_x, ibox[0]+1);
683 snew(axis_y, ibox[1]+1);
684 snew(axis_z, ibox[2]+1);
685 maxbox = max(ibox[0], max(ibox[1], ibox[2]));
686 snew(PP, maxbox*maxbox);
687 snew(WW, maxbox*maxbox);
688 snew(EE, maxbox*maxbox);
689 snew(SS, maxbox*maxbox);
690 for (i = 0; (i < min(neig, 3)); i++)
694 case 0: axis = axis_x; break;
695 case 1: axis = axis_y; break;
696 case 2: axis = axis_z; break;
699 for (j = 0; j <= ibox[i]; j++)
701 axis[j] = min_eig[i] + j/bfac[i];
707 snew(MM, maxbox*maxbox);
708 for (i = 0; (i < ibox[0]); i++)
710 MM[i] = &(M[i*ibox[1]]);
714 for (i = 0; (i < nmap); i++)
716 Mmin = min(Mmin, map[0][i]);
717 Mmax = max(Mmax, map[0][i]);
720 for (i = 0; (i < len); i++)
724 for (i = 0; (i < nmap); i++)
726 index = gmx_nint(mapindex[i]);
729 gmx_fatal(FARGS, "Number of bins in file from -mdata option does not correspond to current analysis");
734 M[index] = map[0][i];
743 pick_minima(logf, ibox, neig, len, W);
748 flags = MAT_SPATIAL_X | MAT_SPATIAL_Y;
751 /* Dump to XPM file */
753 for (i = 0; (i < ibox[0]); i++)
755 snew(PP[i], ibox[1]);
756 for (j = 0; j < ibox[1]; j++)
758 PP[i][j] = P[i*ibox[1]+j];
760 WW[i] = &(W[i*ibox[1]]);
761 EE[i] = &(E[i*ibox[1]]);
762 SS[i] = &(S[i*ibox[1]]);
764 fp = ffopen(xpmP, "w");
765 write_xpm(fp, flags, "Probability Distribution", "", "PC1", "PC2",
766 ibox[0], ibox[1], axis_x, axis_y, PP, 0, Pmax, rlo, rhi, &nlevels);
768 fp = ffopen(xpm, "w");
769 write_xpm(fp, flags, "Gibbs Energy Landscape", "G (kJ/mol)", "PC1", "PC2",
770 ibox[0], ibox[1], axis_x, axis_y, WW, 0, gmax, rlo, rhi, &nlevels);
772 fp = ffopen(xpm2, "w");
773 write_xpm(fp, flags, "Enthalpy Landscape", "H (kJ/mol)", "PC1", "PC2",
774 ibox[0], ibox[1], axis_x, axis_y, EE,
775 emin ? *emin : Emin, emax ? *emax : Einf, rlo, rhi, &nlevels);
777 fp = ffopen(xpm3, "w");
778 write_xpm(fp, flags, "Entropy Landscape", "TDS (kJ/mol)", "PC1", "PC2",
779 ibox[0], ibox[1], axis_x, axis_y, SS, 0, Sinf, rlo, rhi, &nlevels);
783 fp = ffopen(xpm4, "w");
784 write_xpm(fp, flags, "Custom Landscape", mname, "PC1", "PC2",
785 ibox[0], ibox[1], axis_x, axis_y, MM, 0, Minf, rlo, rhi, &nlevels);
791 /* Dump to PDB file */
792 fp = ffopen(pdb, "w");
793 for (i = 0; (i < ibox[0]); i++)
795 xxx[XX] = 3*(i+0.5-ibox[0]/2);
796 for (j = 0; (j < ibox[1]); j++)
798 xxx[YY] = 3*(j+0.5-ibox[1]/2);
799 for (k = 0; (k < ibox[2]); k++)
801 xxx[ZZ] = 3*(k+0.5-ibox[2]/2);
802 index = index3(ibox, i, j, k);
805 fprintf(fp, "%-6s%5u %-4.4s%3.3s %4d %8.3f%8.3f%8.3f%6.2f%6.2f\n",
806 "ATOM", (index+1) %10000, "H", "H", (index+1)%10000,
807 xxx[XX], xxx[YY], xxx[ZZ], 1.0, W[index]);
813 write_xplor("out.xplor", W, ibox, min_eig, max_eig);
816 write_xplor("user.xplor", M, ibox, min_eig, max_eig);
818 nxyz[XX] = imin/(ibox[1]*ibox[2]);
819 nxyz[YY] = (imin-nxyz[XX]*ibox[1]*ibox[2])/ibox[2];
820 nxyz[ZZ] = imin % ibox[2];
821 for (i = 0; (i < ibox[0]); i++)
824 for (j = 0; (j < ibox[1]); j++)
826 WW[i][j] = W[index3(ibox, i, j, nxyz[ZZ])];
829 snew(buf, strlen(xpm)+4);
830 sprintf(buf, "%s", xpm);
831 sprintf(&buf[strlen(xpm)-4], "12.xpm");
832 fp = ffopen(buf, "w");
833 write_xpm(fp, flags, "Gibbs Energy Landscape", "W (kJ/mol)", "PC1", "PC2",
834 ibox[0], ibox[1], axis_x, axis_y, WW, 0, gmax, rlo, rhi, &nlevels);
836 for (i = 0; (i < ibox[0]); i++)
838 for (j = 0; (j < ibox[2]); j++)
840 WW[i][j] = W[index3(ibox, i, nxyz[YY], j)];
843 sprintf(&buf[strlen(xpm)-4], "13.xpm");
844 fp = ffopen(buf, "w");
845 write_xpm(fp, flags, "SHAM Energy Landscape", "kJ/mol", "PC1", "PC3",
846 ibox[0], ibox[2], axis_x, axis_z, WW, 0, gmax, rlo, rhi, &nlevels);
848 for (i = 0; (i < ibox[1]); i++)
850 for (j = 0; (j < ibox[2]); j++)
852 WW[i][j] = W[index3(ibox, nxyz[XX], i, j)];
855 sprintf(&buf[strlen(xpm)-4], "23.xpm");
856 fp = ffopen(buf, "w");
857 write_xpm(fp, flags, "SHAM Energy Landscape", "kJ/mol", "PC2", "PC3",
858 ibox[1], ibox[2], axis_y, axis_z, WW, 0, gmax, rlo, rhi, &nlevels);
869 static void ehisto(const char *fh, int n, real **enerT, const output_env_t oenv)
872 int i, j, k, nbin, blength;
874 real *T, bmin, bmax, bwidth;
882 for (j = 1; (j < n); j++)
884 for (k = 0; (k < nbin); k++)
886 if (T[k] == enerT[1][j])
895 T[nbin] = enerT[1][j];
898 bmin = min(enerT[0][j], bmin);
899 bmax = max(enerT[0][j], bmax);
902 blength = (bmax - bmin)/bwidth + 2;
904 for (i = 0; (i < nbin); i++)
906 snew(histo[i], blength);
908 for (j = 0; (j < n); j++)
910 k = (enerT[0][j]-bmin)/bwidth;
911 histo[bindex[j]][k]++;
913 fp = xvgropen(fh, "Energy distribution", "E (kJ/mol)", "", oenv);
914 for (j = 0; (j < blength); j++)
916 fprintf(fp, "%8.3f", bmin+j*bwidth);
917 for (k = 0; (k < nbin); k++)
919 fprintf(fp, " %6d", histo[k][j]);
926 int gmx_sham(int argc, char *argv[])
928 const char *desc[] = {
929 "[TT]g_sham[tt] makes multi-dimensional free-energy, enthalpy and entropy plots.",
930 "[TT]g_sham[tt] reads one or more [TT].xvg[tt] files and analyzes data sets.",
931 "The basic purpose of [TT]g_sham[tt] is to plot Gibbs free energy landscapes",
932 "(option [TT]-ls[tt])",
933 "by Bolzmann inverting multi-dimensional histograms (option [TT]-lp[tt]),",
935 "make enthalpy (option [TT]-lsh[tt]) and entropy (option [TT]-lss[tt])",
936 "plots. The histograms can be made for any quantities the user supplies.",
937 "A line in the input file may start with a time",
938 "(see option [TT]-time[tt]) and any number of [IT]y[it]-values may follow.",
939 "Multiple sets can also be",
940 "read when they are separated by & (option [TT]-n[tt]),",
941 "in this case only one [IT]y[it]-value is read from each line.",
942 "All lines starting with # and @ are skipped.",
944 "Option [TT]-ge[tt] can be used to supply a file with free energies",
945 "when the ensemble is not a Boltzmann ensemble, but needs to be biased",
946 "by this free energy. One free energy value is required for each",
947 "(multi-dimensional) data point in the [TT]-f[tt] input.",
949 "Option [TT]-ene[tt] can be used to supply a file with energies.",
950 "These energies are used as a weighting function in the single",
951 "histogram analysis method by Kumar et al. When temperatures",
952 "are supplied (as a second column in the file), an experimental",
953 "weighting scheme is applied. In addition the vales",
954 "are used for making enthalpy and entropy plots.",
956 "With option [TT]-dim[tt], dimensions can be gives for distances.",
957 "When a distance is 2- or 3-dimensional, the circumference or surface",
958 "sampled by two particles increases with increasing distance.",
959 "Depending on what one would like to show, one can choose to correct",
960 "the histogram and free-energy for this volume effect.",
961 "The probability is normalized by r and r^2 for dimensions of 2 and 3, ",
963 "A value of -1 is used to indicate an angle in degrees between two",
964 "vectors: a sin(angle) normalization will be applied.",
965 "[BB]Note[bb] that for angles between vectors the inner-product or cosine",
966 "is the natural quantity to use, as it will produce bins of the same",
969 static real tb = -1, te = -1, frac = 0.5, filtlen = 0, binwidth = 0.1;
970 static gmx_bool bHaveT = TRUE, bDer = FALSE, bSubAv = TRUE, bAverCorr = FALSE, bXYdy = FALSE;
971 static gmx_bool bEESEF = FALSE, bEENLC = FALSE, bEeFitAc = FALSE, bPower = FALSE;
972 static gmx_bool bShamEner = TRUE, bSham = TRUE;
973 static real Tref = 298.15, pmin = 0, ttol = 0, pmax = 0, gmax = 0, emin = 0, emax = 0;
974 static rvec nrdim = {1, 1, 1};
975 static rvec nrbox = {32, 32, 32};
976 static rvec xmin = {0, 0, 0}, xmax = {1, 1, 1};
977 static int nsets_in = 1, nb_min = 4, resol = 10, nlevels = 25;
978 static const char *mname = "";
980 { "-time", FALSE, etBOOL, {&bHaveT},
981 "Expect a time in the input" },
982 { "-b", FALSE, etREAL, {&tb},
983 "First time to read from set" },
984 { "-e", FALSE, etREAL, {&te},
985 "Last time to read from set" },
986 { "-ttol", FALSE, etREAL, {&ttol},
987 "Tolerance on time in appropriate units (usually ps)" },
988 { "-n", FALSE, etINT, {&nsets_in},
989 "Read this number of sets separated by lines containing only an ampersand" },
990 { "-d", FALSE, etBOOL, {&bDer},
991 "Use the derivative" },
992 { "-bw", FALSE, etREAL, {&binwidth},
993 "Binwidth for the distribution" },
994 { "-sham", FALSE, etBOOL, {&bSham},
995 "Turn off energy weighting even if energies are given" },
996 { "-tsham", FALSE, etREAL, {&Tref},
997 "Temperature for single histogram analysis" },
998 { "-pmin", FALSE, etREAL, {&pmin},
999 "Minimum probability. Anything lower than this will be set to zero" },
1000 { "-dim", FALSE, etRVEC, {nrdim},
1001 "Dimensions for distances, used for volume correction (max 3 values, dimensions > 3 will get the same value as the last)" },
1002 { "-ngrid", FALSE, etRVEC, {nrbox},
1003 "Number of bins for energy landscapes (max 3 values, dimensions > 3 will get the same value as the last)" },
1004 { "-xmin", FALSE, etRVEC, {xmin},
1005 "Minimum for the axes in energy landscape (see above for > 3 dimensions)" },
1006 { "-xmax", FALSE, etRVEC, {xmax},
1007 "Maximum for the axes in energy landscape (see above for > 3 dimensions)" },
1008 { "-pmax", FALSE, etREAL, {&pmax},
1009 "Maximum probability in output, default is calculate" },
1010 { "-gmax", FALSE, etREAL, {&gmax},
1011 "Maximum free energy in output, default is calculate" },
1012 { "-emin", FALSE, etREAL, {&emin},
1013 "Minimum enthalpy in output, default is calculate" },
1014 { "-emax", FALSE, etREAL, {&emax},
1015 "Maximum enthalpy in output, default is calculate" },
1016 { "-nlevels", FALSE, etINT, {&nlevels},
1017 "Number of levels for energy landscape" },
1018 { "-mname", FALSE, etSTR, {&mname},
1019 "Legend label for the custom landscape" },
1021 #define NPA asize(pa)
1024 int n, e_n, d_n, nlast, s, nset, e_nset, d_nset, i, j = 0, *idim, *ibox;
1025 real **val, **et_val, **dt_val, *t, *e_t, e_dt, d_dt, *d_t, dt, tot, error;
1027 double *av, *sig, cum1, cum2, cum3, cum4, db;
1028 const char *fn_ge, *fn_ene;
1030 gmx_large_int_t num_grid_points;
1033 { efXVG, "-f", "graph", ffREAD },
1034 { efXVG, "-ge", "gibbs", ffOPTRD },
1035 { efXVG, "-ene", "esham", ffOPTRD },
1036 { efXVG, "-dist", "ener", ffOPTWR },
1037 { efXVG, "-histo", "edist", ffOPTWR },
1038 { efNDX, "-bin", "bindex", ffOPTWR },
1039 { efXPM, "-lp", "prob", ffOPTWR },
1040 { efXPM, "-ls", "gibbs", ffOPTWR },
1041 { efXPM, "-lsh", "enthalpy", ffOPTWR },
1042 { efXPM, "-lss", "entropy", ffOPTWR },
1043 { efXPM, "-map", "map", ffOPTWR },
1044 { efPDB, "-ls3", "gibbs3", ffOPTWR },
1045 { efXVG, "-mdata", "mapdata", ffOPTRD },
1046 { efLOG, "-g", "shamlog", ffOPTWR }
1048 #define NFILE asize(fnm)
1053 parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_BE_NICE,
1054 NFILE, fnm, npargs, pa, asize(desc), desc, 0, NULL, &oenv);
1056 val = read_xvg_time(opt2fn("-f", NFILE, fnm), bHaveT,
1057 opt2parg_bSet("-b", npargs, pa), tb-ttol,
1058 opt2parg_bSet("-e", npargs, pa), te+ttol,
1059 nsets_in, &nset, &n, &dt, &t);
1060 printf("Read %d sets of %d points, dt = %g\n\n", nset, n, dt);
1062 fn_ge = opt2fn_null("-ge", NFILE, fnm);
1063 fn_ene = opt2fn_null("-ene", NFILE, fnm);
1065 if (fn_ge && fn_ene)
1067 gmx_fatal(FARGS, "Can not do free energy and energy corrections at the same time");
1070 if (fn_ge || fn_ene)
1072 et_val = read_xvg_time(fn_ge ? fn_ge : fn_ene, bHaveT,
1073 opt2parg_bSet("-b", npargs, pa), tb-ttol,
1074 opt2parg_bSet("-e", npargs, pa), te+ttol,
1075 1, &e_nset, &e_n, &e_dt, &e_t);
1080 gmx_fatal(FARGS, "Can only handle one free energy component in %s",
1086 if (e_nset != 1 && e_nset != 2)
1088 gmx_fatal(FARGS, "Can only handle one energy component or one energy and one T in %s",
1094 gmx_fatal(FARGS, "Number of energies (%d) does not match number of entries (%d) in %s", e_n, n, opt2fn("-f", NFILE, fnm));
1102 if (opt2fn_null("-mdata", NFILE, fnm) != NULL)
1104 dt_val = read_xvg_time(opt2fn("-mdata", NFILE, fnm), bHaveT,
1105 FALSE, tb, FALSE, te,
1106 nsets_in, &d_nset, &d_n, &d_dt, &d_t);
1109 gmx_fatal(FARGS, "Can only handle one mapping data column in %s",
1110 opt2fn("-mdata", NFILE, fnm));
1118 if (fn_ene && et_val)
1120 ehisto(opt2fn("-histo", NFILE, fnm), e_n, et_val, oenv);
1123 snew(idim, max(3, nset));
1124 snew(ibox, max(3, nset));
1125 snew(rmin, max(3, nset));
1126 snew(rmax, max(3, nset));
1127 for (i = 0; (i < min(3, nset)); i++)
1134 for (; (i < nset); i++)
1142 /* Check that the grid size is manageable. */
1143 num_grid_points = ibox[0];
1144 for (i = 1; i < nset; i++)
1146 gmx_large_int_t result;
1147 if (!check_int_multiply_for_overflow(num_grid_points, ibox[i], &result))
1150 "The number of dimensions and grid points is too large for this tool\n"
1151 "to handle with what it knows about the architecture upon which it\n"
1152 "is running. Use a different machine or consult the GROMACS mailing list.");
1154 num_grid_points = result;
1156 /* The number of grid points fits in a gmx_large_int_t. */
1158 do_sham(opt2fn("-dist", NFILE, fnm), opt2fn("-bin", NFILE, fnm),
1159 opt2fn("-lp", NFILE, fnm),
1160 opt2fn("-ls", NFILE, fnm), opt2fn("-lsh", NFILE, fnm),
1161 opt2fn("-lss", NFILE, fnm), opt2fn("-map", NFILE, fnm),
1162 opt2fn("-ls3", NFILE, fnm), opt2fn("-g", NFILE, fnm),
1163 n, nset, val, fn_ge != NULL, e_nset, et_val, d_n, d_t, dt_val, Tref,
1165 opt2parg_bSet("-emin", NPA, pa) ? &emin : NULL,
1166 opt2parg_bSet("-emax", NPA, pa) ? &emax : NULL,
1168 mname, bSham, idim, ibox,
1169 opt2parg_bSet("-xmin", NPA, pa), rmin,
1170 opt2parg_bSet("-xmax", NPA, pa), rmax);