1 %%%%%%%%%%%%%%%%%%%%%%%%% Start LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%
2 % we have our own preamble,
3 % use `noweave -delay` to not print noweb's automatic one
4 % -*- mode: Noweb; noweb-code-mode: c-mode -*-
5 \documentclass[letterpaper, 11pt]{article}
8 \noweboptions{smallcode}
10 \usepackage{url} % needed for natbib for underscores in urls?
11 \usepackage[breaklinks,colorlinks]{hyperref} % for internal links, \phantomsection
12 % breaklinks breaks long links
13 % colorlinks colors link text instead of boxing it
14 \usepackage{amsmath} % for \text{}, \xrightarrow[sub]{super}
15 \usepackage{amsthm} % for theorem and proof environments
16 \newtheorem{theorem}{Theorem}
18 \usepackage{subfig} % Support side-by-side figures
19 \usepackage{pgf} % Fancy graphics
20 \usepackage{tikz} % A nice, inline PGF frontend
21 \usetikzlibrary{automata} % Graph-theory library
23 \usepackage{doc} % BibTeX
24 \usepackage[super,sort&compress]{natbib} % fancy citation extensions
25 % super selects citations in superscript mode
26 % sort&compress automatically sorts and compresses compound citations (\cite{a,b,...})
28 \usepackage[mathletters]{ucs} % use math fonts to allow Greek UTF-8 in the text
29 \usepackage[utf8x]{inputenc} % allow UTF-8 characters
31 %\bibliographystyle{ieeetr} % pick the bibliography style, short and sweet
32 %\bibliographystyle{plain} % pick the bibliography style, includes dates
33 \bibliographystyle{plainnat}
34 \defcitealias{sw:noweb}{{\tt noweb}}
35 \defcitealias{galassi05}{GSL}
36 \defcitealias{sw:check}{{\tt check}}
37 \defcitealias{sw:python}{Python}
42 \textheight 9.5in % leave a bit of extra room for page numbers (at top of page)
45 \setlength{\parindent}{0pt}
46 \setlength{\parskip}{5pt}
48 % For some reason, the \maketitle wants to go on it's own page
49 % so we'll just hardcode the following in our first page.
50 %\title{Sawsim: a sawtooth protein unfolding simulator}
51 %\author{W.~Trevor King}
54 \newcommand{\prog}{[[sawsim]]}
55 \newcommand{\lang}{[[C]]}
56 \newcommand{\p}[3]{\left#1 #2 \right#3} % e.g. \p({complicated expr})
57 \newcommand{\Th}{\ensuremath{\sup{\text{th}}}} % ordinal th
58 \newcommand{\U}[1]{\ensuremath{\text{ #1}}} % units: $1\U{m/s}$
60 % single spaced lists, from Alvin Alexander
61 % http://www.devdaily.com/blog/post/latex/control-line-spacing-in-itemize-enumerate-tags/
62 \newenvironment{packed_item}{
64 \setlength{\itemsep}{1pt}
65 \setlength{\parskip}{0pt}
66 \setlength{\parsep}{0pt}
70 \nwfilename{sawsim.nw}
75 Sawsim: a sawtooth protein unfolding simulator \\
83 %%%%%%%%%%%%%%%%%%%%%%%%%%%% End LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%%
85 \section{Introduction}
87 The unfolding of multi-globular protein strings is a stochastic
88 process. In the \prog\ program, we use Monte Carlo methods to
89 simulate this unfolding through an explicit time-stepping approach.
90 We develop a program in \lang\ to simulate probability of unfolding
91 under a constant extension velocity of a chain of globular domains.
93 In order to extract physical meaning from many mechanical unfolding
94 simulations, it is necessary to compare the experimental results
95 against the results of simulations using different models and
96 parameters. The model and parameters that best match the experimental
97 data provide the ‘best’ model for that experiment.
99 Previous mechanical unfolding studies have rolled their own
100 simulations for modelling their experiment, and there has been little
101 exchange and standardization of these simulations between groups.
102 The \prog\ program attempts to fill this gap, providing a flexible,
103 extensible, \emph{community} program, to make it easier for groups
104 to share the burden and increase the transparency of data analysis.
106 ‘Sawsim’ is blend of ‘sawtooth simulation’.
108 \subsection{Related work}
112 \subsection{About this document}
114 This document was written in \citetalias{sw:noweb} with code blocks in
115 \lang\ and comment blocks in \LaTeX.
117 \subsection{System Requirements and Dependencies}
120 \subsection{Change Log}
122 Version 0.0 used only the unfolded domain WLC to determine the tension
123 and had a fixed timestep $dt$.
125 Version 0.1 added adaptive timesteps, adjusting so that the unfolding
126 probability for a given domain was constant.
128 Version 0.2 added Kramers' model unfolding rate calculations, and a
129 switch to select Bell's or Kramers' model. This induced a major
130 rewrite, introducing the tension group abstraction, and a split to
131 multiple \lang\ source files. Also added a unit-test suites for
132 sanity checking, and switched to SI units throughout.
134 Version 0.3 added integrated Kramers' model unfolding rate
135 calculations. Also replaced some of my hand-coded routines with GNU
136 Scientific Library \citepalias{galassi05} calls.
138 Version 0.4 added Python evaluation for the saddle-point Kramers'
139 model unfolding rate calculations, so that the model functions could
140 be controlled from the command line. Also added interpolating lookup
141 tables to accelerate slow unfolding rate model evaluation, and command
142 line control of tension groups.
144 Version 0.5 added the stiffness environmental parameter and it's
145 associated unfolding models.
147 Version 0.6 generalized from two state unfolding to arbitrary
148 $N$-state rate reactions. This allows simulations of
149 unfolding-refolding, metastable intermediates, etc., but required
150 reorganizing the command line interface.
152 Version 0.7 added tension model inverses, which seems to reduce
153 computation time by about a factor of 2. I was expecting more, but
154 I'll take what I can get.
156 Version 0.8 fixed a minor bug in unfolding probability for
157 multi-domain groups. The probability of at least one domain unfolding
158 had been calculated as $P_N=NP_1$, rather than $P_N=1-(1-P_1)^N$.
159 However, for small $P$ the two are equivalent.
161 Version 0.9 the -P option to sawsim now sets the target probability
162 for the per-state transition $P_N$, not the per-domain transisiton
165 <<version definition>>=
166 #define VERSION "0.9"
172 sawsim - program for simulating single molecule mechanical unfolding.
173 Copyright (C) 2008-2009, William Trevor King
175 This program is free software; you can redistribute it and/or
176 modify it under the terms of the GNU General Public License as
177 published by the Free Software Foundation; either version 3 of the
178 License, or (at your option) any later version.
180 This program is distributed in the hope that it will be useful, but
181 WITHOUT ANY WARRANTY; without even the implied warranty of
182 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
183 See the GNU General Public License for more details.
185 You should have received a copy of the GNU General Public License
186 along with this program; if not, write to the Free Software
187 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
190 The author may be contacted at <wking@drexel.edu> on the Internet, or
191 write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
192 Philadelphia PA 19104, USA.
205 The unfolding system is stored as a chain of \emph{domains} (Figure
206 \ref{fig.domain_chain}). Domains can be folded, globular protein
207 domains, unfolded protein linkers, AFM cantilevers, or other
208 stretchable system components. The system is modeled as graph of
209 possible domain states with transitions between them (Figure
210 \ref{fig.domain_states}).
214 \subfloat[][]{\label{fig.domain_chain}
215 \begin{tikzpicture}[->,node distance=2cm,font=\footnotesize]
216 \tikzstyle{every state}=[fill,draw=red!50,very thick,fill=red!20]
217 \node[state] (A) {domain 1};
218 \node[state] (B) [below of=A] {domain 2};
219 \node[state] (C) [below of=B] {.~.~.};
220 \node[state] (D) [below of=C] {domain $N$};
221 \node (S) [below of=D] {Surface};
222 \node (E) [above of=A] {};
224 \path[-] (A) edge (B)
225 (B) edge node [right] {Tension} (C)
228 \path[->,green] (A) edge node [right,black] {Extension} (E);
232 \subfloat[][]{\label{fig.domain_states}
233 \begin{tikzpicture}[->,node distance=2.5cm,shorten <=1pt,shorten >=1pt,font=\footnotesize]
234 \tikzstyle{every state}=[fill,draw=blue!50,very thick,fill=blue!20]
235 \node[state] (A) {cantilever};
236 \node[state] (C) [below of=A] {transition};
237 \node[state] (B) [left of=C] {folded};
238 \node[state] (D) [right of=C] {unfolded};
240 \path (B) edge [bend left] node [above] {$k_1$} (C)
241 (C) edge [bend left] node [below] {$k_1'$} (B)
242 edge [bend left] node [above] {$k_2$} (D)
243 (D) edge [bend left] node [below] {$k_2'$} (C);
246 \caption{\subref{fig.domain_chain} Extending a chain of domains. One end
247 of the chain is fixed, while the other is extended at a constant
248 speed. The domains are coupled with rigid linkers, so the domains
249 themselves must stretch to accomodate the extension.
250 \subref{fig.domain_states} Each domain exists in a discrete state. At
251 each timestep, it may transition into another state following a
252 user-defined state matrix such as this one, showing a metastable
253 transition state and an explicit ``cantilever'' domain.}
257 Each domain contributes to the total tension, which depends on the
258 chain extension. The particular model for the tension as a function
259 of extension is handled generally with function pointers. So far the
260 following models have been implemented:
262 \item Null (Appendix \ref{sec.null_tension}),
263 \item Constant (Appendix \ref{sec.const_tension}),
264 \item Hookean spring (Appendix \ref{sec.hooke}),
265 \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
266 \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
269 The tension model and parameters used for a particular domain depend
270 on the domain's current state. The transition rates between states
271 are also handled generally with function pointers, with
274 \item Null (Appendix \ref{sec.null_k}),
275 \item Bell's model (Appendix \ref{sec.bell}),
276 \item Spring constant corrected Bell model (Appendix \ref{sec.kbell}),
277 \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
278 \item Kramers' saddle point model (Appendix \ref{sec.kramers}).
281 The unfolding of the chain domains is modeled in two stages. First
282 the tension of the chain is calculated. Then the tension (assumed to
283 be constant over the length of the chain, see Section
284 \ref{sec.timescales}), is applied to each domain, and used to
285 calculate the probability of that domain changing state in the next
286 timestep $dt$. Then the time is stepped forward, and the process is
287 repeated. The simulation is complete when a pre-set time limit
288 [[t_max]] is reached or a pre-selected state [[stop_state]] is empty.
290 int main(int argc, char **argv)
292 <<tension model getopt array definition>>
293 <<k model getopt array definition>>
295 /* variables initialized in get_options() */
296 list_t *state_list=NULL, *transition_list=NULL;
297 environment_t env={0};
298 double P, t_max, dt_max, v, x_step;
299 state_t *stop_state=NULL;
301 /* variables used in the simulation loop */
302 double t=0, x=0, dt, F; /* time, distance, timestep, force */
303 int transition=1; /* boolean marking a transition for tension optimization */
305 get_options(argc, argv, &P, &t_max, &dt_max, &v, &x_step, &stop_state, &env,
306 NUM_TENSION_MODELS, tension_models, NUM_K_MODELS, k_models,
307 &state_list, &transition_list);
310 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Distance (m)\tForce (N)\tStiffness (N/m)\n");
311 if (flags & FLAG_OUTPUT_TRANSITION_FORCES) printf("#Force (N)\tInitial state\tFinal state\n");
312 while ((t_max > 0 && t <= t_max) || (stop_state != NULL && stop_state->num_domains > 0)) {
313 F = find_tension(state_list, &env, x, transition);
315 dt = determine_dt(state_list, transition_list, &env, x, P, dt_max, v);
317 dt = determine_dt(state_list, transition_list, &env, x, P, dt_max, 0);
318 transition=random_transitions(transition_list, F, dt, &env);
319 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness);
324 if (dt == dt_max) { /* step completed */
327 } else { /* still working on this step */
332 destroy_state_list(state_list);
333 destroy_transition_list(transition_list);
337 @ The meat of the simulation is bundled into the three functions
338 [[find_tension]], [[determine_dt]], and [[random_transitions]].
339 [[find_tension]] is discussed in Section \ref{sec.find_tension},
340 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
341 [[random_transitions]] in Section \ref{sec.transition_rate}.
343 The stretched distance is extended in one of two modes depending on
344 [[xstep]]. If [[xstep == 0]], then the pulling is continuous, at
345 least within the limits of the inherent discretization of a time
346 stepped approach. At any rate, the timestep is limited by the maximum
347 allowed timestep [[dt_max]] and the maximum allowed unfolding
348 probability in a given step [[P]]. In the second mode [[xstep > 0]],
349 and the end to end distance increases in discrete steps of that
350 length. The time between steps is chosen to maintain an average
351 unfolding velocity [[v]]. This approach is less work to simulate than
352 smooth pulling and also closer to the reality of many experiments, but
353 it is practically impossible to treat analytically. With both pulling
354 styles implemented, the effects of the discretization can be easily
355 calculated, bridging the gap between theory and experiment.
357 Environmental parameters important in determining reaction rates and
358 tensions (e.g.\ temperature, pH) are stored in a single structure to
359 facilitate extension to more complicated models in the future.
360 <<environment definition>>=
361 typedef struct environment_struct {
371 #define DOUBLE_PRECISION 1e-12
373 <<environment definition>>
374 <<one dimensional function definition>>
375 <<create/destroy definitions>>
377 #endif /* GLOBAL_H */
379 Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
380 included multiple times.
382 \section{Simulation functions}
384 <<simulation functions>>=
388 <<does domain transition>>
389 <<randomly transition domains>>
393 \label{sec.find_tension}
395 Because the stretched system may be made up of several parts (folded
396 domains, unfolded domains, spring-like cantilever, \ldots), we will
397 assign the domains to states with tension models and groups. The
398 models are used to determine the tension of all the domain in a given
399 state following some model (Hookean spring, WLC, \ldots). The domains
400 are assumed to be commutative, so ordering is ignored. The
401 interactions between the groups are assumed to be linear, but the
402 interactions between domains of the same group need not be. Each
403 model has a tension handler function, which gives the tension $F$
404 needed to stretch a given group of domains a distance $x$.
406 Groups represent collections of states which the model should treat as
407 a single entity. To understand the purpose of groups, consider a
408 system with two unfolded domain states (e.g.\ for two different
409 proteins) that were both modeled as WLCs. If both unfolded states had
410 the same persistence length, it would be wise to assign them to the
411 same group. This would reduce both the computational cost of
412 balancing the tension and the error associated with the inter-group
413 interaction linearity. Note that group numbers only matter
414 \emph{within} model classes, since grouping domains with seperate
415 models doesn't make sense. Within designated groups, the tension
416 parameters for each domain are still checked for compatibility, so if
417 you accidentally assign incompatible domains to the same group you
418 should get an appropriate error message.
420 <<find tension definitions>>=
421 #define NUM_TENSION_MODELS 5
425 <<tension handler array global declaration>>=
426 extern one_dim_fn_t *tension_handlers[];
429 <<tension handler array global>>=
430 one_dim_fn_t *tension_handlers[] = {
432 &const_tension_handler,
440 <<tension model global declarations>>=
441 <<tension handler array global declaration>>
444 <<tension handler types>>=
445 typedef struct tension_handler_data_struct {
446 list_t *group_tension_params; /* one item for each domain in the group */
449 } tension_handler_data_t;
453 After sorting the chain into separate groups $G_i$, with tension
454 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
455 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
456 \forall i,j$. Note that there may be several states within each
457 group. [[find_tension]] is basically a wrapper to feed proper input
458 into the [[tension_balance]] function. [[transition]] should be set
459 to 0 if there were no transitions since the previous call to
460 [[find_tension]] to support some optimizations that assume a small
461 increase in tension between steps (only true if no transition has
462 occured). While we're messing with the tension handlers, we also grab
463 a measure of the current linker stiffness for the environmental
464 variable, which is needed by some of the unfolding rate models
465 (e.g.\ the linker-stiffness corrected Bell model, Appendix
468 double find_tension(list_t *state_list, environment_t *env, double x, int transition)
470 static int num_active_groups=0;
471 static one_dim_fn_t **per_group_handlers = NULL;
472 static one_dim_fn_t **per_group_inverse_handlers = NULL;
473 static void **per_group_data = NULL; /* array of pointers to tension_handler_data_t */
474 static double last_x;
475 tension_handler_data_t *tdata;
480 fprintf(stderr, "%s:\tfinding tension at distance %g%s\n", __FUNCTION__, x, transition ? " after a transition" : "");
483 if (transition != 0 || num_active_groups == 0) { /* setup statics */
484 /* get new statics, freeing the old ones if they aren't NULL */
485 get_active_groups(state_list, &num_active_groups, &per_group_handlers, &per_group_inverse_handlers, &per_group_data);
487 /* fill in the group handlers and params */
489 fprintf(stderr, "%s:\t%d groups, tension handler array %p, inverse tension handler array %p, handler data array %p, (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_inverse_handlers, per_group_data, per_group_data[0]);
491 for (i=0; i<num_active_groups; i++) {
492 tdata = (tension_handler_data_t *) per_group_data[i];
494 fprintf(stderr, "%s:\tactive group %d of %d (data %p, param list %p)", __FUNCTION__, i, num_active_groups, tdata, tdata->group_tension_params);
495 list_print(stderr, tdata->group_tension_params, " parameters");
498 /* tdata->persist continues unchanged */
501 } /* else continue with the current statics */
503 F = tension_balance(num_active_groups, per_group_handlers, per_group_inverse_handlers, (void **)per_group_data, last_x, x, &env->stiffness);
509 @ For the moment, we will restrict our group boundaries to a single
510 dimension, so $\sum_i x_i = x$, our total extension (it is this
511 restriction that causes the groups to interact linearly). We'll also
512 restrict our extensions to all be positive. With these restrictions,
513 the problem of balancing the tensions reduces to solving a system of
514 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
515 the number of active groups. In general this can be fairly
516 complicated, but without much loss of practicality we can restrict
517 ourselves to strictly monotonically increasing, non-negative tension
518 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
519 simpler. For example, it guarantees the existence of a unique, real
520 solution for finite forces. See Appendix \ref{sec.tension_balance}
521 for the tension balancing implementation.
523 Each group must define a parameter structure if the tension model is
525 a function to create the parameter structure,
526 a function to destroy the parameter structure, and
527 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as its data argument.
528 The tension-specific parameter structures are contained in the domain
529 group field. For optimization, a group may also define an inverse
530 tension function as an optimization. See the Section
531 \ref{sec.model_selection} for a discussion on the command line
532 framework and inverse function discussion. See the worm-like chain
533 implementation for an example (Section \ref{sec.wlc}).
535 The major limitation of our approach is that all of our tension models
536 are Markovian (memory-less), which is not really a problem for the
537 chains (see \citet{evans99} for WLC thermalization rate calculations).
539 \subsection{Transition rate}
540 \label{sec.transition_rate}
542 Each state transition is modeled as unimolecular, first order reaction
544 \text{State 1} \xrightarrow{k} \text{State 2}
546 With the rate constant $k$ defined by
548 \frac{d[\text{State 2}]}{dt} = -\frac{d[\text{State 1}]}{dt} = k [\text{State 1}].
550 So the probability of a given domain transitioning in a short time
556 For a group of $N$ identical domains, and therefore identical
557 unfolding probabilities $P_1$, the probability of $n$ domains
558 unfolding in a given timestep is given by the binomial distribution
560 P(n) = \frac{N!}{n!(N-n)!}P_1^n(1-P_1)^(N-n).
562 The probability that \emph{none} of these domains unfold is then
566 so the probability that \emph{at least one} domain unfolds is
570 Note that for small $P$,
572 p(n>0) = 1-(1-NP_1+\mathcal{O}(P_1^2)) = NP_1 - \mathcal{O}(P^2)
575 This conversion from $P_1$ to $P_N$ is handled by the [[pN]] macro
577 /* find multi-domain transition rate from the single domain rate */
578 #define PN(P,N) (1.0-pow(1.0-(P), (N)))
582 We can also define a macro to find the approprate $dt$ to achieve a
583 target $P_N$ with a given $k$ and $N$.
585 P_N &= 1-(1-P_1)^N \\
586 1-P_1 &= (1-P_N)^{1/N} \\
587 P_1 &= 1-(1-P_N)^{1/N}
590 #define P1(PN,N) (1.0-pow(1.0-(PN), 1.0/(N)))
593 We take some time to discuss the meaning of $p(n>1)$
594 (i.e. multi-unfolding timesteps) in Section \ref{sec.timescales}.
596 <<does domain transition>>=
597 int domain_transitions(double F, double dt, environment_t *env, transition_t *transition)
598 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to transition structure */
600 k = accel_k(transition->k, F, env, transition->k_params);
601 //(*transition->k)(F, env, domain->k_params);
602 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
603 return happens(PN(k*dt, N_INIT(transition)));
605 @ [[happens]] is a random decision making function defined in Appendix \ref{sec.utils}.
607 <<randomly transition domains>>=
608 int random_transitions(list_t *transition_list, double tension,
609 double dt, environment_t *env)
610 { /* returns 1 if there was a transition and 0 otherwise */
611 while (transition_list != NULL) {
612 if (T(transition_list)->initial_state->num_domains > 0
613 && domain_transitions(tension, dt, env, T(transition_list))) {
614 if (flags & FLAG_OUTPUT_TRANSITION_FORCES)
615 fprintf(stdout, "%g\t%s\t%s\n", tension, T(transition_list)->initial_state->name, T(transition_list)->final_state->name);
616 T(transition_list)->initial_state->num_domains--;
617 T(transition_list)->final_state->num_domains++;
618 return 1; /* our one domain has transitioned, stop looking for others */
620 transition_list = transition_list->next;
626 \subsection{Timescales}
627 \label{sec.timescales}
629 The simulation assumes that chain equilibration is instantaneous
630 relative to the simulation timestep $dt$, so that tension is uniform
631 along the chain. The quality of this assumption depends on your
632 particular chain. For example, a damped spring thermalizes on a
633 timescale of order $\tau = 1/\gamma$ where the damping ratio $\gamma
634 \equiv \omega_0(-\zeta \pm \sqrt{\zeta^2-1}$, the natural angular
635 frequency $omega_0 \equiv \sqrt{k/m}$, $\zeta \equiv b/2\sqrt{km}$,
636 and $b$ sets the drag $F_b = -bv$. For our cantilevers $\tau$ is
637 on the order of milliseconds, which is longer than a timestep.
638 However, the approximation is still reasonable, because there is
639 rarely unfolding during the cantilever return. You could, of course,
640 take cantilever, WLC, etc.\ response times into effect, but that
641 would significantly complicate a simulation that seems to work
642 well enough without bothering :p. For WLC and FJC relaxation times,
643 see the Appendix of \citet{evans99} and \citet{kroy07}.
645 Besides assuming our timestep is much greater than equilibration
646 times, we also force it to be short enough so that only one domain may
647 unfold in a given timestep. For an unfolding timescale of $dt_u =
648 1/k$, we adapt our timesteps to keep $dt \ll dt_u$, so the probability
649 of two domains unfolding in the same timestep is small (see
650 [[determine_dt]] in Section \ref{sec.adaptive_dt} and [[P]] in
651 [[main()]] in Section \ref{sec.main} set by [[-P]] in
652 [[get_options()]] in Section \ref{sec.get_options}). This approach
653 breaks down as the adaptive timestep scheme approaches the transition
654 timestep $dt_t$, but $dt_t \sim 1\U{$\mu$s}$ for Ig27-like proteins
655 \citep{klimov00}, so this shouldn't be a problem. To reassure
656 yourself, you can ask the simulator to print the smallest timestep
657 that was used with TODO.
659 Even if the likelyhood of two domains unfolding in the same timestep
660 is small, if you run long enough simulations it will eventually occur.
661 In this case, we assume that $dt_t \ll dt$, so even if two domains
662 would unfold if we stuck to our unfolding rate $k$ for an entire
663 timestep $dt$, in ``reality'' only one of those domains unfolds.
664 Because we do not know when in the timestep the unfolding took place,
665 we move on to the next timestep without checking for further
666 unfoldings. This ``unchecked timestep portion'' should not contribute
667 significant errors to the calculation, because if $dt$ was small
668 enough that unfolding was unlikely at the pre-unfolding force, it
669 should be even more unlikely at the post-unfolding force, especially
670 over only a fraction of the timestep. The only dangerous cases for
671 the current implementation are those where unfolding intermediates are
672 much less stable than their precursor states, in which case an
673 unfolding event during the remaining timestep may actually be likely.
674 A simple workaround for such cases is to pick the value for [[P]] to
675 be small enough that the $dt \ll dt_u$ assumption survives even under
676 a transition populating the unstable intermediate.
678 \subsection{Adaptive timesteps}
679 \label{sec.adaptive_dt}
681 We'd like to pick $dt$ so the probability of changing state
682 (e.g.\ unfolding) in the next timestep is small. If we don't adapt our
683 timestep, we also risk breaking our approximation that $F(x' \in [x,
684 x+v \cdot dt]) = F(x)$ or that only one domain changes state in a
685 given timestep. The simulation would have been accurate for
686 sufficiently small fixed timesteps, but adaptive timesteps will allow
687 us to move through low tension regions in fewer steps, leading to a
688 more efficient simulation.
690 Consider the two state folded $\rightarrow$ unfolded transition.
691 Because $F(x)$ is monotonically increasing between unfoldings,
692 excessively large timesteps will lead to erroneously small unfolding
693 rates (an thus, higher average unfolding force).
695 The actual adaptive timestep implementation is not particularly
696 interesting, since we are only required to reduce $dt$ to somewhere
697 below a set threshold, so I've removed it to Appendix
698 \ref{sec.adaptive_dt_implementation}.
704 The models provide the physics of the simulation, but the simulation
705 allows interchangeable models, and we are currently focusing on the
706 simulation itself, so we remove the actual model implementation to the
709 The tension models are defined in Appendix \ref{sec.tension_models}
710 and unfolding models are defined in Appendix \ref{sec.k_models}.
713 #define k_B 1.3806503e-23 /* J/K */
717 \section{Command line interface}
719 <<get option functions>>=
721 <<index tension model>>
726 \subsection{Model selection}
727 \label{sec.model_selection}
729 The main difficulty with the command line interface in \prog\ is
730 developing an intuitive method for accessing the possibly large number
731 of available models. We'll treat this generally by defining an array
732 of available models, containing their important parameters.
734 <<get options definitions>>=
735 <<model getopt definitions>>
738 <<create/destroy definitions>>=
739 typedef void *create_data_func_t(char **param_strings);
740 typedef void destroy_data_func_t(void *param_struct);
743 <<model getopt definitions>>=
744 <<tension model getopt definitions>>
745 <<k model getopt definitions>>
748 In order to choose models, we need a method of assembling a reaction
749 graph such as that in Figure \ref{fig.domain_states} and population
750 the graph with a set of domains. First we declare the domain states
753 [[-s <state-name>,<model>[,<group>] -S <model-params> [-N <initial-population]]] \\
757 [[-s unfolded,wlc,1 -S 1e-8,4e-10 -N 5]]
759 This creates the state named [[unfolded]], using the WLC model and one
760 as a group number (defaults to 0, see Section \ref{sec.find_tension}).
761 The tension parameters are then set to [[1e-8,4e-10]], which in the
762 case of the WLC are the contour and persistence lengths respectively.
763 Finally we populate the state by adding five domains to the state.
764 The name of the state is importent for identifying states later.
766 Once the domains have been created and populated, you can added
767 transition rates following
769 [[-k <initial-state>,<final-state>,<model> -K <model-params>]] \\
773 [[-k folded,unfolded,bell -K 3.3e-4,0.25e-9]]
775 This creates a reaction going from the [[folded]] state to the
776 [[unfolded]] state, with the rate constant given by the Bell model
777 (Appendix \ref{sec.bell}). The unfolding parameters are then set to
778 [[3.3e-4,0.25e-9]], which in the case of the Bell model are the
779 unforced rate and transition state distance respectively.
781 \subsubsection{Tension}
783 To access the various tension models from the command line, we define
784 a structure that contains all the neccessary information about the
786 <<tension model getopt definitions>>=
787 typedef struct tension_model_getopt_struct {
788 one_dim_fn_t *handler; /* fn implementing the model on a group */
789 one_dim_fn_t *inverse_handler; /* fn implementing inverse on group */
790 char *name; /* name string identifying the model */
791 char *description; /* a brief description of the model */
792 int num_params; /* number of extra params for the model */
793 char **param_descriptions; /* descriptions of extra params */
794 char *params; /* default values for the extra params */
795 create_data_func_t *creator; /* param-string -> model-data-struct fn */
796 destroy_data_func_t *destructor; /* fn to free the model data structure */
797 } tension_model_getopt_t;
798 @ Set [[inverse_handler]] to [[NULL]] if you wish to invert using
799 bisection. Obviously, this will be slower than computing an
800 analytical inverse and more prone to rounding errors, but it may be
801 easier if a closed-form inverse does not exist.
803 <<tension model getopt array definition>>=
804 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
805 <<null tension model getopt>>,
806 <<constant tension model getopt>>,
807 <<hooke tension model getopt>>,
808 <<worm-like chain tension model getopt>>,
809 <<freely-jointed chain tension model getopt>>
813 \subsubsection{Unfolding rate}
815 <<k model getopt definitions>>=
816 #define NUM_K_MODELS 6
818 typedef struct k_model_getopt_struct {
823 char **param_descriptions;
825 create_data_func_t *creator;
826 destroy_data_func_t *destructor;
830 <<k model getopt array definition>>=
831 k_model_getopt_t k_models[NUM_K_MODELS] = {
832 <<null k model getopt>>,
833 <<const k model getopt>>,
834 <<bell k model getopt>>,
835 <<kbell k model getopt>>,
836 <<kramers k model getopt>>,
837 <<kramers integ k model getopt>>
844 void help(char *prog_name, double P, double t_max, double dt_max, double v, double x_step,
845 state_t *stop_state, environment_t *env,
846 int n_tension_models, tension_model_getopt_t *tension_models,
847 int n_k_models, k_model_getopt_t *k_models,
848 int k_model, list_t *state_list)
853 if (state_list != NULL)
854 state = S(tail(state_list));
856 printf("usage: %s [options]\n", prog_name);
857 printf("Version %s\n\n", VERSION);
858 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
860 printf("Simulation options:\n");
862 printf("-q state\tStop simulation when this state empties (- for no limit, currently %s)\n", stop_state ? stop_state->name : "-" );
863 printf("-t tmax\tMaximum simulated time (-1 for no limit, currently %g s)\n", t_max);
864 printf("-d dt\tMaximum allowed timestep dt (currently %g s)\n", dt_max);
865 printf("-P P\tTarget probability for dt (currently %g)\n", P);
866 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
867 printf("-x dx\tPulling step size (currently %g m)\n", x_step);
868 printf("\tWhen dx = 0, the pulling is continuous.\n");
869 printf("\tWhen dx > 0, the pulling occurs in discrete steps.\n");
871 printf("Environmental options:\n");
872 printf("-T T\tTemperature T (currently %g K)\n", env->T);
873 printf("-C T\tYou can also set the temperature T in Celsius\n");
875 printf("State creation options:\n");
876 printf("-s name,model[[,group],params]\tCreate a new state.\n");
877 printf("Once you've set up the state, you may populate it with domains\n");
878 printf("-N n\tAdd n domains in the current state (%s)\n", state ? state->name : "-");
880 printf("Transition creation options:\n");
881 printf("-k initial_state,final_state,model[,params]\tCreate a new transition\n");
883 printf("To double check your reaction graph, you can produce graphviz dot output\n");
884 printf("describing the current states and transitions.\n");
885 printf("-D\tPrint dot description of states and transitions and exit.\n");
887 printf("Simulation output mode options:\n");
888 printf("There are two output modes. In standard mode, only the transition\n");
889 printf("events are printed. For example:\n");
890 printf(" #Force (N)\tInitial state\tFinal state\n");
891 printf(" 123.456e-12\tfolded\tunfolded\n");
893 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
894 printf(" #Distance (m)\tForce (N)\tStiffness (N/m)\n");
895 printf(" 0.001\t0.0023\t0.002\n");
897 printf("-V\tChange output to verbose mode.\n");
898 printf("-h\tPrint this help and exit.\n");
901 printf("Tension models:\n");
902 for (i=0; i<n_tension_models; i++) {
903 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
904 for (j=0; j < tension_models[i].num_params ; j++)
905 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
906 printf("\t\tdefaults: %s\n", tension_models[i].params);
908 printf("Unfolding rate models:\n");
909 for (i=0; i<n_k_models; i++) {
910 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
911 for (j=0; j < k_models[i].num_params ; j++)
912 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
913 printf("\t\tdefaults: %s\n", k_models[i].params);
916 printf("Examples:\n");
917 printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded. Stop once there are no more folded domains.\n");
918 printf(" $ %s -v1e-6 -s cantilever,hooke,0.05 -N1 -s folded,null -N8 -s 'unfolded,wlc,{0.39e-9,28e-9}' -k 'folded,unfolded,bell,{3.3e-4,0.25e-9}' -q folded\n", prog_name);
919 printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded. Stop after 0.25 seconds\n");
920 printf(" $ %s -v1e-6 -s cantilever,hooke,0.05 -N1 -s 'folded,wlc,{1e-8,4e-10}' -N8 -s 'unfolded,wlc,1,{0.4e-9,28.1e-9}' -k 'folded,unfolded,bell,{5e-5,0.225e-9}' -t0.25\n", prog_name);
924 \subsection{Get options}
925 \label{sec.get_options}
929 void get_options(int argc, char **argv,
930 double *pP, double *pT_max, double *pDt_max, double *pV,
931 double *pX_step, state_t **pStop_state, environment_t *env,
932 int n_tension_models, tension_model_getopt_t *tension_models,
933 int n_k_models, k_model_getopt_t *k_models,
934 list_t **pState_list, list_t **pTransition_list)
936 char *prog_name = NULL;
937 char c, options[] = "q:t:d:P:v:x:T:C:s:N:k:DVh";
938 int tension_model=0, k_model=0, initial_state=-1, final_state=-1;
939 char *state_name=NULL;
940 int i, n, tension_group=0;
941 list_t *temp_list=NULL;
944 void *model_params; /* for INIT_MODEL() */
945 int num_param_args; /* for INIT_MODEL() */
946 char **param_args; /* for INIT_MODEL() */
948 extern int optind, optopt, opterr;
953 for (i=0; i<n_tension_models; i++) {
954 fprintf(stderr, "tension model %i %s:\t model's handler %p\t inverse handler %p\n",
955 i, tension_models[i].name, tension_models[i].handler, tension_models[i].inverse_handler);
956 assert(tension_models[i].handler == tension_handlers[i]);
961 flags = FLAG_OUTPUT_TRANSITION_FORCES;
964 *pT_max = -1; /* s, -1 removes this limit */
965 *pDt_max = 0.001; /* s */
966 *pP = 1e-3; /* % pop per s */
967 *pV = 1e-6; /* m/s */
968 *pX_step = 0.0; /* m */
969 env->T = 300.0; /* K */
971 *pTransition_list = NULL;
973 while ((c=getopt(argc, argv, options)) != -1) {
976 if (STRMATCH(optarg, "-")) {
979 temp_list = head(*pState_list);
980 while (temp_list != NULL) {
981 if (STRMATCH(S(temp_list)->name, optarg)) {
982 *pStop_state = S(temp_list);
985 temp_list = temp_list->next;
989 case 't': *pT_max = safe_strtod(optarg, "-t"); break;
990 case 'd': *pDt_max = safe_strtod(optarg, "-d"); break;
991 case 'P': *pP = safe_strtod(optarg, "-P"); break;
992 case 'v': *pV = safe_strtod(optarg, "-v"); break;
993 case 'x': *pX_step = safe_strtod(optarg, "-x"); break;
994 case 'T': env->T = safe_strtod(optarg, "-T"); break;
995 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
997 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
998 assert(num_strings >= 2 && num_strings <= 4);
1000 state_name = strings[0];
1001 if (state_index(*pState_list, state_name) != -1) {
1002 fprintf(stderr, "Invalid option '%s', repeated state name '%s'\n", optarg, state_name);
1005 tension_model = index_tension_model(n_tension_models, tension_models, strings[1]);
1006 if (num_strings == 4)
1007 tension_group = safe_strtoi(strings[2], optarg);
1010 if (num_strings >= 3)
1011 tension_models[tension_model].params = strings[num_strings-1];
1013 fprintf(stderr, "%s:\t%s -> state %s, model %s, params %s, group %d\n", __FUNCTION__, optarg, state_name, tension_models[tension_model].name, tension_models[tension_model].params, tension_group);
1015 INIT_MODEL(state_name, &tension_models[tension_model], tension_models[tension_model].params, model_params);
1016 push(pState_list, create_state(state_name,
1017 tension_models[tension_model].handler,
1018 tension_models[tension_model].inverse_handler,
1021 tension_models[tension_model].destructor,
1023 *pState_list = tail(*pState_list);
1025 free_parsed_list(num_strings, strings);
1028 n = safe_strtoi(optarg, "-N");
1030 fprintf(stderr, "%s:\tadding %d domains to %s\n", __FUNCTION__, n, S(*pState_list)->name);
1032 S(*pState_list)->num_domains += n;
1035 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
1036 assert(num_strings >= 3 && num_strings <= 4);
1037 initial_state = state_index(*pState_list, strings[0]);
1038 final_state = state_index(*pState_list, strings[1]);
1039 k_model = index_k_model(n_k_models, k_models, strings[2]);
1041 fprintf(stderr, "%s:\t%s -> %s reaction from state %s (%d) -> state %s (%d)\n", __FUNCTION__, optarg, strings[2], strings[0], initial_state, strings[1], final_state);
1043 assert(initial_state != final_state);
1044 if (num_strings == 4)
1045 k_models[k_model].params = strings[3];
1046 INIT_MODEL(state_name, &k_models[k_model], k_models[k_model].params, model_params);
1047 push(pTransition_list, create_transition(list_element(*pState_list, initial_state),
1048 list_element(*pState_list, final_state),
1049 k_models[k_model].k,
1051 k_models[k_model].destructor), 1);
1052 free_parsed_list(num_strings, strings);
1055 print_state_graph(stdout, *pState_list, *pTransition_list);
1059 flags = FLAG_OUTPUT_FULL_CURVE;
1062 fprintf(stderr, "unrecognized option '%c'\n", optopt);
1063 /* fall through to default case */
1065 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);
1069 /* check the arguments */
1070 assert(*pState_list != NULL); /* no domains! */
1071 assert(*pP > 0.0); assert(*pP < 1.0);
1072 assert(*pDt_max > 0.0);
1074 assert(env->T > 0.0);
1076 crosslink_groups(*pState_list, *pTransition_list);
1082 <<index tension model>>=
1083 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
1086 for (i=0; i<n_models; i++)
1087 if (STRMATCH(models[i].name, name))
1089 if (i == n_models) {
1090 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
1098 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
1101 for (i=0; i<n_models; i++)
1102 if (STRMATCH(models[i].name, name))
1104 if (i == n_models) {
1105 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
1112 <<initialize model definition>>=
1113 /* requires int num_param_args and char **param_args in the current scope
1115 * INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
1116 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
1118 #define INIT_MODEL(role, model, param_string, param_pointer) \
1120 parse_list_string((param_string), SEP, '{', '}', \
1121 &num_param_args, ¶m_args); \
1122 if (num_param_args != (model)->num_params) { \
1124 "%s model %s expected %d params,\n", \
1125 (role), (model)->name, (model)->num_params); \
1127 "not the %d params in '%s'\n", \
1128 num_param_args, (param_string)); \
1129 assert(num_param_args == (model)->num_params); \
1131 if ((model)->creator) \
1132 param_pointer = (*(model)->creator)(param_args); \
1134 param_pointer = NULL; \
1135 free_parsed_list(num_param_args, param_args); \
1139 First we define some safe string parsing functions. Currently these
1140 abort on any irregularity, but in the future they could print error
1141 messages, etc. [[static]] because the functions are currently
1142 defined in each [[*.c]] file that uses them, but perhaps they should
1143 be moved to [[utils.h/c]] or some such instead\ldots
1146 static int safe_strtoi(const char *s, const char *description) {
1150 ret = strtol(s, &endp, 10);
1151 if (endp[0] != '\0') { //strlen(endp) > 0) {
1152 fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1153 endp, s, description);
1154 assert(1==0); //strlen(endp) == 0);
1156 fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1159 } else if (errno == ERANGE) {
1160 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1161 assert(errno != ERANGE);
1166 static double safe_strtod(const char *s, const char *description) {
1170 ret = strtod(s, &endp);
1171 if (endp[0] != '\0') { //strlen(endp) > 0) {
1172 fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1173 endp, s, description);
1174 assert(1==0); //strlen(endp) == 0);
1176 fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1179 } else if (errno == ERANGE) {
1180 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1181 assert(errno != ERANGE);
1190 \addcontentsline{toc}{section}{Appendicies}
1192 \section{sawsim.c details}
1196 The general layout of our simulation code is:
1207 We include [[math.h]], so don't forget to link to the libm with `-lm'.
1209 #include <assert.h> /* assert() */
1210 #include <stdlib.h> /* malloc(), free(), rand(), strto*() */
1211 #include <stdio.h> /* fprintf(), stdout */
1212 #include <string.h> /* strlen, strtok() */
1213 #include <math.h> /* exp(), M_PI, sqrt() */
1214 #include <sys/types.h> /* pid_t (returned by getpid()) */
1215 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
1216 #include <errno.h> /* errno, ERANGE (for safe_strto*()) */
1219 #include "tension_balance.h"
1220 #include "k_model.h"
1221 #include "tension_model.h"
1223 #include "accel_k.h"
1227 <<version definition>>
1228 <<flag definitions>>
1229 <<max/min definitions>>
1230 <<string comparison definition and globals>>
1231 <<initialize model definition>>
1232 <<get options definitions>>
1233 <<state definitions>>
1234 <<transition definitions>>
1243 <<create/destroy state>>
1244 <<destroy state list>>
1246 <<create/destroy transition>>
1247 <<destroy transition list>>
1248 <<print state graph>>
1250 <<simulation functions>>
1251 <<boilerplate functions>>
1254 <<boilerplate functions>>=
1256 <<get option functions>>
1262 srand(getpid()*time(NULL)); /* seed rand() */
1266 <<flag definitions>>=
1267 /* in octal b/c of prefixed '0' */
1268 #define FLAG_OUTPUT_FULL_CURVE 01
1269 #define FLAG_OUTPUT_TRANSITION_FORCES 02
1273 static unsigned long int flags = 0;
1276 \subsection{Utilities}
1279 <<max/min definitions>>=
1280 #define MAX(a,b) ((a)>(b) ? (a) : (b))
1281 #define MIN(a,b) ((a)<(b) ? (a) : (b))
1284 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
1285 <<string comparison definition and globals>>=
1286 // Check if two strings match, return 1 if they do
1287 static char *temp_string_A;
1288 static char *temp_string_B;
1289 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
1290 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
1291 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
1292 /* +1 to also compare the '\0' */
1295 We also define a macro for our [[check]] unit testing
1296 <<check relative error>>=
1297 #define CHECK_ERR(max_err, expected, received) \
1299 fail_unless( (received-expected)/expected < max_err, \
1300 "relative error %g >= %g in %s (Expected %g, received %g)", \
1301 (received-expected)/expected, max_err, #received, \
1302 expected, received); \
1303 fail_unless(-(received-expected)/expected < max_err, \
1304 "relative error %g >= %g in %s (Expected %g, received %g)", \
1305 -(received-expected)/expected, max_err, #received, \
1306 expected, received); \
1311 int happens(double probability)
1313 assert(probability >= 0.0); assert(probability <= 1.0);
1314 return (double)rand()/RAND_MAX < probability; /* TODO: replace with GSL rand http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html x*/
1319 \subsection{Adaptive timesteps}
1320 \label{sec.adaptive_dt_implementation}
1322 $F(x)$ increases with $x$, possibly exploding (e.g.\ in the worm-like
1323 chain model), so basing the timestep on the the unfolding probability
1324 at the current tension is dangerous, and we need to search for a $dt$
1325 for which $P_N(F(x+v*dt)) < P_\text{target}$. There are two cases to
1326 consider. In the most common, no domains have unfolded since the last
1327 step, and we expect the next step to be slightly shorter than the
1328 previous one. In the less common, domains did unfold in the last
1329 step, and we expect the next step to be considerably longer than the
1332 double search_dt(transition_t *transition,
1334 environment_t *env, double x,
1335 double target_prob, double max_dt, double v)
1336 { /* Returns a valid timestep dt in seconds for the given transition.
1337 * Takes a list of all states, a pointer env to the environmental data,
1338 * a starting separation x in m, a target_prob between 0 and 1,
1339 * max_dt in s, stretching velocity v in m/s.
1341 double F, k, P, dtCur, dtU, dtUCur, dtL, dt;
1343 /* get upper bound using the current position */
1344 F = find_tension(state_list, env, x, 0); /* BUG. repeated calculation */
1345 //printf("Start with x = %g (F = %g)\n", x, F);
1346 k = accel_k(transition->k, F, env, transition->k_params);
1347 //printf("x %g\tF %g\tk %g\n", x, F, k);
1348 dtU = P1(target_prob, N_INIT(transition)) / k; /* upper bound on dt */
1350 //printf("overshot max_dt\n");
1353 if (v == 0) /* discrete stepping, so force is temporarily constant */
1356 /* set a lower bound on dt too */
1359 /* The dt determined above may produce illegitimate forces or ks.
1360 * Reduce the upper bound until we have valid ks. */
1362 F = find_tension(state_list, env, x+v*dt, 0);
1363 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1366 F = find_tension(state_list, env, x+v*dt, 0);
1368 //printf("Try for dt = %g (F = %g)\n", dt, F);
1369 k = accel_k(transition->k, F, env, transition->k_params);
1370 /* returned k may be -1.0 */
1371 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1372 while (k == -1.0) { /* reduce step until we hit a valid k */
1374 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1375 F = find_tension(state_list, env, x+v*dt, 0);
1376 //printf("Try for dt = %g (F = %g)\n", dt, F);
1377 k = accel_k(transition->k, F, env, transition->k_params);
1378 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1380 assert(dtU > 1e-14); /* timestep to valid k too small */
1381 /* safe timestep back from x+dtU */
1382 dtUCur = P1(target_prob, N_INIT(transition)) / k;
1384 return dt; /* dtU is safe. */
1386 /* dtU wasn't safe, lets see what would be. */
1387 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1388 dt = (dtU + dtL) / 2.0;
1389 F = find_tension(state_list, env, x+v*dt, 0);
1390 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1391 k = accel_k(transition->k, F, env, transition->k_params);
1392 dtCur = P1(target_prob, N_INIT(transition)) / k;
1393 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1394 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1396 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1397 dtU = dt; /* ... stepping out only dtCur would be safe */
1400 break; /* dtCur = dt */
1402 return MAX(dtUCur, dtL);
1407 To determine $dt$ for an array of potentially different folded domains,
1408 we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
1411 double determine_dt(list_t *state_list, list_t *transition_list,
1412 environment_t *env, double x,
1413 double target_prob, double dt_max, double v)
1414 { /* Returns the timestep dt in seconds.
1415 * Takes the state and transition lists, a pointer to the environment,
1416 * an extension x in m, target_prob between 0 and 1, dt_max in s,
1418 double dt=dt_max, new_dt;
1419 assert(target_prob > 0.0); assert(target_prob < 1.0);
1420 assert(dt_max > 0.0);
1421 assert(state_list != NULL);
1422 assert(transition_list != NULL);
1423 transition_list = head(transition_list);
1424 while (transition_list != NULL) {
1425 if (T(transition_list)->initial_state->num_domains > 0) {
1426 new_dt = search_dt(T(transition_list), state_list, env, x, target_prob, dt, v);
1427 dt = MIN(dt, new_dt);
1429 transition_list = transition_list->next;
1436 \subsection{State data}
1437 \label{sec.state_data}
1439 The domains exist in one of an array of possible states. Each state
1440 has a tension model which determines it's force/extension curve
1441 (possibly [[null]] for rigid states, see Appendix
1442 \ref{sec.null_tension}).
1444 Groups are, for example, for two domain states with WLC tensions; one
1445 with $p=0.4\U{nm}$, $L=10\U{nm}$ and the other with $p=0.4\U{nm}$,
1446 $L=28\U{nm}$. Since the persistence lengths are the same, you would
1447 like to add the contour lengths of all the domains in \emph{both}
1448 states, and plug that total contour length into a single WLC
1449 evaluation (see Section \ref{sec.find_tension}).
1450 <<state definitions>>=
1451 typedef struct state_struct {
1452 char *name; /* name string */
1453 one_dim_fn_t *tension_handler; /* tension handler */
1454 one_dim_fn_t *inverse_tension_handler; /* inverse tension handler */
1455 int tension_group; /* for grouping similar tension states */
1456 void *tension_params; /* pointer to tension parameters */
1457 destroy_data_func_t *destroy_tension;
1458 int num_domains; /* number of domains currently in state */
1459 /* cached lookups generated from the above data */
1460 int num_out_trans; /* length of out_trans array */
1461 struct transition_struct **out_trans; /* transitions leaving this state */
1462 int num_grouped_states; /* length of grouped states array */
1463 struct state_struct **grouped_states; /* states grouped with this one (doesn't include self) */
1466 /* get the state data for the current list node */
1467 #define S(list) ((state_t *)(list)->d)
1469 @ [[tension_params]] is a pointer to the parameters used by the
1470 handler function when determining the tension. [[destroy_tension]]
1471 points to a function for freeing [[tension_params]]. We it in the
1472 state structure so that [[destroy_state]] and will have an easier job
1475 [[create_]] and [[destroy_state]] are simple wrappers around
1476 [[malloc]] and [[free]].
1477 <<create/destroy state>>=
1478 state_t *create_state(char *name,
1479 one_dim_fn_t *tension_handler,
1480 one_dim_fn_t *inverse_tension_handler,
1482 void *tension_params,
1483 destroy_data_func_t *destroy_tension,
1486 state_t *ret = (state_t *)malloc(sizeof(state_t));
1487 assert(ret != NULL);
1488 /* make a copy of the name, so we aren't dependent on the original */
1489 ret->name = (char *)malloc(sizeof(char)*(strlen(name)+1));
1490 assert(ret->name != NULL);
1491 strcpy(ret->name, name); /* we know ret->name is long enough */
1493 ret->tension_handler = tension_handler;
1494 ret->inverse_tension_handler = inverse_tension_handler;
1495 ret->tension_group = tension_group;
1496 ret->tension_params = tension_params;
1497 ret->destroy_tension = destroy_tension;
1498 ret->num_domains = num_domains;
1500 ret->num_out_trans = 0;
1501 ret->out_trans = NULL;
1502 ret->num_grouped_states = 0;
1503 ret->grouped_states = NULL;
1506 fprintf(stderr, "generated state %s (%p) with tension handler %p, inverse handler %p, params %p, group %d.\n", ret->name, ret, ret->tension_handler, ret->inverse_tension_handler, ret->tension_params, ret->tension_group);
1511 void destroy_state(state_t *state)
1515 fprintf(stderr, "freeing state %s (%p)\n", state->name, state);
1519 if (state->destroy_tension)
1520 (*state->destroy_tension)(state->tension_params);
1521 if (state->out_trans)
1522 free(state->out_trans);
1523 if (state->grouped_states)
1524 free(state->grouped_states);
1531 We'll be storing the states in a list (see Appendix \ref{sec.lists}),
1532 so we also define a convenience function for cleaning up.
1533 <<destroy state list>>=
1534 void destroy_state_list(list_t *state_list)
1536 state_list = head(state_list);
1537 while (state_list != NULL)
1538 destroy_state((state_t *) pop(&state_list));
1543 We can also define a convenient way to get a state index from it's name.
1545 int state_index(list_t *state_list, char *name)
1548 state_list = head(state_list);
1549 while (state_list != NULL) {
1550 if (STRMATCH(S(state_list)->name, name)) return i;
1551 state_list = state_list->next;
1559 \subsection{Transition data}
1560 \label{sec.transition_data}
1562 Once you've created states (Sections \ref{sec.main},
1563 \ref{sec.model_selection}, and \ref{sec.state_data}), you need to
1564 create transitions between them.
1565 <<transition definitions>>=
1566 typedef struct transition_struct {
1567 state_t *initial_state;
1568 state_t *final_state;
1569 k_func_t *k; /* transition rate model */
1570 void *k_params; /* pointer to k parameters */
1571 destroy_data_func_t *destroy_k;
1574 /* get the transition data for the current list node */
1575 #define T(list) ((transition_t *)(list)->d)
1577 /* get the number of initial state population for a transition state */
1578 #define N_INIT(transition) ((transition)->initial_state->num_domains)
1582 @ [[k]] is a pointer to the function determining the transition rate
1583 for a given tension. [[k_params]] and [[destroy_k]] have similar
1584 roles with regards to [[k]] as [[tension_params]] and
1585 [[destroy_tension]] have with regards to [[tension_handler]] (see
1586 Section \ref{sec.state_data}).
1588 [[create_]] and [[destroy_transition]] are simple wrappers around
1589 [[malloc]] and [[free]].
1590 <<create/destroy transition>>=
1591 transition_t *create_transition(state_t *initial_state, state_t *final_state,
1594 destroy_data_func_t *destroy_k)
1596 transition_t *ret = (transition_t *)malloc(sizeof(transition_t));
1597 assert(ret != NULL);
1598 assert(initial_state != NULL);
1599 assert(final_state != NULL);
1600 ret->initial_state = initial_state;
1601 ret->final_state = final_state;
1603 ret->k_params = k_params;
1604 ret->destroy_k = destroy_k;
1606 fprintf(stderr, "generated transition %p from %s to %s, params %p\n", ret, ret->initial_state->name, ret->final_state->name, ret->k_params);
1611 void destroy_transition(transition_t *transition)
1615 fprintf(stderr, "freeing transition %p (%s to %s)\n", transition, transition->initial_state->name, transition->final_state->name);
1617 if (transition->destroy_k)
1618 (*transition->destroy_k)(transition->k_params);
1625 We'll be storing the transitions in a list (see Appendix
1626 \ref{sec.lists}), so we also define a convenience function for
1628 <<destroy transition list>>=
1629 void destroy_transition_list(list_t *transition_list)
1631 transition_list = head(transition_list);
1632 while (transition_list != NULL)
1633 destroy_transition((transition_t *) pop(&transition_list));
1638 \subsection{Graphviz output}
1639 \label{sec.graphviz_output}
1641 <<print state graph>>=
1642 void print_state_graph(FILE *file, list_t *state_list, list_t *transition_list)
1644 state_list = head(state_list);
1645 transition_list = head(transition_list);
1646 fprintf(file, "digraph states {\n node [shape=record,color=black,fillcolor=white,style=filled];\n\n");
1648 fprintf(file, " n%d [label=\"%s\"]\n", state_index(state_list, S(state_list)->name), S(state_list)->name);
1649 if (state_list->next == NULL) break;
1650 state_list = state_list->next;
1652 fprintf(file, "\n");
1653 while (transition_list != NULL) {
1654 fprintf(file, " n%d -> n%d;\n", state_index(state_list, T(transition_list)->initial_state->name), state_index(state_list, T(transition_list)->final_state->name));
1655 transition_list = transition_list->next;
1657 fprintf(file, "}\n");
1661 \subsection{Domain model and group handling}
1663 <<group functions>>=
1664 <<crosslink groups>>
1665 <<get active groups>>
1668 Scan through a list of states and construct an array of all the
1670 <<get active groups>>=
1671 void get_active_groups(list_t *state_list,
1672 int *pNum_active_groups,
1673 one_dim_fn_t ***pPer_group_handlers,
1674 one_dim_fn_t ***pPer_group_inverse_handlers,
1675 void ***pPer_group_data /* pointer to array of pointers to tension_handler_data_t */)
1677 void **active_groups=NULL;
1678 one_dim_fn_t *handler, *old_handler, *inverse_handler, *old_inverse_handler;
1679 one_dim_fn_t **per_group_handlers=NULL, **per_group_inverse_handlers=NULL;
1680 void **per_group_data=NULL; /* array of pointers to tension_handler_data_t */
1681 int i, j, num_domains, group, old_group, num_active_groups=0;
1682 list_t *active_states_list=NULL;
1683 tension_handler_data_t *tdata=NULL;
1686 /* get all the active groups in a list */
1687 state_list = head(state_list);
1689 fprintf(stderr, "%s:\t", __FUNCTION__);
1690 list_print(stderr, state_list, "states");
1692 while (state_list != NULL) {
1693 state = S(state_list);
1694 handler = state->tension_handler;
1695 inverse_handler = state->inverse_tension_handler;
1696 group = state->tension_group;
1697 num_domains = state->num_domains;
1698 if (list_index(active_states_list, state) == -1) {
1699 /* we haven't added this state (or it's associated group) yet */
1700 if (num_domains > 0 && handler != NULL) { /* add the state */
1701 num_active_groups++;
1702 push(&active_states_list, state, 1);
1704 fprintf(stderr, "%s:\tactive %s state (%p) with %d domain(s)\n", __FUNCTION__, state->name, state, state->num_domains);
1706 for (i=0; i < state->num_grouped_states; i++) {
1707 if (state->grouped_states[i]->num_domains > 0) {
1708 /* add this related (and active) state now */
1709 assert(state->grouped_states[i]->tension_handler == handler);
1710 assert(state->grouped_states[i]->tension_group == group);
1711 push(&active_states_list, state->grouped_states[i], 1);
1713 fprintf(stderr, "%s:\tactive grouped %s state (%p) with %d domain(s)\n", __FUNCTION__, state->grouped_states[i]->name, state->grouped_states[i], state->grouped_states[i]->num_domains);
1717 } /* else empty state or NULL handler */
1718 } /* else state already added */
1719 state_list = state_list->next;
1722 fprintf(stderr, "%s:\t", __FUNCTION__);
1723 list_print(stderr, active_states_list, "active states");
1726 assert(num_active_groups <= list_length(active_states_list));
1727 assert(num_active_groups > 0);
1729 /* allocate the arrays we'll be returning */
1730 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1731 assert(per_group_handlers != NULL);
1732 per_group_inverse_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1733 assert(per_group_inverse_handlers != NULL);
1734 per_group_data = (void **)malloc(sizeof(void *)*num_active_groups);
1735 assert(per_group_data != NULL);
1737 fprintf(stderr, "%s:\tallocating per group data %p:", __FUNCTION__, per_group_data);
1739 for (i=0; i<num_active_groups; i++) { /* we might want to recycle any group data structs allocated in the previous round */
1740 per_group_data[i] = malloc(sizeof(tension_handler_data_t));
1741 assert(per_group_data[i] != NULL);
1743 fprintf(stderr, " %d,%p", i, per_group_data[i]);
1747 fprintf(stderr, "\n");
1750 /* populate the arrays we'll be returning */
1751 active_states_list = head(active_states_list);
1753 old_inverse_handler = NULL;
1754 j = -1; /* j is the current active _group_ index */
1755 while (active_states_list != NULL) {
1756 state = (state_t *) pop(&active_states_list);
1757 handler = state->tension_handler;
1758 inverse_handler = state->inverse_tension_handler;
1759 group = state->tension_group;
1760 num_domains = state->num_domains;
1761 if (handler != old_handler || inverse_handler != old_inverse_handler || group != old_group) {
1762 j++; /* move to the next group */
1763 tdata = (tension_handler_data_t *) per_group_data[j];
1764 per_group_handlers[j] = handler;
1765 per_group_inverse_handlers[j] = inverse_handler;
1766 tdata->group_tension_params = NULL;
1768 tdata->persist = NULL; /* we might want to grab any appropriate data from the old persist here... */
1771 fprintf(stderr, "%s:\tadding active state %s to group %d (%p:%p:%d) with %d domain(s)\n", __FUNCTION__, state->name, j, handler, inverse_handler, group, num_domains);
1773 for (i=0; i<num_domains; i++) { /* add num_domains copies of the tension params */
1774 push(&tdata->group_tension_params, state->tension_params, 1);
1777 fprintf(stderr, "%s:\tpushed %d copies of %p onto data %p's param list %p\n", __FUNCTION__, num_domains, state->tension_params, tdata, tdata->group_tension_params);
1779 old_handler = handler;
1780 old_inverse_handler = inverse_handler;
1784 /* free old groups */
1785 if (*pPer_group_handlers != NULL)
1786 free(*pPer_group_handlers);
1787 if (*pPer_group_inverse_handlers != NULL)
1788 free(*pPer_group_inverse_handlers);
1789 if (*pPer_group_data != NULL) {
1790 for (i=0; i<*pNum_active_groups; i++)
1791 free((*pPer_group_data)[i]);
1792 free(*pPer_group_data);
1795 /* set new groups */
1797 fprintf(stderr, "%s:\t%d groups, tension handler array %p, inverse tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_inverse_handlers, per_group_data, per_group_data[0]);
1799 *pNum_active_groups = num_active_groups;
1800 *pPer_group_handlers = per_group_handlers;
1801 *pPer_group_inverse_handlers = per_group_inverse_handlers;
1802 *pPer_group_data = per_group_data;
1806 <<crosslink groups>>=
1807 void crosslink_groups(list_t *state_groups, list_t *transition_list)
1809 list_t *list, *out_trans=NULL, *associates=NULL;
1810 one_dim_fn_t *handler;
1813 state_groups = head(state_groups);
1814 while (state_groups != NULL) {
1815 /* find transitions leaving this state */
1817 fprintf(stderr, "finding transitions leaving %s\n", S(state_groups)->name);
1819 list = head(transition_list);
1820 while (list != NULL) {
1821 if (T(list)->initial_state == S(state_groups) && list_index(associates, T(list)->initial_state) != -1) {
1822 push(&out_trans, T(list), 1);
1826 n = list_length(out_trans);
1827 S(state_groups)->num_out_trans = n;
1828 S(state_groups)->out_trans = (transition_t **)malloc(sizeof(transition_t *)*n);
1829 assert(S(state_groups)->out_trans != NULL);
1830 for (i=0; i<n; i++) {
1831 S(state_groups)->out_trans[i] = (transition_t *)pop(&out_trans);
1834 /* find states grouped with this state */
1836 fprintf(stderr, "finding states grouped with %s\n", S(state_groups)->name);
1838 handler = S(state_groups)->tension_handler;
1839 group = S(state_groups)->tension_group;
1840 list = head(state_groups);
1841 while (list != NULL) {
1842 if (S(list)->tension_handler == handler && S(list)->tension_group == group && list != state_groups) {
1843 push(&associates, S(list), 1);
1847 n = list_length(associates);
1848 S(state_groups)->num_grouped_states = n;
1849 S(state_groups)->grouped_states = (state_t **)malloc(sizeof(state_t *)*n);
1850 assert(S(state_groups)->grouped_states != NULL);
1851 for (i=0; i<n; i++) {
1852 S(state_groups)->grouped_states[i] = (state_t *)pop(&associates);
1854 state_groups = state_groups->next;
1860 \section{String parsing}
1862 For handling command line arguments, we need parse delimited strings (e.g.\ [["param1,param2,param3..."]]).
1863 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1864 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1865 need to take care of parsing those parameters themselves.
1866 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1872 <<parse definitions>>
1873 <<parse declarations>>
1877 <<parse module makefile lines>>=
1878 NW_SAWSIM_MODS += parse
1879 CHECK_BINS += check_parse
1883 <<parse definitions>>=
1884 #define SEP ',' /* argument separator character */
1887 <<parse declarations>>=
1888 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1889 int *num, char ***string_array);
1890 extern void free_parsed_list(int num, char **string_array);
1893 [[parse_list_string]] allocates memory, don't forget to free it
1894 afterward with [[free_parsed_list]]. It does not alter the original.
1896 The string may start off with a [[deeper]] character
1897 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1898 [[x,y]], where the pointer is one character in on the copied string.
1899 However, when we go to free the memory, we need a pointer to the
1900 beginning of the string. In order to accommodate this for a string
1901 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1902 the first $N$ elements point to the separated arguments, and let the
1903 last element point to the start of the copied string regardless of
1905 <<parse delimited list string functions>>=
1906 /* TODO, split out into parse.hc */
1907 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1910 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1911 if (string[i] == deeper) {depth++;}
1912 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1918 void parse_list_string(char *string, char sep, char deeper, char shallower,
1919 int *num, char ***string_array)
1921 char *str=NULL, **ret=NULL;
1923 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
1925 *string_array = NULL;
1928 /* make a copy of the string, so we don't change the original */
1929 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1930 assert(str != NULL);
1931 strcpy(str, string); /* we know str is long enough */
1932 /* count the number of regions, so we can allocate pointers to them */
1935 n++; i++; /* move on to next argument */
1936 i += next_delim_index(str+i, sep, deeper, shallower);
1937 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1938 } while (str[i] != '\0');
1939 ret = (char **)malloc(sizeof(char *)*(n+1));
1940 assert(ret != NULL);
1941 /* replace the separators with '\0' & assign pointers */
1942 ret[n] = str; /* point to the front of the copied string */
1945 for(i=1; i<n; i++) {
1946 j += next_delim_index(str+j, sep, deeper, shallower);
1948 ret[i] = str+j; /* point to the separated arguments */
1950 /* strip off bounding braces, if any */
1951 for(i=0; i<n; i++) {
1952 if (ret[i][0] == deeper) {
1956 if (ret[i][strlen(ret[i])-1] == shallower)
1957 ret[i][strlen(ret[i])-1] = '\0';
1960 *string_array = ret;
1963 void free_parsed_list(int num, char **string_array)
1966 /* frees all items, since they were allocated as a single string*/
1967 free(string_array[num]);
1968 /* and free the array of pointers */
1974 \subsection{Parsing implementation}
1980 <<parse delimited list string functions>>
1984 #include <assert.h> /* assert() */
1985 #include <stdlib.h> /* NULL */
1986 #include <stdio.h> /* fprintf(), stdout *//*!!*/
1987 #include <string.h> /* strlen() */
1991 \subsection{Parsing unit tests}
1993 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1995 <<parse check includes>>
1996 <<string comparison definition and globals>>
1997 <<check relative error>>
1998 <<parse test suite>>
1999 <<main check program>>
2002 <<parse check includes>>=
2003 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2004 #include <stdio.h> /* printf() */
2005 #include <assert.h> /* assert() */
2006 #include <string.h> /* strlen() */
2011 <<parse test suite>>=
2012 <<parse list string tests>>
2013 <<parse suite definition>>
2016 <<parse suite definition>>=
2017 Suite *test_suite (void)
2019 Suite *s = suite_create ("k model");
2020 <<parse list string test case defs>>
2022 <<parse list string test case add>>
2027 <<parse list string tests>>=
2030 START_TEST(test_next_delim_index)
2032 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
2033 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
2034 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
2035 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
2036 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
2040 START_TEST(test_parse_list_null)
2044 parse_list_string(NULL, SEP, '{', '}',
2045 &num_param_args, ¶m_args);
2046 fail_unless(num_param_args == 0, NULL);
2047 fail_unless(param_args == NULL, NULL);
2050 START_TEST(test_parse_list_single_simple)
2054 parse_list_string("arg", SEP, '{', '}',
2055 &num_param_args, ¶m_args);
2056 fail_unless(num_param_args == 1, NULL);
2057 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
2060 START_TEST(test_parse_list_single_compound)
2064 parse_list_string("{x,y,z}", SEP, '{', '}',
2065 &num_param_args, ¶m_args);
2066 fail_unless(num_param_args == 1, NULL);
2067 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
2070 START_TEST(test_parse_list_double_simple)
2074 parse_list_string("abc,def", SEP, '{', '}',
2075 &num_param_args, ¶m_args);
2076 fail_unless(num_param_args == 2, NULL);
2077 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2078 fail_unless(STRMATCH(param_args[1],"def"), NULL);
2081 START_TEST(test_parse_list_compound)
2085 parse_list_string("abc,{def,ghi}", SEP, '{', '}',
2086 &num_param_args, ¶m_args);
2087 fail_unless(num_param_args == 2, NULL);
2088 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2089 fail_unless(STRMATCH(param_args[1],"def,ghi"), NULL);
2093 <<parse list string test case defs>>=
2094 TCase *tc_parse_list_string = tcase_create("parse list string");
2096 <<parse list string test case add>>=
2097 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
2098 tcase_add_test(tc_parse_list_string, test_parse_list_null);
2099 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
2100 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
2101 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
2102 tcase_add_test(tc_parse_list_string, test_parse_list_compound);
2103 suite_add_tcase(s, tc_parse_list_string);
2107 \section{Unit tests}
2109 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2116 <<check relative error>>
2119 <<main check program>>
2131 <<determine dt tests>>
2133 <<does domain transition tests>>
2134 <<randomly unfold domains tests>>
2135 <<suite definition>>
2138 <<suite definition>>=
2139 Suite *test_suite (void)
2141 Suite *s = suite_create ("sawsim");
2142 <<F test case defs>>
2143 <<determine dt test case defs>>
2144 <<happens test case defs>>
2145 <<does domain transition test case defs>>
2146 <<randomly unfold domains test case defs>>
2149 <<determine dt test case add>>
2150 <<happens test case add>>
2151 <<does domain transition test case add>>
2152 <<randomly unfold domains test case add>>
2155 tcase_add_checked_fixture(tc_strip_address,
2156 setup_strip_address,
2157 teardown_strip_address);
2163 <<main check program>>=
2167 Suite *s = test_suite();
2168 SRunner *sr = srunner_create(s);
2169 srunner_run_all(sr, CK_ENV);
2170 number_failed = srunner_ntests_failed(sr);
2172 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
2176 \subsection{F tests}
2178 <<worm-like chain tests>>
2180 <<F test case defs>>=
2181 <<worm-like chain test case def>>
2183 <<F test case add>>=
2184 <<worm-like chain test case add>>
2187 <<worm-like chain tests>>=
2188 <<worm-like chain function>>
2190 START_TEST(test_wlc_at_zero)
2192 double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
2193 fail_unless(abs(wlc(x, T, p, L)) < lim, \
2194 "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
2195 x, T, p, L, abs(wlc(x, T, p, L)), lim);
2199 START_TEST(test_wlc_at_half)
2201 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
2202 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
2203 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
2205 * linear term = x/L = 0.5
2206 * nonlinear + linear = 0.75 + 0.5 = 1.25
2207 * wlc = 10e21*1.25 = 12.5e21
2209 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
2210 "wlc(%g, %g, %g, %g) = %g != %g",
2211 x, T, p, L, wlc(x, T, p, L), 12.5e21);
2215 <<worm-like chain test case def>>=
2216 TCase *tc_wlc = tcase_create("WLC");
2219 <<worm-like chain test case add>>=
2220 tcase_add_test(tc_wlc, test_wlc_at_zero);
2221 tcase_add_test(tc_wlc, test_wlc_at_half);
2222 suite_add_tcase(s, tc_wlc);
2225 \subsection{Model tests}
2227 Check the searching with [[linear_k]].
2228 Check overwhelming force treatment with the heavyside-esque step [[?]].
2230 Hmm, I'm not really sure what I was doing here...
2231 <<determine dt tests>>=
2232 double linear_k(double F, environment_t *env, void *params)
2234 double Fnot = *(double *)params;
2238 START_TEST(test_determine_dt_linear_k)
2241 transition_t unfold={0};
2242 environment_t env = {0};
2243 double F[]={0,1,2,3,4,5,6};
2244 double dt_max=3.0, Fnot=3.0, x=1.0, targ_p=0.1, v=1.0;
2245 list_t *state_list=NULL, *transition_list=NULL;
2248 folded.tension_handler = &hooke_handler;
2249 folded.num_domains = 1;
2250 unfold.initial_state = &folded;
2251 unfold.k = linear_k;
2252 unfold.k_params = &Fnot;
2253 push(&state_list, &folded, 1);
2254 push(&transition_list, &unfold, 1);
2255 for( i=0; i < sizeof(F)/sizeof(double); i++) {
2256 //fail_unless(determine_dt(state_list, transition_list, &env, x, targ_p, dt_max, v)==1, NULL);
2261 <<determine dt test case defs>>=
2262 TCase *tc_determine_dt = tcase_create("Determine dt");
2264 <<determine dt test case add>>=
2265 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
2266 suite_add_tcase(s, tc_determine_dt);
2271 <<happens test case defs>>=
2273 <<happens test case add>>=
2276 <<does domain transition tests>>=
2278 <<does domain transition test case defs>>=
2280 <<does domain transition test case add>>=
2283 <<randomly unfold domains tests>>=
2285 <<randomly unfold domains test case defs>>=
2287 <<randomly unfold domains test case add>>=
2291 \section{Balancing group extensions}
2292 \label{sec.tension_balance}
2294 For a given total extension $x$ (of the piezo), the various domain
2295 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
2296 amounts, and we need to tweak the portion that each extends to balance
2297 the tension amongst the active groups. Since the tension balancing is
2298 separable from the bulk of the simulation, we break it out into a
2299 separate module. The interface is defined in [[tension_balance.h]],
2300 the implementation in [[tension_balance.c]], and the unit testing in
2301 [[check_tension_balance.c]]
2303 <<tension-balance.h>>=
2304 #ifndef TENSION_BALANCE_H
2305 #define TENSION_BALANCE_H
2307 <<tension balance function declaration>>
2308 #endif /* TENSION_BALANCE_H */
2311 <<tension balance functions>>=
2312 <<one dimensional bracket>>
2313 <<one dimensional solve>>
2315 <<group stiffness function>>
2316 <<chain stiffness function>>
2317 <<tension balance function>>
2320 <<tension balance module makefile lines>>=
2321 NW_SAWSIM_MODS += tension_balance
2322 CHECK_BINS += check_tension_balance
2323 check_tension_balance_MODS =
2326 The entire force balancing problem reduces to a solving two nested
2327 one-dimensional root-finding problems. First we define the one
2328 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
2329 populated group). $x(x_0)$ is also strictly monotonically increasing,
2330 so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
2331 stores the last successful value of $x$, since we expect to be taking
2332 small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
2333 indicates that the groups have changed.
2334 <<tension balance function declaration>>=
2335 double tension_balance(int num_tension_groups,
2336 one_dim_fn_t **tension_handlers,
2337 one_dim_fn_t **inverse_tension_handlers,
2338 void **params, /* array of pointers to tension_handler_data_t */
2343 <<tension balance function>>=
2344 double tension_balance(int num_tension_groups,
2345 one_dim_fn_t **tension_handlers,
2346 one_dim_fn_t **inverse_tension_handlers,
2347 void **params, /* array of pointers to tension_handler_data_t */
2352 static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL,NULL};
2353 double F, xo_guess, xo, lb, ub;
2354 double min_relx=1e-6, min_rely=1e-6;
2355 int max_steps=200, i;
2357 assert(num_tension_groups > 0);
2358 assert(tension_handlers != NULL);
2359 assert(params != NULL);
2360 assert(num_tension_groups > 0);
2362 if (last_x == -1 || x_xo_data.xi == NULL) { /* new group setup, reset x_xo_data */
2363 /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
2364 if (x_xo_data.xi != NULL)
2366 /* construct new x_xo_data */
2367 x_xo_data.n_groups = num_tension_groups;
2368 x_xo_data.tension_handler = tension_handlers;
2369 x_xo_data.inverse_tension_handler = inverse_tension_handlers;
2370 x_xo_data.handler_data = params;
2372 fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p), x_of_xo_data %p\n", __FUNCTION__, x_xo_data.n_groups, x_xo_data.tension_handler, x_xo_data.handler_data, x_xo_data.handler_data[0], &x_xo_data);
2374 x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
2375 assert(x_xo_data.xi != NULL);
2376 for (i=0; i<num_tension_groups; i++)
2377 x_xo_data.xi[i] = last_x;
2380 if (num_tension_groups == 1) { /* only one group, no balancing required */
2382 x_xo_data.xi[0] = xo;
2384 fprintf(stderr, "%s:\tonly one active group, no balancing required\n", __FUNCTION__);
2388 fprintf(stderr, "%s:\tbalancing for x = %g with", __FUNCTION__, x);
2389 for(i=0; i<num_tension_groups; i++) fprintf(stderr, "\tx[%d]=%g", i, x_xo_data.xi[i]);
2390 fprintf(stderr, "\n");
2392 if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
2393 if (x == 0) xo_guess = length_scale;
2394 else xo_guess = x/num_tension_groups;
2396 fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
2398 oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
2399 } else { /* work off the last balanced point */
2401 fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
2404 lb = x_xo_data.xi[0];
2405 ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */
2406 } else if (x < last_x) {
2407 lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
2408 ub = x_xo_data.xi[0];
2409 } else { /* x == last_x */
2410 fprintf(stderr, "not moving\n");
2411 lb= x_xo_data.xi[0];
2412 ub= x_xo_data.xi[0];
2416 fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
2417 __FUNCTION__, x, lb, ub);
2419 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
2420 F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
2422 if (num_tension_groups > 1) {
2423 fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
2424 for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
2425 fprintf(stderr, "\n");
2430 F = (*tension_handlers[0])(xo, params[0]);
2432 if (stiffness != NULL)
2433 *stiffness = chain_stiffness(&x_xo_data, min_relx);
2440 <<tension balance internal definitions>>=
2441 <<x of x0 definitions>>
2444 <<x of x0 definitions>>=
2445 typedef struct x_of_xo_data_struct {
2447 one_dim_fn_t **tension_handler; /* array of fn pointers */
2448 one_dim_fn_t **inverse_tension_handler; /* inverse fn pointers */
2449 void **handler_data; /* array of pointers to tension_handler_data_t */
2450 double *xi; /* array of group extensions */
2455 double x_of_xo(double xo, void *pdata)
2457 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
2458 double F, x=0, xi, xi_guess, lb, ub, slope;
2460 double min_relx=1e-6, min_rely=1e-6;
2462 assert(data->n_groups > 0);
2465 fprintf(stderr, "%s:\tget initial tension at %g with %d groups, handler %p, data %p (x_xo_data %p)\n", __FUNCTION__, xo, data->n_groups, data->tension_handler[0], data->handler_data[0], data);
2468 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
2470 fprintf(stderr, "%s:\tfinding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
2474 if (data->xi) data->xi[0] = xo;
2475 for (i=1; i < data->n_groups; i++) {
2476 if (data->inverse_tension_handler[i] != NULL) {
2477 xi = (*data->inverse_tension_handler[i])(F, data->handler_data[i]);
2478 } else { /* invert numerically */
2479 if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
2480 else xi_guess = data->xi[i];
2482 fprintf(stderr, "%s:\tsearching for proper x[%d] for F[%d](x[%d]) = %g N (handler %p, data %p)\n", __FUNCTION__, i, i, i, F, data->tension_handler[i], data->handler_data[i]);
2484 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
2485 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
2490 fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
2494 fprintf(stderr, "%s:\tfound x(x0=%g) = %g\n", __FUNCTION__, xo, x);
2500 <<group stiffness function>>=
2501 double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double x, double relx)
2503 double dx, xl, F, Fl, stiffness;
2505 assert(relx > 0 && relx < 1);
2506 if (x == 0 || relx == 0) {
2507 dx = 1e-12; /* HACK, how to get length scale? */
2517 F = (*tension_handler)(x, handler_data);
2518 Fl = (*tension_handler)(xl, handler_data);
2519 stiffness = (F-Fl)/dx;
2520 assert(stiffness > 0);
2525 <<chain stiffness function>>=
2526 double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx)
2530 /* add group stiffnesses in series */
2531 for (i=0; i < x_xo_data->n_groups; i++) {
2533 fprintf(stderr, "%s:\tgetting stiffness of active state %d of %d for x[%d]=%g, relx=%g\n", __FUNCTION__, i, x_xo_data->n_groups, i, x_xo_data->xi[i], relx);
2536 stiffness += 1.0/group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx);
2538 stiffness = 1.0 / stiffness;
2544 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
2545 number of steps for monotonically increasing $f(x)$. Simple
2546 bisection, so it's robust and converges fairly quickly. You can set
2547 the maximum number of steps to take, but if you have enough processor
2548 power, it's better to limit your precision with the relative limits
2549 [[min_relx]] and [[y]]. These ensure that the width of the bracket is
2550 small on the length scale of it's larger side. Note that \emph{both}
2551 relative limits must be satisfied for the search to stop.
2552 <<one dimensional function definition>>=
2553 /* equivalent to gsl_func_t */
2554 typedef double one_dim_fn_t(double x, void *params);
2556 <<one dimensional solve declaration>>=
2557 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2558 double min_relx, double min_rely, int max_steps,
2562 <<one dimensional solve>>=
2563 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2564 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2565 double min_relx, double min_rely, int max_steps,
2568 double yx, ylb, yub, x;
2571 ylb = (*f)(lb, params);
2572 yub = (*f)(ub, params);
2574 /* check some simple cases */
2576 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bracketed */
2577 else return lb; /* any x on the range [lb, ub] would work */
2579 if (ylb == y) { x=lb; yx=ylb; goto end; }
2580 if (yub == y) { x=ub; yx=yub; goto end; }
2583 fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2585 assert(ylb < y); assert(yub > y);
2586 for (i=0; i<max_steps; i++) {
2588 yx = (*f)(x, params);
2590 fprintf(stderr, "%s:\tbracket: lb %g, x %g, ub %g\tylb %g, yx %g, yub %g\t(y %g)\n", __FUNCTION__, lb, x, ub, ylb, yx, yub, y);
2592 if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2593 else if (yx > y) { ub=x; yub=yx; }
2594 else /* yx < y */ { lb=x; ylb=yx; }
2599 fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2605 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2606 Generate bracketing $x$ values through bisection or exponential growth.
2607 <<one dimensional bracket declaration>>=
2608 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2611 <<one dimensional bracket>>=
2612 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2614 double yx, step, x=xguess;
2617 fprintf(stderr, "%s:\tbracketing %g (params %p)\n", __FUNCTION__, y, params);
2619 yx = (*f)(x, params);
2621 fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2628 fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2632 if (x == 0) x = length_scale/2.0; /* HACK */
2635 yx = (*f)(x, params);
2637 fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2642 fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2645 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2649 \subsection{Balancing implementation}
2651 <<tension-balance.c>>=
2654 <<tension balance includes>>
2655 #include "tension_balance.h"
2657 double length_scale = 1e-10; /* HACK */
2659 <<tension balance internal definitions>>
2660 <<tension balance functions>>
2663 <<tension balance includes>>=
2664 #include <assert.h> /* assert() */
2665 #include <stdlib.h> /* NULL */
2666 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2667 #include <stdio.h> /* fprintf(), stdout */
2671 \subsection{Balancing unit tests}
2673 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2674 <<check-tension-balance.c>>=
2675 <<tension balance check includes>>
2676 <<tension balance test suite>>
2677 <<main check program>>
2680 <<tension balance check includes>>=
2681 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2682 #include <assert.h> /* assert() */
2685 #include "tension_balance.h"
2688 <<tension balance test suite>>=
2689 <<tension balance function tests>>
2690 <<tension balance suite definition>>
2693 <<tension balance suite definition>>=
2694 Suite *test_suite (void)
2696 Suite *s = suite_create ("tension balance");
2697 <<tension balance function test case defs>>
2699 <<tension balance function test case adds>>
2704 <<tension balance function tests>>=
2705 <<check relative error>>
2707 double hooke(double x, void *pK)
2710 return *((double*)pK) * x;
2713 START_TEST(test_single_function)
2715 double k=5, x=3, last_x=2.0, F;
2716 one_dim_fn_t *handlers[] = {&hooke};
2717 one_dim_fn_t *inverse_handlers[] = {NULL};
2718 void *data[] = {&k};
2719 F = tension_balance(1, handlers, inverse_handlers, data, last_x, x, NULL);
2720 fail_unless(F = k*x, NULL);
2725 We can also test balancing two springs with different spring constants.
2726 The analytic solution for a general number of springs is given in Appendix \ref{sec.math_hooke}.
2727 <<tension balance function tests>>=
2728 START_TEST(test_double_hooke)
2730 double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
2731 one_dim_fn_t *handlers[] = {&hooke, &hooke};
2732 one_dim_fn_t *inverse_handlers[] = {NULL, NULL};
2733 void *data[] = {&k1, &k2};
2734 F = tension_balance(2, handlers, inverse_handlers, data, last_x, x, NULL);
2738 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2739 //CHECK_ERR(1e-6, x1e, xi[0]);
2740 //CHECK_ERR(1e-6, x2e, xi[1]);
2741 CHECK_ERR(1e-6, Fe, F);
2745 TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
2747 <<tension balance function test case defs>>=
2748 TCase *tc_tbfunc = tcase_create("tension balance function");
2751 <<tension balance function test case adds>>=
2752 tcase_add_test(tc_tbfunc, test_single_function);
2753 tcase_add_test(tc_tbfunc, test_double_hooke);
2754 suite_add_tcase(s, tc_tbfunc);
2760 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2761 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2762 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2768 <<list definitions>>
2769 <<list declarations>>
2773 <<list declarations>>=
2774 <<head and tail declarations>>
2775 <<list length declaration>>
2776 <<push declaration>>
2778 <<list destroy declaration>>
2779 <<list to array declaration>>
2780 <<move declaration>>
2781 <<sort declaration>>
2782 <<list index declaration>>
2783 <<list element declaration>>
2784 <<unique declaration>>
2785 <<list print declaration>>
2789 <<create and destroy node>>
2804 <<list module makefile lines>>=
2805 NW_SAWSIM_MODS += list
2806 CHECK_BINS += check_list
2810 <<list definitions>>=
2811 typedef struct list_struct {
2812 struct list_struct *next;
2813 struct list_struct *prev;
2817 <<comparison function typedef>>
2820 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2821 <<head and tail declarations>>=
2822 list_t *head(list_t *list);
2823 list_t *tail(list_t *list);
2826 list_t *head(list_t *list)
2830 while (list->prev != NULL) {
2836 list_t *tail(list_t *list)
2840 while (list->next != NULL) {
2847 <<list length declaration>>=
2848 int list_length(list_t *list);
2851 int list_length(list_t *list)
2858 while (list->next != NULL) {
2866 [[push]] inserts a node relative to the current position in the list
2867 without changing the current position.
2868 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2869 so we set it to that of the pushed domain.
2870 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2871 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2872 <<push declaration>>=
2873 void push(list_t **pList, void *data, int below);
2876 void push(list_t **pList, void *data, int below)
2878 list_t *list, *new_node;
2879 assert(pList != NULL);
2881 new_node = create_node(data);
2884 else if (below==0) { /* insert above */
2885 if (list->prev != NULL)
2886 list->prev->next = new_node;
2887 new_node->prev = list->prev;
2888 list->prev = new_node;
2889 new_node->next = list;
2890 } else { /* insert below */
2891 if (list->next != NULL)
2892 list->next->prev = new_node;
2893 new_node->next = list->next;
2894 list->next = new_node;
2895 new_node->prev = list;
2900 [[pop]] removes the current domain node, moving the current position
2901 to the node after it, unless that node is [[NULL]], in which case move
2902 the current position to the node before it.
2903 <<pop declaration>>=
2904 void *pop(list_t **pList);
2907 void *pop(list_t **pList)
2909 list_t *list, *popped;
2911 assert(pList != NULL);
2913 assert(list != NULL); /* not an empty list */
2915 /* bypass the popped node */
2916 if (list->prev != NULL)
2917 list->prev->next = list->next;
2918 if (list->next != NULL) {
2919 list->next->prev = list->prev;
2920 *pList = list->next;
2922 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
2924 destroy_node(popped);
2929 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
2930 If the cleanup function is [[NULL]], no cleanup function is called.
2931 The nodes are not popped in any particular order.
2932 <<list destroy declaration>>=
2933 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
2936 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
2940 assert(pList != NULL);
2943 assert(list != NULL); /* not an empty list */
2944 while (list != NULL) {
2946 if (destroy != NULL)
2952 [[list_to_array]] converts a list to an array of pointers, preserving order.
2953 <<list to array declaration>>=
2954 void list_to_array(list_t *list, int *length, void ***array);
2957 void list_to_array(list_t *list, int *pLength, void ***pArray)
2961 assert(list != NULL);
2962 assert(pLength != NULL);
2963 assert(pArray != NULL);
2965 length = list_length(list);
2966 array = (void **)malloc(sizeof(void *)*length);
2967 assert(array != NULL);
2970 while (list != NULL) {
2971 array[i++] = list->d;
2979 [[move]] moves the current element along the list in either direction.
2980 <<move declaration>>=
2981 void move(list_t **pList, int downstream);
2984 void move(list_t **pList, int downstream)
2986 list_t *list, *flipper;
2988 assert(pList != NULL);
2989 list = *pList; /* a>B>c>d */
2990 assert(list != NULL); /* not an empty list */
2991 if (downstream == 0)
2992 flipper = list->prev; /* flipper is A */
2994 flipper = list->next; /* flipper is C */
2995 /* ensure there is somewhere to go in the direction we're heading */
2996 assert(flipper != NULL);
2997 /* remove the current element */
2998 data = pop(&list); /* data=B, a>C>d */
2999 /* and put it back in in the correct spot */
3001 if (downstream == 0) {
3002 push(&list, data, 0); /* b>A>c>d */
3003 list = list->prev; /* B>a>c>d */
3005 push(&list, data, 1); /* a>C>b>d */
3006 list = list->next; /* a>c>B>d */
3012 [[sort]] sorts a list from smallest to largest according to a user
3013 specified comparison.
3014 <<comparison function typedef>>=
3015 typedef int (comparison_fn_t)(void *A, void *B);
3018 <<int comparison function>>=
3019 int int_comparison(void *A, void *B)
3024 if (a > b) return 1;
3025 else if (a == b) return 0;
3030 <<sort declaration>>=
3031 void sort(list_t **pList, comparison_fn_t *comp);
3033 Here we use the bubble sort algorithm.
3035 void sort(list_t **pList, comparison_fn_t *comp)
3039 assert(pList != NULL);
3041 assert(list != NULL); /* not an empty list */
3042 while (stable == 0) {
3045 while (list->next != NULL) {
3046 if ((*comp)(list->d, list->next->d) > 0) {
3048 move(&list, 1 /* downstream */);
3058 [[list_index]] finds the location of [[data]] in [[list]]. Returns
3059 [[-1]] if [[data]] is not in [[list]].
3060 <<list index declaration>>=
3061 int list_index(list_t *list, void *data);
3065 int list_index(list_t *list, void *data)
3069 while (list != NULL) {
3070 if (list->d == data) return i;
3079 [[list_element]] returns the data bound to the $i^\text{th}$ node of the list.
3080 <<list element declaration>>=
3081 void *list_element(list_t *list, int i);
3084 void *list_element(list_t *list, int i)
3088 while (list != NULL) {
3089 if (i == j) return list->d;
3099 [[unique]] removes duplicates from a sorted list according to a user
3100 specified comparison. Don't do this unless you have other pointers
3101 any dynamically allocated data in your list, because [[unique]] just
3102 drops any non-unique data without freeing it.
3103 <<unique declaration>>=
3104 void unique(list_t **pList, comparison_fn_t *comp);
3107 void unique(list_t **pList, comparison_fn_t *comp)
3110 assert(pList != NULL);
3112 assert(list != NULL); /* not an empty list */
3114 while (list->next != NULL) {
3115 if ((*comp)(list->d, list->next->d) == 0) {
3124 [[list_print]] prints a list to a [[FILE *]] stream.
3125 <<list print declaration>>=
3126 void list_print(FILE *file, list_t *list, char *list_name);
3129 void list_print(FILE *file, list_t *list, char *list_name)
3131 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
3133 while (list != NULL) {
3134 fprintf(file, " %p", list->d);
3137 fprintf(file, "\n");
3141 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
3142 <<create and destroy node>>=
3143 list_t *create_node(void *data)
3145 list_t *ret = (list_t *)malloc(sizeof(list_t));
3146 assert(ret != NULL);
3153 void destroy_node(list_t *node)
3159 The user must free the data pointed to by the node on their own.
3161 \subsection{List implementation}
3171 #include <assert.h> /* assert() */
3172 #include <stdlib.h> /* malloc(), free(), rand() */
3173 #include <stdio.h> /* fprintf(), stdout, FILE */
3174 #include "global.h" /* destroy_data_func_t */
3177 \subsection{List unit tests}
3179 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3181 <<list check includes>>
3183 <<main check program>>
3186 <<list check includes>>=
3187 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3188 #include <stdio.h> /* FILE */
3194 <<list test suite>>=
3197 <<list suite definition>>
3200 <<list suite definition>>=
3201 Suite *test_suite (void)
3203 Suite *s = suite_create ("list");
3204 <<push test case defs>>
3205 <<pop test case defs>>
3207 <<push test case adds>>
3208 <<pop test case adds>>
3214 START_TEST(test_push)
3219 push(&list, (void *)i, 1);
3220 fail_unless(list == head(list), NULL);
3221 fail_unless( (int)list->d == 0 );
3222 for (i=0; i<n; i++) {
3223 /* we expect 0, n-1, n-2, ..., 1 */
3226 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
3232 <<push test case defs>>=
3233 TCase *tc_push = tcase_create("push");
3236 <<push test case adds>>=
3237 tcase_add_test(tc_push, test_push);
3238 suite_add_tcase(s, tc_push);
3243 <<pop test case defs>>=
3245 <<pop test case adds>>=
3248 \section{Function string evaluation}
3250 For the saddle-point approximated Kramers' model (Section \ref{sec.kramers}) we need the ability to evaluate user-supplied functions ($E(x)$, $x_{ts}(F)$, \ldots).
3251 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
3252 The solution is to run a scripting language as a subprocess accessed via pipes.
3253 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
3255 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
3256 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
3257 Most of this is also POSIX-specific (as far as I know), so you'll have to do some legwork to port it to non-POSIX systems.
3258 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
3260 If you feel the benefits do \emph{not} outweigh the costs, or you are on a non-POSIX system (e.g.\ MS Windows without Cygwin), you should probably hardcode your functions in \lang.
3261 Then you can either statically or dynamically link to those hardcoded functions.
3262 While much less flexible, this approach would be a fairly simple-to-implement fallback.
3264 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
3265 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
3268 #ifndef STRING_EVAL_H
3269 #define STRING_EVAL_H
3271 <<string eval setup declaration>>
3272 <<string eval function declaration>>
3273 <<string eval teardown declaration>>
3274 #endif /* STRING_EVAL_H */
3277 <<string eval module makefile lines>>=
3278 NW_SAWSIM_MODS += string_eval
3279 CHECK_BINS += check_string_eval
3280 check_string_eval_MODS =
3283 For an introduction to POSIX process control, see\\
3284 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
3285 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
3286 and of course, the relavent [[man]] pages.
3288 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
3289 [[execvp]] replaces the calling process' program with a new program.
3290 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
3291 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
3292 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
3293 The new program needs command line arguments, just like it would if you were running it from a shell.
3294 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
3295 with the final array entry being a [[NULL]] pointer.
3297 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
3298 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
3299 <<string eval subprocess definitions>>=
3300 #define SUBPROCESS "python"
3301 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
3302 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
3303 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
3305 The [[i]] option lets Python know that it should run in interactive mode.
3306 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
3307 In interactive mode, python acts on every instruction as soon as it is recieved.
3308 The [[c]] option signals a command line instruction for Python to execute on startup, which in our case simply turns off the default prompting ([[ps1]] and [[ps2]]), since that is mostly for human users, and our program doesn't need it.
3309 %The [[c]] option signals a command line instruction for Python to execute on startup, which in our case simply [[pass]]es so that the silly Python header information is not printed. We leave the prompts in, because we scan for them to determine when the output has completed.
3311 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
3312 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
3313 [[fork]] returns two copies of the same program, executing the original code.
3314 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
3315 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
3317 We communicate with the child (Python) process using \emph{pipes}, with one process writing data into one end of the pipe, and the other process reading the data out of the other end.
3318 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
3319 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
3320 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
3321 <<string eval pipe definitions>>=
3322 #define PIPE_READ 0 /* the end you read from */
3323 #define PIPE_WRITE 1 /* the end you write to */
3325 #define STDIN 0 /* initial index of stdin pair */
3326 #define STDOUT 2 /* initial index of stdout pair */
3329 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
3331 As a finishing touch, we can promote the POSIX file descriptors ([[read]]/[[write]] interface) into the more familiar [[stdio.h]] \emph{streams} ([[fprintf]]/[[fgetc]] interface) using [[fdopen]], which creates a stream from an open file descriptor.
3333 <<string eval setup declaration>>=
3334 extern void string_eval_setup(FILE **pIn, FILE **pOut);
3336 <<string eval setup definition>>=
3337 void string_eval_setup(FILE **pIn, FILE **pOut)
3340 int pfd[4]; /* file descriptors for stdin and stdout */
3342 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
3343 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
3345 pid = fork(); /* split process into two copies */
3347 if (pid == -1) { /* fork error */
3348 perror("fork error");
3350 } else if (pid == 0) { /* child */
3351 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
3352 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
3353 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
3354 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
3355 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
3356 perror("exec error"); /* exec shouldn't return */
3358 } else { /* parent */
3359 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
3360 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
3361 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
3362 if ( *pIn == NULL ) {
3363 perror("fdopen (in)");
3366 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
3367 if ( *pOut == NULL ) {
3368 perror("fdopen (out)");
3375 To use the evaluating subprocess, we just pipe in our command, and read out the results.
3376 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
3377 <<string eval function declaration>>=
3378 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
3380 <<string eval function definition>>=
3381 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
3384 rc = fprintf(in, "%s", input);
3385 assert(rc == strlen(input));
3388 alarm(1); /* set a one second timeout on the read */
3389 assert( fgets(output, buflen, out) != NULL );
3390 alarm(0); /* cancel the timeout */
3391 //fprintf(stderr, "eval: %s ----> %s", input, output);
3394 The [[alarm]] calls set and clear a timeout on the returned output.
3395 If the timeout expires, the process would get a [[SIGALRM]], but it doesn't have a [[SIGALRM]] handler, so it gets a [[SIGKILL]] and dies.
3396 This protects against invalid input for which a line of output is not printed to [[stdout]].
3397 Other invalid inputs (e.g.\ those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
3398 If you are getting strange results, check your python code seperately. TODO, better error handling.
3400 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
3401 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
3402 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
3403 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
3404 <<string eval teardown declaration>>=
3405 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
3408 <<string eval teardown definition>>=
3409 void string_eval_teardown(FILE **pIn, FILE **pOut)
3414 /* redirect Python's stderr */
3415 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
3419 assert( fclose(*pIn) == 0 );
3421 assert( fclose(*pOut) == 0 );
3424 /* wait for python to exit */
3426 pid = wait(&stat_loc);
3433 if (WIFEXITED(stat_loc)) {
3434 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
3435 } else if (WIFSIGNALED(stat_loc)) {
3436 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
3441 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
3443 \subsection{String evaluation implementation}
3447 <<string eval includes>>
3448 #include "string_eval.h"
3449 <<string eval internal definitions>>
3450 <<string eval functions>>
3453 <<string eval includes>>=
3454 #include <assert.h> /* assert() */
3455 #include <stdlib.h> /* NULL */
3456 #include <stdio.h> /* fprintf(), stdout, fdopen() */
3457 #include <string.h> /* strlen() */
3458 #include <sys/types.h> /* pid_t */
3459 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
3460 #include <wait.h> /* wait() */
3463 <<string eval internal definitions>>=
3464 <<string eval subprocess definitions>>
3465 <<string eval pipe definitions>>
3468 <<string eval functions>>=
3469 <<string eval setup definition>>
3470 <<string eval function definition>>
3471 <<string eval teardown definition>>
3474 \subsection{String evaluation unit tests}
3476 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3477 <<check-string-eval.c>>=
3478 <<string eval check includes>>
3479 <<string comparison definition and globals>>
3480 <<string eval test suite>>
3481 <<main check program>>
3484 <<string eval check includes>>=
3485 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3486 #include <stdio.h> /* printf() */
3487 #include <assert.h> /* assert() */
3488 #include <string.h> /* strlen() */
3489 #include <signal.h> /* SIGKILL */
3491 #include "string_eval.h"
3494 <<string eval test suite>>=
3495 <<string eval tests>>
3496 <<string eval suite definition>>
3499 <<string eval suite definition>>=
3500 Suite *test_suite (void)
3502 Suite *s = suite_create ("string eval");
3503 <<string eval test case defs>>
3505 <<string eval test case add>>
3510 <<string eval tests>>=
3511 START_TEST(test_setup_teardown)
3514 string_eval_setup(&in, &out);
3515 string_eval_teardown(&in, &out);
3518 START_TEST(test_invalid_command)
3521 char input[80], output[80]={};
3522 string_eval_setup(&in, &out);
3523 sprintf(input, "print ABCDefg\n");
3524 string_eval(in, out, input, 80, output);
3525 string_eval_teardown(&in, &out);
3528 START_TEST(test_simple_eval)
3531 char input[80], output[80]={};
3532 string_eval_setup(&in, &out);
3533 sprintf(input, "print 3+4*5\n");
3534 string_eval(in, out, input, 80, output);
3535 fail_unless(STRMATCH(output,"23\n"), NULL);
3536 string_eval_teardown(&in, &out);
3539 START_TEST(test_multiple_evals)
3542 char input[80], output[80]={};
3543 string_eval_setup(&in, &out);
3544 sprintf(input, "print 3+4*5\n");
3545 string_eval(in, out, input, 80, output);
3546 fail_unless(STRMATCH(output,"23\n"), NULL);
3547 sprintf(input, "print (3**2 + 4**2)**0.5\n");
3548 string_eval(in, out, input, 80, output);
3549 fail_unless(STRMATCH(output,"5.0\n"), NULL);
3550 string_eval_teardown(&in, &out);
3553 START_TEST(test_eval_with_state)
3556 char input[80], output[80]={};
3557 string_eval_setup(&in, &out);
3558 sprintf(input, "print 3+4*5\n");
3559 fprintf(in, "A = 3\n");
3560 sprintf(input, "print A*3\n");
3561 string_eval(in, out, input, 80, output);
3562 fail_unless(STRMATCH(output,"9\n"), NULL);
3563 string_eval_teardown(&in, &out);
3567 <<string eval test case defs>>=
3568 TCase *tc_string_eval = tcase_create("string_eval");
3570 <<string eval test case add>>=
3571 tcase_add_test(tc_string_eval, test_setup_teardown);
3572 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
3573 tcase_add_test(tc_string_eval, test_simple_eval);
3574 tcase_add_test(tc_string_eval, test_multiple_evals);
3575 tcase_add_test(tc_string_eval, test_eval_with_state);
3576 suite_add_tcase(s, tc_string_eval);
3580 \section{Accelerating function evaluation}
3582 My first version-0.3 code was running very slowly.
3583 With the options suggested in the help
3584 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
3585 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
3586 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
3587 That's only 15 calls per solution, so the search algorithm seems reasonable.
3588 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
3593 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
3595 #endif /* ACCEL_K_H */
3598 <<accel k module makefile lines>>=
3599 NW_SAWSIM_MODS += accel_k
3600 #CHECK_BINS += check_accel_k
3601 check_accel_k_MODS =
3605 #include <assert.h> /* assert() */
3606 #include <stdlib.h> /* realloc(), free(), NULL */
3607 #include "global.h" /* environment_t */
3608 #include "k_model.h" /* k_func_t */
3609 #include "interp.h" /* interp_* */
3610 #include "accel_k.h"
3612 typedef struct accel_k_struct {
3613 interp_table_t *itable;
3619 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3620 static int num_accels = 0;
3621 static accel_k_t *accels=NULL;
3623 /* Wrap k in the standard f(x) acceleration form */
3624 static double k_wrap(double F, void *params)
3626 accel_k_t *p = (accel_k_t *)params;
3627 return (*p->k)(F, p->env, p->params);
3630 static int k_tol(double FA, double kA, double FB, double kB)
3633 if (FB-FA > 1e-12) {
3634 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3635 return 1; /* unnacceptable */
3637 //printf("acceptable tol\n");
3638 return 0; /* acceptable */
3642 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3645 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3646 assert(accels != NULL);
3647 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3649 accels[i].env = env;
3650 accels[i].params = params;
3657 for (i=0; i<num_accels; i++)
3658 interp_table_free(accels[i].itable);
3660 if (accels) free(accels);
3664 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3667 for (i=0; i<num_accels; i++) {
3668 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3669 /* We've already setup this function.
3670 * WARNING: we're assuming param and env are the same. */
3671 return interp_table_eval(accels[i].itable, accels+i, F);
3674 /* set up a new acceleration environment */
3675 i = add_accel_k(k, env, params);
3676 return interp_table_eval(accels[i].itable, accels+i, F);
3680 \section{Tension models}
3681 \label{sec.tension_models}
3683 TODO: tension model intro.
3684 See \url{http://web.mit.edu/course/3/3.11/www/pset03/Rec9.pdf} for a quick introduction to worm-like and freely-jointed chain models.
3686 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3687 The interface is defined in [[tension_model.h]], the implementation in [[tension_model.c]], the unit testing in [[check_k_model.c]], and some other tidbits in [[tension_model_utils.c]].
3689 <<tension-model.h>>=
3690 #ifndef TENSION_MODEL_H
3691 #define TENSION_MODEL_H
3693 <<tension handler types>>
3694 <<constant tension model declarations>>
3695 <<hooke tension model declarations>>
3696 <<worm-like chain tension model declarations>>
3697 <<freely-jointed chain tension model declarations>>
3698 <<find tension definitions>>
3699 <<tension model global declarations>>
3700 #endif /* TENSION_MODEL_H */
3703 <<tension model module makefile lines>>=
3704 NW_SAWSIM_MODS += tension_model
3705 #CHECK_BINS += check_tension_model
3706 check_tension_model_MODS =
3708 <<tension model utils makefile lines>>=
3709 TENSION_MODEL_MODS = tension_model parse list tension_balance
3710 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3711 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3712 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3713 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
3714 notangle -Rtension-model-utils.c $< > $@
3715 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
3716 gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3717 clean_tension_model_utils :
3718 rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3719 @ The pipe symbol [[|]] marks the prerequisits following it (in this
3720 case, the directories) as ‘order-only’ prerequisites. The timestamp
3721 on these prerequisits does not effect whether the rules are executed.
3722 This is appropriate for directories, where we don't need to recompile
3723 after an unrelated has been added to the directory, but only when the
3724 source prerequisites change. See the [[make]] documentation for more
3726 (\url{http://www.gnu.org/software/make/manual/html_node/Prerequisite-Types.html}).
3729 \label{sec.null_tension}
3731 For unstretchable domains.
3733 <<null tension model getopt>>=
3734 {NULL, NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3737 \subsection{Constant}
3738 \label{sec.const_tension}
3740 An infinitely stretchable domain providing a constant tension.
3741 Obviously this cannot be inverted, so you can't put this domain in
3742 series with any other domains. We include it mostly for testing
3745 <<constant tension functions>>=
3746 <<constant tension function>>
3747 <<constant tension parameter create/destroy functions>>
3750 <<constant tension model declarations>>=
3751 <<constant tension function declaration>>
3752 <<constant tension parameter create/destroy function declarations>>
3753 <<constant tension model global declarations>>
3757 <<constant tension function declaration>>=
3758 extern double const_tension_handler(double x, void *pdata);
3760 <<constant tension function>>=
3761 double const_tension_handler(double x, void *pdata)
3763 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3764 list_t *list = data->group_tension_params;
3769 assert(list != NULL); /* empty active group?! */
3770 F = ((const_tension_param_t *)list->d)->F;
3772 fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3774 while (list != NULL) {
3775 assert(((const_tension_param_t *)list->d)->F == F);
3783 <<constant tension structure>>=
3784 typedef struct const_tension_param_struct {
3785 double F; /* tension (force) in N */
3786 } const_tension_param_t;
3790 <<constant tension parameter create/destroy function declarations>>=
3791 extern void *string_create_const_tension_param_t(char **param_strings);
3792 extern void destroy_const_tension_param_t(void *p);
3795 <<constant tension parameter create/destroy functions>>=
3796 const_tension_param_t *create_const_tension_param_t(double F)
3798 const_tension_param_t *ret
3799 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3800 assert(ret != NULL);
3805 void *string_create_const_tension_param_t(char **param_strings)
3807 return create_const_tension_param_t(
3808 safe_strtod(param_strings[0],__FUNCTION__));
3811 void destroy_const_tension_param_t(void *p)
3819 <<constant tension model global declarations>>=
3820 extern int num_const_tension_params;
3821 extern char *const_tension_param_descriptions[];
3822 extern char const_tension_param_string[];
3825 <<constant tension model globals>>=
3826 int num_const_tension_params = 1;
3827 char *const_tension_param_descriptions[] = {"tension F, N"};
3828 char const_tension_param_string[] = "0";
3831 <<constant tension model getopt>>=
3832 {&const_tension_handler, NULL, "const", "an infinitely stretchable domain with constant tension", 1, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t}
3838 <<hooke functions>>=
3839 <<internal hooke functions>>
3841 <<inverse hooke handler>>
3842 <<hooke parameter create/destroy functions>>
3845 <<hooke tension model declarations>>=
3846 <<hooke tension function declaration>>
3847 <<hooke tension parameter create/destroy function declarations>>
3848 <<hooke tension model global declarations>>
3852 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
3853 The behavior of a series of springs $k_i$ in series is given by
3855 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
3857 For a simple proof, see Appendix \ref{sec.math_hooke}.
3859 <<hooke tension function declaration>>=
3860 extern double hooke_handler(double x, void *pdata);
3861 extern double inverse_hooke_handler(double F, void *pdata);
3864 First we define a function that computes the effective spring constant
3865 (stored in a single [[hooke_param_t]]) for the entire group.
3866 <<internal hooke functions>>=
3867 static void hooke_param_grouper(tension_handler_data_t *data,
3868 hooke_param_t *param) {
3869 list_t *list = data->group_tension_params;
3873 assert(list != NULL); /* empty active group?! */
3874 while (list != NULL) {
3875 assert( ((hooke_param_t *)list->d)->k > 0 );
3876 k += 1.0/ ((hooke_param_t *)list->d)->k;
3878 fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
3879 __FUNCTION__, ((hooke_param_t *)list->d)->k);
3885 fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N (from data %p)\n",
3886 __FUNCTION__, k, x, k*x, data);
3893 Everyone knows Hooke's law: $F=kx$.
3895 double hooke_handler(double x, void *pdata)
3897 hooke_param_t param = {0};
3900 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
3906 Inverting Hooke's law gives $x=F/k$.
3907 <<inverse hooke handler>>=
3908 double inverse_hooke_handler(double F, void *pdata)
3910 hooke_param_t param = {0};
3913 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
3919 Not much to keep track of for a single effective spring, just the
3920 spring constant $k$.
3921 <<hooke structure>>=
3922 typedef struct hooke_param_struct {
3923 double k; /* spring constant in N/m */
3928 We wrap up our Hooke implementation with some book-keeping functions.
3929 <<hooke tension parameter create/destroy function declarations>>=
3930 extern void *string_create_hooke_param_t(char **param_strings);
3931 extern void destroy_hooke_param_t(void *p);
3935 <<hooke parameter create/destroy functions>>=
3936 hooke_param_t *create_hooke_param_t(double k)
3938 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
3939 assert(ret != NULL);
3944 void *string_create_hooke_param_t(char **param_strings)
3946 return create_hooke_param_t(safe_strtod(param_strings[0], __FUNCTION__));
3949 void destroy_hooke_param_t(void *p)
3956 <<hooke tension model global declarations>>=
3957 extern int num_hooke_params;
3958 extern char *hooke_param_descriptions[];
3959 extern char hooke_param_string[];
3962 <<hooke tension model globals>>=
3963 int num_hooke_params = 1;
3964 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
3965 char hooke_param_string[]="0.05";
3968 <<hooke tension model getopt>>=
3969 {hooke_handler, inverse_hooke_handler, "hooke", "a Hookean spring (F=kx)", num_hooke_params, hooke_param_descriptions, (char *)hooke_param_string, &string_create_hooke_param_t, &destroy_hooke_param_t}
3972 \subsection{Worm-like chain}
3975 We can model several unfolded domains with a single worm-like chain.
3976 <<worm-like chain functions>>=
3977 <<internal worm-like chain functions>>
3978 <<worm-like chain handler>>
3979 <<inverse worm-like chain handler>>
3980 <<worm-like chain create/destroy functions>>
3983 <<worm-like chain tension model declarations>>=
3985 <<worm-like chain handler declaration>>
3986 <<inverse worm-like chain handler declaration>>
3987 <<worm-like chain create/destroy function declarations>>
3988 <<worm-like chain tension model global declarations>>
3992 <<internal worm-like chain functions>>=
3993 <<worm-like chain function>>
3994 <<inverse worm-like chain function>>
3995 <<worm-like chain parameter grouper>>
3998 The combination of all unfolded domains is modeled as a worm like chain,
3999 with the force $F_{\text{WLC}}$ approximately given by
4001 F_{\text{WLC}} \approx \frac{k_B T}{p}
4002 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
4004 where $x$ is the distance between the fixed ends,
4005 $k_B$ is Boltzmann's constant,
4006 $T$ is the absolute temperature,
4007 $p$ is the persistence length, and
4008 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
4009 <<worm-like chain function>>=
4010 static double wlc(double x, double T, double p, double L)
4012 double xL=0.0; /* uses global k_B */
4014 assert(T > 0); assert(p > 0); assert(L > 0);
4015 if (x >= L) return HUGE_VAL;
4016 xL = x/L; /* unitless */
4017 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
4018 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
4019 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
4024 Inverting the approximate WLC is fairly straightforward, if a bit tedious.
4025 Letting $x_L \equiv x/L$ and $F_T \equiv Fp/k_B T$ we have
4027 F_T &\approx \frac{1}{4} \p({\frac{1}{(1-x_L)^2} - 1}) + x_L \\
4028 F_T (1-x_L)^2 &\approx \frac{1}{4} - \frac{1}{4} (1-x_L)^2 + x_L(1-x_L)^2 \\
4029 0 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) (1-x_L)^2 + x_L(1-x_L)^2 \\
4030 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) \p({1-2x_L+x_L^2})
4031 + x_L - 2x_L^2 + x_L^3 \\
4032 &\approx \p[{\frac{1}{4} - \p({F_T + \frac{1}{4}})}]
4033 + x_L \p[{2F_T + \frac{1}{2} + 1}]
4034 + x_L^2 \p[{-\p({F_T + \frac{1}{4}}) - 2}]
4037 + x_L \p({2F_T + \frac{3}{2}})
4038 - x_L^2 \p({F_T + \frac{9}{4}})
4040 &\approx ax_L^3 + bx_L^2 + cx_L + d
4042 where $a \equiv 1$, $b \equiv -F_T - 9/4$, $c \equiv 2F_T + 3/2$, and
4044 % From \citet{wikipedia_cubic_function} the discriminant
4046 % \Delta &= -4b^3d + b^2c^2 - 4ac^3 + 18abcd - 27a^2d^2 \\
4047 % &= -4F_T\p({F_T + \frac{9}{4}})^3
4048 % + \p({F_T + \frac{9}{4}})^2 \cdot \p({2F_T + \frac{3}{2}})^2
4049 % - 4 \p({2F_T + \frac{3}{2}})^3
4050 % + 18F_T \p({F_T + \frac{9}{4}}) \cdot \p({2F_T + \frac{3}{2}})
4052 % &= -4F_T \p({F_T^3 + \frac{27}{4}F_T^2 +\frac{243}{16}F_T +\frac{729}{64}})
4053 %% a^3 + 3a^2b + 3ab^2 + b^3
4054 % + \p({F_T^2 + \frac{9}{2}F_T + \frac{81}{16}})
4055 % \cdot \p({4F_T^2 + 6F_T + \frac{9}{4}}) \\
4056 % &\qquad - 4 \p({8F_T^3 + 18F_T^2 + \frac{27}{2}F_T + \frac{27}{8}})
4057 %% a^3 + 3a^2b + 3ab^2 + b^3
4058 % + 18F_T \p({2F_T^2 + 6F_T + \frac{27}{8}})
4060 % &= -4F_T^4 - 27F_T^3 - \frac{243}{4}F_T^2 - \frac{729}{16}F_T
4061 % + 4F_T^4 + 24F_T^3 +\frac{99}{2}F_T^2 +\frac{81}{2}F_T +\frac{729}{64}\\
4062 % &\qquad - 32F_T^3 - 72F_T^2 - 54F_T - \frac{27}{2}
4063 % + 36F_T^3 + 108F_T^2 + \frac{243}{4}F_T
4065 % &= F_T^4 \p({-4 + 4}) + F_T^3 \p({-27 + 24 - 32 + 36})
4066 % + F_T^2 \p({-\frac{243}{4} + \frac{99}{2} - 72 + 108 -27}) \\
4067 % &\qquad + F_T \p({-\frac{729}{16} + \frac{81}{2} - 54 + \frac{243}{4}})
4068 % \p({\frac{729}{64} - \frac{27}{2}}) \\
4069 % &= F_T^3 - \frac{9}{4}F_T^2 + \frac{27}{16}F_T - \frac{135}{64}
4071 We can plug these coefficients into the GSL cubic solver to invert the
4072 WLC\citep{galassi05}. The applicable root is always the one which
4073 comes before the singularity, which will be the smallest real root.
4074 <<inverse worm-like chain function>>=
4075 static double inverse_wlc(double F, double T, double p, double L)
4077 double FT = F*p/(k_B*T); /* uses global k_B */
4078 double xL0, xL1, xL2;
4081 assert(T > 0); assert(p > 0); assert(L > 0);
4084 num_roots = gsl_poly_solve_cubic(-(FT+2.25),2*FT+1.5,-FT, &xL0, &xL1, &xL2);
4085 assert(num_roots > 0); assert(xL0 >= -DOUBLE_PRECISION); assert(xL0 < 1);
4086 if (xL0 < 0) xL0 = 0.0;
4092 First we define a function that computes the effective WLC parameters
4093 (stored in a single [[wlc_param_t]]) for the entire group.
4094 <<worm-like chain parameter grouper>>=
4095 static void wlc_param_grouper(tension_handler_data_t *data,
4096 wlc_param_t *param) {
4097 list_t *list = data->group_tension_params;
4098 double p=0.0, L=0.0;
4101 assert(list != NULL); /* empty active group?! */
4102 p = ((wlc_param_t *)list->d)->p;
4103 while (list != NULL) {
4104 /* enforce consistent p */
4105 assert( ((wlc_param_t *)list->d)->p == p);
4106 L += ((wlc_param_t *)list->d)->L;
4108 fprintf(stderr, "%s: WLC section %g m long in series\n",
4109 __FUNCTION__, ((wlc_param_t *)list->d)->L);
4114 fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
4115 __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
4123 <<worm-like chain handler declaration>>=
4124 extern double wlc_handler(double x, void *pdata);
4127 This model requires that each unfolded domain in the group have the
4128 same persistence length.
4129 <<worm-like chain handler>>=
4130 double wlc_handler(double x, void *pdata)
4132 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4133 wlc_param_t param = {0};
4135 wlc_param_grouper(data, ¶m);
4136 return wlc(x, data->env->T, param.p, param.L);
4141 <<inverse worm-like chain handler declaration>>=
4142 extern double inverse_wlc_handler(double F, void *pdata);
4145 <<inverse worm-like chain handler>>=
4146 double inverse_wlc_handler(double F, void *pdata)
4148 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4149 wlc_param_t param = {0};
4151 wlc_param_grouper(data, ¶m);
4152 return inverse_wlc(F, data->env->T, param.p, param.L);
4157 <<worm-like chain structure>>=
4158 typedef struct wlc_param_struct {
4159 double p; /* WLC persistence length (m) */
4160 double L; /* WLC contour length (m) */
4164 <<worm-like chain create/destroy function declarations>>=
4165 extern void *string_create_wlc_param_t(char **param_strings);
4166 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
4169 <<worm-like chain create/destroy functions>>=
4170 wlc_param_t *create_wlc_param_t(double p, double L)
4172 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
4173 assert(ret != NULL);
4179 void *string_create_wlc_param_t(char **param_strings)
4181 return create_wlc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4182 safe_strtod(param_strings[1], __FUNCTION__));
4185 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
4193 We define [[wlc_param_string]] as a global array because hard-coding the values into the tension model getopt structure gives a read-only version.
4194 TODO (now we copy the string before parsing, but still do this for clarity).
4196 <<valid string write code>>=
4197 char string[] = "abc";
4200 <<valid string write code>>=
4201 char *string = "abc";
4205 <<worm-like chain tension model global declarations>>=
4206 extern int num_wlc_params;
4207 extern char *wlc_param_descriptions[];
4208 extern char wlc_param_string[];
4210 <<worm-like chain tension model globals>>=
4211 int num_wlc_params = 2;
4212 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
4213 char wlc_param_string[]="0.39e-9,28.4e-9";
4216 <<worm-like chain tension model getopt>>=
4217 {wlc_handler, inverse_wlc_handler, "wlc", "a worm-like chain (F=kbT/p{.25[(1-x/L)^-2 -1] + x/L})", num_wlc_params, wlc_param_descriptions, (char *)wlc_param_string, &string_create_wlc_param_t, &destroy_wlc_param_t}
4219 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
4221 \subsection{Freely-jointed chain}
4224 <<freely-jointed chain functions>>=
4225 <<freely-jointed chain function>>
4226 <<freely-jointed chain handler>>
4227 <<freely-jointed chain create/destroy functions>>
4230 <<freely-jointed chain tension model declarations>>=
4231 <<freely-jointed chain handler declaration>>
4232 <<freely-jointed chain create/destroy function declarations>>
4233 <<freely-jointed chain tension model global declarations>>
4237 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
4238 The tension of a chain of $N$ such links is given by
4240 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
4242 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}.
4243 <<freely-jointed chain function>>=
4244 double langevin(double x, void *params)
4247 return 1.0/tanh(x) - 1/x;
4249 <<one dimensional bracket declaration>>
4250 <<one dimensional solve declaration>>
4251 double inv_langevin(double x)
4253 double lb=0.0, ub=1.0, ret;
4254 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
4255 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
4259 double fjc(double x, double T, double l, int N)
4261 double L = l*(double)N;
4262 if (x == 0) return 0;
4263 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
4264 return k_B*T/l * inv_langevin(x/L);
4268 <<freely-jointed chain handler declaration>>=
4269 extern double fjc_handler(double x, void *pdata);
4271 <<freely-jointed chain handler>>=
4272 double fjc_handler(double x, void *pdata)
4274 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4275 list_t *list = data->group_tension_params;
4280 assert(list != NULL); /* empty active group?! */
4281 l = ((fjc_param_t *)list->d)->l;
4282 while (list != NULL) {
4283 /* enforce consistent l */
4284 assert( ((fjc_param_t *)list->d)->l == l);
4285 N += ((fjc_param_t *)list->d)->N;
4287 fprintf(stderr, "%s: FJC section %d links long in series\n",
4288 __FUNCTION__, ((fjc_param_t *)list->d)->N);
4293 fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
4294 __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
4296 return fjc(x, data->env->T, l, N);
4300 <<freely-jointed chain structure>>=
4301 typedef struct fjc_param_struct {
4302 double l; /* FJC link length (m) */
4303 int N; /* FJC number of links */
4307 <<freely-jointed chain create/destroy function declarations>>=
4308 extern void *string_create_fjc_param_t(char **param_strings);
4309 extern void destroy_fjc_param_t(void *p);
4312 <<freely-jointed chain create/destroy functions>>=
4313 fjc_param_t *create_fjc_param_t(double l, double N)
4315 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
4316 assert(ret != NULL);
4322 void *string_create_fjc_param_t(char **param_strings)
4324 return create_fjc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4325 safe_strtod(param_strings[1], __FUNCTION__));
4328 void destroy_fjc_param_t(void *p)
4335 <<freely-jointed chain tension model global declarations>>=
4336 extern int num_fjc_params;
4337 extern char *fjc_param_descriptions[];
4338 extern char fjc_param_string[];
4341 <<freely-jointed chain tension model globals>>=
4342 int num_fjc_params = 2;
4343 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
4344 char fjc_param_string[]="0.5e-9,200";
4347 <<freely-jointed chain tension model getopt>>=
4348 {fjc_handler, NULL, "fjc", "a freely-jointed chain [F=kbT/l L^-1(x/L)]", num_fjc_params, fjc_param_descriptions, (char *)fjc_param_string, &string_create_fjc_param_t, &destroy_fjc_param_t}
4351 \subsection{Tension model implementation}
4353 <<tension-model.c>>=
4356 <<tension model includes>>
4357 #include "tension_model.h"
4358 <<tension model internal definitions>>
4359 <<tension model globals>>
4360 <<tension model functions>>
4363 <<tension model includes>>=
4364 #include <assert.h> /* assert() */
4365 #include <stdlib.h> /* NULL, strto*() */
4366 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4367 #include <stdio.h> /* fprintf(), stdout */
4368 #include <errno.h> /* errno, ERANGE */
4371 #include "tension_balance.h" /* oneD_*() */
4374 <<tension model internal definitions>>=
4375 <<constant tension structure>>
4377 <<worm-like chain structure>>
4378 <<freely-jointed chain structure>>
4381 <<tension model globals>>=
4382 <<tension handler array global>>
4383 <<constant tension model globals>>
4384 <<hooke tension model globals>>
4385 <<worm-like chain tension model globals>>
4386 <<freely-jointed chain tension model globals>>
4389 <<tension model functions>>=
4391 <<constant tension functions>>
4393 <<worm-like chain functions>>
4394 <<freely-jointed chain functions>>
4398 \subsection{Utilities}
4400 The tension models can get complicated, and users may want to reassure
4401 themselves that this portion of the input physics are functioning properly. The
4402 stand-alone program [[tension_model_utils]] demonstrates the tension model
4403 interface without the overhead of having to understand a full simulation.
4404 [[tension_model_utils]] takes command line model arguments like the full
4405 simulation, and prints $F(x)$ data to the screen for the selected model over a
4408 <<tension-model-utils.c>>=
4410 <<tension model utility includes>>
4411 <<tension model utility definitions>>
4412 <<tension model utility getopt functions>>
4414 <<tension model utility main>>
4417 <<tension model utility main>>=
4418 int main(int argc, char **argv)
4420 <<tension model getopt array definition>>
4421 tension_model_getopt_t *model = NULL;
4425 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
4426 tension_handler_data_t tdata;
4427 int num_param_args; /* for INIT_MODEL() */
4428 char **param_args; /* for INIT_MODEL() */
4430 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model,
4431 &Fmax, &Xmax, &flags);
4433 if (flags & VFLAG) {
4434 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4436 INIT_MODEL("utility", model, model->params, params);
4437 tdata.group_tension_params = NULL;
4438 push(&tdata.group_tension_params, params, 1);
4440 tdata.persist = NULL;
4441 if (model->handler == NULL) {
4442 printf("No tension function for model %s\n", model->name);
4448 printf("#Distance (m)\tForce (N)\n");
4449 for (i=0; i<=N; i++) {
4450 x = Xmax*i/(double)N;
4451 F = (*model->handler)(x, &tdata);
4452 if (F < 0 || F > Fmax) break;
4453 printf("%g\t%g\n", x, F);
4455 if (flags & VFLAG && i <= N) {
4456 /* explain exit condition */
4458 printf("Impossible force %g\n", F);
4460 printf("Reached large force limit %g > %g\n", F, Fmax);
4463 params = pop(&tdata.group_tension_params);
4464 if (model->destructor)
4465 (*model->destructor)(params);
4470 <<tension model utility includes>>=
4473 #include <string.h> /* strlen() */
4474 #include <assert.h> /* assert() */
4475 #include <errno.h> /* errno, ERANGE */
4479 #include "tension_model.h"
4482 <<tension model utility definitions>>=
4483 <<version definition>>
4484 #define VFLAG 1 /* verbose */
4485 <<string comparison definition and globals>>
4486 <<tension model getopt definitions>>
4487 <<initialize model definition>>
4491 <<tension model utility getopt functions>>=
4493 <<index tension model>>
4494 <<tension model utility help>>
4495 <<tension model utility get options>>
4498 <<tension model utility help>>=
4499 void help(char *prog_name,
4501 int n_tension_models, tension_model_getopt_t *tension_models,
4502 int tension_model, double Fmax, double Xmax)
4505 printf("usage: %s [options]\n", prog_name);
4506 printf("Version %s\n\n", VERSION);
4507 printf("A utility for understanding the available tension models\n\n");
4508 printf("Environmental options:\n");
4509 printf("-T T\tTemperature T (currently %g K)\n", env->T);
4510 printf("-C T\tYou can also set the temperature T in Celsius\n");
4511 printf("Model options:\n");
4512 printf("-m model\tFolded tension model (currently %s)\n",
4513 tension_models[tension_model].name);
4514 printf("-a args\tFolded tension model argument string (currently %s)\n",
4515 tension_models[tension_model].params);
4516 printf("F(x) is calculated for a range of x and printed\n");
4517 printf("For example:\n");
4518 printf(" #Distance (m)\tForce (N)\n");
4519 printf(" 123.456\t7.89\n");
4521 printf("-F\tSet the maximum F value for the standard mode F(x) output (currently %g)\n", Fmax);
4522 printf("-X\tSet the maximum x value for the standard mode F(x) output (currently %g)\n", Xmax);
4523 printf("-V\tChange output to verbose mode\n");
4524 printf("-h\tPrint this help and exit\n");
4526 printf("Tension models:\n");
4527 for (i=0; i<n_tension_models; i++) {
4528 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
4529 for (j=0; j < tension_models[i].num_params ; j++)
4530 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
4531 printf("\t\tdefaults: %s\n", tension_models[i].params);
4536 <<tension model utility get options>>=
4537 void get_options(int argc, char **argv, environment_t *env,
4538 int n_tension_models, tension_model_getopt_t *tension_models,
4539 tension_model_getopt_t **model, double *Fmax, double *Xmax,
4540 unsigned int *flags)
4542 char *prog_name = NULL;
4543 char c, options[] = "T:C:m:a:F:X:Vh";
4544 int tension_model=0, num_strings;
4545 extern char *optarg;
4546 extern int optind, optopt, opterr;
4550 /* setup defaults */
4552 prog_name = argv[0];
4553 env->T = 300.0; /* K */
4554 *Fmax = 1e5; /* N */
4555 *Xmax = 1e-6; /* m */
4557 *model = tension_models;
4559 while ((c=getopt(argc, argv, options)) != -1) {
4561 case 'T': env->T = safe_strtod(optarg, "-T"); break;
4562 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
4564 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
4565 *model = tension_models+tension_model;
4568 tension_models[tension_model].params = optarg;
4570 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
4571 case 'X': *Xmax = safe_strtod(optarg, "-X"); break;
4572 case 'V': *flags |= VFLAG; break;
4574 fprintf(stderr, "unrecognized option '%c'\n", optopt);
4575 /* fall through to default case */
4577 help(prog_name, env, n_tension_models, tension_models, tension_model, *Fmax, *Xmax);
4586 \section{Unfolding rate models}
4587 \label{sec.k_models}
4589 We have two main choices for determining $k$: Bell's model, which assumes the
4590 folded domains exist in a single `folded' state and make instantaneous,
4591 stochastic transitions over a free energy barrier; and Kramers' model, which
4592 treats the unfolding as a thermalized diffusion process.
4593 We deal with these options like we dealt with the different tension models: by bundling all model-specific
4594 parameters into structures, with model specific functions for determining $k$.
4596 <<k func definition>>=
4597 typedef double k_func_t(double F, environment_t *env, void *params);
4600 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
4601 The interface is defined in [[k_model.h]], the implementation in [[k_model.c]], the unit testing in [[check_k_model.c]], and some other tidbits in [[k_model_utils.c]].
4603 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
4604 because \LaTeX\ doesn't like underscores outside math mode.
4605 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
4606 You could use spaces instead (`\verb+ +'),
4607 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
4608 which I don't like as much.
4614 <<k func definition>>
4615 <<null k declarations>>
4616 <<const k declarations>>
4617 <<bell k declarations>>
4618 <<kbell k declarations>>
4619 <<kramers k declarations>>
4620 <<kramers integ k declarations>>
4621 #endif /* K_MODEL_H */
4624 <<k model module makefile lines>>=
4625 NW_SAWSIM_MODS += k_model
4626 CHECK_BINS += check_k_model
4627 check_k_model_MODS = parse string_eval
4629 [[check_k_model]] also depends on the parse module.
4631 <<k model utils makefile lines>>=
4632 K_MODEL_MODS = k_model parse string_eval
4633 K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
4634 $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h)
4635 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
4636 notangle -Rk-model-utils.c $< > $@
4637 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
4638 gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
4639 clean_k_model_utils :
4640 rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
4646 For domains that are never folded.
4648 <<null k declarations>>=
4650 <<null k functions>>=
4655 <<null k model getopt>>=
4656 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
4659 \subsection{Constant rate model}
4662 This is a pretty boring model, but it is usefull for testing the engine.
4666 where $k_0$ is the constant unfolding rate.
4668 <<const k functions>>=
4669 <<const k function>>
4670 <<const k structure create/destroy functions>>
4673 <<const k declarations>>=
4674 <<const k function declaration>>
4675 <<const k structure create/destroy function declarations>>
4676 <<const k global declarations>>
4680 <<const k structure>>=
4681 typedef struct const_k_param_struct {
4682 double knot; /* transition rate k_0 (frac population per s) */
4686 <<const k function declaration>>=
4687 double const_k(double F, environment_t *env, void *const_k_params);
4689 <<const k function>>=
4690 double const_k(double F, environment_t *env, void *const_k_params)
4691 { /* Returns the rate constant k in frac pop per s. */
4692 const_k_param_t *p = (const_k_param_t *) const_k_params;
4694 assert(p->knot > 0);
4699 <<const k structure create/destroy function declarations>>=
4700 void *string_create_const_k_param_t(char **param_strings);
4701 void destroy_const_k_param_t(void *p);
4704 <<const k structure create/destroy functions>>=
4705 const_k_param_t *create_const_k_param_t(double knot)
4707 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
4708 assert(ret != NULL);
4713 void *string_create_const_k_param_t(char **param_strings)
4715 return create_const_k_param_t(safe_strtod(param_strings[0], __FUNCTION__));
4718 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
4725 <<const k global declarations>>=
4726 extern int num_const_k_params;
4727 extern char *const_k_param_descriptions[];
4728 extern char const_k_param_string[];
4731 <<const k globals>>=
4732 int num_const_k_params = 1;
4733 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
4734 char const_k_param_string[]="1";
4737 <<const k model getopt>>=
4738 {"const", "constant rate transitions", &const_k, num_const_k_params, const_k_param_descriptions, (char *)const_k_param_string, &string_create_const_k_param_t, &destroy_const_k_param_t}
4741 \subsection{Bell's model}
4744 Quantitatively, Bell's model gives $k$ as
4746 k = k_0 \cdot e^{F dx / k_B T} \;,
4748 where $k_0$ is the unforced unfolding rate,
4749 $F$ is the applied unfolding force,
4750 $dx$ is the distance to the transition state, and
4751 $k_B T$ is the thermal energy\citep{schlierf06}.
4753 <<bell k functions>>=
4755 <<bell k structure create/destroy functions>>
4758 <<bell k declarations>>=
4759 <<bell k function declaration>>
4760 <<bell k structure create/destroy function declarations>>
4761 <<bell k global declarations>>
4765 <<bell k structure>>=
4766 typedef struct bell_param_struct {
4767 double knot; /* transition rate k_0 (frac population per s) */
4768 double dx; /* `distance to transition state' (nm) */
4772 <<bell k function declaration>>=
4773 double bell_k(double F, environment_t *env, void *bell_params);
4775 <<bell k function>>=
4776 double bell_k(double F, environment_t *env, void *bell_params)
4777 { /* Returns the rate constant k in frac pop per s.
4778 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4779 * uses global k_B in J/K */
4780 bell_param_t *p = (bell_param_t *) bell_params;
4781 assert(F >= 0); assert(env->T > 0);
4783 assert(p->knot > 0); //assert(p->dx >= 0); /* p < 0 for refolding reactions, etc. */
4785 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
4786 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
4787 p->knot * exp(F*p->dx / (k_B*env->T) ));
4788 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
4790 return p->knot * exp(F*p->dx / (k_B*env->T) );
4794 <<bell k structure create/destroy function declarations>>=
4795 void *string_create_bell_param_t(char **param_strings);
4796 void destroy_bell_param_t(void *p);
4799 <<bell k structure create/destroy functions>>=
4800 bell_param_t *create_bell_param_t(double knot, double dx)
4802 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
4803 assert(ret != NULL);
4809 void *string_create_bell_param_t(char **param_strings)
4811 return create_bell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4812 safe_strtod(param_strings[1], __FUNCTION__));
4815 void destroy_bell_param_t(void *p /* bell_param_t * */)
4822 <<bell k global declarations>>=
4823 extern int num_bell_params;
4824 extern char *bell_param_descriptions[];
4825 extern char bell_param_string[];
4829 int num_bell_params = 2;
4830 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4831 char bell_param_string[]="3.3e-4,0.25e-9";
4834 <<bell k model getopt>>=
4835 {"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}
4837 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4838 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4841 \subsection{Linker stiffness corrected Bell model}
4844 Walton et. al showed that the Bell model constant force approximation
4845 doesn't hold for biotin-streptavadin binding, and I extended their
4846 results to I27 for the 2009 Biophysical Society Annual
4847 Meeting\cite{walton08,king09}. More details to follow when I get done
4848 with the conference\ldots
4850 We adjust Bell's model to
4852 k = k_0 \cdot e^{\frac{F dx - \frac{1}{2} \kappa dx^2}{k_B T}} \;,
4854 where $k_0$ is the unforced unfolding rate,
4855 $F$ is the applied unfolding force,
4856 $\kappa$ is the effective spring constant,
4857 $dx$ is the distance to the transition state, and
4858 $k_B T$ is the thermal energy\citep{schlierf06}.
4860 <<kbell k functions>>=
4861 <<kbell k function>>
4862 <<kbell k structure create/destroy functions>>
4865 <<kbell k declarations>>=
4866 <<kbell k function declaration>>
4867 <<kbell k structure create/destroy function declarations>>
4868 <<kbell k global declarations>>
4872 <<kbell k structure>>=
4873 typedef struct kbell_param_struct {
4874 double knot; /* transition rate k_0 (frac population per s) */
4875 double dx; /* `distance to transition state' (nm) */
4879 <<kbell k function declaration>>=
4880 double kbell_k(double F, environment_t *env, void *kbell_params);
4882 <<kbell k function>>=
4883 double kbell_k(double F, environment_t *env, void *kbell_params)
4884 { /* Returns the rate constant k in frac pop per s.
4885 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4886 * uses global k_B in J/K */
4887 kbell_param_t *p = (kbell_param_t *) kbell_params;
4888 assert(F >= 0); assert(env->T > 0); assert(env->stiffness >= 0);
4890 assert(p->knot > 0); assert(p->dx >= 0);
4891 return p->knot * exp((F*p->dx - 0.5*env->stiffness*(p->dx)*(p->dx))/(k_B*env->T) );
4895 <<kbell k structure create/destroy function declarations>>=
4896 void *string_create_kbell_param_t(char **param_strings);
4897 void destroy_kbell_param_t(void *p);
4900 <<kbell k structure create/destroy functions>>=
4901 kbell_param_t *create_kbell_param_t(double knot, double dx)
4903 kbell_param_t *ret = (kbell_param_t *) malloc(sizeof(kbell_param_t));
4904 assert(ret != NULL);
4910 void *string_create_kbell_param_t(char **param_strings)
4912 return create_kbell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4913 safe_strtod(param_strings[1], __FUNCTION__));
4916 void destroy_kbell_param_t(void *p /* kbell_param_t * */)
4923 <<kbell k global declarations>>=
4924 extern int num_kbell_params;
4925 extern char *kbell_param_descriptions[];
4926 extern char kbell_param_string[];
4929 <<kbell k globals>>=
4930 int num_kbell_params = 2;
4931 char *kbell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4932 char kbell_param_string[]="3.3e-4,0.25e-9";
4935 <<kbell k model getopt>>=
4936 {"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}
4938 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4939 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4942 \subsection{Kramer's model}
4945 Kramer's model gives $k$ as
4947 % \frac{1}{k} = \frac{1}{D}
4948 % \int_{x_\text{min}}^{x_\text{max}}
4949 % e^{\frac{-U_F(x)}{k_B T}}
4950 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4953 %where $D$ is the diffusion constant,
4954 %$U_F(x)$ is the free energy along the unfolding coordinate, and
4955 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4956 % before the minimum and after the maximum, respectively \citep{schlierf06}.
4958 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
4960 where $D$ is the diffusion constant,
4962 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
4964 are length scales where
4965 $x_c(F)$ is the location of the bound state and
4966 $x_{ts}(F)$ is the location of the transition state,
4967 $E(F,x)$ is the free energy, and
4968 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
4969 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
4970 \citet{evans97} Eqn.~3).
4972 <<kramers k functions>>=
4973 <<kramers k function>>
4974 <<kramers k structure create/destroy functions>>
4977 <<kramers k declarations>>=
4978 <<kramers k function declaration>>
4979 <<kramers k structure create/destroy function declarations>>
4980 <<kramers k global declarations>>
4984 <<kramers k structure>>=
4985 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
4986 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
4988 typedef struct kramers_param_struct {
4989 double D; /* diffusion rate (um/s) */
4996 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
4997 //kramers_x_func_t *xts; /* function returning transition x (nm) */
4998 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
4999 //kramers_E_func_t *E; /* function returning E (J) */
5000 //double *E_params; /* parametrize the energy landscape */
5001 //int n_E_params; /* length of E_params array */
5005 <<kramers k function declaration>>=
5006 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
5007 extern double kramers_k(double F, environment_t *env, void *kramers_params);
5009 <<kramers k function>>=
5010 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
5012 static char input[160]={0};
5013 static char output[80]={0};
5014 /* setup the environment */
5015 fprintf(in, "F = %g; T = %g\n", F, T);
5016 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
5017 string_eval(in, out, input, 80, output);
5018 return safe_strtod(output, __FUNCTION__);
5021 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
5023 static char input[160]={0};
5024 static char output[80]={0};
5025 /* setup the environment */
5026 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
5027 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
5028 string_eval(in, out, input, 80, output);
5029 return safe_strtod(output, __FUNCTION__);
5032 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
5034 kramers_param_t *p = (kramers_param_t *) kramers_params;
5035 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
5038 double kramers_k(double F, environment_t *env, void *kramers_params)
5039 { /* Returns the rate constant k in frac pop per s.
5040 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5041 * uses global k_B in J/K */
5042 kramers_param_t *p = (kramers_param_t *) kramers_params;
5043 double kbT, xc, xts, lc, lts, Eb;
5044 assert(F >= 0); assert(env->T > 0);
5047 assert(p->xc != NULL); assert(p->xts != NULL);
5048 assert(p->ddEdxx != NULL); assert(p->E != NULL);
5049 assert(p->in != NULL); assert(p->out != NULL);
5051 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
5052 if (xc == -1.0) return -1.0; /* error (Too much force) */
5053 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
5054 if (xc == -1.0) return -1.0; /* error (Too much force) */
5055 lc = sqrt(2.0*M_PI*kbT /
5056 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
5057 lts = sqrt(-2.0*M_PI*kbT/
5058 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
5059 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
5060 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
5061 //fprintf(stderr,"D = %g, lc = %g, lts = %g, Eb = %g, kbT = %g, k = %g\n", p->D, lc, lts, Eb, kbT, p->D/(lc*lts) * exp(-Eb/kbT));
5062 return p->D/(lc*lts) * exp(-Eb/kbT);
5066 <<kramers k structure create/destroy function declarations>>=
5067 void *string_create_kramers_param_t(char **param_strings);
5068 void destroy_kramers_param_t(void *p);
5071 <<kramers k structure create/destroy functions>>=
5072 kramers_param_t *create_kramers_param_t(double D,
5073 char *xc_fn, char *xts_fn,
5074 char *E_fn, char *ddEdxx_fn,
5075 char *global_define_string)
5076 // kramers_x_func_t *xc, kramers_x_func_t *xts,
5077 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
5078 // double *E_params, int n_E_params)
5080 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
5081 assert(ret != NULL);
5082 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5083 assert(ret->xc != NULL);
5084 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5085 assert(ret->xts != NULL);
5086 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
5087 assert(ret->E != NULL);
5088 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
5089 assert(ret->ddEdxx != NULL);
5091 strcpy(ret->xc, xc_fn);
5092 strcpy(ret->xts, xts_fn);
5093 strcpy(ret->E, E_fn);
5094 strcpy(ret->ddEdxx, ddEdxx_fn);
5095 //ret->E_params = E_params;
5096 //ret->n_E_params = n_E_params;
5097 string_eval_setup(&ret->in, &ret->out);
5098 fprintf(ret->in, "kB = %g\n", k_B);
5099 fprintf(ret->in, "%s\n", global_define_string);
5103 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
5104 void *string_create_kramers_param_t(char **param_strings)
5106 return create_kramers_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5114 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
5116 kramers_param_t *pk = (kramers_param_t *)p;
5118 string_eval_teardown(&pk->in, &pk->out);
5120 // free(pk->E_params);
5126 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
5127 Schlierf and Rief used a cubic-spline interpolation routine and the double integral form of Kramers' theory, so we get to pick an actual function to approximate the energy landscape.
5128 We follow \citet{shillcock98} and use
5130 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
5131 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
5134 where TODO, variable meanings.%\citep{shillcock98}.
5135 The first and second derivatives of $E(F,E_b,x,x_s)$ are
5137 E'(F,E_b,x,x_s) &=\frac{27 E_b}{4 x_s}\frac{x}{x_s}\p({4-9\frac{x}{x_s}})-F\\
5138 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
5140 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
5141 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
5142 The roots are therefor at
5144 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
5145 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
5146 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
5149 <<kramers k global declarations>>=
5150 extern int num_kramers_params;
5151 extern char *kramers_param_descriptions[];
5152 extern char kramers_param_string[];
5155 <<kramers k globals>>=
5156 int num_kramers_params = 6;
5157 char *kramers_param_descriptions[] = {"Diffusion constant D, m^2/s", "constant parameter declarations", "bound position xc, m", "transition state xts, m","energy E, J", "d^2E/dx^2, J/m^2"};
5158 char kramers_param_string[]="0.2e-12,{Eb = 20*300*kB;xs = 1.5e-9},{(F*xs > 3*Eb) and (-1) or (2*xs/9.0*(1 - (1 - F*xs/(3*Eb))**0.5))},{2*xs/9.0*(1 + (1 - F*xs/(3*Eb))**0.5)},{27.0/4 * Eb * (x/xs)**2 * (2-3*x/xs) - F*x},{27.0*Eb/(2*xs**2) * (2-9*x/xs)}";
5160 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
5161 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
5162 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
5163 Returning an $x_c$ position of $-1\U{m}$ lets the simulation know that the rate is poorly defined for that force, and the timestepping algorithm in Section \ref{sec.adaptive_dt} reduces it's timestep appropriately.
5164 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
5165 It works so long as [[val_1]] is not [[false]].
5167 <<kramers k model getopt>>=
5168 {"kramers", "thermalized, diffusion-limited transitions (function parameters are parsed by Python. For distances, the environment variables force 'F' and temperature 'T' are defined, as well as the Boltzmann constant 'kB'. For energies, the position 'x' is also defined. Functions must evaluate to a float representing their output and produce no output on stdout.", &kramers_k, num_kramers_params, kramers_param_descriptions, (char *)kramers_param_string, &string_create_kramers_param_t, &destroy_kramers_param_t}
5171 \subsection{Kramer's model (natural cubic splines)}
5172 \label{sec.kramers_integ}
5174 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
5175 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
5176 \citet{schlierf06} Eqn.~2)
5178 \frac{1}{k} = \frac{1}{D}
5179 \int_{x_\text{min}}^{x_\text{max}}
5180 e^{\frac{U_F(x)}{k_B T}}
5181 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5184 where $D$ is the diffusion constant,
5185 $U_F(x)$ is the free energy along the unfolding coordinate, and
5186 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
5187 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
5189 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
5190 interpolating between them with cubic splines.
5191 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
5193 <<kramers integ k functions>>=
5194 <<spline functions>>
5195 <<kramers integ A k functions>>
5196 <<kramers integ B k functions>>
5197 <<kramers integ k function>>
5198 <<kramers integ k structure create/destroy functions>>
5201 <<kramers integ k declarations>>=
5202 <<kramers integ k function declaration>>
5203 <<kramers integ k structure create/destroy function declarations>>
5204 <<kramers integ k global declarations>>
5208 <<kramers integ k structure>>=
5209 typedef double func_t(double x, void *params);
5210 typedef struct kramers_integ_param_struct {
5211 double D; /* diffusion rate (um/s) */
5212 double x_min; /* integration bounds */
5214 func_t *E; /* function returning E (J) */
5215 void *E_params; /* parametrize the energy landscape */
5216 destroy_data_func_t *destroy_E_params;
5218 double F; /* for passing into GSL evaluations */
5220 } kramers_integ_param_t;
5223 <<spline param structure>>=
5224 typedef struct spline_param_struct {
5225 int n_knots; /* natural cubic spline knots for U(x) */
5228 gsl_spline *spline; /* GSL spline parameters */
5229 gsl_interp_accel *acc; /* GSL spline acceleration data */
5233 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
5235 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5237 <<kramers integ A k functions>>=
5238 double kramers_integ_k_integrandA(double x, void *params)
5240 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5241 double U = (*p->E)(x, p->E_params);
5242 double Fx = p->F * x;
5243 double kbT = k_B * p->env->T;
5244 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
5245 return exp(-(U-Fx)/kbT);
5247 double kramers_integ_k_integralA(double x, void *params)
5249 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5251 double result, abserr;
5252 size_t neval; /* for qng */
5253 static gsl_integration_workspace *w = NULL;
5254 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5255 f.function = &kramers_integ_k_integrandA;
5257 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5258 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5259 //fprintf(stderr, "integralA = %g\n", result);
5265 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
5266 \int_{x_\text{min}}^{x_\text{max}}
5267 e^{\frac{U_F(x)}{k_B T}}
5268 \text{Integral}_A(x)
5271 <<kramers integ B k functions>>=
5272 double kramers_integ_k_integrandB(double x, void *params)
5274 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5275 double U = (*p->E)(x, p->E_params);
5276 double Fx = p->F * x;
5277 double kbT = k_B * p->env->T;
5278 double intA = kramers_integ_k_integralA(x, params);
5279 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
5280 return exp((U-Fx)/kbT)*intA;
5282 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
5284 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5286 double result, abserr;
5287 size_t neval; /* for qng */
5288 static gsl_integration_workspace *w = NULL;
5289 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5290 f.function = &kramers_integ_k_integrandB;
5294 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5295 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5296 //fprintf(stderr, "integralB = %g\n", result);
5301 With the integrals taken care of, evaluating $k$ is simply
5303 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
5305 <<kramers integ k function declaration>>=
5306 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
5308 <<kramers integ k function>>=
5309 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
5310 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
5311 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
5312 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
5316 <<kramers integ k structure create/destroy function declarations>>=
5317 void *string_create_kramers_integ_param_t(char **param_strings);
5318 void destroy_kramers_integ_param_t(void *p);
5321 <<kramers integ k structure create/destroy functions>>=
5322 kramers_integ_param_t *create_kramers_integ_param_t(double D,
5323 double xmin, double xmax, func_t *E,
5325 destroy_data_func_t *destroy_E_params)
5327 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
5328 assert(ret != NULL);
5333 ret->E_params = E_params;
5334 ret->destroy_E_params = destroy_E_params;
5338 void *string_create_kramers_integ_param_t(char **param_strings)
5340 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
5341 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
5342 void *E_params = string_create_spline_param_t(param_strings+1);
5343 return create_kramers_integ_param_t(
5344 safe_strtod(param_strings[0], __FUNCTION__),
5345 safe_strtod(param_strings[2], __FUNCTION__),
5346 safe_strtod(param_strings[3], __FUNCTION__),
5347 &spline_eval, E_params, destroy_spline_param_t);
5350 void destroy_kramers_integ_param_t(void *params)
5352 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5355 (*p->destroy_E_params)(p->E_params);
5361 Finally we have the GSL spline wrappers:
5362 <<spline functions>>=
5364 <<create/destroy spline>>
5367 <<spline function>>=
5368 double spline_eval(double x, void *spline_params)
5370 spline_param_t *p = (spline_param_t *)(spline_params);
5371 return gsl_spline_eval(p->spline, x, p->acc);
5375 <<create/destroy spline>>=
5376 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
5378 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
5379 assert(ret != NULL);
5380 ret->n_knots = n_knots;
5383 ret->acc = gsl_interp_accel_alloc();
5384 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
5385 gsl_spline_init(ret->spline, x, y, n_knots);
5389 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
5390 void *string_create_spline_param_t(char **param_strings)
5394 double *x=NULL, *y=NULL;
5395 /* break into ordered pairs */
5396 parse_list_string(param_strings[0], SEP, '(', ')',
5397 &num_knots, &knot_coords);
5398 x = (double *)malloc(sizeof(double)*num_knots);
5400 y = (double *)malloc(sizeof(double)*num_knots);
5402 for (i=0; i<num_knots; i++) {
5403 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
5404 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
5407 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
5409 free_parsed_list(num_knots, knot_coords);
5410 return create_spline_param_t(num_knots, x, y);
5413 void destroy_spline_param_t(void *params)
5415 spline_param_t *p = (spline_param_t *)params;
5418 gsl_spline_free(p->spline);
5420 gsl_interp_accel_free(p->acc);
5430 <<kramers integ k global declarations>>=
5431 extern int num_kramers_integ_params;
5432 extern char *kramers_integ_param_descriptions[];
5433 extern char kramers_integ_param_string[];
5436 <<kramers integ k globals>>=
5437 int num_kramers_integ_params = 4;
5438 char *kramers_integ_param_descriptions[] = {
5439 "diffusion D, m^2 / s",
5440 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
5441 "starting integration bound x_min, m",
5442 "ending integration bound x_max, m"};
5443 //char kramers_integ_param_string[]="0.2e-12,{(3.5e-9,-5*k_B*300),(4e-9,-20*k_B*300),(4.25e-9,-8*k_B*300),(5e-9,0),(5.5e-9,-10*k_B*300)},3e-9,6e-9";
5444 //char kramers_integ_param_string[]="0.2e-12,{(3.5e-9,-2.1e-20),(4e-9,-8.3e-20),(4.25e-9,-3.3e-20),(5e-9,0),(5.5e-9,-4.1e-20)},3e-9,6e-9";
5445 char kramers_integ_param_string[]="0.2e-12,{(2e-9,0),(2.1e-9,-3.7e-20),(2.2e-9,-4.4e-20),(2.35e-9,-2.5e-20),(2.41e-9,-7e-21),(2.45e-9,8e-22),(2.51e-9,6.9e-21),(2.64e-9,1.39e-20),(2.9e-9,2.55e-20),(3.52e-9,2.9e-20),(3.7e-9,1.45e-20),(4e-9,-1.7e-20),(6e-9,-2e-20),(9e-9,-3e-20),(1e-8,-3e-20)},2e-9,6e-9";
5446 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
5447 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
5448 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
5451 <<kramers integ k model getopt>>=
5452 {"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}
5454 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
5455 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5457 \subsection{Unfolding model implementation}
5461 <<k model includes>>
5462 #include "k_model.h"
5463 <<k model internal definitions>>
5465 <<k model functions>>
5468 <<k model includes>>=
5469 #include <assert.h> /* assert() */
5470 #include <stdlib.h> /* NULL, malloc(), strto*() */
5471 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
5472 #include <stdio.h> /* fprintf(), stdout */
5473 #include <string.h> /* strlen(), strcpy() */
5474 #include <errno.h> /* errno, ERANGE */
5475 #include <gsl/gsl_integration.h>
5476 #include <gsl/gsl_spline.h>
5481 <<k model internal definitions>>=
5482 <<const k structure>>
5483 <<bell k structure>>
5484 <<kbell k structure>>
5485 <<kramers k structure>>
5486 <<kramers integ k structure>>
5487 <<spline param structure>>
5490 <<k model globals>>=
5495 <<kramers k globals>>
5496 <<kramers integ k globals>>
5499 <<k model functions>>=
5501 <<null k functions>>
5502 <<const k functions>>
5503 <<bell k functions>>
5504 <<kbell k functions>>
5505 <<kramers k functions>>
5506 <<kramers integ k functions>>
5509 \subsection{Unfolding model unit tests}
5511 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
5512 <<check-k-model.c>>=
5513 <<k model check includes>>
5514 <<check relative error>>
5516 <<k model test suite>>
5517 <<main check program>>
5520 <<k model check includes>>=
5521 #include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
5522 #include <stdio.h> /* sprintf() */
5523 #include <assert.h> /* assert() */
5524 #include <math.h> /* exp() */
5525 #include <errno.h> /* errno, ERANGE */
5528 #include "k_model.h"
5531 <<k model test suite>>=
5535 <<k model suite definition>>
5538 <<k model suite definition>>=
5539 Suite *test_suite (void)
5541 Suite *s = suite_create ("k model");
5542 <<const k test case defs>>
5543 <<bell k test case defs>>
5545 <<const k test case adds>>
5546 <<bell k test case adds>>
5551 \subsubsection{Constant}
5553 <<const k test case defs>>=
5554 TCase *tc_const_k = tcase_create("Constant k");
5557 <<const k test case adds>>=
5558 tcase_add_test(tc_const_k, test_const_k_create_destroy);
5559 tcase_add_test(tc_const_k, test_const_k_over_range);
5560 suite_add_tcase(s, tc_const_k);
5565 START_TEST(test_const_k_create_destroy)
5567 char *knot[] = {"1","2","3","4","5","6"};
5568 char *params[] = {knot[0]};
5571 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5572 params[0] = knot[i];
5573 p = string_create_const_k_param_t(params);
5574 destroy_const_k_param_t(p);
5579 START_TEST(test_const_k_over_range)
5581 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
5582 char *knot[] = {"1","2","3","4","5","6"};
5583 char *params[] = {knot[0]};
5586 char param_string[80];
5589 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5590 params[0] = knot[i];
5591 p = string_create_const_k_param_t(params);
5592 for ( F=Fm; F<FM; F+=dF ) {
5593 fail_unless(const_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
5594 "const_k(%g, %g, &{%s}) = %g != %s",
5595 F, T, knot[i], const_k(F, &env, p), knot[i]);
5597 destroy_const_k_param_t(p);
5603 \subsubsection{Bell}
5605 <<bell k test case defs>>=
5606 TCase *tc_bell = tcase_create("Bell's k");
5609 <<bell k test case adds>>=
5610 tcase_add_test(tc_bell, test_bell_k_create_destroy);
5611 tcase_add_test(tc_bell, test_bell_k_at_zero);
5612 tcase_add_test(tc_bell, test_bell_k_at_one);
5613 suite_add_tcase(s, tc_bell);
5617 START_TEST(test_bell_k_create_destroy)
5619 char *knot[] = {"1","2","3","4","5","6"};
5621 char *params[] = {knot[0], dx};
5624 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5625 params[0] = knot[i];
5626 p = string_create_bell_param_t(params);
5627 destroy_bell_param_t(p);
5632 START_TEST(test_bell_k_at_zero)
5634 double F=0.0, T=300.0;
5635 char *knot[] = {"1","2","3","4","5","6"};
5637 char *params[] = {knot[0], dx};
5640 char param_string[80];
5643 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5644 params[0] = knot[i];
5645 p = string_create_bell_param_t(params);
5646 fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
5647 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
5648 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
5649 destroy_bell_param_t(p);
5654 START_TEST(test_bell_k_at_one)
5657 char *knot[] = {"1","2","3","4","5","6"};
5659 char *params[] = {knot[0], dx};
5660 double F= k_B*T/safe_strtod(dx, __FUNCTION__);
5663 char param_string[80];
5666 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5667 params[0] = knot[i];
5668 p = string_create_bell_param_t(params);
5669 CHECK_ERR(1e-6, safe_strtod(knot[i], __FUNCTION__)*exp(1.0), bell_k(F, &env, p));
5670 //fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__)*exp(1.0),
5671 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
5672 // F, T, knot[i], dx, bell_k(F, &env, p), safe_strtod(knot[i], __FUNCTION__)*exp(1.0));
5673 destroy_bell_param_t(p);
5679 <<kramers k tests>>=
5682 <<kramers k test case def>>=
5685 <<kramers k test case add>>=
5688 <<k model function tests>>=
5689 <<check relative error>>
5691 START_TEST(test_something)
5693 double k=5, x=3, last_x=2.0, F;
5694 one_dim_fn_t *handlers[] = {&hooke};
5695 void *data[] = {&k};
5696 double xi[] = {0.0};
5698 int new_active_groups = 1;
5699 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
5700 fail_unless(F = k*x, NULL);
5705 \subsection{Utilities}
5707 The unfolding models can get complicated, and are sometimes hard to explain to others.
5708 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
5709 the overhead of having to understand a full simulation.
5710 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
5711 or special mode, where the operation depends on the specific model selected.
5713 <<k-model-utils.c>>=
5715 <<k model utility includes>>
5716 <<k model utility definitions>>
5717 <<k model utility getopt functions>>
5718 <<k model utility multi model E>>
5719 <<k model utility main>>
5722 <<k model utility main>>=
5723 int main(int argc, char **argv)
5725 <<k model getopt array definition>>
5726 k_model_getopt_t *model = NULL;
5731 int num_param_args; /* for INIT_MODEL() */
5732 char **param_args; /* for INIT_MODEL() */
5734 double special_xmin;
5735 double special_xmax;
5736 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
5737 &Fmax, &special_xmin, &special_xmax, &flags);
5738 if (flags & VFLAG) {
5739 printf("#initializing model %s with parameters %s\n", model->name, model->params);
5741 INIT_MODEL("utility", model, model->params, params);
5744 if (model->k == NULL) {
5745 printf("No k function for model %s\n", model->name);
5749 printf("Fmax = %g <= 0\n", Fmax);
5755 printf("#F (N)\tk (%% pop. per s)\n");
5756 for (i=0; i<=N; i++) {
5757 F = Fmax*i/(double)N;
5758 k = (*model->k)(F, &env, params);
5760 printf("%g\t%g\n", F, k);
5765 if (model->k == NULL || model->k == &bell_k) {
5766 printf("No special mode for model %s\n", model->name);
5769 if (special_xmax <= special_xmin) {
5770 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
5774 double Xrng=(special_xmax-special_xmin), x, E;
5776 printf("#x (m)\tE (J)\n");
5777 for (i=0; i<=N; i++) {
5778 x = special_xmin + Xrng*i/(double)N;
5779 E = multi_model_E(model->k, params, &env, x);
5780 printf("%g\t%g\n", x, E);
5787 if (model->destructor)
5788 (*model->destructor)(params);
5793 <<k model utility multi model E>>=
5794 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
5796 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
5798 if (k == kramers_integ_k)
5799 E = (*p->E)(x, p->E_params);
5801 E = kramers_E(0, env, model_params, x);
5807 <<k model utility includes>>=
5810 #include <string.h> /* strlen() */
5811 #include <assert.h> /* assert() */
5812 #include <errno.h> /* errno, ERANGE */
5815 #include "string_eval.h"
5816 #include "k_model.h"
5819 <<k model utility definitions>>=
5820 <<version definition>>
5821 #define VFLAG 1 /* verbose */
5822 enum mode_t {M_K_OF_F, M_SPECIAL};
5823 <<string comparison definition and globals>>
5824 <<k model getopt definitions>>
5825 <<initialize model definition>>
5826 <<kramers k structure>>
5827 <<kramers integ k structure>>
5831 <<k model utility getopt functions>>=
5834 <<k model utility help>>
5835 <<k model utility get options>>
5838 <<k model utility help>>=
5839 void help(char *prog_name,
5841 int n_k_models, k_model_getopt_t *k_models,
5842 int k_model, double Fmax, double special_xmin, double special_xmax)
5845 printf("usage: %s [options]\n", prog_name);
5846 printf("Version %s\n\n", VERSION);
5847 printf("A utility for understanding the available unfolding models\n\n");
5848 printf("Environmental options:\n");
5849 printf("-T T\tTemperature T (currently %g K)\n", env->T);
5850 printf("-C T\tYou can also set the temperature T in Celsius\n");
5851 printf("Model options:\n");
5852 printf("-k model\tTransition rate model (currently %s)\n",
5853 k_models[k_model].name);
5854 printf("-K args\tTransition rate model argument string (currently %s)\n",
5855 k_models[k_model].params);
5856 printf("There are two output modes. In standard mode, k(F) is printed\n");
5857 printf("For example:\n");
5858 printf(" #Force (N)\tk (% pop. per s)\n");
5859 printf(" 123.456\t7.89\n");
5861 printf("In special mode, the output depends on the model.\n");
5862 printf("For models defining an energy function E(x), that function is printed\n");
5863 printf(" #Distance (m)\tE (J)\n");
5864 printf(" 0.0012\t0.0034\n");
5866 printf("-m\tChange output to standard mode\n");
5867 printf("-M\tChange output to special mode\n");
5868 printf("-F\tSet the maximum F value for the standard mode k(F) output (currently %g)\n", Fmax);
5869 printf("-x\tSet the minimum x value for the special mode E(x) output (currently %g)\n", special_xmin);
5870 printf("-X\tSet the maximum x value for the special mode E(x) output (currently %g)\n", special_xmax);
5871 printf("-V\tChange output to verbose mode\n");
5872 printf("-h\tPrint this help and exit\n");
5874 printf("Unfolding rate models:\n");
5875 for (i=0; i<n_k_models; i++) {
5876 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
5877 for (j=0; j < k_models[i].num_params ; j++)
5878 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
5879 printf("\t\tdefaults: %s\n", k_models[i].params);
5884 <<k model utility get options>>=
5885 void get_options(int argc, char **argv, environment_t *env,
5886 int n_k_models, k_model_getopt_t *k_models,
5887 enum mode_t *mode, k_model_getopt_t **model,
5888 double *Fmax, double *special_xmin, double *special_xmax,
5889 unsigned int *flags)
5891 char *prog_name = NULL;
5892 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
5894 extern char *optarg;
5895 extern int optind, optopt, opterr;
5899 /* setup defaults */
5901 prog_name = argv[0];
5902 env->T = 300.0; /* K */
5906 *Fmax = 1e-9; /* N */
5907 *special_xmax = 1e-8;
5908 *special_xmin = 0.0;
5910 while ((c=getopt(argc, argv, options)) != -1) {
5912 case 'T': env->T = safe_strtod(optarg, "-T"); break;
5913 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
5915 k_model = index_k_model(n_k_models, k_models, optarg);
5916 *model = k_models+k_model;
5919 k_models[k_model].params = optarg;
5921 case 'm': *mode = M_K_OF_F; break;
5922 case 'M': *mode = M_SPECIAL; break;
5923 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
5924 case 'x': *special_xmin = safe_strtod(optarg, "-x"); break;
5925 case 'X': *special_xmax = safe_strtod(optarg, "-X"); break;
5926 case 'V': *flags |= VFLAG; break;
5928 fprintf(stderr, "unrecognized option '%c'\n", optopt);
5929 /* fall through to default case */
5931 help(prog_name, env, n_k_models, k_models, k_model, *Fmax, *special_xmin, *special_xmax);
5940 \section{\LaTeX\ documentation}
5942 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
5943 The comment blocks are extracted (with nicely formatted code blocks), using
5944 <<latex makefile lines>>=
5945 $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5946 noweave -latex -index -delay $< > $@
5947 $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR)
5951 <<latex makefile lines>>=
5952 $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
5954 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5955 $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
5956 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5957 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5958 mv $(BUILD_DIR)/sawsim.pdf $@
5960 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
5961 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
5962 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
5964 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
5965 <<latex makefile lines>>=
5967 rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
5968 $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
5969 $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.toc \
5970 $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib
5971 clean_latex : semi-clean_latex
5972 rm -f $(DOC_DIR)/sawsim.pdf
5977 [[make]] is a common utility on *nix systems for managing dependencies while building software.
5978 In this case, we have several \emph{targets} that we'd like to build:
5979 the main simulation program \prog;
5980 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
5981 the documentation [[sawsim.pdf]];
5982 and the [[Makefile]] itself.
5983 Besides the generated files, there is the utility target
5984 [[clean]] for removing all generated files except [[Makefile]].
5986 % [[dist]] for generating a distributable tar file.
5988 Extract the makefile with
5989 `[[notangle -Rmakefile src/sawsim.nw | sed 's/ /\t/' > Makefile]]'.
5990 The sed is needed because notangle replaces tabs with eight spaces,
5991 but [[make]] needs tabs.
5993 # Decide what directories to put things in
5998 # Note: Cannot use, for example, `./build', because make eats the `./'
5999 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
6001 # Modules (X.c, X.h) defined in the noweb file
6004 # Modules defined outside the noweb file
6005 FREE_SAWSIM_MODS = interp tavl
6007 <<list module makefile lines>>
6008 <<tension balance module makefile lines>>
6009 <<tension model module makefile lines>>
6010 <<k model module makefile lines>>
6011 <<parse module makefile lines>>
6012 <<string eval module makefile lines>>
6013 <<accel k module makefile lines>>
6015 SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS)
6017 # Everything needed for sawsim under one roof. sawsim.c must come first
6018 SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \
6019 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h)
6020 # Libraries needed to compile sawsim
6021 LIBS = gsl gslcblas m
6022 CHECK_LIBS = gsl gslcblas m check
6023 # The non-check binaries generated
6024 BINS = sawsim tension_model_utils k_model_utils sawsim_profile
6027 # Define the major targets
6028 all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
6030 view : $(DOC_DIR)/sawsim.pdf
6032 profile : $(BIN_DIR)/sawsim_profile | $(BIN_DIR)
6033 $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \
6034 -mnull -Mwlc -A0.39e-9,28e-9 -F8
6035 gprof $(BIN_DIR)/sawsim_profile gmon.out > $@
6036 check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim
6037 $(SHELL) -e -c 'for B in $^; do ./$$B; done'
6038 clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \
6039 clean_tension_model_utils \
6040 clean_k_model_utils clean_latex clean_check_sawsim
6041 rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \
6042 $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \
6043 $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
6044 $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
6045 $(BUILD_DIR)/global.h ./gmon.out
6046 $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
6048 # Various builds of sawsim
6049 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
6050 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6051 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
6052 gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6053 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
6054 gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6056 # Create the directories
6057 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
6060 # Copy over the source external to sawsim
6061 # Note: Cannot use, for example, `./build', because make eats the `./'
6062 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
6064 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6068 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6072 ## Basic source generated with noweb
6073 # The central files sawsim.c and global.h...
6074 $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6076 $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6077 notangle -Rglobal.h $< > $@
6078 # ... and the modules
6079 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6080 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6081 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6082 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6083 $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6084 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6085 # Note: I use `_' as a space in my file names, but noweb doesn't like
6086 # that so we use `-' to name the noweb roots and substitute here.
6088 ## Building the unit-test programs
6090 $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
6091 notangle -Rchecks $< > $@
6092 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
6093 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
6094 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR)
6095 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
6096 clean_check_sawsim :
6097 rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
6098 # ... and the modules
6100 $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
6101 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \
6102 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \
6103 $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6104 $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6105 $$(BUILD_DIR)/global.h | $(BIN_DIR)
6106 gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
6107 $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
6109 # todo: clean up the required modules to
6110 $(CHECK_BINS:%=clean_%) :
6111 rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
6113 # Cleaning up the modules
6115 $(SAWSIM_MODS:%=clean_%) :
6116 rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h)
6118 <<tension model utils makefile lines>>
6119 <<k model utils makefile lines>>
6120 <<latex makefile lines>>
6122 Makefile : $(SRC_DIR)/sawsim.nw
6123 notangle -Rmakefile $< | sed 's/ /\t/' > $@
6125 This is fairly self-explanatory. We compile a dynamically linked
6126 version ([[sawsim]]) and a statically linked version
6127 ([[sawsim_static]]). The static version is about 9 times as big, but
6128 it runs without needing dynamic GSL libraries to link to, so it's a
6129 better format for distributing to a cluster for parallel evaluation.
6133 \subsection{Hookean springs in series}
6134 \label{sec.math_hooke}
6137 The effective spring constant for $n$ springs $k_i$ in series is given by
6139 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6145 F &= k_i x_i &&\forall i \le n \\
6146 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
6147 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
6148 F &= k_1 x_1 = k_\text{eff} x \\
6149 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
6150 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6155 \addcontentsline{toc}{section}{References}
6156 \bibliography{sawsim}