mde_delta_h_coll_init(md->dhc, ir);
}
md->fp_dhdl = NULL;
+ snew(md->dE,ir->fepvals->n_lambda);
}
else
{
md->fp_dhdl = fp_dhdl;
+ snew(md->dE,ir->fepvals->n_lambda);
}
if (ir->bSimTemp) {
int i;
return md;
}
+/* print a lambda vector to a string
+ fep = the inputrec's FEP input data
+ i = the index of the lambda vector
+ get_native_lambda = whether to print the native lambda
+ get_names = whether to print the names rather than the values
+ str = the pre-allocated string buffer to print to. */
+static void print_lambda_vector(t_lambda *fep, int i,
+ gmx_bool get_native_lambda, gmx_bool get_names,
+ char *str)
+{
+ size_t nps=0, np;
+ int j,k=0;
+ int Nsep=0;
+
+ for (j=0;j<efptNR;j++)
+ {
+ if (fep->separate_dvdl[j])
+ Nsep ++;
+ }
+ str[0]=0; /* reset the string */
+ if (Nsep > 1)
+ {
+ str += sprintf(str, "("); /* set the opening parenthesis*/
+ }
+ for (j=0;j<efptNR;j++)
+ {
+ if (fep->separate_dvdl[j])
+ {
+ double lam;
+ if (!get_names)
+ {
+ if (get_native_lambda && fep->init_lambda >= 0)
+ {
+ str += sprintf(str,"%.4f", fep->init_lambda);
+ }
+ else
+ {
+ str += sprintf(str,"%.4f", fep->all_lambda[j][i]);
+ }
+ }
+ else
+ {
+ str += sprintf(str,"%s", efpt_singular_names[j]);
+ }
+ /* print comma for the next item */
+ if (k<Nsep-1)
+ {
+ str += sprintf(str,", ");
+ }
+ k++;
+ }
+ }
+ if (Nsep > 1)
+ {
+ /* and add the closing parenthesis */
+ str += sprintf(str, ")");
+ }
+}
+
+
extern FILE *open_dhdl(const char *filename,const t_inputrec *ir,
const output_env_t oenv)
{
*lambdastate="\\lambda state",*remain="remaining";
char title[STRLEN],label_x[STRLEN],label_y[STRLEN];
int i,np,nps,nsets,nsets_de,nsetsbegin;
- t_lambda *fep;
+ int n_lambda_terms=0;
+ t_lambda *fep=ir->fepvals; /* for simplicity */
+ t_expanded *expand=ir->expandedvals;
char **setname;
- char buf[STRLEN];
+ char buf[STRLEN], lambda_vec_str[STRLEN], lambda_name_str[STRLEN];
int bufplace=0;
int nsets_dhdl = 0;
int s = 0;
int nsetsextend;
+ gmx_bool write_pV = FALSE;
- /* for simplicity */
- fep = ir->fepvals;
+ /* count the number of different lambda terms */
+ for (i=0;i<efptNR;i++)
+ {
+ if (fep->separate_dvdl[i])
+ {
+ n_lambda_terms++;
+ }
+ }
if (fep->n_lambda == 0)
{
}
if (ir->efep != efepSLOWGROWTH)
{
- if (fep->n_lambda == 0)
+ if ( (fep->init_lambda >= 0) && (n_lambda_terms == 1 ))
{
- sprintf(&(buf[bufplace]),"%s = %g",
- lambda,fep->init_lambda);
+ /* compatibility output */
+ sprintf(&(buf[bufplace]),"%s = %.4f", lambda,fep->init_lambda);
}
else
{
- sprintf(&(buf[bufplace]),"%s = %d",
- lambdastate,fep->init_fep_state);
+ print_lambda_vector(fep, fep->init_fep_state, TRUE, FALSE,
+ lambda_vec_str);
+ print_lambda_vector(fep, fep->init_fep_state, TRUE, TRUE,
+ lambda_name_str);
+ sprintf(&(buf[bufplace]),"%s %d: %s = %s",
+ lambdastate,fep->init_fep_state,
+ lambda_name_str, lambda_vec_str);
}
}
xvgr_subtitle(fp,buf,oenv);
- for (i=0;i<efptNR;i++)
+
+ nsets_dhdl=0;
+ if (fep->dhdl_derivatives == edhdlderivativesYES)
{
- if (fep->separate_dvdl[i]) {nsets_dhdl++;}
+ nsets_dhdl = n_lambda_terms;
}
-
/* count the number of delta_g states */
nsets_de = fep->n_lambda;
nsets = nsets_dhdl + nsets_de; /* dhdl + fep differences */
- if (fep->n_lambda>0 && ir->bExpanded)
+ if (fep->n_lambda>0 && (expand->elmcmove > elmcmoveNO))
{
nsets += 1; /*add fep state for expanded ensemble */
}
}
nsetsextend = nsets;
- if ((ir->epc!=epcNO) && (fep->n_lambda>0))
+ if ((ir->epc!=epcNO) && (fep->n_lambda>0) && (fep->init_lambda < 0))
{
- nsetsextend += 1; /* for PV term, other terms possible if required for the reduced potential (only needed with foreign lambda) */
+ nsetsextend += 1; /* for PV term, other terms possible if required for
+ the reduced potential (only needed with foreign
+ lambda, and only output when init_lambda is not
+ set in order to maintain compatibility of the
+ dhdl.xvg file) */
+ write_pV = TRUE;
}
snew(setname,nsetsextend);
- if (ir->bExpanded)
+ if (expand->elmcmove > elmcmoveNO)
{
/* state for the fep_vals, if we have alchemical sampling */
sprintf(buf,"%s","Thermodynamic state");
s+=1;
}
- for (i=0;i<efptNR;i++)
+ if (fep->dhdl_derivatives == edhdlderivativesYES)
{
- if (fep->separate_dvdl[i]) {
- sprintf(buf,"%s (%s)",dhdl,efpt_names[i]);
- setname[s] = strdup(buf);
- s+=1;
+ for (i=0;i<efptNR;i++)
+ {
+ if (fep->separate_dvdl[i]) {
+
+ if ( (fep->init_lambda >= 0) && (n_lambda_terms == 1 ))
+ {
+ /* compatibility output */
+ sprintf(buf,"%s %s %.4f",dhdl,lambda, fep->init_lambda);
+ }
+ else
+ {
+ double lam=fep->init_lambda;
+ if (fep->init_lambda < 0)
+ {
+ lam=fep->all_lambda[i][fep->init_fep_state];
+ }
+ sprintf(buf,"%s %s = %.4f",dhdl, efpt_singular_names[i],
+ lam);
+ }
+ setname[s] = strdup(buf);
+ s+=1;
+ }
}
}
* from this xvg legend.
*/
- if (ir->bExpanded) {
+ if (expand->elmcmove > elmcmoveNO) {
nsetsbegin = 1; /* for including the expanded ensemble */
} else {
nsetsbegin = 0;
}
nsetsbegin += nsets_dhdl;
- for(s=nsetsbegin; s<nsets; s++)
+ for(i=0; i<fep->n_lambda; i++)
{
- nps = sprintf(buf,"%s %s (",deltag,lambda);
- for (i=0;i<efptNR;i++)
+ print_lambda_vector(fep, i, FALSE, FALSE, lambda_vec_str);
+ if ( (fep->init_lambda >= 0) && (n_lambda_terms == 1 ))
{
- if (fep->separate_dvdl[i])
- {
- np = sprintf(&buf[nps],"%g,",fep->all_lambda[i][s-(nsetsbegin)]);
- nps += np;
- }
+ /* for compatible dhdl.xvg files */
+ nps = sprintf(buf,"%s %s %s",deltag,lambda, lambda_vec_str);
}
- if (ir->bSimTemp)
+ else
{
- /* print the temperature for this state if doing simulated annealing */
- sprintf(&buf[nps],"T = %g (%s))",ir->simtempvals->temperatures[s-(nsetsbegin)],unit_temp_K);
+ nps = sprintf(buf,"%s %s to %s",deltag,lambda, lambda_vec_str);
}
- else
+
+ if (ir->bSimTemp)
{
- sprintf(&buf[nps-1],")"); /* -1 to overwrite the last comma */
+ /* print the temperature for this state if doing simulated annealing */
+ sprintf(&buf[nps],"T = %g (%s)",
+ ir->simtempvals->temperatures[s-(nsetsbegin)],
+ unit_temp_K);
}
setname[s] = strdup(buf);
+ s++;
}
- if (ir->epc!=epcNO) {
+ if (write_pV) {
np = sprintf(buf,"pV (%s)",unit_energy);
- setname[nsetsextend-1] = strdup(buf); /* the first entry after nsets */
+ setname[nsetsextend-1] = strdup(buf); /* the first entry after
+ nsets */
}
xvgr_legend(fp,nsetsextend,(const char **)setname,oenv);
real eee[egNR];
real ecopy[F_NRE];
double store_dhdl[efptNR];
- double *dE=NULL;
real store_energy=0;
real tmp;
ebin_increase_count(md->ebin,bSum);
/* BAR + thermodynamic integration values */
- if ((md->fp_dhdl || md->dhc) && bDoDHDL && (enerd->n_lambda > 0))
+ if ((md->fp_dhdl || md->dhc) && bDoDHDL)
{
- snew(dE,enerd->n_lambda-1);
for(i=0; i<enerd->n_lambda-1; i++) {
- dE[i] = enerd->enerpart_lambda[i+1]-enerd->enerpart_lambda[0]; /* zero for simulated tempering */
+ /* zero for simulated tempering */
+ md->dE[i] = enerd->enerpart_lambda[i+1]-enerd->enerpart_lambda[0];
if (md->temperatures!=NULL)
{
/* MRS: is this right, given the way we have defined the exchange probabilities? */
/* is this even useful to have at all? */
- dE[i] += (md->temperatures[i]/md->temperatures[state->fep_state]-1.0)*enerd->term[F_EKIN];
+ md->dE[i] += (md->temperatures[i]/
+ md->temperatures[state->fep_state]-1.0)*
+ enerd->term[F_EKIN];
}
}
- }
-
- if (md->fp_dhdl && bDoDHDL)
- {
- fprintf(md->fp_dhdl,"%.4f",time);
- /* the current free energy state */
- /* print the current state if we are doing expanded ensemble */
- if (expand->elmcmove > elmcmoveNO) {
- fprintf(md->fp_dhdl," %4d",state->fep_state);
- }
- /* total energy (for if the temperature changes */
- if (fep->bPrintEnergy)
+ if (md->fp_dhdl)
{
- store_energy = enerd->term[F_ETOT];
- fprintf(md->fp_dhdl," %#.8g",store_energy);
- }
+ fprintf(md->fp_dhdl,"%.4f",time);
+ /* the current free energy state */
- for (i=0;i<efptNR;i++)
- {
- if (fep->separate_dvdl[i])
+ /* print the current state if we are doing expanded ensemble */
+ if (expand->elmcmove > elmcmoveNO) {
+ fprintf(md->fp_dhdl," %4d",state->fep_state);
+ }
+ /* total energy (for if the temperature changes */
+ if (fep->bPrintEnergy)
{
- fprintf(md->fp_dhdl," %#.8g",enerd->term[F_DVDL+i]); /* assumes F_DVDL is first */
+ store_energy = enerd->term[F_ETOT];
+ fprintf(md->fp_dhdl," %#.8g",store_energy);
}
- }
- for(i=1; i<enerd->n_lambda; i++)
- {
- fprintf(md->fp_dhdl," %#.8g",dE[i-1]);
+ if (fep->dhdl_derivatives == edhdlderivativesYES)
+ {
+ for (i=0;i<efptNR;i++)
+ {
+ if (fep->separate_dvdl[i])
+ {
+ /* assumes F_DVDL is first */
+ fprintf(md->fp_dhdl," %#.8g",enerd->term[F_DVDL+i]);
+ }
+ }
+ }
+ for(i=0;i<fep->n_lambda;i++)
+ {
+ fprintf(md->fp_dhdl," %#.8g",md->dE[i]);
+ }
+ if ((md->epc!=epcNO) &&
+ (enerd->n_lambda > 0) &&
+ (fep->init_lambda<0))
+ {
+ fprintf(md->fp_dhdl," %#.8g",pv); /* PV term only needed when
+ there are alternate state
+ lambda and we're not in
+ compatibility mode */
+ }
+ fprintf(md->fp_dhdl,"\n");
+ /* and the binary free energy output */
}
- if ((md->epc!=epcNO) && (enerd->n_lambda > 0))
- {
- fprintf(md->fp_dhdl," %#.8g",pv); /* PV term only needed when there are alternate state lambda */
- }
- fprintf(md->fp_dhdl,"\n");
- /* and the binary free energy output */
- }
- if (md->dhc && bDoDHDL)
- {
- int idhdl = 0;
- for (i=0;i<efptNR;i++)
+ if (md->dhc && bDoDHDL)
{
- if (fep->separate_dvdl[i])
+ int idhdl = 0;
+ for (i=0;i<efptNR;i++)
{
- store_dhdl[idhdl] = enerd->term[F_DVDL+i]; /* assumes F_DVDL is first */
- idhdl+=1;
+ if (fep->separate_dvdl[i])
+ {
+ /* assumes F_DVDL is first */
+ store_dhdl[idhdl] = enerd->term[F_DVDL+i];
+ idhdl+=1;
+ }
}
+ store_energy = enerd->term[F_ETOT];
+ /* store_dh is dE */
+ mde_delta_h_coll_add_dh(md->dhc,
+ (double)state->fep_state,
+ store_energy,
+ pv,
+ store_dhdl,
+ md->dE,
+ time);
}
- /* store_dh is dE */
- mde_delta_h_coll_add_dh(md->dhc,
- (double)state->fep_state,
- store_energy,
- pv,
- (expand->elamstats>elamstatsNO),
- (fep->bPrintEnergy),
- (md->epc!=epcNO),
- idhdl,
- fep->n_lambda,
- store_dhdl,
- dE,
- time);
- }
- if ((md->fp_dhdl || md->dhc) && bDoDHDL && (enerd->n_lambda >0))
- {
- sfree(dE);
}
}
/* initialize the delta_h list */
static void mde_delta_h_init(t_mde_delta_h *dh, int nbins,
- double dx, unsigned int ndhmax)
+ double dx, unsigned int ndhmax,
+ int type, int derivative, int nlambda,
+ double *lambda)
{
int i;
+ dh->type=type;
+ dh->derivative=derivative;
+ dh->lambda=lambda;
+ dh->nlambda=nlambda;
+
+ snew(dh->lambda, nlambda);
+ for(i=0;i<nlambda;i++)
+ {
+ dh->lambda[i] = lambda[i];
+ }
+
+
+ snew(dh->subblock_meta_d, dh->nlambda+1);
+
dh->ndhmax=ndhmax+2;
for(i=0;i<2;i++)
{
/* first check which type we should use: histogram or raw data */
if (dh->nhist == 0)
{
- unsigned int i;
+ int i;
/* We write raw data.
- Raw data consists of 3 subblocks: a block with the
- the foreign lambda, and the data itself */
+ Raw data consists of 3 subblocks: an int metadata block
+ with type and derivative index, a foreign lambda block
+ and and the data itself */
add_subblocks_enxblock(blk, 3);
blk->id=enxDH;
/* subblock 1 */
- dh->subblock_i[0]=dh->derivative ? 1 : 0; /* derivative type */
- blk->sub[0].nr=1;
+ dh->subblock_meta_i[0]=dh->type; /* block data type */
+ dh->subblock_meta_i[1]=dh->derivative; /* derivative direction if
+ applicable (in indices
+ starting from first coord in
+ the main delta_h_coll) */
+ blk->sub[0].nr=2;
blk->sub[0].type=xdr_datatype_int;
- blk->sub[0].ival=dh->subblock_i;
+ blk->sub[0].ival=dh->subblock_meta_i;
/* subblock 2 */
- dh->subblock_d[0]=dh->lambda;
- blk->sub[1].nr=1;
+ for(i=0;i<dh->nlambda;i++)
+ {
+ dh->subblock_meta_d[i]=dh->lambda[i];
+ }
+ blk->sub[1].nr=dh->nlambda;
blk->sub[1].type=xdr_datatype_double;
- blk->sub[1].dval=dh->subblock_d;
+ blk->sub[1].dval=dh->subblock_meta_d;
/* subblock 3 */
/* check if there's actual data to be written. */
/*if (dh->ndh > 1)*/
if (dh->ndh > 0)
{
+ unsigned int i;
+
blk->sub[2].nr=dh->ndh;
/* For F@H for now. */
#undef GMX_DOUBLE
blk->sub[2].fval=NULL;
#else
blk->sub[2].type=xdr_datatype_double;
-#endif
blk->sub[2].dval=NULL;
+#endif
}
}
else
{
int nhist_written=0;
int i;
+ int k;
+ /* TODO histogram metadata */
/* check if there's actual data to be written. */
if (dh->ndh > 1)
{
add_subblocks_enxblock(blk, nhist_written+2);
blk->id=enxDHHIST;
- /* subblock 1: the foreign lambda value + the histogram spacing */
- dh->subblock_d[0]=dh->lambda;
- dh->subblock_d[1]=dh->dx;
- blk->sub[0].nr=2;
+ /* subblock 1: the lambda value + the histogram spacing */
+ if (dh->nlambda == 1)
+ {
+ /* for backward compatibility */
+ dh->subblock_meta_d[0]=dh->lambda[0];
+ }
+ else
+ {
+ dh->subblock_meta_d[0]=-1;
+ for(i=0;i<dh->nlambda;i++)
+ {
+ dh->subblock_meta_d[2+i]=dh->lambda[i];
+ }
+ }
+ dh->subblock_meta_d[1]=dh->dx;
+ blk->sub[0].nr = 2+ ((dh->nlambda>1) ? dh->nlambda : 0);
blk->sub[0].type=xdr_datatype_double;
- blk->sub[0].dval=dh->subblock_d;
+ blk->sub[0].dval=dh->subblock_meta_d;
/* subblock 2: the starting point(s) as a long integer */
- dh->subblock_l[0]=nhist_written;
- dh->subblock_l[1]=dh->derivative ? 1 : 0;
+ dh->subblock_meta_l[0]=nhist_written;
+ dh->subblock_meta_l[1]=dh->type; /*dh->derivative ? 1 : 0;*/
+ k=2;
for(i=0;i<nhist_written;i++)
- dh->subblock_l[2+i]=dh->x0[i];
+ dh->subblock_meta_l[k++]=dh->x0[i];
+ /* append the derivative data */
+ dh->subblock_meta_l[k++]=dh->derivative;
- blk->sub[1].nr=nhist_written+2;
+ blk->sub[1].nr=nhist_written+3;
blk->sub[1].type=xdr_datatype_large_int;
- blk->sub[1].lval=dh->subblock_l;
+ blk->sub[1].lval=dh->subblock_meta_l;
/* subblock 3 + 4 : the histogram data */
for(i=0;i<nhist_written;i++)
/* initialize the collection*/
void mde_delta_h_coll_init(t_mde_delta_h_coll *dhc, const t_inputrec *ir)
{
- int i;
+ int i,j,n;
double lambda;
+ double *lambda_vec;
int ndhmax=ir->nstenergy/ir->nstcalcenergy;
+ t_lambda *fep=ir->fepvals;
dhc->temperature=ir->opts.ref_t[0]; /* only store system temperature */
dhc->start_time=0.;
dhc->delta_time=ir->delta_t*ir->fepvals->nstdhdl;
dhc->start_time_set=FALSE;
- /* for continuous change of lambda values */
+ /* this is the compatibility lambda value. If it is >=0, it is valid,
+ and there is either an old-style lambda or a slow growth simulation. */
dhc->start_lambda=ir->fepvals->init_lambda;
+ /* for continuous change of lambda values */
dhc->delta_lambda=ir->fepvals->delta_lambda*ir->fepvals->nstdhdl;
- /* total number of raw data points in the sample */
- dhc->ndh = 0;
-
- /* include one more for the specification of the state, by lambda or fep_state, store as double for now*/
- if (ir->expandedvals->elamstats > elamstatsNO) {
- dhc->ndh +=1;
+ if (dhc->start_lambda < 0)
+ {
+ /* create the native lambda vectors */
+ dhc->lambda_index=fep->init_fep_state;
+ dhc->n_lambda_vec=0;
+ for(i=0;i<efptNR;i++)
+ {
+ if (fep->separate_dvdl[i])
+ {
+ dhc->n_lambda_vec++;
+ }
+ }
+ snew(dhc->native_lambda_vec, dhc->n_lambda_vec);
+ snew(dhc->native_lambda_components, dhc->n_lambda_vec);
+ j=0;
+ for(i=0;i<efptNR;i++)
+ {
+ if (fep->separate_dvdl[i])
+ {
+ dhc->native_lambda_components[j]=i;
+ if (fep->init_fep_state >=0 &&
+ fep->init_fep_state < fep->n_lambda)
+ {
+ dhc->native_lambda_vec[j]=
+ fep->all_lambda[i][fep->init_fep_state];
+ }
+ else
+ {
+ dhc->native_lambda_vec[j]=-1;
+ }
+ j++;
+ }
+ }
}
-
- /* whether to print energies */
- if (ir->fepvals->bPrintEnergy) {
- dhc->ndh += 1;
+ else
+ {
+ /* don't allocate the meta-data subblocks for lambda vectors */
+ dhc->native_lambda_vec=NULL;
+ dhc->n_lambda_vec=0;
+ dhc->native_lambda_components=0;
+ dhc->lambda_index=-1;
}
+ /* allocate metadata subblocks */
+ snew(dhc->subblock_d, 5 + dhc->n_lambda_vec);
+ snew(dhc->subblock_i, 1 + dhc->n_lambda_vec);
+
+ /* now decide which data to write out */
+ dhc->nlambda=0;
+ dhc->ndhdl=0;
+ dhc->dh_expanded=NULL;
+ dhc->dh_energy=NULL;
+ dhc->dh_pv=NULL;
+
+ /* total number of raw data point collections in the sample */
+ dhc->ndh = 0;
- /* add the dhdl's */
- for (i=0;i<efptNR;i++)
{
- if (ir->fepvals->separate_dvdl[i])
+ gmx_bool bExpanded=FALSE;
+ gmx_bool bEnergy=FALSE;
+ gmx_bool bPV=FALSE;
+ int n_lambda_components=0;
+
+ /* first count the number of states */
+
+ /* add the dhdl's */
+ if (fep->dhdl_derivatives == edhdlderivativesYES)
{
- dhc->ndh+=1;
+ for (i=0;i<efptNR;i++)
+ {
+ if (ir->fepvals->separate_dvdl[i])
+ {
+ dhc->ndh+=1;
+ dhc->ndhdl+=1;
+ }
+ }
}
- }
+ /* add the lambdas */
+ dhc->nlambda = ir->fepvals->n_lambda;
+ dhc->ndh += dhc->nlambda;
+ /* another compatibility check */
+ if (dhc->start_lambda < 0)
+ {
+ /* include one more for the specification of the state, by lambda or
+ fep_state*/
+ if (ir->expandedvals->elmcmove > elmcmoveNO) {
+ dhc->ndh +=1;
+ bExpanded=TRUE;
+ }
+ /* whether to print energies */
+ if (ir->fepvals->bPrintEnergy) {
+ dhc->ndh += 1;
+ bEnergy=TRUE;
+ }
+ if (ir->epc > epcNO) {
+ dhc->ndh += 1; /* include pressure-volume work */
+ bPV=TRUE;
+ }
+ }
+ /* allocate them */
+ snew(dhc->dh, dhc->ndh);
+
+ /* now initialize them */
+ /* the order, for now, must match that of the dhdl.xvg file because of
+ how g_energy -odh is implemented */
+ n=0;
+ if (bExpanded)
+ {
+ dhc->dh_expanded=dhc->dh+n;
+ mde_delta_h_init(dhc->dh+n, ir->fepvals->dh_hist_size,
+ ir->fepvals->dh_hist_spacing, ndhmax,
+ dhbtEXPANDED, 0, 0, NULL);
+ n++;
+ }
+ if (bEnergy)
+ {
+ dhc->dh_energy=dhc->dh+n;
+ mde_delta_h_init(dhc->dh+n, ir->fepvals->dh_hist_size,
+ ir->fepvals->dh_hist_spacing, ndhmax,
+ dhbtEN, 0, 0, NULL);
+ n++;
+ }
+ /* add the dhdl's */
+ n_lambda_components=0;
+ if (fep->dhdl_derivatives == edhdlderivativesYES)
+ {
+ dhc->dh_dhdl = dhc->dh + n;
+ for (i=0;i<efptNR;i++)
+ {
+ if (ir->fepvals->separate_dvdl[i])
+ {
+ /* we give it init_lambda for compatibility */
+ mde_delta_h_init(dhc->dh+n, ir->fepvals->dh_hist_size,
+ ir->fepvals->dh_hist_spacing, ndhmax,
+ dhbtDHDL, n_lambda_components, 1,
+ &(fep->init_lambda));
+ n++;
+ n_lambda_components++;
+ }
+ }
+ }
+ else
+ {
+ for (i=0;i<efptNR;i++)
+ {
+ if (ir->fepvals->separate_dvdl[i])
+ {
+ n_lambda_components++; /* count the components */
+ }
+ }
- /* add the lambdas */
- dhc->ndh += ir->fepvals->n_lambda;
+ }
+ /* add the lambdas */
+ dhc->dh_du = dhc->dh + n;
+ snew(lambda_vec, n_lambda_components);
+ for(i=0;i<ir->fepvals->n_lambda;i++)
+ {
+ int k=0;
- if (ir->epc > epcNO) {
- dhc->ndh += 1; /* include pressure-volume work */
- }
+ for(j=0;j<efptNR;j++)
+ {
+ if (ir->fepvals->separate_dvdl[j])
+ {
+ lambda_vec[k++] = fep->all_lambda[j][i];
+ }
+ }
- snew(dhc->dh, dhc->ndh);
- for(i=0;i<dhc->ndh;i++)
- {
- mde_delta_h_init(dhc->dh+i, ir->fepvals->dh_hist_size,
- ir->fepvals->dh_hist_spacing, ndhmax);
+ mde_delta_h_init(dhc->dh+n, ir->fepvals->dh_hist_size,
+ ir->fepvals->dh_hist_spacing, ndhmax,
+ dhbtDH, 0, n_lambda_components, lambda_vec);
+ n++;
+ }
+ sfree(lambda_vec);
+ if (bPV)
+ {
+ dhc->dh_pv=dhc->dh+n;
+ mde_delta_h_init(dhc->dh+n, ir->fepvals->dh_hist_size,
+ ir->fepvals->dh_hist_spacing, ndhmax,
+ dhbtPV, 0, 0, NULL);
+ n++;
+ }
}
}
double fep_state,
double energy,
double pV,
- int bExpanded,
- int bPrintEnergy,
- int bPressure,
- int ndhdl,
- int nlambda,
double *dhdl,
double *foreign_dU,
double time)
{
- int i,n;
+ int i;
if (!dhc->start_time_set)
{
dhc->start_time=time;
}
- n = 0;
- if (bExpanded)
+ for (i=0;i<dhc->ndhdl;i++)
{
- mde_delta_h_add_dh(dhc->dh+n,fep_state,time);
- n++;
+ mde_delta_h_add_dh(dhc->dh_dhdl+i, dhdl[i], time);
}
- if (bPrintEnergy)
+ for (i=0;i<dhc->nlambda;i++)
{
- mde_delta_h_add_dh(dhc->dh+n,energy,time);
- n++;
+ mde_delta_h_add_dh(dhc->dh_du+i, foreign_dU[i], time);
}
- for (i=0;i<ndhdl;i++)
+ if (dhc->dh_pv != NULL)
{
- mde_delta_h_add_dh(dhc->dh+n, dhdl[i], time);
- n++;
+ mde_delta_h_add_dh(dhc->dh_pv, pV, time);
}
- for (i=0;i<nlambda;i++)
+ if (dhc->dh_energy != NULL)
{
- mde_delta_h_add_dh(dhc->dh+n, foreign_dU[i], time);
- n++;
+ mde_delta_h_add_dh(dhc->dh_energy,energy,time);
}
- if (bPressure)
+ if (dhc->dh_expanded != NULL)
{
- mde_delta_h_add_dh(dhc->dh+n, pV, time);
- n++;
+ mde_delta_h_add_dh(dhc->dh_expanded,fep_state,time);
}
-}
-/* write the data associated with all the du blocks, but not the blocks
- themselves. Essentially, the metadata. Or -- is this generated every time?*/
+}
+/* write the metadata associated with all the du blocks, and call
+ handle_block to write out all the du blocks */
void mde_delta_h_coll_handle_block(t_mde_delta_h_coll *dhc,
t_enxframe *fr, int nblock)
{
add_blocks_enxframe(fr, nblock);
blk=fr->block + (nblock-1);
- add_subblocks_enxblock(blk, 1);
+ /* only allocate lambda vector component blocks if they must be written out
+ for backward compatibility */
+ if (dhc->native_lambda_components!=NULL)
+ {
+ add_subblocks_enxblock(blk, 2);
+ }
+ else
+ {
+ add_subblocks_enxblock(blk, 1);
+ }
dhc->subblock_d[0] = dhc->temperature; /* temperature */
dhc->subblock_d[1] = dhc->start_time; /* time of first sample */
dhc->subblock_d[2] = dhc->delta_time; /* time difference between samples */
- dhc->subblock_d[3] = dhc->start_lambda; /* lambda at starttime */
+ dhc->subblock_d[3] = dhc->start_lambda; /* old-style lambda at starttime */
dhc->subblock_d[4] = dhc->delta_lambda; /* lambda diff. between samples */
+ /* set the lambda vector components if they exist */
+ if (dhc->native_lambda_components!=NULL)
+ {
+ for(i=0;i<dhc->n_lambda_vec;i++)
+ {
+ dhc->subblock_d[5+i] = dhc->native_lambda_vec[i];
+ }
+ }
blk->id=enxDHCOLL;
- blk->sub[0].nr=5;
+ blk->sub[0].nr=5 + dhc->n_lambda_vec;
blk->sub[0].type=xdr_datatype_double;
blk->sub[0].dval=dhc->subblock_d;
+ if (dhc->native_lambda_components != NULL)
+ {
+ dhc->subblock_i[0] = dhc->lambda_index;
+ /* set the lambda vector component IDs if they exist */
+ dhc->subblock_i[1] = dhc->n_lambda_vec;
+ for(i=0;i<dhc->n_lambda_vec;i++)
+ {
+ dhc->subblock_i[i+2] = dhc->native_lambda_components[i];
+ }
+ blk->sub[1].nr=2 + dhc->n_lambda_vec;
+ blk->sub[1].type=xdr_datatype_int;
+ blk->sub[1].ival=dhc->subblock_i;
+ }
+
for(i=0;i<dhc->ndh;i++)
{
nblock++;
#include "xvgr.h"
#include "gmx_ana.h"
#include "maths.h"
+#include "names.h"
+#include "mdebin.h"
/* Suppress Cygwin compiler warnings from using newlib version of
* ctype.h */
#undef isdigit
#endif
-/* the dhdl.xvg data from a simulation (actually obsolete, but still
- here for reading the dhdl.xvg file*/
+
+/* Structure for the names of lambda vector components */
+typedef struct lambda_components_t
+{
+ char **names; /* Array of strings with names for the lambda vector
+ components */
+ int N; /* The number of components */
+ int Nalloc; /* The number of allocated components */
+} lambda_components_t;
+
+/* Structure for a lambda vector or a dhdl derivative direction */
+typedef struct lambda_vec_t
+{
+ double *val; /* The lambda vector component values. Only valid if
+ dhdl == -1 */
+ int dhdl; /* The coordinate index for the derivative described by this
+ structure, or -1 */
+ const lambda_components_t *lc; /* the associated lambda_components
+ structure */
+ int index; /* The state number (init-lambda-state) of this lambda
+ vector, if known. If not, it is set to -1 */
+} lambda_vec_t;
+
+/* the dhdl.xvg data from a simulation */
typedef struct xvg_t
{
- char *filename;
+ const char *filename;
int ftp; /* file type */
int nset; /* number of lambdas, including dhdl */
int *np; /* number of data points (du or hists) per lambda */
int np_alloc; /* number of points (du or hists) allocated */
double temp; /* temperature */
- double *lambda; /* the lambdas (of first index for y). */
+ lambda_vec_t *lambda; /* the lambdas (of first index for y). */
double *t; /* the times (of second index for y) */
double **y; /* the dU values. y[0] holds the derivative, while
further ones contain the energy differences between
the native lambda and the 'foreign' lambdas. */
-
- double native_lambda; /* the native lambda */
+ lambda_vec_t native_lambda; /* the native lambda */
struct xvg_t *next, *prev; /*location in the global linked list of xvg_ts*/
} xvg_t;
-
typedef struct hist_t
{
unsigned int *bin[2]; /* the (forward + reverse) histogram values */
/* an aggregate of samples for partial free energy calculation */
typedef struct samples_t
{
- double native_lambda;
- double foreign_lambda;
+ lambda_vec_t *native_lambda; /* pointer to native lambda vector */
+ lambda_vec_t *foreign_lambda; /* pointer to foreign lambda vector */
double temp; /* the temperature */
gmx_bool derivative; /* whether this sample is a derivative */
foreign lambda) */
typedef struct sample_coll_t
{
- double native_lambda; /* these should be the same for all samples in the histogram?*/
- double foreign_lambda; /* collection */
+ lambda_vec_t *native_lambda; /* these should be the same for all samples
+ in the histogram */
+ lambda_vec_t *foreign_lambda; /* collection */
double temp; /* the temperature */
int nsamples; /* the number of samples */
} sample_coll_t;
/* all the samples associated with a lambda point */
-typedef struct lambda_t
+typedef struct lambda_data_t
{
- double lambda; /* the native lambda (at start time if dynamic) */
+ lambda_vec_t *lambda; /* the native lambda (at start time if dynamic) */
double temp; /* temperature */
sample_coll_t *sc; /* the samples */
sample_coll_t sc_head; /*the pre-allocated list head for the linked list.*/
- struct lambda_t *next, *prev; /* the next and prev in the list */
-} lambda_t;
+ struct lambda_data_t *next, *prev; /* the next and prev in the list */
+} lambda_data_t;
+
+/* Top-level data structure of simulation data */
+typedef struct sim_data_t
+{
+ lambda_data_t *lb; /* a lambda data linked list */
+ lambda_data_t lb_head; /* The head element of the linked list */
+ lambda_components_t lc; /* the allowed components of the lambda
+ vectors */
+} sim_data_t;
-/* calculated values. */
+/* Top-level data structure with calculated values. */
typedef struct {
sample_coll_t *a, *b; /* the simulation data */
} barres_t;
+/* Initialize a lambda_components structure */
+static void lambda_components_init(lambda_components_t *lc)
+{
+ lc->N=0;
+ lc->Nalloc=2;
+ snew(lc->names, lc->Nalloc);
+}
+
+/* Add a component to a lambda_components structure */
+static void lambda_components_add(lambda_components_t *lc,
+ const char *name, size_t name_length)
+{
+ while (lc->N + 1 > lc->Nalloc)
+ {
+ lc->Nalloc = (lc->Nalloc == 0) ? 2 : 2*lc->Nalloc;
+ srealloc( lc->names, lc->Nalloc );
+ }
+ snew(lc->names[lc->N], name_length+1);
+ strncpy(lc->names[lc->N], name, name_length);
+ lc->N++;
+}
+
+/* check whether a component with index 'index' matches the given name, or
+ is also NULL. Returns TRUE if this is the case.
+ the string name does not need to end */
+static gmx_bool lambda_components_check(const lambda_components_t *lc,
+ int index,
+ const char *name,
+ size_t name_length)
+{
+ size_t len;
+ if (index >= lc->N)
+ return FALSE;
+ if (name == NULL && lc->names[index] == NULL)
+ return TRUE;
+ if ((name == NULL) != (lc->names[index] == NULL))
+ return FALSE;
+ len = strlen(lc->names[index]);
+ if (len != name_length)
+ return FALSE;
+ if (strncmp(lc->names[index], name, name_length)== 0)
+ return TRUE;
+ return FALSE;
+}
+
+/* Find the index of a given lambda component name, or -1 if not found */
+static int lambda_components_find(const lambda_components_t *lc,
+ const char *name,
+ size_t name_length)
+{
+ int i;
+
+ for(i=0;i<lc->N;i++)
+ {
+ if (strncmp(lc->names[i], name, name_length) == 0)
+ return i;
+ }
+ return -1;
+}
+
+
+
+/* initialize a lambda vector */
+static void lambda_vec_init(lambda_vec_t *lv, const lambda_components_t *lc)
+{
+ snew(lv->val, lc->N);
+ lv->index = -1;
+ lv->dhdl = -1;
+ lv->lc = lc;
+}
+
+static void lambda_vec_destroy(lambda_vec_t *lv)
+{
+ sfree(lv->val);
+}
+
+static void lambda_vec_copy(lambda_vec_t *lv, const lambda_vec_t *orig)
+{
+ int i;
+
+ lambda_vec_init(lv, orig->lc);
+ lv->dhdl=orig->dhdl;
+ lv->index=orig->index;
+ for(i=0;i<lv->lc->N;i++)
+ {
+ lv->val[i] = orig->val[i];
+ }
+}
+
+/* write a lambda vec to a preallocated string */
+static void lambda_vec_print(const lambda_vec_t *lv, char *str, gmx_bool named)
+{
+ int i;
+ size_t np;
+
+ str[0]=0; /* reset the string */
+ if (lv -> dhdl < 0)
+ {
+ if (named)
+ {
+ str += sprintf(str, "delta H to ");
+ }
+ if (lv->lc->N > 1)
+ {
+ str += sprintf(str, "(");
+ }
+ for(i=0;i<lv->lc->N;i++)
+ {
+ str += sprintf(str, "%g", lv->val[i]);
+ if (i < lv->lc->N-1)
+ {
+ str += sprintf(str, ", ");
+ }
+ }
+ if (lv->lc->N > 1)
+ {
+ str += sprintf(str, ")");
+ }
+ }
+ else
+ {
+ /* this lambda vector describes a derivative */
+ str += sprintf(str, "dH/dl");
+ if (strlen(lv->lc->names[lv->dhdl]) > 0)
+ {
+ str += sprintf(str, " (%s)", lv->lc->names[lv->dhdl]);
+ }
+ }
+}
+
+/* write a shortened version of the lambda vec to a preallocated string */
+static void lambda_vec_print_short(const lambda_vec_t *lv, char *str)
+{
+ int i;
+ size_t np;
+
+ if (lv->index >= 0)
+ {
+ sprintf(str, "%6d", lv->index);
+ }
+ else
+ {
+ if (lv->dhdl < 0)
+ {
+ sprintf(str, "%6.3f", lv->val[0]);
+ }
+ else
+ {
+ sprintf(str, "dH/dl[%d]", lv->dhdl);
+ }
+ }
+}
+
+/* write an intermediate version of two lambda vecs to a preallocated string */
+static void lambda_vec_print_intermediate(const lambda_vec_t *a,
+ const lambda_vec_t *b, char *str)
+{
+ int i;
+ size_t np;
+
+ str[0]=0;
+ if ( (a->index >= 0) && (b->index>=0) )
+ {
+ sprintf(str, "%6.3f", ((double)a->index+(double)b->index)/2.);
+ }
+ else
+ {
+ if ( (a->dhdl < 0) && (b->dhdl < 0) )
+ {
+ sprintf(str, "%6.3f", (a->val[0]+b->val[0])/2.);
+ }
+ }
+}
+
+
+
+/* calculate the difference in lambda vectors: c = a-b.
+ c must be initialized already, and a and b must describe non-derivative
+ lambda points */
+static void lambda_vec_diff(const lambda_vec_t *a, const lambda_vec_t *b,
+ lambda_vec_t *c)
+{
+ int i;
+
+ if ( (a->dhdl > 0) || (b->dhdl > 0) )
+ {
+ gmx_fatal(FARGS,
+ "Trying to calculate the difference between derivatives instead of lambda points");
+ }
+ if ((a->lc != b->lc) || (a->lc != c->lc) )
+ {
+ gmx_fatal(FARGS,
+ "Trying to calculate the difference lambdas with differing basis set");
+ }
+ for(i=0;i<a->lc->N;i++)
+ {
+ c->val[i] = a->val[i] - b->val[i];
+ }
+}
+
+/* calculate and return the absolute difference in lambda vectors: c = |a-b|.
+ a and b must describe non-derivative lambda points */
+static double lambda_vec_abs_diff(const lambda_vec_t *a, const lambda_vec_t *b)
+{
+ int i;
+ double ret=0.;
+
+ if ( (a->dhdl > 0) || (b->dhdl > 0) )
+ {
+ gmx_fatal(FARGS,
+ "Trying to calculate the difference between derivatives instead of lambda points");
+ }
+ if (a->lc != b->lc)
+ {
+ gmx_fatal(FARGS,
+ "Trying to calculate the difference lambdas with differing basis set");
+ }
+ for(i=0;i<a->lc->N;i++)
+ {
+ double df=a->val[i] - b->val[i];
+ ret += df*df;
+ }
+ return sqrt(ret);
+}
+
+
+/* check whether two lambda vectors are the same */
+static gmx_bool lambda_vec_same(const lambda_vec_t *a, const lambda_vec_t *b)
+{
+ int i;
+
+ if (a->lc != b->lc)
+ return FALSE;
+ if (a->dhdl < 0)
+ {
+ for(i=0;i<a->lc->N;i++)
+ {
+ if (!gmx_within_tol(a->val[i], b->val[i], 10*GMX_REAL_EPS))
+ {
+ return FALSE;
+ }
+ }
+ return TRUE;
+ }
+ else
+ {
+ /* they're derivatives, so we check whether the indices match */
+ return (a->dhdl == b->dhdl);
+ }
+}
+
+/* Compare the sort order of two foreign lambda vectors
+
+ returns 1 if a is 'bigger' than b,
+ returns 0 if they're the same,
+ returns -1 if a is 'smaller' than b.*/
+static gmx_bool lambda_vec_cmp_foreign(const lambda_vec_t *a,
+ const lambda_vec_t *b)
+{
+ int i;
+ double norm_a=0, norm_b=0;
+ gmx_bool different=FALSE;
+
+ if (a->lc != b->lc)
+ {
+ gmx_fatal(FARGS, "Can't compare lambdas with differing basis sets");
+ }
+ /* if either one has an index we sort based on that */
+ if ((a->index >= 0) || (b->index >= 0))
+ {
+ if (a->index == b->index)
+ return 0;
+ return (a->index > b->index) ? 1 : -1;
+ }
+ if (a->dhdl >= 0 || b->dhdl >= 0)
+ {
+ /* lambda vectors that are derivatives always sort higher than those
+ without derivatives */
+ if ((a->dhdl >= 0) != (b->dhdl >= 0) )
+ {
+ return (a->dhdl >= 0) ? 1 : -1 ;
+ }
+ return a->dhdl > b->dhdl;
+ }
+
+ /* neither has an index, so we can only sort on the lambda components,
+ which is only valid if there is one component */
+ for(i=0;i<a->lc->N;i++)
+ {
+ if (!gmx_within_tol(a->val[i], b->val[i], 10*GMX_REAL_EPS))
+ {
+ different=TRUE;
+ }
+ norm_a += a->val[i]*a->val[i];
+ norm_b += b->val[i]*b->val[i];
+ }
+ if (!different)
+ {
+ return 0;
+ }
+ return norm_a > norm_b;
+}
+
+/* Compare the sort order of two native lambda vectors
+
+ returns 1 if a is 'bigger' than b,
+ returns 0 if they're the same,
+ returns -1 if a is 'smaller' than b.*/
+static gmx_bool lambda_vec_cmp_native(const lambda_vec_t *a,
+ const lambda_vec_t *b)
+{
+ int i;
+
+ if (a->lc != b->lc)
+ {
+ gmx_fatal(FARGS, "Can't compare lambdas with differing basis sets");
+ }
+ /* if either one has an index we sort based on that */
+ if ((a->index >= 0) || (b->index >= 0))
+ {
+ if (a->index == b->index)
+ return 0;
+ return (a->index > b->index) ? 1 : -1;
+ }
+ /* neither has an index, so we can only sort on the lambda components,
+ which is only valid if there is one component */
+ if (a->lc->N > 1)
+ {
+ gmx_fatal(FARGS,
+ "Can't compare lambdas with no index and > 1 component");
+ }
+ if (a->dhdl >= 0 || b->dhdl >= 0)
+ {
+ gmx_fatal(FARGS,
+ "Can't compare native lambdas that are derivatives");
+ }
+ if (gmx_within_tol(a->val[0], b->val[0], 10*GMX_REAL_EPS))
+ {
+ return 0;
+ }
+ return a->val[0] > b->val[0] ? 1 : -1;
+}
+
+
static void hist_init(hist_t *h, int nhist, int *nbin)
ba->y=NULL;
}
-static void samples_init(samples_t *s, double native_lambda,
- double foreign_lambda, double temp,
+static void samples_init(samples_t *s, lambda_vec_t *native_lambda,
+ lambda_vec_t *foreign_lambda, double temp,
gmx_bool derivative, const char *filename)
{
s->native_lambda=native_lambda;
r->s=NULL;
}
-static void sample_coll_init(sample_coll_t *sc, double native_lambda,
- double foreign_lambda, double temp)
+static void sample_coll_init(sample_coll_t *sc, lambda_vec_t *native_lambda,
+ lambda_vec_t *foreign_lambda, double temp)
{
sc->native_lambda = native_lambda;
sc->foreign_lambda = foreign_lambda;
}
-static void lambda_init(lambda_t *l, double native_lambda, double temp)
+static void lambda_data_init(lambda_data_t *l, lambda_vec_t *native_lambda,
+ double temp)
{
l->lambda=native_lambda;
l->temp=temp;
l->sc=&(l->sc_head);
- sample_coll_init(l->sc, native_lambda, 0., 0.);
+ sample_coll_init(l->sc, native_lambda, NULL, 0.);
l->sc->next=l->sc;
l->sc->prev=l->sc;
}
}
-
-
-static gmx_bool lambda_same(double lambda1, double lambda2)
-{
- return gmx_within_tol(lambda1, lambda2, 10*GMX_REAL_EPS);
-}
-
/* calculate the total number of samples in a sample collection */
static void sample_coll_calc_ntot(sample_coll_t *sc)
{
/* find the barsamples_t associated with a lambda that corresponds to
a specific foreign lambda */
-static sample_coll_t *lambda_find_sample_coll(lambda_t *l,
- double foreign_lambda)
+static sample_coll_t *lambda_data_find_sample_coll(lambda_data_t *l,
+ lambda_vec_t *foreign_lambda)
{
sample_coll_t *sc=l->sc->next;
while(sc != l->sc)
{
- if (lambda_same(sc->foreign_lambda,foreign_lambda))
+ if (lambda_vec_same(sc->foreign_lambda,foreign_lambda))
{
return sc;
}
}
/* insert li into an ordered list of lambda_colls */
-static void lambda_insert_sample_coll(lambda_t *l, sample_coll_t *sc)
+static void lambda_data_insert_sample_coll(lambda_data_t *l, sample_coll_t *sc)
{
sample_coll_t *scn=l->sc->next;
while ( (scn!=l->sc) )
{
- if (scn->foreign_lambda > sc->foreign_lambda)
+ if (lambda_vec_cmp_foreign(scn->foreign_lambda, sc->foreign_lambda) > 0)
break;
scn=scn->next;
}
}
/* insert li into an ordered list of lambdas */
-static void lambda_insert_lambda(lambda_t *head, lambda_t *li)
+static void lambda_data_insert_lambda(lambda_data_t *head, lambda_data_t *li)
{
- lambda_t *lc=head->next;
+ lambda_data_t *lc=head->next;
while (lc!=head)
{
- if (lc->lambda > li->lambda)
+ if (lambda_vec_cmp_native(lc->lambda, li->lambda) > 0)
break;
lc=lc->next;
}
gmx_fatal(FARGS, "Temperatures in files %s and %s are not the same!",
s->filename, sc->next->s[0]->filename);
}
- if (sc->native_lambda != s->native_lambda)
+ if (!lambda_vec_same(sc->native_lambda, s->native_lambda))
{
gmx_fatal(FARGS, "Native lambda in files %s and %s are not the same (and they should be)!",
s->filename, sc->next->s[0]->filename);
}
- if (sc->foreign_lambda != s->foreign_lambda)
+ if (!lambda_vec_same(sc->foreign_lambda,s->foreign_lambda))
{
gmx_fatal(FARGS, "Foreign lambda in files %s and %s are not the same (and they should be)!",
s->filename, sc->next->s[0]->filename);
/* insert a sample into a lambda_list, creating the right sample_coll if
neccesary */
-static void lambda_list_insert_sample(lambda_t *head, samples_t *s)
+static void lambda_data_list_insert_sample(lambda_data_t *head, samples_t *s)
{
gmx_bool found=FALSE;
sample_coll_t *sc;
sample_range_t r;
- lambda_t *l=head->next;
+ lambda_data_t *l=head->next;
- /* first search for the right lambda_t */
+ /* first search for the right lambda_data_t */
while(l != head)
{
- if (lambda_same(l->lambda, s->native_lambda) )
+ if (lambda_vec_same(l->lambda, s->native_lambda) )
{
found=TRUE;
break;
if (!found)
{
snew(l, 1); /* allocate a new one */
- lambda_init(l, s->native_lambda, s->temp); /* initialize it */
- lambda_insert_lambda(head, l); /* add it to the list */
+ lambda_data_init(l, s->native_lambda, s->temp); /* initialize it */
+ lambda_data_insert_lambda(head, l); /* add it to the list */
}
/* now look for a sample collection */
- sc=lambda_find_sample_coll(l, s->foreign_lambda);
+ sc=lambda_data_find_sample_coll(l, s->foreign_lambda);
if (!sc)
{
snew(sc, 1); /* allocate a new one */
sample_coll_init(sc, s->native_lambda, s->foreign_lambda, s->temp);
- lambda_insert_sample_coll(l, sc);
+ lambda_data_insert_sample_coll(l, sc);
}
/* now insert the samples into the sample coll */
}
/* write a collection of histograms to a file */
-void lambdas_histogram(lambda_t *bl_head, const char *filename,
- int nbin_default, const output_env_t oenv)
+void sim_data_histogram(sim_data_t *sd, const char *filename,
+ int nbin_default, const output_env_t oenv)
{
char label_x[STRLEN];
const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda";
const char *title="N(\\DeltaH)";
const char *label_y="Samples";
FILE *fp;
- lambda_t *bl;
+ lambda_data_t *bl;
int nsets=0;
char **setnames=NULL;
gmx_bool first_set=FALSE;
double dx=0;
double min=0;
int i;
+ lambda_data_t *bl_head = sd->lb;
printf("\nWriting histogram to %s\n", filename);
sprintf(label_x, "\\DeltaH (%s)", unit_energy);
/* iterate over all samples */
while(sc!=bl->sc)
{
+ char buf[STRLEN], buf2[STRLEN];
+
nsets++;
srenew(setnames, nsets);
snew(setnames[nsets-1], STRLEN);
- if (!lambda_same(sc->foreign_lambda, sc->native_lambda))
+ if (sc->foreign_lambda->dhdl < 0)
{
- sprintf(setnames[nsets-1], "N(%s(%s=%g) | %s=%g)",
- deltag, lambda, sc->foreign_lambda, lambda,
- sc->native_lambda);
+ lambda_vec_print(sc->native_lambda, buf, FALSE);
+ lambda_vec_print(sc->foreign_lambda, buf2, FALSE);
+ sprintf(setnames[nsets-1], "N(%s(%s=%s) | %s=%s)",
+ deltag, lambda, buf2, lambda, buf);
}
else
{
- sprintf(setnames[nsets-1], "N(%s | %s=%g)",
- dhdl, lambda, sc->native_lambda);
+ lambda_vec_print(sc->native_lambda, buf, FALSE);
+ sprintf(setnames[nsets-1], "N(%s | %s=%s)",
+ dhdl, lambda, buf);
}
sc=sc->next;
}
/* create a collection (array) of barres_t object given a ordered linked list
of barlamda_t sample collections */
-static barres_t *barres_list_create(lambda_t *bl_head, int *nres,
+static barres_t *barres_list_create(sim_data_t *sd, int *nres,
gmx_bool use_dhdl)
{
- lambda_t *bl;
+ lambda_data_t *bl;
int nlambda=0;
barres_t *res;
int i;
gmx_bool dhdl=FALSE;
gmx_bool first=TRUE;
+ lambda_data_t *bl_head=sd->lb;
/* first count the lambdas */
bl=bl_head->next;
barres_t *br=&(res[*nres]);
/* there is always a previous one. we search for that as a foreign
lambda: */
- scprev=lambda_find_sample_coll(bl->prev, bl->lambda);
- sc=lambda_find_sample_coll(bl, bl->prev->lambda);
+ scprev=lambda_data_find_sample_coll(bl->prev, bl->lambda);
+ sc=lambda_data_find_sample_coll(bl, bl->prev->lambda);
barres_init(br);
{
/* we use dhdl */
- scprev=lambda_find_sample_coll(bl->prev, bl->prev->lambda);
- sc=lambda_find_sample_coll(bl, bl->lambda);
+ scprev=lambda_data_find_sample_coll(bl->prev, bl->prev->lambda);
+ sc=lambda_data_find_sample_coll(bl, bl->lambda);
if (first)
{
{
barres_t *br=&(res[i]);
- delta_lambda=fabs(br->b->native_lambda-br->a->native_lambda);
+ delta_lambda=lambda_vec_abs_diff(br->b->native_lambda,
+ br->a->native_lambda);
for(j=0;j<br->a->nsamples;j++)
{
sample_coll_calc_ntot(sc);
}
-static void lambdas_impose_times(lambda_t *head, double begin, double end)
+static void sim_data_impose_times(sim_data_t *sd, double begin, double end)
{
double first_t, last_t;
double begin_t, end_t;
- lambda_t *lc;
+ lambda_data_t *lc;
+ lambda_data_t *head=sd->lb;
int j;
if (begin<=0 && end<0)
}
}
+/* Initialize a sim_data structure */
+static void sim_data_init(sim_data_t *sd)
+{
+ /* make linked list */
+ sd->lb = &(sd->lb_head);
+ sd->lb->next = sd->lb;
+ sd->lb->prev = sd->lb;
+
+ lambda_components_init(&(sd->lc));
+}
+
static double calc_bar_sum(int n,const double *W,double Wfac,double sbMmDG)
{
M = log(n1/n2);
- if (!lambda_same(ca->native_lambda, ca->foreign_lambda))
+ /*if (!lambda_vec_same(ca->native_lambda, ca->foreign_lambda))*/
+ if (ca->foreign_lambda->dhdl<0)
{
/* this is the case when the delta U were calculated directly
(i.e. we're not scaling dhdl) */
{
/* we're using dhdl, so delta_lambda needs to be a
multiplication factor. */
- double delta_lambda=cb->native_lambda-ca->native_lambda;
+ /*double delta_lambda=cb->native_lambda-ca->native_lambda;*/
+ double delta_lambda=lambda_vec_abs_diff(cb->native_lambda,
+ ca->native_lambda);
+ if (cb->native_lambda->lc->N > 1)
+ {
+ gmx_fatal(FARGS,
+ "Can't (yet) do multi-component dhdl interpolation");
+ }
+
Wfac1 = beta*delta_lambda;
Wfac2 = -beta*delta_lambda;
}
n2 = cb->ntot;
/* to ensure the work values are the same as during the delta_G */
- if (!lambda_same(ca->native_lambda, ca->foreign_lambda))
+ /*if (!lambda_vec_same(ca->native_lambda, ca->foreign_lambda))*/
+ if (ca->foreign_lambda->dhdl<0)
{
/* this is the case when the delta U were calculated directly
(i.e. we're not scaling dhdl) */
{
/* we're using dhdl, so delta_lambda needs to be a
multiplication factor. */
- double delta_lambda=cb->native_lambda-ca->native_lambda;
+ double delta_lambda=lambda_vec_abs_diff(cb->native_lambda,
+ ca->native_lambda);
Wfac1 = beta*delta_lambda;
Wfac2 = -beta*delta_lambda;
}
n2 = cb->ntot;
/* to ensure the work values are the same as during the delta_G */
- if (!lambda_same(ca->native_lambda, ca->foreign_lambda))
+ /*if (!lambda_vec_same(ca->native_lambda, ca->foreign_lambda))*/
+ if (ca->foreign_lambda->dhdl<0)
{
/* this is the case when the delta U were calculated directly
(i.e. we're not scaling dhdl) */
{
/* we're using dhdl, so delta_lambda needs to be a
multiplication factor. */
- double delta_lambda=cb->native_lambda-ca->native_lambda;
+ double delta_lambda=lambda_vec_abs_diff(cb->native_lambda,
+ ca->native_lambda);
Wfac1 = beta*delta_lambda;
Wfac2 = -beta*delta_lambda;
}
return sqrt(svar/(nbmax + 1 - nbmin));
}
+
+/* Seek the end of an identifier (consecutive non-spaces), followed by
+ an optional number of spaces or '='-signs. Returns a pointer to the
+ first non-space value found after that. Returns NULL if the string
+ ends before that.
+ */
+static const char *find_value(const char *str)
+{
+ gmx_bool name_end_found=FALSE;
+
+ /* if the string is a NULL pointer, return a NULL pointer. */
+ if (str==NULL)
+ {
+ return NULL;
+ }
+ while(*str != '\0')
+ {
+ /* first find the end of the name */
+ if (! name_end_found)
+ {
+ if ( isspace(*str) || (*str == '=') )
+ {
+ name_end_found=TRUE;
+ }
+ }
+ else
+ {
+ if (! ( isspace(*str) || (*str == '=') ))
+ {
+ return str;
+ }
+ }
+ str++;
+ }
+ return NULL;
+}
+
+
+
+/* read a vector-notation description of a lambda vector */
+static gmx_bool read_lambda_compvec(const char *str,
+ lambda_vec_t *lv,
+ const lambda_components_t *lc_in,
+ lambda_components_t *lc_out,
+ const char **end,
+ const char *fn)
+{
+ gmx_bool initialize_lc = FALSE; /* whether to initialize the lambda
+ components, or to check them */
+ gmx_bool start_reached=FALSE; /* whether the start of component names
+ has been reached */
+ gmx_bool vector=FALSE; /* whether there are multiple components */
+ int n = 0; /* current component number */
+ const char *val_start=NULL; /* start of the component name, or NULL
+ if not in a value */
+ char *strtod_end;
+ gmx_bool OK=TRUE;
+
+ if (end)
+ *end=str;
+
+
+ if (lc_out && lc_out->N == 0)
+ {
+ initialize_lc = TRUE;
+ }
+
+ if (lc_in == NULL)
+ lc_in = lc_out;
+
+ while(1)
+ {
+ if (!start_reached)
+ {
+ if (isalnum(*str))
+ {
+ vector=FALSE;
+ start_reached=TRUE;
+ val_start=str;
+ }
+ else if (*str=='(')
+ {
+ vector=TRUE;
+ start_reached=TRUE;
+ }
+ else if (! isspace(*str))
+ {
+ gmx_fatal(FARGS, "Error in lambda components in %s", fn);
+ }
+ }
+ else
+ {
+ if (val_start)
+ {
+ if (isspace(*str) || *str==')' || *str==',' || *str=='\0')
+ {
+ /* end of value */
+ if (lv == NULL)
+ {
+ if (initialize_lc)
+ {
+ lambda_components_add(lc_out, val_start,
+ (str-val_start));
+ }
+ else
+ {
+ if (!lambda_components_check(lc_out, n, val_start,
+ (str-val_start)))
+ {
+ return FALSE;
+ }
+ }
+ }
+ else
+ {
+ /* add a vector component to lv */
+ lv->val[n] = strtod(val_start, &strtod_end);
+ if (val_start == strtod_end)
+ {
+ gmx_fatal(FARGS,
+ "Error reading lambda vector in %s", fn);
+ }
+ }
+ /* reset for the next identifier */
+ val_start=NULL;
+ n++;
+ if (!vector)
+ {
+ return OK;
+ }
+ }
+ }
+ else if (isalnum(*str))
+ {
+ val_start=str;
+ }
+ if (*str == ')')
+ {
+ str++;
+ if (end)
+ *end=str;
+ if (!vector)
+ {
+ gmx_fatal(FARGS, "Error in lambda components in %s", fn);
+ }
+ else
+ {
+ if (n == lc_in->N)
+ {
+ return OK;
+ }
+ else if (lv == NULL)
+ {
+ return FALSE;
+ }
+ else
+ {
+ gmx_fatal(FARGS, "Incomplete lambda vector data in %s",
+ fn);
+ return FALSE;
+ }
+
+ }
+ }
+ }
+ if (*str == '\0')
+ {
+ break;
+ }
+ str++;
+ if (end)
+ {
+ *end=str;
+ }
+ }
+ if (vector)
+ {
+ gmx_fatal(FARGS, "Incomplete lambda components data in %s", fn);
+ return FALSE;
+ }
+ return OK;
+}
+
+/* read and check the component names from a string */
+static gmx_bool read_lambda_components(const char *str,
+ lambda_components_t *lc,
+ const char **end,
+ const char *fn)
+{
+ return read_lambda_compvec(str, NULL, NULL, lc, end, fn);
+}
+
+/* read an initialized lambda vector from a string */
+static gmx_bool read_lambda_vector(const char *str,
+ lambda_vec_t *lv,
+ const char **end,
+ const char *fn)
+{
+ return read_lambda_compvec(str, lv, lv->lc, NULL, end, fn);
+}
+
+
+
/* deduce lambda value from legend.
-input:
- bdhdl = if true, value may be a derivative.
-output:
- bdhdl = whether the legend was for a derivative.
+ fn = the file name
+ legend = the legend string
+ ba = the xvg data
+ lam = the initialized lambda vector
+returns whether to use the data in this set.
*/
-static double legend2lambda(char *fn,const char *legend,gmx_bool *bdhdl)
+static gmx_bool legend2lambda(const char *fn,
+ const char *legend,
+ xvg_t *ba,
+ lambda_vec_t *lam)
{
double lambda=0;
- const char *ptr;
+ const char *ptr=NULL, *ptr2=NULL;
gmx_bool ok=FALSE;
+ gmx_bool bdhdl=FALSE;
+ const char *tostr=" to ";
if (legend == NULL)
{
gmx_fatal(FARGS,"There is no legend in file '%s', can not deduce lambda",fn);
}
- ptr = strrchr(legend,' ');
- if (strstr(legend,"dH"))
+ /* look for the last 'to': */
+ ptr2=legend;
+ do
{
- if (! (*bdhdl))
+ ptr2 = strstr(ptr2, tostr);
+ if (ptr2!=NULL)
{
- ok=FALSE;
- }
- else
- {
- ok=TRUE;
+ ptr=ptr2;
+ ptr2++;
}
}
+ while(ptr2!=NULL && *ptr2!= '\0' );
+
+ if (ptr)
+ {
+ ptr += strlen(tostr)-1; /* and advance past that 'to' */
+ }
else
{
- if (strchr(legend,'D') != NULL && strchr(legend,'H') != NULL)
+ /* look for the = sign */
+ ptr = strrchr(legend,'=');
+ if (!ptr)
{
- ok=TRUE;
- *bdhdl=FALSE;
+ /* otherwise look for the last space */
+ ptr = strrchr(legend,' ');
}
}
+
+ if (strstr(legend,"dH"))
+ {
+ ok=TRUE;
+ bdhdl=TRUE;
+ }
+ else if (strchr(legend,'D') != NULL && strchr(legend,'H') != NULL)
+ {
+ ok=TRUE;
+ bdhdl=FALSE;
+ }
+ else /*if (strstr(legend, "pV"))*/
+ {
+ return FALSE;
+ }
if (!ptr)
{
ok=FALSE;
if (!ok)
{
- printf("%s\n", legend);
gmx_fatal(FARGS,"There is no proper lambda legend in file '%s', can not deduce lambda",fn);
}
- if (sscanf(ptr,"%lf",&lambda) != 1)
+ if (!bdhdl)
{
- gmx_fatal(FARGS,"There is no proper lambda legend in file '%s', can not deduce lambda",fn);
+ ptr=find_value(ptr);
+ if (!ptr || !read_lambda_vector(ptr, lam, NULL, fn))
+ {
+ gmx_fatal(FARGS, "lambda vector '%s' %s faulty", legend, fn);
+ }
}
+ else
+ {
+ int dhdl_index;
+ const char *end;
+ char buf[STRLEN];
- return lambda;
+ ptr = strrchr(legend, '=');
+ end=ptr;
+ if (ptr)
+ {
+ /* there must be a component name */
+ ptr--;
+ if (ptr < legend)
+ {
+ gmx_fatal(FARGS, "dhdl legend '%s' %s faulty", legend, fn);
+ }
+ /* now backtrack to the start of the identifier */
+ while(isspace(*ptr))
+ {
+ end=ptr;
+ ptr--;
+ if (ptr < legend)
+ {
+ gmx_fatal(FARGS, "dhdl legend '%s' %s faulty", legend, fn);
+ }
+ }
+ while(!isspace(*ptr))
+ {
+ ptr--;
+ if (ptr < legend)
+ {
+ gmx_fatal(FARGS, "dhdl legend '%s' %s faulty", legend, fn);
+ }
+ }
+ ptr++;
+ strncpy(buf, ptr, (end-ptr));
+ buf[(end-ptr)]='\0';
+ dhdl_index=lambda_components_find(lam->lc, ptr, (end-ptr));
+ if (dhdl_index < 0)
+ {
+ char buf[STRLEN];
+ strncpy(buf, ptr, (end-ptr));
+ buf[(end-ptr)]='\0';
+ gmx_fatal(FARGS,
+ "Did not find lambda component for '%s' in %s",
+ buf, fn);
+ }
+ }
+ else
+ {
+ if (lam->lc->N > 1)
+ {
+ gmx_fatal(FARGS,
+ "dhdl without component name with >1 lambda component in %s",
+ fn);
+ }
+ dhdl_index=0;
+ }
+ lam->dhdl = dhdl_index;
+ }
+ return TRUE;
}
-static gmx_bool subtitle2lambda(const char *subtitle,double *lambda)
+static gmx_bool subtitle2lambda(const char *subtitle, xvg_t *ba, const char *fn,
+ lambda_components_t *lc)
{
gmx_bool bFound;
- char *ptr;
+ const char *ptr;
+ char *end;
+ double native_lambda;
bFound = FALSE;
- /* plain text lambda string */
- ptr = strstr(subtitle,"lambda");
- if (ptr == NULL)
- {
- /* xmgrace formatted lambda string */
- ptr = strstr(subtitle,"\\xl\\f{}");
- }
- if (ptr == NULL)
- {
- /* xmgr formatted lambda string */
- ptr = strstr(subtitle,"\\8l\\4");
- }
- if (ptr != NULL)
+ /* first check for a state string */
+ ptr = strstr(subtitle, "state");
+ if (ptr)
{
- ptr = strstr(ptr,"=");
+ int index=-1;
+ const char *val_end;
+
+ /* the new 4.6 style lambda vectors */
+ ptr = find_value(ptr);
+ if (ptr)
+ {
+ index = strtol(ptr, &end, 10);
+ if (ptr == end)
+ {
+ gmx_fatal(FARGS,"Incomplete state data in %s", fn);
+ return FALSE;
+ }
+ ptr=end;
+ }
+ else
+ {
+ gmx_fatal(FARGS,"Incomplete state data in %s", fn);
+ return FALSE;
+ }
+ /* now find the lambda vector component names */
+ while(*ptr!='(' && ! isalnum(*ptr))
+ {
+ ptr++;
+ if (*ptr == '\0')
+ {
+ gmx_fatal(FARGS,
+ "Incomplete lambda vector component data in %s", fn);
+ return FALSE;
+ }
+ }
+ val_end=ptr;
+ if (!read_lambda_components(ptr, lc, &val_end, fn))
+ {
+ gmx_fatal(FARGS,
+ "lambda vector components in %s don't match those previously read",
+ fn);
+ }
+ ptr=find_value(val_end);
+ if (!ptr)
+ {
+ gmx_fatal(FARGS,"Incomplete state data in %s", fn);
+ return FALSE;
+ }
+ lambda_vec_init(&(ba->native_lambda), lc);
+ if (!read_lambda_vector(ptr, &(ba->native_lambda), NULL, fn))
+ {
+ gmx_fatal(FARGS, "lambda vector in %s faulty", fn);
+ }
+ ba->native_lambda.index=index;
+ bFound=TRUE;
}
- if (ptr != NULL)
+ else
{
- bFound = (sscanf(ptr+1,"%lf",lambda) == 1);
+ /* compatibility mode: check for lambda in other ways. */
+ /* plain text lambda string */
+ ptr = strstr(subtitle,"lambda");
+ if (ptr == NULL)
+ {
+ /* xmgrace formatted lambda string */
+ ptr = strstr(subtitle,"\\xl\\f{}");
+ }
+ if (ptr == NULL)
+ {
+ /* xmgr formatted lambda string */
+ ptr = strstr(subtitle,"\\8l\\4");
+ }
+ if (ptr != NULL)
+ {
+ ptr = strstr(ptr,"=");
+ }
+ if (ptr != NULL)
+ {
+ bFound = (sscanf(ptr+1,"%lf",&(native_lambda)) == 1);
+ /* add the lambda component name as an empty string */
+ if (lc->N > 0)
+ {
+ if (!lambda_components_check(lc, 0, "", 0))
+ {
+ gmx_fatal(FARGS,
+ "lambda vector components in %s don't match those previously read",
+ fn);
+ }
+ }
+ else
+ {
+ lambda_components_add(lc, "", 0);
+ }
+ lambda_vec_init(&(ba->native_lambda), lc);
+ ba->native_lambda.val[0]=native_lambda;
+ }
}
return bFound;
}
-static double filename2lambda(char *fn)
+static void filename2lambda(const char *fn, xvg_t *ba)
{
double lambda;
- char *ptr,*endptr,*digitptr;
+ const char *ptr,*digitptr;
+ char *endptr;
int dirsep;
ptr = fn;
/* go to the end of the path string and search backward to find the last
{
gmx_fatal(FARGS,"Malformed number in file path '%s'",fn);
}
-
- return lambda;
}
-static void read_bar_xvg_lowlevel(char *fn, real *temp, xvg_t *ba)
+static void read_bar_xvg_lowlevel(const char *fn, real *temp, xvg_t *ba,
+ lambda_components_t *lc)
{
int i;
char *subtitle,**legend,*ptr;
int np;
gmx_bool native_lambda_read=FALSE;
+ char buf[STRLEN];
+ lambda_vec_t lv;
xvg_init(ba);
{
gmx_fatal(FARGS,"File %s contains no usable data.",fn);
}
+ /* Reorder the data */
ba->t = ba->y[0];
+ for(i=1; i<ba->nset; i++)
+ {
+ ba->y[i-1] = ba->y[i];
+ }
+ ba->nset--;
snew(ba->np,ba->nset);
for(i=0;i<ba->nset;i++)
ba->np[i]=np;
-
ba->temp = -1;
if (subtitle != NULL)
{
+ /* try to extract temperature */
ptr = strstr(subtitle,"T =");
if (ptr != NULL)
{
/* Try to deduce lambda from the subtitle */
if (subtitle)
{
- if (subtitle2lambda(subtitle,&(ba->native_lambda)))
+ if (subtitle2lambda(subtitle,ba, fn, lc))
{
native_lambda_read=TRUE;
}
}
- snew(ba->lambda,ba->nset-1);
+ snew(ba->lambda,ba->nset);
if (legend == NULL)
{
- /* Check if we have a single set, no legend, nset=2 means t and dH/dl */
- if (ba->nset == 2)
+ /* Check if we have a single set, no legend, nset=1 means t and dH/dl */
+ if (ba->nset == 1)
{
if (!native_lambda_read)
{
/* Deduce lambda from the file name */
- ba->native_lambda = filename2lambda(fn);
+ filename2lambda(fn, ba);
native_lambda_read=TRUE;
}
ba->lambda[0] = ba->native_lambda;
}
else
{
- for(i=0; i<ba->nset-1; i++)
+ for(i=0; i<ba->nset; )
{
- gmx_bool is_dhdl=(i==0);
+ gmx_bool use=FALSE;
/* Read lambda from the legend */
- ba->lambda[i] = legend2lambda(fn,legend[i], &is_dhdl);
-
- if (is_dhdl && !native_lambda_read)
+ lambda_vec_init( &(ba->lambda[i]), lc );
+ lambda_vec_copy( &(ba->lambda[i]), &(ba->native_lambda));
+ use=legend2lambda(fn,legend[i], ba, &(ba->lambda[i]));
+ if (use)
{
- ba->native_lambda = ba->lambda[i];
- native_lambda_read=TRUE;
+ lambda_vec_print(&(ba->lambda[i]), buf, FALSE);
+ i++;
+ }
+ else
+ {
+ int j;
+ printf("%s: Ignoring set '%s'.\n", fn, legend[i]);
+ for(j=i+1; j<ba->nset; j++)
+ {
+ ba->y[j-1] = ba->y[j];
+ legend[j-1] = legend[j];
+ }
+ ba->nset--;
}
}
}
gmx_fatal(FARGS,"File %s contains multiple sets but no indication of the native lambda",fn);
}
- /* Reorder the data */
- for(i=1; i<ba->nset; i++)
- {
- ba->y[i-1] = ba->y[i];
- }
if (legend != NULL)
{
for(i=0; i<ba->nset-1; i++)
}
sfree(legend);
}
- ba->nset--;
}
-static void read_bar_xvg(char *fn, real *temp, lambda_t *lambda_head)
+static void read_bar_xvg(char *fn, real *temp, sim_data_t *sd)
{
xvg_t *barsim;
samples_t *s;
snew(barsim,1);
- read_bar_xvg_lowlevel(fn, temp, barsim);
+ read_bar_xvg_lowlevel(fn, temp, barsim, &(sd->lc));
if (barsim->nset <1 )
{
snew(s, barsim->nset);
for(i=0;i<barsim->nset;i++)
{
- samples_init(s+i, barsim->native_lambda, barsim->lambda[i],
- barsim->temp, lambda_same(barsim->native_lambda,
- barsim->lambda[i]),
+ samples_init(s+i, &(barsim->native_lambda), &(barsim->lambda[i]),
+ barsim->temp, lambda_vec_same(&(barsim->native_lambda),
+ &(barsim->lambda[i])),
fn);
s[i].du=barsim->y[i];
s[i].ndu=barsim->np[i];
s[i].t=barsim->t;
- lambda_list_insert_sample(lambda_head, s+i);
+ lambda_data_list_insert_sample(sd->lb, s+i);
}
- printf("%s: %.1f - %.1f; lambda = %.3f\n foreign lambdas:",
- fn, s[0].t[0], s[0].t[s[0].ndu-1], s[0].native_lambda);
- for(i=0;i<barsim->nset;i++)
{
- printf(" %.3f (%d pts)", s[i].foreign_lambda, s[i].ndu);
+ char buf[STRLEN];
+
+ lambda_vec_print(s[0].native_lambda, buf, FALSE);
+ printf("%s: %.1f - %.1f; lambda = %s\n dH/dl & foreign lambdas:\n",
+ fn, s[0].t[0], s[0].t[s[0].ndu-1], buf);
+ for(i=0;i<barsim->nset;i++)
+ {
+ lambda_vec_print(s[i].foreign_lambda, buf, TRUE);
+ printf(" %s (%d pts)\n", buf, s[i].ndu);
+ }
}
printf("\n\n");
}
static void read_edr_rawdh_block(samples_t **smp, int *ndu, t_enxblock *blk,
double start_time, double delta_time,
- double native_lambda, double temp,
+ lambda_vec_t *native_lambda, double temp,
double *last_t, const char *filename)
{
- int j;
+ int i,j;
gmx_bool allocated;
- double foreign_lambda;
- int derivative;
+ double old_foreign_lambda;
+ lambda_vec_t *foreign_lambda;
+ int type;
samples_t *s; /* convenience pointer */
int startj;
"Unexpected/corrupted block data in file %s around time %g.",
filename, start_time);
}
-
- derivative = blk->sub[0].ival[0];
- foreign_lambda = blk->sub[1].dval[0];
+
+ snew(foreign_lambda, 1);
+ lambda_vec_init(foreign_lambda, native_lambda->lc);
+ lambda_vec_copy(foreign_lambda, native_lambda);
+ type = blk->sub[0].ival[0];
+ if (type == dhbtDH)
+ {
+ for(i=0;i<native_lambda->lc->N;i++)
+ {
+ foreign_lambda->val[i] = blk->sub[1].dval[i];
+ }
+ }
+ else
+ {
+ if (blk->sub[0].nr > 1)
+ {
+ foreign_lambda->dhdl = blk->sub[0].ival[1];
+ }
+ else
+ foreign_lambda->dhdl = 0;
+ }
if (! *smp)
{
/* initialize the samples structure if it's empty. */
snew(*smp, 1);
samples_init(*smp, native_lambda, foreign_lambda, temp,
- derivative!=0, filename);
+ type==dhbtDHDL, filename);
(*smp)->start_time=start_time;
(*smp)->delta_time=delta_time;
}
s=*smp;
/* now double check */
- if ( ! lambda_same(s->foreign_lambda, foreign_lambda) ||
- ( (derivative!=0) != (s->derivative!=0) ) )
+ if ( ! lambda_vec_same(s->foreign_lambda, foreign_lambda) )
{
- fprintf(stderr, "Got foreign lambda=%g, expected: %g\n",
- foreign_lambda, s->foreign_lambda);
- fprintf(stderr, "Got derivative=%d, expected: %d\n",
- derivative, s->derivative);
+ char buf[STRLEN], buf2[STRLEN];
+ lambda_vec_print(foreign_lambda, buf, FALSE);
+ lambda_vec_print(s->foreign_lambda, buf2, FALSE);
+ fprintf(stderr, "Got foreign lambda=%s, expected: %s\n", buf, buf2);
gmx_fatal(FARGS, "Corrupted data in file %s around t=%g.",
filename, start_time);
}
static samples_t *read_edr_hist_block(int *nsamples, t_enxblock *blk,
double start_time, double delta_time,
- double native_lambda, double temp,
+ lambda_vec_t *native_lambda, double temp,
double *last_t, const char *filename)
{
int i,j;
samples_t *s;
int nhist;
- double foreign_lambda;
- int derivative;
+ double old_foreign_lambda;
+ lambda_vec_t *foreign_lambda;
+ int type;
int nbins[2];
/* check the block types etc. */
snew(s, 1);
*nsamples=1;
- foreign_lambda=blk->sub[0].dval[0];
- derivative=(int)(blk->sub[1].lval[1]);
- if (derivative)
- foreign_lambda=native_lambda;
+ snew(foreign_lambda, 1);
+ lambda_vec_init(foreign_lambda, native_lambda->lc);
+ lambda_vec_copy(foreign_lambda, native_lambda);
+ type = (int)(blk->sub[1].lval[1]);
+ if (type == dhbtDH)
+ {
+ double old_foreign_lambda;
+
+ old_foreign_lambda=blk->sub[0].dval[0];
+ if (old_foreign_lambda >= 0)
+ {
+ foreign_lambda->val[0]=old_foreign_lambda;
+ if (foreign_lambda->lc->N > 1)
+ {
+ gmx_fatal(FARGS,
+ "Single-component lambda in multi-component file %s",
+ filename);
+ }
+ }
+ else
+ {
+ for(i=0;i<native_lambda->lc->N;i++)
+ {
+ foreign_lambda->val[i] = blk->sub[0].dval[i+2];
+ }
+ }
+ }
+ else
+ {
+ if (foreign_lambda->lc->N > 1)
+ {
+ if (blk->sub[1].nr < 3 + nhist)
+ {
+ gmx_fatal(FARGS,
+ "Missing derivative coord in multi-component file %s",
+ filename);
+ }
+ foreign_lambda->dhdl = blk->sub[1].lval[2 + nhist];
+ }
+ else
+ {
+ foreign_lambda->dhdl = 0;
+ }
+ }
- samples_init(s, native_lambda, foreign_lambda, temp,
- derivative!=0, filename);
+ samples_init(s, native_lambda, foreign_lambda, temp, type==dhbtDHDL,
+ filename);
snew(s->hist, 1);
for(i=0;i<nhist;i++)
for(i=0;i<nhist;i++)
{
- s->hist->x0[i]=blk->sub[1].lval[2+i];
+ s->hist->x0[i] = blk->sub[1].lval[2+i];
s->hist->dx[i] = blk->sub[0].dval[1];
if (i==1)
+ {
s->hist->dx[i] = - s->hist->dx[i];
+ }
}
s->hist->start_time = start_time;
}
-static void read_barsim_edr(char *fn, real *temp, lambda_t *lambda_head)
+static void read_barsim_edr(char *fn, real *temp, sim_data_t *sd)
{
- int i;
+ int i,j;
ener_file_t fp;
t_enxframe *fr;
int nre;
samples_t **samples_rawdh=NULL; /* contains samples for raw delta_h */
int *nhists=NULL; /* array to keep count & print at end */
int *npts=NULL; /* array to keep count & print at end */
- double *lambdas=NULL; /* array to keep count & print at end */
- double native_lambda=-1;
+ lambda_vec_t **lambdas=NULL; /* array to keep count & print at end */
+ lambda_vec_t *native_lambda;
double end_time; /* the end time of the last batch of samples */
int nsamples=0;
+ lambda_vec_t start_lambda;
fp = open_enx(fn,"r");
do_enxnms(fp,&nre,&enm);
snew(fr, 1);
+ snew(native_lambda, 1);
+ start_lambda.lc = NULL;
+
while(do_enx(fp, fr))
{
/* count the data blocks */
int nlam=0;
int k;
/* DHCOLL block information: */
- double start_time=0, delta_time=0, start_lambda=0, delta_lambda=0;
+ double start_time=0, delta_time=0, old_start_lambda=0, delta_lambda=0;
double rtemp=0;
/* count the blocks and handle collection information: */
rtemp = 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];
+ old_start_lambda = fr->block[i].sub[0].dval[3];
delta_lambda = fr->block[i].sub[0].dval[4];
if (delta_lambda>0)
}
*temp=rtemp;
+ if (old_start_lambda >= 0)
+ {
+ if (sd->lc.N > 0)
+ {
+ if (!lambda_components_check(&(sd->lc), 0, "", 0))
+ {
+ gmx_fatal(FARGS,
+ "lambda vector components in %s don't match those previously read",
+ fn);
+ }
+ }
+ else
+ {
+ lambda_components_add(&(sd->lc), "", 0);
+ }
+ if (!start_lambda.lc)
+ {
+ lambda_vec_init(&start_lambda, &(sd->lc));
+ }
+ start_lambda.val[0]=old_start_lambda;
+ }
+ else
+ {
+ /* read lambda vector */
+ int n_lambda_vec;
+ gmx_bool check=(sd->lc.N > 0);
+ if (fr->block[i].nsub < 2)
+ {
+ gmx_fatal(FARGS,
+ "No lambda vector, but start_lambda=%g\n",
+ old_start_lambda);
+ }
+ n_lambda_vec=fr->block[i].sub[1].ival[1];
+ for(j=0;j<n_lambda_vec;j++)
+ {
+ const char *name=
+ efpt_singular_names[fr->block[i].sub[1].ival[1+j]];
+ if (check)
+ {
+ /* check the components */
+ lambda_components_check(&(sd->lc), j, name,
+ strlen(name));
+ }
+ else
+ {
+ lambda_components_add(&(sd->lc), name,
+ strlen(name));
+ }
+ }
+ lambda_vec_init(&start_lambda, &(sd->lc));
+ start_lambda.index=fr->block[i].sub[1].ival[0];
+ for(j=0;j<n_lambda_vec;j++)
+ {
+ start_lambda.val[j]=fr->block[i].sub[0].dval[5+j];
+ }
+ }
if (first_t < 0)
first_t=start_time;
}
if (nlam != 1)
{
- gmx_fatal(FARGS, "Did not find a delta h information in file %s" , fn);
+ gmx_fatal(FARGS, "Did not find delta H information in file %s", fn);
}
if (nblocks_raw > 0 && nblocks_hist > 0 )
{
if (nsamples > 0)
{
/* check the native lambda */
- if (!lambda_same(start_lambda, native_lambda) )
+ if (!lambda_vec_same(&start_lambda, native_lambda) )
{
gmx_fatal(FARGS, "Native lambda not constant in file %s: started at %g, and becomes %g at time %g",
fn, native_lambda, start_lambda, start_time);
if (samples_rawdh[i])
{
/* insert it into the existing list */
- lambda_list_insert_sample(lambda_head,
- samples_rawdh[i]);
+ lambda_data_list_insert_sample(sd->lb,
+ samples_rawdh[i]);
/* and make sure we'll allocate a new one this time
around */
samples_rawdh[i]=NULL;
{
/* this is the first round; allocate the associated data
structures */
- native_lambda=start_lambda;
+ /*native_lambda=start_lambda;*/
+ lambda_vec_init(native_lambda, &(sd->lc));
+ lambda_vec_copy(native_lambda, &start_lambda);
nsamples=nblocks_raw+nblocks_hist;
snew(nhists, nsamples);
snew(npts, nsamples);
{
nhists[i]=0;
npts[i]=0;
- lambdas[i]=-1;
+ lambdas[i]=NULL;
samples_rawdh[i]=NULL; /* init to NULL so we know which
ones contain values */
}
k=0; /* counter for the lambdas, etc. arrays */
for(i=0;i<fr->nblock;i++)
{
- if (fr->block[i].id == enxDH)
- {
- int ndu;
- read_edr_rawdh_block(&(samples_rawdh[k]),
- &ndu,
- &(fr->block[i]),
- start_time, delta_time,
- start_lambda, rtemp,
- &last_t, fn);
- npts[k] += ndu;
- if (samples_rawdh[k])
+ if (fr->block[i].id == enxDH)
+ {
+ int type=(fr->block[i].sub[0].ival[0]);
+ if (type == dhbtDH || type == dhbtDHDL)
{
- lambdas[k]=samples_rawdh[k]->foreign_lambda;
+ int ndu;
+ read_edr_rawdh_block(&(samples_rawdh[k]),
+ &ndu,
+ &(fr->block[i]),
+ start_time, delta_time,
+ native_lambda, rtemp,
+ &last_t, fn);
+ npts[k] += ndu;
+ if (samples_rawdh[k])
+ {
+ lambdas[k]=samples_rawdh[k]->foreign_lambda;
+ }
+ k++;
}
- k++;
}
else if (fr->block[i].id == enxDHHIST)
{
- int j;
- int nb=0;
- samples_t *s; /* this is where the data will go */
- s=read_edr_hist_block(&nb, &(fr->block[i]),
- start_time, delta_time,
- start_lambda, rtemp,
- &last_t, fn);
- nhists[k] += nb;
- if (nb>0)
+ int type=(int)(fr->block[i].sub[1].lval[1]);
+ if (type == dhbtDH || type == dhbtDHDL)
{
- lambdas[k]= s->foreign_lambda;
- }
- k++;
- /* and insert the new sample immediately */
- for(j=0;j<nb;j++)
- {
- lambda_list_insert_sample(lambda_head, s+j);
+ int j;
+ int nb=0;
+ samples_t *s; /* this is where the data will go */
+ s=read_edr_hist_block(&nb, &(fr->block[i]),
+ start_time, delta_time,
+ native_lambda, rtemp,
+ &last_t, fn);
+ nhists[k] += nb;
+ if (nb>0)
+ {
+ lambdas[k]=s->foreign_lambda;
+ }
+ k++;
+ /* and insert the new sample immediately */
+ for(j=0;j<nb;j++)
+ {
+ lambda_data_list_insert_sample(sd->lb, s+j);
+ }
}
}
}
if (samples_rawdh[i])
{
/* insert it into the existing list */
- lambda_list_insert_sample(lambda_head, samples_rawdh[i]);
+ lambda_data_list_insert_sample(sd->lb, samples_rawdh[i]);
}
}
- fprintf(stderr, "\n");
- printf("%s: %.1f - %.1f; lambda = %.3f\n foreign lambdas:",
- fn, first_t, last_t, native_lambda);
- for(i=0;i<nsamples;i++)
{
- if (nhists[i] > 0)
- {
- printf(" %.3f (%d hists)", lambdas[i], nhists[i]);
- }
- else
+ char buf[STRLEN];
+ printf("\n");
+ lambda_vec_print(native_lambda, buf, FALSE);
+ printf("%s: %.1f - %.1f; lambda = %s\n foreign lambdas:\n",
+ fn, first_t, last_t, buf);
+ for(i=0;i<nsamples;i++)
{
- printf(" %.3f (%d pts)", lambdas[i], npts[i]);
+ if (lambdas[i])
+ {
+ lambda_vec_print(lambdas[i], buf, TRUE);
+ if (nhists[i] > 0)
+ {
+ printf(" %s (%d hists)\n", buf, nhists[i]);
+ }
+ else
+ {
+ printf(" %s (%d pts)\n", buf, npts[i]);
+ }
+ }
}
}
printf("\n\n");
int nedrfile=0;
char **fxvgnms;
char **fedrnms;
- lambda_t *lb; /* the pre-processed lambda data (linked list head) */
- lambda_t lambda_head; /* the head element */
+ sim_data_t sim_data; /* the simulation data */
barres_t *results; /* the results */
int nresults; /* number of results in results array */
double *partsum;
double prec,dg_tot,dg,sig, dg_tot_max, dg_tot_min;
FILE *fpb,*fpi;
- char lamformat[20];
- char dgformat[20],xvg2format[STRLEN],xvg3format[STRLEN],buf[STRLEN];
+ char dgformat[20],xvg2format[STRLEN],xvg3format[STRLEN];
+ char buf[STRLEN], buf2[STRLEN];
char ktformat[STRLEN], sktformat[STRLEN];
char kteformat[STRLEN], skteformat[STRLEN];
output_env_t oenv;
nedrfile = opt2fns(&fedrnms,"-g",NFILE,fnm);
}
+ sim_data_init(&sim_data);
+#if 0
/* make linked list */
lb=&lambda_head;
- lambda_init(lb, 0, 0);
+ lambda_data_init(lb, 0, 0);
lb->next=lb;
lb->prev=lb;
+#endif
nfile_tot = nxvgfile + nedrfile;
/* read in all files. First xvg files */
for(f=0; f<nxvgfile; f++)
{
- read_bar_xvg(fxvgnms[f],&temp,lb);
+ read_bar_xvg(fxvgnms[f],&temp,&sim_data);
nf++;
}
/* then .edr files */
for(f=0; f<nedrfile; f++)
{
- read_barsim_edr(fedrnms[f],&temp,lb);;
+ read_barsim_edr(fedrnms[f],&temp,&sim_data);;
nf++;
}
/* fix the times to allow for equilibration */
- lambdas_impose_times(lb, begin, end);
+ sim_data_impose_times(&sim_data, begin, end);
if (opt2bSet("-oh",NFILE,fnm))
{
- lambdas_histogram(lb, opt2fn("-oh",NFILE,fnm), nbin, oenv);
+ sim_data_histogram(&sim_data, opt2fn("-oh",NFILE,fnm), nbin, oenv);
}
/* assemble the output structures from the lambdas */
- results=barres_list_create(lb, &nresults, use_dhdl);
+ results=barres_list_create(&sim_data, &nresults, use_dhdl);
sum_disc_err=barres_list_max_disc_err(results, nresults);
}
- sprintf(lamformat,"%%6.3f");
+ /*sprintf(lamformat,"%%6.3f");*/
sprintf( dgformat,"%%%d.%df",3+nd,nd);
/* the format strings of the results in kT */
sprintf( ktformat,"%%%d.%df",5+nd,nd);
/* the format strings of the errors in kT */
sprintf( kteformat,"%%%d.%df",3+nd,nd);
sprintf( skteformat,"%%%ds",4+nd);
- sprintf(xvg2format,"%s %s\n","%g",dgformat);
- sprintf(xvg3format,"%s %s %s\n","%g",dgformat,dgformat);
+ sprintf(xvg2format,"%s %s\n","%s",dgformat);
+ sprintf(xvg3format,"%s %s %s\n","%s",dgformat,dgformat);
printf("\n");
for(f=0; f<nresults; f++)
{
- printf(lamformat, results[f].a->native_lambda);
- printf(" ");
- printf(lamformat, results[f].b->native_lambda);
- printf(" ");
+ lambda_vec_print_short(results[f].a->native_lambda, buf);
+ printf("%s ", buf);
+ lambda_vec_print_short(results[f].b->native_lambda, buf);
+ printf("%s ", buf);
printf(ktformat, results[f].dg);
printf(" ");
if (bEE)
if (fpi != NULL)
{
- fprintf(fpi, xvg2format, results[f].a->native_lambda, dg_tot);
+ lambda_vec_print_short(results[f].a->native_lambda, buf);
+ fprintf(fpi, xvg2format, buf, dg_tot);
}
if (fpb != NULL)
{
- fprintf(fpb, xvg3format,
- 0.5*(results[f].a->native_lambda +
- results[f].b->native_lambda),
- results[f].dg,results[f].dg_err);
+ lambda_vec_print_intermediate(results[f].a->native_lambda,
+ results[f].b->native_lambda,
+ buf);
+
+ fprintf(fpb, xvg3format, buf, results[f].dg,results[f].dg_err);
}
- /*printf("lambda %4.2f - %4.2f, DG ", results[f].lambda_a,
- results[f].lambda_b);*/
- printf("lambda ");
- printf(lamformat, results[f].a->native_lambda);
- printf(" - ");
- printf(lamformat, results[f].b->native_lambda);
+ printf("point ");
+ lambda_vec_print_short(results[f].a->native_lambda, buf);
+ lambda_vec_print_short(results[f].b->native_lambda, buf2);
+ printf("%s - %s", buf, buf2);
printf(", DG ");
printf(dgformat,results[f].dg*kT);
dg_tot += results[f].dg;
}
printf("\n");
- printf("total ");
- printf(lamformat, results[0].a->native_lambda);
- printf(" - ");
- printf(lamformat, results[nresults-1].b->native_lambda);
+ printf("total ");
+ lambda_vec_print_short(results[0].a->native_lambda, buf);
+ lambda_vec_print_short(results[nresults-1].b->native_lambda, buf2);
+ printf("%s - %s", buf, buf2);
printf(", DG ");
printf(dgformat,dg_tot*kT);
if (fpi != NULL)
{
- fprintf(fpi, xvg2format,
- results[nresults-1].b->native_lambda, dg_tot);
+ lambda_vec_print_short(results[nresults-1].b->native_lambda, buf);
+ fprintf(fpi, xvg2format, buf, dg_tot);
ffclose(fpi);
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff