From 62eff8ac09814757097844c6eb242e38ae7a5e77 Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Fri, 5 Sep 2008 21:48:59 -0400 Subject: [PATCH] Added the option to make discrete steps in x. --- TODO | 2 +- src/sawsim.nw | 54 +++++++++++++++++++++++++++++++++++++++++---------- 2 files changed, 45 insertions(+), 11 deletions(-) diff --git a/TODO b/TODO index 9eaf623..a19485f 100644 --- a/TODO +++ b/TODO @@ -1 +1 @@ -fix arbitrary length scales and solution tolerances somhow +maybe histogram without dice rolling for identical domains? diff --git a/src/sawsim.nw b/src/sawsim.nw index c43a9b3..bb09ab4 100644 --- a/src/sawsim.nw +++ b/src/sawsim.nw @@ -213,20 +213,34 @@ int main(int argc, char **argv) <> list_t *domain_list=NULL; /* variables initialized in get_options() */ environment_t env={0}; - double P, dt_max, v; + double P, dt_max, v, xstep; int num_folded=0, unfolding; double x=0, dt, F; /* variables used in the simulation loop */ - get_options(argc, argv, &P, &dt_max, &v, &env, NUM_TENSION_MODELS, + get_options(argc, argv, &P, &dt_max, &v, &xstep, &env, NUM_TENSION_MODELS, tension_models, NUM_K_MODELS, k_models, &domain_list, &num_folded); setup(); if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Position (m)\tForce (N)\n"); if (flags & FLAG_OUTPUT_UNFOLDING_FORCES) printf("#Unfolding Force (N)\n"); while (num_folded > 0) { F = find_tension(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, unfolding); - dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, v); + if (xstep == 0) + dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, v); + else + dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, 0); unfolding=random_unfoldings(domain_list, F, dt, &env, &num_folded); if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\n", x, F); - x += v*dt; + if (xstep == 0) { + x += v*dt; + } else { + if (dt == dt_max) { /* step completed */ + fprintf(stderr, "step completed\n"); + x += xstep; + dt_max = xstep / v; + } else { /* still working on this step */ + fprintf(stderr, "partial step %g of %g\n", dt, dt_max); + dt_max -= dt; + } + } } destroy_domain_list(domain_list); free_accels(); @@ -238,6 +252,20 @@ int main(int argc, char **argv) [[determine_dt]] in Section \ref{sec.adaptive_dt}, and [[random_unfoldings]] in Section \ref{sec.unfolding_rate}. +The stretched distance is extended in one of two modes depending on +[[xstep]]. If [[xstep == 0]], then the pulling is continuous, at +least within the limits of the inherent discretization of a time +stepped approach. At any rate, the timestep is limited by the maximum +allowed timestep [[dt_max]] and the maximum allowed unfolding +probability in a given step [[P]]. In the second mode [[xstep > 0]], +and the end to end distance increases in discrete steps of that +length. The time between steps is chosen to maintain an average +unfolding velocity [[v]]. This approach is less work to simulate than +smooth pulling and also closer to the reality of many experiments, but +it is practically impossible to treat analytically. With both pulling +styles implemented, the effects of the discretization can be easily +calculated, bridging the gap between theory and experiment. + Environmental parameters important in determining reaction rates and tensions (e.g. temperature, pH) are stored in a single structure to facilitate extension to more complicated models in the future. @@ -609,7 +637,7 @@ k_model_getopt_t k_models[NUM_K_MODELS] = { \subsection{help} <>= -void help(char *prog_name, double P, double dt_max, double v, +void help(char *prog_name, double P, double dt_max, double v, double xstep, environment_t *env, int n_tension_models, tension_model_getopt_t *tension_models, int folded_tension_model, int unfolded_tension_model, @@ -624,6 +652,9 @@ void help(char *prog_name, double P, double dt_max, double v, printf("-P P\tTarget probability for dt (currently %g)\n", P); printf("-t dt\tMaximum allowed timestep dt (currently %g)\n", dt_max); printf("-v v\tPulling velocity v (currently %g nm/s)\n", v); + printf("-x dx\tPulling step size (currently %g m)\n", xstep); + printf("\tWhen dx = 0, the pulling is continuous\n"); + printf("\tWhen dx > 0, the pulling occurs in discrete steps\n"); printf("Environmental options:\n"); printf("-T T\tTemperature T (currently %g K)\n", env->T); printf("-C T\tYou can also set the temperature T in Celsius\n"); @@ -686,14 +717,14 @@ void help(char *prog_name, double P, double dt_max, double v, <>= void get_options(int argc, char **argv, - double *pP, double *pDt_max, double *pV, + double *pP, double *pDt_max, double *pV, double *pXstep, 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 **pDomain_list, int *pNum_folded) { char *prog_name = NULL; - char c, options[] = "P:t:v:T:C:m:a:M:A:k:K:F:U:Vh"; + char c, options[] = "P:t:v:x:T:C:m:a:M:A:k:K:F:U:Vh"; int ftension_model=0, utension_model=0, k_model=0; char *ftension_params=NULL, *utension_params=NULL; int i, n, ftension_group=0, utension_group=0; @@ -718,6 +749,7 @@ void get_options(int argc, char **argv, *pP = 1e-3; /* % pop per s */ *pDt_max = 0.001; /* s */ *pV = 1e-6; /* m/s */ + *pXstep = 0.0; /* m */ env->T = 300.0; /* K */ while ((c=getopt(argc, argv, options)) != -1) { @@ -725,6 +757,7 @@ void get_options(int argc, char **argv, case 'P': *pP = atof(optarg); break; case 't': *pDt_max = atof(optarg); break; case 'v': *pV = atof(optarg); break; + case 'x': *pXstep = atof(optarg); break; case 'T': env->T = atof(optarg); break; case 'C': env->T = atof(optarg)+273.15; break; case 'm': @@ -797,7 +830,7 @@ void get_options(int argc, char **argv, fprintf(stderr, "unrecognized option '%c'\n", optopt); /* fall through to default case */ default: - help(prog_name, *pP, *pDt_max, *pV, env, n_tension_models, tension_models, ftension_model, utension_model, n_k_models, k_models, k_model); + help(prog_name, *pP, *pDt_max, *pV, *pXstep, env, n_tension_models, tension_models, ftension_model, utension_model, n_k_models, k_models, k_model); exit(1); } } @@ -1059,6 +1092,9 @@ double search_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, //printf("overshot max_dt\n"); dtU = max_dt; } + if (v == 0) /* discrete stepping, so force is temporatily constant */ + return dtU; + /* set a lower bound on dt too */ dtL = 0.0; @@ -1122,8 +1158,6 @@ double determine_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, l assert(target_prob > 0.0); assert(target_prob < 1.0); assert(dt_max > 0.0); - /* .5 nm steps = v * dt */ - //return 0.5e-9/v; while (domain_list != NULL) { if (D(domain_list)->state == DS_FOLDED) { new_dt = search_dt(num_tension_handlers, tension_handlers, domain_list, env, x, target_prob, dt, v); -- 2.26.2