/* 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->lambda_stop_n - ir->fepvals->lambda_start_n;
+ 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=ir->fepvals->lambda_start_n;i<ir->fepvals->lambda_stop_n;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++;
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff