static void do_dhdl(t_enxframe *fr, FILE **fp_dhdl, const char *filename,
+ int *blocks, int *hists, int *samples, int *nlambdas,
const output_env_t oenv)
{
const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda";
{
gmx_fatal(FARGS, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do.");
}
-
if (! *fp_dhdl)
{
if (nblock_dh>0)
+ (*hists)+=nblock_hist;
+ (*blocks)+=nblock_dh;
+ (*nlambdas) = nblock_hist+nblock_dh;
+
+
/* write the data */
if (nblock_hist > 0)
{
+ gmx_large_int_t sum=0;
/* histograms */
for(i=0;i<fr->nblock;i++)
{
xmax=(x0+k+1)*dx;
fprintf(*fp_dhdl,"%g %d\n%g %d\n", xmin, hist,
xmax, hist);
+ sum+=hist;
}
/* multiple histogram data blocks in one histogram
mean that the second one is the reverse of the first one:
}
}
}
+
+ (*samples) += (int)(sum/nblock_hist);
}
else
{
}
}
+
if (first)
{
xvgr_legend(*fp_dhdl, nblock_dh, (const char**)setnames, oenv);
sfree(setnames);
}
+ (*samples) += len;
for(i=0;i<len;i++)
{
double time=start_time + delta_time*i;
"Volume", "Pressure"
};
- FILE *out,*fp_pairs=NULL,*fort=NULL,*fodt=NULL,*foten=NULL;
+ FILE *out=NULL,*fp_pairs=NULL,*fort=NULL,*fodt=NULL,*foten=NULL;
FILE *fp_dhdl=NULL;
FILE **drout;
ener_file_t fp;
t_enxblock *blk=NULL;
t_enxblock *blk_disre=NULL;
int ndisre=0;
+ int dh_blocks=0, dh_hists=0, dh_samples=0, dh_lambdas;
t_filenm fnm[] = {
{ efEDR, "-f", NULL, ffREAD },
bVisco = opt2bSet("-vis",NFILE,fnm);
- if (!bDisRe) {
- if (bVisco) {
- nset=asize(setnm);
- snew(set,nset);
- /* This is nasty code... To extract Pres tensor, Volume and Temperature */
- for(j=0; j<nset; j++) {
- for(i=0; i<nre; i++) {
- if (strstr(enm[i].name,setnm[j])) {
- set[j]=i;
- break;
- }
- }
- if (i == nre) {
- if (gmx_strcasecmp(setnm[j],"Volume")==0) {
- printf("Enter the box volume (" unit_volume "): ");
- if(1 != scanf("%lf",&dbl))
- {
- gmx_fatal(FARGS,"Error reading user input");
- }
- Vaver = dbl;
- } else
- gmx_fatal(FARGS,"Could not find term %s for viscosity calculation",
- setnm[j]);
- }
+ if (!bDisRe && !bDHDL)
+ {
+ if (bVisco) {
+ nset=asize(setnm);
+ snew(set,nset);
+ /* This is nasty code... To extract Pres tensor, Volume and Temperature */
+ for(j=0; j<nset; j++) {
+ for(i=0; i<nre; i++) {
+ if (strstr(enm[i].name,setnm[j])) {
+ set[j]=i;
+ break;
+ }
+ }
+ if (i == nre) {
+ if (gmx_strcasecmp(setnm[j],"Volume")==0) {
+ printf("Enter the box volume (" unit_volume "): ");
+ if(1 != scanf("%lf",&dbl))
+ {
+ gmx_fatal(FARGS,"Error reading user input");
+ }
+ Vaver = dbl;
+ } else
+ gmx_fatal(FARGS,"Could not find term %s for viscosity calculation",
+ setnm[j]);
+ }
+ }
}
- }
- else
- {
- set=select_by_name(nre,enm,&nset);
- }
- /* Print all the different units once */
- sprintf(buf,"(%s)",enm[set[0]].unit);
- for(i=1; i<nset; i++) {
- for(j=0; j<i; j++) {
- if (strcmp(enm[set[i]].unit,enm[set[j]].unit) == 0) {
- break;
- }
+ else
+ {
+ set=select_by_name(nre,enm,&nset);
}
- if (j == i) {
- strcat(buf,", (");
- strcat(buf,enm[set[i]].unit);
- strcat(buf,")");
+ /* Print all the different units once */
+ sprintf(buf,"(%s)",enm[set[0]].unit);
+ for(i=1; i<nset; i++) {
+ for(j=0; j<i; j++) {
+ if (strcmp(enm[set[i]].unit,enm[set[j]].unit) == 0) {
+ break;
+ }
+ }
+ if (j == i) {
+ strcat(buf,", (");
+ strcat(buf,enm[set[i]].unit);
+ strcat(buf,")");
+ }
}
- }
- out=xvgropen(opt2fn("-o",NFILE,fnm),"Gromacs Energies","Time (ps)",buf,
- oenv);
-
- snew(leg,nset+1);
- for(i=0; (i<nset); i++)
- leg[i] = enm[set[i]].name;
- if (bSum) {
- leg[nset]=strdup("Sum");
- xvgr_legend(out,nset+1,(const char**)leg,oenv);
- }
- else
- xvgr_legend(out,nset,(const char**)leg,oenv);
+ out=xvgropen(opt2fn("-o",NFILE,fnm),"Gromacs Energies","Time (ps)",buf,
+ oenv);
+
+ snew(leg,nset+1);
+ for(i=0; (i<nset); i++)
+ leg[i] = enm[set[i]].name;
+ if (bSum) {
+ leg[nset]=strdup("Sum");
+ xvgr_legend(out,nset+1,(const char**)leg,oenv);
+ }
+ else
+ xvgr_legend(out,nset,(const char**)leg,oenv);
- snew(bIsEner,nset);
- for(i=0; (i<nset); i++) {
- bIsEner[i] = FALSE;
- for (j=0; (j <= F_ETOT); j++)
- bIsEner[i] = bIsEner[i] ||
- (gmx_strcasecmp(interaction_function[j].longname,leg[i]) == 0);
- }
-
- if (bPrAll && nset > 1) {
- gmx_fatal(FARGS,"Printing averages can only be done when a single set is selected");
- }
+ snew(bIsEner,nset);
+ for(i=0; (i<nset); i++) {
+ bIsEner[i] = FALSE;
+ for (j=0; (j <= F_ETOT); j++)
+ bIsEner[i] = bIsEner[i] ||
+ (gmx_strcasecmp(interaction_function[j].longname,leg[i]) == 0);
+ }
- time = NULL;
+ if (bPrAll && nset > 1) {
+ gmx_fatal(FARGS,"Printing averages can only be done when a single set is selected");
+ }
- if (bORIRE || bOTEN)
- get_orires_parms(ftp2fn(efTPX,NFILE,fnm),&nor,&nex,&or_label,&oobs);
-
- if (bORIRE) {
- if (bOrinst)
- enx_i = enxORI;
- else
- enx_i = enxOR;
-
- if (bORA || bODA)
- snew(orient,nor);
- if (bODR)
- snew(odrms,nor);
- if (bORT || bODT) {
- fprintf(stderr,"Select the orientation restraint labels you want (-1 is all)\n");
- fprintf(stderr,"End your selection with 0\n");
- j = -1;
- orsel = NULL;
- do {
- j++;
- srenew(orsel,j+1);
- if(1 != scanf("%d",&(orsel[j])))
- {
- gmx_fatal(FARGS,"Error reading user input");
- }
- } while (orsel[j] > 0);
- if (orsel[0] == -1) {
- fprintf(stderr,"Selecting all %d orientation restraints\n",nor);
- norsel = nor;
- srenew(orsel,nor);
- for(i=0; i<nor; i++)
- orsel[i] = i;
- } else {
- /* Build the selection */
- norsel=0;
- for(i=0; i<j; i++) {
- for(k=0; k<nor; k++)
- if (or_label[k] == orsel[i]) {
- orsel[norsel] = k;
- norsel++;
- break;
- }
- if (k == nor)
- fprintf(stderr,"Orientation restraint label %d not found\n",
- orsel[i]);
- }
- }
- snew(odtleg,norsel);
- for(i=0; i<norsel; i++) {
- snew(odtleg[i],256);
- sprintf(odtleg[i],"%d",or_label[orsel[i]]);
- }
- if (bORT) {
- fort=xvgropen(opt2fn("-ort",NFILE,fnm), "Calculated orientations",
- "Time (ps)","",oenv);
- if (bOrinst)
- fprintf(fort,"%s",orinst_sub);
- xvgr_legend(fort,norsel,(const char**)odtleg,oenv);
- }
- if (bODT) {
- fodt=xvgropen(opt2fn("-odt",NFILE,fnm),
- "Orientation restraint deviation",
- "Time (ps)","",oenv);
- if (bOrinst)
- fprintf(fodt,"%s",orinst_sub);
- xvgr_legend(fodt,norsel,(const char**)odtleg,oenv);
- }
+ time = NULL;
+
+ if (bORIRE || bOTEN)
+ get_orires_parms(ftp2fn(efTPX,NFILE,fnm),&nor,&nex,&or_label,&oobs);
+
+ if (bORIRE) {
+ if (bOrinst)
+ enx_i = enxORI;
+ else
+ enx_i = enxOR;
+
+ if (bORA || bODA)
+ snew(orient,nor);
+ if (bODR)
+ snew(odrms,nor);
+ if (bORT || bODT) {
+ fprintf(stderr,"Select the orientation restraint labels you want (-1 is all)\n");
+ fprintf(stderr,"End your selection with 0\n");
+ j = -1;
+ orsel = NULL;
+ do {
+ j++;
+ srenew(orsel,j+1);
+ if(1 != scanf("%d",&(orsel[j])))
+ {
+ gmx_fatal(FARGS,"Error reading user input");
+ }
+ } while (orsel[j] > 0);
+ if (orsel[0] == -1) {
+ fprintf(stderr,"Selecting all %d orientation restraints\n",nor);
+ norsel = nor;
+ srenew(orsel,nor);
+ for(i=0; i<nor; i++)
+ orsel[i] = i;
+ } else {
+ /* Build the selection */
+ norsel=0;
+ for(i=0; i<j; i++) {
+ for(k=0; k<nor; k++)
+ if (or_label[k] == orsel[i]) {
+ orsel[norsel] = k;
+ norsel++;
+ break;
+ }
+ if (k == nor)
+ fprintf(stderr,"Orientation restraint label %d not found\n",
+ orsel[i]);
+ }
+ }
+ snew(odtleg,norsel);
+ for(i=0; i<norsel; i++) {
+ snew(odtleg[i],256);
+ sprintf(odtleg[i],"%d",or_label[orsel[i]]);
+ }
+ if (bORT) {
+ fort=xvgropen(opt2fn("-ort",NFILE,fnm), "Calculated orientations",
+ "Time (ps)","",oenv);
+ if (bOrinst)
+ fprintf(fort,"%s",orinst_sub);
+ xvgr_legend(fort,norsel,(const char**)odtleg,oenv);
+ }
+ if (bODT) {
+ fodt=xvgropen(opt2fn("-odt",NFILE,fnm),
+ "Orientation restraint deviation",
+ "Time (ps)","",oenv);
+ if (bOrinst)
+ fprintf(fodt,"%s",orinst_sub);
+ xvgr_legend(fodt,norsel,(const char**)odtleg,oenv);
+ }
+ }
}
- }
- if (bOTEN) {
- foten=xvgropen(opt2fn("-oten",NFILE,fnm),
- "Order tensor","Time (ps)","",oenv);
- snew(otenleg,bOvec ? nex*12 : nex*3);
- for(i=0; i<nex; i++) {
- for(j=0; j<3; j++) {
- sprintf(buf,"eig%d",j+1);
- otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
- }
- if (bOvec) {
- for(j=0; j<9; j++) {
- sprintf(buf,"vec%d%s",j/3+1,j%3==0 ? "x" : (j%3==1 ? "y" : "z"));
- otenleg[12*i+3+j] = strdup(buf);
- }
- }
+ if (bOTEN) {
+ foten=xvgropen(opt2fn("-oten",NFILE,fnm),
+ "Order tensor","Time (ps)","",oenv);
+ snew(otenleg,bOvec ? nex*12 : nex*3);
+ for(i=0; i<nex; i++) {
+ for(j=0; j<3; j++) {
+ sprintf(buf,"eig%d",j+1);
+ otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
+ }
+ if (bOvec) {
+ for(j=0; j<9; j++) {
+ sprintf(buf,"vec%d%s",j/3+1,j%3==0 ? "x" : (j%3==1 ? "y" : "z"));
+ otenleg[12*i+3+j] = strdup(buf);
+ }
+ }
+ }
+ xvgr_legend(foten,bOvec ? nex*12 : nex*3,(const char**)otenleg,oenv);
}
- xvgr_legend(foten,bOvec ? nex*12 : nex*3,(const char**)otenleg,oenv);
- }
}
- else {
- nbounds=get_bounds(ftp2fn(efTPX,NFILE,fnm),&bounds,&index,&pair,&npairs,
- &mtop,&top,&ir);
- snew(violaver,npairs);
- out=xvgropen(opt2fn("-o",NFILE,fnm),"Sum of Violations",
- "Time (ps)","nm",oenv);
- xvgr_legend(out,2,drleg,oenv);
- if (bDRAll) {
- fp_pairs=xvgropen(opt2fn("-pairs",NFILE,fnm),"Pair Distances",
- "Time (ps)","Distance (nm)",oenv);
- if (output_env_get_print_xvgr_codes(oenv))
- fprintf(fp_pairs,"@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
- ir.dr_tau);
- }
+ else if (bDisRe)
+ {
+ nbounds=get_bounds(ftp2fn(efTPX,NFILE,fnm),&bounds,&index,&pair,&npairs,
+ &mtop,&top,&ir);
+ snew(violaver,npairs);
+ out=xvgropen(opt2fn("-o",NFILE,fnm),"Sum of Violations",
+ "Time (ps)","nm",oenv);
+ xvgr_legend(out,2,drleg,oenv);
+ if (bDRAll) {
+ fp_pairs=xvgropen(opt2fn("-pairs",NFILE,fnm),"Pair Distances",
+ "Time (ps)","Distance (nm)",oenv);
+ if (output_env_get_print_xvgr_codes(oenv))
+ fprintf(fp_pairs,"@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
+ ir.dr_tau);
+ }
}
/*
* Store energies for analysis afterwards...
*/
- if (!bDisRe && (fr->nre > 0)) {
+ if (!bDisRe && !bDHDL && (fr->nre > 0)) {
if (edat.nframes % 1000 == 0) {
srenew(time,edat.nframes+1000);
}
teller_disre++;
}
}
+ else if (bDHDL)
+ {
+ do_dhdl(fr, &fp_dhdl, opt2fn("-odh",NFILE,fnm),
+ &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas,
+ oenv);
+ }
/*******************************************
* E N E R G I E S
*******************************************/
fprintf(foten," %g",vals[i*12+j]);
fprintf(foten,"\n");
}
- if (bDHDL)
- {
- do_dhdl(fr, &fp_dhdl, opt2fn("-odh",NFILE,fnm), oenv);
- }
}
}
}
fprintf(stderr,"\n");
close_enx(fp);
-
- ffclose(out);
+ if (out)
+ ffclose(out);
- if (bDHDL)
- ffclose(fp_dhdl);
if (bDRAll)
- ffclose(fp_pairs);
+ ffclose(fp_pairs);
if (bORT)
- ffclose(fort);
+ ffclose(fort);
if (bODT)
- ffclose(fodt);
- if (bORA) {
- out = xvgropen(opt2fn("-ora",NFILE,fnm),
- "Average calculated orientations",
- "Restraint label","",oenv);
- if (bOrinst)
- fprintf(out,"%s",orinst_sub);
- for(i=0; i<nor; i++)
- fprintf(out,"%5d %g\n",or_label[i],orient[i]/norfr);
- ffclose(out);
+ ffclose(fodt);
+ if (bORA)
+ {
+ out = xvgropen(opt2fn("-ora",NFILE,fnm),
+ "Average calculated orientations",
+ "Restraint label","",oenv);
+ if (bOrinst)
+ fprintf(out,"%s",orinst_sub);
+ for(i=0; i<nor; i++)
+ fprintf(out,"%5d %g\n",or_label[i],orient[i]/norfr);
+ ffclose(out);
}
if (bODA) {
- out = xvgropen(opt2fn("-oda",NFILE,fnm),
- "Average restraint deviation",
- "Restraint label","",oenv);
- if (bOrinst)
- fprintf(out,"%s",orinst_sub);
- for(i=0; i<nor; i++)
- fprintf(out,"%5d %g\n",or_label[i],orient[i]/norfr-oobs[i]);
- ffclose(out);
+ out = xvgropen(opt2fn("-oda",NFILE,fnm),
+ "Average restraint deviation",
+ "Restraint label","",oenv);
+ if (bOrinst)
+ fprintf(out,"%s",orinst_sub);
+ for(i=0; i<nor; i++)
+ fprintf(out,"%5d %g\n",or_label[i],orient[i]/norfr-oobs[i]);
+ ffclose(out);
}
if (bODR) {
- out = xvgropen(opt2fn("-odr",NFILE,fnm),
- "RMS orientation restraint deviations",
- "Restraint label","",oenv);
- if (bOrinst)
- fprintf(out,"%s",orinst_sub);
- for(i=0; i<nor; i++)
- fprintf(out,"%5d %g\n",or_label[i],sqrt(odrms[i]/norfr));
- ffclose(out);
+ out = xvgropen(opt2fn("-odr",NFILE,fnm),
+ "RMS orientation restraint deviations",
+ "Restraint label","",oenv);
+ if (bOrinst)
+ fprintf(out,"%s",orinst_sub);
+ for(i=0; i<nor; i++)
+ fprintf(out,"%5d %g\n",or_label[i],sqrt(odrms[i]/norfr));
+ ffclose(out);
}
if (bOTEN)
- ffclose(foten);
-
- if (bDisRe) {
- analyse_disre(opt2fn("-viol",NFILE,fnm),
- teller_disre,violaver,bounds,index,pair,nbounds,oenv);
- } else {
- analyse_ener(opt2bSet("-corr",NFILE,fnm),opt2fn("-corr",NFILE,fnm),
- bFee,bSum,bFluct,opt2parg_bSet("-nmol",npargs,ppa),
- bVisco,opt2fn("-vis",NFILE,fnm),
- nmol,nconstr,start_step,start_t,frame[cur].step,frame[cur].t,
- time,reftemp,&edat,
- nset,set,bIsEner,leg,enm,Vaver,ezero,nbmin,nbmax,
- oenv);
+ ffclose(foten);
+
+ if (bDisRe)
+ {
+ analyse_disre(opt2fn("-viol",NFILE,fnm),
+ teller_disre,violaver,bounds,index,pair,nbounds,oenv);
+ }
+ else if (bDHDL)
+ {
+ if (fp_dhdl)
+ {
+ ffclose(fp_dhdl);
+ }
+
+ printf("\nWritten %d lambda values with %d samples as ",
+ dh_lambdas, dh_samples);
+ if (dh_hists > 0)
+ {
+ printf("%d dH histograms ", dh_hists);
+ }
+ if (dh_blocks> 0)
+ {
+ printf("%d dH data blocks ", dh_blocks);
+ }
+ printf("to %s\n", opt2fn("-odh",NFILE,fnm));
+ }
+ else
+ {
+ analyse_ener(opt2bSet("-corr",NFILE,fnm),opt2fn("-corr",NFILE,fnm),
+ bFee,bSum,bFluct,opt2parg_bSet("-nmol",npargs,ppa),
+ bVisco,opt2fn("-vis",NFILE,fnm),
+ nmol,nconstr,start_step,start_t,frame[cur].step,frame[cur].t,
+ time,reftemp,&edat,
+ nset,set,bIsEner,leg,enm,Vaver,ezero,nbmin,nbmax,
+ oenv);
}
if (opt2bSet("-f2",NFILE,fnm)) {
- fec(opt2fn("-f2",NFILE,fnm), opt2fn("-ravg",NFILE,fnm),
- reftemp, nset, set, leg, &edat, time ,oenv);
+ fec(opt2fn("-f2",NFILE,fnm), opt2fn("-ravg",NFILE,fnm),
+ reftemp, nset, set, leg, &edat, time ,oenv);
}
-
+
{
- const char *nxy = "-nxy";
-
- do_view(oenv,opt2fn("-o",NFILE,fnm),nxy);
- do_view(oenv,opt2fn_null("-ravg",NFILE,fnm),nxy);
- do_view(oenv,opt2fn_null("-ora",NFILE,fnm),nxy);
- do_view(oenv,opt2fn_null("-ort",NFILE,fnm),nxy);
- do_view(oenv,opt2fn_null("-oda",NFILE,fnm),nxy);
- do_view(oenv,opt2fn_null("-odr",NFILE,fnm),nxy);
- do_view(oenv,opt2fn_null("-odt",NFILE,fnm),nxy);
- do_view(oenv,opt2fn_null("-oten",NFILE,fnm),nxy);
- do_view(oenv,opt2fn_null("-odh",NFILE,fnm),nxy);
+ const char *nxy = "-nxy";
+
+ do_view(oenv,opt2fn("-o",NFILE,fnm),nxy);
+ do_view(oenv,opt2fn_null("-ravg",NFILE,fnm),nxy);
+ do_view(oenv,opt2fn_null("-ora",NFILE,fnm),nxy);
+ do_view(oenv,opt2fn_null("-ort",NFILE,fnm),nxy);
+ do_view(oenv,opt2fn_null("-oda",NFILE,fnm),nxy);
+ do_view(oenv,opt2fn_null("-odr",NFILE,fnm),nxy);
+ do_view(oenv,opt2fn_null("-odt",NFILE,fnm),nxy);
+ do_view(oenv,opt2fn_null("-oten",NFILE,fnm),nxy);
+ do_view(oenv,opt2fn_null("-odh",NFILE,fnm),nxy);
}
thanx(stderr);
-
+
return 0;
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff