+% sawsim - simulating single molecule mechanical unfolding.
+% Copyright (C) 2008-2010 William Trevor King
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+% The author may be contacted at <wking@drexel.edu> on the Internet, or
+% write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
+% Philadelphia PA 19104, USA.
+
%%%%%%%%%%%%%%%%%%%%%%%%% Start LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%
% we have our own preamble,
% use `noweave -delay` to not print noweb's automatic one
\subsection{Related work}
-TODO References
+Sawim has been published\citep{king10}! It is, as far as I know, the
+only open-source Monte Carlo force spectroscopy simulator, but similar
+closed source simulators have been around since the seminal paper by
+\citet{rief97a}. In chronological order, major contributions have
+come from
+
+\begin{packed_item}
+ \item \citet{rief97a} \citeyear{rief97a}:
+ \item \citet{rief97b} \citeyear{rief97b}:
+ \item \citet{kellermayer97} \citeyear{kellermayer97}:
+ \item \citet{rief98} \citeyear{rief98}:
+ first paper to focus mainly on the simulation
+ \item \citet{oberhauser98} \citeyear{oberhauser98}:
+ \item \citet{carrion-vazquez99a} \citeyear{carrion-vazquez99a}:
+ \item \citet{best02} \citeyear{best02}:
+ pointed out large fit valley
+ \item \citet{zinober02} \citeyear{zinober02}:
+ introduced the scaffold effect
+ \item \citet{jollymore09} \citeyear{jollymore09}:
+ \item \citet{king10} \citeyear{king10}:
+ introduced sawsim
+\end{packed_item}
+
\subsection{About this document}
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}
<<license>>=
- sawsim - program for simulating single molecule mechanical unfolding.
- Copyright (C) 2008-2009, William Trevor King
+ sawsim - simulating single molecule mechanical unfolding.
+ Copyright (C) 2008-2010, William Trevor King
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License
- along with this program; if not, write to the Free Software
- Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
- 02111-1307, USA.
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
The author may be contacted at <wking@drexel.edu> on the Internet, or
write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
/* 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
$$
particular chain. For example, a damped spring thermalizes on a
timescale of order $\tau = 1/\gamma$ where the damping ratio $\gamma
\equiv \omega_0(-\zeta \pm \sqrt{\zeta^2-1}$, the natural angular
-frequency $omega_0 \equiv \sqrt{k/m}$, $\zeta \equiv b/2\sqrt{km}$,
+frequency $\omega_0 \equiv \sqrt{k/m}$, $\zeta \equiv b/2\sqrt{km}$,
and $b$ sets the drag $F_b = -bv$. For our cantilevers $\tau$ is
on the order of milliseconds, which is longer than a timestep.
However, the approximation is still reasonable, because there is
simulation itself, so we remove the actual model implementation to the
appendices.
-The tension models are defined in Appendix \ref{sec.tension_models}
+The tension models are defined in Appendix \ref{sec.tension_models},
and unfolding models are defined in Appendix \ref{sec.k_models}.
<<model globals>>=
#define k_B 1.3806503e-23 /* J/K */
+
@
\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);
}
}
#include "tension_model.h"
#include "parse.h"
#include "accel_k.h"
+
@
<<definitions>>=
-(received-expected)/expected, max_err, #received, \
expected, received); \
} while(0)
+
@
<<happens>>=
\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;
Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
<<check-parse.c>>=
+<<license comment>>
<<parse check includes>>
<<string comparison definition and globals>>
<<check relative error>>
<<parse suite definition>>=
Suite *test_suite (void)
{
- Suite *s = suite_create ("k model");
+ Suite *s = suite_create ("parse");
<<parse list string test case defs>>
- <<parse list string test case add>>
+ <<parse list string test case adds>>
return s;
}
@
<<parse list string test case defs>>=
TCase *tc_parse_list_string = tcase_create("parse list string");
@
-<<parse list string test case add>>=
+<<parse list string test case adds>>=
//tcase_add_test(tc_parse_list_string, test_next_delim_index);
tcase_add_test(tc_parse_list_string, test_parse_list_null);
tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
@
<<test suite>>=
-<<F tests>>
<<determine dt tests>>
<<happens tests>>
<<does domain transition tests>>
Suite *test_suite (void)
{
Suite *s = suite_create ("sawsim");
- <<F test case defs>>
<<determine dt test case defs>>
<<happens test case defs>>
<<does domain transition test case defs>>
<<randomly unfold domains test case defs>>
- <<F test case add>>
- <<determine dt test case add>>
- <<happens test case add>>
- <<does domain transition test case add>>
- <<randomly unfold domains test case add>>
+ <<determine dt test case adds>>
+ <<happens test case adds>>
+ <<does domain transition test case adds>>
+ <<randomly unfold domains test case adds>>
/*
tcase_add_checked_fixture(tc_strip_address,
}
@
-\subsection{F tests}
-<<F tests>>=
-<<worm-like chain tests>>
-@
-<<F test case defs>>=
-<<worm-like chain test case def>>
-@
-<<F test case add>>=
-<<worm-like chain test case add>>
-@
-
-<<worm-like chain tests>>=
-<<worm-like chain function>>
-
-START_TEST(test_wlc_at_zero)
-{
- double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
- fail_unless(abs(wlc(x, T, p, L)) < lim, \
- "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
- x, T, p, L, abs(wlc(x, T, p, L)), lim);
-}
-END_TEST
-
-START_TEST(test_wlc_at_half)
-{
- double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
- /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
- * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
- * = 0.25 * 3 = 3/4
- * linear term = x/L = 0.5
- * nonlinear + linear = 0.75 + 0.5 = 1.25
- * wlc = 10e21*1.25 = 12.5e21
- */
- fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
- "wlc(%g, %g, %g, %g) = %g != %g",
- x, T, p, L, wlc(x, T, p, L), 12.5e21);
-}
-END_TEST
-@
-<<worm-like chain test case def>>=
-TCase *tc_wlc = tcase_create("WLC");
-@
-
-<<worm-like chain test case add>>=
-tcase_add_test(tc_wlc, test_wlc_at_zero);
-tcase_add_test(tc_wlc, test_wlc_at_half);
-suite_add_tcase(s, tc_wlc);
-@
\subsection{Model tests}
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
<<determine dt test case defs>>=
TCase *tc_determine_dt = tcase_create("Determine dt");
@
-<<determine dt test case add>>=
+<<determine dt test case adds>>=
tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
suite_add_tcase(s, tc_determine_dt);
@
@
<<happens test case defs>>=
@
-<<happens test case add>>=
+<<happens test case adds>>=
@
<<does domain transition tests>>=
@
<<does domain transition test case defs>>=
@
-<<does domain transition test case add>>=
+<<does domain transition test case adds>>=
@
<<randomly unfold domains tests>>=
@
<<randomly unfold domains test case defs>>=
@
-<<randomly unfold domains test case add>>=
+<<randomly unfold domains test case adds>>=
@
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
Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
<<check-tension-balance.c>>=
+<<license comment>>
<<tension balance check includes>>
<<tension balance test suite>>
<<main check program>>
Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
<<check-list.c>>=
+<<license comment>>
<<list check includes>>
<<list test suite>>
<<main check program>>
Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
<<check-string-eval.c>>=
+<<license comment>>
<<string eval check includes>>
<<string comparison definition and globals>>
<<string eval test suite>>
Suite *s = suite_create ("string eval");
<<string eval test case defs>>
- <<string eval test case add>>
+ <<string eval test case adds>>
return s;
}
@
FILE *in, *out;
char input[80], output[80]={};
string_eval_setup(&in, &out);
- sprintf(input, "print ABCDefg\n");
+ sprintf(input, "print undefined_name__traceback_should_be_printed_to_stderr\n");
string_eval(in, out, input, 80, output);
string_eval_teardown(&in, &out);
}
<<string eval test case defs>>=
TCase *tc_string_eval = tcase_create("string_eval");
@
-<<string eval test case add>>=
+<<string eval test case adds>>=
tcase_add_test(tc_string_eval, test_setup_teardown);
tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
tcase_add_test(tc_string_eval, test_simple_eval);
<<accel-k.h>>=
#ifndef ACCEL_K_H
#define ACCEL_K_H
+<<license comment>>
double accel_k(k_func_t *k, double F, environment_t *env, void *params);
void free_accels();
#endif /* ACCEL_K_H */
<<accel k module makefile lines>>=
NW_SAWSIM_MODS += accel_k
-#CHECK_BINS += check_accel_k
-check_accel_k_MODS =
+CHECK_BINS += check_accel_k
+check_accel_k_MODS = interp k_model parse string_eval tavl
@
<<accel-k.c>>=
+<<license comment>>
#include <assert.h> /* assert() */
#include <stdlib.h> /* realloc(), free(), NULL */
#include "global.h" /* environment_t */
}
@
+\subsection{Unit tests}
+
+Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
+<<check-accel-k.c>>=
+<<license comment>>
+<<accel k check includes>>
+<<check relative error>>
+<<model globals>>
+<<accel k test suite>>
+<<main check program>>
+@
+
+<<accel k check includes>>=
+#include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
+#include <stdio.h> /* sprintf() */
+#include <assert.h> /* assert() */
+#include <math.h> /* exp() */
+#include <errno.h> /* errno, ERANGE */
+<<check includes>>
+#include "global.h"
+#include "k_model.h"
+#include "accel_k.h"
+
+@
+
+<<accel k test suite>>=
+<<accel k suite tests>>
+<<accel k suite definition>>
+@
+
+<<accel k suite definition>>=
+Suite *test_suite (void)
+{
+ Suite *s = suite_create ("accelerated k model");
+ <<accel k test case defs>>
+
+ <<accel k test case adds>>
+ return s;
+}
+@
+
+<<accel k test case defs>>=
+TCase *tc_accel_k = tcase_create("accelerated k model");
+@
+
+<<accel k test case adds>>=
+tcase_add_test(tc_accel_k, test_accel_k_accuracy);
+suite_add_tcase(s, tc_accel_k);
+@
+
+<<accel k suite tests>>=
+START_TEST(test_accel_k_accuracy)
+{
+ double k, k_accel, F, Fm=0.0, FM=300e-12, dF=0.5e-12;
+ k_func_t *k_func = bell_k;
+ char *param_strings[] = {"3.3e-4", "0.25e-9"};
+ void *params = string_create_bell_param_t(param_strings);
+ environment_t env = {};
+ env.T = 300.0;
+
+ for(F=Fm; F <= FM; F+=dF) {
+ k = k_func(F, &env, params);
+ k_accel = accel_k(k_func, F, &env, params);
+ CHECK_ERR(1e-8, k, k_accel);
+ }
+ destroy_bell_param_t(params);
+}
+END_TEST
+
+@
+
\section{Tension models}
\label{sec.tension_models}
<<tension model module makefile lines>>=
NW_SAWSIM_MODS += tension_model
-#CHECK_BINS += check_tension_model
-check_tension_model_MODS =
+CHECK_BINS += check_tension_model
+check_tension_model_MODS = list tension_balance
@
<<tension model utils makefile lines>>=
TENSION_MODEL_MODS = tension_model parse list tension_balance
$(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
notangle -Rtension-model-utils.c $< > $@
$(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
- gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
+ $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
clean_tension_model_utils :
rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
@ The pipe symbol [[|]] marks the prerequisits following it (in this
$$
where $L=Nl$ is the total length of the chain, and $\mathcal{L}(\alpha) \equiv \coth{\alpha} - \frac{1}{\alpha}$ is the Langevin function\citep{hatfield99}.
<<freely-jointed chain function>>=
-double langevin(double x, void *params)
+static double langevin(double x, void *params)
{
if (x==0) return 0;
return 1.0/tanh(x) - 1/x;
}
-<<one dimensional bracket declaration>>
-<<one dimensional solve declaration>>
-double inv_langevin(double x)
+
+static double inv_langevin(double x)
{
double lb=0.0, ub=1.0, ret;
oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
return ret;
}
-double fjc(double x, double T, double l, int N)
+static double fjc(double x, double T, double l, int N)
{
double L = l*(double)N;
if (x == 0) return 0;
//printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
return k_B*T/l * inv_langevin(x/L);
}
+
@
<<freely-jointed chain handler declaration>>=
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
@
+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);
@
}
@
+\subsection{Tension model unit tests}
+
+Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
+<<check-tension-model.c>>=
+<<license comment>>
+<<tension model check includes>>
+<<check relative error>>
+<<model globals>>
+<<tension model test suite>>
+<<main check program>>
+@
+
+<<tension model check includes>>=
+#include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
+#include <stdio.h> /* sprintf() */
+#include <assert.h> /* assert() */
+#include <math.h> /* exp() */
+#include <errno.h> /* errno, ERANGE */
+<<check includes>>
+#include "global.h"
+#include "list.h"
+#include "tension_balance.h" /* oneD_*() */
+#include "tension_model.h"
+@
+
+<<tension model test suite>>=
+<<safe strto*>>
+<<const tension tests>>
+<<hooke tests>>
+<<worm-like chain tests>>
+<<freely-jointed chain tests>>
+<<piston tests>>
+<<tension model suite definition>>
+@
+
+<<tension model suite definition>>=
+Suite *test_suite (void)
+{
+ Suite *s = suite_create ("tension model");
+ <<const tension test case defs>>
+ <<hooke test case defs>>
+ <<worm-like chain test case defs>>
+ <<freely-jointed chain test case defs>>
+ <<piston test case defs>>
+
+ <<const tension test case adds>>
+ <<hooke test case adds>>
+ <<worm-like chain test case adds>>
+ <<freely-jointed chain test case adds>>
+ <<piston test case adds>>
+ return s;
+}
+@
+
+\subsubsection{Constant}
+
+<<const tension test case defs>>=
+TCase *tc_const_tension = tcase_create("Constant tension");
+@
+
+<<const tension test case adds>>=
+tcase_add_test(tc_const_tension, test_const_tension_create_destroy);
+suite_add_tcase(s, tc_const_tension);
+@
+
+<<const tension tests>>=
+START_TEST(test_const_tension_create_destroy)
+{
+ char *tension[] = {"1","2.2", "3"};
+ char *params[] = {tension[0]};
+ void *p = NULL;
+ int i;
+ for( i=0; i < sizeof(tension)/sizeof(char *); i++) {
+ params[0] = tension[i];
+ p = string_create_const_tension_param_t(params);
+ destroy_const_tension_param_t(p);
+ }
+}
+END_TEST
+
+@ TODO: In order to test the constant tension handler itself, we'd
+have to construct a group.
+
+
+\subsubsection{Hooke}
+
+<<hooke test case defs>>=
+TCase *tc_hooke = tcase_create("Hooke");
+@
+
+<<hooke test case adds>>=
+tcase_add_test(tc_hooke, test_hooke_create_destroy);
+suite_add_tcase(s, tc_hooke);
+
+@
+
+<<hooke tests>>=
+START_TEST(test_hooke_create_destroy)
+{
+ char *k[] = {"1","2.2", "3"};
+ char *params[] = {k[0]};
+ void *p = NULL;
+ int i;
+ for( i=0; i < sizeof(k)/sizeof(char *); i++) {
+ params[0] = k[i];
+ p = string_create_hooke_param_t(params);
+ destroy_hooke_param_t(p);
+ }
+}
+END_TEST
+
+@ TODO: In order to test the Hooke tension handler itself, we'd
+have to construct a group.
+
+
+\subsubsection{Worm-like chain}
+
+<<worm-like chain test case defs>>=
+TCase *tc_wlc = tcase_create("WLC");
+@
+
+<<worm-like chain test case adds>>=
+tcase_add_test(tc_wlc, test_wlc_at_zero);
+tcase_add_test(tc_wlc, test_wlc_at_half);
+suite_add_tcase(s, tc_wlc);
+
+@
+
+<<worm-like chain tests>>=
+<<worm-like chain function>>
+START_TEST(test_wlc_at_zero)
+{
+ double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
+ fail_unless(abs(wlc(x, T, p, L)) < lim, \
+ "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
+ x, T, p, L, abs(wlc(x, T, p, L)), lim);
+}
+END_TEST
+
+START_TEST(test_wlc_at_half)
+{
+ double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
+ /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
+ * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
+ * = 0.25 * 3 = 3/4
+ * linear term = x/L = 0.5
+ * nonlinear + linear = 0.75 + 0.5 = 1.25
+ * wlc = 10e21*1.25 = 12.5e21
+ */
+ fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
+ "wlc(%g, %g, %g, %g) = %g != %g",
+ x, T, p, L, wlc(x, T, p, L), 12.5e21);
+}
+END_TEST
+
+@
+
+
+\subsubsection{Freely-jointed chain}
+
+<<freely-jointed chain test case defs>>=
+TCase *tc_fjc = tcase_create("FJC");
+@
+
+<<freely-jointed chain test case adds>>=
+tcase_add_test(tc_fjc, test_fjc_at_zero);
+tcase_add_test(tc_fjc, test_fjc_at_half);
+suite_add_tcase(s, tc_fjc);
+
+@
+
+<<freely-jointed chain tests>>=
+<<freely-jointed chain function>>
+START_TEST(test_fjc_at_zero)
+{
+ int N=10;
+ double T=1.0, l=1.0, p=0.1, x=0.0, lim=1e-30;
+ fail_unless(abs(fjc(x, T, l, N)) < lim, \
+ "abs(fjc(x=%g,T=%g,l=%g,N=%d)) = %g >= %g",
+ x, T, l, N, abs(fjc(x, T, l, N)), lim);
+}
+END_TEST
+
+START_TEST(test_fjc_at_half)
+{
+ int N = 10;
+ double T=1.0/k_B, l=1.0, x=5.0;
+ /* prefactor = k_B T / l = k_B (1.0/k_B) / 1.0 = 1.0 J/nm = 1.0e21 pN
+ * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
+ * = 0.25 * 3 = 3/4
+ * linear term = x/L = 0.5
+ * nonlinear + linear = 0.75 + 0.5 = 1.25
+ * fjc = 10e21*1.25 = 12.5e21
+ */
+ fail_unless(fjc(x, T, l, N)-12.5e21 < 1e16,
+ "fjc(%g, %g, %g, %d) = %g != %g",
+ x, T, l, N, fjc(x, T, l, N), 12.5e21);
+}
+END_TEST
+
+@
+
+
+\subsubsection{Piston}
+
+<<piston test case defs>>=
+TCase *tc_piston = tcase_create("Piston");
+@
+
+<<piston test case adds>>=
+tcase_add_test(tc_piston, test_piston_create_destroy);
+suite_add_tcase(s, tc_piston);
+
+@
+
+<<piston tests>>=
+START_TEST(test_piston_create_destroy)
+{
+ char *L[] = {"1","2.2", "3"};
+ char *params[] = {L[0]};
+ void *p = NULL;
+ int i;
+ for( i=0; i < sizeof(L)/sizeof(char *); i++) {
+ params[0] = L[i];
+ p = string_create_piston_tension_param_t(params);
+ destroy_piston_tension_param_t(p);
+ }
+}
+END_TEST
+
+@ TODO: In order to test the piston tension handler itself, we'd
+have to construct a group.
+
\section{Unfolding rate models}
\label{sec.k_models}
$(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
notangle -Rk-model-utils.c $< > $@
$(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
- gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
+ $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
clean_k_model_utils :
rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
@
<<bell k model getopt>>=
{"bell", "thermalized, two-state transitions", &bell_k, num_bell_params, bell_param_descriptions, (char *)bell_param_string, &string_create_bell_param_t, &destroy_bell_param_t}
@
-Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
+Initialized with Titin I27 parameters\citep{carrion-vazquez99a}.
% Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
<<kbell k model getopt>>=
{"kbell", "thermalized, two-state transitions, corrected for effective linker stiffness", &kbell_k, num_kbell_params, kbell_param_descriptions, (char *)kbell_param_string, &string_create_kbell_param_t, &destroy_kbell_param_t}
@
-Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
+Initialized with Titin I27 parameters\citep{carrion-vazquez99a}.
% Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
<<kramers integ k model getopt>>=
{"kramers_integ", "thermalized, diffusion-limited transitions", &kramers_integ_k, num_kramers_integ_params, kramers_integ_param_descriptions, (char *)kramers_integ_param_string, &string_create_kramers_integ_param_t, &destroy_kramers_integ_param_t}
@
-Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
+Initialized with Titin I27 parameters\citep{carrion-vazquez99a}.
% Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
\subsection{Unfolding model implementation}
Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
<<check-k-model.c>>=
+<<license comment>>
<<k model check includes>>
<<check relative error>>
<<model globals>>
<<kramers k tests>>=
@
-<<kramers k test case def>>=
+<<kramers k test case defs>>=
@
-<<kramers k test case add>>=
+<<kramers k test case adds>>=
@
<<k model function tests>>=
k_models[k_model].params);
printf("There are two output modes. In standard mode, k(F) is printed\n");
printf("For example:\n");
- printf(" #Force (N)\tk (% pop. per s)\n");
+ printf(" #Force (N)\tk (%% pop. per s)\n");
printf(" 123.456\t7.89\n");
printf(" ...\n");
printf("In special mode, the output depends on the model.\n");
$(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
$(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
$(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
+ $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
+ $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
$(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
mv $(BUILD_DIR)/sawsim.pdf $@
@
The sed is needed because notangle replaces tabs with eight spaces,
but [[make]] needs tabs.
<<makefile>>=
+# Customize compilation
+CC = gcc
+CFLAGS =
+LDFLAGS =
+
# Decide what directories to put things in
SRC_DIR = src
BUILD_DIR = build
# Modules (X.c, X.h) defined in the noweb file
NW_SAWSIM_MODS =
+# Check binaries (check_*) defined in the noweb file
CHECK_BINS =
# Modules defined outside the noweb file
FREE_SAWSIM_MODS = interp tavl
DOCS = sawsim.pdf
# Define the major targets
-all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
+all : ./Makefile all_bin all_doc
+
+.PHONY: all_bin all_doc
+all_bin : $(BINS:%=$(BIN_DIR)/%)
+all_doc : $(DOCS:%=$(DOC_DIR)/%)
view : $(DOC_DIR)/sawsim.pdf
xpdf $< &
$(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
$(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
$(BUILD_DIR)/global.h ./gmon.out
- $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
+ $(SHELL) -e -c "rm -rf $(BUILD_DIR) $(DOC_DIR)"
# Various builds of sawsim
$(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
- gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
+ $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
$(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
- gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
+ $(CC) $(CFLAGS) $(LDFLAGS) -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
$(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
- gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
+ $(CC) $(CFLAGS) $(LDFLAGS) -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
# Create the directories
$(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
$(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
$(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
$(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR)
- gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
+ $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
clean_check_sawsim :
rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
# ... and the modules
$$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
$$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
$$(BUILD_DIR)/global.h | $(BIN_DIR)
- gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
+ $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
$(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
$(CHECK_LIBS:%=-l%)
-# todo: clean up the required modules to
+# TODO: clean up the required modules too
$(CHECK_BINS:%=clean_%) :
rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
projects / sawsim.git / blobdiff