\item Hookean spring (Appendix \ref{sec.hooke}),
\item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
\item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
+ \item Piston (Appendix \ref{sec.piston_tension}),
\end{packed_item}
The tension model and parameters used for a particular domain depend
should get an appropriate error message.
<<find tension definitions>>=
-#define NUM_TENSION_MODELS 5
+#define NUM_TENSION_MODELS 6
@
<<constant tension model getopt>>,
<<hooke tension model getopt>>,
<<worm-like chain tension model getopt>>,
- <<freely-jointed chain tension model getopt>>
+ <<freely-jointed chain tension model getopt>>,
+ <<piston tension model getopt>>
};
@
extern char *optarg;
extern int optind, optopt, opterr;
- assert (argc > 0);
+ assert(argc > 0);
#ifdef DEBUG
for (i=0; i<n_tension_models; i++) {
<<x of x0>>
<<group stiffness function>>
<<chain stiffness function>>
+<<full chain stiffness function>>
<<tension balance function>>
@
F = (*tension_handlers[0])(xo, params[0]);
- if (stiffness != NULL)
+ if (stiffness != NULL) {
*stiffness = chain_stiffness(&x_xo_data, min_relx);
-
+ if (*stiffness == 0.0) { /* re-calculate based on full chain */
+ *stiffness = full_chain_stiffness(num_tension_groups, tension_handlers,
+ inverse_tension_handlers, params,
+ last_x, x, min_relx, F);
+ }
+ }
return F;
}
}
@
+Determine the stiffness of a single tension group by taking a small
+step [[dx]] back from the position [[x]] for which we want the
+stiffness. The touchy part is determining a good [[dx]] and ensuring
+the whole interval is on [[x>=0]].
<<group stiffness function>>=
double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double x, double relx)
{
F = (*tension_handler)(x, handler_data);
Fl = (*tension_handler)(xl, handler_data);
stiffness = (F-Fl)/dx;
- assert(stiffness > 0);
+ assert(stiffness >= 0);
return stiffness;
}
@
+Attempt to calculate the chain stiffness by adding group stiffnesses
+as springs in series. In the case where a group has undefined
+stiffness (e.g. piston tension, Section \ref{sec.piston_tension}),
+this algorithm breaks down. In that case, [[tension_balance()]]
+(\ref{sec.tension_balance}) should use [[full_chain_stiffness()]]
+which uses the full chain's $F(x)$ rather than that of the individual
+domains'. We attempt the series approach first because, lacking the
+numerical inversion inside [[tension_balance()]], it is both faster
+and more accurate than the full-chain derivative.
<<chain stiffness function>>=
double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx)
{
- double stiffness;
+ double stiffness, gstiff;
int i;
/* add group stiffnesses in series */
for (i=0; i < x_xo_data->n_groups; i++) {
fprintf(stderr, "%s:\tgetting stiffness of active state %d of %d for x[%d]=%g, relx=%g\n", __FUNCTION__, i, x_xo_data->n_groups, i, x_xo_data->xi[i], relx);
fflush(stderr);
#endif
- stiffness += 1.0/group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx);
+ gstiff = group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx);
+ if (gstiff == 0.0);
+ return 0.0;
+ stiffness += 1.0/gstiff;
}
stiffness = 1.0 / stiffness;
return stiffness;
@
+Determine the chain stiffness using only [[tension_balance()]]. This
+works around difficulties with tension models that have undefined
+stiffness (see discussion for [[chain_stiffness()]] above).
+<<full chain stiffness function>>=
+double full_chain_stiffness(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 relx,
+ double F /* already evaluated F(x) */)
+{
+ double dx, xl, Fl, stiffness;
+
+ assert(x >= 0);
+ assert(relx > 0 && relx < 1);
+
+ /* Other option for dx that we ignore because of our poor tension_balance()
+ * resolution (i.e. bad slopes if you zoom in too close):
+ * if (last_x != -1.0 && last_x != x) // last_x limited
+ * dx fabs(x-last_x);
+ * else
+ * dx = HUGE_VAL;
+ * if (1==1) { // constant max-value limited
+ * dx_p = 1e-12;
+ * dx = MIN(dx, dx_p);
+ * }
+ * if (x != 0 && relx != 0) { // relx limited
+ * dx_p = x*relx;
+ * dx = MIN(dx, dx_p);
+ * }
+ * TODO, determine which of (min_relx, min_rely, max_steps) is actually
+ * limiting tension_balance.
+ */
+ dx = 1e-12; /* HACK, how to get length scale? */
+
+ xl = x-dx;
+ if (xl >= 0) {
+ Fl = tension_balance(num_tension_groups, tension_handlers,
+ inverse_tension_handlers, params,
+ last_x, xl, NULL);
+ } else {
+ xl = x;
+ Fl = F;
+ x += dx;
+ F = tension_balance(num_tension_groups, tension_handlers,
+ inverse_tension_handlers, params,
+ last_x, x, NULL);
+ }
+ stiffness = (F-Fl)/dx;
+ assert(stiffness >= 0);
+ return stiffness;
+}
+
+@
+
Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
number of steps for monotonically increasing $f(x)$. Simple
bisection, so it's robust and converges fairly quickly. You can set
<<hooke tension model declarations>>
<<worm-like chain tension model declarations>>
<<freely-jointed chain tension model declarations>>
+<<piston tension model declarations>>
<<find tension definitions>>
<<tension model global declarations>>
#endif /* TENSION_MODEL_H */
list_t *list = data->group_tension_params;
double F;
- assert (x >= 0.0);
+ assert(x >= 0.0);
list = head(list);
assert(list != NULL); /* empty active group?! */
F = ((const_tension_param_t *)list->d)->F;
@
<<constant tension model getopt>>=
-{&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}
+{&const_tension_handler, NULL, "const", "an infinitely stretchable domain with constant tension", num_const_tension_params, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t}
@
\subsection{Hooke}
{
hooke_param_t param = {0};
- assert (x >= 0.0);
+ assert(x >= 0.0);
hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
return param.k*x;
}
{
hooke_param_t param = {0};
- assert (F >= 0.0);
+ assert(F >= 0.0);
hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
return F/param.k;
}
tension_handler_data_t *data = (tension_handler_data_t *)pdata;
wlc_param_t param = {0};
+ assert(x >= 0.0);
wlc_param_grouper(data, ¶m);
return wlc(x, data->env->T, param.p, param.L);
}
@
<<worm-like chain tension model getopt>>=
-{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}
+{&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.
{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{Piston}
+\label{sec.piston_tension}
+
+An infinitely stretchable domain with no tension for extensions less
+than a particular threshold $L$, and infinite tension for $x>L$. The
+tension at $x=L$ is undefined, but may be any positive value. The
+model is called the ``piston'' model because it resembles a
+frictionless piston sliding in a rigid cylinder of length $L$.
+
+Note that because of it's infinte stiffness at $x=L$, fully extended
+domains with this tension model will be identical to completely rigid
+domains. However, there is a distinction between this model and rigid
+domains for $x<L$, because any reactions that occur at $F=0$
+(e.g. refolding) will have more time to occur.
+
+Because the tension at $x=L$ is undefined, the user must make sure
+domains with this tension model are never the initial active group,
+because it would break [[tension_balance()]] in [[find_tension()]]
+(see Section \ref{sec.tension_balance}).
+
+<<piston tension functions>>=
+<<piston tension parameter grouper>>
+<<piston tension handler>>
+<<inverse piston tension handler>>
+<<piston tension parameter create/destroy functions>>
+@
+
+<<piston tension model declarations>>=
+<<piston tension handler declaration>>
+<<inverse piston tension handler declaration>>
+<<piston tension parameter create/destroy function declarations>>
+<<piston tension model global declarations>>
+
+@
+
+First we define a function that computes the effective piston parameters
+(stored in a single [[piston_tension_param_t]]) for the entire group.
+<<piston tension parameter grouper>>=
+static void piston_tension_param_grouper(tension_handler_data_t *data,
+ piston_tension_param_t *param) {
+ list_t *list = data->group_tension_params;
+ double L=0;
+
+ list = head(list);
+ assert(list != NULL); /* empty active group?! */
+ while (list != NULL) {
+ L += ((piston_tension_param_t *)list->d)->L;
+ list = list->next;
+ }
+ param->L = L;
+}
+
+<<piston tension handler declaration>>=
+extern double piston_tension_handler(double x, void *pdata);
+@
+<<piston tension handler>>=
+double piston_tension_handler(double x, void *pdata)
+{
+ tension_handler_data_t *data = (tension_handler_data_t *)pdata;
+ piston_tension_param_t param = {0};
+ double F;
+
+ piston_tension_param_grouper(data, ¶m);
+ if (x <= param.L) F = 0;
+ else F = HUGE_VAL;
+#ifdef DEBUG
+ fprintf(stderr, "%s: x %g, L %g, tension %g N\n", __FUNCTION__, x, param.L, F);
+#endif
+ return F;
+}
+
+@
+
+<<inverse piston tension handler declaration>>=
+extern double inverse_piston_tension_handler(double x, void *pdata);
+@
+<<inverse piston tension handler>>=
+double inverse_piston_tension_handler(double F, void *pdata)
+{
+ tension_handler_data_t *data = (tension_handler_data_t *)pdata;
+ piston_tension_param_t param = {0};
+
+ assert(F >= 0.0);
+ piston_tension_param_grouper(data, ¶m);
+ if (F == 0.0) return 0;
+ return param.L;
+}
+
+@
+
+<<piston tension structure>>=
+typedef struct piston_tension_param_struct {
+ double L; /* length of piston in m */
+} piston_tension_param_t;
+
+@
+
+<<piston tension parameter create/destroy function declarations>>=
+extern void *string_create_piston_tension_param_t(char **param_strings);
+extern void destroy_piston_tension_param_t(void *p);
+
+@
+
+<<piston tension parameter create/destroy functions>>=
+piston_tension_param_t *create_piston_tension_param_t(double L)
+{
+ piston_tension_param_t *ret
+ = (piston_tension_param_t *) malloc(sizeof(piston_tension_param_t));
+ assert(ret != NULL);
+ ret->L = L;
+ return ret;
+}
+
+void *string_create_piston_tension_param_t(char **param_strings)
+{
+ return create_piston_tension_param_t(
+ safe_strtod(param_strings[0],__FUNCTION__));
+}
+
+void destroy_piston_tension_param_t(void *p)
+{
+ if (p)
+ free(p);
+}
+
+@
+
+<<piston tension model global declarations>>=
+extern int num_piston_tension_params;
+extern char *piston_tension_param_descriptions[];
+extern char piston_tension_param_string[];
+
+@
+
+<<piston tension model globals>>=
+int num_piston_tension_params = 1;
+char *piston_tension_param_descriptions[] = {"length L, m"};
+char piston_tension_param_string[] = "0";
+
+@
+
+<<piston tension model getopt>>=
+{&piston_tension_handler, &inverse_piston_tension_handler, "piston", "a domain that slides without tension for x<L, but has infinite tension for x>L.", num_piston_tension_params, piston_tension_param_descriptions, piston_tension_param_string, &string_create_piston_tension_param_t, &destroy_piston_tension_param_t}
+@
+
\subsection{Tension model implementation}
<<tension-model.c>>=
<<hooke structure>>
<<worm-like chain structure>>
<<freely-jointed chain structure>>
+<<piston tension structure>>
@
<<tension model globals>>=
<<hooke tension model globals>>
<<worm-like chain tension model globals>>
<<freely-jointed chain tension model globals>>
+<<piston tension model globals>>
@
<<tension model functions>>=
<<hooke functions>>
<<worm-like chain functions>>
<<freely-jointed chain functions>>
+<<piston tension functions>>
@
-
\subsection{Utilities}
The tension models can get complicated, and users may want to reassure
extern char *optarg;
extern int optind, optopt, opterr;
- assert (argc > 0);
+ assert(argc > 0);
/* setup defaults */
extern char *optarg;
extern int optind, optopt, opterr;
- assert (argc > 0);
+ assert(argc > 0);
/* setup defaults */
projects / sawsim.git / blobdiff