/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2010,2011,2012,2013,2014,2015,2017,2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2010-2018, The GROMACS development team.
+ * Copyright (c) 2019, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
/* 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 */
+ 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 */
+ 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
{
- 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 */
- 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. */
- lambda_vec_t native_lambda; /* the native lambda */
-
- struct xvg_t *next, *prev; /*location in the global linked list of xvg_ts*/
+ 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 */
+ 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. */
+ 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 */
- double dx[2]; /* the histogram spacing. The reverse
- dx is the negative of the forward dx.*/
- int64_t x0[2]; /* the (forward + reverse) histogram start
- point(s) as int */
+ unsigned int* bin[2]; /* the (forward + reverse) histogram values */
+ double dx[2]; /* the histogram spacing. The reverse
+ dx is the negative of the forward dx.*/
+ int64_t x0[2]; /* the (forward + reverse) histogram start
+ point(s) as int */
- int nbin[2]; /* the (forward+reverse) number of bins */
- int64_t sum; /* the total number of counts. Must be
- the same for forward + reverse. */
- int nhist; /* number of hist datas (forward or reverse) */
+ int nbin[2]; /* the (forward+reverse) number of bins */
+ int64_t sum; /* the total number of counts. Must be
+ the same for forward + reverse. */
+ int nhist; /* number of hist datas (forward or reverse) */
- double start_time, delta_time; /* start time, end time of histogram */
+ double start_time, delta_time; /* start time, end time of histogram */
} hist_t;
/* an aggregate of samples for partial free energy calculation */
typedef struct samples_t
{
- lambda_vec_t *native_lambda; /* pointer to native lambda vector */
- lambda_vec_t *foreign_lambda; /* pointer to foreign lambda vector */
+ 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 */
/* The samples come either as either delta U lists: */
int ndu; /* the number of delta U samples */
- double *du; /* the delta u's */
- double *t; /* the times associated with those samples, or: */
+ double* du; /* the delta u's */
+ double* t; /* the times associated with those samples, or: */
double start_time, delta_time; /*start time and delta time for linear time*/
/* or as histograms: */
- hist_t *hist; /* a histogram */
+ hist_t* hist; /* a histogram */
/* allocation data: (not NULL for data 'owned' by this struct) */
- double *du_alloc, *t_alloc; /* allocated delta u arrays */
- size_t ndu_alloc, nt_alloc; /* pre-allocated sizes */
- hist_t *hist_alloc; /* allocated hist */
+ double *du_alloc, *t_alloc; /* allocated delta u arrays */
+ size_t ndu_alloc, nt_alloc; /* pre-allocated sizes */
+ hist_t* hist_alloc; /* allocated hist */
- int64_t ntot; /* total number of samples */
- const char *filename; /* the file name this sample comes from */
+ int64_t ntot; /* total number of samples */
+ const char* filename; /* the file name this sample comes from */
} samples_t;
/* a sample range (start to end for du-style data, or boolean
for both du-style data and histograms */
typedef struct sample_range_t
{
- int start, end; /* start and end index for du style data */
- gmx_bool use; /* whether to use this sample */
+ int start, end; /* start and end index for du style data */
+ gmx_bool use; /* whether to use this sample */
- samples_t *s; /* the samples this range belongs to */
+ samples_t* s; /* the samples this range belongs to */
} sample_range_t;
foreign lambda) */
typedef struct sample_coll_t
{
- 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 */
+ 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 */
- samples_t **s; /* the samples themselves */
- sample_range_t *r; /* the sample ranges */
- int nsamples_alloc; /* number of allocated samples */
+ int nsamples; /* the number of samples */
+ samples_t** s; /* the samples themselves */
+ sample_range_t* r; /* the sample ranges */
+ int nsamples_alloc; /* number of allocated samples */
- int64_t ntot; /* total number of samples in the ranges of
- this collection */
+ int64_t ntot; /* total number of samples in the ranges of
+ this collection */
struct sample_coll_t *next, *prev; /* next and previous in the list */
} sample_coll_t;
/* all the samples associated with a lambda point */
typedef struct lambda_data_t
{
- lambda_vec_t *lambda; /* the native lambda (at start time if dynamic) */
- double temp; /* temperature */
+ 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; /* the samples */
- sample_coll_t sc_head; /*the pre-allocated list head for the linked list.*/
+ sample_coll_t sc_head; /*the pre-allocated list head for the linked list.*/
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_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 */
+ lambda_components_t lc; /* the allowed components of the lambda
+ vectors */
} sim_data_t;
/* Top-level data structure with calculated values. */
-typedef struct {
- sample_coll_t *a, *b; /* the simulation data */
+typedef struct
+{
+ sample_coll_t *a, *b; /* the simulation data */
- double dg; /* the free energy difference */
- double dg_err; /* the free energy difference */
+ double dg; /* the free energy difference */
+ double dg_err; /* the free energy difference */
- double dg_disc_err; /* discretization error */
- double dg_histrange_err; /* histogram range error */
+ double dg_disc_err; /* discretization error */
+ double dg_histrange_err; /* histogram range error */
- double sa; /* relative entropy of b in state a */
- double sa_err; /* error in sa */
- double sb; /* relative entropy of a in state b */
- double sb_err; /* error in sb */
+ double sa; /* relative entropy of b in state a */
+ double sa_err; /* error in sa */
+ double sb; /* relative entropy of a in state b */
+ double sb_err; /* error in sb */
- double dg_stddev; /* expected dg stddev per sample */
- double dg_stddev_err; /* error in dg_stddev */
+ double dg_stddev; /* expected dg stddev per sample */
+ double dg_stddev_err; /* error in dg_stddev */
} barres_t;
/* Initialize a lambda_components structure */
-static void lambda_components_init(lambda_components_t *lc)
+static void lambda_components_init(lambda_components_t* lc)
{
lc->N = 0;
lc->Nalloc = 2;
}
/* Add a component to a lambda_components structure */
-static void lambda_components_add(lambda_components_t *lc,
- const char *name, size_t name_length)
+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;
+ lc->Nalloc = (lc->Nalloc == 0) ? 2 : 2 * lc->Nalloc;
srenew(lc->names, lc->Nalloc);
}
- snew(lc->names[lc->N], name_length+1);
+ snew(lc->names[lc->N], name_length + 1);
std::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)
+static gmx_bool lambda_components_check(const lambda_components_t* lc, int index, const char* name, size_t name_length)
{
size_t len;
if (!lc || index >= lc->N)
{
return TRUE;
}
- if (((name != nullptr) && (lc->names[index] == nullptr)) ||
- ((name == nullptr) && (lc->names[index] != nullptr)))
+ if (((name != nullptr) && (lc->names[index] == nullptr))
+ || ((name == nullptr) && (lc->names[index] != nullptr)))
{
return FALSE;
}
- GMX_RELEASE_ASSERT((name != nullptr) || (name_length == 0),
- "If name is empty, the length of the substring to examine within it must be zero");
+ GMX_RELEASE_ASSERT(
+ (name != nullptr) || (name_length == 0),
+ "If name is empty, the length of the substring to examine within it must be zero");
len = std::strlen(lc->names[index]);
if (len != name_length)
{
}
/* 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)
+static int lambda_components_find(const lambda_components_t* lc, const char* name, size_t name_length)
{
int i;
}
-
/* initialize a lambda vector */
-static void lambda_vec_init(lambda_vec_t *lv, const lambda_components_t *lc)
+static void lambda_vec_init(lambda_vec_t* lv, const lambda_components_t* lc)
{
snew(lv->val, lc->N);
lv->index = -1;
lv->lc = lc;
}
-static void lambda_vec_copy(lambda_vec_t *lv, const lambda_vec_t *orig)
+static void lambda_vec_copy(lambda_vec_t* lv, const lambda_vec_t* orig)
{
int i;
}
/* write a lambda vec to a preallocated string */
-static void lambda_vec_print(const lambda_vec_t *lv, char *str, gmx_bool named)
+static void lambda_vec_print(const lambda_vec_t* lv, char* str, gmx_bool named)
{
- int i;
+ int i;
str[0] = 0; /* reset the string */
if (lv->dhdl < 0)
for (i = 0; i < lv->lc->N; i++)
{
str += sprintf(str, "%g", lv->val[i]);
- if (i < lv->lc->N-1)
+ if (i < lv->lc->N - 1)
{
str += sprintf(str, ", ");
}
}
/* 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)
+static void lambda_vec_print_short(const lambda_vec_t* lv, char* str)
{
if (lv->index >= 0)
{
}
/* 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)
+static void lambda_vec_print_intermediate(const lambda_vec_t* a, const lambda_vec_t* b, char* str)
{
str[0] = 0;
- if ( (a->index >= 0) && (b->index >= 0) )
+ if ((a->index >= 0) && (b->index >= 0))
{
- sprintf(str, "%6.3f", (a->index+b->index)/2.0);
+ sprintf(str, "%6.3f", (a->index + b->index) / 2.0);
}
else
{
- if ( (a->dhdl < 0) && (b->dhdl < 0) )
+ if ((a->dhdl < 0) && (b->dhdl < 0))
{
- sprintf(str, "%6.3f", (a->val[0]+b->val[0])/2.0);
+ sprintf(str, "%6.3f", (a->val[0] + b->val[0]) / 2.0);
}
}
}
/* 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)
+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) )
+ if ((a->dhdl > 0) || (b->dhdl > 0))
{
- gmx_fatal(FARGS,
- "Trying to calculate the difference between derivatives instead of lambda points");
+ 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");
+ 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;
+ ret += df * df;
}
return std::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)
+static gmx_bool lambda_vec_same(const lambda_vec_t* a, const lambda_vec_t* b)
{
int i;
{
for (i = 0; i < a->lc->N; i++)
{
- if (!gmx_within_tol(a->val[i], b->val[i], 10*GMX_REAL_EPS))
+ if (!gmx_within_tol(a->val[i], b->val[i], 10 * GMX_REAL_EPS))
{
return FALSE;
}
returns 1 if a is 'bigger' than b,
returns 0 if they're the same,
returns -1 if a is 'smaller' than b.*/
-static int lambda_vec_cmp_foreign(const lambda_vec_t *a,
- const lambda_vec_t *b)
+static int lambda_vec_cmp_foreign(const lambda_vec_t* a, const lambda_vec_t* b)
{
int i;
- double norm_a = 0, norm_b = 0;
+ double norm_a = 0, norm_b = 0;
gmx_bool different = FALSE;
if (a->lc != b->lc)
{
/* lambda vectors that are derivatives always sort higher than those
without derivatives */
- if ((a->dhdl >= 0) != (b->dhdl >= 0) )
+ if ((a->dhdl >= 0) != (b->dhdl >= 0))
{
return (a->dhdl >= 0) ? 1 : -1;
}
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))
+ 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];
+ norm_a += a->val[i] * a->val[i];
+ norm_b += b->val[i] * b->val[i];
}
if (!different)
{
returns 1 if a is 'bigger' than b,
returns 0 if they're the same,
returns -1 if a is 'smaller' than b.*/
-static int lambda_vec_cmp_native(const lambda_vec_t *a,
- const lambda_vec_t *b)
+static int lambda_vec_cmp_native(const lambda_vec_t* a, const lambda_vec_t* b)
{
if (a->lc != b->lc)
{
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");
+ 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");
+ 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))
+ if (gmx_within_tol(a->val[0], b->val[0], 10 * GMX_REAL_EPS))
{
return 0;
}
}
-
-
-static void hist_init(hist_t *h, int nhist, int *nbin)
+static void hist_init(hist_t* h, int nhist, int* nbin)
{
int i;
if (nhist > 2)
h->x0[i] = 0;
h->nbin[i] = nbin[i];
h->start_time = h->delta_time = 0;
- h->dx[i] = 0;
+ h->dx[i] = 0;
}
h->sum = 0;
h->nhist = nhist;
}
-static void xvg_init(xvg_t *ba)
+static void xvg_init(xvg_t* ba)
{
ba->filename = nullptr;
ba->nset = 0;
ba->y = nullptr;
}
-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)
+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;
s->foreign_lambda = foreign_lambda;
s->du = nullptr;
s->t = nullptr;
s->start_time = s->delta_time = 0;
- s->hist = nullptr;
- s->du_alloc = nullptr;
- s->t_alloc = nullptr;
- s->hist_alloc = nullptr;
- s->ndu_alloc = 0;
- s->nt_alloc = 0;
+ s->hist = nullptr;
+ s->du_alloc = nullptr;
+ s->t_alloc = nullptr;
+ s->hist_alloc = nullptr;
+ s->ndu_alloc = 0;
+ s->nt_alloc = 0;
s->ntot = 0;
s->filename = filename;
}
-static void sample_range_init(sample_range_t *r, samples_t *s)
+static void sample_range_init(sample_range_t* r, samples_t* s)
{
r->start = 0;
r->end = s->ndu;
r->s = nullptr;
}
-static void sample_coll_init(sample_coll_t *sc, lambda_vec_t *native_lambda,
- lambda_vec_t *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;
sc->next = sc->prev = nullptr;
}
-static void sample_coll_destroy(sample_coll_t *sc)
+static void sample_coll_destroy(sample_coll_t* sc)
{
/* don't free the samples themselves */
sfree(sc->r);
}
-static void lambda_data_init(lambda_data_t *l, lambda_vec_t *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->prev = l->sc;
}
-static void barres_init(barres_t *br)
+static void barres_init(barres_t* br)
{
br->dg = 0;
br->dg_err = 0;
/* calculate the total number of samples in a sample collection */
-static void sample_coll_calc_ntot(sample_coll_t *sc)
+static void sample_coll_calc_ntot(sample_coll_t* sc)
{
int i;
/* find the barsamples_t associated with a lambda that corresponds to
a specific foreign lambda */
-static sample_coll_t *lambda_data_find_sample_coll(lambda_data_t *l,
- lambda_vec_t *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;
+ sample_coll_t* sc = l->sc->next;
while (sc != l->sc)
{
}
/* insert li into an ordered list of lambda_colls */
-static void lambda_data_insert_sample_coll(lambda_data_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) )
+ sample_coll_t* scn = l->sc->next;
+ while ((scn != l->sc))
{
if (lambda_vec_cmp_foreign(scn->foreign_lambda, sc->foreign_lambda) > 0)
{
}
/* insert li into an ordered list of lambdas */
-static void lambda_data_insert_lambda(lambda_data_t *head, lambda_data_t *li)
+static void lambda_data_insert_lambda(lambda_data_t* head, lambda_data_t* li)
{
- lambda_data_t *lc = head->next;
+ lambda_data_t* lc = head->next;
while (lc != head)
{
if (lambda_vec_cmp_native(lc->lambda, li->lambda) > 0)
/* insert a sample and a sample_range into a sample_coll. The
samples are stored as a pointer, the range is copied. */
-static void sample_coll_insert_sample(sample_coll_t *sc, samples_t *s,
- sample_range_t *r)
+static void sample_coll_insert_sample(sample_coll_t* sc, samples_t* s, sample_range_t* r)
{
/* first check if it belongs here */
GMX_ASSERT(sc->next->s, "Next not properly initialized!");
if (sc->temp != s->temp)
{
- gmx_fatal(FARGS, "Temperatures in files %s and %s are not the same!",
- s->filename, sc->next->s[0]->filename);
+ gmx_fatal(FARGS, "Temperatures in files %s and %s are not the same!", s->filename,
+ sc->next->s[0]->filename);
}
if (!lambda_vec_same(sc->native_lambda, s->native_lambda))
{
}
/* check if there's room */
- if ( (sc->nsamples + 1) > sc->nsamples_alloc)
+ if ((sc->nsamples + 1) > sc->nsamples_alloc)
{
- sc->nsamples_alloc = std::max(2*sc->nsamples_alloc, 2);
+ sc->nsamples_alloc = std::max(2 * sc->nsamples_alloc, 2);
srenew(sc->s, sc->nsamples_alloc);
srenew(sc->r, sc->nsamples_alloc);
}
/* insert a sample into a lambda_list, creating the right sample_coll if
neccesary */
-static void lambda_data_list_insert_sample(lambda_data_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_coll_t* sc;
sample_range_t r;
- lambda_data_t *l = head->next;
+ lambda_data_t* l = head->next;
/* first search for the right lambda_data_t */
while (l != head)
{
- if (lambda_vec_same(l->lambda, s->native_lambda) )
+ if (lambda_vec_same(l->lambda, s->native_lambda))
{
found = TRUE;
break;
/* make a histogram out of a sample collection */
-static void sample_coll_make_hist(sample_coll_t *sc, int **bin,
- int *nbin_alloc, int *nbin,
- double *dx, double *xmin, int nbin_default)
+static void
+sample_coll_make_hist(sample_coll_t* sc, int** bin, int* nbin_alloc, int* nbin, double* dx, double* xmin, int nbin_default)
{
int i, j, k;
gmx_bool dx_set = FALSE;
gmx_bool xmax_set = FALSE;
gmx_bool xmax_set_hard = FALSE; /* whether the xmax is bounded by the
limits of a histogram */
- double xmax = -1;
+ double xmax = -1;
/* first determine dx and xmin; try the histograms */
for (i = 0; i < sc->nsamples; i++)
{
if (sc->s[i]->hist)
{
- hist_t *hist = sc->s[i]->hist;
+ hist_t* hist = sc->s[i]->hist;
for (k = 0; k < hist->nhist; k++)
{
double hdx = hist->dx[k];
- double xmax_now = (hist->x0[k]+hist->nbin[k])*hdx;
+ double xmax_now = (hist->x0[k] + hist->nbin[k]) * hdx;
/* we use the biggest dx*/
- if ( (!dx_set) || hist->dx[0] > *dx)
+ if ((!dx_set) || hist->dx[0] > *dx)
{
dx_set = TRUE;
*dx = hist->dx[0];
}
- if ( (!xmin_set) || (hist->x0[k]*hdx) < *xmin)
+ if ((!xmin_set) || (hist->x0[k] * hdx) < *xmin)
{
xmin_set = TRUE;
- *xmin = (hist->x0[k]*hdx);
+ *xmin = (hist->x0[k] * hdx);
}
- if ( (!xmax_set) || (xmax_now > xmax && !xmax_set_hard) )
+ if ((!xmax_set) || (xmax_now > xmax && !xmax_set_hard))
{
xmax_set = TRUE;
xmax = xmax_now;
- if (hist->bin[k][hist->nbin[k]-1] != 0)
+ if (hist->bin[k][hist->nbin[k] - 1] != 0)
{
xmax_set_hard = TRUE;
}
}
- if (hist->bin[k][hist->nbin[k]-1] != 0 && (xmax_now < xmax) )
+ if (hist->bin[k][hist->nbin[k] - 1] != 0 && (xmax_now < xmax))
{
xmax_set_hard = TRUE;
xmax = xmax_now;
int endi = sc->r[i].end;
double du_xmin = sc->s[i]->du[starti];
double du_xmax = sc->s[i]->du[starti];
- for (j = starti+1; j < endi; j++)
+ for (j = starti + 1; j < endi; j++)
{
if (sc->s[i]->du[j] < du_xmin)
{
}
/* and now change the limits */
- if ( (!xmin_set) || (du_xmin < *xmin) )
+ if ((!xmin_set) || (du_xmin < *xmin))
{
xmin_set = TRUE;
*xmin = du_xmin;
}
- if ( (!xmax_set) || ((du_xmax > xmax) && !xmax_set_hard) )
+ if ((!xmax_set) || ((du_xmax > xmax) && !xmax_set_hard))
{
xmax_set = TRUE;
xmax = du_xmax;
if (!dx_set)
{
*nbin = nbin_default;
- *dx = (xmax-(*xmin))/((*nbin)-2); /* -2 because we want the last bin to
- be 0, and we count from 0 */
+ *dx = (xmax - (*xmin)) / ((*nbin) - 2); /* -2 because we want the last bin to
+ be 0, and we count from 0 */
}
else
{
- *nbin = static_cast<int>((xmax-(*xmin))/(*dx));
+ *nbin = static_cast<int>((xmax - (*xmin)) / (*dx));
}
if (*nbin > *nbin_alloc)
{
if (sc->s[i]->hist)
{
- hist_t *hist = sc->s[i]->hist;
+ hist_t* hist = sc->s[i]->hist;
for (k = 0; k < hist->nhist; k++)
{
double hdx = hist->dx[k];
- double xmin_hist = hist->x0[k]*hdx;
+ double xmin_hist = hist->x0[k] * hdx;
for (j = 0; j < hist->nbin[k]; j++)
{
/* calculate the bin corresponding to the middle of the
original bin */
- double x = hdx*(j+0.5) + xmin_hist;
- int binnr = static_cast<int>((x-(*xmin))/(*dx));
+ double x = hdx * (j + 0.5) + xmin_hist;
+ int binnr = static_cast<int>((x - (*xmin)) / (*dx));
if (binnr >= *nbin || binnr < 0)
{
- binnr = (*nbin)-1;
+ binnr = (*nbin) - 1;
}
(*bin)[binnr] += hist->bin[k][j];
int endi = sc->r[i].end;
for (j = starti; j < endi; j++)
{
- int binnr = static_cast<int>((sc->s[i]->du[j]-(*xmin))/(*dx));
+ int binnr = static_cast<int>((sc->s[i]->du[j] - (*xmin)) / (*dx));
if (binnr >= *nbin || binnr < 0)
{
- binnr = (*nbin)-1;
+ binnr = (*nbin) - 1;
}
(*bin)[binnr]++;
}
/* write a collection of histograms to a file */
-static void sim_data_histogram(sim_data_t *sd, const char *filename,
- int nbin_default, const gmx_output_env_t *oenv)
+static void sim_data_histogram(sim_data_t* sd, const char* filename, int nbin_default, const gmx_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_data_t *bl;
+ const char * dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
+ const char* title = "N(\\DeltaH)";
+ const char* label_y = "Samples";
+ FILE* fp;
+ lambda_data_t* bl;
int nsets = 0;
- char **setnames = nullptr;
+ char** setnames = nullptr;
gmx_bool first_set = FALSE;
/* histogram data: */
- int *hist = nullptr;
+ int* hist = nullptr;
int nbin = 0;
int nbin_alloc = 0;
double dx = 0;
double minval = 0;
int i;
- lambda_data_t *bl_head = sd->lb;
+ lambda_data_t* bl_head = sd->lb;
printf("\nWriting histogram to %s\n", filename);
sprintf(label_x, "\\DeltaH (%s)", unit_energy);
/* iterate over all lambdas */
while (bl != bl_head)
{
- sample_coll_t *sc = bl->sc->next;
+ sample_coll_t* sc = bl->sc->next;
/* iterate over all samples */
while (sc != bl->sc)
nsets++;
srenew(setnames, nsets);
- snew(setnames[nsets-1], STRLEN);
+ snew(setnames[nsets - 1], STRLEN);
if (sc->foreign_lambda->dhdl < 0)
{
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);
+ sprintf(setnames[nsets - 1], "N(%s(%s=%s) | %s=%s)", deltag, lambda, buf2, lambda, buf);
}
else
{
lambda_vec_print(sc->native_lambda, buf, FALSE);
- sprintf(setnames[nsets-1], "N(%s | %s=%s)",
- dhdl, lambda, buf);
+ sprintf(setnames[nsets - 1], "N(%s | %s=%s)", dhdl, lambda, buf);
}
sc = sc->next;
}
/* iterate over all lambdas */
while (bl != bl_head)
{
- sample_coll_t *sc = bl->sc->next;
+ sample_coll_t* sc = bl->sc->next;
/* iterate over all samples */
while (sc != bl->sc)
xvgr_new_dataset(fp, 0, 0, nullptr, oenv);
}
- sample_coll_make_hist(sc, &hist, &nbin_alloc, &nbin, &dx, &minval,
- nbin_default);
+ sample_coll_make_hist(sc, &hist, &nbin_alloc, &nbin, &dx, &minval, nbin_default);
for (i = 0; i < nbin; i++)
{
- double xmin = i*dx + minval;
- double xmax = (i+1)*dx + minval;
+ double xmin = i * dx + minval;
+ double xmax = (i + 1) * dx + minval;
fprintf(fp, "%g %d\n%g %d\n", xmin, hist[i], xmax, hist[i]);
}
xvgrclose(fp);
}
-static int
-snprint_lambda_vec(char *str, int sz, const char *label, lambda_vec_t *lambda)
+static int snprint_lambda_vec(char* str, int sz, const char* label, lambda_vec_t* lambda)
{
int n = 0;
- n += snprintf(str+n, sz-n, "lambda vector [%s]: ", label);
+ n += snprintf(str + n, sz - n, "lambda vector [%s]: ", label);
if (lambda->index >= 0)
{
- n += snprintf(str+n, sz-n, " init-lambda-state=%d", lambda->index);
+ n += snprintf(str + n, sz - n, " init-lambda-state=%d", lambda->index);
}
if (lambda->dhdl >= 0)
{
- n += snprintf(str+n, sz-n, " dhdl index=%d", lambda->dhdl);
+ n += snprintf(str + n, sz - n, " dhdl index=%d", lambda->dhdl);
}
else
{
int i;
for (i = 0; i < lambda->lc->N; i++)
{
- n += snprintf(str+n, sz-n, " (%s) l=%g", lambda->lc->names[i], lambda->val[i]);
+ n += snprintf(str + n, sz - n, " (%s) l=%g", lambda->lc->names[i], lambda->val[i]);
}
}
return n;
/* create a collection (array) of barres_t object given a ordered linked list
of barlamda_t sample collections */
-static barres_t *barres_list_create(sim_data_t *sd, int *nres,
- gmx_bool use_dhdl)
+static barres_t* barres_list_create(sim_data_t* sd, int* nres, gmx_bool use_dhdl)
{
- lambda_data_t *bl;
+ lambda_data_t* bl;
int nlambda = 0;
- barres_t *res;
+ barres_t* res;
gmx_bool dhdl = FALSE;
gmx_bool first = TRUE;
- lambda_data_t *bl_head = sd->lb;
+ lambda_data_t* bl_head = sd->lb;
/* first count the lambdas */
bl = bl_head->next;
nlambda++;
bl = bl->next;
}
- snew(res, nlambda-1);
+ snew(res, nlambda - 1);
/* next put the right samples in the res */
*nres = 0;
while (bl != bl_head)
{
sample_coll_t *sc, *scprev;
- barres_t *br = &(res[*nres]);
+ barres_t* br = &(res[*nres]);
/* there is always a previous one. we search for that as a foreign
lambda: */
scprev = lambda_data_find_sample_coll(bl->prev, bl->lambda);
if (first)
{
- printf("\nWARNING: Using the derivative data (dH/dlambda) to extrapolate delta H values.\nThis will only work if the Hamiltonian is linear in lambda.\n");
+ printf("\nWARNING: Using the derivative data (dH/dlambda) to extrapolate delta H "
+ "values.\nThis will only work if the Hamiltonian is linear in lambda.\n");
dhdl = TRUE;
}
if (!dhdl)
{
- gmx_fatal(FARGS, "Some dhdl files contain only one value (dH/dl), while others \ncontain multiple values (dH/dl and/or Delta H), will not proceed \nbecause of possible inconsistencies.\n");
+ gmx_fatal(FARGS,
+ "Some dhdl files contain only one value (dH/dl), while others \ncontain "
+ "multiple values (dH/dl and/or Delta H), will not proceed \nbecause of "
+ "possible inconsistencies.\n");
}
}
else if (!scprev && !sc)
snprint_lambda_vec(descX, STRLEN, "X", bl->prev->lambda);
snprint_lambda_vec(descY, STRLEN, "Y", bl->lambda);
- gmx_fatal(FARGS, "There is no path between the states X & Y below that is covered by foreign lambdas:\ncannot proceed with BAR.\nUse thermodynamic integration of dH/dl by calculating the averages of dH/dl\nwith g_analyze and integrating them.\nAlternatively, use the -extp option if (and only if) the Hamiltonian\ndepends linearly on lambda, which is NOT normally the case.\n\n%s\n%s\n", descX, descY);
+ gmx_fatal(FARGS,
+ "There is no path between the states X & Y below that is covered by foreign "
+ "lambdas:\ncannot proceed with BAR.\nUse thermodynamic integration of dH/dl "
+ "by calculating the averages of dH/dl\nwith g_analyze and integrating "
+ "them.\nAlternatively, use the -extp option if (and only if) the "
+ "Hamiltonian\ndepends linearly on lambda, which is NOT normally the "
+ "case.\n\n%s\n%s\n",
+ descX, descY);
}
/* normal delta H */
char descX[STRLEN], descY[STRLEN];
snprint_lambda_vec(descX, STRLEN, "X", bl->lambda);
snprint_lambda_vec(descY, STRLEN, "Y", bl->prev->lambda);
- gmx_fatal(FARGS, "Could not find a set for foreign lambda (state X below)\nin the files for main lambda (state Y below)\n\n%s\n%s\n", descX, descY);
+ gmx_fatal(FARGS,
+ "Could not find a set for foreign lambda (state X below)\nin the files for "
+ "main lambda (state Y below)\n\n%s\n%s\n",
+ descX, descY);
}
if (!sc)
{
char descX[STRLEN], descY[STRLEN];
snprint_lambda_vec(descX, STRLEN, "X", bl->prev->lambda);
snprint_lambda_vec(descY, STRLEN, "Y", bl->lambda);
- gmx_fatal(FARGS, "Could not find a set for foreign lambda (state X below)\nin the files for main lambda (state Y below)\n\n%s\n%s\n", descX, descY);
+ gmx_fatal(FARGS,
+ "Could not find a set for foreign lambda (state X below)\nin the files for "
+ "main lambda (state Y below)\n\n%s\n%s\n",
+ descX, descY);
}
br->a = scprev;
br->b = sc;
}
/* estimate the maximum discretization error */
-static double barres_list_max_disc_err(barres_t *res, int nres)
+static double barres_list_max_disc_err(barres_t* res, int nres)
{
int i, j;
double disc_err = 0.;
for (i = 0; i < nres; i++)
{
- barres_t *br = &(res[i]);
+ barres_t* br = &(res[i]);
- delta_lambda = lambda_vec_abs_diff(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++)
{
double Wfac = 1.;
if (br->a->s[j]->derivative)
{
- Wfac = delta_lambda;
+ Wfac = delta_lambda;
}
- disc_err = std::max(disc_err, Wfac*br->a->s[j]->hist->dx[0]);
+ disc_err = std::max(disc_err, Wfac * br->a->s[j]->hist->dx[0]);
}
}
for (j = 0; j < br->b->nsamples; j++)
double Wfac = 1.;
if (br->b->s[j]->derivative)
{
- Wfac = delta_lambda;
+ Wfac = delta_lambda;
}
- disc_err = std::max(disc_err, Wfac*br->b->s[j]->hist->dx[0]);
+ disc_err = std::max(disc_err, Wfac * br->b->s[j]->hist->dx[0]);
}
}
}
/* impose start and end times on a sample collection, updating sample_ranges */
-static void sample_coll_impose_times(sample_coll_t *sc, double begin_t,
- double end_t)
+static void sample_coll_impose_times(sample_coll_t* sc, double begin_t, double end_t)
{
int i;
for (i = 0; i < sc->nsamples; i++)
{
- samples_t *s = sc->s[i];
- sample_range_t *r = &(sc->r[i]);
+ samples_t* s = sc->s[i];
+ sample_range_t* r = &(sc->r[i]);
if (s->hist)
{
- double end_time = s->hist->delta_time*s->hist->sum +
- s->hist->start_time;
+ double end_time = s->hist->delta_time * s->hist->sum + s->hist->start_time;
if (s->hist->start_time < begin_t || end_time > end_t)
{
r->use = FALSE;
double end_time;
if (s->start_time < begin_t)
{
- r->start = static_cast<int>((begin_t - s->start_time)/s->delta_time);
+ r->start = static_cast<int>((begin_t - s->start_time) / s->delta_time);
}
- end_time = s->delta_time*s->ndu + s->start_time;
+ end_time = s->delta_time * s->ndu + s->start_time;
if (end_time > end_t)
{
- r->end = static_cast<int>((end_t - s->start_time)/s->delta_time);
+ r->end = static_cast<int>((end_t - s->start_time) / s->delta_time);
}
}
else
sample_coll_calc_ntot(sc);
}
-static void sim_data_impose_times(sim_data_t *sd, 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_data_t *lc;
- lambda_data_t *head = sd->lb;
+ lambda_data_t* lc;
+ lambda_data_t* head = sd->lb;
int j;
if (begin <= 0 && end < 0)
lc = head->next;
while (lc != head)
{
- sample_coll_t *sc = lc->sc->next;
+ sample_coll_t* sc = lc->sc->next;
while (sc != lc->sc)
{
for (j = 0; j < sc->nsamples; j++)
double start_t, end_t;
start_t = sc->s[j]->start_time;
- end_t = sc->s[j]->start_time;
+ end_t = sc->s[j]->start_time;
if (sc->s[j]->hist)
{
- end_t += sc->s[j]->delta_time*sc->s[j]->hist->sum;
+ end_t += sc->s[j]->delta_time * sc->s[j]->hist->sum;
}
else
{
if (sc->s[j]->t)
{
- end_t = sc->s[j]->t[sc->s[j]->ndu-1];
+ end_t = sc->s[j]->t[sc->s[j]->ndu - 1];
}
else
{
- end_t += sc->s[j]->delta_time*sc->s[j]->ndu;
+ end_t += sc->s[j]->delta_time * sc->s[j]->ndu;
}
}
lc = head->next;
while (lc != head)
{
- sample_coll_t *sc = lc->sc->next;
+ sample_coll_t* sc = lc->sc->next;
while (sc != lc->sc)
{
sample_coll_impose_times(sc, begin_t, end_t);
/* create subsample i out of ni from an existing sample_coll */
-static gmx_bool sample_coll_create_subsample(sample_coll_t *sc,
- sample_coll_t *sc_orig,
- int i, int ni)
+static gmx_bool sample_coll_create_subsample(sample_coll_t* sc, sample_coll_t* sc_orig, int i, int ni)
{
- int j;
+ int j;
- int64_t ntot_start;
- int64_t ntot_end;
- int64_t ntot_so_far;
+ int64_t ntot_start;
+ int64_t ntot_end;
+ int64_t ntot_so_far;
*sc = *sc_orig; /* just copy all fields */
/* now fix start and end fields */
/* the casts avoid possible overflows */
- ntot_start = static_cast<int64_t>(sc_orig->ntot*static_cast<double>(i)/static_cast<double>(ni));
- ntot_end = static_cast<int64_t>(sc_orig->ntot*static_cast<double>(i+1)/static_cast<double>(ni));
+ ntot_start = static_cast<int64_t>(sc_orig->ntot * static_cast<double>(i) / static_cast<double>(ni));
+ ntot_end = static_cast<int64_t>(sc_orig->ntot * static_cast<double>(i + 1) / static_cast<double>(ni));
ntot_so_far = 0;
for (j = 0; j < sc->nsamples; j++)
{
}
/* check if we're in range at all */
- if ( (new_end < new_start) || (new_start > sc->r[j].end) )
+ if ((new_end < new_start) || (new_start > sc->r[j].end))
{
new_start = 0;
new_end = 0;
}
/* and write the new range */
- GMX_RELEASE_ASSERT(new_start <= std::numeric_limits<int>::max(), "Value of 'new_start' too large for int converstion");
- GMX_RELEASE_ASSERT(new_end <= std::numeric_limits<int>::max(), "Value of 'new_end' too large for int converstion");
+ GMX_RELEASE_ASSERT(new_start <= std::numeric_limits<int>::max(),
+ "Value of 'new_start' too large for int converstion");
+ GMX_RELEASE_ASSERT(new_end <= std::numeric_limits<int>::max(),
+ "Value of 'new_end' too large for int converstion");
sc->r[j].start = static_cast<int>(new_start);
sc->r[j].end = static_cast<int>(new_end);
}
/* first calculate normalized bounds
(where 0 is the start of the hist range, and 1 the end) */
- ntot_start_norm = (ntot_start-ntot_so_far)/static_cast<double>(ntot_add);
- ntot_end_norm = (ntot_end-ntot_so_far)/static_cast<double>(ntot_add);
+ ntot_start_norm = (ntot_start - ntot_so_far) / static_cast<double>(ntot_add);
+ ntot_end_norm = (ntot_end - ntot_so_far) / static_cast<double>(ntot_add);
/* now fix the boundaries */
ntot_start_norm = std::min(1.0, std::max(0.0, ntot_start_norm));
}
/* calculate minimum and maximum work values in sample collection */
-static void sample_coll_min_max(sample_coll_t *sc, double Wfac,
- double *Wmin, double *Wmax)
+static void sample_coll_min_max(sample_coll_t* sc, double Wfac, double* Wmin, double* Wmax)
{
int i, j;
- *Wmin = std::numeric_limits<float>::max();
+ *Wmin = std::numeric_limits<float>::max();
*Wmax = -std::numeric_limits<float>::max();
for (i = 0; i < sc->nsamples; i++)
{
- samples_t *s = sc->s[i];
- sample_range_t *r = &(sc->r[i]);
+ samples_t* s = sc->s[i];
+ sample_range_t* r = &(sc->r[i]);
if (r->use)
{
if (!s->hist)
{
for (j = r->start; j < r->end; j++)
{
- *Wmin = std::min(*Wmin, s->du[j]*Wfac);
- *Wmax = std::max(*Wmax, s->du[j]*Wfac);
+ *Wmin = std::min(*Wmin, s->du[j] * Wfac);
+ *Wmax = std::max(*Wmax, s->du[j] * Wfac);
}
}
else
{
int hd = 0; /* determine the histogram direction: */
double dx;
- if ( (s->hist->nhist > 1) && (Wfac < 0) )
+ if ((s->hist->nhist > 1) && (Wfac < 0))
{
hd = 1;
}
dx = s->hist->dx[hd];
- for (j = s->hist->nbin[hd]-1; j >= 0; j--)
+ for (j = s->hist->nbin[hd] - 1; j >= 0; j--)
{
- *Wmin = std::min(*Wmin, Wfac*(s->hist->x0[hd])*dx);
- *Wmax = std::max(*Wmax, Wfac*(s->hist->x0[hd])*dx);
+ *Wmin = std::min(*Wmin, Wfac * (s->hist->x0[hd]) * dx);
+ *Wmax = std::max(*Wmax, Wfac * (s->hist->x0[hd]) * dx);
/* look for the highest value bin with values */
if (s->hist->bin[hd][j] > 0)
{
- *Wmin = std::min(*Wmin, Wfac*(j+s->hist->x0[hd]+1)*dx);
- *Wmax = std::max(*Wmax, Wfac*(j+s->hist->x0[hd]+1)*dx);
+ *Wmin = std::min(*Wmin, Wfac * (j + s->hist->x0[hd] + 1) * dx);
+ *Wmax = std::max(*Wmax, Wfac * (j + s->hist->x0[hd] + 1) * dx);
break;
}
}
}
/* Initialize a sim_data structure */
-static void sim_data_init(sim_data_t *sd)
+static void sim_data_init(sim_data_t* sd)
{
/* make linked list */
sd->lb = &(sd->lb_head);
}
-static double calc_bar_sum(int n, const double *W, double Wfac, double sbMmDG)
+static double calc_bar_sum(int n, const double* W, double Wfac, double sbMmDG)
{
int i;
double sum;
for (i = 0; i < n; i++)
{
- sum += 1./(1. + std::exp(Wfac*W[i] + sbMmDG));
+ sum += 1. / (1. + std::exp(Wfac * W[i] + sbMmDG));
}
return sum;
if type== 0, calculate the best estimate for the average,
if type==-1, calculate the minimum possible value given the histogram
if type== 1, calculate the maximum possible value given the histogram */
-static double calc_bar_sum_hist(const hist_t *hist, double Wfac, double sbMmDG,
- int type)
+static double calc_bar_sum_hist(const hist_t* hist, double Wfac, double sbMmDG, int type)
{
double sum = 0.;
int i;
int hd = 0; /* determine the histogram direction: */
double dx;
- if ( (hist->nhist > 1) && (Wfac < 0) )
+ if ((hist->nhist > 1) && (Wfac < 0))
{
hd = 1;
}
dx = hist->dx[hd];
- maxbin = hist->nbin[hd]-1;
+ maxbin = hist->nbin[hd] - 1;
if (type == 1)
{
maxbin = hist->nbin[hd]; /* we also add whatever was out of range */
for (i = 0; i < maxbin; i++)
{
- double x = Wfac*((i+hist->x0[hd])+0.5)*dx; /* bin middle */
- double pxdx = hist->bin[0][i]*normdx; /* p(x)dx */
+ double x = Wfac * ((i + hist->x0[hd]) + 0.5) * dx; /* bin middle */
+ double pxdx = hist->bin[0][i] * normdx; /* p(x)dx */
- sum += pxdx/(1. + std::exp(x + sbMmDG));
+ sum += pxdx / (1. + std::exp(x + sbMmDG));
}
return sum;
}
-static double calc_bar_lowlevel(sample_coll_t *ca, sample_coll_t *cb,
- double temp, double tol, int type)
+static double calc_bar_lowlevel(sample_coll_t* ca, sample_coll_t* cb, double temp, double tol, int type)
{
double kT, beta, M;
int i;
double DG0, DG1, DG2, dDG1;
double n1, n2; /* numbers of samples as doubles */
- kT = BOLTZ*temp;
- beta = 1/kT;
+ kT = BOLTZ * temp;
+ beta = 1 / kT;
/* count the numbers of samples */
n1 = ca->ntot;
n2 = cb->ntot;
- M = std::log(n1/n2);
+ M = std::log(n1 / n2);
/*if (!lambda_vec_same(ca->native_lambda, ca->foreign_lambda))*/
if (ca->foreign_lambda->dhdl < 0)
/* 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);
+ 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");
+ gmx_fatal(FARGS, "Can't (yet) do multi-component dhdl interpolation");
}
- Wfac1 = beta*delta_lambda;
- Wfac2 = -beta*delta_lambda;
+ Wfac1 = beta * delta_lambda;
+ Wfac2 = -beta * delta_lambda;
}
if (beta < 1)
smaller than what we get out of the BAR averages.
For the comparison we can use twice the tolerance. */
- while (DG2 - DG0 > 2*tol)
+ while (DG2 - DG0 > 2 * tol)
{
- DG1 = 0.5*(DG0 + DG2);
+ DG1 = 0.5 * (DG0 + DG2);
/* calculate the BAR averages */
dDG1 = 0.;
for (i = 0; i < ca->nsamples; i++)
{
- samples_t *s = ca->s[i];
- sample_range_t *r = &(ca->r[i]);
+ samples_t* s = ca->s[i];
+ sample_range_t* r = &(ca->r[i]);
if (r->use)
{
if (s->hist)
{
- dDG1 += calc_bar_sum_hist(s->hist, Wfac1, (M-DG1), type);
+ dDG1 += calc_bar_sum_hist(s->hist, Wfac1, (M - DG1), type);
}
else
{
- dDG1 += calc_bar_sum(r->end - r->start, s->du + r->start,
- Wfac1, (M-DG1));
+ dDG1 += calc_bar_sum(r->end - r->start, s->du + r->start, Wfac1, (M - DG1));
}
}
}
for (i = 0; i < cb->nsamples; i++)
{
- samples_t *s = cb->s[i];
- sample_range_t *r = &(cb->r[i]);
+ samples_t* s = cb->s[i];
+ sample_range_t* r = &(cb->r[i]);
if (r->use)
{
if (s->hist)
{
- dDG1 -= calc_bar_sum_hist(s->hist, Wfac2, -(M-DG1), type);
+ dDG1 -= calc_bar_sum_hist(s->hist, Wfac2, -(M - DG1), type);
}
else
{
- dDG1 -= calc_bar_sum(r->end - r->start, s->du + r->start,
- Wfac2, -(M-DG1));
+ dDG1 -= calc_bar_sum(r->end - r->start, s->du + r->start, Wfac2, -(M - DG1));
}
}
}
}
}
- return 0.5*(DG0 + DG2);
+ return 0.5 * (DG0 + DG2);
}
-static void calc_rel_entropy(sample_coll_t *ca, sample_coll_t *cb,
- double temp, double dg, double *sa, double *sb)
+static void calc_rel_entropy(sample_coll_t* ca, sample_coll_t* cb, double temp, double dg, double* sa, double* sb)
{
int i, j;
double W_ab = 0.;
double Wfac1, Wfac2;
double n1, n2;
- kT = BOLTZ*temp;
- beta = 1/kT;
+ kT = BOLTZ * temp;
+ beta = 1 / kT;
/* count the numbers of samples */
n1 = ca->ntot;
{
/* we're using dhdl, so delta_lambda needs to be a
multiplication factor. */
- double delta_lambda = lambda_vec_abs_diff(cb->native_lambda,
- ca->native_lambda);
- Wfac1 = beta*delta_lambda;
- Wfac2 = -beta*delta_lambda;
+ double delta_lambda = lambda_vec_abs_diff(cb->native_lambda, ca->native_lambda);
+ Wfac1 = beta * delta_lambda;
+ Wfac2 = -beta * delta_lambda;
}
/* first calculate the average work in both directions */
for (i = 0; i < ca->nsamples; i++)
{
- samples_t *s = ca->s[i];
- sample_range_t *r = &(ca->r[i]);
+ samples_t* s = ca->s[i];
+ sample_range_t* r = &(ca->r[i]);
if (r->use)
{
if (!s->hist)
{
for (j = r->start; j < r->end; j++)
{
- W_ab += Wfac1*s->du[j];
+ W_ab += Wfac1 * s->du[j];
}
}
else
double normdx = 1.;
double dx;
int hd = 0; /* histogram direction */
- if ( (s->hist->nhist > 1) && (Wfac1 < 0) )
+ if ((s->hist->nhist > 1) && (Wfac1 < 0))
{
hd = 1;
}
for (j = 0; j < s->hist->nbin[0]; j++)
{
- double x = Wfac1*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
- double pxdx = s->hist->bin[0][j]*normdx; /* p(x)dx */
- W_ab += pxdx*x;
+ double x = Wfac1 * ((j + s->hist->x0[0]) + 0.5) * dx; /*bin ctr*/
+ double pxdx = s->hist->bin[0][j] * normdx; /* p(x)dx */
+ W_ab += pxdx * x;
}
}
}
for (i = 0; i < cb->nsamples; i++)
{
- samples_t *s = cb->s[i];
- sample_range_t *r = &(cb->r[i]);
+ samples_t* s = cb->s[i];
+ sample_range_t* r = &(cb->r[i]);
if (r->use)
{
if (!s->hist)
{
for (j = r->start; j < r->end; j++)
{
- W_ba += Wfac1*s->du[j];
+ W_ba += Wfac1 * s->du[j];
}
}
else
double normdx = 1.;
double dx;
int hd = 0; /* histogram direction */
- if ( (s->hist->nhist > 1) && (Wfac2 < 0) )
+ if ((s->hist->nhist > 1) && (Wfac2 < 0))
{
hd = 1;
}
for (j = 0; j < s->hist->nbin[0]; j++)
{
- double x = Wfac1*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
- double pxdx = s->hist->bin[0][j]*normdx; /* p(x)dx */
- W_ba += pxdx*x;
+ double x = Wfac1 * ((j + s->hist->x0[0]) + 0.5) * dx; /*bin ctr*/
+ double pxdx = s->hist->bin[0][j] * normdx; /* p(x)dx */
+ W_ba += pxdx * x;
}
}
}
*sb = (W_ba + dg);
}
-static void calc_dg_stddev(sample_coll_t *ca, sample_coll_t *cb,
- double temp, double dg, double *stddev)
+static void calc_dg_stddev(sample_coll_t* ca, sample_coll_t* cb, double temp, double dg, double* stddev)
{
int i, j;
double M;
double Wfac1, Wfac2;
double n1, n2;
- kT = BOLTZ*temp;
- beta = 1/kT;
+ kT = BOLTZ * temp;
+ beta = 1 / kT;
/* count the numbers of samples */
n1 = ca->ntot;
{
/* we're using dhdl, so delta_lambda needs to be a
multiplication factor. */
- double delta_lambda = lambda_vec_abs_diff(cb->native_lambda,
- ca->native_lambda);
- Wfac1 = beta*delta_lambda;
- Wfac2 = -beta*delta_lambda;
+ double delta_lambda = lambda_vec_abs_diff(cb->native_lambda, ca->native_lambda);
+ Wfac1 = beta * delta_lambda;
+ Wfac2 = -beta * delta_lambda;
}
- M = std::log(n1/n2);
+ M = std::log(n1 / n2);
/* calculate average in both directions */
for (i = 0; i < ca->nsamples; i++)
{
- samples_t *s = ca->s[i];
- sample_range_t *r = &(ca->r[i]);
+ samples_t* s = ca->s[i];
+ sample_range_t* r = &(ca->r[i]);
if (r->use)
{
if (!s->hist)
{
for (j = r->start; j < r->end; j++)
{
- sigmafact += 1./(2. + 2.*std::cosh((M + Wfac1*s->du[j] - dg)));
+ sigmafact += 1. / (2. + 2. * std::cosh((M + Wfac1 * s->du[j] - dg)));
}
}
else
double normdx = 1.;
double dx;
int hd = 0; /* histogram direction */
- if ( (s->hist->nhist > 1) && (Wfac1 < 0) )
+ if ((s->hist->nhist > 1) && (Wfac1 < 0))
{
hd = 1;
}
for (j = 0; j < s->hist->nbin[0]; j++)
{
- double x = Wfac1*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
- double pxdx = s->hist->bin[0][j]*normdx; /* p(x)dx */
+ double x = Wfac1 * ((j + s->hist->x0[0]) + 0.5) * dx; /*bin ctr*/
+ double pxdx = s->hist->bin[0][j] * normdx; /* p(x)dx */
- sigmafact += pxdx/(2. + 2.*std::cosh((M + x - dg)));
+ sigmafact += pxdx / (2. + 2. * std::cosh((M + x - dg)));
}
}
}
}
for (i = 0; i < cb->nsamples; i++)
{
- samples_t *s = cb->s[i];
- sample_range_t *r = &(cb->r[i]);
+ samples_t* s = cb->s[i];
+ sample_range_t* r = &(cb->r[i]);
if (r->use)
{
if (!s->hist)
{
for (j = r->start; j < r->end; j++)
{
- sigmafact += 1./(2. + 2.*std::cosh((M - Wfac2*s->du[j] - dg)));
+ sigmafact += 1. / (2. + 2. * std::cosh((M - Wfac2 * s->du[j] - dg)));
}
}
else
double normdx = 1.;
double dx;
int hd = 0; /* histogram direction */
- if ( (s->hist->nhist > 1) && (Wfac2 < 0) )
+ if ((s->hist->nhist > 1) && (Wfac2 < 0))
{
hd = 1;
}
for (j = 0; j < s->hist->nbin[0]; j++)
{
- double x = Wfac2*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
- double pxdx = s->hist->bin[0][j]*normdx; /* p(x)dx */
+ double x = Wfac2 * ((j + s->hist->x0[0]) + 0.5) * dx; /*bin ctr*/
+ double pxdx = s->hist->bin[0][j] * normdx; /* p(x)dx */
- sigmafact += pxdx/(2. + 2.*std::cosh((M - x - dg)));
+ sigmafact += pxdx / (2. + 2. * std::cosh((M - x - dg)));
}
}
}
/* Eq. 10 from
Shirts, Bair, Hooker & Pande, Phys. Rev. Lett 91, 140601 (2003): */
- *stddev = std::sqrt(((1.0/sigmafact) - ( (n1+n2)/n1 + (n1+n2)/n2 )));
+ *stddev = std::sqrt(((1.0 / sigmafact) - ((n1 + n2) / n1 + (n1 + n2) / n2)));
}
-
-static void calc_bar(barres_t *br, double tol,
- int npee_min, int npee_max, gmx_bool *bEE,
- double *partsum)
+static void calc_bar(barres_t* br, double tol, int npee_min, int npee_max, gmx_bool* bEE, double* partsum)
{
- int npee, p;
- double dg_sig2, sa_sig2, sb_sig2, stddev_sig2; /* intermediate variance values
- for calculated quantities */
+ int npee, p;
+ double dg_sig2, sa_sig2, sb_sig2, stddev_sig2; /* intermediate variance values
+ for calculated quantities */
double temp = br->a->temp;
int i;
double dg_min, dg_max;
dg_min = calc_bar_lowlevel(br->a, br->b, temp, tol, -1);
dg_max = calc_bar_lowlevel(br->a, br->b, temp, tol, 1);
- if (std::abs(dg_max - dg_min) > GMX_REAL_EPS*10)
+ if (std::abs(dg_max - dg_min) > GMX_REAL_EPS * 10)
{
/* the histogram range error is the biggest of the differences
between the best estimate and the extremes */
}
calc_rel_entropy(br->a, br->b, temp, br->dg, &(br->sa), &(br->sb));
- calc_dg_stddev(br->a, br->b, temp, br->dg, &(br->dg_stddev) );
+ calc_dg_stddev(br->a, br->b, temp, br->dg, &(br->dg_stddev));
dg_sig2 = 0;
sa_sig2 = 0;
sample_coll_t ca, cb;
/* initialize the samples */
- sample_coll_init(&ca, br->a->native_lambda, br->a->foreign_lambda,
- br->a->temp);
- sample_coll_init(&cb, br->b->native_lambda, br->b->foreign_lambda,
- br->b->temp);
+ sample_coll_init(&ca, br->a->native_lambda, br->a->foreign_lambda, br->a->temp);
+ sample_coll_init(&cb, br->b->native_lambda, br->b->foreign_lambda, br->b->temp);
for (npee = npee_min; npee <= npee_max; npee++)
{
if (!cac || !cbc)
{
- printf("WARNING: histogram number incompatible with block number for averaging: can't do error estimate\n");
+ printf("WARNING: histogram number incompatible with block number for "
+ "averaging: can't do error estimate\n");
*bEE = FALSE;
if (cac)
{
return;
}
- dgp = calc_bar_lowlevel(&ca, &cb, temp, tol, 0);
- dgs += dgp;
- dgs2 += dgp*dgp;
+ dgp = calc_bar_lowlevel(&ca, &cb, temp, tol, 0);
+ dgs += dgp;
+ dgs2 += dgp * dgp;
- partsum[npee*(npee_max+1)+p] += dgp;
+ partsum[npee * (npee_max + 1) + p] += dgp;
calc_rel_entropy(&ca, &cb, temp, dgp, &sac, &sbc);
- dsa += sac;
- dsa2 += sac*sac;
- dsb += sbc;
- dsb2 += sbc*sbc;
- calc_dg_stddev(&ca, &cb, temp, dgp, &stddevc );
+ dsa += sac;
+ dsa2 += sac * sac;
+ dsb += sbc;
+ dsb2 += sbc * sbc;
+ calc_dg_stddev(&ca, &cb, temp, dgp, &stddevc);
- dstddev += stddevc;
- dstddev2 += stddevc*stddevc;
+ dstddev += stddevc;
+ dstddev2 += stddevc * stddevc;
sample_coll_destroy(&ca);
sample_coll_destroy(&cb);
}
- dgs /= npee;
- dgs2 /= npee;
- dg_sig2 += (dgs2-dgs*dgs)/(npee-1);
+ dgs /= npee;
+ dgs2 /= npee;
+ dg_sig2 += (dgs2 - dgs * dgs) / (npee - 1);
- dsa /= npee;
- dsa2 /= npee;
- dsb /= npee;
- dsb2 /= npee;
- sa_sig2 += (dsa2-dsa*dsa)/(npee-1);
- sb_sig2 += (dsb2-dsb*dsb)/(npee-1);
+ dsa /= npee;
+ dsa2 /= npee;
+ dsb /= npee;
+ dsb2 /= npee;
+ sa_sig2 += (dsa2 - dsa * dsa) / (npee - 1);
+ sb_sig2 += (dsb2 - dsb * dsb) / (npee - 1);
- dstddev /= npee;
- dstddev2 /= npee;
- stddev_sig2 += (dstddev2-dstddev*dstddev)/(npee-1);
+ dstddev /= npee;
+ dstddev2 /= npee;
+ stddev_sig2 += (dstddev2 - dstddev * dstddev) / (npee - 1);
}
- br->dg_err = std::sqrt(dg_sig2/(npee_max - npee_min + 1));
- br->sa_err = std::sqrt(sa_sig2/(npee_max - npee_min + 1));
- br->sb_err = std::sqrt(sb_sig2/(npee_max - npee_min + 1));
- br->dg_stddev_err = std::sqrt(stddev_sig2/(npee_max - npee_min + 1));
+ br->dg_err = std::sqrt(dg_sig2 / (npee_max - npee_min + 1));
+ br->sa_err = std::sqrt(sa_sig2 / (npee_max - npee_min + 1));
+ br->sb_err = std::sqrt(sb_sig2 / (npee_max - npee_min + 1));
+ br->dg_stddev_err = std::sqrt(stddev_sig2 / (npee_max - npee_min + 1));
}
}
-static double bar_err(int nbmin, int nbmax, const double *partsum)
+static double bar_err(int nbmin, int nbmax, const double* partsum)
{
int nb, b;
double svar, s, s2, dg;
s2 = 0;
for (b = 0; b < nb; b++)
{
- dg = partsum[nb*(nbmax+1)+b];
- s += dg;
- s2 += dg*dg;
+ dg = partsum[nb * (nbmax + 1) + b];
+ s += dg;
+ s2 += dg * dg;
}
- s /= nb;
- s2 /= nb;
- svar += (s2 - s*s)/(nb - 1);
+ s /= nb;
+ s2 /= nb;
+ svar += (s2 - s * s) / (nb - 1);
}
- return std::sqrt(svar/(nbmax + 1 - nbmin));
+ return std::sqrt(svar / (nbmax + 1 - nbmin));
}
first non-space value found after that. Returns NULL if the string
ends before that.
*/
-static const char *find_value(const char *str)
+static const char* find_value(const char* str)
{
gmx_bool name_end_found = FALSE;
/* first find the end of the name */
if (!name_end_found)
{
- if (std::isspace(*str) || (*str == '=') )
+ if (std::isspace(*str) || (*str == '='))
{
name_end_found = TRUE;
}
}
else
{
- if (!( std::isspace(*str) || (*str == '=') ))
+ if (!(std::isspace(*str) || (*str == '=')))
{
return str;
}
/* 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 = nullptr; /* start of the component name, or NULL
- if not in a value */
- char *strtod_end;
- gmx_bool OK = TRUE;
+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 = nullptr; /* start of the component name, or NULL
+ if not in a value */
+ char* strtod_end;
+ gmx_bool OK = TRUE;
if (end)
{
{
if (initialize_lc)
{
- lambda_components_add(lc_out, val_start,
- (str-val_start));
+ lambda_components_add(lc_out, val_start, (str - val_start));
}
else
{
- if (!lambda_components_check(lc_out, n, val_start,
- (str-val_start)))
+ if (!lambda_components_check(lc_out, n, val_start, (str - val_start)))
{
return FALSE;
}
lv->val[n] = strtod(val_start, &strtod_end);
if (val_start == strtod_end)
{
- gmx_fatal(FARGS,
- "Error reading lambda vector in %s", fn);
+ gmx_fatal(FARGS, "Error reading lambda vector in %s", fn);
}
}
/* reset for the next identifier */
}
else
{
- gmx_fatal(FARGS, "Incomplete lambda vector data in %s",
- fn);
+ gmx_fatal(FARGS, "Incomplete lambda vector data in %s", fn);
return FALSE;
}
-
}
}
}
}
/* 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)
+static gmx_bool read_lambda_components(const char* str, lambda_components_t* lc, const char** end, const char* fn)
{
return read_lambda_compvec(str, nullptr, nullptr, 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)
+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, nullptr, end, fn);
}
-
/* deduce lambda value from legend.
fn = the file name
legend = the legend string
lam = the initialized lambda vector
returns whether to use the data in this set.
*/
-static gmx_bool legend2lambda(const char *fn,
- const char *legend,
- lambda_vec_t *lam)
+static gmx_bool legend2lambda(const char* fn, const char* legend, lambda_vec_t* lam)
{
- const char *ptr = nullptr, *ptr2 = nullptr;
- gmx_bool ok = FALSE;
- gmx_bool bdhdl = FALSE;
- const char *tostr = " to ";
+ const char *ptr = nullptr, *ptr2 = nullptr;
+ gmx_bool ok = FALSE;
+ gmx_bool bdhdl = FALSE;
+ const char* tostr = " to ";
if (legend == nullptr)
{
ptr = ptr2;
ptr2++;
}
- }
- while (ptr2 != nullptr && *ptr2 != '\0');
+ } while (ptr2 != nullptr && *ptr2 != '\0');
if (ptr)
{
- ptr += std::strlen(tostr)-1; /* and advance past that 'to' */
+ ptr += std::strlen(tostr) - 1; /* and advance past that 'to' */
}
else
{
else
{
int dhdl_index;
- const char *end;
+ const char* end;
ptr = std::strrchr(legend, '=');
end = ptr;
}
}
ptr++;
- dhdl_index = lambda_components_find(lam->lc, ptr, (end-ptr));
+ dhdl_index = lambda_components_find(lam->lc, ptr, (end - ptr));
if (dhdl_index < 0)
{
char buf[STRLEN];
- std::strncpy(buf, ptr, (end-ptr));
- buf[(end-ptr)] = '\0';
- gmx_fatal(FARGS,
- "Did not find lambda component for '%s' in %s",
- buf, fn);
+ std::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);
+ gmx_fatal(FARGS, "dhdl without component name with >1 lambda component in %s", fn);
}
dhdl_index = 0;
}
return TRUE;
}
-static gmx_bool subtitle2lambda(const char *subtitle, xvg_t *ba, const char *fn,
- lambda_components_t *lc)
+static gmx_bool subtitle2lambda(const char* subtitle, xvg_t* ba, const char* fn, lambda_components_t* lc)
{
gmx_bool bFound;
- const char *ptr;
- char *end;
+ const char* ptr;
+ char* end;
double native_lambda;
bFound = FALSE;
if (ptr)
{
int index = -1;
- const char *val_end;
+ const char* val_end;
/* the new 4.6 style lambda vectors */
ptr = find_value(ptr);
ptr++;
if (*ptr == '\0')
{
- gmx_fatal(FARGS,
- "Incomplete lambda vector component data in %s", fn);
+ 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);
+ gmx_fatal(FARGS, "lambda vector components in %s don't match those previously read", fn);
}
ptr = find_value(val_end);
if (!ptr)
}
if (ptr != nullptr)
{
- bFound = (sscanf(ptr+1, "%lf", &(native_lambda)) == 1);
+ 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);
+ "lambda vector components in %s don't match those previously read", fn);
}
}
else
return bFound;
}
-static void read_bar_xvg_lowlevel(const char *fn, const real *temp, xvg_t *ba,
- lambda_components_t *lc)
+static void read_bar_xvg_lowlevel(const char* fn, const 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];
+ int i;
+ char * subtitle, **legend, *ptr;
+ int np;
+ gmx_bool native_lambda_read = FALSE;
+ char buf[STRLEN];
xvg_init(ba);
gmx_fatal(FARGS, "File %s contains no usable data.", fn);
}
/* Reorder the data */
- ba->t = ba->y[0];
+ ba->t = ba->y[0];
for (i = 1; i < ba->nset; i++)
{
- ba->y[i-1] = ba->y[i];
+ ba->y[i - 1] = ba->y[i];
}
ba->nset--;
{
if (ba->temp <= 0)
{
- gmx_fatal(FARGS, "Found temperature of %f in file '%s'",
- ba->temp, fn);
+ gmx_fatal(FARGS, "Found temperature of %f in file '%s'", ba->temp, fn);
}
}
}
{
if (*temp <= 0)
{
- gmx_fatal(FARGS, "Did not find a temperature in the subtitle in file '%s', use the -temp option of [TT]gmx bar[tt]", fn);
+ gmx_fatal(FARGS,
+ "Did not find a temperature in the subtitle in file '%s', use the -temp "
+ "option of [TT]gmx bar[tt]",
+ fn);
}
ba->temp = *temp;
}
}
else
{
- gmx_fatal(FARGS, "File %s contains multiple sets but no legends, can not determine the lambda values", fn);
+ gmx_fatal(FARGS,
+ "File %s contains multiple sets but no legends, can not determine the lambda "
+ "values",
+ fn);
}
}
else
{
- for (i = 0; i < ba->nset; )
+ for (i = 0; i < ba->nset;)
{
/* Read lambda from the legend */
- lambda_vec_init( &(ba->lambda[i]), lc );
- lambda_vec_copy( &(ba->lambda[i]), &(ba->native_lambda));
+ lambda_vec_init(&(ba->lambda[i]), lc);
+ lambda_vec_copy(&(ba->lambda[i]), &(ba->native_lambda));
gmx_bool use = legend2lambda(fn, legend[i], &(ba->lambda[i]));
if (use)
{
{
int j;
printf("%s: Ignoring set '%s'.\n", fn, legend[i]);
- for (j = i+1; j < ba->nset; j++)
+ for (j = i + 1; j < ba->nset; j++)
{
- ba->y[j-1] = ba->y[j];
- legend[j-1] = legend[j];
+ ba->y[j - 1] = ba->y[j];
+ legend[j - 1] = legend[j];
}
ba->nset--;
}
if (legend != nullptr)
{
- for (i = 0; i < ba->nset-1; i++)
+ for (i = 0; i < ba->nset - 1; i++)
{
sfree(legend[i]);
}
}
}
-static void read_bar_xvg(const char *fn, real *temp, sim_data_t *sd)
+static void read_bar_xvg(const char* fn, real* temp, sim_data_t* sd)
{
- xvg_t *barsim;
- samples_t *s;
+ xvg_t* barsim;
+ samples_t* s;
int i;
snew(barsim, 1);
gmx_fatal(FARGS, "File '%s' contains fewer than two columns", fn);
}
- if (!gmx_within_tol(*temp, barsim->temp, GMX_FLOAT_EPS) && (*temp > 0) )
+ if (!gmx_within_tol(*temp, barsim->temp, GMX_FLOAT_EPS) && (*temp > 0))
{
gmx_fatal(FARGS, "Temperature in file %s different from earlier files or setting\n", fn);
}
snew(s, barsim->nset);
for (i = 0; i < barsim->nset; i++)
{
- samples_init(s+i, &(barsim->native_lambda), &(barsim->lambda[i]),
- barsim->temp, lambda_vec_same(&(barsim->native_lambda),
- &(barsim->lambda[i])),
- fn);
+ 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_data_list_insert_sample(sd->lb, s+i);
+ lambda_data_list_insert_sample(sd->lb, s + i);
}
{
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);
+ 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("\n\n");
}
-static void read_edr_rawdh_block(samples_t **smp, int *ndu, t_enxblock *blk,
- double start_time, double delta_time,
- lambda_vec_t *native_lambda, double temp,
- double *last_t, const char *filename)
+static void read_edr_rawdh_block(samples_t** smp,
+ int* ndu,
+ t_enxblock* blk,
+ double start_time,
+ double delta_time,
+ lambda_vec_t* native_lambda,
+ double temp,
+ double* last_t,
+ const char* filename)
{
int i, j;
- lambda_vec_t *foreign_lambda;
+ lambda_vec_t* foreign_lambda;
int type;
- samples_t *s; /* convenience pointer */
+ samples_t* s; /* convenience pointer */
int startj;
/* check the block types etc. */
- if ( (blk->nsub < 3) ||
- (blk->sub[0].type != xdr_datatype_int) ||
- (blk->sub[1].type != xdr_datatype_double) ||
- (
- (blk->sub[2].type != xdr_datatype_float) &&
- (blk->sub[2].type != xdr_datatype_double)
- ) ||
- (blk->sub[0].nr < 1) ||
- (blk->sub[1].nr < 1) )
+ if ((blk->nsub < 3) || (blk->sub[0].type != xdr_datatype_int) || (blk->sub[1].type != xdr_datatype_double)
+ || ((blk->sub[2].type != xdr_datatype_float) && (blk->sub[2].type != xdr_datatype_double))
+ || (blk->sub[0].nr < 1) || (blk->sub[1].nr < 1))
{
- gmx_fatal(FARGS,
- "Unexpected/corrupted block data in file %s around time %f.",
- filename, start_time);
+ gmx_fatal(FARGS, "Unexpected/corrupted block data in file %s around time %f.", filename, start_time);
}
snew(foreign_lambda, 1);
{
/* initialize the samples structure if it's empty. */
snew(*smp, 1);
- samples_init(*smp, native_lambda, foreign_lambda, temp,
- type == dhbtDHDL, filename);
+ samples_init(*smp, native_lambda, foreign_lambda, temp, type == dhbtDHDL, filename);
(*smp)->start_time = start_time;
(*smp)->delta_time = delta_time;
}
s = *smp;
/* now double check */
- if (!lambda_vec_same(s->foreign_lambda, foreign_lambda) )
+ if (!lambda_vec_same(s->foreign_lambda, foreign_lambda))
{
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=%f.",
- filename, start_time);
+ gmx_fatal(FARGS, "Corrupted data in file %s around t=%f.", filename, start_time);
}
/* make room for the data */
if (gmx::index(s->ndu_alloc) < s->ndu + blk->sub[2].nr)
{
- s->ndu_alloc += (s->ndu_alloc < static_cast<size_t>(blk->sub[2].nr)) ?
- blk->sub[2].nr*2 : s->ndu_alloc;
+ s->ndu_alloc += (s->ndu_alloc < static_cast<size_t>(blk->sub[2].nr)) ? blk->sub[2].nr * 2
+ : s->ndu_alloc;
srenew(s->du_alloc, s->ndu_alloc);
s->du = s->du_alloc;
}
- startj = s->ndu;
- s->ndu += blk->sub[2].nr;
+ startj = s->ndu;
+ s->ndu += blk->sub[2].nr;
s->ntot += blk->sub[2].nr;
- *ndu = blk->sub[2].nr;
+ *ndu = blk->sub[2].nr;
/* and copy the data*/
for (j = 0; j < blk->sub[2].nr; j++)
{
if (blk->sub[2].type == xdr_datatype_float)
{
- s->du[startj+j] = blk->sub[2].fval[j];
+ s->du[startj + j] = blk->sub[2].fval[j];
}
else
{
- s->du[startj+j] = blk->sub[2].dval[j];
+ s->du[startj + j] = blk->sub[2].dval[j];
}
}
- if (start_time + blk->sub[2].nr*delta_time > *last_t)
+ if (start_time + blk->sub[2].nr * delta_time > *last_t)
{
- *last_t = start_time + blk->sub[2].nr*delta_time;
+ *last_t = start_time + blk->sub[2].nr * delta_time;
}
}
-static samples_t *read_edr_hist_block(int *nsamples, t_enxblock *blk,
- double start_time, double delta_time,
- lambda_vec_t *native_lambda, double temp,
- double *last_t, const char *filename)
+static samples_t* read_edr_hist_block(int* nsamples,
+ t_enxblock* blk,
+ double start_time,
+ double delta_time,
+ lambda_vec_t* native_lambda,
+ double temp,
+ double* last_t,
+ const char* filename)
{
int i, j;
- samples_t *s;
+ samples_t* s;
int nhist;
- lambda_vec_t *foreign_lambda;
+ lambda_vec_t* foreign_lambda;
int type;
int nbins[2];
/* check the block types etc. */
- if ( (blk->nsub < 2) ||
- (blk->sub[0].type != xdr_datatype_double) ||
- (blk->sub[1].type != xdr_datatype_int64) ||
- (blk->sub[0].nr < 2) ||
- (blk->sub[1].nr < 2) )
+ if ((blk->nsub < 2) || (blk->sub[0].type != xdr_datatype_double)
+ || (blk->sub[1].type != xdr_datatype_int64) || (blk->sub[0].nr < 2) || (blk->sub[1].nr < 2))
{
- gmx_fatal(FARGS,
- "Unexpected/corrupted block data in file %s around time %f",
- filename, start_time);
+ gmx_fatal(FARGS, "Unexpected/corrupted block data in file %s around time %f", filename, start_time);
}
- nhist = blk->nsub-2;
+ nhist = blk->nsub - 2;
if (nhist == 0)
{
return nullptr;
}
if (nhist > 2)
{
- gmx_fatal(FARGS,
- "Unexpected/corrupted block data in file %s around time %f",
- filename, start_time);
+ gmx_fatal(FARGS, "Unexpected/corrupted block data in file %s around time %f", filename, start_time);
}
snew(s, 1);
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);
+ 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];
+ foreign_lambda->val[i] = blk->sub[0].dval[i + 2];
}
}
}
{
if (blk->sub[1].nr < 3 + nhist)
{
- gmx_fatal(FARGS,
- "Missing derivative coord in multi-component file %s",
- filename);
+ gmx_fatal(FARGS, "Missing derivative coord in multi-component file %s", filename);
}
foreign_lambda->dhdl = blk->sub[1].lval[2 + nhist];
}
}
}
- samples_init(s, native_lambda, foreign_lambda, temp, type == dhbtDHDL,
- filename);
+ samples_init(s, native_lambda, foreign_lambda, temp, type == dhbtDHDL, filename);
snew(s->hist, 1);
for (i = 0; i < nhist; i++)
{
- nbins[i] = blk->sub[i+2].nr;
+ nbins[i] = blk->sub[i + 2].nr;
}
hist_init(s->hist, nhist, nbins);
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)
{
for (i = 0; i < nhist; i++)
{
- int64_t sum = 0;
+ int64_t sum = 0;
for (j = 0; j < s->hist->nbin[i]; j++)
{
- int binv = static_cast<int>(blk->sub[i+2].ival[j]);
+ int binv = static_cast<int>(blk->sub[i + 2].ival[j]);
s->hist->bin[i][j] = binv;
- sum += binv;
-
+ sum += binv;
}
if (i == 0)
{
{
if (s->ntot != sum)
{
- gmx_fatal(FARGS, "Histogram counts don't match in %s",
- filename);
+ gmx_fatal(FARGS, "Histogram counts don't match in %s", filename);
}
}
}
- if (start_time + s->hist->sum*delta_time > *last_t)
+ if (start_time + s->hist->sum * delta_time > *last_t)
{
- *last_t = start_time + s->hist->sum*delta_time;
+ *last_t = start_time + s->hist->sum * delta_time;
}
return s;
}
-static void read_barsim_edr(const char *fn, real *temp, sim_data_t *sd)
+static void read_barsim_edr(const char* fn, real* temp, sim_data_t* sd)
{
int i, j;
ener_file_t fp;
- t_enxframe *fr;
+ t_enxframe* fr;
int nre;
- gmx_enxnm_t *enm = nullptr;
+ gmx_enxnm_t* enm = nullptr;
double first_t = -1;
double last_t = -1;
- samples_t **samples_rawdh = nullptr; /* contains samples for raw delta_h */
- int *nhists = nullptr; /* array to keep count & print at end */
- int *npts = nullptr; /* array to keep count & print at end */
- lambda_vec_t **lambdas = nullptr; /* array to keep count & print at end */
- lambda_vec_t *native_lambda;
+ samples_t** samples_rawdh = nullptr; /* contains samples for raw delta_h */
+ int* nhists = nullptr; /* array to keep count & print at end */
+ int* npts = nullptr; /* array to keep count & print at end */
+ lambda_vec_t** lambdas = nullptr; /* array to keep count & print at end */
+ lambda_vec_t* native_lambda;
int nsamples = 0;
lambda_vec_t start_lambda;
while (do_enx(fp, fr))
{
/* count the data blocks */
- int nblocks_raw = 0;
- int nblocks_hist = 0;
- int nlam = 0;
- int k;
+ int nblocks_raw = 0;
+ int nblocks_hist = 0;
+ int nlam = 0;
+ int k;
/* DHCOLL block information: */
double start_time = 0, delta_time = 0, old_start_lambda = 0, delta_lambda = 0;
- double rtemp = 0;
+ double rtemp = 0;
/* count the blocks and handle collection information: */
for (i = 0; i < fr->nblock; i++)
if (fr->block[i].id == enxDHCOLL)
{
nlam++;
- if ( (fr->block[i].nsub < 1) ||
- (fr->block[i].sub[0].type != xdr_datatype_double) ||
- (fr->block[i].sub[0].nr < 5))
+ 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 in file %s", fn);
}
{
gmx_fatal(FARGS, "Lambda values not constant in %s: can't apply BAR method", fn);
}
- if ( ( *temp != rtemp) && (*temp > 0) )
+ if ((*temp != rtemp) && (*temp > 0))
{
- gmx_fatal(FARGS, "Temperature in file %s different from earlier files or setting\n", fn);
+ gmx_fatal(FARGS,
+ "Temperature in file %s different from earlier files or setting\n", fn);
}
*temp = rtemp;
if (!lambda_components_check(&(sd->lc), 0, "", 0))
{
gmx_fatal(FARGS,
- "lambda vector components in %s don't match those previously read",
+ "lambda vector components in %s don't match those previously "
+ "read",
fn);
}
}
gmx_bool check = (sd->lc.N > 0);
if (fr->block[i].nsub < 2)
{
- gmx_fatal(FARGS,
- "No lambda vector, but start_lambda=%f\n",
- old_start_lambda);
+ gmx_fatal(FARGS, "No lambda vector, but start_lambda=%f\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]];
+ 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,
- std::strlen(name));
+ lambda_components_check(&(sd->lc), j, name, std::strlen(name));
}
else
{
- lambda_components_add(&(sd->lc), name,
- std::strlen(name));
+ lambda_components_add(&(sd->lc), name, std::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];
+ start_lambda.val[j] = fr->block[i].sub[0].dval[5 + j];
}
}
if (first_t < 0)
/*native_lambda=start_lambda;*/
lambda_vec_init(native_lambda, &(sd->lc));
lambda_vec_copy(native_lambda, &start_lambda);
- nsamples = nblocks_raw+nblocks_hist;
+ nsamples = nblocks_raw + nblocks_hist;
snew(nhists, nsamples);
snew(npts, nsamples);
snew(lambdas, nsamples);
// nsamples > 0 means this is NOT the first iteration
/* check the native lambda */
- if (!lambda_vec_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 %f, and becomes %f at time %f",
+ gmx_fatal(FARGS,
+ "Native lambda not constant in file %s: started at %f, and becomes %f at "
+ "time %f",
fn, native_lambda->val[0], start_lambda.val[0], start_time);
}
/* check the number of samples against the previous number */
- if ( ((nblocks_raw+nblocks_hist) != nsamples) || (nlam != 1 ) )
+ if (((nblocks_raw + nblocks_hist) != nsamples) || (nlam != 1))
{
- gmx_fatal(FARGS, "Unexpected block count in %s: was %d, now %d\n",
- fn, nsamples+1, nblocks_raw+nblocks_hist+nlam);
+ gmx_fatal(FARGS, "Unexpected block count in %s: was %d, now %d\n", fn, nsamples + 1,
+ nblocks_raw + nblocks_hist + nlam);
}
/* check whether last iterations's end time matches with
the currrent start time */
- if ( (std::abs(last_t - start_time) > 2*delta_time) && last_t >= 0)
+ if ((std::abs(last_t - start_time) > 2 * delta_time) && last_t >= 0)
{
/* it didn't. We need to store our samples and reallocate */
for (i = 0; i < nsamples; i++)
if (samples_rawdh[i])
{
/* insert it into the existing list */
- lambda_data_list_insert_sample(sd->lb,
- 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] = nullptr;
if (type == dhbtDH || type == dhbtDHDL)
{
int ndu;
- read_edr_rawdh_block(&(samples_rawdh[k]),
- &ndu,
- &(fr->block[i]),
- start_time, delta_time,
- native_lambda, rtemp,
- &last_t, fn);
+ 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])
{
{
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);
+ 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)
{
/* and insert the new sample immediately */
for (j = 0; j < nb; j++)
{
- lambda_data_list_insert_sample(sd->lb, s+j);
+ lambda_data_list_insert_sample(sd->lb, s + j);
}
}
}
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);
+ printf("%s: %.1f - %.1f; lambda = %s\n foreign lambdas:\n", fn, first_t, last_t, buf);
for (i = 0; i < nsamples; i++)
{
if (lambdas[i])
}
-int gmx_bar(int argc, char *argv[])
+int gmx_bar(int argc, char* argv[])
{
- static const char *desc[] = {
+ static const char* desc[] = {
"[THISMODULE] calculates free energy difference estimates through ",
"Bennett's acceptance ratio method (BAR). It also automatically",
"adds series of individual free energies obtained with BAR into",
"These can contain either lists of energy differences (see the ",
"[REF].mdp[ref] option [TT]separate_dhdl_file[tt]), or a series of ",
"histograms (see the [REF].mdp[ref] options [TT]dh_hist_size[tt] and ",
- "[TT]dh_hist_spacing[tt]).", "The temperature and [GRK]lambda[grk] ",
+ "[TT]dh_hist_spacing[tt]).",
+ "The temperature and [GRK]lambda[grk] ",
"values are automatically deduced from the [TT]ener.edr[tt] file.[PAR]",
"In addition to the [REF].mdp[ref] option [TT]foreign_lambda[tt], ",
"[TT]-oh[tt] option to write series of histograms, together with the ",
"[TT]-nbin[tt] option.[PAR]"
};
- static real begin = 0, end = -1, temp = -1;
- int nd = 2, nbmin = 5, nbmax = 5;
- int nbin = 100;
- gmx_bool use_dhdl = FALSE;
- t_pargs pa[] = {
- { "-b", FALSE, etREAL, {&begin}, "Begin time for BAR" },
- { "-e", FALSE, etREAL, {&end}, "End time for BAR" },
- { "-temp", FALSE, etREAL, {&temp}, "Temperature (K)" },
- { "-prec", FALSE, etINT, {&nd}, "The number of digits after the decimal point" },
- { "-nbmin", FALSE, etINT, {&nbmin}, "Minimum number of blocks for error estimation" },
- { "-nbmax", FALSE, etINT, {&nbmax}, "Maximum number of blocks for error estimation" },
- { "-nbin", FALSE, etINT, {&nbin}, "Number of bins for histogram output"},
- { "-extp", FALSE, etBOOL, {&use_dhdl}, "Whether to linearly extrapolate dH/dl values to use as energies"}
+ static real begin = 0, end = -1, temp = -1;
+ int nd = 2, nbmin = 5, nbmax = 5;
+ int nbin = 100;
+ gmx_bool use_dhdl = FALSE;
+ t_pargs pa[] = {
+ { "-b", FALSE, etREAL, { &begin }, "Begin time for BAR" },
+ { "-e", FALSE, etREAL, { &end }, "End time for BAR" },
+ { "-temp", FALSE, etREAL, { &temp }, "Temperature (K)" },
+ { "-prec", FALSE, etINT, { &nd }, "The number of digits after the decimal point" },
+ { "-nbmin", FALSE, etINT, { &nbmin }, "Minimum number of blocks for error estimation" },
+ { "-nbmax", FALSE, etINT, { &nbmax }, "Maximum number of blocks for error estimation" },
+ { "-nbin", FALSE, etINT, { &nbin }, "Number of bins for histogram output" },
+ { "-extp",
+ FALSE,
+ etBOOL,
+ { &use_dhdl },
+ "Whether to linearly extrapolate dH/dl values to use as energies" }
};
- t_filenm fnm[] = {
- { efXVG, "-f", "dhdl", ffOPTRDMULT },
- { efEDR, "-g", "ener", ffOPTRDMULT },
- { efXVG, "-o", "bar", ffOPTWR },
- { efXVG, "-oi", "barint", ffOPTWR },
- { efXVG, "-oh", "histogram", ffOPTWR }
- };
+ t_filenm fnm[] = { { efXVG, "-f", "dhdl", ffOPTRDMULT },
+ { efEDR, "-g", "ener", ffOPTRDMULT },
+ { efXVG, "-o", "bar", ffOPTWR },
+ { efXVG, "-oi", "barint", ffOPTWR },
+ { efXVG, "-oh", "histogram", ffOPTWR } };
#define NFILE asize(fnm)
- int f;
- int nf = 0; /* file counter */
- int nfile_tot; /* total number of input files */
- sim_data_t sim_data; /* the simulation data */
- barres_t *results; /* the results */
- int nresults; /* number of results in results array */
+ int f;
+ int nf = 0; /* file counter */
+ int nfile_tot; /* total number of input files */
+ sim_data_t sim_data; /* the simulation data */
+ barres_t* results; /* the results */
+ int nresults; /* number of results in results array */
- double *partsum;
+ double* partsum;
double prec, dg_tot;
- FILE *fpb, *fpi;
+ FILE * fpb, *fpi;
char dgformat[20], xvg2format[STRLEN], xvg3format[STRLEN];
char buf[STRLEN], buf2[STRLEN];
char ktformat[STRLEN], sktformat[STRLEN];
char kteformat[STRLEN], skteformat[STRLEN];
- gmx_output_env_t *oenv;
+ gmx_output_env_t* oenv;
double kT;
gmx_bool result_OK = TRUE, bEE = TRUE;
- gmx_bool disc_err = FALSE;
- double sum_disc_err = 0.; /* discretization error */
- gmx_bool histrange_err = FALSE;
- double sum_histrange_err = 0.; /* histogram range error */
- double stat_err = 0.; /* statistical error */
+ gmx_bool disc_err = FALSE;
+ double sum_disc_err = 0.; /* discretization error */
+ gmx_bool histrange_err = FALSE;
+ double sum_histrange_err = 0.; /* histogram range error */
+ double stat_err = 0.; /* statistical error */
- if (!parse_common_args(&argc, argv,
- PCA_CAN_VIEW,
- NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
+ if (!parse_common_args(&argc, argv, PCA_CAN_VIEW, NFILE, fnm, asize(pa), pa, asize(desc), desc,
+ 0, nullptr, &oenv))
{
return 0;
}
}
prec = std::pow(10.0, static_cast<double>(-nd));
- snew(partsum, (nbmax+1)*(nbmax+1));
+ snew(partsum, (nbmax + 1) * (nbmax + 1));
nf = 0;
/* read in all files. First xvg files */
- for (const std::string &filenm : xvgFiles)
+ for (const std::string& filenm : xvgFiles)
{
read_bar_xvg(filenm.c_str(), &temp, &sim_data);
nf++;
}
/* then .edr files */
- for (const std::string &filenm : edrFiles)
+ for (const std::string& filenm : edrFiles)
{
- read_barsim_edr(filenm.c_str(), &temp, &sim_data);;
+ read_barsim_edr(filenm.c_str(), &temp, &sim_data);
+
nf++;
}
{
prec = sum_disc_err;
nd = static_cast<int>(std::ceil(-std::log10(prec)));
- printf("WARNING: setting the precision to %g because that is the minimum\n reasonable number, given the expected discretization error.\n", prec);
+ printf("WARNING: setting the precision to %g because that is the minimum\n "
+ "reasonable number, given the expected discretization error.\n",
+ prec);
}
/*sprintf(lamformat,"%%6.3f");*/
- sprintf( dgformat, "%%%d.%df", 3+nd, nd);
+ sprintf(dgformat, "%%%d.%df", 3 + nd, nd);
/* the format strings of the results in kT */
- sprintf( ktformat, "%%%d.%df", 5+nd, nd);
- sprintf( sktformat, "%%%ds", 6+nd);
+ sprintf(ktformat, "%%%d.%df", 5 + nd, nd);
+ sprintf(sktformat, "%%%ds", 6 + nd);
/* the format strings of the errors in kT */
- sprintf( kteformat, "%%%d.%df", 3+nd, nd);
- sprintf( skteformat, "%%%ds", 4+nd);
+ sprintf(kteformat, "%%%d.%df", 3 + nd, nd);
+ sprintf(skteformat, "%%%ds", 4 + nd);
sprintf(xvg2format, "%s %s\n", "%s", dgformat);
sprintf(xvg3format, "%s %s %s\n", "%s", dgformat, dgformat);
-
fpb = nullptr;
if (opt2bSet("-o", NFILE, fnm))
{
sprintf(buf, "%s (%s)", "\\DeltaG", "kT");
- fpb = xvgropen_type(opt2fn("-o", NFILE, fnm), "Free energy differences",
- "\\lambda", buf, exvggtXYDY, oenv);
+ fpb = xvgropen_type(opt2fn("-o", NFILE, fnm), "Free energy differences", "\\lambda", buf,
+ exvggtXYDY, oenv);
}
fpi = nullptr;
if (opt2bSet("-oi", NFILE, fnm))
{
sprintf(buf, "%s (%s)", "\\DeltaG", "kT");
- fpi = xvgropen(opt2fn("-oi", NFILE, fnm), "Free energy integral",
- "\\lambda", buf, oenv);
+ fpi = xvgropen(opt2fn("-oi", NFILE, fnm), "Free energy integral", "\\lambda", buf, oenv);
}
-
if (nbmin > nbmax)
{
nbmin = nbmax;
/* Determine the free energy difference with a factor of 10
* more accuracy than requested for printing.
*/
- calc_bar(&(results[f]), 0.1*prec, nbmin, nbmax,
- &bEE, partsum);
+ calc_bar(&(results[f]), 0.1 * prec, nbmin, nbmax, &bEE, partsum);
- if (results[f].dg_disc_err > prec/10.)
+ if (results[f].dg_disc_err > prec / 10.)
{
disc_err = TRUE;
}
- if (results[f].dg_histrange_err > prec/10.)
+ if (results[f].dg_histrange_err > prec / 10.)
{
histrange_err = TRUE;
}
}
/* print results in kT */
- kT = BOLTZ*temp;
+ kT = BOLTZ * temp;
printf("\nTemperature: %g K\n", temp);
printf("\nDetailed results in kT (see help for explanation):\n\n");
printf("%6s ", " lam_A");
printf("%6s ", " lam_B");
- printf(sktformat, "DG ");
+ printf(sktformat, "DG ");
if (bEE)
{
printf(skteformat, "+/- ");
{
printf(skteformat, "range ");
}
- printf(sktformat, "s_A ");
+ printf(sktformat, "s_A ");
if (bEE)
{
- printf(skteformat, "+/- " );
+ printf(skteformat, "+/- ");
}
- printf(sktformat, "s_B ");
+ printf(sktformat, "s_B ");
if (bEE)
{
- printf(skteformat, "+/- " );
+ printf(skteformat, "+/- ");
}
- printf(sktformat, "stdev ");
+ printf(sktformat, "stdev ");
if (bEE)
{
printf(skteformat, "+/- ");
printf("%s ", buf);
lambda_vec_print_short(results[f].b->native_lambda, buf);
printf("%s ", buf);
- printf(ktformat, results[f].dg);
+ printf(ktformat, results[f].dg);
printf(" ");
if (bEE)
{
printf(kteformat, results[f].dg_histrange_err);
printf(" ");
}
- printf(ktformat, results[f].sa);
+ printf(ktformat, results[f].sa);
printf(" ");
if (bEE)
{
printf(kteformat, results[f].sa_err);
printf(" ");
}
- printf(ktformat, results[f].sb);
+ printf(ktformat, results[f].sb);
printf(" ");
if (bEE)
{
printf(kteformat, results[f].sb_err);
printf(" ");
}
- printf(ktformat, results[f].dg_stddev);
+ printf(ktformat, results[f].dg_stddev);
printf(" ");
if (bEE)
{
printf("\n");
/* Check for negative relative entropy with a 95% certainty. */
- if (results[f].sa < -2*results[f].sa_err ||
- results[f].sb < -2*results[f].sb_err)
+ if (results[f].sa < -2 * results[f].sa_err || results[f].sb < -2 * results[f].sb_err)
{
result_OK = FALSE;
}
/* final results in kJ/mol */
printf("\n\nFinal results in kJ/mol:\n\n");
- dg_tot = 0;
+ dg_tot = 0;
for (f = 0; f < nresults; f++)
{
if (fpb != nullptr)
{
- lambda_vec_print_intermediate(results[f].a->native_lambda,
- results[f].b->native_lambda,
- buf);
+ 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("%s - %s", buf, buf2);
printf(", DG ");
- printf(dgformat, results[f].dg*kT);
+ printf(dgformat, results[f].dg * kT);
if (bEE)
{
printf(" +/- ");
- printf(dgformat, results[f].dg_err*kT);
+ printf(dgformat, results[f].dg_err * kT);
}
if (histrange_err)
{
printf(" (max. range err. = ");
- printf(dgformat, results[f].dg_histrange_err*kT);
+ printf(dgformat, results[f].dg_histrange_err * kT);
printf(")");
- sum_histrange_err += results[f].dg_histrange_err*kT;
+ sum_histrange_err += results[f].dg_histrange_err * kT;
}
printf("\n");
printf("\n");
printf("total ");
lambda_vec_print_short(results[0].a->native_lambda, buf);
- lambda_vec_print_short(results[nresults-1].b->native_lambda, buf2);
+ lambda_vec_print_short(results[nresults - 1].b->native_lambda, buf2);
printf("%s - %s", buf, buf2);
printf(", DG ");
- printf(dgformat, dg_tot*kT);
+ printf(dgformat, dg_tot * kT);
if (bEE)
{
- stat_err = bar_err(nbmin, nbmax, partsum)*kT;
+ stat_err = bar_err(nbmin, nbmax, partsum) * kT;
printf(" +/- ");
printf(dgformat, std::max(std::max(stat_err, sum_disc_err), sum_histrange_err));
}
printf(dgformat, sum_disc_err);
if (bEE && stat_err < sum_disc_err)
{
- printf("WARNING: discretization error (%g) is larger than statistical error.\n Decrease histogram spacing for more accurate results\n", stat_err);
+ printf("WARNING: discretization error (%g) is larger than statistical error.\n "
+ "Decrease histogram spacing for more accurate results\n",
+ stat_err);
}
}
if (histrange_err)
printf(dgformat, sum_histrange_err);
if (bEE && stat_err < sum_histrange_err)
{
- printf("WARNING: histogram range error (%g) is larger than statistical error.\n Increase histogram range for more accurate results\n", stat_err);
+ printf("WARNING: histogram range error (%g) is larger than statistical error.\n "
+ "Increase histogram range for more accurate results\n",
+ stat_err);
}
-
}
printf("\n");
if (fpi != nullptr)
{
- lambda_vec_print_short(results[nresults-1].b->native_lambda, buf);
+ lambda_vec_print_short(results[nresults - 1].b->native_lambda, buf);
fprintf(fpi, xvg2format, buf, dg_tot);
xvgrclose(fpi);
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff