<<k model getopt array definition>>
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();
[[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.
\subsection{help}
<<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,
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");
<<get options>>=
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;
*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) {
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':
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);
}
}
//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;
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);
projects / sawsim.git / commitdiff