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 0.9 added the [[-P]] option to sawsim, setting the target
+probability for the per-state transition $P_N$, not the per-domain
+transisiton $P_1$.
+
+Version 0.10 added the [[-F]] option to sawsim, setting a limit on the
+force change $dF$ in a single timestep for continuous pulling
+simulations. It also added the piston tension model (Section
+\ref{sec.piston_tension}), and adjusted the stiffness calculations to
+handle domains with undefined stiffness.
<<version definition>>=
-#define VERSION "0.9"
+#define VERSION "0.10"
@
\subsection{License}
/* variables initialized in get_options() */
list_t *state_list=NULL, *transition_list=NULL;
environment_t env={0};
- double P, t_max, dt_max, v, x_step;
+ double P, dF, t_max, dt_max, v, x_step;
state_t *stop_state=NULL;
/* variables used in the simulation loop */
double t=0, x=0, dt, F; /* time, distance, timestep, force */
- int transition=1; /* boolean marking a transition for tension optimization */
+ int transition=1; /* boolean marking a transition for tension and timestep optimization */
- get_options(argc, argv, &P, &t_max, &dt_max, &v, &x_step, &stop_state, &env,
- NUM_TENSION_MODELS, tension_models, NUM_K_MODELS, k_models,
- &state_list, &transition_list);
+ get_options(argc, argv, &P, &dF, &t_max, &dt_max, &v, &x_step,
+ &stop_state, &env, NUM_TENSION_MODELS, tension_models,
+ NUM_K_MODELS, k_models, &state_list, &transition_list);
setup();
if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Distance (m)\tForce (N)\tStiffness (N/m)\n");
if (flags & FLAG_OUTPUT_TRANSITION_FORCES) printf("#Force (N)\tInitial state\tFinal state\n");
while ((t_max > 0 && t <= t_max) || (stop_state != NULL && stop_state->num_domains > 0)) {
- F = find_tension(state_list, &env, x, transition);
+ F = find_tension(state_list, &env, x, transition, 0);
if (x_step == 0)
- dt = determine_dt(state_list, transition_list, &env, x, P, dt_max, v);
+ dt = determine_dt(state_list, transition_list, &env, F, x, dF, P, dt_max, v, transition);
else
- dt = determine_dt(state_list, transition_list, &env, x, P, dt_max, 0);
+ dt = determine_dt(state_list, transition_list, &env, F, x, dF, P, dt_max, 0, transition);
transition=random_transitions(transition_list, F, dt, &env);
if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness);
t += dt;
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 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
-\ref{sec.kbell}).
+occured). While we're messing with the tension handlers and if
+[[const_env==0]], 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 \ref{sec.kbell}).
<<find tension>>=
-double find_tension(list_t *state_list, environment_t *env, double x, int transition)
-{
+double find_tension(list_t *state_list, environment_t *env, double x,
+ int transition, int const_env)
+{ // TODO: !cont_env -> step_call, only alter env or statics if step_call==1
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;
- double F;
+ double F, *pStiffness;
int i;
#ifdef DEBUG
last_x = -1.0;
} /* else continue with the current statics */
- F = tension_balance(num_active_groups, per_group_handlers, per_group_inverse_handlers, (void **)per_group_data, last_x, x, &env->stiffness);
+ if (const_env == 1)
+ pStiffness = NULL;
+ else
+ pStiffness = &env->stiffness;
+
+ F = tension_balance(num_active_groups, per_group_handlers, per_group_inverse_handlers, (void **)per_group_data, last_x, x, pStiffness);
last_x = x;
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).
+ 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
$$
\subsection{help}
<<help>>=
-void help(char *prog_name, double P, double t_max, double dt_max, double v, double x_step,
+void help(char *prog_name, double P, double dF,
+ double t_max, double dt_max, double v, double x_step,
state_t *stop_state, environment_t *env,
int n_tension_models, tension_model_getopt_t *tension_models,
int n_k_models, k_model_getopt_t *k_models,
printf("-q state\tStop simulation when this state empties (- for no limit, currently %s)\n", stop_state ? stop_state->name : "-" );
printf("-t tmax\tMaximum simulated time (-1 for no limit, currently %g s)\n", t_max);
printf("-d dt\tMaximum allowed timestep dt (currently %g s)\n", dt_max);
- printf("-P P\tTarget probability for dt (currently %g)\n", P);
+ printf("-P P\tMaximum transition probability for dt (currently %g)\n", P);
+ printf("-F dF\tMaximum change in force over dt. Only relevant if dx > 0. (currently %g N)\n", dF);
printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
printf("-x dx\tPulling step size (currently %g m)\n", x_step);
printf("\tWhen dx = 0, the pulling is continuous.\n");
<<get options>>=
<<safe strto*>>
-void get_options(int argc, char **argv,
- double *pP, double *pT_max, double *pDt_max, double *pV,
+void get_options(int argc, char **argv, double *pP, double *pDF,
+ double *pT_max, double *pDt_max, double *pV,
double *pX_step, state_t **pStop_state, environment_t *env,
int n_tension_models, tension_model_getopt_t *tension_models,
int n_k_models, k_model_getopt_t *k_models,
list_t **pState_list, list_t **pTransition_list)
{
char *prog_name = NULL;
- char c, options[] = "q:t:d:P:v:x:T:C:s:N:k:DVh";
+ char c, options[] = "q:t:d:P:F:v:x:T:C:s:N:k:DVh";
int tension_model=0, k_model=0, initial_state=-1, final_state=-1;
char *state_name=NULL;
int i, n, tension_group=0;
*pT_max = -1; /* s, -1 removes this limit */
*pDt_max = 0.001; /* s */
*pP = 1e-3; /* % pop per s */
+ *pDF = 1e-12; /* N */
*pV = 1e-6; /* m/s */
*pX_step = 0.0; /* m */
env->T = 300.0; /* K */
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 'F': *pDF = safe_strtod(optarg, "-F"); 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;
fprintf(stderr, "unrecognized option '%c'\n", optopt);
/* fall through to default case */
default:
- help(prog_name, *pP, *pT_max, *pDt_max, *pV, *pX_step, *pStop_state, env, n_tension_models, tension_models, n_k_models, k_models, k_model, *pState_list);
+ help(prog_name, *pP, *pDF, *pT_max, *pDt_max, *pV, *pX_step,
+ *pStop_state, env, n_tension_models, tension_models,
+ n_k_models, k_models, k_model, *pState_list);
exit(1);
}
}
\label{sec.adaptive_dt_implementation}
$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_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
+chain model), so basing the timestep on the unfolding probability at
+the current tension is dangerous, and we need to search for a $dt$ for
+which $P_N(F(x+v*dt))<P_\text{max}$ (See Section
+\ref{sec.transition_rate} for a discussion of $P_N$). 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
step, and we expect the next step to be considerably longer than the
previous one.
<<search dt>>=
double search_dt(transition_t *transition,
list_t *state_list,
environment_t *env, double x,
- double target_prob, double max_dt, double v)
+ double max_prob, double max_dt, double v)
{ /* Returns a valid timestep dt in seconds for the given transition.
* Takes a list of all states, a pointer env to the environmental data,
- * a starting separation x in m, a target_prob between 0 and 1,
+ * a starting separation x in m, a max_prob between 0 and 1,
* max_dt in s, stretching velocity v in m/s.
*/
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 */
+ F = find_tension(state_list, env, x, 0, 1); /* 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 = P1(target_prob, N_INIT(transition)) / k; /* upper bound on dt */
+ dtU = P1(max_prob, N_INIT(transition)) / k; /* upper bound on dt */
if (dtU > max_dt) {
//printf("overshot max_dt\n");
dtU = max_dt;
/* The dt determined above may produce illegitimate forces or ks.
* Reduce the upper bound until we have valid ks. */
dt = dtU;
- F = find_tension(state_list, env, x+v*dt, 0);
+ F = find_tension(state_list, env, x+v*dt, 0, 1);
while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
dtU /= 2.0;
dt = dtU;
- F = find_tension(state_list, env, x+v*dt, 0);
+ F = find_tension(state_list, env, x+v*dt, 0, 1);
}
//printf("Try for dt = %g (F = %g)\n", dt, F);
k = accel_k(transition->k, F, env, transition->k_params);
while (k == -1.0) { /* reduce step until we hit a valid k */
dtU /= 2.0;
dt = dtU; /* hopefully, we can use the max dt, see if it works */
- F = find_tension(state_list, env, x+v*dt, 0);
+ F = find_tension(state_list, env, x+v*dt, 0, 1);
//printf("Try for dt = %g (F = %g)\n", dt, F);
k = accel_k(transition->k, F, env, transition->k_params);
//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 */
/* safe timestep back from x+dtU */
- dtUCur = P1(target_prob, N_INIT(transition)) / k;
+ dtUCur = P1(max_prob, N_INIT(transition)) / k;
if (dtUCur >= dt)
return dt; /* dtU is safe. */
/* dtU wasn't safe, lets see what would be. */
while (dtU > 1.1*dtL) { /* until the range is reasonably small */
dt = (dtU + dtL) / 2.0;
- F = find_tension(state_list, env, x+v*dt, 0);
+ F = find_tension(state_list, env, x+v*dt, 0, 1);
//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 = P1(target_prob, N_INIT(transition)) / k;
+ dtCur = P1(max_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;
@
-To determine $dt$ for an array of potentially different folded domains,
-we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
+Determine $dt$ for an array of potentially different folded domains.
+We apply [[search_dt()]] to each populated domain to find the maximum
+$dt$ the domain can handle. The final $dt$ is less than the
+individual domain $dt$s (to satisfy [[search_dt()]], see above). If
+$v>0$ (continuous pulling), $dt$ is also reduced to satisfy
+$|F(x+v*dt)-F(x)|<dF_\text{max}$, which limits errors due to our
+transition rate assumption that the force does not change appreciably
+over a single timestep.
<<determine dt>>=
<<search dt>>
double determine_dt(list_t *state_list, list_t *transition_list,
- environment_t *env, double x,
- double target_prob, double dt_max, double v)
+ environment_t *env, double F, double x, double max_dF,
+ double max_prob, double dt_max, double v, int transition)
{ /* Returns the timestep dt in seconds.
* Takes the state and transition lists, a pointer to the environment,
- * an extension x in m, target_prob between 0 and 1, dt_max in s,
- * and v in m/s. */
- double dt=dt_max, new_dt;
- assert(target_prob > 0.0); assert(target_prob < 1.0);
+ * the force F(x) in N, an extension x in m, max_dF in N,
+ * max_prob between 0 and 1, dt_max in s, v in m/s, and transition in {0,1}*/
+ double dt=dt_max, new_dt, F_new;
+ assert(max_dF > 0.0);
+ assert(max_prob > 0.0); assert(max_prob < 1.0);
assert(dt_max > 0.0);
assert(state_list != NULL);
assert(transition_list != NULL);
transition_list = head(transition_list);
+
+ if (v > 0) {
+ /* Ensure |dF| < max_dF */
+ F_new = find_tension(state_list, env, x+v*dt, transition, 1);
+ while (fabs(F_new - F) > max_dF) {
+ dt /= 2.0;
+ F_new = find_tension(state_list, env, x+v*dt, transition, 1);
+ }
+ }
+
+ /* Ensure all individual domains pass search_dt() */
while (transition_list != NULL) {
if (T(transition_list)->initial_state->num_domains > 0) {
- new_dt = search_dt(T(transition_list), state_list, env, x, target_prob, dt, v);
+ new_dt = search_dt(T(transition_list),state_list,env,x,max_prob,dt,v);
dt = MIN(dt, new_dt);
}
transition_list = transition_list->next;
transition_t unfold={0};
environment_t env = {0};
double F[]={0,1,2,3,4,5,6};
- double dt_max=3.0, Fnot=3.0, x=1.0, targ_p=0.1, v=1.0;
+ double dt_max=3.0, Fnot=3.0, x=1.0, max_p=0.1, max_dF=0.1, v=1.0;
list_t *state_list=NULL, *transition_list=NULL;
int i;
env.T = 300.0;
push(&state_list, &folded, 1);
push(&transition_list, &unfold, 1);
for( i=0; i < sizeof(F)/sizeof(double); i++) {
- //fail_unless(determine_dt(state_list, transition_list, &env, x, targ_p, dt_max, v)==1, NULL);
+ //fail_unless(determine_dt(state_list, transition_list, &env, x, max_p, max_dF, dt_max, v)==1, NULL);
}
}
END_TEST
<<x of x0>>
<<group stiffness function>>
<<chain stiffness function>>
+<<full chain stiffness function>>
<<tension balance function>>
@
lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
ub = x_xo_data.xi[0];
} else { /* x == last_x */
+#ifdef DEBUG
fprintf(stderr, "not moving\n");
- lb= x_xo_data.xi[0];
- ub= x_xo_data.xi[0];
+#endif
+ lb = x_xo_data.xi[0];
+ ub = x_xo_data.xi[0];
}
}
#ifdef DEBUG
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)
{
}
@
+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
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$.
+$$
+ F = \begin{cases}
+ 0 & \text{if $x<L$}, \\
+ \infty & \text{if $x>L$}.
+ \end{cases}
+$$
+
+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,
@
+We abuse [[x_of_xo()]] and [[oneD_solve()]] with this inverse
+definition (see Section \ref{sec.tension_balance}), because the
+$x(F=0)<=L$, but our function returns $x(F=0)=0$. This is fine when
+the total extension \emph{is} zero, but cheats a bit when there is
+some total extension. For any non-zero extension, the initial active
+group produces some tension (we hope. True for all our other tension
+models). This tension fully extends the piston group (of which there
+should be only one, since all piston states can get lumped into the
+same tension group.). If the total extension is $\le$ the target
+extension, the full extension is accurate, so the inverse was valid
+after all. If not, [[oneD_solve()]] will halve the extension of the
+first group, reducing the overall tension. For total extension less
+than the piston extension, this bisection continues for [[max_steps]],
+leaving $x_0$ reduced by a factor of $2^\text{max\_steps}$, a very
+small number. Because of this, the returned tension $F_0(x_0)$ is
+very small, even though the total extension found by [[x_of_xo()]]
+is still $>$ the piston length.
+
+While this works, it is clearly not the most elegant or robust
+solution. TODO: think of (and implement) a better procedure.
+
<<inverse piston tension handler declaration>>=
extern double inverse_piston_tension_handler(double x, void *pdata);
@
projects / sawsim.git / blobdiff