/* 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;
\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.transtion_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
projects / sawsim.git / commitdiff