Introduce SimulatorBuilder
[gromacs.git] / src / gromacs / gmxana / gmx_sham.cpp
blobb8550687a2bd52e933c377fd94dc95dfd97b5906
1 /*
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,2015,2016,2017,2018, 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.
37 #include "gmxpre.h"
39 #include <cassert>
40 #include <cmath>
41 #include <cstdlib>
42 #include <cstring>
44 #include <algorithm>
46 #include "gromacs/commandline/pargs.h"
47 #include "gromacs/fileio/matio.h"
48 #include "gromacs/fileio/pdbio.h"
49 #include "gromacs/fileio/xvgr.h"
50 #include "gromacs/gmxana/gmx_ana.h"
51 #include "gromacs/gmxana/gstat.h"
52 #include "gromacs/math/functions.h"
53 #include "gromacs/math/units.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/topology/topology.h"
56 #include "gromacs/utility/arraysize.h"
57 #include "gromacs/utility/fatalerror.h"
58 #include "gromacs/utility/futil.h"
59 #include "gromacs/utility/smalloc.h"
62 static int index2(const int *ibox, int x, int y)
64 return (ibox[1]*x+y);
67 static int index3(const int *ibox, int x, int y, int z)
69 return (ibox[2]*(ibox[1]*x+y)+z);
72 static int64_t indexn(int ndim, const int *ibox, const int *nxyz)
74 int64_t d, dd;
75 int k, kk;
77 /* Compute index in 1-D array */
78 d = 0;
79 for (k = 0; (k < ndim); k++)
81 dd = nxyz[k];
82 for (kk = k+1; (kk < ndim); kk++)
84 dd = dd*ibox[kk];
86 d += dd;
88 return d;
91 typedef struct {
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 */
98 real * ed; /* data */
99 } XplorMap;
101 static void lo_write_xplor(XplorMap * map, const char * file)
103 FILE * fp;
104 int z, i, j, n;
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");
122 z = map->dmin[2];
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]);
135 fprintf(fp, "\n");
138 fprintf(fp, " -9999\n");
139 gmx_ffclose(fp);
142 static void write_xplor(const char *file, const real *data, int *ibox, const real dmin[], const real dmax[])
144 XplorMap *xm;
145 int i, j, k, n;
147 snew(xm, 1);
148 xm->Nx = ibox[XX];
149 xm->Ny = ibox[YY];
150 xm->Nz = ibox[ZZ];
151 snew(xm->ed, xm->Nx*xm->Ny*xm->Nz);
152 n = 0;
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);
175 sfree(xm->ed);
176 sfree(xm);
179 static void normalize_p_e(int len, double *P, const int *nbin, real *E, real pmin)
181 int i;
182 double Ptot = 0;
184 for (i = 0; (i < len); i++)
186 Ptot += P[i];
187 if (nbin[i] > 0)
189 E[i] = E[i]/nbin[i];
192 printf("Ptot = %g\n", Ptot);
193 for (i = 0; (i < len); i++)
195 P[i] = P[i]/Ptot;
196 /* Have to check for pmin after normalizing to prevent "stretching"
197 * the energies.
199 if (P[i] < pmin)
201 P[i] = 0;
206 typedef struct {
207 int64_t index;
208 real ener;
209 } t_minimum;
211 static int comp_minima(const void *a, const void *b)
213 const t_minimum *ma = reinterpret_cast<const t_minimum*>(a);
214 const t_minimum *mb = reinterpret_cast<const t_minimum*>(b);;
216 if (ma->ener < mb->ener)
218 return -1;
220 else if (ma->ener > mb->ener)
222 return 1;
224 else
226 return 0;
230 static inline
231 void print_minimum(FILE *fp, int num, const t_minimum *min)
233 fprintf(fp,
234 "Minimum %d at index " "%" PRId64 " energy %10.3f\n",
235 num, min->index, min->ener);
238 static inline
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;
246 static inline
247 gmx_bool is_local_minimum_from_below(const t_minimum *this_min,
248 int dimension_index,
249 int dimension_min,
250 int neighbour_index,
251 const real *W)
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. */
259 static inline
260 gmx_bool is_local_minimum_from_above(const t_minimum *this_min,
261 int dimension_index,
262 int dimension_max,
263 int neighbour_index,
264 const real *W)
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[])
274 FILE *fp;
275 int i, j, k, nmin;
276 t_minimum *mm, this_min;
277 int *this_point;
278 int loopmax, loopcounter;
280 snew(mm, len);
281 nmin = 0;
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. */
286 switch (ndim)
288 case 0:
289 /* This is probably impossible to reach anyway. */
290 break;
291 case 2:
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);
305 nmin++;
309 break;
310 case 3:
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);
328 nmin++;
333 break;
334 default:
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
339 * start. */
340 snew(this_point, ndim);
342 /* Determine the number of points of the ndim-dimensional
343 * grid. */
344 loopmax = ibox[0];
345 for (i = 1; i < ndim; i++)
347 loopmax *= ibox[i];
350 loopcounter = 0;
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
360 * dimension? */
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];
367 this_point[i]--;
368 bMin = bMin &&
369 is_local_minimum_from_below(&this_min, index, 0, indexn(ndim, ibox, this_point), W);
370 this_point[i] += 2;
371 bMin = bMin &&
372 is_local_minimum_from_above(&this_min, index, ibox[i]-1, indexn(ndim, ibox, this_point), W);
373 this_point[i]--;
375 if (bMin)
377 add_minimum(fp, nmin, &this_min, mm);
378 nmin++;
381 /* update global loop counter */
382 loopcounter++;
384 /* Avoid underflow of this_point[i] */
385 if (loopmax > loopcounter)
387 /* update this_point non-recursively */
388 i = ndim-1;
389 this_point[i]++;
390 while (ibox[i] == this_point[i])
392 this_point[i] = 0;
393 i--;
394 /* this_point[i] cannot underflow because
395 * loopmax > loopcounter. */
396 this_point[i]++;
401 sfree(this_point);
402 break;
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]);
410 gmx_ffclose(fp);
411 sfree(mm);
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,
417 const char *logf,
418 int n, int neig, real **eig,
419 gmx_bool bGE, int nenerT, real **enerT,
420 real Tref,
421 real pmax, real gmax,
422 const real *emin, const real *emax, int nlevels, real pmin,
423 const int *idim, int *ibox,
424 gmx_bool bXmin, real *xmin, gmx_bool bXmax, real *xmax)
426 FILE *fp;
427 real *min_eig, *max_eig;
428 real *axis_x, *axis_y, *axis_z, *axis = nullptr;
429 double *P;
430 real **PP, *W, *E, **WW, **EE, *S, **SS, *M, *bE;
431 rvec xxx;
432 char *buf;
433 double *bfac, efac, bref, Pmax, Wmin, Wmax, Winf, Emin, Emax, Einf, Smin, Smax, Sinf;
434 real *delta;
435 int i, j, k, imin, len, index, *nbin, *bindex, bi;
436 int *nxyz, maxbox;
437 t_blocka *b;
438 gmx_bool bOutside;
439 unsigned int flags;
440 t_rgb rlo = { 0, 0, 0 };
441 t_rgb rhi = { 1, 1, 1 };
443 /* Determine extremes for the eigenvectors */
444 snew(min_eig, neig);
445 snew(max_eig, neig);
446 snew(nxyz, neig);
447 snew(bfac, neig);
448 snew(delta, neig);
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] = std::min(min_eig[i], eig[i][j]);
457 max_eig[i] = std::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.
463 if (bXmax)
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];
471 else
473 max_eig[i] += delta[i];
476 if (bXmin)
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];
484 else
486 min_eig[i] -= delta[i];
488 bfac[i] = ibox[i]/(max_eig[i]-min_eig[i]);
490 /* Do the binning */
491 bref = 1/(BOLTZ*Tref);
492 snew(bE, n);
493 if (bGE || nenerT == 2)
495 assert(enerT);
496 Emin = 1e8;
497 for (j = 0; (j < n); j++)
499 if (bGE)
501 bE[j] = bref*enerT[0][j];
503 else
505 bE[j] = (bref - 1/(BOLTZ*enerT[1][j]))*enerT[0][j];
507 Emin = std::min(Emin, static_cast<double>(bE[j]));
510 else
512 Emin = 0;
514 len = 1;
515 for (i = 0; (i < neig); i++)
517 len = len*ibox[i];
519 printf("There are %d bins in the %d-dimensional histogram. Beta-Emin = %g\n",
520 len, neig, Emin);
521 snew(P, len);
522 snew(W, len);
523 snew(E, len);
524 snew(S, len);
525 snew(M, len);
526 snew(nbin, len);
527 snew(bindex, n);
530 /* Loop over projections */
531 for (j = 0; (j < n); j++)
533 /* Loop over dimensions */
534 bOutside = FALSE;
535 for (i = 0; (i < neig); i++)
537 nxyz[i] = static_cast<int>(bfac[i]*(eig[i][j]-min_eig[i]));
538 if (nxyz[i] < 0 || nxyz[i] >= ibox[i])
540 bOutside = TRUE;
543 if (!bOutside)
545 index = indexn(neig, ibox, nxyz);
546 range_check(index, 0, len);
547 /* Compute the exponential factor */
548 if (enerT)
550 efac = std::exp(-bE[j]+Emin);
552 else
554 efac = 1;
556 /* Apply the bin volume correction for a multi-dimensional distance */
557 for (i = 0; i < neig; i++)
559 if (idim[i] == 2)
561 efac /= eig[i][j];
563 else if (idim[i] == 3)
565 efac /= gmx::square(eig[i][j]);
567 else if (idim[i] == -1)
569 efac /= std::sin(DEG2RAD*eig[i][j]);
572 /* Update the probability */
573 P[index] += efac;
574 /* Update the energy */
575 if (enerT)
577 E[index] += enerT[0][j];
579 /* Statistics: which "structure" in which bin */
580 nbin[index]++;
581 bindex[j] = index;
584 /* Normalize probability */
585 normalize_p_e(len, P, nbin, E, pmin);
586 Pmax = 0;
587 /* Compute boundaries for the Free energy */
588 Wmin = 1e8;
589 imin = -1;
590 Wmax = -1e8;
591 /* Recompute Emin: it may have changed due to averaging */
592 Emin = 1e8;
593 Emax = -1e8;
594 for (i = 0; (i < len); i++)
596 if (P[i] != 0)
598 Pmax = std::max(P[i], Pmax);
599 W[i] = -BOLTZ*Tref*std::log(P[i]);
600 if (W[i] < Wmin)
602 Wmin = W[i];
603 imin = i;
605 Emin = std::min(static_cast<double>(E[i]), Emin);
606 Emax = std::max(static_cast<double>(E[i]), Emax);
607 Wmax = std::max(static_cast<double>(W[i]), Wmax);
610 if (pmax > 0)
612 Pmax = pmax;
614 if (gmax > 0)
616 Wmax = gmax;
618 else
620 Wmax -= Wmin;
622 Winf = Wmax+1;
623 Einf = Emax+1;
624 Smin = Emin-Wmax;
625 Smax = Emax-Smin;
626 Sinf = Smax+1;
627 /* Write out the free energy as a function of bin index */
628 fp = gmx_ffopen(fn, "w");
629 for (i = 0; (i < len); i++)
631 if (P[i] != 0)
633 W[i] -= Wmin;
634 S[i] = E[i]-W[i]-Smin;
635 fprintf(fp, "%5d %10.5e %10.5e %10.5e\n", i, W[i], E[i], S[i]);
637 else
639 W[i] = Winf;
640 E[i] = Einf;
641 S[i] = Sinf;
644 gmx_ffclose(fp);
645 /* Organize the structures in the bins */
646 snew(b, 1);
647 snew(b->index, len+1);
648 snew(b->a, n);
649 b->index[0] = 0;
650 for (i = 0; (i < len); i++)
652 b->index[i+1] = b->index[i]+nbin[i];
653 nbin[i] = 0;
655 for (i = 0; (i < n); i++)
657 bi = bindex[i];
658 b->a[b->index[bi]+nbin[bi]] = i;
659 nbin[bi]++;
661 /* Consistency check */
662 /* This no longer applies when we allow the plot to be smaller
663 than the sampled space.
664 for(i=0; (i<len); i++) {
665 if (nbin[i] != (b->index[i+1] - b->index[i]))
666 gmx_fatal(FARGS,"nbin[%d] = %d, should be %d",i,nbin[i],
667 b->index[i+1] - b->index[i]);
670 /* Write the index file */
671 fp = gmx_ffopen(ndx, "w");
672 for (i = 0; (i < len); i++)
674 if (nbin[i] > 0)
676 fprintf(fp, "[ %d ]\n", i);
677 for (j = b->index[i]; (j < b->index[i+1]); j++)
679 fprintf(fp, "%d\n", b->a[j]+1);
683 gmx_ffclose(fp);
684 snew(axis_x, ibox[0]+1);
685 snew(axis_y, ibox[1]+1);
686 snew(axis_z, ibox[2]+1);
687 maxbox = std::max(ibox[0], std::max(ibox[1], ibox[2]));
688 snew(PP, maxbox*maxbox);
689 snew(WW, maxbox*maxbox);
690 snew(EE, maxbox*maxbox);
691 snew(SS, maxbox*maxbox);
692 for (i = 0; (i < std::min(neig, 3)); i++)
694 switch (i)
696 case 0: axis = axis_x; break;
697 case 1: axis = axis_y; break;
698 case 2: axis = axis_z; break;
699 default: break;
701 for (j = 0; j <= ibox[i]; j++)
703 axis[j] = min_eig[i] + j/bfac[i];
707 pick_minima(logf, ibox, neig, len, W);
708 if (gmax <= 0)
710 gmax = Winf;
712 flags = MAT_SPATIAL_X | MAT_SPATIAL_Y;
713 if (neig == 2)
715 /* Dump to XPM file */
716 snew(PP, ibox[0]);
717 for (i = 0; (i < ibox[0]); i++)
719 snew(PP[i], ibox[1]);
720 for (j = 0; j < ibox[1]; j++)
722 PP[i][j] = P[i*ibox[1]+j];
724 WW[i] = &(W[i*ibox[1]]);
725 EE[i] = &(E[i*ibox[1]]);
726 SS[i] = &(S[i*ibox[1]]);
728 fp = gmx_ffopen(xpmP, "w");
729 write_xpm(fp, flags, "Probability Distribution", "", "PC1", "PC2",
730 ibox[0], ibox[1], axis_x, axis_y, PP, 0, Pmax, rlo, rhi, &nlevels);
731 gmx_ffclose(fp);
732 fp = gmx_ffopen(xpm, "w");
733 write_xpm(fp, flags, "Gibbs Energy Landscape", "G (kJ/mol)", "PC1", "PC2",
734 ibox[0], ibox[1], axis_x, axis_y, WW, 0, gmax, rlo, rhi, &nlevels);
735 gmx_ffclose(fp);
736 fp = gmx_ffopen(xpm2, "w");
737 write_xpm(fp, flags, "Enthalpy Landscape", "H (kJ/mol)", "PC1", "PC2",
738 ibox[0], ibox[1], axis_x, axis_y, EE,
739 emin ? *emin : Emin, emax ? *emax : Einf, rlo, rhi, &nlevels);
740 gmx_ffclose(fp);
741 fp = gmx_ffopen(xpm3, "w");
742 write_xpm(fp, flags, "Entropy Landscape", "TDS (kJ/mol)", "PC1", "PC2",
743 ibox[0], ibox[1], axis_x, axis_y, SS, 0, Sinf, rlo, rhi, &nlevels);
744 gmx_ffclose(fp);
746 else if (neig == 3)
748 /* Dump to PDB file */
749 fp = gmx_ffopen(pdb, "w");
750 for (i = 0; (i < ibox[0]); i++)
752 xxx[XX] = 3*i+1.5*(1-ibox[0]);
753 for (j = 0; (j < ibox[1]); j++)
755 xxx[YY] = 3*j+1.5*(1-ibox[1]);
756 for (k = 0; (k < ibox[2]); k++)
758 xxx[ZZ] = 3*k+1.5*(1-ibox[2]);
759 index = index3(ibox, i, j, k);
760 if (P[index] > 0)
762 fprintf(fp, "%-6s%5d %-4.4s%3.3s %4d %8.3f%8.3f%8.3f%6.2f%6.2f\n",
763 "ATOM", (index+1) %10000, "H", "H", (index+1)%10000,
764 xxx[XX], xxx[YY], xxx[ZZ], 1.0, W[index]);
769 gmx_ffclose(fp);
770 write_xplor("out.xplor", W, ibox, min_eig, max_eig);
771 nxyz[XX] = imin/(ibox[1]*ibox[2]);
772 nxyz[YY] = (imin-nxyz[XX]*ibox[1]*ibox[2])/ibox[2];
773 nxyz[ZZ] = imin % ibox[2];
774 for (i = 0; (i < ibox[0]); i++)
776 snew(WW[i], maxbox);
777 for (j = 0; (j < ibox[1]); j++)
779 WW[i][j] = W[index3(ibox, i, j, nxyz[ZZ])];
782 snew(buf, std::strlen(xpm)+4);
783 sprintf(buf, "%s", xpm);
784 sprintf(&buf[std::strlen(xpm)-4], "12.xpm");
785 fp = gmx_ffopen(buf, "w");
786 write_xpm(fp, flags, "Gibbs Energy Landscape", "W (kJ/mol)", "PC1", "PC2",
787 ibox[0], ibox[1], axis_x, axis_y, WW, 0, gmax, rlo, rhi, &nlevels);
788 gmx_ffclose(fp);
789 for (i = 0; (i < ibox[0]); i++)
791 for (j = 0; (j < ibox[2]); j++)
793 WW[i][j] = W[index3(ibox, i, nxyz[YY], j)];
796 sprintf(&buf[std::strlen(xpm)-4], "13.xpm");
797 fp = gmx_ffopen(buf, "w");
798 write_xpm(fp, flags, "SHAM Energy Landscape", "kJ/mol", "PC1", "PC3",
799 ibox[0], ibox[2], axis_x, axis_z, WW, 0, gmax, rlo, rhi, &nlevels);
800 gmx_ffclose(fp);
801 for (i = 0; (i < ibox[1]); i++)
803 for (j = 0; (j < ibox[2]); j++)
805 WW[i][j] = W[index3(ibox, nxyz[XX], i, j)];
808 sprintf(&buf[std::strlen(xpm)-4], "23.xpm");
809 fp = gmx_ffopen(buf, "w");
810 write_xpm(fp, flags, "SHAM Energy Landscape", "kJ/mol", "PC2", "PC3",
811 ibox[1], ibox[2], axis_y, axis_z, WW, 0, gmax, rlo, rhi, &nlevels);
812 gmx_ffclose(fp);
813 sfree(buf);
817 static void ehisto(const char *fh, int n, real **enerT, const gmx_output_env_t *oenv)
819 FILE *fp;
820 int i, j, k, nbin, blength;
821 int *bindex;
822 real *T, bmin, bmax, bwidth;
823 int **histo;
825 bmin = 1e8;
826 bmax = -1e8;
827 snew(bindex, n);
828 snew(T, n);
829 nbin = 0;
830 for (j = 1; (j < n); j++)
832 for (k = 0; (k < nbin); k++)
834 if (T[k] == enerT[1][j])
836 bindex[j] = k;
837 break;
840 if (k == nbin)
842 bindex[j] = nbin;
843 T[nbin] = enerT[1][j];
844 nbin++;
846 bmin = std::min(enerT[0][j], bmin);
847 bmax = std::max(enerT[0][j], bmax);
849 bwidth = 1.0;
850 blength = static_cast<int>((bmax - bmin)/bwidth + 2);
851 snew(histo, nbin);
852 for (i = 0; (i < nbin); i++)
854 snew(histo[i], blength);
856 for (j = 0; (j < n); j++)
858 k = static_cast<int>((enerT[0][j]-bmin)/bwidth);
859 histo[bindex[j]][k]++;
861 fp = xvgropen(fh, "Energy distribution", "E (kJ/mol)", "", oenv);
862 for (j = 0; (j < blength); j++)
864 fprintf(fp, "%8.3f", bmin+j*bwidth);
865 for (k = 0; (k < nbin); k++)
867 fprintf(fp, " %6d", histo[k][j]);
869 fprintf(fp, "\n");
871 xvgrclose(fp);
874 int gmx_sham(int argc, char *argv[])
876 const char *desc[] = {
877 "[THISMODULE] makes multi-dimensional free-energy, enthalpy and entropy plots.",
878 "[THISMODULE] reads one or more [REF].xvg[ref] files and analyzes data sets.",
879 "The basic purpose of [THISMODULE] is to plot Gibbs free energy landscapes",
880 "(option [TT]-ls[tt])",
881 "by Bolzmann inverting multi-dimensional histograms (option [TT]-lp[tt]),",
882 "but it can also",
883 "make enthalpy (option [TT]-lsh[tt]) and entropy (option [TT]-lss[tt])",
884 "plots. The histograms can be made for any quantities the user supplies.",
885 "A line in the input file may start with a time",
886 "(see option [TT]-time[tt]) and any number of [IT]y[it]-values may follow.",
887 "Multiple sets can also be",
888 "read when they are separated by & (option [TT]-n[tt]),",
889 "in this case only one [IT]y[it]-value is read from each line.",
890 "All lines starting with # and @ are skipped.",
891 "[PAR]",
892 "Option [TT]-ge[tt] can be used to supply a file with free energies",
893 "when the ensemble is not a Boltzmann ensemble, but needs to be biased",
894 "by this free energy. One free energy value is required for each",
895 "(multi-dimensional) data point in the [TT]-f[tt] input.",
896 "[PAR]",
897 "Option [TT]-ene[tt] can be used to supply a file with energies.",
898 "These energies are used as a weighting function in the single",
899 "histogram analysis method by Kumar et al. When temperatures",
900 "are supplied (as a second column in the file), an experimental",
901 "weighting scheme is applied. In addition the vales",
902 "are used for making enthalpy and entropy plots.",
903 "[PAR]",
904 "With option [TT]-dim[tt], dimensions can be gives for distances.",
905 "When a distance is 2- or 3-dimensional, the circumference or surface",
906 "sampled by two particles increases with increasing distance.",
907 "Depending on what one would like to show, one can choose to correct",
908 "the histogram and free-energy for this volume effect.",
909 "The probability is normalized by r and r^2 for dimensions of 2 and 3, ",
910 "respectively.",
911 "A value of -1 is used to indicate an angle in degrees between two",
912 "vectors: a sin(angle) normalization will be applied.",
913 "[BB]Note[bb] that for angles between vectors the inner-product or cosine",
914 "is the natural quantity to use, as it will produce bins of the same",
915 "volume."
917 static real tb = -1, te = -1;
918 static gmx_bool bHaveT = TRUE, bDer = FALSE;
919 static gmx_bool 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, nlevels = 25;
925 t_pargs pa[] = {
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)
965 int n, e_n, nset, e_nset = 0, i, *idim, *ibox;
966 real **val, **et_val, *t, *e_t, e_dt, dt;
967 real *rmin, *rmax;
968 const char *fn_ge, *fn_ene;
969 gmx_output_env_t *oenv;
970 int64_t num_grid_points;
972 t_filenm fnm[] = {
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)
988 int npargs;
990 npargs = asize(pa);
991 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW,
992 NFILE, fnm, npargs, pa, asize(desc), desc, 0, nullptr, &oenv))
994 return 0;
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);
1017 if (fn_ge)
1019 if (e_nset != 1)
1021 gmx_fatal(FARGS, "Can only handle one free energy component in %s",
1022 fn_ge);
1025 else
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",
1030 fn_ene);
1033 if (e_n != n)
1035 gmx_fatal(FARGS, "Number of energies (%d) does not match number of entries (%d) in %s", e_n, n, opt2fn("-f", NFILE, fnm));
1038 else
1040 et_val = nullptr;
1043 if (fn_ene && et_val)
1045 ehisto(opt2fn("-histo", NFILE, fnm), e_n, et_val, oenv);
1048 snew(idim, std::max(3, nset));
1049 snew(ibox, std::max(3, nset));
1050 snew(rmin, std::max(3, nset));
1051 snew(rmax, std::max(3, nset));
1052 for (i = 0; (i < std::min(3, nset)); i++)
1054 idim[i] = static_cast<int>(nrdim[i]);
1055 ibox[i] = static_cast<int>(nrbox[i]);
1056 rmin[i] = xmin[i];
1057 rmax[i] = xmax[i];
1059 for (; (i < nset); i++)
1061 idim[i] = static_cast<int>(nrdim[2]);
1062 ibox[i] = static_cast<int>(nrbox[2]);
1063 rmin[i] = xmin[2];
1064 rmax[i] = xmax[2];
1067 /* Check that the grid size is manageable. */
1068 num_grid_points = ibox[0];
1069 for (i = 1; i < nset; i++)
1071 int64_t result;
1072 if (!check_int_multiply_for_overflow(num_grid_points, ibox[i], &result))
1074 gmx_fatal(FARGS,
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 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 != nullptr, e_nset, et_val, Tref,
1087 pmax, gmax,
1088 opt2parg_bSet("-emin", NPA, pa) ? &emin : nullptr,
1089 opt2parg_bSet("-emax", NPA, pa) ? &emax : nullptr,
1090 nlevels, pmin,
1091 idim, ibox,
1092 opt2parg_bSet("-xmin", NPA, pa), rmin,
1093 opt2parg_bSet("-xmax", NPA, pa), rmax);
1095 return 0;