From 585d002200daaecbdf4d2681c2d0b68b63ba0a10 Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Sun, 13 Sep 2009 11:47:58 -0400 Subject: [PATCH] Added max dF "-F" option to sawsim. Also added const_env argument to find_tension(). As I state in a comment in find_tension(), it should really be a step_call/new_step/whatever option. --- src/sawsim.nw | 131 +++++++++++++++++++++++++++++++------------------- 1 file changed, 81 insertions(+), 50 deletions(-) diff --git a/src/sawsim.nw b/src/sawsim.nw index 87ed220..90ab17b 100644 --- a/src/sawsim.nw +++ b/src/sawsim.nw @@ -296,26 +296,26 @@ int main(int argc, char **argv) /* 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; @@ -460,21 +460,22 @@ into the [[tension_balance]] function. [[transition]] should be set 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}). <>= -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 @@ -501,7 +502,12 @@ double find_tension(list_t *state_list, environment_t *env, double x, int transi 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; @@ -843,7 +849,8 @@ k_model_getopt_t k_models[NUM_K_MODELS] = { \subsection{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, @@ -864,7 +871,8 @@ void help(char *prog_name, double P, double t_max, double dt_max, double v, doub 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"); @@ -928,15 +936,15 @@ void help(char *prog_name, double P, double t_max, double dt_max, double v, doub <>= <> -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; @@ -966,6 +974,7 @@ void get_options(int argc, char **argv, *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 */ @@ -991,6 +1000,7 @@ void get_options(int argc, char **argv, 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; @@ -1064,7 +1074,9 @@ void get_options(int argc, char **argv, 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); } } @@ -1322,32 +1334,33 @@ int happens(double probability) \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))>= 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; @@ -1361,11 +1374,11 @@ double search_dt(transition_t *transition, /* 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); @@ -1374,24 +1387,24 @@ double search_dt(transition_t *transition, 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; @@ -1406,26 +1419,44 @@ double search_dt(transition_t *transition, @ -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)|>= <> 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; @@ -2243,7 +2274,7 @@ START_TEST(test_determine_dt_linear_k) 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; @@ -2255,7 +2286,7 @@ START_TEST(test_determine_dt_linear_k) 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 -- 2.26.2