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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
48 #include "gmx_fatal.h"
63 static int index2(int *ibox,int x,int y)
68 static int index3(int *ibox,int x,int y,int z)
70 return (ibox[2]*(ibox[1]*x+y)+z);
73 static gmx_large_int_t indexn(int ndim, const int *ibox, const int *nxyz)
75 gmx_large_int_t d, dd;
78 /* Compute index in 1-D array */
80 for(k=0; (k<ndim); k++) {
82 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++) {
122 fprintf(fp, "%8d\n", z) ;
123 for (i = 0; i < map->Nx*map->Ny; i += 6) {
124 for (j = 0; j < 6; j++)
125 if (i+j < map->Nx*map->Ny)
126 fprintf(fp, "%12.5E", map->ed[n*map->Nx*map->Ny+i+j]);
130 fprintf(fp, " -9999\n") ;
134 static void write_xplor(const char *file,real *data,int *ibox,real dmin[],real dmax[])
143 snew(xm->ed,xm->Nx*xm->Ny*xm->Nz);
145 for(k=0; (k<xm->Nz); k++)
146 for(j=0; (j<xm->Ny); j++)
147 for(i=0; (i<xm->Nx); i++)
148 xm->ed[n++] = data[index3(ibox,i,j,k)];
149 xm->cell[0] = dmax[XX]-dmin[XX];
150 xm->cell[1] = dmax[YY]-dmin[YY];
151 xm->cell[2] = dmax[ZZ]-dmin[ZZ];
152 xm->cell[3] = xm->cell[4] = xm->cell[5] = 90;
154 clear_ivec(xm->dmin);
155 xm->dmax[XX] = ibox[XX]-1;
156 xm->dmax[YY] = ibox[YY]-1;
157 xm->dmax[ZZ] = ibox[ZZ]-1;
159 lo_write_xplor(xm,file);
165 static void normalize_p_e(int len,double *P,int *nbin,real *E,real pmin)
170 for(i=0; (i<len); i++) {
175 printf("Ptot = %g\n",Ptot);
176 for(i=0; (i<len); i++) {
178 /* Have to check for pmin after normalizing to prevent "stretching"
187 gmx_large_int_t index;
191 static int comp_minima(const void *a,const void *b)
193 t_minimum *ma = (t_minimum *) a;
194 t_minimum *mb = (t_minimum *) b;
196 if (ma->ener < mb->ener)
198 else if (ma->ener > mb->ener)
205 void print_minimum(FILE *fp, int num, const t_minimum *min)
208 "Minimum %d at index " gmx_large_int_pfmt " energy %10.3f\n",
209 num, min->index, min->ener);
213 void add_minimum(FILE *fp, int num, const t_minimum *min, t_minimum *mm)
215 print_minimum(fp, num, min);
216 mm[num].index = min->index;
217 mm[num].ener = min->ener;
221 gmx_bool is_local_minimum_from_below(const t_minimum *this_min,
227 return ((dimension_index == dimension_min) ||
228 ((dimension_index > dimension_min) &&
229 (this_min->ener < W[neighbour_index])));
230 /* Note over/underflow within W cannot occur. */
234 gmx_bool is_local_minimum_from_above(const t_minimum *this_min,
240 return ((dimension_index == dimension_max) ||
241 ((dimension_index < dimension_max) &&
242 (this_min->ener < W[neighbour_index])));
243 /* Note over/underflow within W cannot occur. */
246 static void pick_minima(const char *logfile,int *ibox,int ndim,int len,real W[])
250 t_minimum *mm, this_min;
252 int loopmax, loopcounter;
256 fp = ffopen(logfile,"w");
257 /* Loop over each element in the array of dimenion ndim seeking
258 * minima with respect to every dimension. Specialized loops for
259 * speed with ndim == 2 and ndim == 3. */
263 /* This is probably impossible to reach anyway. */
266 for(i=0; (i<ibox[0]); i++) {
267 for(j=0; (j<ibox[1]); j++) {
268 /* Get the index of this point in the flat array */
269 this_min.index = index2(ibox,i,j);
270 this_min.ener = W[this_min.index];
271 if (is_local_minimum_from_below(&this_min, i, 0, index2(ibox,i-1,j ), W) &&
272 is_local_minimum_from_above(&this_min, i, ibox[0]-1, index2(ibox,i+1,j ), W) &&
273 is_local_minimum_from_below(&this_min, j, 0, index2(ibox,i ,j-1), W) &&
274 is_local_minimum_from_above(&this_min, j, ibox[1]-1, index2(ibox,i ,j+1), W))
276 add_minimum(fp, nmin, &this_min, mm);
283 for(i=0; (i<ibox[0]); i++) {
284 for(j=0; (j<ibox[1]); j++) {
285 for(k=0; (k<ibox[2]); k++) {
286 /* Get the index of this point in the flat array */
287 this_min.index = index3(ibox,i,j,k);
288 this_min.ener = W[this_min.index];
289 if (is_local_minimum_from_below(&this_min, i, 0, index3(ibox,i-1,j ,k ), W) &&
290 is_local_minimum_from_above(&this_min, i, ibox[0]-1, index3(ibox,i+1,j ,k ), W) &&
291 is_local_minimum_from_below(&this_min, j, 0, index3(ibox,i ,j-1,k ), W) &&
292 is_local_minimum_from_above(&this_min, j, ibox[1]-1, index3(ibox,i ,j+1,k ), W) &&
293 is_local_minimum_from_below(&this_min, k, 0, index3(ibox,i ,j ,k-1), W) &&
294 is_local_minimum_from_above(&this_min, k, ibox[2]-1, index3(ibox,i ,j ,k+1), W))
296 add_minimum(fp, nmin, &this_min, mm);
304 /* Note this treats ndim == 1 and ndim > 3 */
306 /* Set up an ndim-dimensional vector to loop over the points
307 * on the grid. (0,0,0, ... 0) is an acceptable place to
309 snew(this_point, ndim);
311 /* Determine the number of points of the ndim-dimensional
314 for (i = 1; i < ndim; i++)
320 while (loopmax > loopcounter)
322 gmx_bool bMin = TRUE;
324 /* Get the index of this_point in the flat array */
325 this_min.index = indexn(ndim, ibox, this_point);
326 this_min.ener = W[this_min.index];
328 /* Is this_point a minimum from above and below in each
330 for (i = 0; bMin && (i < ndim); i++)
332 /* Save the index of this_point within the curent
333 * dimension so we can change that index in the
334 * this_point array for use with indexn(). */
335 int index = this_point[i];
338 is_local_minimum_from_below(&this_min, index, 0, indexn(ndim, ibox, this_point), W);
341 is_local_minimum_from_above(&this_min, index, ibox[i]-1, indexn(ndim, ibox, this_point), W);
346 add_minimum(fp, nmin, &this_min, mm);
350 /* update global loop counter */
353 /* Avoid underflow of this_point[i] */
354 if (loopmax > loopcounter)
356 /* update this_point non-recursively */
359 while (ibox[i] == this_point[i])
363 /* this_point[i] cannot underflow because
364 * loopmax > loopcounter. */
373 qsort(mm,nmin,sizeof(mm[0]),comp_minima);
374 fprintf(fp,"Minima sorted after energy\n");
375 for(i=0; (i<nmin); i++)
377 print_minimum(fp, i, &mm[i]);
383 static void do_sham(const char *fn,const char *ndx,
384 const char *xpmP,const char *xpm,const char *xpm2,
385 const char *xpm3,const char *xpm4,const char *pdb,
387 int n,int neig,real **eig,
388 gmx_bool bGE,int nenerT,real **enerT,
389 int nmap,real *mapindex,real **map,
392 real *emin,real *emax,int nlevels,real pmin,
393 const char *mname,gmx_bool bSham,int *idim,int *ibox,
394 gmx_bool bXmin,real *xmin,gmx_bool bXmax,real *xmax)
397 real *min_eig,*max_eig;
398 real *axis_x,*axis_y,*axis_z,*axis=NULL;
400 real **PP,*W,*E,**WW,**EE,*S,**SS,*M,**MM,*bE;
403 double *bfac,efac,bref,Pmax,Wmin,Wmax,Winf,Emin,Emax,Einf,Smin,Smax,Sinf,Mmin,Mmax,Minf;
405 int i,j,k,imin,len,index,d,*nbin,*bindex,bi;
410 t_rgb rlo = { 0, 0, 0 };
411 t_rgb rhi = { 1, 1, 1 };
413 /* Determine extremes for the eigenvectors */
420 for(i=0; (i<neig); i++) {
421 /* Check for input constraints */
422 min_eig[i] = max_eig[i] = eig[i][0];
423 for(j=0; (j<n); j++) {
424 min_eig[i] = min(min_eig[i],eig[i][j]);
425 max_eig[i] = max(max_eig[i],eig[i][j]);
426 delta[i] = (max_eig[i]-min_eig[i])/(2.0*ibox[i]);
428 /* Add some extra space, half a bin on each side, unless the
429 * user has set the limits.
432 if (max_eig[i] > xmax[i]) {
433 gmx_warning("Your xmax[%d] value %f is smaller than the largest data point %f",i,xmax[i],max_eig[i]);
435 max_eig[i] = xmax[i];
438 max_eig[i] += delta[i];
441 if (min_eig[i] < xmin[i]) {
442 gmx_warning("Your xmin[%d] value %f is larger than the smallest data point %f",i,xmin[i],min_eig[i]);
444 min_eig[i] = xmin[i];
447 min_eig[i] -= delta[i];
448 bfac[i] = ibox[i]/(max_eig[i]-min_eig[i]);
451 bref = 1/(BOLTZ*Tref);
453 if (bGE || nenerT==2) {
455 for(j=0; (j<n); j++) {
457 bE[j] = bref*enerT[0][j];
459 bE[j] = (bref - 1/(BOLTZ*enerT[1][j]))*enerT[0][j];
460 Emin = min(Emin,bE[j]);
466 for(i=0; (i<neig); i++)
468 printf("There are %d bins in the %d-dimensional histogram. Beta-Emin = %g\n",
479 /* Loop over projections */
480 for(j=0; (j<n); j++) {
481 /* Loop over dimensions */
483 for(i=0; (i<neig); i++) {
484 nxyz[i] = bfac[i]*(eig[i][j]-min_eig[i]);
485 if (nxyz[i] < 0 || nxyz[i] >= ibox[i])
489 index = indexn(neig,ibox,nxyz);
490 range_check(index,0,len);
491 /* Compute the exponential factor */
493 efac = exp(-bE[j]+Emin);
496 /* Apply the bin volume correction for a multi-dimensional distance */
497 for(i=0; i<neig; i++) {
500 else if (idim[i] == 3)
501 efac /= sqr(eig[i][j]);
502 else if (idim[i] == -1)
503 efac /= sin(DEG2RAD*eig[i][j]);
505 /* Update the probability */
507 /* Update the energy */
509 E[index] += enerT[0][j];
510 /* Statistics: which "structure" in which bin */
515 /* Normalize probability */
516 normalize_p_e(len,P,nbin,E,pmin);
518 /* Compute boundaries for the Free energy */
522 /* Recompute Emin: it may have changed due to averaging */
525 for(i=0; (i<len); i++) {
527 Pmax = max(P[i],Pmax);
528 W[i] = -BOLTZ*Tref*log(P[i]);
533 Emin = min(E[i],Emin);
534 Emax = max(E[i],Emax);
535 Wmax = max(W[i],Wmax);
551 /* Write out the free energy as a function of bin index */
553 for(i=0; (i<len); i++)
556 S[i] = E[i]-W[i]-Smin;
557 fprintf(fp,"%5d %10.5e %10.5e %10.5e\n",i,W[i],E[i],S[i]);
565 /* Organize the structures in the bins */
567 snew(b->index,len+1);
570 for(i=0; (i<len); i++) {
571 b->index[i+1] = b->index[i]+nbin[i];
574 for(i=0; (i<n); i++) {
576 b->a[b->index[bi]+nbin[bi]] = i;
579 /* Consistency check */
580 /* This no longer applies when we allow the plot to be smaller
581 than the sampled space.
582 for(i=0; (i<len); i++) {
583 if (nbin[i] != (b->index[i+1] - b->index[i]))
584 gmx_fatal(FARGS,"nbin[%d] = %d, should be %d",i,nbin[i],
585 b->index[i+1] - b->index[i]);
588 /* Write the index file */
589 fp = ffopen(ndx,"w");
590 for(i=0; (i<len); i++) {
592 fprintf(fp,"[ %d ]\n",i);
593 for(j=b->index[i]; (j<b->index[i+1]); j++)
594 fprintf(fp,"%d\n",b->a[j]+1);
598 snew(axis_x,ibox[0]+1);
599 snew(axis_y,ibox[1]+1);
600 snew(axis_z,ibox[2]+1);
601 maxbox = max(ibox[0],max(ibox[1],ibox[2]));
602 snew(PP,maxbox*maxbox);
603 snew(WW,maxbox*maxbox);
604 snew(EE,maxbox*maxbox);
605 snew(SS,maxbox*maxbox);
606 for(i=0; (i<min(neig,3)); i++) {
608 case 0: axis = axis_x; break;
609 case 1: axis = axis_y; break;
610 case 2: axis = axis_z; break;
613 for(j=0; j<=ibox[i]; j++)
614 axis[j] = min_eig[i] + j/bfac[i];
618 snew(MM,maxbox*maxbox);
619 for(i=0; (i<ibox[0]); i++)
620 MM[i] = &(M[i*ibox[1]]);
623 for(i=0; (i<nmap); i++) {
624 Mmin = min(Mmin,map[0][i]);
625 Mmax = max(Mmax,map[0][i]);
628 for(i=0; (i<len); i++)
630 for(i=0; (i<nmap); i++) {
631 index = gmx_nint(mapindex[i]);
633 gmx_fatal(FARGS,"Number of bins in file from -mdata option does not correspond to current analysis");
636 M[index] = map[0][i];
643 pick_minima(logf,ibox,neig,len,W);
646 flags = MAT_SPATIAL_X | MAT_SPATIAL_Y;
648 /* Dump to XPM file */
650 for(i=0; (i<ibox[0]); i++) {
652 for(j=0; j<ibox[1]; j++) {
653 PP[i][j] = P[i*ibox[1]+j];
655 WW[i] = &(W[i*ibox[1]]);
656 EE[i] = &(E[i*ibox[1]]);
657 SS[i] = &(S[i*ibox[1]]);
659 fp = ffopen(xpmP,"w");
660 write_xpm(fp,flags,"Probability Distribution","","PC1","PC2",
661 ibox[0],ibox[1],axis_x,axis_y,PP,0,Pmax,rlo,rhi,&nlevels);
663 fp = ffopen(xpm,"w");
664 write_xpm(fp,flags,"Gibbs Energy Landscape","G (kJ/mol)","PC1","PC2",
665 ibox[0],ibox[1],axis_x,axis_y,WW,0,gmax,rlo,rhi,&nlevels);
667 fp = ffopen(xpm2,"w");
668 write_xpm(fp,flags,"Enthalpy Landscape","H (kJ/mol)","PC1","PC2",
669 ibox[0],ibox[1],axis_x,axis_y,EE,
670 emin ? *emin : Emin,emax ? *emax : Einf,rlo,rhi,&nlevels);
672 fp = ffopen(xpm3,"w");
673 write_xpm(fp,flags,"Entropy Landscape","TDS (kJ/mol)","PC1","PC2",
674 ibox[0],ibox[1],axis_x,axis_y,SS,0,Sinf,rlo,rhi,&nlevels);
677 fp = ffopen(xpm4,"w");
678 write_xpm(fp,flags,"Custom Landscape",mname,"PC1","PC2",
679 ibox[0],ibox[1],axis_x,axis_y,MM,0,Minf,rlo,rhi,&nlevels);
683 else if (neig == 3) {
684 /* Dump to PDB file */
685 fp = ffopen(pdb,"w");
686 for(i=0; (i<ibox[0]); i++) {
687 xxx[XX] = 3*(i+0.5-ibox[0]/2);
688 for(j=0; (j<ibox[1]); j++) {
689 xxx[YY] = 3*(j+0.5-ibox[1]/2);
690 for(k=0; (k<ibox[2]); k++) {
691 xxx[ZZ] = 3*(k+0.5-ibox[2]/2);
692 index = index3(ibox,i,j,k);
694 fprintf(fp,"%-6s%5u %-4.4s%3.3s %4d %8.3f%8.3f%8.3f%6.2f%6.2f\n",
695 "ATOM",(index+1) %10000,"H","H",(index+1)%10000,
696 xxx[XX],xxx[YY],xxx[ZZ],1.0,W[index]);
701 write_xplor("out.xplor",W,ibox,min_eig,max_eig);
703 write_xplor("user.xplor",M,ibox,min_eig,max_eig);
704 nxyz[XX] = imin/(ibox[1]*ibox[2]);
705 nxyz[YY] = (imin-nxyz[XX]*ibox[1]*ibox[2])/ibox[2];
706 nxyz[ZZ] = imin % ibox[2];
707 for(i=0; (i<ibox[0]); i++) {
709 for(j=0; (j<ibox[1]); j++)
710 WW[i][j] = W[index3(ibox,i,j,nxyz[ZZ])];
712 snew(buf,strlen(xpm)+4);
713 sprintf(buf,"%s",xpm);
714 sprintf(&buf[strlen(xpm)-4],"12.xpm");
715 fp = ffopen(buf,"w");
716 write_xpm(fp,flags,"Gibbs Energy Landscape","W (kJ/mol)","PC1","PC2",
717 ibox[0],ibox[1],axis_x,axis_y,WW,0,gmax,rlo,rhi,&nlevels);
719 for(i=0; (i<ibox[0]); i++) {
720 for(j=0; (j<ibox[2]); j++)
721 WW[i][j] = W[index3(ibox,i,nxyz[YY],j)];
723 sprintf(&buf[strlen(xpm)-4],"13.xpm");
724 fp = ffopen(buf,"w");
725 write_xpm(fp,flags,"SHAM Energy Landscape","kJ/mol","PC1","PC3",
726 ibox[0],ibox[2],axis_x,axis_z,WW,0,gmax,rlo,rhi,&nlevels);
728 for(i=0; (i<ibox[1]); i++) {
729 for(j=0; (j<ibox[2]); j++)
730 WW[i][j] = W[index3(ibox,nxyz[XX],i,j)];
732 sprintf(&buf[strlen(xpm)-4],"23.xpm");
733 fp = ffopen(buf,"w");
734 write_xpm(fp,flags,"SHAM Energy Landscape","kJ/mol","PC2","PC3",
735 ibox[1],ibox[2],axis_y,axis_z,WW,0,gmax,rlo,rhi,&nlevels);
745 static void ehisto(const char *fh,int n,real **enerT, const output_env_t oenv)
748 int i,j,k,nbin,blength;
750 real *T,bmin,bmax,bwidth;
758 for(j=1; (j<n); j++) {
759 for(k=0; (k<nbin); k++) {
760 if (T[k] == enerT[1][j]) {
767 T[nbin] = enerT[1][j];
770 bmin = min(enerT[0][j],bmin);
771 bmax = max(enerT[0][j],bmax);
774 blength = (bmax - bmin)/bwidth + 2;
776 for(i=0; (i<nbin); i++) {
777 snew(histo[i],blength);
779 for(j=0; (j<n); j++) {
780 k = (enerT[0][j]-bmin)/bwidth;
781 histo[bindex[j]][k]++;
783 fp = xvgropen(fh,"Energy distribution","E (kJ/mol)","",oenv);
784 for(j=0; (j<blength); j++) {
785 fprintf(fp,"%8.3f",bmin+j*bwidth);
786 for(k=0; (k<nbin); k++) {
787 fprintf(fp," %6d",histo[k][j]);
794 int gmx_sham(int argc,char *argv[])
796 const char *desc[] = {
797 "[TT]g_sham[tt] makes multi-dimensional free-energy, enthalpy and entropy plots.",
798 "[TT]g_sham[tt] reads one or more [TT].xvg[tt] files and analyzes data sets.",
799 "The basic purpose of [TT]g_sham[tt] is to plot Gibbs free energy landscapes",
800 "(option [TT]-ls[tt])",
801 "by Bolzmann inverting multi-dimensional histograms (option [TT]-lp[tt]),",
803 "make enthalpy (option [TT]-lsh[tt]) and entropy (option [TT]-lss[tt])",
804 "plots. The histograms can be made for any quantities the user supplies.",
805 "A line in the input file may start with a time",
806 "(see option [TT]-time[tt]) and any number of [IT]y[it]-values may follow.",
807 "Multiple sets can also be",
808 "read when they are separated by & (option [TT]-n[tt]),",
809 "in this case only one [IT]y[it]-value is read from each line.",
810 "All lines starting with # and @ are skipped.",
812 "Option [TT]-ge[tt] can be used to supply a file with free energies",
813 "when the ensemble is not a Boltzmann ensemble, but needs to be biased",
814 "by this free energy. One free energy value is required for each",
815 "(multi-dimensional) data point in the [TT]-f[tt] input.",
817 "Option [TT]-ene[tt] can be used to supply a file with energies.",
818 "These energies are used as a weighting function in the single",
819 "histogram analysis method by Kumar et al. When temperatures",
820 "are supplied (as a second column in the file), an experimental",
821 "weighting scheme is applied. In addition the vales",
822 "are used for making enthalpy and entropy plots.",
824 "With option [TT]-dim[tt], dimensions can be gives for distances.",
825 "When a distance is 2- or 3-dimensional, the circumference or surface",
826 "sampled by two particles increases with increasing distance.",
827 "Depending on what one would like to show, one can choose to correct",
828 "the histogram and free-energy for this volume effect.",
829 "The probability is normalized by r and r^2 for dimensions of 2 and 3, ",
831 "A value of -1 is used to indicate an angle in degrees between two",
832 "vectors: a sin(angle) normalization will be applied.",
833 "[BB]Note[bb] that for angles between vectors the inner-product or cosine",
834 "is the natural quantity to use, as it will produce bins of the same",
837 static real tb=-1,te=-1,frac=0.5,filtlen=0,binwidth=0.1;
838 static gmx_bool bHaveT=TRUE,bDer=FALSE,bSubAv=TRUE,bAverCorr=FALSE,bXYdy=FALSE;
839 static gmx_bool bEESEF=FALSE,bEENLC=FALSE,bEeFitAc=FALSE,bPower=FALSE;
840 static gmx_bool bShamEner=TRUE,bSham=TRUE;
841 static real Tref=298.15,pmin=0,ttol=0,pmax=0,gmax=0,emin=0,emax=0;
842 static rvec nrdim = {1,1,1};
843 static rvec nrbox = {32,32,32};
844 static rvec xmin = {0,0,0}, xmax={1,1,1};
845 static int nsets_in=1,nb_min=4,resol=10,nlevels=25;
846 static const char *mname="";
848 { "-time", FALSE, etBOOL, {&bHaveT},
849 "Expect a time in the input" },
850 { "-b", FALSE, etREAL, {&tb},
851 "First time to read from set" },
852 { "-e", FALSE, etREAL, {&te},
853 "Last time to read from set" },
854 { "-ttol", FALSE, etREAL, {&ttol},
855 "Tolerance on time in appropriate units (usually ps)" },
856 { "-n", FALSE, etINT, {&nsets_in},
857 "Read this number of sets separated by lines containing only an ampersand" },
858 { "-d", FALSE, etBOOL, {&bDer},
859 "Use the derivative" },
860 { "-bw", FALSE, etREAL, {&binwidth},
861 "Binwidth for the distribution" },
862 { "-sham", FALSE, etBOOL, {&bSham},
863 "Turn off energy weighting even if energies are given" },
864 { "-tsham", FALSE, etREAL, {&Tref},
865 "Temperature for single histogram analysis" },
866 { "-pmin", FALSE, etREAL, {&pmin},
867 "Minimum probability. Anything lower than this will be set to zero" },
868 { "-dim", FALSE, etRVEC, {nrdim},
869 "Dimensions for distances, used for volume correction (max 3 values, dimensions > 3 will get the same value as the last)" },
870 { "-ngrid", FALSE, etRVEC, {nrbox},
871 "Number of bins for energy landscapes (max 3 values, dimensions > 3 will get the same value as the last)" },
872 { "-xmin", FALSE, etRVEC, {xmin},
873 "Minimum for the axes in energy landscape (see above for > 3 dimensions)" },
874 { "-xmax", FALSE, etRVEC, {xmax},
875 "Maximum for the axes in energy landscape (see above for > 3 dimensions)" },
876 { "-pmax", FALSE, etREAL, {&pmax},
877 "Maximum probability in output, default is calculate" },
878 { "-gmax", FALSE, etREAL, {&gmax},
879 "Maximum free energy in output, default is calculate" },
880 { "-emin", FALSE, etREAL, {&emin},
881 "Minimum enthalpy in output, default is calculate" },
882 { "-emax", FALSE, etREAL, {&emax},
883 "Maximum enthalpy in output, default is calculate" },
884 { "-nlevels", FALSE, etINT, {&nlevels},
885 "Number of levels for energy landscape" },
886 { "-mname", FALSE, etSTR, {&mname},
887 "Legend label for the custom landscape" },
889 #define NPA asize(pa)
892 int n,e_n,d_n,nlast,s,nset,e_nset,d_nset,i,j=0,*idim,*ibox;
893 real **val,**et_val,**dt_val,*t,*e_t,e_dt,d_dt,*d_t,dt,tot,error;
895 double *av,*sig,cum1,cum2,cum3,cum4,db;
896 const char *fn_ge,*fn_ene;
898 gmx_large_int_t num_grid_points;
901 { efXVG, "-f", "graph", ffREAD },
902 { efXVG, "-ge", "gibbs", ffOPTRD },
903 { efXVG, "-ene", "esham", ffOPTRD },
904 { efXVG, "-dist", "ener", ffOPTWR },
905 { efXVG, "-histo","edist", ffOPTWR },
906 { efNDX, "-bin", "bindex", ffOPTWR },
907 { efXPM, "-lp", "prob", ffOPTWR },
908 { efXPM, "-ls", "gibbs", ffOPTWR },
909 { efXPM, "-lsh", "enthalpy", ffOPTWR },
910 { efXPM, "-lss", "entropy", ffOPTWR },
911 { efXPM, "-map", "map", ffOPTWR },
912 { efPDB, "-ls3", "gibbs3", ffOPTWR },
913 { efXVG, "-mdata","mapdata", ffOPTRD },
914 { efLOG, "-g", "shamlog", ffOPTWR }
916 #define NFILE asize(fnm)
922 CopyRight(stderr,argv[0]);
923 parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_BE_NICE ,
924 NFILE,fnm,npargs,pa,asize(desc),desc,0,NULL,&oenv);
926 val=read_xvg_time(opt2fn("-f",NFILE,fnm),bHaveT,
927 opt2parg_bSet("-b",npargs,pa),tb-ttol,
928 opt2parg_bSet("-e",npargs,pa),te+ttol,
929 nsets_in,&nset,&n,&dt,&t);
930 printf("Read %d sets of %d points, dt = %g\n\n",nset,n,dt);
932 fn_ge = opt2fn_null("-ge",NFILE,fnm);
933 fn_ene = opt2fn_null("-ene",NFILE,fnm);
936 gmx_fatal(FARGS,"Can not do free energy and energy corrections at the same time");
938 if (fn_ge || fn_ene) {
939 et_val=read_xvg_time(fn_ge ? fn_ge : fn_ene,bHaveT,
940 opt2parg_bSet("-b",npargs,pa),tb-ttol,
941 opt2parg_bSet("-e",npargs,pa),te+ttol,
942 1,&e_nset,&e_n,&e_dt,&e_t);
945 gmx_fatal(FARGS,"Can only handle one free energy component in %s",
948 if (e_nset!=1 && e_nset!=2)
949 gmx_fatal(FARGS,"Can only handle one energy component or one energy and one T in %s",
953 gmx_fatal(FARGS,"Number of energies (%d) does not match number of entries (%d) in %s",e_n,n,opt2fn("-f",NFILE,fnm));
958 if (opt2fn_null("-mdata",NFILE,fnm) != NULL) {
959 dt_val=read_xvg_time(opt2fn("-mdata",NFILE,fnm),bHaveT,
961 nsets_in,&d_nset,&d_n,&d_dt,&d_t);
963 gmx_fatal(FARGS,"Can only handle one mapping data column in %s",
964 opt2fn("-mdata",NFILE,fnm));
969 if (fn_ene && et_val)
970 ehisto(opt2fn("-histo",NFILE,fnm),e_n,et_val,oenv);
972 snew(idim,max(3,nset));
973 snew(ibox,max(3,nset));
974 snew(rmin,max(3,nset));
975 snew(rmax,max(3,nset));
976 for(i=0; (i<min(3,nset)); i++) {
982 for(; (i<nset); i++) {
989 /* Check that the grid size is manageable. */
990 num_grid_points = ibox[0];
991 for(i = 1; i < nset; i++)
993 gmx_large_int_t result;
994 if (!check_int_multiply_for_overflow(num_grid_points, ibox[i], &result))
997 "The number of dimensions and grid points is too large for this tool\n"
998 "to handle with what it knows about the architecture upon which it\n"
999 "is running. Use a different machine or consult the GROMACS mailing list.");
1001 num_grid_points = result;
1003 /* The number of grid points fits in a gmx_large_int_t. */
1005 do_sham(opt2fn("-dist",NFILE,fnm),opt2fn("-bin",NFILE,fnm),
1006 opt2fn("-lp",NFILE,fnm),
1007 opt2fn("-ls",NFILE,fnm),opt2fn("-lsh",NFILE,fnm),
1008 opt2fn("-lss",NFILE,fnm),opt2fn("-map",NFILE,fnm),
1009 opt2fn("-ls3",NFILE,fnm),opt2fn("-g",NFILE,fnm),
1010 n,nset,val,fn_ge!=NULL,e_nset,et_val,d_n,d_t,dt_val,Tref,
1012 opt2parg_bSet("-emin",NPA,pa) ? &emin : NULL,
1013 opt2parg_bSet("-emax",NPA,pa) ? &emax : NULL,
1015 mname,bSham,idim,ibox,
1016 opt2parg_bSet("-xmin",NPA,pa),rmin,
1017 opt2parg_bSet("-xmax",NPA,pa),rmax);