computation time by about a factor of 2. I was expecting more, but
I'll take what I can get.
+Version 0.8 fixed a minor bug in unfolding probability for
+multi-domain groups. The probability of at least one domain unfolding
+had been calculated as $P_N=NP_1$, rather than $P_N=1-(1-P_1)^N$.
+However, for small $P$ the two are equivalent.
+
+Version 0.9 the -P option to sawsim now sets the target probability
+for the per-state transition $P_N$, not the per-domain transisiton
+$P_1$.
+
<<version definition>>=
-#define VERSION "0.7"
+#define VERSION "0.9"
@
\subsection{License}
\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
-on whether the domain's current state. The transition rates between
-states are also handled generally with function pointers, with
+on the domain's current state. The transition rates between states
+are also handled generally with function pointers, with
implementations for
\begin{packed_item}
\item Null (Appendix \ref{sec.null_k}),
model has a tension handler function, which gives the tension $F$
needed to stretch a given group of domains a distance $x$.
-GROUPS ARE CURRENTLY DISABLED.
-
Groups represent collections of states which the model should treat as
a single entity. To understand the purpose of groups, consider a
system with two unfolded domain states (e.g.\ for two different
balancing the tension and the error associated with the inter-group
interaction linearity. Note that group numbers only matter
\emph{within} model classes, since grouping domains with seperate
-models doesn't make sense.
+models doesn't make sense. Within designated groups, the tension
+parameters for each domain are still checked for compatibility, so if
+you accidentally assign incompatible domains to the same group you
+should get an appropriate error message.
<<find tension definitions>>=
-#define NUM_TENSION_MODELS 5
+#define NUM_TENSION_MODELS 6
@
to 0 if there were no transitions since the previous call to
[[find_tension]] to support some optimizations that assume a small
increase in tension between steps (only true if no transition has
-occured). While were messing with the tension handlers, we also grab
+occured). While we're messing with the tension handlers, we also grab
a measure of the current linker stiffness for the environmental
variable, which is needed by some of the unfolding rate models
(e.g.\ the linker-stiffness corrected Bell model, Appendix
parametrized,
a function to create the parameter structure,
a function to destroy the parameter structure, and
- a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as it's data argument.
+ a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as its data argument.
The tension-specific parameter structures are contained in the domain
-group field. See the Section \ref{sec.model_selection} for a
-discussion on the command line framework. See the worm-like chain
+group field. For optimization, a group may also define an inverse
+tension function as an optimization. See the Section
+\ref{sec.model_selection} for a discussion on the command line
+framework and inverse function discussion. See the worm-like chain
implementation for an example (Section \ref{sec.wlc}).
The major limitation of our approach is that all of our tension models
So the probability of a given domain transitioning in a short time
$dt$ is given by
$$
- P = k \cdot dt.
+ P_1 = k \cdot dt.
$$
+For a group of $N$ identical domains, and therefore identical
+unfolding probabilities $P_1$, the probability of $n$ domains
+unfolding in a given timestep is given by the binomial distribution
+$$
+ P(n) = \frac{N!}{n!(N-n)!}P_1^n(1-P_1)^(N-n).
+$$
+The probability that \emph{none} of these domains unfold is then
+$$
+ P(0) = (1-P_1)^N,
+$$
+so the probability that \emph{at least one} domain unfolds is
+$$
+ P_N = 1-(1-P_1)^N.
+$$
+Note that for small $P$,
+$$
+ p(n>0) = 1-(1-NP_1+\mathcal{O}(P_1^2)) = NP_1 - \mathcal{O}(P^2)
+ \approx NP_1.
+$$
+This conversion from $P_1$ to $P_N$ is handled by the [[pN]] macro
+<<PN definition>>=
+/* find multi-domain transition rate from the single domain rate */
+#define PN(P,N) (1.0-pow(1.0-(P), (N)))
+
+@
+
+We can also define a macro to find the approprate $dt$ to achieve a
+target $P_N$ with a given $k$ and $N$.
+\begin{align}
+ P_N &= 1-(1-P_1)^N \\
+ 1-P_1 &= (1-P_N)^{1/N} \\
+ P_1 &= 1-(1-P_N)^{1/N}
+\end{align}
+<<P1 definition>>=
+#define P1(PN,N) (1.0-pow(1.0-(PN), 1.0/(N)))
+@
+
+We take some time to discuss the meaning of $p(n>1)$
+(i.e. multi-unfolding timesteps) in Section \ref{sec.timescales}.
+
<<does domain transition>>=
-int domain_transitions(double F, double dt, environment_t *env, int num_domains, transition_t *transition)
-{ /* returns 1 or 0, F in N, dt in s, pointer to env. data, number of 'reactant' domains, pointer to transition structure */
+int domain_transitions(double F, double dt, environment_t *env, transition_t *transition)
+{ /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to transition structure */
double k;
k = accel_k(transition->k, F, env, transition->k_params);
//(*transition->k)(F, env, domain->k_params);
//printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
- return happens(k*dt*num_domains); /* N dice rolls for prob. k*dt event */
+ return happens(PN(k*dt, N_INIT(transition)));
}
@ [[happens]] is a random decision making function defined in Appendix \ref{sec.utils}.
{ /* returns 1 if there was a transition and 0 otherwise */
while (transition_list != NULL) {
if (T(transition_list)->initial_state->num_domains > 0
- && domain_transitions(tension, dt, env, T(transition_list)->initial_state->num_domains, T(transition_list))) {
+ && domain_transitions(tension, dt, env, T(transition_list))) {
if (flags & FLAG_OUTPUT_TRANSITION_FORCES)
fprintf(stdout, "%g\t%s\t%s\n", tension, T(transition_list)->initial_state->name, T(transition_list)->final_state->name);
T(transition_list)->initial_state->num_domains--;
Besides assuming our timestep is much greater than equilibration
times, we also force it to be short enough so that only one domain may
unfold in a given timestep. For an unfolding timescale of $dt_u =
-1/k$ we adapt our timesteps to keep $dt \ll dt_u$, so the probability
-of two domains unfolding in the same timestep is negligible. This
-approach breaks down as the adaptive timestep scheme approaches $dt
-\sim dt_u$, but $dt_u \sim 1\U{$\mu$s}$ for Ig27-like proteins
-\citep{klimov00}, so this shouldn't be a problem. To reassure
+1/k$, we adapt our timesteps to keep $dt \ll dt_u$, so the probability
+of two domains unfolding in the same timestep is small (see
+[[determine_dt]] in Section \ref{sec.adaptive_dt} and [[P]] in
+[[main()]] in Section \ref{sec.main} set by [[-P]] in
+[[get_options()]] in Section \ref{sec.get_options}). This approach
+breaks down as the adaptive timestep scheme approaches the transition
+timestep $dt_t$, but $dt_t \sim 1\U{$\mu$s}$ for Ig27-like proteins
+\citep{klimov00}, so this shouldn't be a problem. To reassure
yourself, you can ask the simulator to print the smallest timestep
that was used with TODO.
+Even if the likelyhood of two domains unfolding in the same timestep
+is small, if you run long enough simulations it will eventually occur.
+In this case, we assume that $dt_t \ll dt$, so even if two domains
+would unfold if we stuck to our unfolding rate $k$ for an entire
+timestep $dt$, in ``reality'' only one of those domains unfolds.
+Because we do not know when in the timestep the unfolding took place,
+we move on to the next timestep without checking for further
+unfoldings. This ``unchecked timestep portion'' should not contribute
+significant errors to the calculation, because if $dt$ was small
+enough that unfolding was unlikely at the pre-unfolding force, it
+should be even more unlikely at the post-unfolding force, especially
+over only a fraction of the timestep. The only dangerous cases for
+the current implementation are those where unfolding intermediates are
+much less stable than their precursor states, in which case an
+unfolding event during the remaining timestep may actually be likely.
+A simple workaround for such cases is to pick the value for [[P]] to
+be small enough that the $dt \ll dt_u$ assumption survives even under
+a transition populating the unstable intermediate.
+
\subsection{Adaptive timesteps}
\label{sec.adaptive_dt}
<<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>>
};
@
@
\subsection{Get options}
+\label{sec.get_options}
<<get options>>=
+<<safe strto*>>
void get_options(int argc, char **argv,
double *pP, double *pT_max, double *pDt_max, double *pV,
double *pX_step, state_t **pStop_state, environment_t *env,
extern char *optarg;
extern int optind, optopt, opterr;
- assert (argc > 0);
+ assert(argc > 0);
#ifdef DEBUG
for (i=0; i<n_tension_models; i++) {
}
}
break;
- case 't': *pT_max = atof(optarg); break;
- case 'd': *pDt_max = atof(optarg); break;
- case 'P': *pP = atof(optarg); break;
- case 'v': *pV = atof(optarg); break;
- case 'x': *pX_step = atof(optarg); break;
- case 'T': env->T = atof(optarg); break;
- case 'C': env->T = atof(optarg)+273.15; break;
+ case 't': *pT_max = safe_strtod(optarg, "-t"); break;
+ case 'd': *pDt_max = safe_strtod(optarg, "-d"); break;
+ case 'P': *pP = safe_strtod(optarg, "-P"); break;
+ case 'v': *pV = safe_strtod(optarg, "-v"); break;
+ case 'x': *pX_step = safe_strtod(optarg, "-x"); break;
+ case 'T': env->T = safe_strtod(optarg, "-T"); break;
+ case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
case 's':
parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
assert(num_strings >= 2 && num_strings <= 4);
}
tension_model = index_tension_model(n_tension_models, tension_models, strings[1]);
if (num_strings == 4)
- tension_group = atoi(strings[2]);
+ tension_group = safe_strtoi(strings[2], optarg);
else
tension_group = 0;
if (num_strings >= 3)
free_parsed_list(num_strings, strings);
break;
case 'N':
- n = atoi(optarg);
+ n = safe_strtoi(optarg, "-N");
#ifdef DEBUG
fprintf(stderr, "%s:\tadding %d domains to %s\n", __FUNCTION__, n, S(*pState_list)->name);
#endif
} while (0);
@
+First we define some safe string parsing functions. Currently these
+abort on any irregularity, but in the future they could print error
+messages, etc. [[static]] because the functions are currently
+defined in each [[*.c]] file that uses them, but perhaps they should
+be moved to [[utils.h/c]] or some such instead\ldots
+
+<<safe strto*>>=
+static int safe_strtoi(const char *s, const char *description) {
+ char *endp = NULL;
+ long int ret;
+ assert(s != NULL);
+ ret = strtol(s, &endp, 10);
+ if (endp[0] != '\0') { //strlen(endp) > 0) {
+ fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
+ endp, s, description);
+ assert(1==0); //strlen(endp) == 0);
+ } if (endp == s) {
+ fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
+ s, description);
+ assert(endp != s);
+ } else if (errno == ERANGE) {
+ fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
+ assert(errno != ERANGE);
+ }
+ return (int) ret;
+}
+
+static double safe_strtod(const char *s, const char *description) {
+ char *endp = NULL;
+ double ret;
+ assert(s != NULL);
+ ret = strtod(s, &endp);
+ if (endp[0] != '\0') { //strlen(endp) > 0) {
+ fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
+ endp, s, description);
+ assert(1==0); //strlen(endp) == 0);
+ } if (endp == s) {
+ fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
+ s, description);
+ assert(endp != s);
+ } else if (errno == ERANGE) {
+ fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
+ assert(errno != ERANGE);
+ }
+ return ret;
+}
+
+@
+
\phantomsection
\appendix
\addcontentsline{toc}{section}{Appendicies}
We include [[math.h]], so don't forget to link to the libm with `-lm'.
<<includes>>=
#include <assert.h> /* assert() */
-#include <stdlib.h> /* malloc(), free(), rand() */
+#include <stdlib.h> /* malloc(), free(), rand(), strto*() */
#include <stdio.h> /* fprintf(), stdout */
#include <string.h> /* strlen, strtok() */
#include <math.h> /* exp(), M_PI, sqrt() */
#include <sys/types.h> /* pid_t (returned by getpid()) */
#include <unistd.h> /* getpid() (for seeding rand()), getopt() */
+#include <errno.h> /* errno, ERANGE (for safe_strto*()) */
#include "global.h"
#include "list.h"
#include "tension_balance.h"
$F(x)$ increases with $x$, possibly exploding (e.g.\ in the worm-like
chain model), so basing the timestep on the the unfolding probability
at the current tension is dangerous, and we need to search for a $dt$
-for which $P(F(x+v*dt)) < P_\text{target}$. There are two cases to
+for which $P_N(F(x+v*dt)) < P_\text{target}$. There are two cases to
consider. In the most common, no domains have unfolded since the last
step, and we expect the next step to be slightly shorter than the
previous one. In the less common, domains did unfold in the last
* a starting separation x in m, a target_prob between 0 and 1,
* max_dt in s, stretching velocity v in m/s.
*/
- double F, k, dtCur, dtU, dtUCur, dtL, dt;
+ double F, k, P, dtCur, dtU, dtUCur, dtL, dt;
/* get upper bound using the current position */
F = find_tension(state_list, env, x, 0); /* BUG. repeated calculation */
//printf("Start with x = %g (F = %g)\n", x, F);
k = accel_k(transition->k, F, env, transition->k_params);
//printf("x %g\tF %g\tk %g\n", x, F, k);
- dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */
+ dtU = P1(target_prob, N_INIT(transition)) / k; /* upper bound on dt */
if (dtU > max_dt) {
//printf("overshot max_dt\n");
dtU = max_dt;
//printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
}
assert(dtU > 1e-14); /* timestep to valid k too small */
- dtUCur = target_prob / k; /* safe timestep back from x+dtU */
+ /* safe timestep back from x+dtU */
+ dtUCur = P1(target_prob, N_INIT(transition)) / k;
if (dtUCur >= dt)
return dt; /* dtU is safe. */
F = find_tension(state_list, env, x+v*dt, 0);
//printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
k = accel_k(transition->k, F, env, transition->k_params);
- dtCur = target_prob / k;
+ dtCur = P1(target_prob, N_INIT(transition)) / k;
//printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
if (dtCur > dt) /* safe timestep back from x+dt covers dt */
dtL = dt;
/* get the transition data for the current list node */
#define T(list) ((transition_t *)(list)->d)
+
+/* get the number of initial state population for a transition state */
+#define N_INIT(transition) ((transition)->initial_state->num_domains)
+
+<<PN definition>>
+<<P1 definition>>
@ [[k]] is a pointer to the function determining the transition rate
for a given tension. [[k_params]] and [[destroy_k]] have similar
roles with regards to [[k]] as [[tension_params]] and
@
<<parse check includes>>=
-#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
-#include <stdio.h> /* printf() */
-#include <assert.h> /* assert() */
-#include <string.h> /* strlen() */
+#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
+#include <stdio.h> /* printf() */
+#include <assert.h> /* assert() */
+#include <string.h> /* strlen() */
<<check includes>>
#include "parse.h"
@
<<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
@
<<string eval check includes>>=
-#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
-#include <stdio.h> /* printf() */
-#include <assert.h> /* assert() */
-#include <string.h> /* strlen() */
-#include <signal.h> /* SIGKILL */
+#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
+#include <stdio.h> /* printf() */
+#include <assert.h> /* assert() */
+#include <string.h> /* strlen() */
+#include <signal.h> /* SIGKILL */
<<check includes>>
#include "string_eval.h"
@
<<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;
void *string_create_const_tension_param_t(char **param_strings)
{
- return create_const_tension_param_t(atof(param_strings[0]));
+ return create_const_tension_param_t(
+ safe_strtod(param_strings[0],__FUNCTION__));
}
void destroy_const_tension_param_t(void *p)
@
<<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;
}
void *string_create_hooke_param_t(char **param_strings)
{
- return create_hooke_param_t(atof(param_strings[0]));
+ return create_hooke_param_t(safe_strtod(param_strings[0], __FUNCTION__));
}
void destroy_hooke_param_t(void *p)
int num_roots;
assert(F >= 0);
assert(T > 0); assert(p > 0); assert(L > 0);
+ if (F == HUGE_VAL)
+ return L;
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;
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);
}
void *string_create_wlc_param_t(char **param_strings)
{
- return create_wlc_param_t(atof(param_strings[0]), atof(param_strings[1]));
+ return create_wlc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
+ safe_strtod(param_strings[1], __FUNCTION__));
}
void destroy_wlc_param_t(void *p /* wlc_param_t * */)
@
<<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.
void *string_create_fjc_param_t(char **param_strings)
{
- return create_fjc_param_t(atof(param_strings[0]), atof(param_strings[1]));
+ return create_fjc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
+ safe_strtod(param_strings[1], __FUNCTION__));
}
void destroy_fjc_param_t(void *p)
{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>>=
<<tension model includes>>=
#include <assert.h> /* assert() */
-#include <stdlib.h> /* NULL */
+#include <stdlib.h> /* NULL, strto*() */
#include <math.h> /* HUGE_VAL, sqrt(), exp() */
#include <stdio.h> /* fprintf(), stdout */
+#include <errno.h> /* errno, ERANGE */
#include "global.h"
#include "list.h"
#include "tension_balance.h" /* oneD_*() */
<<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>>=
+<<safe strto*>>
<<constant tension 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
<<tension model utility includes>>=
#include <stdlib.h>
#include <stdio.h>
-#include <string.h> /* strlen() */
-#include <assert.h> /* assert() */
+#include <string.h> /* strlen() */
+#include <assert.h> /* assert() */
+#include <errno.h> /* errno, ERANGE */
#include "global.h"
#include "parse.h"
#include "list.h"
@
<<tension model utility getopt functions>>=
+<<safe strto*>>
<<index tension model>>
<<tension model utility help>>
<<tension model utility get options>>
extern char *optarg;
extern int optind, optopt, opterr;
- assert (argc > 0);
+ assert(argc > 0);
/* setup defaults */
prog_name = argv[0];
- env->T = 300.0; /* K */
- *Fmax = 1e5;
- *Xmax = 1e-6;
+ env->T = 300.0; /* K */
+ *Fmax = 1e5; /* N */
+ *Xmax = 1e-6; /* m */
*flags = 0;
*model = tension_models;
while ((c=getopt(argc, argv, options)) != -1) {
switch(c) {
- case 'T': env->T = atof(optarg); break;
- case 'C': env->T = atof(optarg)+273.15; break;
+ case 'T': env->T = safe_strtod(optarg, "-T"); break;
+ case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
case 'm':
tension_model = index_tension_model(n_tension_models, tension_models, optarg);
*model = tension_models+tension_model;
case 'a':
tension_models[tension_model].params = optarg;
break;
- case 'F': *Fmax = atof(optarg); break;
- case 'X': *Xmax = atof(optarg); break;
- case 'V': *flags |= VFLAG; break;
+ case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
+ case 'X': *Xmax = safe_strtod(optarg, "-X"); break;
+ case 'V': *flags |= VFLAG; break;
case '?':
fprintf(stderr, "unrecognized option '%c'\n", optopt);
/* fall through to default case */
void *string_create_const_k_param_t(char **param_strings)
{
- return create_const_k_param_t(atof(param_strings[0]));
+ return create_const_k_param_t(safe_strtod(param_strings[0], __FUNCTION__));
}
void destroy_const_k_param_t(void *p /* const_k_param_t * */)
void *string_create_bell_param_t(char **param_strings)
{
- return create_bell_param_t(atof(param_strings[0]), atof(param_strings[1]));
+ return create_bell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
+ safe_strtod(param_strings[1], __FUNCTION__));
}
void destroy_bell_param_t(void *p /* bell_param_t * */)
void *string_create_kbell_param_t(char **param_strings)
{
- return create_kbell_param_t(atof(param_strings[0]), atof(param_strings[1]));
+ return create_kbell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
+ safe_strtod(param_strings[1], __FUNCTION__));
}
void destroy_kbell_param_t(void *p /* kbell_param_t * */)
fprintf(in, "F = %g; T = %g\n", F, T);
sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
string_eval(in, out, input, 80, output);
- return atof(output);
+ return safe_strtod(output, __FUNCTION__);
}
double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
string_eval(in, out, input, 80, output);
- return atof(output);
+ return safe_strtod(output, __FUNCTION__);
}
double kramers_E(double F, environment_t *env, void *kramers_params, double x)
/* takes an array of 4 functions: {"(1,2),(2,0),..."} */
void *string_create_kramers_param_t(char **param_strings)
{
- return create_kramers_param_t(atof(param_strings[0]),
+ return create_kramers_param_t(safe_strtod(param_strings[0], __FUNCTION__),
param_strings[2],
param_strings[3],
param_strings[4],
//printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
// param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
void *E_params = string_create_spline_param_t(param_strings+1);
- return create_kramers_integ_param_t(atof(param_strings[0]),
- atof(param_strings[2]),
- atof(param_strings[3]),
- &spline_eval, E_params,
- destroy_spline_param_t);
+ return create_kramers_integ_param_t(
+ safe_strtod(param_strings[0], __FUNCTION__),
+ safe_strtod(param_strings[2], __FUNCTION__),
+ safe_strtod(param_strings[3], __FUNCTION__),
+ &spline_eval, E_params, destroy_spline_param_t);
}
void destroy_kramers_integ_param_t(void *params)
@
<<k model includes>>=
-#include <assert.h> /* assert() */
-#include <stdlib.h> /* NULL, malloc() */
-#include <math.h> /* HUGE_VAL, sqrt(), exp() */
-#include <stdio.h> /* fprintf(), stdout */
-#include <string.h> /* strlen(), strcpy() */
+#include <assert.h> /* assert() */
+#include <stdlib.h> /* NULL, malloc(), strto*() */
+#include <math.h> /* HUGE_VAL, sqrt(), exp() */
+#include <stdio.h> /* fprintf(), stdout */
+#include <string.h> /* strlen(), strcpy() */
+#include <errno.h> /* errno, ERANGE */
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include "global.h"
@
<<k model functions>>=
+<<safe strto*>>
<<null k functions>>
<<const k functions>>
<<bell k functions>>
@
<<k model check includes>>=
-#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
-#include <stdio.h> /* sprintf() */
-#include <assert.h> /* assert() */
-#include <math.h> /* exp() */
+#include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
+#include <stdio.h> /* sprintf() */
+#include <assert.h> /* assert() */
+#include <math.h> /* exp() */
+#include <errno.h> /* errno, ERANGE */
<<check includes>>
#include "global.h"
#include "k_model.h"
@
<<k model test suite>>=
+<<safe strto*>>
<<const k tests>>
<<bell k tests>>
<<k model suite definition>>
params[0] = knot[i];
p = string_create_const_k_param_t(params);
for ( F=Fm; F<FM; F+=dF ) {
- fail_unless(const_k(F, &env, p)==atof(knot[i]),
+ fail_unless(const_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
"const_k(%g, %g, &{%s}) = %g != %s",
F, T, knot[i], const_k(F, &env, p), knot[i]);
}
for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
params[0] = knot[i];
p = string_create_bell_param_t(params);
- fail_unless(bell_k(F, &env, p)==atof(knot[i]),
+ fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
"bell_k(%g, %g, &{%s,%s}) = %g != %s",
F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
destroy_bell_param_t(p);
char *knot[] = {"1","2","3","4","5","6"};
char *dx="1";
char *params[] = {knot[0], dx};
- double F= k_B*T/atof(dx);
+ double F= k_B*T/safe_strtod(dx, __FUNCTION__);
void *p = NULL;
environment_t env;
char param_string[80];
for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
params[0] = knot[i];
p = string_create_bell_param_t(params);
- CHECK_ERR(1e-6, atof(knot[i])*exp(1.0), bell_k(F, &env, p));
- //fail_unless(bell_k(F, &env, p)==atof(knot[i])*exp(1.0),
+ CHECK_ERR(1e-6, safe_strtod(knot[i], __FUNCTION__)*exp(1.0), bell_k(F, &env, p));
+ //fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__)*exp(1.0),
// "bell_k(%g, %g, &{%s,%s}) = %g != %g",
- // F, T, knot[i], dx, bell_k(F, &env, p), atof(knot[i])*exp(1.0));
+ // F, T, knot[i], dx, bell_k(F, &env, p), safe_strtod(knot[i], __FUNCTION__)*exp(1.0));
destroy_bell_param_t(p);
}
}
<<k model utility includes>>=
#include <stdlib.h>
#include <stdio.h>
-#include <string.h> /* strlen() */
-#include <assert.h> /* assert() */
+#include <string.h> /* strlen() */
+#include <assert.h> /* assert() */
+#include <errno.h> /* errno, ERANGE */
#include "global.h"
#include "parse.h"
#include "string_eval.h"
@
<<k model utility getopt functions>>=
+<<safe strto*>>
<<index k model>>
<<k model utility help>>
<<k model utility get options>>
extern char *optarg;
extern int optind, optopt, opterr;
- assert (argc > 0);
+ assert(argc > 0);
/* setup defaults */
prog_name = argv[0];
- env->T = 300.0; /* K */
+ env->T = 300.0; /* K */
*mode = M_K_OF_F;
*flags = 0;
*model = k_models;
- *Fmax = 1e-9;
+ *Fmax = 1e-9; /* N */
*special_xmax = 1e-8;
*special_xmin = 0.0;
while ((c=getopt(argc, argv, options)) != -1) {
switch(c) {
- case 'T': env->T = atof(optarg); break;
- case 'C': env->T = atof(optarg)+273.15; break;
+ case 'T': env->T = safe_strtod(optarg, "-T"); break;
+ case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
case 'k':
k_model = index_k_model(n_k_models, k_models, optarg);
*model = k_models+k_model;
case 'K':
k_models[k_model].params = optarg;
break;
- case 'm': *mode = M_K_OF_F; break;
- case 'M': *mode = M_SPECIAL; break;
- case 'F': *Fmax = atof(optarg); break;
- case 'x': *special_xmin = atof(optarg); break;
- case 'X': *special_xmax = atof(optarg); break;
- case 'V': *flags |= VFLAG; break;
+ case 'm': *mode = M_K_OF_F; break;
+ case 'M': *mode = M_SPECIAL; break;
+ case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
+ case 'x': *special_xmin = safe_strtod(optarg, "-x"); break;
+ case 'X': *special_xmax = safe_strtod(optarg, "-X"); break;
+ case 'V': *flags |= VFLAG; break;
case '?':
fprintf(stderr, "unrecognized option '%c'\n", optopt);
/* fall through to default case */
projects / sawsim.git / blobdiff