unfolding-refolding, metastable intermediates, etc., but required
reorganizing the command line interface.
+Version 0.7 added tension model inverses, which seems to reduce
+computation time by about a factor of 2. I was expecting more, but
+I'll take what I can get.
+
<<version definition>>=
-#define VERSION "0.6"
+#define VERSION "0.7"
@
\subsection{License}
<<global.h>>=
#ifndef GLOBAL_H
#define GLOBAL_H
+
+#define DOUBLE_PRECISION 1e-12
+
<<environment definition>>
<<one dimensional function definition>>
<<create/destroy definitions>>
{
static int num_active_groups=0;
static one_dim_fn_t **per_group_handlers = NULL;
+ static one_dim_fn_t **per_group_inverse_handlers = NULL;
static void **per_group_data = NULL; /* array of pointers to tension_handler_data_t */
static double last_x;
tension_handler_data_t *tdata;
if (transition != 0 || num_active_groups == 0) { /* setup statics */
/* get new statics, freeing the old ones if they aren't NULL */
- get_active_groups(state_list, &num_active_groups, &per_group_handlers, &per_group_data);
+ get_active_groups(state_list, &num_active_groups, &per_group_handlers, &per_group_inverse_handlers, &per_group_data);
/* fill in the group handlers and params */
#ifdef DEBUG
- fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_data, per_group_data[0]);
+ fprintf(stderr, "%s:\t%d groups, tension handler array %p, inverse tension handler array %p, handler data array %p, (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_inverse_handlers, per_group_data, per_group_data[0]);
#endif
for (i=0; i<num_active_groups; i++) {
tdata = (tension_handler_data_t *) per_group_data[i];
last_x = -1.0;
} /* else continue with the current statics */
- F = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_data, last_x, x, &env->stiffness);
+ F = tension_balance(num_active_groups, per_group_handlers, per_group_inverse_handlers, (void **)per_group_data, last_x, x, &env->stiffness);
last_x = x;
<<tension model getopt definitions>>=
typedef struct tension_model_getopt_struct {
one_dim_fn_t *handler; /* fn implementing the model on a group */
+ one_dim_fn_t *inverse_handler; /* fn implementing inverse on group */
char *name; /* name string identifying the model */
char *description; /* a brief description of the model */
int num_params; /* number of extra params for the model */
create_data_func_t *creator; /* param-string -> model-data-struct fn */
destroy_data_func_t *destructor; /* fn to free the model data structure */
} tension_model_getopt_t;
-@
+@ Set [[inverse_handler]] to [[NULL]] if you wish to invert using
+bisection. Obviously, this will be slower than computing an
+analytical inverse and more prone to rounding errors, but it may be
+easier if a closed-form inverse does not exist.
<<tension model getopt array definition>>=
tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
#ifdef DEBUG
for (i=0; i<n_tension_models; i++) {
- fprintf(stderr, "tension model %i %s:\t model's handler %p\n",
- i, tension_models[i].name, tension_models[i].handler);
+ fprintf(stderr, "tension model %i %s:\t model's handler %p\t inverse handler %p\n",
+ i, tension_models[i].name, tension_models[i].handler, tension_models[i].inverse_handler);
assert(tension_models[i].handler == tension_handlers[i]);
}
#endif
if (num_strings >= 3)
tension_models[tension_model].params = strings[num_strings-1];
#ifdef DEBUG
- fprintf(stderr, "%s:\t%s -> state %s, model %s, params %s (%s), group %g\n", __FUNCTION__, optarg, state_name, tension_models[tension_model].name, tension_models[tension_model].params, tension_group);
+ fprintf(stderr, "%s:\t%s -> state %s, model %s, params %s, group %d\n", __FUNCTION__, optarg, state_name, tension_models[tension_model].name, tension_models[tension_model].params, tension_group);
#endif
INIT_MODEL(state_name, &tension_models[tension_model], tension_models[tension_model].params, model_params);
push(pState_list, create_state(state_name,
tension_models[tension_model].handler,
+ tension_models[tension_model].inverse_handler,
tension_group,
model_params,
tension_models[tension_model].destructor,
typedef struct state_struct {
char *name; /* name string */
one_dim_fn_t *tension_handler; /* tension handler */
+ one_dim_fn_t *inverse_tension_handler; /* inverse tension handler */
int tension_group; /* for grouping similar tension states */
void *tension_params; /* pointer to tension parameters */
destroy_data_func_t *destroy_tension;
<<create/destroy state>>=
state_t *create_state(char *name,
one_dim_fn_t *tension_handler,
+ one_dim_fn_t *inverse_tension_handler,
int tension_group,
void *tension_params,
destroy_data_func_t *destroy_tension,
strcpy(ret->name, name); /* we know ret->name is long enough */
ret->tension_handler = tension_handler;
+ ret->inverse_tension_handler = inverse_tension_handler;
ret->tension_group = tension_group;
ret->tension_params = tension_params;
ret->destroy_tension = destroy_tension;
ret->grouped_states = NULL;
#ifdef DEBUG
- fprintf(stderr, "generated state %s (%p) with tension handler %p, params %p, group %d.\n", ret->name, ret, ret->tension_handler, ret->tension_params, ret->tension_group);
+ fprintf(stderr, "generated state %s (%p) with tension handler %p, inverse handler %p, params %p, group %d.\n", ret->name, ret, ret->tension_handler, ret->inverse_tension_handler, ret->tension_params, ret->tension_group);
#endif
return ret;
}
void get_active_groups(list_t *state_list,
int *pNum_active_groups,
one_dim_fn_t ***pPer_group_handlers,
+ one_dim_fn_t ***pPer_group_inverse_handlers,
void ***pPer_group_data /* pointer to array of pointers to tension_handler_data_t */)
{
void **active_groups=NULL;
- one_dim_fn_t *handler, *old_handler, **per_group_handlers=NULL;
+ one_dim_fn_t *handler, *old_handler, *inverse_handler, *old_inverse_handler;
+ one_dim_fn_t **per_group_handlers=NULL, **per_group_inverse_handlers=NULL;
void **per_group_data=NULL; /* array of pointers to tension_handler_data_t */
int i, j, num_domains, group, old_group, num_active_groups=0;
list_t *active_states_list=NULL;
while (state_list != NULL) {
state = S(state_list);
handler = state->tension_handler;
+ inverse_handler = state->inverse_tension_handler;
group = state->tension_group;
num_domains = state->num_domains;
if (list_index(active_states_list, state) == -1) {
/* allocate the arrays we'll be returning */
per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
assert(per_group_handlers != NULL);
+ per_group_inverse_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
+ assert(per_group_inverse_handlers != NULL);
per_group_data = (void **)malloc(sizeof(void *)*num_active_groups);
assert(per_group_data != NULL);
#ifdef DEBUG
/* populate the arrays we'll be returning */
active_states_list = head(active_states_list);
old_handler = NULL;
+ old_inverse_handler = NULL;
j = -1; /* j is the current active _group_ index */
while (active_states_list != NULL) {
state = (state_t *) pop(&active_states_list);
handler = state->tension_handler;
+ inverse_handler = state->inverse_tension_handler;
group = state->tension_group;
num_domains = state->num_domains;
- if (handler != old_handler || group != old_group) {
+ if (handler != old_handler || inverse_handler != old_inverse_handler || group != old_group) {
j++; /* move to the next group */
tdata = (tension_handler_data_t *) per_group_data[j];
per_group_handlers[j] = handler;
+ per_group_inverse_handlers[j] = inverse_handler;
tdata->group_tension_params = NULL;
tdata->env = NULL;
tdata->persist = NULL; /* we might want to grab any appropriate data from the old persist here... */
}
#ifdef DEBUG
- fprintf(stderr, "%s:\tadding active state %s to group %d (%p:%d) with %d domain(s)\n", __FUNCTION__, state->name, j, handler, group, num_domains);
+ fprintf(stderr, "%s:\tadding active state %s to group %d (%p:%p:%d) with %d domain(s)\n", __FUNCTION__, state->name, j, handler, inverse_handler, group, num_domains);
#endif
for (i=0; i<num_domains; i++) { /* add num_domains copies of the tension params */
push(&tdata->group_tension_params, state->tension_params, 1);
fprintf(stderr, "%s:\tpushed %d copies of %p onto data %p's param list %p\n", __FUNCTION__, num_domains, state->tension_params, tdata, tdata->group_tension_params);
#endif
old_handler = handler;
+ old_inverse_handler = inverse_handler;
old_group = group;
}
/* free old groups */
if (*pPer_group_handlers != NULL)
free(*pPer_group_handlers);
+ if (*pPer_group_inverse_handlers != NULL)
+ free(*pPer_group_inverse_handlers);
if (*pPer_group_data != NULL) {
for (i=0; i<*pNum_active_groups; i++)
free((*pPer_group_data)[i]);
/* set new groups */
#ifdef DEBUG
- fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_data, per_group_data[0]);
+ fprintf(stderr, "%s:\t%d groups, tension handler array %p, inverse tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_inverse_handlers, per_group_data, per_group_data[0]);
#endif
*pNum_active_groups = num_active_groups;
*pPer_group_handlers = per_group_handlers;
+ *pPer_group_inverse_handlers = per_group_inverse_handlers;
*pPer_group_data = per_group_data;
}
@
@
<<worm-like chain tests>>=
-extern double wlc(double x, double T, double p, double L);
+<<worm-like chain function>>
+
START_TEST(test_wlc_at_zero)
{
double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
<<tension balance function declaration>>=
double tension_balance(int num_tension_groups,
one_dim_fn_t **tension_handlers,
+ one_dim_fn_t **inverse_tension_handlers,
void **params, /* array of pointers to tension_handler_data_t */
double last_x,
double x,
<<tension balance function>>=
double tension_balance(int num_tension_groups,
one_dim_fn_t **tension_handlers,
+ one_dim_fn_t **inverse_tension_handlers,
void **params, /* array of pointers to tension_handler_data_t */
double last_x,
double x,
double *stiffness)
{
- static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL};
+ static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL,NULL};
double F, xo_guess, xo, lb, ub;
double min_relx=1e-6, min_rely=1e-6;
int max_steps=200, i;
/* construct new x_xo_data */
x_xo_data.n_groups = num_tension_groups;
x_xo_data.tension_handler = tension_handlers;
+ x_xo_data.inverse_tension_handler = inverse_tension_handlers;
x_xo_data.handler_data = params;
#ifdef DEBUG
- fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, x_xo_data.n_groups, x_xo_data.tension_handler, x_xo_data.handler_data, x_xo_data.handler_data[0]);
+ fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p), x_of_xo_data %p\n", __FUNCTION__, x_xo_data.n_groups, x_xo_data.tension_handler, x_xo_data.handler_data, x_xo_data.handler_data[0], &x_xo_data);
#endif
x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
assert(x_xo_data.xi != NULL);
typedef struct x_of_xo_data_struct {
int n_groups;
one_dim_fn_t **tension_handler; /* array of fn pointers */
+ one_dim_fn_t **inverse_tension_handler; /* inverse fn pointers */
void **handler_data; /* array of pointers to tension_handler_data_t */
double *xi; /* array of group extensions */
} x_of_xo_data_t;
assert(data->n_groups > 0);
data->xi[0] = xo;
#ifdef DEBUG
- fprintf(stderr, "%s:\tget initial tension at %g with %d groups, handler %p, data %p\n", __FUNCTION__, xo, data->n_groups, data->tension_handler[0], data->handler_data[0]);
+ fprintf(stderr, "%s:\tget initial tension at %g with %d groups, handler %p, data %p (x_xo_data %p)\n", __FUNCTION__, xo, data->n_groups, data->tension_handler[0], data->handler_data[0], data);
fflush(stderr);
#endif
F = (*data->tension_handler[0])(xo, data->handler_data[0]);
x += xo;
if (data->xi) data->xi[0] = xo;
for (i=1; i < data->n_groups; i++) {
- if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
- else xi_guess = data->xi[i];
+ if (data->inverse_tension_handler[i] != NULL) {
+ xi = (*data->inverse_tension_handler[i])(F, data->handler_data[i]);
+ } else { /* invert numerically */
+ if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
+ else xi_guess = data->xi[i];
#ifdef DEBUG
- fprintf(stderr, "%s:\tsearching for proper x[%d] for F[%d](x[%d]) = %g N (handler %p, data %p)\n", __FUNCTION__, i, i, i, F, data->tension_handler[i], data->handler_data[i]);
+ fprintf(stderr, "%s:\tsearching for proper x[%d] for F[%d](x[%d]) = %g N (handler %p, data %p)\n", __FUNCTION__, i, i, i, F, data->tension_handler[i], data->handler_data[i]);
#endif
- oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
- xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
+ oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
+ xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
+ }
+ x += xi;
+ data->xi[i] = xi;
#ifdef DEBUG
- fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
+ fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
#endif
- data->xi[i] = xi;
- x += xi;
- if (data->xi) data->xi[i] = xi;
}
#ifdef DEBUG
fprintf(stderr, "%s:\tfound x(x0=%g) = %g\n", __FUNCTION__, xo, x);
{
double yx, step, x=xguess;
int i=0;
+#ifdef DEBUG
+ fprintf(stderr, "%s:\tbracketing %g (params %p)\n", __FUNCTION__, y, params);
+#endif
yx = (*f)(x, params);
#ifdef DEBUG
fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
{
double k=5, x=3, last_x=2.0, F;
one_dim_fn_t *handlers[] = {&hooke};
+ one_dim_fn_t *inverse_handlers[] = {NULL};
void *data[] = {&k};
- F = tension_balance(1, handlers, data, last_x, x, NULL);
+ F = tension_balance(1, handlers, inverse_handlers, data, last_x, x, NULL);
fail_unless(F = k*x, NULL);
}
END_TEST
{
double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
one_dim_fn_t *handlers[] = {&hooke, &hooke};
+ one_dim_fn_t *inverse_handlers[] = {NULL, NULL};
void *data[] = {&k1, &k2};
- F = tension_balance(2, handlers, data, last_x, x, NULL);
+ F = tension_balance(2, handlers, inverse_handlers, data, last_x, x, NULL);
x1e = x*k2/(k1+k2);
Fe = k1*x1e;
x2e = x1e*k1/k2;
For unstretchable domains.
<<null tension model getopt>>=
-{NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
+{NULL, NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
@
\subsection{Constant}
\label{sec.const_tension}
+An infinitely stretchable domain providing a constant tension.
+Obviously this cannot be inverted, so you can't put this domain in
+series with any other domains. We include it mostly for testing
+purposes.
+
<<constant tension functions>>=
<<constant tension function>>
<<constant tension parameter create/destroy functions>>
@
-An infinitely stretchable domain providing a constant tension.
-
<<constant tension function declaration>>=
extern double const_tension_handler(double x, void *pdata);
@
@
<<constant tension model getopt>>=
-{&const_tension_handler, "const", "an infinitely stretchable domain with constant tension", 1, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t}
+{&const_tension_handler, NULL, "const", "an infinitely stretchable domain with constant tension", 1, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t}
@
\subsection{Hooke}
\label{sec.hooke}
<<hooke functions>>=
-<<hooke function>>
+<<internal hooke functions>>
+<<hooke handler>>
+<<inverse hooke handler>>
<<hooke parameter create/destroy functions>>
@
<<hooke tension function declaration>>=
extern double hooke_handler(double x, void *pdata);
+extern double inverse_hooke_handler(double F, void *pdata);
@
-<<hooke function>>=
-double hooke_handler(double x, void *pdata)
-{
- tension_handler_data_t *data = (tension_handler_data_t *)pdata;
+First we define a function that computes the effective spring constant
+(stored in a single [[hooke_param_t]]) for the entire group.
+<<internal hooke functions>>=
+static void hooke_param_grouper(tension_handler_data_t *data,
+ hooke_param_t *param) {
list_t *list = data->group_tension_params;
double k=0.0;
- assert (x >= 0.0);
list = head(list);
assert(list != NULL); /* empty active group?! */
while (list != NULL) {
fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N (from data %p)\n",
__FUNCTION__, k, x, k*x, data);
#endif
- return k*x;
+ param->k = k;
+}
+
+@
+
+Everyone knows Hooke's law: $F=kx$.
+<<hooke handler>>=
+double hooke_handler(double x, void *pdata)
+{
+ hooke_param_t param = {0};
+
+ assert (x >= 0.0);
+ hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
+ return param.k*x;
+}
+
+@
+
+Inverting Hooke's law gives $x=F/k$.
+<<inverse hooke handler>>=
+double inverse_hooke_handler(double F, void *pdata)
+{
+ hooke_param_t param = {0};
+
+ assert (F >= 0.0);
+ hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
+ return F/param.k;
}
+
@
+Not much to keep track of for a single effective spring, just the
+spring constant $k$.
<<hooke structure>>=
typedef struct hooke_param_struct {
double k; /* spring constant in N/m */
} hooke_param_t;
+
@
+We wrap up our Hooke implementation with some book-keeping functions.
<<hooke tension parameter create/destroy function declarations>>=
extern void *string_create_hooke_param_t(char **param_strings);
extern void destroy_hooke_param_t(void *p);
+
@
<<hooke parameter create/destroy functions>>=
@
<<hooke tension model getopt>>=
-{hooke_handler, "hooke", "a Hookean spring (F=kx)", num_hooke_params, hooke_param_descriptions, (char *)hooke_param_string, &string_create_hooke_param_t, &destroy_hooke_param_t}
+{hooke_handler, inverse_hooke_handler, "hooke", "a Hookean spring (F=kx)", num_hooke_params, hooke_param_descriptions, (char *)hooke_param_string, &string_create_hooke_param_t, &destroy_hooke_param_t}
@
\subsection{Worm-like chain}
We can model several unfolded domains with a single worm-like chain.
<<worm-like chain functions>>=
-<<worm-like chain function>>
+<<internal worm-like chain functions>>
<<worm-like chain handler>>
+<<inverse worm-like chain handler>>
<<worm-like chain create/destroy functions>>
@
<<worm-like chain tension model declarations>>=
+
<<worm-like chain handler declaration>>
+<<inverse worm-like chain handler declaration>>
<<worm-like chain create/destroy function declarations>>
<<worm-like chain tension model global declarations>>
@
+<<internal worm-like chain functions>>=
+<<worm-like chain function>>
+<<inverse worm-like chain function>>
+<<worm-like chain parameter grouper>>
+@
+
The combination of all unfolded domains is modeled as a worm like chain,
with the force $F_{\text{WLC}}$ approximately given by
$$
$p$ is the persistence length, and
$L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
<<worm-like chain function>>=
-double wlc(double x, double T, double p, double L)
+static double wlc(double x, double T, double p, double L)
{/* N m K m m */
double xL=0.0; /* uses global k_B */
assert(x >= 0);
// k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
}
+
@
-This model assumes that each unfolded domain has the same persistence length.
-<<worm-like chain handler declaration>>=
-extern double wlc_handler(double x, void *pdata);
+Inverting the approximate WLC is fairly straightforward, if a bit tedious.
+Letting $x_L \equiv x/L$ and $F_T \equiv Fp/k_B T$ we have
+\begin{align}
+ F_T &\approx \frac{1}{4} \p({\frac{1}{(1-x_L)^2} - 1}) + x_L \\
+ F_T (1-x_L)^2 &\approx \frac{1}{4} - \frac{1}{4} (1-x_L)^2 + x_L(1-x_L)^2 \\
+ 0 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) (1-x_L)^2 + x_L(1-x_L)^2 \\
+ &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) \p({1-2x_L+x_L^2})
+ + x_L - 2x_L^2 + x_L^3 \\
+ &\approx \p[{\frac{1}{4} - \p({F_T + \frac{1}{4}})}]
+ + x_L \p[{2F_T + \frac{1}{2} + 1}]
+ + x_L^2 \p[{-\p({F_T + \frac{1}{4}}) - 2}]
+ + x_L^3 \\
+ &\approx -F_T
+ + x_L \p({2F_T + \frac{3}{2}})
+ - x_L^2 \p({F_T + \frac{9}{4}})
+ + x_L^3 \\
+ &\approx ax_L^3 + bx_L^2 + cx_L + d
+\end{align}
+where $a \equiv 1$, $b \equiv -F_T - 9/4$, $c \equiv 2F_T + 3/2$, and
+$d \equiv -F_T$.
+% From \citet{wikipedia_cubic_function} the discriminant
+%\begin{align}
+% \Delta &= -4b^3d + b^2c^2 - 4ac^3 + 18abcd - 27a^2d^2 \\
+% &= -4F_T\p({F_T + \frac{9}{4}})^3
+% + \p({F_T + \frac{9}{4}})^2 \cdot \p({2F_T + \frac{3}{2}})^2
+% - 4 \p({2F_T + \frac{3}{2}})^3
+% + 18F_T \p({F_T + \frac{9}{4}}) \cdot \p({2F_T + \frac{3}{2}})
+% - 27F_T^2 \\
+% &= -4F_T \p({F_T^3 + \frac{27}{4}F_T^2 +\frac{243}{16}F_T +\frac{729}{64}})
+%% a^3 + 3a^2b + 3ab^2 + b^3
+% + \p({F_T^2 + \frac{9}{2}F_T + \frac{81}{16}})
+% \cdot \p({4F_T^2 + 6F_T + \frac{9}{4}}) \\
+% &\qquad - 4 \p({8F_T^3 + 18F_T^2 + \frac{27}{2}F_T + \frac{27}{8}})
+%% a^3 + 3a^2b + 3ab^2 + b^3
+% + 18F_T \p({2F_T^2 + 6F_T + \frac{27}{8}})
+% - 27F_T^2 \\
+% &= -4F_T^4 - 27F_T^3 - \frac{243}{4}F_T^2 - \frac{729}{16}F_T
+% + 4F_T^4 + 24F_T^3 +\frac{99}{2}F_T^2 +\frac{81}{2}F_T +\frac{729}{64}\\
+% &\qquad - 32F_T^3 - 72F_T^2 - 54F_T - \frac{27}{2}
+% + 36F_T^3 + 108F_T^2 + \frac{243}{4}F_T
+% - 27F_T^2 \\
+% &= F_T^4 \p({-4 + 4}) + F_T^3 \p({-27 + 24 - 32 + 36})
+% + F_T^2 \p({-\frac{243}{4} + \frac{99}{2} - 72 + 108 -27}) \\
+% &\qquad + F_T \p({-\frac{729}{16} + \frac{81}{2} - 54 + \frac{243}{4}})
+% \p({\frac{729}{64} - \frac{27}{2}}) \\
+% &= F_T^3 - \frac{9}{4}F_T^2 + \frac{27}{16}F_T - \frac{135}{64}
+%\end{align}
+We can plug these coefficients into the GSL cubic solver to invert the
+WLC\citep{galassi05}. The applicable root is always the one which
+comes before the singularity, which will be the smallest real root.
+<<inverse worm-like chain function>>=
+static double inverse_wlc(double F, double T, double p, double L)
+{/* m N K m m */
+ double FT = F*p/(k_B*T); /* uses global k_B */
+ double xL0, xL1, xL2;
+ int num_roots;
+ assert(F >= 0);
+ assert(T > 0); assert(p > 0); assert(L > 0);
+ num_roots = gsl_poly_solve_cubic(-(FT+2.25),2*FT+1.5,-FT, &xL0, &xL1, &xL2);
+ assert(num_roots > 0); assert(xL0 >= -DOUBLE_PRECISION); assert(xL0 < 1);
+ if (xL0 < 0) xL0 = 0.0;
+ return xL0*L;
+}
+
@
-<<worm-like chain handler>>=
-double wlc_handler(double x, void *pdata)
-{
- tension_handler_data_t *data = (tension_handler_data_t *)pdata;
+First we define a function that computes the effective WLC parameters
+(stored in a single [[wlc_param_t]]) for the entire group.
+<<worm-like chain parameter grouper>>=
+static void wlc_param_grouper(tension_handler_data_t *data,
+ wlc_param_t *param) {
list_t *list = data->group_tension_params;
- double p, L=0.0;
+ double p=0.0, L=0.0;
list = head(list);
assert(list != NULL); /* empty active group?! */
fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
__FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
#endif
- return wlc(x, data->env->T, p, L);
+ param->p = p;
+ param->L = L;
}
+
+@
+
+<<worm-like chain handler declaration>>=
+extern double wlc_handler(double x, void *pdata);
+@
+
+This model requires that each unfolded domain in the group have the
+same persistence length.
+<<worm-like chain handler>>=
+double wlc_handler(double x, void *pdata)
+{
+ tension_handler_data_t *data = (tension_handler_data_t *)pdata;
+ wlc_param_t param = {0};
+
+ wlc_param_grouper(data, ¶m);
+ return wlc(x, data->env->T, param.p, param.L);
+}
+
+@
+
+<<inverse worm-like chain handler declaration>>=
+extern double inverse_wlc_handler(double F, void *pdata);
+@
+
+<<inverse worm-like chain handler>>=
+double inverse_wlc_handler(double F, void *pdata)
+{
+ tension_handler_data_t *data = (tension_handler_data_t *)pdata;
+ wlc_param_t param = {0};
+
+ wlc_param_grouper(data, ¶m);
+ return inverse_wlc(F, data->env->T, param.p, param.L);
+}
+
@
<<worm-like chain structure>>=
if (p)
free(p);
}
+
@
We define [[wlc_param_string]] as a global array because hard-coding the values into the tension model getopt structure gives a read-only version.
@
<<worm-like chain tension model getopt>>=
-{wlc_handler, "wlc", "a worm-like chain (F=kbT/p{.25[(1-x/L)^-2 -1] + x/L})", num_wlc_params, wlc_param_descriptions, (char *)wlc_param_string, &string_create_wlc_param_t, &destroy_wlc_param_t}
+{wlc_handler, inverse_wlc_handler, "wlc", "a worm-like chain (F=kbT/p{.25[(1-x/L)^-2 -1] + x/L})", num_wlc_params, wlc_param_descriptions, (char *)wlc_param_string, &string_create_wlc_param_t, &destroy_wlc_param_t}
@
Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
@
<<freely-jointed chain tension model getopt>>=
-{fjc_handler, "fjc", "a freely-jointed chain [F=kbT/l L^-1(x/L)]", num_fjc_params, fjc_param_descriptions, (char *)fjc_param_string, &string_create_fjc_param_t, &destroy_fjc_param_t}
+{fjc_handler, NULL, "fjc", "a freely-jointed chain [F=kbT/l L^-1(x/L)]", num_fjc_params, fjc_param_descriptions, (char *)fjc_param_string, &string_create_fjc_param_t, &destroy_fjc_param_t}
@
\subsection{Tension model implementation}
semi-clean_latex :
rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
$(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
- $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.tex \
- $(BUILD_DIR)/sawsim.bib
+ $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.toc \
+ $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib
clean_latex : semi-clean_latex
rm -f $(DOC_DIR)/sawsim.pdf
@
projects / sawsim.git / commitdiff