#include "gmx_ana.h"
#include "mdebin.h"
-static real minthird=-1.0/3.0,minsixth=-1.0/6.0;
+static real minthird = -1.0/3.0, minsixth = -1.0/6.0;
typedef struct {
real sum;
} exactsum_t;
typedef struct {
- real *ener;
+ real *ener;
exactsum_t *es;
- gmx_bool bExactStat;
- double av;
- double rmsd;
- double ee;
- double slope;
+ gmx_bool bExactStat;
+ double av;
+ double rmsd;
+ double ee;
+ double slope;
} enerdat_t;
typedef struct {
- gmx_large_int_t nsteps;
- gmx_large_int_t npoints;
- int nframes;
+ gmx_large_int_t nsteps;
+ gmx_large_int_t npoints;
+ int nframes;
int *step;
int *steps;
int *points;
enerdat_t *s;
} enerdata_t;
-static double mypow(double x,double y)
+static double mypow(double x, double y)
{
- if (x > 0)
- return pow(x,y);
- else
- return 0.0;
+ if (x > 0)
+ {
+ return pow(x, y);
+ }
+ else
+ {
+ return 0.0;
+ }
}
-static int *select_it(int nre,char *nm[],int *nset)
+static int *select_it(int nre, char *nm[], int *nset)
{
- gmx_bool *bE;
- int n,k,j,i;
- int *set;
- gmx_bool bVerbose = TRUE;
+ gmx_bool *bE;
+ int n, k, j, i;
+ int *set;
+ gmx_bool bVerbose = TRUE;
- if ((getenv("VERBOSE")) != NULL)
- bVerbose = FALSE;
+ if ((getenv("VERBOSE")) != NULL)
+ {
+ bVerbose = FALSE;
+ }
- fprintf(stderr,"Select the terms you want from the following list\n");
- fprintf(stderr,"End your selection with 0\n");
+ fprintf(stderr, "Select the terms you want from the following list\n");
+ fprintf(stderr, "End your selection with 0\n");
- if ( bVerbose ) {
- for(k=0; (k<nre); ) {
- for(j=0; (j<4) && (k<nre); j++,k++)
- fprintf(stderr," %3d=%14s",k+1,nm[k]);
- fprintf(stderr,"\n");
+ if (bVerbose)
+ {
+ for (k = 0; (k < nre); )
+ {
+ for (j = 0; (j < 4) && (k < nre); j++, k++)
+ {
+ fprintf(stderr, " %3d=%14s", k+1, nm[k]);
+ }
+ fprintf(stderr, "\n");
+ }
}
- }
- snew(bE,nre);
- do {
- if(1 != scanf("%d",&n))
+ snew(bE, nre);
+ do
{
- gmx_fatal(FARGS,"Error reading user input");
+ if (1 != scanf("%d", &n))
+ {
+ gmx_fatal(FARGS, "Error reading user input");
+ }
+ if ((n > 0) && (n <= nre))
+ {
+ bE[n-1] = TRUE;
+ }
}
- if ((n>0) && (n<=nre))
- bE[n-1]=TRUE;
- } while (n != 0);
+ while (n != 0);
- snew(set,nre);
- for(i=(*nset)=0; (i<nre); i++)
- if (bE[i])
- set[(*nset)++]=i;
+ snew(set, nre);
+ for (i = (*nset) = 0; (i < nre); i++)
+ {
+ if (bE[i])
+ {
+ set[(*nset)++] = i;
+ }
+ }
- sfree(bE);
+ sfree(bE);
- return set;
+ return set;
}
static void chomp(char *buf)
{
- int len = strlen(buf);
+ int len = strlen(buf);
- while ((len > 0) && (buf[len-1] == '\n')) {
- buf[len-1] = '\0';
- len--;
- }
+ while ((len > 0) && (buf[len-1] == '\n'))
+ {
+ buf[len-1] = '\0';
+ len--;
+ }
}
-static int *select_by_name(int nre,gmx_enxnm_t *nm,int *nset)
+static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset)
{
- gmx_bool *bE;
- int n,k,kk,j,i,nmatch,nind,nss;
- int *set;
- gmx_bool bEOF,bVerbose = TRUE,bLong=FALSE;
- char *ptr,buf[STRLEN];
- const char *fm4="%3d %-14s";
- const char *fm2="%3d %-34s";
- char **newnm=NULL;
-
- if ((getenv("VERBOSE")) != NULL)
- bVerbose = FALSE;
-
- fprintf(stderr,"\n");
- fprintf(stderr,"Select the terms you want from the following list by\n");
- fprintf(stderr,"selecting either (part of) the name or the number or a combination.\n");
- fprintf(stderr,"End your selection with an empty line or a zero.\n");
- fprintf(stderr,"-------------------------------------------------------------------\n");
-
- snew(newnm,nre);
- j = 0;
- for(k=0; k<nre; k++) {
- newnm[k] = strdup(nm[k].name);
- /* Insert dashes in all the names */
- while ((ptr = strchr(newnm[k],' ')) != NULL) {
- *ptr = '-';
- }
- if ( bVerbose ) {
- if (j == 0) {
- if (k > 0) {
- fprintf(stderr,"\n");
- }
- bLong = FALSE;
- for(kk=k; kk<k+4; kk++) {
- if (kk < nre && strlen(nm[kk].name) > 14) {
- bLong = TRUE;
- }
- }
- } else {
- fprintf(stderr," ");
- }
- if (!bLong) {
- fprintf(stderr,fm4,k+1,newnm[k]);
- j++;
- if (j == 4) {
- j = 0;
- }
- } else {
- fprintf(stderr,fm2,k+1,newnm[k]);
- j++;
- if (j == 2) {
- j = 0;
- }
- }
- }
- }
- if ( bVerbose ) {
- fprintf(stderr,"\n\n");
- }
-
- snew(bE,nre);
-
- bEOF = FALSE;
- while (!bEOF && (fgets2(buf,STRLEN-1,stdin))) {
- /* Remove newlines */
- chomp(buf);
-
- /* Remove spaces */
- trim(buf);
-
- /* Empty line means end of input */
- bEOF = (strlen(buf) == 0);
- if (!bEOF) {
- ptr = buf;
- do {
- if (!bEOF) {
- /* First try to read an integer */
- nss = sscanf(ptr,"%d",&nind);
- if (nss == 1) {
- /* Zero means end of input */
- if (nind == 0) {
- bEOF = TRUE;
- } else if ((1<=nind) && (nind<=nre)) {
- bE[nind-1] = TRUE;
- } else {
- fprintf(stderr,"number %d is out of range\n",nind);
- }
- }
- else {
- /* Now try to read a string */
- i = strlen(ptr);
- nmatch = 0;
- for(nind=0; nind<nre; nind++) {
- if (gmx_strcasecmp(newnm[nind],ptr) == 0) {
- bE[nind] = TRUE;
- nmatch++;
- }
- }
- if (nmatch == 0) {
- i = strlen(ptr);
- nmatch = 0;
- for(nind=0; nind<nre; nind++) {
- if (gmx_strncasecmp(newnm[nind],ptr,i) == 0) {
- bE[nind] = TRUE;
- nmatch++;
- }
- }
- if (nmatch == 0) {
- fprintf(stderr,"String '%s' does not match anything\n",ptr);
- }
- }
- }
- }
- /* Look for the first space, and remove spaces from there */
- if ((ptr = strchr(ptr,' ')) != NULL) {
- trim(ptr);
- }
- } while (!bEOF && (ptr && (strlen(ptr) > 0)));
- }
- }
-
- snew(set,nre);
- for(i=(*nset)=0; (i<nre); i++)
- if (bE[i])
- set[(*nset)++]=i;
-
- sfree(bE);
-
- if (*nset == 0)
- gmx_fatal(FARGS,"No energy terms selected");
-
- for(i=0; (i<nre); i++)
- sfree(newnm[i]);
- sfree(newnm);
-
- return set;
+ gmx_bool *bE;
+ int n, k, kk, j, i, nmatch, nind, nss;
+ int *set;
+ gmx_bool bEOF, bVerbose = TRUE, bLong = FALSE;
+ char *ptr, buf[STRLEN];
+ const char *fm4 = "%3d %-14s";
+ const char *fm2 = "%3d %-34s";
+ char **newnm = NULL;
+
+ if ((getenv("VERBOSE")) != NULL)
+ {
+ bVerbose = FALSE;
+ }
+
+ fprintf(stderr, "\n");
+ fprintf(stderr, "Select the terms you want from the following list by\n");
+ fprintf(stderr, "selecting either (part of) the name or the number or a combination.\n");
+ fprintf(stderr, "End your selection with an empty line or a zero.\n");
+ fprintf(stderr, "-------------------------------------------------------------------\n");
+
+ snew(newnm, nre);
+ j = 0;
+ for (k = 0; k < nre; k++)
+ {
+ newnm[k] = strdup(nm[k].name);
+ /* Insert dashes in all the names */
+ while ((ptr = strchr(newnm[k], ' ')) != NULL)
+ {
+ *ptr = '-';
+ }
+ if (bVerbose)
+ {
+ if (j == 0)
+ {
+ if (k > 0)
+ {
+ fprintf(stderr, "\n");
+ }
+ bLong = FALSE;
+ for (kk = k; kk < k+4; kk++)
+ {
+ if (kk < nre && strlen(nm[kk].name) > 14)
+ {
+ bLong = TRUE;
+ }
+ }
+ }
+ else
+ {
+ fprintf(stderr, " ");
+ }
+ if (!bLong)
+ {
+ fprintf(stderr, fm4, k+1, newnm[k]);
+ j++;
+ if (j == 4)
+ {
+ j = 0;
+ }
+ }
+ else
+ {
+ fprintf(stderr, fm2, k+1, newnm[k]);
+ j++;
+ if (j == 2)
+ {
+ j = 0;
+ }
+ }
+ }
+ }
+ if (bVerbose)
+ {
+ fprintf(stderr, "\n\n");
+ }
+
+ snew(bE, nre);
+
+ bEOF = FALSE;
+ while (!bEOF && (fgets2(buf, STRLEN-1, stdin)))
+ {
+ /* Remove newlines */
+ chomp(buf);
+
+ /* Remove spaces */
+ trim(buf);
+
+ /* Empty line means end of input */
+ bEOF = (strlen(buf) == 0);
+ if (!bEOF)
+ {
+ ptr = buf;
+ do
+ {
+ if (!bEOF)
+ {
+ /* First try to read an integer */
+ nss = sscanf(ptr, "%d", &nind);
+ if (nss == 1)
+ {
+ /* Zero means end of input */
+ if (nind == 0)
+ {
+ bEOF = TRUE;
+ }
+ else if ((1 <= nind) && (nind <= nre))
+ {
+ bE[nind-1] = TRUE;
+ }
+ else
+ {
+ fprintf(stderr, "number %d is out of range\n", nind);
+ }
+ }
+ else
+ {
+ /* Now try to read a string */
+ i = strlen(ptr);
+ nmatch = 0;
+ for (nind = 0; nind < nre; nind++)
+ {
+ if (gmx_strcasecmp(newnm[nind], ptr) == 0)
+ {
+ bE[nind] = TRUE;
+ nmatch++;
+ }
+ }
+ if (nmatch == 0)
+ {
+ i = strlen(ptr);
+ nmatch = 0;
+ for (nind = 0; nind < nre; nind++)
+ {
+ if (gmx_strncasecmp(newnm[nind], ptr, i) == 0)
+ {
+ bE[nind] = TRUE;
+ nmatch++;
+ }
+ }
+ if (nmatch == 0)
+ {
+ fprintf(stderr, "String '%s' does not match anything\n", ptr);
+ }
+ }
+ }
+ }
+ /* Look for the first space, and remove spaces from there */
+ if ((ptr = strchr(ptr, ' ')) != NULL)
+ {
+ trim(ptr);
+ }
+ }
+ while (!bEOF && (ptr && (strlen(ptr) > 0)));
+ }
+ }
+
+ snew(set, nre);
+ for (i = (*nset) = 0; (i < nre); i++)
+ {
+ if (bE[i])
+ {
+ set[(*nset)++] = i;
+ }
+ }
+
+ sfree(bE);
+
+ if (*nset == 0)
+ {
+ gmx_fatal(FARGS, "No energy terms selected");
+ }
+
+ for (i = 0; (i < nre); i++)
+ {
+ sfree(newnm[i]);
+ }
+ sfree(newnm);
+
+ return set;
}
static void get_dhdl_parms(const char *topnm, t_inputrec *ir)
{
- gmx_mtop_t mtop;
- int natoms;
+ gmx_mtop_t mtop;
+ int natoms;
t_iatom *iatom;
- matrix box;
+ matrix box;
/* all we need is the ir to be able to write the label */
- read_tpx(topnm,ir,box,&natoms,NULL,NULL,NULL,&mtop);
+ read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, &mtop);
}
static void get_orires_parms(const char *topnm,
- int *nor,int *nex,int **label,real **obs)
+ int *nor, int *nex, int **label, real **obs)
{
- gmx_mtop_t mtop;
- gmx_localtop_t *top;
- t_inputrec ir;
- t_iparams *ip;
- int natoms,i;
- t_iatom *iatom;
- int nb;
- matrix box;
-
- read_tpx(topnm,&ir,box,&natoms,NULL,NULL,NULL,&mtop);
- top = gmx_mtop_generate_local_top(&mtop,&ir);
-
- ip = top->idef.iparams;
- iatom = top->idef.il[F_ORIRES].iatoms;
-
- /* Count how many distance restraint there are... */
- nb = top->idef.il[F_ORIRES].nr;
- if (nb == 0)
- gmx_fatal(FARGS,"No orientation restraints in topology!\n");
-
- *nor = nb/3;
- *nex = 0;
- snew(*label,*nor);
- snew(*obs,*nor);
- for(i=0; i<nb; i+=3) {
- (*label)[i/3] = ip[iatom[i]].orires.label;
- (*obs)[i/3] = ip[iatom[i]].orires.obs;
- if (ip[iatom[i]].orires.ex >= *nex)
- *nex = ip[iatom[i]].orires.ex+1;
- }
- fprintf(stderr,"Found %d orientation restraints with %d experiments",
- *nor,*nex);
+ gmx_mtop_t mtop;
+ gmx_localtop_t *top;
+ t_inputrec ir;
+ t_iparams *ip;
+ int natoms, i;
+ t_iatom *iatom;
+ int nb;
+ matrix box;
+
+ read_tpx(topnm, &ir, box, &natoms, NULL, NULL, NULL, &mtop);
+ top = gmx_mtop_generate_local_top(&mtop, &ir);
+
+ ip = top->idef.iparams;
+ iatom = top->idef.il[F_ORIRES].iatoms;
+
+ /* Count how many distance restraint there are... */
+ nb = top->idef.il[F_ORIRES].nr;
+ if (nb == 0)
+ {
+ gmx_fatal(FARGS, "No orientation restraints in topology!\n");
+ }
+
+ *nor = nb/3;
+ *nex = 0;
+ snew(*label, *nor);
+ snew(*obs, *nor);
+ for (i = 0; i < nb; i += 3)
+ {
+ (*label)[i/3] = ip[iatom[i]].orires.label;
+ (*obs)[i/3] = ip[iatom[i]].orires.obs;
+ if (ip[iatom[i]].orires.ex >= *nex)
+ {
+ *nex = ip[iatom[i]].orires.ex+1;
+ }
+ }
+ fprintf(stderr, "Found %d orientation restraints with %d experiments",
+ *nor, *nex);
}
static int get_bounds(const char *topnm,
- real **bounds,int **index,int **dr_pair,int *npairs,
- gmx_mtop_t *mtop,gmx_localtop_t **ltop,t_inputrec *ir)
+ real **bounds, int **index, int **dr_pair, int *npairs,
+ gmx_mtop_t *mtop, gmx_localtop_t **ltop, t_inputrec *ir)
{
- gmx_localtop_t *top;
- t_functype *functype;
- t_iparams *ip;
- int natoms,i,j,k,type,ftype,natom;
- t_ilist *disres;
- t_iatom *iatom;
- real *b;
- int *ind,*pair;
- int nb,label1;
- matrix box;
-
- read_tpx(topnm,ir,box,&natoms,NULL,NULL,NULL,mtop);
- snew(*ltop,1);
- top = gmx_mtop_generate_local_top(mtop,ir);
- *ltop = top;
-
- functype = top->idef.functype;
- ip = top->idef.iparams;
-
- /* Count how many distance restraint there are... */
- nb=top->idef.il[F_DISRES].nr;
- if (nb == 0)
- gmx_fatal(FARGS,"No distance restraints in topology!\n");
-
- /* Allocate memory */
- snew(b,nb);
- snew(ind,nb);
- snew(pair,nb+1);
-
- /* Fill the bound array */
- nb=0;
- for(i=0; (i<top->idef.ntypes); i++) {
- ftype = functype[i];
- if (ftype == F_DISRES) {
-
- label1 = ip[i].disres.label;
- b[nb] = ip[i].disres.up1;
- ind[nb] = label1;
- nb++;
- }
- }
- *bounds = b;
-
- /* Fill the index array */
- label1 = -1;
- disres = &(top->idef.il[F_DISRES]);
- iatom = disres->iatoms;
- for(i=j=k=0; (i<disres->nr); ) {
- type = iatom[i];
- ftype = top->idef.functype[type];
- natom = interaction_function[ftype].nratoms+1;
- if (label1 != top->idef.iparams[type].disres.label) {
- pair[j] = k;
- label1 = top->idef.iparams[type].disres.label;
- j ++;
- }
- k++;
- i += natom;
- }
- pair[j] = k;
- *npairs = k;
- if (j != nb)
- gmx_incons("get_bounds for distance restraints");
-
- *index = ind;
- *dr_pair = pair;
-
- return nb;
+ gmx_localtop_t *top;
+ t_functype *functype;
+ t_iparams *ip;
+ int natoms, i, j, k, type, ftype, natom;
+ t_ilist *disres;
+ t_iatom *iatom;
+ real *b;
+ int *ind, *pair;
+ int nb, label1;
+ matrix box;
+
+ read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, mtop);
+ snew(*ltop, 1);
+ top = gmx_mtop_generate_local_top(mtop, ir);
+ *ltop = top;
+
+ functype = top->idef.functype;
+ ip = top->idef.iparams;
+
+ /* Count how many distance restraint there are... */
+ nb = top->idef.il[F_DISRES].nr;
+ if (nb == 0)
+ {
+ gmx_fatal(FARGS, "No distance restraints in topology!\n");
+ }
+
+ /* Allocate memory */
+ snew(b, nb);
+ snew(ind, nb);
+ snew(pair, nb+1);
+
+ /* Fill the bound array */
+ nb = 0;
+ for (i = 0; (i < top->idef.ntypes); i++)
+ {
+ ftype = functype[i];
+ if (ftype == F_DISRES)
+ {
+
+ label1 = ip[i].disres.label;
+ b[nb] = ip[i].disres.up1;
+ ind[nb] = label1;
+ nb++;
+ }
+ }
+ *bounds = b;
+
+ /* Fill the index array */
+ label1 = -1;
+ disres = &(top->idef.il[F_DISRES]);
+ iatom = disres->iatoms;
+ for (i = j = k = 0; (i < disres->nr); )
+ {
+ type = iatom[i];
+ ftype = top->idef.functype[type];
+ natom = interaction_function[ftype].nratoms+1;
+ if (label1 != top->idef.iparams[type].disres.label)
+ {
+ pair[j] = k;
+ label1 = top->idef.iparams[type].disres.label;
+ j++;
+ }
+ k++;
+ i += natom;
+ }
+ pair[j] = k;
+ *npairs = k;
+ if (j != nb)
+ {
+ gmx_incons("get_bounds for distance restraints");
+ }
+
+ *index = ind;
+ *dr_pair = pair;
+
+ return nb;
}
-static void calc_violations(real rt[],real rav3[],int nb,int index[],
- real bounds[],real *viol,double *st,double *sa)
+static void calc_violations(real rt[], real rav3[], int nb, int index[],
+ real bounds[], real *viol, double *st, double *sa)
{
- const real sixth=1.0/6.0;
- int i,j;
- double rsum,rav,sumaver,sumt;
-
- sumaver = 0;
- sumt = 0;
- for(i=0; (i<nb); i++) {
- rsum = 0.0;
- rav = 0.0;
- for(j=index[i]; (j<index[i+1]); j++) {
- if (viol)
- viol[j] += mypow(rt[j],-3.0);
- rav += sqr(rav3[j]);
- rsum += mypow(rt[j],-6);
- }
- rsum = max(0.0,mypow(rsum,-sixth)-bounds[i]);
- rav = max(0.0,mypow(rav, -sixth)-bounds[i]);
-
- sumt += rsum;
- sumaver += rav;
- }
- *st = sumt;
- *sa = sumaver;
+ const real sixth = 1.0/6.0;
+ int i, j;
+ double rsum, rav, sumaver, sumt;
+
+ sumaver = 0;
+ sumt = 0;
+ for (i = 0; (i < nb); i++)
+ {
+ rsum = 0.0;
+ rav = 0.0;
+ for (j = index[i]; (j < index[i+1]); j++)
+ {
+ if (viol)
+ {
+ viol[j] += mypow(rt[j], -3.0);
+ }
+ rav += sqr(rav3[j]);
+ rsum += mypow(rt[j], -6);
+ }
+ rsum = max(0.0, mypow(rsum, -sixth)-bounds[i]);
+ rav = max(0.0, mypow(rav, -sixth)-bounds[i]);
+
+ sumt += rsum;
+ sumaver += rav;
+ }
+ *st = sumt;
+ *sa = sumaver;
}
static void analyse_disre(const char *voutfn, int nframes,
- real violaver[], real bounds[], int index[],
- int pair[], int nbounds,
+ real violaver[], real bounds[], int index[],
+ int pair[], int nbounds,
const output_env_t oenv)
{
FILE *vout;
- double sum,sumt,sumaver;
- int i,j;
+ double sum, sumt, sumaver;
+ int i, j;
/* Subtract bounds from distances, to calculate violations */
- calc_violations(violaver,violaver,
- nbounds,pair,bounds,NULL,&sumt,&sumaver);
+ calc_violations(violaver, violaver,
+ nbounds, pair, bounds, NULL, &sumt, &sumaver);
#ifdef DEBUG
- fprintf(stdout,"\nSum of violations averaged over simulation: %g nm\n",
+ fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
sumaver);
- fprintf(stdout,"Largest violation averaged over simulation: %g nm\n\n",
+ fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n",
sumt);
#endif
- vout=xvgropen(voutfn,"r\\S-3\\N average violations","DR Index","nm",
- oenv);
+ vout = xvgropen(voutfn, "r\\S-3\\N average violations", "DR Index", "nm",
+ oenv);
sum = 0.0;
sumt = 0.0;
- for(i=0; (i<nbounds); i++) {
+ for (i = 0; (i < nbounds); i++)
+ {
/* Do ensemble averaging */
sumaver = 0;
- for(j=pair[i]; (j<pair[i+1]); j++)
+ for (j = pair[i]; (j < pair[i+1]); j++)
+ {
sumaver += sqr(violaver[j]/nframes);
- sumaver = max(0.0,mypow(sumaver,minsixth)-bounds[i]);
+ }
+ sumaver = max(0.0, mypow(sumaver, minsixth)-bounds[i]);
sumt += sumaver;
- sum = max(sum,sumaver);
- fprintf(vout,"%10d %10.5e\n",index[i],sumaver);
+ sum = max(sum, sumaver);
+ fprintf(vout, "%10d %10.5e\n", index[i], sumaver);
}
#ifdef DEBUG
- for(j=0; (j<dr.ndr); j++)
- fprintf(vout,"%10d %10.5e\n",j,mypow(violaver[j]/nframes,minthird));
+ for (j = 0; (j < dr.ndr); j++)
+ {
+ fprintf(vout, "%10d %10.5e\n", j, mypow(violaver[j]/nframes, minthird));
+ }
#endif
ffclose(vout);
- fprintf(stdout,"\nSum of violations averaged over simulation: %g nm\n",
+ fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
sumt);
- fprintf(stdout,"Largest violation averaged over simulation: %g nm\n\n",sum);
+ fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n", sum);
- do_view(oenv,voutfn,"-graphtype bar");
+ do_view(oenv, voutfn, "-graphtype bar");
}
-static void einstein_visco(const char *fn,const char *fni,int nsets,
- int nframes,real **sum,
- real V,real T,int nsteps,double time[],
+static void einstein_visco(const char *fn, const char *fni, int nsets,
+ int nframes, real **sum,
+ real V, real T, int nsteps, double time[],
const output_env_t oenv)
{
- FILE *fp0,*fp1;
- double av[4],avold[4];
- double fac,dt,di;
- int i,j,m,nf4;
+ FILE *fp0, *fp1;
+ double av[4], avold[4];
+ double fac, dt, di;
+ int i, j, m, nf4;
if (nframes < 1)
+ {
return;
+ }
dt = (time[1]-time[0]);
nf4 = nframes/4+1;
- for(i=0; i<=nsets; i++)
+ for (i = 0; i <= nsets; i++)
+ {
avold[i] = 0;
- fp0=xvgropen(fni,"Shear viscosity integral",
- "Time (ps)","(kg m\\S-1\\N s\\S-1\\N ps)",oenv);
- fp1=xvgropen(fn,"Shear viscosity using Einstein relation",
- "Time (ps)","(kg m\\S-1\\N s\\S-1\\N)",oenv);
- for(i=1; i<nf4; i++) {
+ }
+ fp0 = xvgropen(fni, "Shear viscosity integral",
+ "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv);
+ fp1 = xvgropen(fn, "Shear viscosity using Einstein relation",
+ "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv);
+ for (i = 1; i < nf4; i++)
+ {
fac = dt*nframes/nsteps;
- for(m=0; m<=nsets; m++)
+ for (m = 0; m <= nsets; m++)
+ {
av[m] = 0;
- for(j=0; j<nframes-i; j++) {
- for(m=0; m<nsets; m++) {
+ }
+ for (j = 0; j < nframes-i; j++)
+ {
+ for (m = 0; m < nsets; m++)
+ {
di = sqr(fac*(sum[m][j+i]-sum[m][j]));
av[m] += di;
}
/* Convert to SI for the viscosity */
fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nframes-i);
- fprintf(fp0,"%10g",time[i]-time[0]);
- for(m=0; (m<=nsets); m++) {
+ fprintf(fp0, "%10g", time[i]-time[0]);
+ for (m = 0; (m <= nsets); m++)
+ {
av[m] = fac*av[m];
- fprintf(fp0," %10g",av[m]);
+ fprintf(fp0, " %10g", av[m]);
}
- fprintf(fp0,"\n");
- fprintf(fp1,"%10g",0.5*(time[i]+time[i-1])-time[0]);
- for(m=0; (m<=nsets); m++) {
- fprintf(fp1," %10g",(av[m]-avold[m])/dt);
+ fprintf(fp0, "\n");
+ fprintf(fp1, "%10g", 0.5*(time[i]+time[i-1])-time[0]);
+ for (m = 0; (m <= nsets); m++)
+ {
+ fprintf(fp1, " %10g", (av[m]-avold[m])/dt);
avold[m] = av[m];
}
- fprintf(fp1,"\n");
+ fprintf(fp1, "\n");
}
ffclose(fp0);
ffclose(fp1);
typedef struct {
gmx_large_int_t np;
- double sum;
- double sav;
- double sav2;
+ double sum;
+ double sav;
+ double sav2;
} ee_sum_t;
typedef struct {
- int b;
- ee_sum_t sum;
+ int b;
+ ee_sum_t sum;
gmx_large_int_t nst;
gmx_large_int_t nst_min;
} ener_ee_t;
ees->sum = 0;
}
-static void add_ee_sum(ee_sum_t *ees,double sum,int np)
+static void add_ee_sum(ee_sum_t *ees, double sum, int np)
{
ees->np += np;
ees->sum += sum;
{
double av;
- av = ees->sum/ees->np;
+ av = ees->sum/ees->np;
ees->sav += av;
ees->sav2 += av*av;
ees->np = 0;
ees->sum = 0;
}
-static double calc_ee2(int nb,ee_sum_t *ees)
+static double calc_ee2(int nb, ee_sum_t *ees)
{
return (ees->sav2/nb - dsqr(ees->sav/nb))/(nb - 1);
}
if (debug)
{
char buf[STEPSTRSIZE];
- fprintf(debug,"Storing average for err.est.: %s steps\n",
- gmx_step_str(eee->nst,buf));
+ fprintf(debug, "Storing average for err.est.: %s steps\n",
+ gmx_step_str(eee->nst, buf));
}
add_ee_av(&eee->sum);
eee->b++;
eee->nst = 0;
}
-static void calc_averages(int nset,enerdata_t *edat,int nbmin,int nbmax)
+static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax)
{
- int nb,i,f,nee;
- double sum,sum2,sump,see2;
- gmx_large_int_t steps,np,p,bound_nb;
- enerdat_t *ed;
- exactsum_t *es;
- gmx_bool bAllZero;
- double x,sx,sy,sxx,sxy;
- ener_ee_t *eee;
+ int nb, i, f, nee;
+ double sum, sum2, sump, see2;
+ gmx_large_int_t steps, np, p, bound_nb;
+ enerdat_t *ed;
+ exactsum_t *es;
+ gmx_bool bAllZero;
+ double x, sx, sy, sxx, sxy;
+ ener_ee_t *eee;
/* Check if we have exact statistics over all points */
- for(i=0; i<nset; i++)
+ for (i = 0; i < nset; i++)
{
- ed = &edat->s[i];
+ ed = &edat->s[i];
ed->bExactStat = FALSE;
if (edat->npoints > 0)
{
* But if all energy values are 0, we still have exact sums.
*/
bAllZero = TRUE;
- for(f=0; f<edat->nframes && !ed->bExactStat; f++)
+ for (f = 0; f < edat->nframes && !ed->bExactStat; f++)
{
if (ed->ener[i] != 0)
{
}
}
- snew(eee,nbmax+1);
- for(i=0; i<nset; i++)
+ snew(eee, nbmax+1);
+ for (i = 0; i < nset; i++)
{
ed = &edat->s[i];
sy = 0;
sxx = 0;
sxy = 0;
- for(nb=nbmin; nb<=nbmax; nb++)
+ for (nb = nbmin; nb <= nbmax; nb++)
{
eee[nb].b = 0;
clear_ee_sum(&eee[nb].sum);
- eee[nb].nst = 0;
+ eee[nb].nst = 0;
eee[nb].nst_min = 0;
}
- for(f=0; f<edat->nframes; f++)
+ for (f = 0; f < edat->nframes; f++)
{
es = &ed->es[f];
sxx += p*x*x;
sxy += x*sump;
- for(nb=nbmin; nb<=nbmax; nb++)
+ for (nb = nbmin; nb <= nbmax; nb++)
{
/* Check if the current end step is closer to the desired
* block boundary than the next end step.
}
if (ed->bExactStat)
{
- add_ee_sum(&eee[nb].sum,es->sum,edat->points[f]);
+ add_ee_sum(&eee[nb].sum, es->sum, edat->points[f]);
}
else
{
- add_ee_sum(&eee[nb].sum,edat->s[i].ener[f],1);
+ add_ee_sum(&eee[nb].sum, edat->s[i].ener[f], 1);
}
bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
if (edat->step[f]*nb >= bound_nb)
nee = 0;
see2 = 0;
- for(nb=nbmin; nb<=nbmax; nb++)
+ for (nb = nbmin; nb <= nbmax; nb++)
{
/* Check if we actually got nb blocks and if the smallest
* block is not shorter than 80% of the average.
*/
if (debug)
{
- char buf1[STEPSTRSIZE],buf2[STEPSTRSIZE];
- fprintf(debug,"Requested %d blocks, we have %d blocks, min %s nsteps %s\n",
- nb,eee[nb].b,
- gmx_step_str(eee[nb].nst_min,buf1),
- gmx_step_str(edat->nsteps,buf2));
+ char buf1[STEPSTRSIZE], buf2[STEPSTRSIZE];
+ fprintf(debug, "Requested %d blocks, we have %d blocks, min %s nsteps %s\n",
+ nb, eee[nb].b,
+ gmx_step_str(eee[nb].nst_min, buf1),
+ gmx_step_str(edat->nsteps, buf2));
}
if (eee[nb].b == nb && 5*nb*eee[nb].nst_min >= 4*edat->nsteps)
{
- see2 += calc_ee2(nb,&eee[nb].sum);
+ see2 += calc_ee2(nb, &eee[nb].sum);
nee++;
}
}
sfree(eee);
}
-static enerdata_t *calc_sum(int nset,enerdata_t *edat,int nbmin,int nbmax)
+static enerdata_t *calc_sum(int nset, enerdata_t *edat, int nbmin, int nbmax)
{
enerdata_t *esum;
- enerdat_t *s;
- int f,i;
- double sum;
+ enerdat_t *s;
+ int f, i;
+ double sum;
- snew(esum,1);
+ snew(esum, 1);
*esum = *edat;
- snew(esum->s,1);
+ snew(esum->s, 1);
s = &esum->s[0];
- snew(s->ener,esum->nframes);
- snew(s->es ,esum->nframes);
+ snew(s->ener, esum->nframes);
+ snew(s->es, esum->nframes);
s->bExactStat = TRUE;
s->slope = 0;
- for(i=0; i<nset; i++)
+ for (i = 0; i < nset; i++)
{
if (!edat->s[i].bExactStat)
{
s->slope += edat->s[i].slope;
}
- for(f=0; f<edat->nframes; f++)
+ for (f = 0; f < edat->nframes; f++)
{
sum = 0;
- for(i=0; i<nset; i++)
+ for (i = 0; i < nset; i++)
{
sum += edat->s[i].ener[f];
}
s->ener[f] = sum;
- sum = 0;
- for(i=0; i<nset; i++)
+ sum = 0;
+ for (i = 0; i < nset; i++)
{
sum += edat->s[i].es[f].sum;
}
s->es[f].sum2 = 0;
}
- calc_averages(1,esum,nbmin,nbmax);
+ calc_averages(1, esum, nbmin, nbmax);
return esum;
}
-static char *ee_pr(double ee,char *buf)
+static char *ee_pr(double ee, char *buf)
{
char tmp[100];
double rnd;
if (ee < 0)
{
- sprintf(buf,"%s","--");
+ sprintf(buf, "%s", "--");
}
else
{
/* Round to two decimals by printing. */
- sprintf(tmp,"%.1e",ee);
- sscanf(tmp,"%lf",&rnd);
- sprintf(buf,"%g",rnd);
+ sprintf(tmp, "%.1e", ee);
+ sscanf(tmp, "%lf", &rnd);
+ sprintf(buf, "%g", rnd);
}
return buf;
}
-static void remove_drift(int nset,int nbmin,int nbmax,real dt,enerdata_t *edat)
+static void remove_drift(int nset, int nbmin, int nbmax, real dt, enerdata_t *edat)
{
/* Remove the drift by performing a fit to y = ax+b.
Uses 5 iterations. */
- int i,j,k;
- double delta,d,sd,sd2;
+ int i, j, k;
+ double delta, d, sd, sd2;
edat->npoints = edat->nframes;
- edat->nsteps = edat->nframes;
+ edat->nsteps = edat->nframes;
- for(k=0; (k<5); k++)
+ for (k = 0; (k < 5); k++)
{
- for(i=0; (i<nset); i++)
+ for (i = 0; (i < nset); i++)
{
delta = edat->s[i].slope*dt;
if (NULL != debug)
- fprintf(debug,"slope for set %d is %g\n",i,edat->s[i].slope);
+ {
+ fprintf(debug, "slope for set %d is %g\n", i, edat->s[i].slope);
+ }
- for(j=0; (j<edat->nframes); j++)
+ for (j = 0; (j < edat->nframes); j++)
{
edat->s[i].ener[j] -= j*delta;
edat->s[i].es[j].sum = 0;
edat->s[i].es[j].sum2 = 0;
}
}
- calc_averages(nset,edat,nbmin,nbmax);
+ calc_averages(nset, edat, nbmin, nbmax);
}
}
static void calc_fluctuation_props(FILE *fp,
- gmx_bool bDriftCorr,real dt,
- int nset,int set[],int nmol,
- char **leg,enerdata_t *edat,
- int nbmin,int nbmax)
+ gmx_bool bDriftCorr, real dt,
+ int nset, int set[], int nmol,
+ char **leg, enerdata_t *edat,
+ int nbmin, int nbmax)
{
- int i,j;
- double vvhh,vv,v,h,hh2,vv2,varv,hh,varh,tt,cv,cp,alpha,kappa,dcp,et,varet;
+ int i, j;
+ double vvhh, vv, v, h, hh2, vv2, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet;
double NANO3;
- enum { eVol, eEnth, eTemp, eEtot, eNR };
+ enum {
+ eVol, eEnth, eTemp, eEtot, eNR
+ };
const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
- int ii[eNR];
+ int ii[eNR];
NANO3 = NANO*NANO*NANO;
if (!bDriftCorr)
{
- fprintf(fp,"\nYou may want to use the -driftcorr flag in order to correct\n"
+ fprintf(fp, "\nYou may want to use the -driftcorr flag in order to correct\n"
"for spurious drift in the graphs. Note that this is not\n"
"a substitute for proper equilibration and sampling!\n");
}
else
{
- remove_drift(nset,nbmin,nbmax,dt,edat);
+ remove_drift(nset, nbmin, nbmax, dt, edat);
}
- for(i=0; (i<eNR); i++)
+ for (i = 0; (i < eNR); i++)
{
- for(ii[i]=0; (ii[i]<nset &&
- (gmx_strcasecmp(leg[ii[i]],my_ener[i]) != 0)); ii[i]++)
+ for (ii[i] = 0; (ii[i] < nset &&
+ (gmx_strcasecmp(leg[ii[i]], my_ener[i]) != 0)); ii[i]++)
+ {
;
+ }
/* if (ii[i] < nset)
fprintf(fp,"Found %s data.\n",my_ener[i]);
-*/ }
+ */ }
/* Compute it all! */
vvhh = alpha = kappa = cp = dcp = cv = NOTSET;
{
/* Only compute cv in constant volume runs, which we can test
by checking whether the enthalpy was computed.
- */
+ */
et = edat->s[ii[eEtot]].av;
varet = sqr(edat->s[ii[eEtot]].rmsd);
cv = KILO*((varet/nmol)/(BOLTZ*tt*tt));
if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
{
vvhh = 0;
- for(j=0; (j<edat->nframes); j++)
+ for (j = 0; (j < edat->nframes); j++)
{
- v = edat->s[ii[eVol]].ener[j]*NANO3;
- h = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO;
+ v = edat->s[ii[eVol]].ener[j]*NANO3;
+ h = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO;
vvhh += (v*h);
}
vvhh /= edat->nframes;
if (tt != NOTSET)
{
if (nmol < 2)
- fprintf(fp,"\nWARNING: nmol = %d, this may not be what you want.\n",
+ {
+ fprintf(fp, "\nWARNING: nmol = %d, this may not be what you want.\n",
nmol);
- fprintf(fp,"\nTemperature dependent fluctuation properties at T = %g.\n",tt);
- fprintf(fp,"\nHeat capacities obtained from fluctuations do *not* include\n");
- fprintf(fp,"quantum corrections. If you want to get a more accurate estimate\n");
- fprintf(fp,"please use the g_dos program.\n\n");
- fprintf(fp,"WARNING: Please verify that your simulations are converged and perform\n"
+ }
+ fprintf(fp, "\nTemperature dependent fluctuation properties at T = %g.\n", tt);
+ fprintf(fp, "\nHeat capacities obtained from fluctuations do *not* include\n");
+ fprintf(fp, "quantum corrections. If you want to get a more accurate estimate\n");
+ fprintf(fp, "please use the g_dos program.\n\n");
+ fprintf(fp, "WARNING: Please verify that your simulations are converged and perform\n"
"a block-averaging error analysis (not implemented in g_energy yet)\n");
if (debug != NULL)
{
if (varv != NOTSET)
- fprintf(fp,"varv = %10g (m^6)\n",varv*AVOGADRO/nmol);
+ {
+ fprintf(fp, "varv = %10g (m^6)\n", varv*AVOGADRO/nmol);
+ }
if (vvhh != NOTSET)
- fprintf(fp,"vvhh = %10g (m^3 J)\n",vvhh);
+ {
+ fprintf(fp, "vvhh = %10g (m^3 J)\n", vvhh);
+ }
}
if (vv != NOTSET)
- fprintf(fp,"Volume = %10g m^3/mol\n",
+ {
+ fprintf(fp, "Volume = %10g m^3/mol\n",
vv*AVOGADRO/nmol);
+ }
if (varh != NOTSET)
- fprintf(fp,"Enthalpy = %10g kJ/mol\n",
- hh*AVOGADRO/(KILO*nmol));
+ {
+ fprintf(fp, "Enthalpy = %10g kJ/mol\n",
+ hh*AVOGADRO/(KILO*nmol));
+ }
if (alpha != NOTSET)
- fprintf(fp,"Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n",
+ {
+ fprintf(fp, "Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n",
alpha);
+ }
if (kappa != NOTSET)
{
- fprintf(fp,"Isothermal Compressibility Kappa = %10g (J/m^3)\n",
+ fprintf(fp, "Isothermal Compressibility Kappa = %10g (J/m^3)\n",
kappa);
- fprintf(fp,"Adiabatic bulk modulus = %10g (m^3/J)\n",
+ fprintf(fp, "Adiabatic bulk modulus = %10g (m^3/J)\n",
1.0/kappa);
}
if (cp != NOTSET)
- fprintf(fp,"Heat capacity at constant pressure Cp = %10g J/mol K\n",
+ {
+ fprintf(fp, "Heat capacity at constant pressure Cp = %10g J/mol K\n",
cp);
+ }
if (cv != NOTSET)
- fprintf(fp,"Heat capacity at constant volume Cv = %10g J/mol K\n",
+ {
+ fprintf(fp, "Heat capacity at constant volume Cv = %10g J/mol K\n",
cv);
+ }
if (dcp != NOTSET)
- fprintf(fp,"Cp-Cv = %10g J/mol K\n",
+ {
+ fprintf(fp, "Cp-Cv = %10g J/mol K\n",
dcp);
- please_cite(fp,"Allen1987a");
+ }
+ please_cite(fp, "Allen1987a");
}
else
{
- fprintf(fp,"You should select the temperature in order to obtain fluctuation properties.\n");
+ fprintf(fp, "You should select the temperature in order to obtain fluctuation properties.\n");
}
}
-static void analyse_ener(gmx_bool bCorr,const char *corrfn,
- gmx_bool bFee,gmx_bool bSum,gmx_bool bFluct,gmx_bool bTempFluct,
- gmx_bool bVisco,const char *visfn,int nmol,
- gmx_large_int_t start_step,double start_t,
- gmx_large_int_t step,double t,
+static void analyse_ener(gmx_bool bCorr, const char *corrfn,
+ gmx_bool bFee, gmx_bool bSum, gmx_bool bFluct, gmx_bool bTempFluct,
+ gmx_bool bVisco, const char *visfn, int nmol,
+ gmx_large_int_t start_step, double start_t,
+ gmx_large_int_t step, double t,
double time[], real reftemp,
enerdata_t *edat,
- int nset,int set[],gmx_bool *bIsEner,
- char **leg,gmx_enxnm_t *enm,
- real Vaver,real ezero,
- int nbmin,int nbmax,
+ int nset, int set[], gmx_bool *bIsEner,
+ char **leg, gmx_enxnm_t *enm,
+ real Vaver, real ezero,
+ int nbmin, int nbmax,
const output_env_t oenv)
{
- FILE *fp;
- /* Check out the printed manual for equations! */
- double Dt,aver,stddev,errest,delta_t,totaldrift;
- enerdata_t *esum=NULL;
- real xxx,integral,intBulk,Temp=0,Pres=0;
- real sfrac,oldfrac,diffsum,diffav,fstep,pr_aver,pr_stddev,pr_errest;
- double beta=0,expE,expEtot,*fee=NULL;
- gmx_large_int_t nsteps;
- int nexact,nnotexact;
- double x1m,x1mk;
- int i,j,k,nout;
- real chi2;
- char buf[256],eebuf[100];
-
- nsteps = step - start_step + 1;
- if (nsteps < 1) {
- fprintf(stdout,"Not enough steps (%s) for statistics\n",
- gmx_step_str(nsteps,buf));
- }
- else {
- /* Calculate the time difference */
- delta_t = t - start_t;
-
- fprintf(stdout,"\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
- gmx_step_str(nsteps,buf),start_t,t,nset);
-
- calc_averages(nset,edat,nbmin,nbmax);
-
- if (bSum) {
- esum = calc_sum(nset,edat,nbmin,nbmax);
- }
-
- if (edat->npoints == 0) {
- nexact = 0;
- nnotexact = nset;
- } else {
- nexact = 0;
- nnotexact = 0;
- for(i=0; (i<nset); i++) {
- if (edat->s[i].bExactStat) {
- nexact++;
- } else {
- nnotexact++;
- }
- }
- }
-
- if (nnotexact == 0) {
- fprintf(stdout,"All statistics are over %s points\n",
- gmx_step_str(edat->npoints,buf));
- } else if (nexact == 0 || edat->npoints == edat->nframes) {
- fprintf(stdout,"All statistics are over %d points (frames)\n",
- edat->nframes);
- } else {
- fprintf(stdout,"The term%s",nnotexact==1 ? "" : "s");
- for(i=0; (i<nset); i++) {
- if (!edat->s[i].bExactStat) {
- fprintf(stdout," '%s'",leg[i]);
- }
- }
- fprintf(stdout," %s has statistics over %d points (frames)\n",
- nnotexact==1 ? "is" : "are",edat->nframes);
- fprintf(stdout,"All other statistics are over %s points\n",
- gmx_step_str(edat->npoints,buf));
- }
- fprintf(stdout,"\n");
-
- fprintf(stdout,"%-24s %10s %10s %10s %10s",
- "Energy","Average","Err.Est.","RMSD","Tot-Drift");
- if (bFee)
- fprintf(stdout," %10s\n","-kT ln<e^(E/kT)>");
- else
- fprintf(stdout,"\n");
- fprintf(stdout,"-------------------------------------------------------------------------------\n");
-
- /* Initiate locals, only used with -sum */
- expEtot=0;
- if (bFee) {
- beta = 1.0/(BOLTZ*reftemp);
- snew(fee,nset);
- }
- for(i=0; (i<nset); i++) {
- aver = edat->s[i].av;
- stddev = edat->s[i].rmsd;
- errest = edat->s[i].ee;
-
- if (bFee) {
- expE = 0;
- for(j=0; (j<edat->nframes); j++) {
- expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol);
- }
- if (bSum)
- expEtot+=expE/edat->nframes;
-
- fee[i] = log(expE/edat->nframes)/beta + aver/nmol;
- }
- if (strstr(leg[i],"empera") != NULL) {
- Temp = aver;
- } else if (strstr(leg[i],"olum") != NULL) {
- Vaver= aver;
- } else if (strstr(leg[i],"essure") != NULL) {
- Pres = aver;
- }
- if (bIsEner[i]) {
- pr_aver = aver/nmol-ezero;
- pr_stddev = stddev/nmol;
- pr_errest = errest/nmol;
- }
- else {
- pr_aver = aver;
- pr_stddev = stddev;
- pr_errest = errest;
- }
-
- /* Multiply the slope in steps with the number of steps taken */
- totaldrift = (edat->nsteps - 1)*edat->s[i].slope;
- if (bIsEner[i])
- {
- totaldrift /= nmol;
- }
-
- fprintf(stdout,"%-24s %10g %10s %10g %10g",
- leg[i],pr_aver,ee_pr(pr_errest,eebuf),pr_stddev,totaldrift);
- if (bFee)
- fprintf(stdout," %10g",fee[i]);
-
- fprintf(stdout," (%s)\n",enm[set[i]].unit);
-
- if (bFluct) {
- for(j=0; (j<edat->nframes); j++)
- edat->s[i].ener[j] -= aver;
- }
- }
- if (bSum) {
- totaldrift = (edat->nsteps - 1)*esum->s[0].slope;
- fprintf(stdout,"%-24s %10g %10s %10s %10g (%s)",
- "Total",esum->s[0].av/nmol,ee_pr(esum->s[0].ee/nmol,eebuf),
- "--",totaldrift/nmol,enm[set[0]].unit);
- /* pr_aver,pr_stddev,a,totaldrift */
- if (bFee)
- fprintf(stdout," %10g %10g\n",
- log(expEtot)/beta + esum->s[0].av/nmol,log(expEtot)/beta);
- else
- fprintf(stdout,"\n");
- }
-
- /* Do correlation function */
- if (edat->nframes > 1)
- {
- Dt = delta_t/(edat->nframes - 1);
+ FILE *fp;
+ /* Check out the printed manual for equations! */
+ double Dt, aver, stddev, errest, delta_t, totaldrift;
+ enerdata_t *esum = NULL;
+ real xxx, integral, intBulk, Temp = 0, Pres = 0;
+ real sfrac, oldfrac, diffsum, diffav, fstep, pr_aver, pr_stddev, pr_errest;
+ double beta = 0, expE, expEtot, *fee = NULL;
+ gmx_large_int_t nsteps;
+ int nexact, nnotexact;
+ double x1m, x1mk;
+ int i, j, k, nout;
+ real chi2;
+ char buf[256], eebuf[100];
+
+ nsteps = step - start_step + 1;
+ if (nsteps < 1)
+ {
+ fprintf(stdout, "Not enough steps (%s) for statistics\n",
+ gmx_step_str(nsteps, buf));
}
else
{
- Dt = 0;
- }
- if (bVisco) {
- const char* leg[] = { "Shear", "Bulk" };
- real factor;
- real **eneset;
- real **enesum;
-
- /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
-
- /* Symmetrise tensor! (and store in first three elements)
- * And subtract average pressure!
- */
- snew(eneset,12);
- for(i=0; i<12; i++) {
- snew(eneset[i],edat->nframes);
- }
- snew(enesum,3);
- for(i=0; i<3; i++) {
- snew(enesum[i],edat->nframes);
- }
- for(i=0; (i<edat->nframes); i++) {
- eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]);
- eneset[1][i] = 0.5*(edat->s[2].ener[i]+edat->s[6].ener[i]);
- eneset[2][i] = 0.5*(edat->s[5].ener[i]+edat->s[7].ener[i]);
- for(j=3; j<=11; j++) {
- eneset[j][i] = edat->s[j].ener[i];
- }
- eneset[11][i] -= Pres;
- enesum[0][i] = 0.5*(edat->s[1].es[i].sum+edat->s[3].es[i].sum);
- enesum[1][i] = 0.5*(edat->s[2].es[i].sum+edat->s[6].es[i].sum);
- enesum[2][i] = 0.5*(edat->s[5].es[i].sum+edat->s[7].es[i].sum);
- }
-
- einstein_visco("evisco.xvg","eviscoi.xvg",
- 3,edat->nframes,enesum,Vaver,Temp,nsteps,time,oenv);
-
- /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
- /* Do it for shear viscosity */
- strcpy(buf,"Shear Viscosity");
- low_do_autocorr(corrfn,oenv,buf,edat->nframes,3,
- (edat->nframes+1)/2,eneset,Dt,
- eacNormal,1,TRUE,FALSE,FALSE,0.0,0.0,0,1);
-
- /* Now for bulk viscosity */
- strcpy(buf,"Bulk Viscosity");
- low_do_autocorr(corrfn,oenv,buf,edat->nframes,1,
- (edat->nframes+1)/2,&(eneset[11]),Dt,
- eacNormal,1,TRUE,FALSE,FALSE,0.0,0.0,0,1);
-
- factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
- fp=xvgropen(visfn,buf,"Time (ps)","\\8h\\4 (cp)",oenv);
- xvgr_legend(fp,asize(leg),leg,oenv);
-
- /* Use trapezium rule for integration */
- integral = 0;
- intBulk = 0;
- nout = get_acfnout();
- if ((nout < 2) || (nout >= edat->nframes/2))
- nout = edat->nframes/2;
- for(i=1; (i<nout); i++)
- {
- integral += 0.5*(eneset[0][i-1] + eneset[0][i])*factor;
- intBulk += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
- fprintf(fp,"%10g %10g %10g\n",(i*Dt),integral,intBulk);
- }
- ffclose(fp);
- }
- else if (bCorr) {
- if (bFluct)
- strcpy(buf,"Autocorrelation of Energy Fluctuations");
- else
- strcpy(buf,"Energy Autocorrelation");
-#if 0
- do_autocorr(corrfn,oenv,buf,edat->nframes,
- bSum ? 1 : nset,
- bSum ? &edat->s[nset-1].ener : eneset,
- (delta_t/edat->nframes),eacNormal,FALSE);
-#endif
- }
- }
-}
-
-static void print_time(FILE *fp,double t)
-{
- fprintf(fp,"%12.6f",t);
-}
-
-static void print1(FILE *fp,gmx_bool bDp,real e)
-{
- if (bDp)
- fprintf(fp," %16.12f",e);
- else
- fprintf(fp," %10.6f",e);
-}
-
-static void fec(const char *ene2fn, const char *runavgfn,
- real reftemp, int nset, int set[], char *leg[],
- enerdata_t *edat, double time[],
- const output_env_t oenv)
-{
- const char* ravgleg[] = { "\\8D\\4E = E\\sB\\N-E\\sA\\N",
- "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N" };
- FILE *fp;
- ener_file_t enx;
- int nre,timecheck,step,nenergy,nenergy2,maxenergy;
- int i,j;
- gmx_bool bCont;
- real aver, beta;
- real **eneset2;
- double dE, sum;
- gmx_enxnm_t *enm=NULL;
- t_enxframe *fr;
- char buf[22];
-
- /* read second energy file */
- snew(fr,1);
- enm = NULL;
- enx = open_enx(ene2fn,"r");
- do_enxnms(enx,&(fr->nre),&enm);
-
- snew(eneset2,nset+1);
- nenergy2=0;
- maxenergy=0;
- timecheck=0;
- do {
- /* This loop searches for the first frame (when -b option is given),
- * or when this has been found it reads just one energy frame
- */
- do {
- bCont = do_enx(enx,fr);
-
- if (bCont)
- timecheck = check_times(fr->t);
-
- } while (bCont && (timecheck < 0));
-
- /* Store energies for analysis afterwards... */
- if ((timecheck == 0) && bCont) {
- if (fr->nre > 0) {
- if ( nenergy2 >= maxenergy ) {
- maxenergy += 1000;
- for(i=0; i<=nset; i++)
- srenew(eneset2[i],maxenergy);
- }
- if (fr->t != time[nenergy2])
- fprintf(stderr,"\nWARNING time mismatch %g!=%g at frame %s\n",
- fr->t, time[nenergy2], gmx_step_str(fr->step,buf));
- for(i=0; i<nset; i++)
- eneset2[i][nenergy2] = fr->ener[set[i]].e;
- nenergy2++;
- }
- }
- } while (bCont && (timecheck == 0));
-
- /* check */
- if (edat->nframes != nenergy2) {
- fprintf(stderr,"\nWARNING file length mismatch %d!=%d\n",
- edat->nframes,nenergy2);
- }
- nenergy = min(edat->nframes,nenergy2);
-
- /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
- fp=NULL;
- if (runavgfn) {
- fp=xvgropen(runavgfn,"Running average free energy difference",
- "Time (" unit_time ")","\\8D\\4E (" unit_energy ")",oenv);
- xvgr_legend(fp,asize(ravgleg),ravgleg,oenv);
- }
- fprintf(stdout,"\n%-24s %10s\n",
- "Energy","dF = -kT ln < exp(-(EB-EA)/kT) >A");
- sum=0;
- beta = 1.0/(BOLTZ*reftemp);
- for(i=0; i<nset; i++) {
- if (gmx_strcasecmp(leg[i],enm[set[i]].name)!=0)
- fprintf(stderr,"\nWARNING energy set name mismatch %s!=%s\n",
- leg[i],enm[set[i]].name);
- for(j=0; j<nenergy; j++) {
- dE = eneset2[i][j] - edat->s[i].ener[j];
- sum += exp(-dE*beta);
- if (fp)
- fprintf(fp,"%10g %10g %10g\n",
- time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) );
- }
- aver = -BOLTZ*reftemp*log(sum/nenergy);
- fprintf(stdout,"%-24s %10g\n",leg[i],aver);
- }
- if(fp) ffclose(fp);
- sfree(fr);
-}
+ /* Calculate the time difference */
+ delta_t = t - start_t;
+ fprintf(stdout, "\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
+ gmx_step_str(nsteps, buf), start_t, t, nset);
-static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl,
- const char *filename, gmx_bool bDp,
- int *blocks, int *hists, int *samples, int *nlambdas,
- const output_env_t oenv)
-{
- const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda";
- char title[STRLEN],label_x[STRLEN],label_y[STRLEN], legend[STRLEN];
- char buf[STRLEN];
- gmx_bool first=FALSE;
- int nblock_hist=0, nblock_dh=0, nblock_dhcoll=0;
- int i,j,k;
- /* coll data */
- double temp=0, start_time=0, delta_time=0, start_lambda=0, delta_lambda=0;
- static int setnr=0;
- double *native_lambda_vec=NULL;
- const char **lambda_components=NULL;
- int n_lambda_vec=0;
- gmx_bool changing_lambda=FALSE;
- int lambda_fep_state;
+ calc_averages(nset, edat, nbmin, nbmax);
- /* now count the blocks & handle the global dh data */
- for(i=0;i<fr->nblock;i++)
- {
- if (fr->block[i].id == enxDHHIST)
+ if (bSum)
{
- nblock_hist++;
+ esum = calc_sum(nset, edat, nbmin, nbmax);
}
- else if (fr->block[i].id == enxDH)
+
+ if (edat->npoints == 0)
{
- nblock_dh++;
+ nexact = 0;
+ nnotexact = nset;
}
- else if (fr->block[i].id == enxDHCOLL)
+ else
{
- nblock_dhcoll++;
- if ( (fr->block[i].nsub < 1) ||
- (fr->block[i].sub[0].type != xdr_datatype_double) ||
- (fr->block[i].sub[0].nr < 5))
- {
- gmx_fatal(FARGS, "Unexpected block data");
- }
-
- /* read the data from the DHCOLL block */
- temp = fr->block[i].sub[0].dval[0];
- start_time = fr->block[i].sub[0].dval[1];
- delta_time = fr->block[i].sub[0].dval[2];
- start_lambda = fr->block[i].sub[0].dval[3];
- delta_lambda = fr->block[i].sub[0].dval[4];
- changing_lambda = (delta_lambda != 0);
- if (fr->block[i].nsub > 1)
+ nexact = 0;
+ nnotexact = 0;
+ for (i = 0; (i < nset); i++)
{
- lambda_fep_state=fr->block[i].sub[1].ival[0];
- if (n_lambda_vec==0)
+ if (edat->s[i].bExactStat)
{
- n_lambda_vec=fr->block[i].sub[1].ival[1];
+ nexact++;
}
else
{
- if (n_lambda_vec!=fr->block[i].sub[1].ival[1])
- {
- gmx_fatal(FARGS,
- "Unexpected change of basis set in lambda");
- }
+ nnotexact++;
}
- if (lambda_components == NULL)
- snew(lambda_components, n_lambda_vec);
- if (native_lambda_vec == NULL)
- snew(native_lambda_vec, n_lambda_vec);
- for(j=0;j<n_lambda_vec;j++)
+ }
+ }
+
+ if (nnotexact == 0)
+ {
+ fprintf(stdout, "All statistics are over %s points\n",
+ gmx_step_str(edat->npoints, buf));
+ }
+ else if (nexact == 0 || edat->npoints == edat->nframes)
+ {
+ fprintf(stdout, "All statistics are over %d points (frames)\n",
+ edat->nframes);
+ }
+ else
+ {
+ fprintf(stdout, "The term%s", nnotexact == 1 ? "" : "s");
+ for (i = 0; (i < nset); i++)
+ {
+ if (!edat->s[i].bExactStat)
{
- native_lambda_vec[j] = fr->block[i].sub[0].dval[5+j];
- lambda_components[j] =
- efpt_singular_names[fr->block[i].sub[1].ival[2+j]];
+ fprintf(stdout, " '%s'", leg[i]);
}
}
+ fprintf(stdout, " %s has statistics over %d points (frames)\n",
+ nnotexact == 1 ? "is" : "are", edat->nframes);
+ fprintf(stdout, "All other statistics are over %s points\n",
+ gmx_step_str(edat->npoints, buf));
}
- }
+ fprintf(stdout, "\n");
- if (nblock_hist == 0 && nblock_dh == 0)
- {
- /* don't do anything */
- return;
- }
- if (nblock_hist>0 && nblock_dh>0)
- {
- 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)
+ fprintf(stdout, "%-24s %10s %10s %10s %10s",
+ "Energy", "Average", "Err.Est.", "RMSD", "Tot-Drift");
+ if (bFee)
{
- /* we have standard, non-histogram data --
- call open_dhdl to open the file */
- /* TODO this is an ugly hack that needs to be fixed: this will only
- work if the order of data is always the same and if we're
- only using the g_energy compiled with the mdrun that produced
- the ener.edr. */
- *fp_dhdl=open_dhdl(filename,ir,oenv);
+ fprintf(stdout, " %10s\n", "-kT ln<e^(E/kT)>");
}
else
{
- sprintf(title,"N(%s)",deltag);
- sprintf(label_x,"%s (%s)",deltag,unit_energy);
- sprintf(label_y,"Samples");
- *fp_dhdl=xvgropen_type(filename, title, label_x, label_y, exvggtXNY,oenv);
- sprintf(buf,"T = %g (K), %s = %g", temp, lambda, start_lambda);
- xvgr_subtitle(*fp_dhdl,buf,oenv);
+ fprintf(stdout, "\n");
}
- }
-
- (*hists)+=nblock_hist;
- (*blocks)+=nblock_dh;
- (*nlambdas) = nblock_hist+nblock_dh;
+ fprintf(stdout, "-------------------------------------------------------------------------------\n");
+
+ /* Initiate locals, only used with -sum */
+ expEtot = 0;
+ if (bFee)
+ {
+ beta = 1.0/(BOLTZ*reftemp);
+ snew(fee, nset);
+ }
+ for (i = 0; (i < nset); i++)
+ {
+ aver = edat->s[i].av;
+ stddev = edat->s[i].rmsd;
+ errest = edat->s[i].ee;
+
+ if (bFee)
+ {
+ expE = 0;
+ for (j = 0; (j < edat->nframes); j++)
+ {
+ expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol);
+ }
+ if (bSum)
+ {
+ expEtot += expE/edat->nframes;
+ }
+
+ fee[i] = log(expE/edat->nframes)/beta + aver/nmol;
+ }
+ if (strstr(leg[i], "empera") != NULL)
+ {
+ Temp = aver;
+ }
+ else if (strstr(leg[i], "olum") != NULL)
+ {
+ Vaver = aver;
+ }
+ else if (strstr(leg[i], "essure") != NULL)
+ {
+ Pres = aver;
+ }
+ if (bIsEner[i])
+ {
+ pr_aver = aver/nmol-ezero;
+ pr_stddev = stddev/nmol;
+ pr_errest = errest/nmol;
+ }
+ else
+ {
+ pr_aver = aver;
+ pr_stddev = stddev;
+ pr_errest = errest;
+ }
+
+ /* Multiply the slope in steps with the number of steps taken */
+ totaldrift = (edat->nsteps - 1)*edat->s[i].slope;
+ if (bIsEner[i])
+ {
+ totaldrift /= nmol;
+ }
+
+ fprintf(stdout, "%-24s %10g %10s %10g %10g",
+ leg[i], pr_aver, ee_pr(pr_errest, eebuf), pr_stddev, totaldrift);
+ if (bFee)
+ {
+ fprintf(stdout, " %10g", fee[i]);
+ }
+
+ fprintf(stdout, " (%s)\n", enm[set[i]].unit);
+
+ if (bFluct)
+ {
+ for (j = 0; (j < edat->nframes); j++)
+ {
+ edat->s[i].ener[j] -= aver;
+ }
+ }
+ }
+ if (bSum)
+ {
+ totaldrift = (edat->nsteps - 1)*esum->s[0].slope;
+ fprintf(stdout, "%-24s %10g %10s %10s %10g (%s)",
+ "Total", esum->s[0].av/nmol, ee_pr(esum->s[0].ee/nmol, eebuf),
+ "--", totaldrift/nmol, enm[set[0]].unit);
+ /* pr_aver,pr_stddev,a,totaldrift */
+ if (bFee)
+ {
+ fprintf(stdout, " %10g %10g\n",
+ log(expEtot)/beta + esum->s[0].av/nmol, log(expEtot)/beta);
+ }
+ else
+ {
+ fprintf(stdout, "\n");
+ }
+ }
+
+ /* Do correlation function */
+ if (edat->nframes > 1)
+ {
+ Dt = delta_t/(edat->nframes - 1);
+ }
+ else
+ {
+ Dt = 0;
+ }
+ if (bVisco)
+ {
+ const char* leg[] = { "Shear", "Bulk" };
+ real factor;
+ real **eneset;
+ real **enesum;
+
+ /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
+
+ /* Symmetrise tensor! (and store in first three elements)
+ * And subtract average pressure!
+ */
+ snew(eneset, 12);
+ for (i = 0; i < 12; i++)
+ {
+ snew(eneset[i], edat->nframes);
+ }
+ snew(enesum, 3);
+ for (i = 0; i < 3; i++)
+ {
+ snew(enesum[i], edat->nframes);
+ }
+ for (i = 0; (i < edat->nframes); i++)
+ {
+ eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]);
+ eneset[1][i] = 0.5*(edat->s[2].ener[i]+edat->s[6].ener[i]);
+ eneset[2][i] = 0.5*(edat->s[5].ener[i]+edat->s[7].ener[i]);
+ for (j = 3; j <= 11; j++)
+ {
+ eneset[j][i] = edat->s[j].ener[i];
+ }
+ eneset[11][i] -= Pres;
+ enesum[0][i] = 0.5*(edat->s[1].es[i].sum+edat->s[3].es[i].sum);
+ enesum[1][i] = 0.5*(edat->s[2].es[i].sum+edat->s[6].es[i].sum);
+ enesum[2][i] = 0.5*(edat->s[5].es[i].sum+edat->s[7].es[i].sum);
+ }
+
+ einstein_visco("evisco.xvg", "eviscoi.xvg",
+ 3, edat->nframes, enesum, Vaver, Temp, nsteps, time, oenv);
+
+ /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
+ /* Do it for shear viscosity */
+ strcpy(buf, "Shear Viscosity");
+ low_do_autocorr(corrfn, oenv, buf, edat->nframes, 3,
+ (edat->nframes+1)/2, eneset, Dt,
+ eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0, 1);
+
+ /* Now for bulk viscosity */
+ strcpy(buf, "Bulk Viscosity");
+ low_do_autocorr(corrfn, oenv, buf, edat->nframes, 1,
+ (edat->nframes+1)/2, &(eneset[11]), Dt,
+ eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0, 1);
+
+ factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
+ fp = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv);
+ xvgr_legend(fp, asize(leg), leg, oenv);
+
+ /* Use trapezium rule for integration */
+ integral = 0;
+ intBulk = 0;
+ nout = get_acfnout();
+ if ((nout < 2) || (nout >= edat->nframes/2))
+ {
+ nout = edat->nframes/2;
+ }
+ for (i = 1; (i < nout); i++)
+ {
+ integral += 0.5*(eneset[0][i-1] + eneset[0][i])*factor;
+ intBulk += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
+ fprintf(fp, "%10g %10g %10g\n", (i*Dt), integral, intBulk);
+ }
+ ffclose(fp);
+ }
+ else if (bCorr)
+ {
+ if (bFluct)
+ {
+ strcpy(buf, "Autocorrelation of Energy Fluctuations");
+ }
+ else
+ {
+ strcpy(buf, "Energy Autocorrelation");
+ }
+#if 0
+ do_autocorr(corrfn, oenv, buf, edat->nframes,
+ bSum ? 1 : nset,
+ bSum ? &edat->s[nset-1].ener : eneset,
+ (delta_t/edat->nframes), eacNormal, FALSE);
+#endif
+ }
+ }
+}
+
+static void print_time(FILE *fp, double t)
+{
+ fprintf(fp, "%12.6f", t);
+}
+
+static void print1(FILE *fp, gmx_bool bDp, real e)
+{
+ if (bDp)
+ {
+ fprintf(fp, " %16.12f", e);
+ }
+ else
+ {
+ fprintf(fp, " %10.6f", e);
+ }
+}
+
+static void fec(const char *ene2fn, const char *runavgfn,
+ real reftemp, int nset, int set[], char *leg[],
+ enerdata_t *edat, double time[],
+ const output_env_t oenv)
+{
+ const char * ravgleg[] = {
+ "\\8D\\4E = E\\sB\\N-E\\sA\\N",
+ "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N"
+ };
+ FILE *fp;
+ ener_file_t enx;
+ int nre, timecheck, step, nenergy, nenergy2, maxenergy;
+ int i, j;
+ gmx_bool bCont;
+ real aver, beta;
+ real **eneset2;
+ double dE, sum;
+ gmx_enxnm_t *enm = NULL;
+ t_enxframe *fr;
+ char buf[22];
+
+ /* read second energy file */
+ snew(fr, 1);
+ enm = NULL;
+ enx = open_enx(ene2fn, "r");
+ do_enxnms(enx, &(fr->nre), &enm);
+
+ snew(eneset2, nset+1);
+ nenergy2 = 0;
+ maxenergy = 0;
+ timecheck = 0;
+ do
+ {
+ /* This loop searches for the first frame (when -b option is given),
+ * or when this has been found it reads just one energy frame
+ */
+ do
+ {
+ bCont = do_enx(enx, fr);
+
+ if (bCont)
+ {
+ timecheck = check_times(fr->t);
+ }
+
+ }
+ while (bCont && (timecheck < 0));
+
+ /* Store energies for analysis afterwards... */
+ if ((timecheck == 0) && bCont)
+ {
+ if (fr->nre > 0)
+ {
+ if (nenergy2 >= maxenergy)
+ {
+ maxenergy += 1000;
+ for (i = 0; i <= nset; i++)
+ {
+ srenew(eneset2[i], maxenergy);
+ }
+ }
+ if (fr->t != time[nenergy2])
+ {
+ fprintf(stderr, "\nWARNING time mismatch %g!=%g at frame %s\n",
+ fr->t, time[nenergy2], gmx_step_str(fr->step, buf));
+ }
+ for (i = 0; i < nset; i++)
+ {
+ eneset2[i][nenergy2] = fr->ener[set[i]].e;
+ }
+ nenergy2++;
+ }
+ }
+ }
+ while (bCont && (timecheck == 0));
+
+ /* check */
+ if (edat->nframes != nenergy2)
+ {
+ fprintf(stderr, "\nWARNING file length mismatch %d!=%d\n",
+ edat->nframes, nenergy2);
+ }
+ nenergy = min(edat->nframes, nenergy2);
+
+ /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
+ fp = NULL;
+ if (runavgfn)
+ {
+ fp = xvgropen(runavgfn, "Running average free energy difference",
+ "Time (" unit_time ")", "\\8D\\4E (" unit_energy ")", oenv);
+ xvgr_legend(fp, asize(ravgleg), ravgleg, oenv);
+ }
+ fprintf(stdout, "\n%-24s %10s\n",
+ "Energy", "dF = -kT ln < exp(-(EB-EA)/kT) >A");
+ sum = 0;
+ beta = 1.0/(BOLTZ*reftemp);
+ for (i = 0; i < nset; i++)
+ {
+ if (gmx_strcasecmp(leg[i], enm[set[i]].name) != 0)
+ {
+ fprintf(stderr, "\nWARNING energy set name mismatch %s!=%s\n",
+ leg[i], enm[set[i]].name);
+ }
+ for (j = 0; j < nenergy; j++)
+ {
+ dE = eneset2[i][j] - edat->s[i].ener[j];
+ sum += exp(-dE*beta);
+ if (fp)
+ {
+ fprintf(fp, "%10g %10g %10g\n",
+ time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) );
+ }
+ }
+ aver = -BOLTZ*reftemp*log(sum/nenergy);
+ fprintf(stdout, "%-24s %10g\n", leg[i], aver);
+ }
+ if (fp)
+ {
+ ffclose(fp);
+ }
+ sfree(fr);
+}
+
+
+static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl,
+ const char *filename, gmx_bool bDp,
+ int *blocks, int *hists, int *samples, int *nlambdas,
+ const output_env_t oenv)
+{
+ const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
+ char title[STRLEN], label_x[STRLEN], label_y[STRLEN], legend[STRLEN];
+ char buf[STRLEN];
+ gmx_bool first = FALSE;
+ int nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0;
+ int i, j, k;
+ /* coll data */
+ double temp = 0, start_time = 0, delta_time = 0, start_lambda = 0, delta_lambda = 0;
+ static int setnr = 0;
+ double *native_lambda_vec = NULL;
+ const char **lambda_components = NULL;
+ int n_lambda_vec = 0;
+ gmx_bool changing_lambda = FALSE;
+ int lambda_fep_state;
+
+ /* now count the blocks & handle the global dh data */
+ for (i = 0; i < fr->nblock; i++)
+ {
+ if (fr->block[i].id == enxDHHIST)
+ {
+ nblock_hist++;
+ }
+ else if (fr->block[i].id == enxDH)
+ {
+ nblock_dh++;
+ }
+ else if (fr->block[i].id == enxDHCOLL)
+ {
+ nblock_dhcoll++;
+ if ( (fr->block[i].nsub < 1) ||
+ (fr->block[i].sub[0].type != xdr_datatype_double) ||
+ (fr->block[i].sub[0].nr < 5))
+ {
+ gmx_fatal(FARGS, "Unexpected block data");
+ }
+
+ /* read the data from the DHCOLL block */
+ temp = fr->block[i].sub[0].dval[0];
+ start_time = fr->block[i].sub[0].dval[1];
+ delta_time = fr->block[i].sub[0].dval[2];
+ start_lambda = fr->block[i].sub[0].dval[3];
+ delta_lambda = fr->block[i].sub[0].dval[4];
+ changing_lambda = (delta_lambda != 0);
+ if (fr->block[i].nsub > 1)
+ {
+ lambda_fep_state = fr->block[i].sub[1].ival[0];
+ if (n_lambda_vec == 0)
+ {
+ n_lambda_vec = fr->block[i].sub[1].ival[1];
+ }
+ else
+ {
+ if (n_lambda_vec != fr->block[i].sub[1].ival[1])
+ {
+ gmx_fatal(FARGS,
+ "Unexpected change of basis set in lambda");
+ }
+ }
+ if (lambda_components == NULL)
+ {
+ snew(lambda_components, n_lambda_vec);
+ }
+ if (native_lambda_vec == NULL)
+ {
+ snew(native_lambda_vec, n_lambda_vec);
+ }
+ for (j = 0; j < n_lambda_vec; j++)
+ {
+ native_lambda_vec[j] = fr->block[i].sub[0].dval[5+j];
+ lambda_components[j] =
+ efpt_singular_names[fr->block[i].sub[1].ival[2+j]];
+ }
+ }
+ }
+ }
+
+ if (nblock_hist == 0 && nblock_dh == 0)
+ {
+ /* don't do anything */
+ return;
+ }
+ if (nblock_hist > 0 && nblock_dh > 0)
+ {
+ 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)
+ {
+ /* we have standard, non-histogram data --
+ call open_dhdl to open the file */
+ /* TODO this is an ugly hack that needs to be fixed: this will only
+ work if the order of data is always the same and if we're
+ only using the g_energy compiled with the mdrun that produced
+ the ener.edr. */
+ *fp_dhdl = open_dhdl(filename, ir, oenv);
+ }
+ else
+ {
+ sprintf(title, "N(%s)", deltag);
+ sprintf(label_x, "%s (%s)", deltag, unit_energy);
+ sprintf(label_y, "Samples");
+ *fp_dhdl = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
+ sprintf(buf, "T = %g (K), %s = %g", temp, lambda, start_lambda);
+ xvgr_subtitle(*fp_dhdl, buf, oenv);
+ }
+ }
+
+ (*hists) += nblock_hist;
+ (*blocks) += nblock_dh;
+ (*nlambdas) = nblock_hist+nblock_dh;
/* write the data */
if (nblock_hist > 0)
{
- gmx_large_int_t sum=0;
+ gmx_large_int_t sum = 0;
/* histograms */
- for(i=0;i<fr->nblock;i++)
+ for (i = 0; i < fr->nblock; i++)
{
- t_enxblock *blk=&(fr->block[i]);
- if (blk->id==enxDHHIST)
+ t_enxblock *blk = &(fr->block[i]);
+ if (blk->id == enxDHHIST)
{
- double foreign_lambda, dx;
+ double foreign_lambda, dx;
gmx_large_int_t x0;
- int nhist, derivative;
+ int nhist, derivative;
/* check the block types etc. */
if ( (blk->nsub < 2) ||
{
gmx_fatal(FARGS, "Unexpected block data in file");
}
- foreign_lambda=blk->sub[0].dval[0];
- dx=blk->sub[0].dval[1];
- nhist=blk->sub[1].lval[0];
- derivative=blk->sub[1].lval[1];
- for(j=0;j<nhist;j++)
+ foreign_lambda = blk->sub[0].dval[0];
+ dx = blk->sub[0].dval[1];
+ nhist = blk->sub[1].lval[0];
+ derivative = blk->sub[1].lval[1];
+ for (j = 0; j < nhist; j++)
{
const char *lg[1];
- x0=blk->sub[1].lval[2+j];
+ x0 = blk->sub[1].lval[2+j];
if (!derivative)
{
dhdl, lambda, start_lambda);
}
- lg[0]=legend;
+ lg[0] = legend;
xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
setnr++;
- for(k=0;k<blk->sub[j+2].nr;k++)
+ for (k = 0; k < blk->sub[j+2].nr; k++)
{
- int hist;
+ int hist;
double xmin, xmax;
- hist=blk->sub[j+2].ival[k];
- xmin=(x0+k)*dx;
- xmax=(x0+k+1)*dx;
- fprintf(*fp_dhdl,"%g %d\n%g %d\n", xmin, hist,
+ hist = blk->sub[j+2].ival[k];
+ xmin = (x0+k)*dx;
+ xmax = (x0+k+1)*dx;
+ fprintf(*fp_dhdl, "%g %d\n%g %d\n", xmin, hist,
xmax, hist);
- sum+=hist;
+ sum += hist;
}
/* multiple histogram data blocks in one histogram
mean that the second one is the reverse of the first one:
for dhdl derivatives, it's important to know both the
maximum and minimum values */
- dx=-dx;
+ dx = -dx;
}
}
}
else
{
/* raw dh */
- int len=0;
- char **setnames=NULL;
- int nnames=nblock_dh;
+ int len = 0;
+ char **setnames = NULL;
+ int nnames = nblock_dh;
- for(i=0;i<fr->nblock;i++)
+ for (i = 0; i < fr->nblock; i++)
{
- t_enxblock *blk=&(fr->block[i]);
+ t_enxblock *blk = &(fr->block[i]);
if (blk->id == enxDH)
{
if (len == 0)
{
- len=blk->sub[2].nr;
+ len = blk->sub[2].nr;
}
else
{
- if (len!=blk->sub[2].nr)
+ if (len != blk->sub[2].nr)
{
gmx_fatal(FARGS, "Length inconsistency in dhdl data");
}
}
(*samples) += len;
- for(i=0;i<len;i++)
+ for (i = 0; i < len; i++)
{
- double time=start_time + delta_time*i;
+ double time = start_time + delta_time*i;
- fprintf(*fp_dhdl,"%.4f ", time);
+ fprintf(*fp_dhdl, "%.4f ", time);
- for(j=0;j<fr->nblock;j++)
+ for (j = 0; j < fr->nblock; j++)
{
- t_enxblock *blk=&(fr->block[j]);
+ t_enxblock *blk = &(fr->block[j]);
if (blk->id == enxDH)
{
double value;
if (blk->sub[2].type == xdr_datatype_float)
{
- value=blk->sub[2].fval[i];
+ value = blk->sub[2].fval[i];
}
else
{
- value=blk->sub[2].dval[i];
+ value = blk->sub[2].dval[i];
}
/* we need to decide which data type it is based on the count*/
- if (j==1 && ir->bExpanded)
+ if (j == 1 && ir->bExpanded)
{
- fprintf(*fp_dhdl,"%4d", (int)value); /* if expanded ensembles and zero, this is a state value, it's an integer. We need a cleaner conditional than if j==1! */
- } else {
- if (bDp) {
- fprintf(*fp_dhdl," %#.12g", value); /* print normal precision */
+ fprintf(*fp_dhdl, "%4d", (int)value); /* if expanded ensembles and zero, this is a state value, it's an integer. We need a cleaner conditional than if j==1! */
+ }
+ else
+ {
+ if (bDp)
+ {
+ fprintf(*fp_dhdl, " %#.12g", value); /* print normal precision */
}
else
{
- fprintf(*fp_dhdl," %#.8g", value); /* print normal precision */
+ fprintf(*fp_dhdl, " %#.8g", value); /* print normal precision */
}
}
}
}
-int gmx_energy(int argc,char *argv[])
+int gmx_energy(int argc, char *argv[])
{
- const char *desc[] = {
-
- "[TT]g_energy[tt] extracts energy components or distance restraint",
- "data from an energy file. The user is prompted to interactively",
- "select the desired energy terms.[PAR]",
-
- "Average, RMSD, and drift are calculated with full precision from the",
- "simulation (see printed manual). Drift is calculated by performing",
- "a least-squares fit of the data to a straight line. The reported total drift",
- "is the difference of the fit at the first and last point.",
- "An error estimate of the average is given based on a block averages",
- "over 5 blocks using the full-precision averages. The error estimate",
- "can be performed over multiple block lengths with the options",
- "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
- "[BB]Note[bb] that in most cases the energy files contains averages over all",
- "MD steps, or over many more points than the number of frames in",
- "energy file. This makes the [TT]g_energy[tt] statistics output more accurate",
- "than the [TT].xvg[tt] output. When exact averages are not present in the energy",
- "file, the statistics mentioned above are simply over the single, per-frame",
- "energy values.[PAR]",
-
- "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
-
- "Some fluctuation-dependent properties can be calculated provided",
- "the correct energy terms are selected, and that the command line option",
- "[TT]-fluct_props[tt] is given. The following properties",
- "will be computed:[BR]",
- "Property Energy terms needed[BR]",
- "---------------------------------------------------[BR]",
- "Heat capacity C[SUB]p[sub] (NPT sims): Enthalpy, Temp [BR]",
- "Heat capacity C[SUB]v[sub] (NVT sims): Etot, Temp [BR]",
- "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
- "Isothermal compressibility: Vol, Temp [BR]",
- "Adiabatic bulk modulus: Vol, Temp [BR]",
- "---------------------------------------------------[BR]",
- "You always need to set the number of molecules [TT]-nmol[tt].",
- "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
- "for quantum effects. Use the [TT]g_dos[tt] program if you need that (and you do).[PAR]"
- "When the [TT]-viol[tt] option is set, the time averaged",
- "violations are plotted and the running time-averaged and",
- "instantaneous sum of violations are recalculated. Additionally",
- "running time-averaged and instantaneous distances between",
- "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
-
- "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
- "[TT]-odt[tt] are used for analyzing orientation restraint data.",
- "The first two options plot the orientation, the last three the",
- "deviations of the orientations from the experimental values.",
- "The options that end on an 'a' plot the average over time",
- "as a function of restraint. The options that end on a 't'",
- "prompt the user for restraint label numbers and plot the data",
- "as a function of time. Option [TT]-odr[tt] plots the RMS",
- "deviation as a function of restraint.",
- "When the run used time or ensemble averaged orientation restraints,",
- "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
- "not ensemble-averaged orientations and deviations instead of",
- "the time and ensemble averages.[PAR]",
-
- "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
- "tensor for each orientation restraint experiment. With option",
- "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
-
- "Option [TT]-odh[tt] extracts and plots the free energy data",
- "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
- "from the [TT]ener.edr[tt] file.[PAR]",
-
- "With [TT]-fee[tt] an estimate is calculated for the free-energy",
- "difference with an ideal gas state: [BR]",
- " [GRK]Delta[grk] A = A(N,V,T) - A[SUB]idealgas[sub](N,V,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln][BR]",
- " [GRK]Delta[grk] G = G(N,p,T) - G[SUB]idealgas[sub](N,p,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln][BR]",
- "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
- "the average is over the ensemble (or time in a trajectory).",
- "Note that this is in principle",
- "only correct when averaging over the whole (Boltzmann) ensemble",
- "and using the potential energy. This also allows for an entropy",
- "estimate using:[BR]",
- " [GRK]Delta[grk] S(N,V,T) = S(N,V,T) - S[SUB]idealgas[sub](N,V,T) = ([CHEVRON]U[SUB]pot[sub][chevron] - [GRK]Delta[grk] A)/T[BR]",
- " [GRK]Delta[grk] S(N,p,T) = S(N,p,T) - S[SUB]idealgas[sub](N,p,T) = ([CHEVRON]U[SUB]pot[sub][chevron] + pV - [GRK]Delta[grk] G)/T",
- "[PAR]",
-
- "When a second energy file is specified ([TT]-f2[tt]), a free energy",
- "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
- "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
- "files, and the average is over the ensemble A. The running average",
- "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
- "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
-
- };
- static gmx_bool bSum=FALSE,bFee=FALSE,bPrAll=FALSE,bFluct=FALSE,bDriftCorr=FALSE;
- static gmx_bool bDp=FALSE,bMutot=FALSE,bOrinst=FALSE,bOvec=FALSE,bFluctProps=FALSE;
- static int skip=0,nmol=1,nbmin=5,nbmax=5;
- static real reftemp=300.0,ezero=0;
- t_pargs pa[] = {
- { "-fee", FALSE, etBOOL, {&bFee},
- "Do a free energy estimate" },
- { "-fetemp", FALSE, etREAL,{&reftemp},
- "Reference temperature for free energy calculation" },
- { "-zero", FALSE, etREAL, {&ezero},
- "Subtract a zero-point energy" },
- { "-sum", FALSE, etBOOL, {&bSum},
- "Sum the energy terms selected rather than display them all" },
- { "-dp", FALSE, etBOOL, {&bDp},
- "Print energies in high precision" },
- { "-nbmin", FALSE, etINT, {&nbmin},
- "Minimum number of blocks for error estimate" },
- { "-nbmax", FALSE, etINT, {&nbmax},
- "Maximum number of blocks for error estimate" },
- { "-mutot",FALSE, etBOOL, {&bMutot},
- "Compute the total dipole moment from the components" },
- { "-skip", FALSE, etINT, {&skip},
- "Skip number of frames between data points" },
- { "-aver", FALSE, etBOOL, {&bPrAll},
- "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
- { "-nmol", FALSE, etINT, {&nmol},
- "Number of molecules in your sample: the energies are divided by this number" },
- { "-fluct_props", FALSE, etBOOL, {&bFluctProps},
- "Compute properties based on energy fluctuations, like heat capacity" },
- { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
- "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
- { "-fluc", FALSE, etBOOL, {&bFluct},
- "Calculate autocorrelation of energy fluctuations rather than energy itself" },
- { "-orinst", FALSE, etBOOL, {&bOrinst},
- "Analyse instantaneous orientation data" },
- { "-ovec", FALSE, etBOOL, {&bOvec},
- "Also plot the eigenvectors with [TT]-oten[tt]" }
- };
- const char* drleg[] = {
- "Running average",
- "Instantaneous"
- };
- static const char *setnm[] = {
- "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
- "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
- "Volume", "Pressure"
- };
-
- FILE *out=NULL,*fp_pairs=NULL,*fort=NULL,*fodt=NULL,*foten=NULL;
- FILE *fp_dhdl=NULL;
- FILE **drout;
- ener_file_t fp;
- int timecheck=0;
- gmx_mtop_t mtop;
- gmx_localtop_t *top=NULL;
- t_inputrec ir;
- t_energy **ee;
- enerdata_t edat;
- gmx_enxnm_t *enm=NULL;
- t_enxframe *frame,*fr=NULL;
- int cur=0;
+ const char *desc[] = {
+
+ "[TT]g_energy[tt] extracts energy components or distance restraint",
+ "data from an energy file. The user is prompted to interactively",
+ "select the desired energy terms.[PAR]",
+
+ "Average, RMSD, and drift are calculated with full precision from the",
+ "simulation (see printed manual). Drift is calculated by performing",
+ "a least-squares fit of the data to a straight line. The reported total drift",
+ "is the difference of the fit at the first and last point.",
+ "An error estimate of the average is given based on a block averages",
+ "over 5 blocks using the full-precision averages. The error estimate",
+ "can be performed over multiple block lengths with the options",
+ "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
+ "[BB]Note[bb] that in most cases the energy files contains averages over all",
+ "MD steps, or over many more points than the number of frames in",
+ "energy file. This makes the [TT]g_energy[tt] statistics output more accurate",
+ "than the [TT].xvg[tt] output. When exact averages are not present in the energy",
+ "file, the statistics mentioned above are simply over the single, per-frame",
+ "energy values.[PAR]",
+
+ "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
+
+ "Some fluctuation-dependent properties can be calculated provided",
+ "the correct energy terms are selected, and that the command line option",
+ "[TT]-fluct_props[tt] is given. The following properties",
+ "will be computed:[BR]",
+ "Property Energy terms needed[BR]",
+ "---------------------------------------------------[BR]",
+ "Heat capacity C[SUB]p[sub] (NPT sims): Enthalpy, Temp [BR]",
+ "Heat capacity C[SUB]v[sub] (NVT sims): Etot, Temp [BR]",
+ "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
+ "Isothermal compressibility: Vol, Temp [BR]",
+ "Adiabatic bulk modulus: Vol, Temp [BR]",
+ "---------------------------------------------------[BR]",
+ "You always need to set the number of molecules [TT]-nmol[tt].",
+ "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
+ "for quantum effects. Use the [TT]g_dos[tt] program if you need that (and you do).[PAR]"
+ "When the [TT]-viol[tt] option is set, the time averaged",
+ "violations are plotted and the running time-averaged and",
+ "instantaneous sum of violations are recalculated. Additionally",
+ "running time-averaged and instantaneous distances between",
+ "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
+
+ "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
+ "[TT]-odt[tt] are used for analyzing orientation restraint data.",
+ "The first two options plot the orientation, the last three the",
+ "deviations of the orientations from the experimental values.",
+ "The options that end on an 'a' plot the average over time",
+ "as a function of restraint. The options that end on a 't'",
+ "prompt the user for restraint label numbers and plot the data",
+ "as a function of time. Option [TT]-odr[tt] plots the RMS",
+ "deviation as a function of restraint.",
+ "When the run used time or ensemble averaged orientation restraints,",
+ "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
+ "not ensemble-averaged orientations and deviations instead of",
+ "the time and ensemble averages.[PAR]",
+
+ "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
+ "tensor for each orientation restraint experiment. With option",
+ "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
+
+ "Option [TT]-odh[tt] extracts and plots the free energy data",
+ "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
+ "from the [TT]ener.edr[tt] file.[PAR]",
+
+ "With [TT]-fee[tt] an estimate is calculated for the free-energy",
+ "difference with an ideal gas state: [BR]",
+ " [GRK]Delta[grk] A = A(N,V,T) - A[SUB]idealgas[sub](N,V,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln][BR]",
+ " [GRK]Delta[grk] G = G(N,p,T) - G[SUB]idealgas[sub](N,p,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln][BR]",
+ "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
+ "the average is over the ensemble (or time in a trajectory).",
+ "Note that this is in principle",
+ "only correct when averaging over the whole (Boltzmann) ensemble",
+ "and using the potential energy. This also allows for an entropy",
+ "estimate using:[BR]",
+ " [GRK]Delta[grk] S(N,V,T) = S(N,V,T) - S[SUB]idealgas[sub](N,V,T) = ([CHEVRON]U[SUB]pot[sub][chevron] - [GRK]Delta[grk] A)/T[BR]",
+ " [GRK]Delta[grk] S(N,p,T) = S(N,p,T) - S[SUB]idealgas[sub](N,p,T) = ([CHEVRON]U[SUB]pot[sub][chevron] + pV - [GRK]Delta[grk] G)/T",
+ "[PAR]",
+
+ "When a second energy file is specified ([TT]-f2[tt]), a free energy",
+ "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
+ "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
+ "files, and the average is over the ensemble A. The running average",
+ "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
+ "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
+
+ };
+ static gmx_bool bSum = FALSE, bFee = FALSE, bPrAll = FALSE, bFluct = FALSE, bDriftCorr = FALSE;
+ static gmx_bool bDp = FALSE, bMutot = FALSE, bOrinst = FALSE, bOvec = FALSE, bFluctProps = FALSE;
+ static int skip = 0, nmol = 1, nbmin = 5, nbmax = 5;
+ static real reftemp = 300.0, ezero = 0;
+ t_pargs pa[] = {
+ { "-fee", FALSE, etBOOL, {&bFee},
+ "Do a free energy estimate" },
+ { "-fetemp", FALSE, etREAL, {&reftemp},
+ "Reference temperature for free energy calculation" },
+ { "-zero", FALSE, etREAL, {&ezero},
+ "Subtract a zero-point energy" },
+ { "-sum", FALSE, etBOOL, {&bSum},
+ "Sum the energy terms selected rather than display them all" },
+ { "-dp", FALSE, etBOOL, {&bDp},
+ "Print energies in high precision" },
+ { "-nbmin", FALSE, etINT, {&nbmin},
+ "Minimum number of blocks for error estimate" },
+ { "-nbmax", FALSE, etINT, {&nbmax},
+ "Maximum number of blocks for error estimate" },
+ { "-mutot", FALSE, etBOOL, {&bMutot},
+ "Compute the total dipole moment from the components" },
+ { "-skip", FALSE, etINT, {&skip},
+ "Skip number of frames between data points" },
+ { "-aver", FALSE, etBOOL, {&bPrAll},
+ "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
+ { "-nmol", FALSE, etINT, {&nmol},
+ "Number of molecules in your sample: the energies are divided by this number" },
+ { "-fluct_props", FALSE, etBOOL, {&bFluctProps},
+ "Compute properties based on energy fluctuations, like heat capacity" },
+ { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
+ "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
+ { "-fluc", FALSE, etBOOL, {&bFluct},
+ "Calculate autocorrelation of energy fluctuations rather than energy itself" },
+ { "-orinst", FALSE, etBOOL, {&bOrinst},
+ "Analyse instantaneous orientation data" },
+ { "-ovec", FALSE, etBOOL, {&bOvec},
+ "Also plot the eigenvectors with [TT]-oten[tt]" }
+ };
+ const char * drleg[] = {
+ "Running average",
+ "Instantaneous"
+ };
+ static const char *setnm[] = {
+ "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
+ "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
+ "Volume", "Pressure"
+ };
+
+ FILE *out = NULL, *fp_pairs = NULL, *fort = NULL, *fodt = NULL, *foten = NULL;
+ FILE *fp_dhdl = NULL;
+ FILE **drout;
+ ener_file_t fp;
+ int timecheck = 0;
+ gmx_mtop_t mtop;
+ gmx_localtop_t *top = NULL;
+ t_inputrec ir;
+ t_energy **ee;
+ enerdata_t edat;
+ gmx_enxnm_t *enm = NULL;
+ t_enxframe *frame, *fr = NULL;
+ int cur = 0;
#define NEXT (1-cur)
- int nre,teller,teller_disre,nfr;
- gmx_large_int_t start_step;
- int nor=0,nex=0,norfr=0,enx_i=0;
- real start_t;
- real *bounds=NULL,*violaver=NULL,*oobs=NULL,*orient=NULL,*odrms=NULL;
- int *index=NULL,*pair=NULL,norsel=0,*orsel=NULL,*or_label=NULL;
- int nbounds=0,npairs;
- gmx_bool bDisRe,bDRAll,bORA,bORT,bODA,bODR,bODT,bORIRE,bOTEN,bDHDL;
- gmx_bool bFoundStart,bCont,bEDR,bVisco;
- double sum,sumaver,sumt,ener,dbl;
- double *time=NULL;
- real Vaver;
- int *set=NULL,i,j,k,nset,sss;
- gmx_bool *bIsEner=NULL;
- char **pairleg,**odtleg,**otenleg;
- char **leg=NULL;
- char **nms;
- char *anm_j,*anm_k,*resnm_j,*resnm_k;
- int resnr_j,resnr_k;
- const char *orinst_sub = "@ subtitle \"instantaneous\"\n";
- char buf[256];
- output_env_t oenv;
- t_enxblock *blk=NULL;
- t_enxblock *blk_disre=NULL;
- int ndisre=0;
- int dh_blocks=0, dh_hists=0, dh_samples=0, dh_lambdas=0;
-
- t_filenm fnm[] = {
- { efEDR, "-f", NULL, ffREAD },
- { efEDR, "-f2", NULL, ffOPTRD },
- { efTPX, "-s", NULL, ffOPTRD },
- { efXVG, "-o", "energy", ffWRITE },
- { efXVG, "-viol", "violaver",ffOPTWR },
- { efXVG, "-pairs","pairs", ffOPTWR },
- { efXVG, "-ora", "orienta", ffOPTWR },
- { efXVG, "-ort", "orientt", ffOPTWR },
- { efXVG, "-oda", "orideva", ffOPTWR },
- { efXVG, "-odr", "oridevr", ffOPTWR },
- { efXVG, "-odt", "oridevt", ffOPTWR },
- { efXVG, "-oten", "oriten", ffOPTWR },
- { efXVG, "-corr", "enecorr", ffOPTWR },
- { efXVG, "-vis", "visco", ffOPTWR },
- { efXVG, "-ravg", "runavgdf",ffOPTWR },
- { efXVG, "-odh", "dhdl" ,ffOPTWR }
- };
+ int nre, teller, teller_disre, nfr;
+ gmx_large_int_t start_step;
+ int nor = 0, nex = 0, norfr = 0, enx_i = 0;
+ real start_t;
+ real *bounds = NULL, *violaver = NULL, *oobs = NULL, *orient = NULL, *odrms = NULL;
+ int *index = NULL, *pair = NULL, norsel = 0, *orsel = NULL, *or_label = NULL;
+ int nbounds = 0, npairs;
+ gmx_bool bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN, bDHDL;
+ gmx_bool bFoundStart, bCont, bEDR, bVisco;
+ double sum, sumaver, sumt, ener, dbl;
+ double *time = NULL;
+ real Vaver;
+ int *set = NULL, i, j, k, nset, sss;
+ gmx_bool *bIsEner = NULL;
+ char **pairleg, **odtleg, **otenleg;
+ char **leg = NULL;
+ char **nms;
+ char *anm_j, *anm_k, *resnm_j, *resnm_k;
+ int resnr_j, resnr_k;
+ const char *orinst_sub = "@ subtitle \"instantaneous\"\n";
+ char buf[256];
+ output_env_t oenv;
+ t_enxblock *blk = NULL;
+ t_enxblock *blk_disre = NULL;
+ int ndisre = 0;
+ int dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0;
+
+ t_filenm fnm[] = {
+ { efEDR, "-f", NULL, ffREAD },
+ { efEDR, "-f2", NULL, ffOPTRD },
+ { efTPX, "-s", NULL, ffOPTRD },
+ { efXVG, "-o", "energy", ffWRITE },
+ { efXVG, "-viol", "violaver", ffOPTWR },
+ { efXVG, "-pairs", "pairs", ffOPTWR },
+ { efXVG, "-ora", "orienta", ffOPTWR },
+ { efXVG, "-ort", "orientt", ffOPTWR },
+ { efXVG, "-oda", "orideva", ffOPTWR },
+ { efXVG, "-odr", "oridevr", ffOPTWR },
+ { efXVG, "-odt", "oridevt", ffOPTWR },
+ { efXVG, "-oten", "oriten", ffOPTWR },
+ { efXVG, "-corr", "enecorr", ffOPTWR },
+ { efXVG, "-vis", "visco", ffOPTWR },
+ { efXVG, "-ravg", "runavgdf", ffOPTWR },
+ { efXVG, "-odh", "dhdl", ffOPTWR }
+ };
#define NFILE asize(fnm)
- int npargs;
- t_pargs *ppa;
-
- npargs = asize(pa);
- ppa = add_acf_pargs(&npargs,pa);
- parse_common_args(&argc,argv,
- PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
- NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
-
- bDRAll = opt2bSet("-pairs",NFILE,fnm);
- bDisRe = opt2bSet("-viol",NFILE,fnm) || bDRAll;
- bORA = opt2bSet("-ora",NFILE,fnm);
- bORT = opt2bSet("-ort",NFILE,fnm);
- bODA = opt2bSet("-oda",NFILE,fnm);
- bODR = opt2bSet("-odr",NFILE,fnm);
- bODT = opt2bSet("-odt",NFILE,fnm);
- bORIRE = bORA || bORT || bODA || bODR || bODT;
- bOTEN = opt2bSet("-oten",NFILE,fnm);
- bDHDL = opt2bSet("-odh",NFILE,fnm);
-
- nset = 0;
-
- snew(frame,2);
- fp = open_enx(ftp2fn(efEDR,NFILE,fnm),"r");
- do_enxnms(fp,&nre,&enm);
-
- Vaver = -1;
-
- bVisco = opt2bSet("-vis",NFILE,fnm);
-
- 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;
- }
- }
- 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);
-
- 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");
- }
-
- 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);
- }
- }
- }
- xvgr_legend(foten,bOvec ? nex*12 : nex*3,(const char**)otenleg,oenv);
- }
- }
- 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);
- }
- } else if (bDHDL) {
- get_dhdl_parms(ftp2fn(efTPX,NFILE,fnm),&ir);
- }
-
- /* Initiate energies and set them to zero */
- edat.nsteps = 0;
- edat.npoints = 0;
- edat.nframes = 0;
- edat.step = NULL;
- edat.steps = NULL;
- edat.points = NULL;
- snew(edat.s,nset);
-
- /* Initiate counters */
- teller = 0;
- teller_disre = 0;
- bFoundStart = FALSE;
- start_step = 0;
- start_t = 0;
- do {
- /* This loop searches for the first frame (when -b option is given),
- * or when this has been found it reads just one energy frame
- */
- do {
- bCont = do_enx(fp,&(frame[NEXT]));
- if (bCont) {
- timecheck = check_times(frame[NEXT].t);
- }
- } while (bCont && (timecheck < 0));
-
- if ((timecheck == 0) && bCont) {
- /* We read a valid frame, so we can use it */
- fr = &(frame[NEXT]);
-
- if (fr->nre > 0) {
- /* The frame contains energies, so update cur */
- cur = NEXT;
-
- if (edat.nframes % 1000 == 0)
- {
- srenew(edat.step,edat.nframes+1000);
- memset(&(edat.step[edat.nframes]),0,1000*sizeof(edat.step[0]));
- srenew(edat.steps,edat.nframes+1000);
- memset(&(edat.steps[edat.nframes]),0,1000*sizeof(edat.steps[0]));
- srenew(edat.points,edat.nframes+1000);
- memset(&(edat.points[edat.nframes]),0,1000*sizeof(edat.points[0]));
-
- for(i=0; i<nset; i++)
- {
- srenew(edat.s[i].ener,edat.nframes+1000);
- memset(&(edat.s[i].ener[edat.nframes]),0,
- 1000*sizeof(edat.s[i].ener[0]));
- srenew(edat.s[i].es ,edat.nframes+1000);
- memset(&(edat.s[i].es[edat.nframes]),0,
- 1000*sizeof(edat.s[i].es[0]));
- }
- }
-
- nfr = edat.nframes;
- edat.step[nfr] = fr->step;
-
- if (!bFoundStart)
- {
- bFoundStart = TRUE;
- /* Initiate the previous step data */
- start_step = fr->step;
- start_t = fr->t;
- /* Initiate the energy sums */
- edat.steps[nfr] = 1;
- edat.points[nfr] = 1;
- for(i=0; i<nset; i++)
- {
- sss = set[i];
- edat.s[i].es[nfr].sum = fr->ener[sss].e;
- edat.s[i].es[nfr].sum2 = 0;
- }
- edat.nsteps = 1;
- edat.npoints = 1;
- }
- else
- {
- edat.steps[nfr] = fr->nsteps;
- {
- if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
- {
- if (fr->nsum <= 1)
- {
- edat.points[nfr] = 1;
- for(i=0; i<nset; i++)
- {
- sss = set[i];
- edat.s[i].es[nfr].sum = fr->ener[sss].e;
- edat.s[i].es[nfr].sum2 = 0;
- }
- edat.npoints += 1;
- }
- else
- {
- edat.points[nfr] = fr->nsum;
- for(i=0; i<nset; i++)
- {
- sss = set[i];
- edat.s[i].es[nfr].sum = fr->ener[sss].esum;
- edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
- }
- edat.npoints += fr->nsum;
- }
- }
- else
- {
- /* The interval does not match fr->nsteps:
- * can not do exact averages.
- */
- edat.npoints = 0;
- }
- edat.nsteps = fr->step - start_step + 1;
- }
- }
- for(i=0; i<nset; i++)
- {
- edat.s[i].ener[nfr] = fr->ener[set[i]].e;
- }
- }
- /*
- * Define distance restraint legends. Can only be done after
- * the first frame has been read... (Then we know how many there are)
- */
- blk_disre=find_block_id_enxframe(fr, enxDISRE, NULL);
- if (bDisRe && bDRAll && !leg && blk_disre)
- {
- t_iatom *fa;
- t_iparams *ip;
-
- fa = top->idef.il[F_DISRES].iatoms;
- ip = top->idef.iparams;
- if (blk_disre->nsub != 2 ||
- (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
- {
- gmx_incons("Number of disre sub-blocks not equal to 2");
- }
-
- ndisre=blk_disre->sub[0].nr ;
- if (ndisre != top->idef.il[F_DISRES].nr/3)
- {
- gmx_fatal(FARGS,"Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
- ndisre,top->idef.il[F_DISRES].nr/3);
- }
- snew(pairleg,ndisre);
- for(i=0; i<ndisre; i++)
- {
- snew(pairleg[i],30);
- j=fa[3*i+1];
- k=fa[3*i+2];
- gmx_mtop_atominfo_global(&mtop,j,&anm_j,&resnr_j,&resnm_j);
- gmx_mtop_atominfo_global(&mtop,k,&anm_k,&resnr_k,&resnm_k);
- sprintf(pairleg[i],"%d %s %d %s (%d)",
- resnr_j,anm_j,resnr_k,anm_k,
- ip[fa[3*i]].disres.label);
- }
- set=select_it(ndisre,pairleg,&nset);
- snew(leg,2*nset);
- for(i=0; (i<nset); i++)
- {
- snew(leg[2*i],32);
- sprintf(leg[2*i], "a %s",pairleg[set[i]]);
- snew(leg[2*i+1],32);
- sprintf(leg[2*i+1],"i %s",pairleg[set[i]]);
- }
- xvgr_legend(fp_pairs,2*nset,(const char**)leg,oenv);
- }
-
- /*
- * Store energies for analysis afterwards...
- */
- if (!bDisRe && !bDHDL && (fr->nre > 0)) {
- if (edat.nframes % 1000 == 0) {
- srenew(time,edat.nframes+1000);
- }
- time[edat.nframes] = fr->t;
- edat.nframes++;
- }
- /*
- * Printing time, only when we do not want to skip frames
- */
- if (!skip || teller % skip == 0) {
- if (bDisRe) {
- /*******************************************
- * D I S T A N C E R E S T R A I N T S
- *******************************************/
- if (ndisre > 0)
- {
+ int npargs;
+ t_pargs *ppa;
+
+ npargs = asize(pa);
+ ppa = add_acf_pargs(&npargs, pa);
+ parse_common_args(&argc, argv,
+ PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
+ NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv);
+
+ bDRAll = opt2bSet("-pairs", NFILE, fnm);
+ bDisRe = opt2bSet("-viol", NFILE, fnm) || bDRAll;
+ bORA = opt2bSet("-ora", NFILE, fnm);
+ bORT = opt2bSet("-ort", NFILE, fnm);
+ bODA = opt2bSet("-oda", NFILE, fnm);
+ bODR = opt2bSet("-odr", NFILE, fnm);
+ bODT = opt2bSet("-odt", NFILE, fnm);
+ bORIRE = bORA || bORT || bODA || bODR || bODT;
+ bOTEN = opt2bSet("-oten", NFILE, fnm);
+ bDHDL = opt2bSet("-odh", NFILE, fnm);
+
+ nset = 0;
+
+ snew(frame, 2);
+ fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
+ do_enxnms(fp, &nre, &enm);
+
+ Vaver = -1;
+
+ bVisco = opt2bSet("-vis", NFILE, fnm);
+
+ 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;
+ }
+ }
+ 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);
+ }
+
+ 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");
+ }
+
+ 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);
+ }
+ }
+ }
+ xvgr_legend(foten, bOvec ? nex*12 : nex*3, (const char**)otenleg, oenv);
+ }
+ }
+ 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);
+ }
+ }
+ }
+ else if (bDHDL)
+ {
+ get_dhdl_parms(ftp2fn(efTPX, NFILE, fnm), &ir);
+ }
+
+ /* Initiate energies and set them to zero */
+ edat.nsteps = 0;
+ edat.npoints = 0;
+ edat.nframes = 0;
+ edat.step = NULL;
+ edat.steps = NULL;
+ edat.points = NULL;
+ snew(edat.s, nset);
+
+ /* Initiate counters */
+ teller = 0;
+ teller_disre = 0;
+ bFoundStart = FALSE;
+ start_step = 0;
+ start_t = 0;
+ do
+ {
+ /* This loop searches for the first frame (when -b option is given),
+ * or when this has been found it reads just one energy frame
+ */
+ do
+ {
+ bCont = do_enx(fp, &(frame[NEXT]));
+ if (bCont)
+ {
+ timecheck = check_times(frame[NEXT].t);
+ }
+ }
+ while (bCont && (timecheck < 0));
+
+ if ((timecheck == 0) && bCont)
+ {
+ /* We read a valid frame, so we can use it */
+ fr = &(frame[NEXT]);
+
+ if (fr->nre > 0)
+ {
+ /* The frame contains energies, so update cur */
+ cur = NEXT;
+
+ if (edat.nframes % 1000 == 0)
+ {
+ srenew(edat.step, edat.nframes+1000);
+ memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0]));
+ srenew(edat.steps, edat.nframes+1000);
+ memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0]));
+ srenew(edat.points, edat.nframes+1000);
+ memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0]));
+
+ for (i = 0; i < nset; i++)
+ {
+ srenew(edat.s[i].ener, edat.nframes+1000);
+ memset(&(edat.s[i].ener[edat.nframes]), 0,
+ 1000*sizeof(edat.s[i].ener[0]));
+ srenew(edat.s[i].es, edat.nframes+1000);
+ memset(&(edat.s[i].es[edat.nframes]), 0,
+ 1000*sizeof(edat.s[i].es[0]));
+ }
+ }
+
+ nfr = edat.nframes;
+ edat.step[nfr] = fr->step;
+
+ if (!bFoundStart)
+ {
+ bFoundStart = TRUE;
+ /* Initiate the previous step data */
+ start_step = fr->step;
+ start_t = fr->t;
+ /* Initiate the energy sums */
+ edat.steps[nfr] = 1;
+ edat.points[nfr] = 1;
+ for (i = 0; i < nset; i++)
+ {
+ sss = set[i];
+ edat.s[i].es[nfr].sum = fr->ener[sss].e;
+ edat.s[i].es[nfr].sum2 = 0;
+ }
+ edat.nsteps = 1;
+ edat.npoints = 1;
+ }
+ else
+ {
+ edat.steps[nfr] = fr->nsteps;
+ {
+ if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
+ {
+ if (fr->nsum <= 1)
+ {
+ edat.points[nfr] = 1;
+ for (i = 0; i < nset; i++)
+ {
+ sss = set[i];
+ edat.s[i].es[nfr].sum = fr->ener[sss].e;
+ edat.s[i].es[nfr].sum2 = 0;
+ }
+ edat.npoints += 1;
+ }
+ else
+ {
+ edat.points[nfr] = fr->nsum;
+ for (i = 0; i < nset; i++)
+ {
+ sss = set[i];
+ edat.s[i].es[nfr].sum = fr->ener[sss].esum;
+ edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
+ }
+ edat.npoints += fr->nsum;
+ }
+ }
+ else
+ {
+ /* The interval does not match fr->nsteps:
+ * can not do exact averages.
+ */
+ edat.npoints = 0;
+ }
+ edat.nsteps = fr->step - start_step + 1;
+ }
+ }
+ for (i = 0; i < nset; i++)
+ {
+ edat.s[i].ener[nfr] = fr->ener[set[i]].e;
+ }
+ }
+ /*
+ * Define distance restraint legends. Can only be done after
+ * the first frame has been read... (Then we know how many there are)
+ */
+ blk_disre = find_block_id_enxframe(fr, enxDISRE, NULL);
+ if (bDisRe && bDRAll && !leg && blk_disre)
+ {
+ t_iatom *fa;
+ t_iparams *ip;
+
+ fa = top->idef.il[F_DISRES].iatoms;
+ ip = top->idef.iparams;
+ if (blk_disre->nsub != 2 ||
+ (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
+ {
+ gmx_incons("Number of disre sub-blocks not equal to 2");
+ }
+
+ ndisre = blk_disre->sub[0].nr;
+ if (ndisre != top->idef.il[F_DISRES].nr/3)
+ {
+ gmx_fatal(FARGS, "Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
+ ndisre, top->idef.il[F_DISRES].nr/3);
+ }
+ snew(pairleg, ndisre);
+ for (i = 0; i < ndisre; i++)
+ {
+ snew(pairleg[i], 30);
+ j = fa[3*i+1];
+ k = fa[3*i+2];
+ gmx_mtop_atominfo_global(&mtop, j, &anm_j, &resnr_j, &resnm_j);
+ gmx_mtop_atominfo_global(&mtop, k, &anm_k, &resnr_k, &resnm_k);
+ sprintf(pairleg[i], "%d %s %d %s (%d)",
+ resnr_j, anm_j, resnr_k, anm_k,
+ ip[fa[3*i]].disres.label);
+ }
+ set = select_it(ndisre, pairleg, &nset);
+ snew(leg, 2*nset);
+ for (i = 0; (i < nset); i++)
+ {
+ snew(leg[2*i], 32);
+ sprintf(leg[2*i], "a %s", pairleg[set[i]]);
+ snew(leg[2*i+1], 32);
+ sprintf(leg[2*i+1], "i %s", pairleg[set[i]]);
+ }
+ xvgr_legend(fp_pairs, 2*nset, (const char**)leg, oenv);
+ }
+
+ /*
+ * Store energies for analysis afterwards...
+ */
+ if (!bDisRe && !bDHDL && (fr->nre > 0))
+ {
+ if (edat.nframes % 1000 == 0)
+ {
+ srenew(time, edat.nframes+1000);
+ }
+ time[edat.nframes] = fr->t;
+ edat.nframes++;
+ }
+ /*
+ * Printing time, only when we do not want to skip frames
+ */
+ if (!skip || teller % skip == 0)
+ {
+ if (bDisRe)
+ {
+ /*******************************************
+ * D I S T A N C E R E S T R A I N T S
+ *******************************************/
+ if (ndisre > 0)
+ {
#ifndef GMX_DOUBLE
- float *disre_rt = blk_disre->sub[0].fval;
- float *disre_rm3tav = blk_disre->sub[1].fval;
+ float *disre_rt = blk_disre->sub[0].fval;
+ float *disre_rm3tav = blk_disre->sub[1].fval;
#else
- double *disre_rt = blk_disre->sub[0].dval;
- double *disre_rm3tav = blk_disre->sub[1].dval;
+ double *disre_rt = blk_disre->sub[0].dval;
+ double *disre_rm3tav = blk_disre->sub[1].dval;
#endif
- print_time(out,fr->t);
- if (violaver == NULL)
- snew(violaver,ndisre);
-
- /* Subtract bounds from distances, to calculate violations */
- calc_violations(disre_rt, disre_rm3tav,
- nbounds,pair,bounds,violaver,&sumt,&sumaver);
-
- fprintf(out," %8.4f %8.4f\n",sumaver,sumt);
- if (bDRAll) {
- print_time(fp_pairs,fr->t);
- for(i=0; (i<nset); i++) {
- sss=set[i];
- fprintf(fp_pairs," %8.4f", mypow(disre_rm3tav[sss],minthird));
- fprintf(fp_pairs," %8.4f", disre_rt[sss]);
- }
- fprintf(fp_pairs,"\n");
- }
- teller_disre++;
- }
- }
- else if (bDHDL)
- {
- do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh",NFILE,fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
- }
-
- /*******************************************
- * E N E R G I E S
- *******************************************/
- else {
- if (fr->nre > 0) {
- if (bPrAll)
- {
- /* We skip frames with single points (usually only the first frame),
- * since they would result in an average plot with outliers.
- */
- if (fr->nsum > 1) {
- print_time(out,fr->t);
- print1(out,bDp,fr->ener[set[0]].e);
- print1(out,bDp,fr->ener[set[0]].esum/fr->nsum);
- print1(out,bDp,sqrt(fr->ener[set[0]].eav/fr->nsum));
- fprintf(out,"\n");
- }
- }
- else
- {
- print_time(out,fr->t);
- if (bSum)
- {
- sum = 0;
- for(i=0; i<nset; i++)
- {
- sum += fr->ener[set[i]].e;
- }
- print1(out,bDp,sum/nmol-ezero);
- }
- else
- {
- for(i=0; (i<nset); i++)
- {
- if (bIsEner[i])
- {
- print1(out,bDp,(fr->ener[set[i]].e)/nmol-ezero);
- }
- else
- {
- print1(out,bDp,fr->ener[set[i]].e);
- }
- }
- }
- fprintf(out,"\n");
- }
- }
- blk = find_block_id_enxframe(fr, enx_i, NULL);
- if (bORIRE && blk)
- {
+ print_time(out, fr->t);
+ if (violaver == NULL)
+ {
+ snew(violaver, ndisre);
+ }
+
+ /* Subtract bounds from distances, to calculate violations */
+ calc_violations(disre_rt, disre_rm3tav,
+ nbounds, pair, bounds, violaver, &sumt, &sumaver);
+
+ fprintf(out, " %8.4f %8.4f\n", sumaver, sumt);
+ if (bDRAll)
+ {
+ print_time(fp_pairs, fr->t);
+ for (i = 0; (i < nset); i++)
+ {
+ sss = set[i];
+ fprintf(fp_pairs, " %8.4f", mypow(disre_rm3tav[sss], minthird));
+ fprintf(fp_pairs, " %8.4f", disre_rt[sss]);
+ }
+ fprintf(fp_pairs, "\n");
+ }
+ teller_disre++;
+ }
+ }
+ else if (bDHDL)
+ {
+ do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh", NFILE, fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
+ }
+
+ /*******************************************
+ * E N E R G I E S
+ *******************************************/
+ else
+ {
+ if (fr->nre > 0)
+ {
+ if (bPrAll)
+ {
+ /* We skip frames with single points (usually only the first frame),
+ * since they would result in an average plot with outliers.
+ */
+ if (fr->nsum > 1)
+ {
+ print_time(out, fr->t);
+ print1(out, bDp, fr->ener[set[0]].e);
+ print1(out, bDp, fr->ener[set[0]].esum/fr->nsum);
+ print1(out, bDp, sqrt(fr->ener[set[0]].eav/fr->nsum));
+ fprintf(out, "\n");
+ }
+ }
+ else
+ {
+ print_time(out, fr->t);
+ if (bSum)
+ {
+ sum = 0;
+ for (i = 0; i < nset; i++)
+ {
+ sum += fr->ener[set[i]].e;
+ }
+ print1(out, bDp, sum/nmol-ezero);
+ }
+ else
+ {
+ for (i = 0; (i < nset); i++)
+ {
+ if (bIsEner[i])
+ {
+ print1(out, bDp, (fr->ener[set[i]].e)/nmol-ezero);
+ }
+ else
+ {
+ print1(out, bDp, fr->ener[set[i]].e);
+ }
+ }
+ }
+ fprintf(out, "\n");
+ }
+ }
+ blk = find_block_id_enxframe(fr, enx_i, NULL);
+ if (bORIRE && blk)
+ {
#ifndef GMX_DOUBLE
- xdr_datatype dt=xdr_datatype_float;
+ xdr_datatype dt = xdr_datatype_float;
#else
- xdr_datatype dt=xdr_datatype_double;
+ xdr_datatype dt = xdr_datatype_double;
#endif
- real *vals;
+ real *vals;
- if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
- {
- gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
- }
+ if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
+ {
+ gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
+ }
#ifndef GMX_DOUBLE
- vals=blk->sub[0].fval;
+ vals = blk->sub[0].fval;
#else
- vals=blk->sub[0].dval;
+ vals = blk->sub[0].dval;
#endif
- if (blk->sub[0].nr != (size_t)nor)
- gmx_fatal(FARGS,"Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
- if (bORA || bODA)
- {
- for(i=0; i<nor; i++)
- orient[i] += vals[i];
- }
- if (bODR)
- {
- for(i=0; i<nor; i++)
- odrms[i] += sqr(vals[i]-oobs[i]);
- }
- if (bORT)
- {
- fprintf(fort," %10f",fr->t);
- for(i=0; i<norsel; i++)
- fprintf(fort," %g",vals[orsel[i]]);
- fprintf(fort,"\n");
- }
- if (bODT)
- {
- fprintf(fodt," %10f",fr->t);
- for(i=0; i<norsel; i++)
- fprintf(fodt," %g", vals[orsel[i]]-oobs[orsel[i]]);
- fprintf(fodt,"\n");
- }
- norfr++;
- }
- blk = find_block_id_enxframe(fr, enxORT, NULL);
- if (bOTEN && blk)
- {
+ if (blk->sub[0].nr != (size_t)nor)
+ {
+ gmx_fatal(FARGS, "Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
+ }
+ if (bORA || bODA)
+ {
+ for (i = 0; i < nor; i++)
+ {
+ orient[i] += vals[i];
+ }
+ }
+ if (bODR)
+ {
+ for (i = 0; i < nor; i++)
+ {
+ odrms[i] += sqr(vals[i]-oobs[i]);
+ }
+ }
+ if (bORT)
+ {
+ fprintf(fort, " %10f", fr->t);
+ for (i = 0; i < norsel; i++)
+ {
+ fprintf(fort, " %g", vals[orsel[i]]);
+ }
+ fprintf(fort, "\n");
+ }
+ if (bODT)
+ {
+ fprintf(fodt, " %10f", fr->t);
+ for (i = 0; i < norsel; i++)
+ {
+ fprintf(fodt, " %g", vals[orsel[i]]-oobs[orsel[i]]);
+ }
+ fprintf(fodt, "\n");
+ }
+ norfr++;
+ }
+ blk = find_block_id_enxframe(fr, enxORT, NULL);
+ if (bOTEN && blk)
+ {
#ifndef GMX_DOUBLE
- xdr_datatype dt=xdr_datatype_float;
+ xdr_datatype dt = xdr_datatype_float;
#else
- xdr_datatype dt=xdr_datatype_double;
+ xdr_datatype dt = xdr_datatype_double;
#endif
- real *vals;
+ real *vals;
- if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
- gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
+ if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
+ {
+ gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
+ }
#ifndef GMX_DOUBLE
- vals=blk->sub[0].fval;
+ vals = blk->sub[0].fval;
#else
- vals=blk->sub[0].dval;
+ vals = blk->sub[0].dval;
#endif
- if (blk->sub[0].nr != (size_t)(nex*12))
- gmx_fatal(FARGS,"Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
- blk->sub[0].nr/12, nex);
- fprintf(foten," %10f",fr->t);
- for(i=0; i<nex; i++)
- for(j=0; j<(bOvec?12:3); j++)
- fprintf(foten," %g",vals[i*12+j]);
- fprintf(foten,"\n");
- }
- }
- }
- teller++;
- }
- } while (bCont && (timecheck == 0));
-
- fprintf(stderr,"\n");
- close_enx(fp);
- if (out)
- ffclose(out);
-
- if (bDRAll)
- ffclose(fp_pairs);
-
- if (bORT)
- 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);
- }
- 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);
- }
- 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);
- }
- if (bOTEN)
- 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("\n\nWrote %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
- {
- gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f",NFILE,fnm));
- }
-
- }
- else
- {
- double dt = (frame[cur].t-start_t)/(edat.nframes-1);
- analyse_ener(opt2bSet("-corr",NFILE,fnm),opt2fn("-corr",NFILE,fnm),
- bFee,bSum,bFluct,opt2parg_bSet("-nmol",npargs,ppa),
- bVisco,opt2fn("-vis",NFILE,fnm),
- nmol,
- start_step,start_t,frame[cur].step,frame[cur].t,
- time,reftemp,&edat,
- nset,set,bIsEner,leg,enm,Vaver,ezero,nbmin,nbmax,
- oenv);
- if (bFluctProps)
- calc_fluctuation_props(stdout,bDriftCorr,dt,nset,set,nmol,leg,&edat,
- nbmin,nbmax);
- }
- if (opt2bSet("-f2",NFILE,fnm)) {
- 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);
- }
- thanx(stderr);
-
- return 0;
+ if (blk->sub[0].nr != (size_t)(nex*12))
+ {
+ gmx_fatal(FARGS, "Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
+ blk->sub[0].nr/12, nex);
+ }
+ fprintf(foten, " %10f", fr->t);
+ for (i = 0; i < nex; i++)
+ {
+ for (j = 0; j < (bOvec ? 12 : 3); j++)
+ {
+ fprintf(foten, " %g", vals[i*12+j]);
+ }
+ }
+ fprintf(foten, "\n");
+ }
+ }
+ }
+ teller++;
+ }
+ }
+ while (bCont && (timecheck == 0));
+
+ fprintf(stderr, "\n");
+ close_enx(fp);
+ if (out)
+ {
+ ffclose(out);
+ }
+
+ if (bDRAll)
+ {
+ ffclose(fp_pairs);
+ }
+
+ if (bORT)
+ {
+ 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);
+ }
+ 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);
+ }
+ 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);
+ }
+ if (bOTEN)
+ {
+ 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("\n\nWrote %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
+ {
+ gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f", NFILE, fnm));
+ }
+
+ }
+ else
+ {
+ double dt = (frame[cur].t-start_t)/(edat.nframes-1);
+ analyse_ener(opt2bSet("-corr", NFILE, fnm), opt2fn("-corr", NFILE, fnm),
+ bFee, bSum, bFluct, opt2parg_bSet("-nmol", npargs, ppa),
+ bVisco, opt2fn("-vis", NFILE, fnm),
+ nmol,
+ start_step, start_t, frame[cur].step, frame[cur].t,
+ time, reftemp, &edat,
+ nset, set, bIsEner, leg, enm, Vaver, ezero, nbmin, nbmax,
+ oenv);
+ if (bFluctProps)
+ {
+ calc_fluctuation_props(stdout, bDriftCorr, dt, nset, set, nmol, leg, &edat,
+ nbmin, nbmax);
+ }
+ }
+ if (opt2bSet("-f2", NFILE, fnm))
+ {
+ 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);
+ }
+ thanx(stderr);
+
+ return 0;
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff