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 definition>>=
162 #define VERSION "0.8"
168 sawsim - program for simulating single molecule mechanical unfolding.
169 Copyright (C) 2008-2009, William Trevor King
171 This program is free software; you can redistribute it and/or
172 modify it under the terms of the GNU General Public License as
173 published by the Free Software Foundation; either version 3 of the
174 License, or (at your option) any later version.
176 This program is distributed in the hope that it will be useful, but
177 WITHOUT ANY WARRANTY; without even the implied warranty of
178 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
179 See the GNU General Public License for more details.
181 You should have received a copy of the GNU General Public License
182 along with this program; if not, write to the Free Software
183 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
186 The author may be contacted at <wking@drexel.edu> on the Internet, or
187 write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
188 Philadelphia PA 19104, USA.
201 The unfolding system is stored as a chain of \emph{domains} (Figure
202 \ref{fig.domain_chain}). Domains can be folded, globular protein
203 domains, unfolded protein linkers, AFM cantilevers, or other
204 stretchable system components. The system is modeled as graph of
205 possible domain states with transitions between them (Figure
206 \ref{fig.domain_states}).
210 \subfloat[][]{\label{fig.domain_chain}
211 \begin{tikzpicture}[->,node distance=2cm,font=\footnotesize]
212 \tikzstyle{every state}=[fill,draw=red!50,very thick,fill=red!20]
213 \node[state] (A) {domain 1};
214 \node[state] (B) [below of=A] {domain 2};
215 \node[state] (C) [below of=B] {.~.~.};
216 \node[state] (D) [below of=C] {domain $N$};
217 \node (S) [below of=D] {Surface};
218 \node (E) [above of=A] {};
220 \path[-] (A) edge (B)
221 (B) edge node [right] {Tension} (C)
224 \path[->,green] (A) edge node [right,black] {Extension} (E);
228 \subfloat[][]{\label{fig.domain_states}
229 \begin{tikzpicture}[->,node distance=2.5cm,shorten <=1pt,shorten >=1pt,font=\footnotesize]
230 \tikzstyle{every state}=[fill,draw=blue!50,very thick,fill=blue!20]
231 \node[state] (A) {cantilever};
232 \node[state] (C) [below of=A] {transition};
233 \node[state] (B) [left of=C] {folded};
234 \node[state] (D) [right of=C] {unfolded};
236 \path (B) edge [bend left] node [above] {$k_1$} (C)
237 (C) edge [bend left] node [below] {$k_1'$} (B)
238 edge [bend left] node [above] {$k_2$} (D)
239 (D) edge [bend left] node [below] {$k_2'$} (C);
242 \caption{\subref{fig.domain_chain} Extending a chain of domains. One end
243 of the chain is fixed, while the other is extended at a constant
244 speed. The domains are coupled with rigid linkers, so the domains
245 themselves must stretch to accomodate the extension.
246 \subref{fig.domain_states} Each domain exists in a discrete state. At
247 each timestep, it may transition into another state following a
248 user-defined state matrix such as this one, showing a metastable
249 transition state and an explicit ``cantilever'' domain.}
253 Each domain contributes to the total tension, which depends on the
254 chain extension. The particular model for the tension as a function
255 of extension is handled generally with function pointers. So far the
256 following models have been implemented:
258 \item Null (Appendix \ref{sec.null_tension}),
259 \item Constant (Appendix \ref{sec.const_tension}),
260 \item Hookean spring (Appendix \ref{sec.hooke}),
261 \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
262 \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
265 The tension model and parameters used for a particular domain depend
266 on the domain's current state. The transition rates between states
267 are also handled generally with function pointers, with
270 \item Null (Appendix \ref{sec.null_k}),
271 \item Bell's model (Appendix \ref{sec.bell}),
272 \item Spring constant corrected Bell model (Appendix \ref{sec.kbell}),
273 \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
274 \item Kramers' saddle point model (Appendix \ref{sec.kramers}).
277 The unfolding of the chain domains is modeled in two stages. First
278 the tension of the chain is calculated. Then the tension (assumed to
279 be constant over the length of the chain, see Section
280 \ref{sec.timescales}), is applied to each domain, and used to
281 calculate the probability of that domain changing state in the next
282 timestep $dt$. Then the time is stepped forward, and the process is
283 repeated. The simulation is complete when a pre-set time limit
284 [[t_max]] is reached or a pre-selected state [[stop_state]] is empty.
286 int main(int argc, char **argv)
288 <<tension model getopt array definition>>
289 <<k model getopt array definition>>
291 /* variables initialized in get_options() */
292 list_t *state_list=NULL, *transition_list=NULL;
293 environment_t env={0};
294 double P, t_max, dt_max, v, x_step;
295 state_t *stop_state=NULL;
297 /* variables used in the simulation loop */
298 double t=0, x=0, dt, F; /* time, distance, timestep, force */
299 int transition=1; /* boolean marking a transition for tension optimization */
301 get_options(argc, argv, &P, &t_max, &dt_max, &v, &x_step, &stop_state, &env,
302 NUM_TENSION_MODELS, tension_models, NUM_K_MODELS, k_models,
303 &state_list, &transition_list);
306 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Distance (m)\tForce (N)\tStiffness (N/m)\n");
307 if (flags & FLAG_OUTPUT_TRANSITION_FORCES) printf("#Force (N)\tInitial state\tFinal state\n");
308 while ((t_max > 0 && t <= t_max) || (stop_state != NULL && stop_state->num_domains > 0)) {
309 F = find_tension(state_list, &env, x, transition);
311 dt = determine_dt(state_list, transition_list, &env, x, P, dt_max, v);
313 dt = determine_dt(state_list, transition_list, &env, x, P, dt_max, 0);
314 transition=random_transitions(transition_list, F, dt, &env);
315 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness);
320 if (dt == dt_max) { /* step completed */
323 } else { /* still working on this step */
328 destroy_state_list(state_list);
329 destroy_transition_list(transition_list);
333 @ The meat of the simulation is bundled into the three functions
334 [[find_tension]], [[determine_dt]], and [[random_transitions]].
335 [[find_tension]] is discussed in Section \ref{sec.find_tension},
336 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
337 [[random_transitions]] in Section \ref{sec.transition_rate}.
339 The stretched distance is extended in one of two modes depending on
340 [[xstep]]. If [[xstep == 0]], then the pulling is continuous, at
341 least within the limits of the inherent discretization of a time
342 stepped approach. At any rate, the timestep is limited by the maximum
343 allowed timestep [[dt_max]] and the maximum allowed unfolding
344 probability in a given step [[P]]. In the second mode [[xstep > 0]],
345 and the end to end distance increases in discrete steps of that
346 length. The time between steps is chosen to maintain an average
347 unfolding velocity [[v]]. This approach is less work to simulate than
348 smooth pulling and also closer to the reality of many experiments, but
349 it is practically impossible to treat analytically. With both pulling
350 styles implemented, the effects of the discretization can be easily
351 calculated, bridging the gap between theory and experiment.
353 Environmental parameters important in determining reaction rates and
354 tensions (e.g.\ temperature, pH) are stored in a single structure to
355 facilitate extension to more complicated models in the future.
356 <<environment definition>>=
357 typedef struct environment_struct {
367 #define DOUBLE_PRECISION 1e-12
369 <<environment definition>>
370 <<one dimensional function definition>>
371 <<create/destroy definitions>>
373 #endif /* GLOBAL_H */
375 Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
376 included multiple times.
378 \section{Simulation functions}
380 <<simulation functions>>=
384 <<does domain transition>>
385 <<randomly transition domains>>
389 \label{sec.find_tension}
391 Because the stretched system may be made up of several parts (folded
392 domains, unfolded domains, spring-like cantilever, \ldots), we will
393 assign the domains to states with tension models and groups. The
394 models are used to determine the tension of all the domain in a given
395 state following some model (Hookean spring, WLC, \ldots). The domains
396 are assumed to be commutative, so ordering is ignored. The
397 interactions between the groups are assumed to be linear, but the
398 interactions between domains of the same group need not be. Each
399 model has a tension handler function, which gives the tension $F$
400 needed to stretch a given group of domains a distance $x$.
402 Groups represent collections of states which the model should treat as
403 a single entity. To understand the purpose of groups, consider a
404 system with two unfolded domain states (e.g.\ for two different
405 proteins) that were both modeled as WLCs. If both unfolded states had
406 the same persistence length, it would be wise to assign them to the
407 same group. This would reduce both the computational cost of
408 balancing the tension and the error associated with the inter-group
409 interaction linearity. Note that group numbers only matter
410 \emph{within} model classes, since grouping domains with seperate
411 models doesn't make sense. Within designated groups, the tension
412 parameters for each domain are still checked for compatibility, so if
413 you accidentally assign incompatible domains to the same group you
414 should get an appropriate error message.
416 <<find tension definitions>>=
417 #define NUM_TENSION_MODELS 5
421 <<tension handler array global declaration>>=
422 extern one_dim_fn_t *tension_handlers[];
425 <<tension handler array global>>=
426 one_dim_fn_t *tension_handlers[] = {
428 &const_tension_handler,
436 <<tension model global declarations>>=
437 <<tension handler array global declaration>>
440 <<tension handler types>>=
441 typedef struct tension_handler_data_struct {
442 list_t *group_tension_params; /* one item for each domain in the group */
445 } tension_handler_data_t;
449 After sorting the chain into separate groups $G_i$, with tension
450 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
451 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
452 \forall i,j$. Note that there may be several states within each
453 group. [[find_tension]] is basically a wrapper to feed proper input
454 into the [[tension_balance]] function. [[transition]] should be set
455 to 0 if there were no transitions since the previous call to
456 [[find_tension]] to support some optimizations that assume a small
457 increase in tension between steps (only true if no transition has
458 occured). While we're messing with the tension handlers, we also grab
459 a measure of the current linker stiffness for the environmental
460 variable, which is needed by some of the unfolding rate models
461 (e.g.\ the linker-stiffness corrected Bell model, Appendix
464 double find_tension(list_t *state_list, environment_t *env, double x, int transition)
466 static int num_active_groups=0;
467 static one_dim_fn_t **per_group_handlers = NULL;
468 static one_dim_fn_t **per_group_inverse_handlers = NULL;
469 static void **per_group_data = NULL; /* array of pointers to tension_handler_data_t */
470 static double last_x;
471 tension_handler_data_t *tdata;
476 fprintf(stderr, "%s:\tfinding tension at distance %g%s\n", __FUNCTION__, x, transition ? " after a transition" : "");
479 if (transition != 0 || num_active_groups == 0) { /* setup statics */
480 /* get new statics, freeing the old ones if they aren't NULL */
481 get_active_groups(state_list, &num_active_groups, &per_group_handlers, &per_group_inverse_handlers, &per_group_data);
483 /* fill in the group handlers and params */
485 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]);
487 for (i=0; i<num_active_groups; i++) {
488 tdata = (tension_handler_data_t *) per_group_data[i];
490 fprintf(stderr, "%s:\tactive group %d of %d (data %p, param list %p)", __FUNCTION__, i, num_active_groups, tdata, tdata->group_tension_params);
491 list_print(stderr, tdata->group_tension_params, " parameters");
494 /* tdata->persist continues unchanged */
497 } /* else continue with the current statics */
499 F = tension_balance(num_active_groups, per_group_handlers, per_group_inverse_handlers, (void **)per_group_data, last_x, x, &env->stiffness);
505 @ For the moment, we will restrict our group boundaries to a single
506 dimension, so $\sum_i x_i = x$, our total extension (it is this
507 restriction that causes the groups to interact linearly). We'll also
508 restrict our extensions to all be positive. With these restrictions,
509 the problem of balancing the tensions reduces to solving a system of
510 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
511 the number of active groups. In general this can be fairly
512 complicated, but without much loss of practicality we can restrict
513 ourselves to strictly monotonically increasing, non-negative tension
514 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
515 simpler. For example, it guarantees the existence of a unique, real
516 solution for finite forces. See Appendix \ref{sec.tension_balance}
517 for the tension balancing implementation.
519 Each group must define a parameter structure if the tension model is
521 a function to create the parameter structure,
522 a function to destroy the parameter structure, and
523 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as its data argument.
524 The tension-specific parameter structures are contained in the domain
525 group field. For optimization, a group may also define an inverse
526 tension function as an optimization. See the Section
527 \ref{sec.model_selection} for a discussion on the command line
528 framework and inverse function discussion. See the worm-like chain
529 implementation for an example (Section \ref{sec.wlc}).
531 The major limitation of our approach is that all of our tension models
532 are Markovian (memory-less), which is not really a problem for the
533 chains (see \citet{evans99} for WLC thermalization rate calculations).
535 \subsection{Transition rate}
536 \label{sec.transition_rate}
538 Each state transition is modeled as unimolecular, first order reaction
540 \text{State 1} \xrightarrow{k} \text{State 2}
542 With the rate constant $k$ defined by
544 \frac{d[\text{State 2}]}{dt} = -\frac{d[\text{State 1}]}{dt} = k [\text{State 1}].
546 So the probability of a given domain transitioning in a short time
552 For a group of $N$ identical domains, and therefore identicalunfolding
553 probabilities $P$, the probability of $n$ domains unfolding in a given
554 timestep is given by the binomial distribution
556 p(n) = \frac{N!}{n!(N-n)!}P^n(1-P)^(N-n).
558 The probability that \emph{none} of these domains unfold is then
562 so the probability that \emph{at least one} domain unfolds is
566 Note that for small $P$,
568 p(n>0) = 1-(1-NP+\mathcal{O}(P^2)) = NP - \mathcal{O}(P^2) \approx NP.
570 We take some time to discuss the meaning of $p(n>1)$
571 (i.e. multi-unfolding timesteps) in Section \ref{sec.timescales}.
573 <<does domain transition>>=
574 int domain_transitions(double F, double dt, environment_t *env, int num_domains, transition_t *transition)
575 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, number of 'reactant' domains, pointer to transition structure */
577 k = accel_k(transition->k, F, env, transition->k_params);
578 //(*transition->k)(F, env, domain->k_params);
579 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
580 return happens(1-pow((1.0-k*dt), num_domains)); /* N dice rolls for prob. k*dt event */
582 @ [[happens]] is a random decision making function defined in Appendix \ref{sec.utils}.
584 <<randomly transition domains>>=
585 int random_transitions(list_t *transition_list, double tension,
586 double dt, environment_t *env)
587 { /* returns 1 if there was a transition and 0 otherwise */
588 while (transition_list != NULL) {
589 if (T(transition_list)->initial_state->num_domains > 0
590 && domain_transitions(tension, dt, env, T(transition_list)->initial_state->num_domains, T(transition_list))) {
591 if (flags & FLAG_OUTPUT_TRANSITION_FORCES)
592 fprintf(stdout, "%g\t%s\t%s\n", tension, T(transition_list)->initial_state->name, T(transition_list)->final_state->name);
593 T(transition_list)->initial_state->num_domains--;
594 T(transition_list)->final_state->num_domains++;
595 return 1; /* our one domain has transitioned, stop looking for others */
597 transition_list = transition_list->next;
603 \subsection{Timescales}
604 \label{sec.timescales}
606 The simulation assumes that chain equilibration is instantaneous
607 relative to the simulation timestep $dt$, so that tension is uniform
608 along the chain. The quality of this assumption depends on your
609 particular chain. For example, a damped spring thermalizes on a
610 timescale of order $\tau = 1/\gamma$ where the damping ratio $\gamma
611 \equiv \omega_0(-\zeta \pm \sqrt{\zeta^2-1}$, the natural angular
612 frequency $omega_0 \equiv \sqrt{k/m}$, $\zeta \equiv b/2\sqrt{km}$,
613 and $b$ sets the drag $F_b = -bv$. For our cantilevers $\tau$ is
614 on the order of milliseconds, which is longer than a timestep.
615 However, the approximation is still reasonable, because there is
616 rarely unfolding during the cantilever return. You could, of course,
617 take cantilever, WLC, etc.\ response times into effect, but that
618 would significantly complicate a simulation that seems to work
619 well enough without bothering :p. For WLC and FJC relaxation times,
620 see the Appendix of \citet{evans99} and \citet{kroy07}.
622 Besides assuming our timestep is much greater than equilibration
623 times, we also force it to be short enough so that only one domain may
624 unfold in a given timestep. For an unfolding timescale of $dt_u =
625 1/k$, we adapt our timesteps to keep $dt \ll dt_u$, so the probability
626 of two domains unfolding in the same timestep is small (see
627 [[determine_dt]] in Section \ref{sec.adaptive_dt} and [[P]] in
628 [[main()]] in Section \ref{sec.main} set by [[-P]] in
629 [[get_options()]] in Section \ref{sec.get_options}). This approach
630 breaks down as the adaptive timestep scheme approaches the transition
631 timestep $dt_t$, but $dt_t \sim 1\U{$\mu$s}$ for Ig27-like proteins
632 \citep{klimov00}, so this shouldn't be a problem. To reassure
633 yourself, you can ask the simulator to print the smallest timestep
634 that was used with TODO.
636 Even if the likelyhood of two domains unfolding in the same timestep
637 is small, if you run long enough simulations it will eventually occur.
638 In this case, we assume that $dt_t \ll dt$, so even if two domains
639 would unfold if we stuck to our unfolding rate $k$ for an entire
640 timestep $dt$, in ``reality'' only one of those domains unfolds.
641 Because we do not know when in the timestep the unfolding took place,
642 we move on to the next timestep without checking for further
643 unfoldings. This ``unchecked timestep portion'' should not contribute
644 significant errors to the calculation, because if $dt$ was small
645 enough that unfolding was unlikely at the pre-unfolding force, it
646 should be even more unlikely at the post-unfolding force, especially
647 over only a fraction of the timestep. The only dangerous cases for
648 the current implementation are those where unfolding intermediates are
649 much less stable than their precursor states, in which case an
650 unfolding event during the remaining timestep may actually be likely.
651 A simple workaround for such cases is to pick the value for [[P]] to
652 be small enough that the $dt \ll dt_u$ assumption survives even under
653 a transition populating the unstable intermediate.
655 \subsection{Adaptive timesteps}
656 \label{sec.adaptive_dt}
658 We'd like to pick $dt$ so the probability of changing state
659 (e.g.\ unfolding) in the next timestep is small. If we don't adapt our
660 timestep, we also risk breaking our approximation that $F(x' \in [x,
661 x+v \cdot dt]) = F(x)$ or that only one domain changes state in a
662 given timestep. The simulation would have been accurate for
663 sufficiently small fixed timesteps, but adaptive timesteps will allow
664 us to move through low tension regions in fewer steps, leading to a
665 more efficient simulation.
667 Consider the two state folded $\rightarrow$ unfolded transition.
668 Because $F(x)$ is monotonically increasing between unfoldings,
669 excessively large timesteps will lead to erroneously small unfolding
670 rates (an thus, higher average unfolding force).
672 The actual adaptive timestep implementation is not particularly
673 interesting, since we are only required to reduce $dt$ to somewhere
674 below a set threshold, so I've removed it to Appendix
675 \ref{sec.adaptive_dt_implementation}.
681 The models provide the physics of the simulation, but the simulation
682 allows interchangeable models, and we are currently focusing on the
683 simulation itself, so we remove the actual model implementation to the
686 The tension models are defined in Appendix \ref{sec.tension_models}
687 and unfolding models are defined in Appendix \ref{sec.k_models}.
690 #define k_B 1.3806503e-23 /* J/K */
694 \section{Command line interface}
696 <<get option functions>>=
698 <<index tension model>>
703 \subsection{Model selection}
704 \label{sec.model_selection}
706 The main difficulty with the command line interface in \prog\ is
707 developing an intuitive method for accessing the possibly large number
708 of available models. We'll treat this generally by defining an array
709 of available models, containing their important parameters.
711 <<get options definitions>>=
712 <<model getopt definitions>>
715 <<create/destroy definitions>>=
716 typedef void *create_data_func_t(char **param_strings);
717 typedef void destroy_data_func_t(void *param_struct);
720 <<model getopt definitions>>=
721 <<tension model getopt definitions>>
722 <<k model getopt definitions>>
725 In order to choose models, we need a method of assembling a reaction
726 graph such as that in Figure \ref{fig.domain_states} and population
727 the graph with a set of domains. First we declare the domain states
730 [[-s <state-name>,<model>[,<group>] -S <model-params> [-N <initial-population]]] \\
734 [[-s unfolded,wlc,1 -S 1e-8,4e-10 -N 5]]
736 This creates the state named [[unfolded]], using the WLC model and one
737 as a group number (defaults to 0, see Section \ref{sec.find_tension}).
738 The tension parameters are then set to [[1e-8,4e-10]], which in the
739 case of the WLC are the contour and persistence lengths respectively.
740 Finally we populate the state by adding five domains to the state.
741 The name of the state is importent for identifying states later.
743 Once the domains have been created and populated, you can added
744 transition rates following
746 [[-k <initial-state>,<final-state>,<model> -K <model-params>]] \\
750 [[-k folded,unfolded,bell -K 3.3e-4,0.25e-9]]
752 This creates a reaction going from the [[folded]] state to the
753 [[unfolded]] state, with the rate constant given by the Bell model
754 (Appendix \ref{sec.bell}). The unfolding parameters are then set to
755 [[3.3e-4,0.25e-9]], which in the case of the Bell model are the
756 unforced rate and transition state distance respectively.
758 \subsubsection{Tension}
760 To access the various tension models from the command line, we define
761 a structure that contains all the neccessary information about the
763 <<tension model getopt definitions>>=
764 typedef struct tension_model_getopt_struct {
765 one_dim_fn_t *handler; /* fn implementing the model on a group */
766 one_dim_fn_t *inverse_handler; /* fn implementing inverse on group */
767 char *name; /* name string identifying the model */
768 char *description; /* a brief description of the model */
769 int num_params; /* number of extra params for the model */
770 char **param_descriptions; /* descriptions of extra params */
771 char *params; /* default values for the extra params */
772 create_data_func_t *creator; /* param-string -> model-data-struct fn */
773 destroy_data_func_t *destructor; /* fn to free the model data structure */
774 } tension_model_getopt_t;
775 @ Set [[inverse_handler]] to [[NULL]] if you wish to invert using
776 bisection. Obviously, this will be slower than computing an
777 analytical inverse and more prone to rounding errors, but it may be
778 easier if a closed-form inverse does not exist.
780 <<tension model getopt array definition>>=
781 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
782 <<null tension model getopt>>,
783 <<constant tension model getopt>>,
784 <<hooke tension model getopt>>,
785 <<worm-like chain tension model getopt>>,
786 <<freely-jointed chain tension model getopt>>
790 \subsubsection{Unfolding rate}
792 <<k model getopt definitions>>=
793 #define NUM_K_MODELS 6
795 typedef struct k_model_getopt_struct {
800 char **param_descriptions;
802 create_data_func_t *creator;
803 destroy_data_func_t *destructor;
807 <<k model getopt array definition>>=
808 k_model_getopt_t k_models[NUM_K_MODELS] = {
809 <<null k model getopt>>,
810 <<const k model getopt>>,
811 <<bell k model getopt>>,
812 <<kbell k model getopt>>,
813 <<kramers k model getopt>>,
814 <<kramers integ k model getopt>>
821 void help(char *prog_name, double P, double t_max, double dt_max, double v, double x_step,
822 state_t *stop_state, environment_t *env,
823 int n_tension_models, tension_model_getopt_t *tension_models,
824 int n_k_models, k_model_getopt_t *k_models,
825 int k_model, list_t *state_list)
830 if (state_list != NULL)
831 state = S(tail(state_list));
833 printf("usage: %s [options]\n", prog_name);
834 printf("Version %s\n\n", VERSION);
835 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
837 printf("Simulation options:\n");
839 printf("-q state\tStop simulation when this state empties (- for no limit, currently %s)\n", stop_state ? stop_state->name : "-" );
840 printf("-t tmax\tMaximum simulated time (-1 for no limit, currently %g s)\n", t_max);
841 printf("-d dt\tMaximum allowed timestep dt (currently %g s)\n", dt_max);
842 printf("-P P\tTarget probability for dt (currently %g)\n", P);
843 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
844 printf("-x dx\tPulling step size (currently %g m)\n", x_step);
845 printf("\tWhen dx = 0, the pulling is continuous.\n");
846 printf("\tWhen dx > 0, the pulling occurs in discrete steps.\n");
848 printf("Environmental options:\n");
849 printf("-T T\tTemperature T (currently %g K)\n", env->T);
850 printf("-C T\tYou can also set the temperature T in Celsius\n");
852 printf("State creation options:\n");
853 printf("-s name,model[[,group],params]\tCreate a new state.\n");
854 printf("Once you've set up the state, you may populate it with domains\n");
855 printf("-N n\tAdd n domains in the current state (%s)\n", state ? state->name : "-");
857 printf("Transition creation options:\n");
858 printf("-k initial_state,final_state,model[,params]\tCreate a new transition\n");
860 printf("To double check your reaction graph, you can produce graphviz dot output\n");
861 printf("describing the current states and transitions.\n");
862 printf("-D\tPrint dot description of states and transitions and exit.\n");
864 printf("Simulation output mode options:\n");
865 printf("There are two output modes. In standard mode, only the transition\n");
866 printf("events are printed. For example:\n");
867 printf(" #Force (N)\tInitial state\tFinal state\n");
868 printf(" 123.456e-12\tfolded\tunfolded\n");
870 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
871 printf(" #Distance (m)\tForce (N)\tStiffness (N/m)\n");
872 printf(" 0.001\t0.0023\t0.002\n");
874 printf("-V\tChange output to verbose mode.\n");
875 printf("-h\tPrint this help and exit.\n");
878 printf("Tension models:\n");
879 for (i=0; i<n_tension_models; i++) {
880 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
881 for (j=0; j < tension_models[i].num_params ; j++)
882 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
883 printf("\t\tdefaults: %s\n", tension_models[i].params);
885 printf("Unfolding rate models:\n");
886 for (i=0; i<n_k_models; i++) {
887 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
888 for (j=0; j < k_models[i].num_params ; j++)
889 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
890 printf("\t\tdefaults: %s\n", k_models[i].params);
893 printf("Examples:\n");
894 printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded. Stop once there are no more folded domains.\n");
895 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);
896 printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded. Stop after 0.25 seconds\n");
897 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);
901 \subsection{Get options}
902 \label{sec.get_options}
906 void get_options(int argc, char **argv,
907 double *pP, double *pT_max, double *pDt_max, double *pV,
908 double *pX_step, state_t **pStop_state, environment_t *env,
909 int n_tension_models, tension_model_getopt_t *tension_models,
910 int n_k_models, k_model_getopt_t *k_models,
911 list_t **pState_list, list_t **pTransition_list)
913 char *prog_name = NULL;
914 char c, options[] = "q:t:d:P:v:x:T:C:s:N:k:DVh";
915 int tension_model=0, k_model=0, initial_state=-1, final_state=-1;
916 char *state_name=NULL;
917 int i, n, tension_group=0;
918 list_t *temp_list=NULL;
921 void *model_params; /* for INIT_MODEL() */
922 int num_param_args; /* for INIT_MODEL() */
923 char **param_args; /* for INIT_MODEL() */
925 extern int optind, optopt, opterr;
930 for (i=0; i<n_tension_models; i++) {
931 fprintf(stderr, "tension model %i %s:\t model's handler %p\t inverse handler %p\n",
932 i, tension_models[i].name, tension_models[i].handler, tension_models[i].inverse_handler);
933 assert(tension_models[i].handler == tension_handlers[i]);
938 flags = FLAG_OUTPUT_TRANSITION_FORCES;
941 *pT_max = -1; /* s, -1 removes this limit */
942 *pDt_max = 0.001; /* s */
943 *pP = 1e-3; /* % pop per s */
944 *pV = 1e-6; /* m/s */
945 *pX_step = 0.0; /* m */
946 env->T = 300.0; /* K */
948 *pTransition_list = NULL;
950 while ((c=getopt(argc, argv, options)) != -1) {
953 if (STRMATCH(optarg, "-")) {
956 temp_list = head(*pState_list);
957 while (temp_list != NULL) {
958 if (STRMATCH(S(temp_list)->name, optarg)) {
959 *pStop_state = S(temp_list);
962 temp_list = temp_list->next;
966 case 't': *pT_max = safe_strtod(optarg, "-t"); break;
967 case 'd': *pDt_max = safe_strtod(optarg, "-d"); break;
968 case 'P': *pP = safe_strtod(optarg, "-P"); break;
969 case 'v': *pV = safe_strtod(optarg, "-v"); break;
970 case 'x': *pX_step = safe_strtod(optarg, "-x"); break;
971 case 'T': env->T = safe_strtod(optarg, "-T"); break;
972 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
974 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
975 assert(num_strings >= 2 && num_strings <= 4);
977 state_name = strings[0];
978 if (state_index(*pState_list, state_name) != -1) {
979 fprintf(stderr, "Invalid option '%s', repeated state name '%s'\n", optarg, state_name);
982 tension_model = index_tension_model(n_tension_models, tension_models, strings[1]);
983 if (num_strings == 4)
984 tension_group = safe_strtoi(strings[2], optarg);
987 if (num_strings >= 3)
988 tension_models[tension_model].params = strings[num_strings-1];
990 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);
992 INIT_MODEL(state_name, &tension_models[tension_model], tension_models[tension_model].params, model_params);
993 push(pState_list, create_state(state_name,
994 tension_models[tension_model].handler,
995 tension_models[tension_model].inverse_handler,
998 tension_models[tension_model].destructor,
1000 *pState_list = tail(*pState_list);
1002 free_parsed_list(num_strings, strings);
1005 n = safe_strtoi(optarg, "-N");
1007 fprintf(stderr, "%s:\tadding %d domains to %s\n", __FUNCTION__, n, S(*pState_list)->name);
1009 S(*pState_list)->num_domains += n;
1012 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
1013 assert(num_strings >= 3 && num_strings <= 4);
1014 initial_state = state_index(*pState_list, strings[0]);
1015 final_state = state_index(*pState_list, strings[1]);
1016 k_model = index_k_model(n_k_models, k_models, strings[2]);
1018 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);
1020 assert(initial_state != final_state);
1021 if (num_strings == 4)
1022 k_models[k_model].params = strings[3];
1023 INIT_MODEL(state_name, &k_models[k_model], k_models[k_model].params, model_params);
1024 push(pTransition_list, create_transition(list_element(*pState_list, initial_state),
1025 list_element(*pState_list, final_state),
1026 k_models[k_model].k,
1028 k_models[k_model].destructor), 1);
1029 free_parsed_list(num_strings, strings);
1032 print_state_graph(stdout, *pState_list, *pTransition_list);
1036 flags = FLAG_OUTPUT_FULL_CURVE;
1039 fprintf(stderr, "unrecognized option '%c'\n", optopt);
1040 /* fall through to default case */
1042 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);
1046 /* check the arguments */
1047 assert(*pState_list != NULL); /* no domains! */
1048 assert(*pP > 0.0); assert(*pP < 1.0);
1049 assert(*pDt_max > 0.0);
1051 assert(env->T > 0.0);
1053 crosslink_groups(*pState_list, *pTransition_list);
1059 <<index tension model>>=
1060 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
1063 for (i=0; i<n_models; i++)
1064 if (STRMATCH(models[i].name, name))
1066 if (i == n_models) {
1067 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
1075 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
1078 for (i=0; i<n_models; i++)
1079 if (STRMATCH(models[i].name, name))
1081 if (i == n_models) {
1082 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
1089 <<initialize model definition>>=
1090 /* requires int num_param_args and char **param_args in the current scope
1092 * INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
1093 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
1095 #define INIT_MODEL(role, model, param_string, param_pointer) \
1097 parse_list_string((param_string), SEP, '{', '}', \
1098 &num_param_args, ¶m_args); \
1099 if (num_param_args != (model)->num_params) { \
1101 "%s model %s expected %d params,\n", \
1102 (role), (model)->name, (model)->num_params); \
1104 "not the %d params in '%s'\n", \
1105 num_param_args, (param_string)); \
1106 assert(num_param_args == (model)->num_params); \
1108 if ((model)->creator) \
1109 param_pointer = (*(model)->creator)(param_args); \
1111 param_pointer = NULL; \
1112 free_parsed_list(num_param_args, param_args); \
1116 First we define some safe string parsing functions. Currently these
1117 abort on any irregularity, but in the future they could print error
1118 messages, etc. [[static]] because the functions are currently
1119 defined in each [[*.c]] file that uses them, but perhaps they should
1120 be moved to [[utils.h/c]] or some such instead\ldots
1123 static int safe_strtoi(const char *s, const char *description) {
1127 ret = strtol(s, &endp, 10);
1128 if (endp[0] != '\0') { //strlen(endp) > 0) {
1129 fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1130 endp, s, description);
1131 assert(1==0); //strlen(endp) == 0);
1133 fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1136 } else if (errno == ERANGE) {
1137 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1138 assert(errno != ERANGE);
1143 static double safe_strtod(const char *s, const char *description) {
1147 ret = strtod(s, &endp);
1148 if (endp[0] != '\0') { //strlen(endp) > 0) {
1149 fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1150 endp, s, description);
1151 assert(1==0); //strlen(endp) == 0);
1153 fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1156 } else if (errno == ERANGE) {
1157 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1158 assert(errno != ERANGE);
1167 \addcontentsline{toc}{section}{Appendicies}
1169 \section{sawsim.c details}
1173 The general layout of our simulation code is:
1184 We include [[math.h]], so don't forget to link to the libm with `-lm'.
1186 #include <assert.h> /* assert() */
1187 #include <stdlib.h> /* malloc(), free(), rand(), strto*() */
1188 #include <stdio.h> /* fprintf(), stdout */
1189 #include <string.h> /* strlen, strtok() */
1190 #include <math.h> /* exp(), M_PI, sqrt() */
1191 #include <sys/types.h> /* pid_t (returned by getpid()) */
1192 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
1193 #include <errno.h> /* errno, ERANGE (for safe_strto*()) */
1196 #include "tension_balance.h"
1197 #include "k_model.h"
1198 #include "tension_model.h"
1200 #include "accel_k.h"
1204 <<version definition>>
1205 <<flag definitions>>
1206 <<max/min definitions>>
1207 <<string comparison definition and globals>>
1208 <<initialize model definition>>
1209 <<get options definitions>>
1210 <<state definitions>>
1211 <<transition definitions>>
1220 <<create/destroy state>>
1221 <<destroy state list>>
1223 <<create/destroy transition>>
1224 <<destroy transition list>>
1225 <<print state graph>>
1227 <<simulation functions>>
1228 <<boilerplate functions>>
1231 <<boilerplate functions>>=
1233 <<get option functions>>
1239 srand(getpid()*time(NULL)); /* seed rand() */
1243 <<flag definitions>>=
1244 /* in octal b/c of prefixed '0' */
1245 #define FLAG_OUTPUT_FULL_CURVE 01
1246 #define FLAG_OUTPUT_TRANSITION_FORCES 02
1250 static unsigned long int flags = 0;
1253 \subsection{Utilities}
1256 <<max/min definitions>>=
1257 #define MAX(a,b) ((a)>(b) ? (a) : (b))
1258 #define MIN(a,b) ((a)<(b) ? (a) : (b))
1261 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
1262 <<string comparison definition and globals>>=
1263 // Check if two strings match, return 1 if they do
1264 static char *temp_string_A;
1265 static char *temp_string_B;
1266 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
1267 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
1268 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
1269 /* +1 to also compare the '\0' */
1272 We also define a macro for our [[check]] unit testing
1273 <<check relative error>>=
1274 #define CHECK_ERR(max_err, expected, received) \
1276 fail_unless( (received-expected)/expected < max_err, \
1277 "relative error %g >= %g in %s (Expected %g, received %g)", \
1278 (received-expected)/expected, max_err, #received, \
1279 expected, received); \
1280 fail_unless(-(received-expected)/expected < max_err, \
1281 "relative error %g >= %g in %s (Expected %g, received %g)", \
1282 -(received-expected)/expected, max_err, #received, \
1283 expected, received); \
1288 int happens(double probability)
1290 assert(probability >= 0.0); assert(probability <= 1.0);
1291 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*/
1296 \subsection{Adaptive timesteps}
1297 \label{sec.adaptive_dt_implementation}
1299 $F(x)$ increases with $x$, possibly exploding (e.g.\ in the worm-like
1300 chain model), so basing the timestep on the the unfolding probability
1301 at the current tension is dangerous, and we need to search for a $dt$
1302 for which $P(F(x+v*dt)) < P_\text{target}$. There are two cases to
1303 consider. In the most common, no domains have unfolded since the last
1304 step, and we expect the next step to be slightly shorter than the
1305 previous one. In the less common, domains did unfold in the last
1306 step, and we expect the next step to be considerably longer than the
1309 double search_dt(transition_t *transition,
1311 environment_t *env, double x,
1312 double target_prob, double max_dt, double v)
1313 { /* Returns a valid timestep dt in seconds for the given transition.
1314 * Takes a list of all states, a pointer env to the environmental data,
1315 * a starting separation x in m, a target_prob between 0 and 1,
1316 * max_dt in s, stretching velocity v in m/s.
1318 double F, k, dtCur, dtU, dtUCur, dtL, dt;
1320 /* get upper bound using the current position */
1321 F = find_tension(state_list, env, x, 0); /* BUG. repeated calculation */
1322 //printf("Start with x = %g (F = %g)\n", x, F);
1323 k = accel_k(transition->k, F, env, transition->k_params);
1324 //printf("x %g\tF %g\tk %g\n", x, F, k);
1325 dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */
1327 //printf("overshot max_dt\n");
1330 if (v == 0) /* discrete stepping, so force is temporarily constant */
1333 /* set a lower bound on dt too */
1336 /* The dt determined above may produce illegitimate forces or ks.
1337 * Reduce the upper bound until we have valid ks. */
1339 F = find_tension(state_list, env, x+v*dt, 0);
1340 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1343 F = find_tension(state_list, env, x+v*dt, 0);
1345 //printf("Try for dt = %g (F = %g)\n", dt, F);
1346 k = accel_k(transition->k, F, env, transition->k_params);
1347 /* returned k may be -1.0 */
1348 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1349 while (k == -1.0) { /* reduce step until we hit a valid k */
1351 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1352 F = find_tension(state_list, env, x+v*dt, 0);
1353 //printf("Try for dt = %g (F = %g)\n", dt, F);
1354 k = accel_k(transition->k, F, env, transition->k_params);
1355 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1357 assert(dtU > 1e-14); /* timestep to valid k too small */
1358 dtUCur = target_prob / k; /* safe timestep back from x+dtU */
1360 return dt; /* dtU is safe. */
1362 /* dtU wasn't safe, lets see what would be. */
1363 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1364 dt = (dtU + dtL) / 2.0;
1365 F = find_tension(state_list, env, x+v*dt, 0);
1366 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1367 k = accel_k(transition->k, F, env, transition->k_params);
1368 dtCur = target_prob / k;
1369 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1370 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1372 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1373 dtU = dt; /* ... stepping out only dtCur would be safe */
1376 break; /* dtCur = dt */
1378 return MAX(dtUCur, dtL);
1383 To determine $dt$ for an array of potentially different folded domains,
1384 we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
1387 double determine_dt(list_t *state_list, list_t *transition_list,
1388 environment_t *env, double x,
1389 double target_prob, double dt_max, double v)
1390 { /* Returns the timestep dt in seconds.
1391 * Takes the state and transition lists, a pointer to the environment,
1392 * an extension x in m, target_prob between 0 and 1, dt_max in s,
1394 double dt=dt_max, new_dt;
1395 assert(target_prob > 0.0); assert(target_prob < 1.0);
1396 assert(dt_max > 0.0);
1397 assert(state_list != NULL);
1398 assert(transition_list != NULL);
1399 transition_list = head(transition_list);
1400 while (transition_list != NULL) {
1401 if (T(transition_list)->initial_state->num_domains > 0) {
1402 new_dt = search_dt(T(transition_list), state_list, env, x, target_prob, dt, v);
1403 dt = MIN(dt, new_dt);
1405 transition_list = transition_list->next;
1412 \subsection{State data}
1413 \label{sec.state_data}
1415 The domains exist in one of an array of possible states. Each state
1416 has a tension model which determines it's force/extension curve
1417 (possibly [[null]] for rigid states, see Appendix
1418 \ref{sec.null_tension}).
1420 Groups are, for example, for two domain states with WLC tensions; one
1421 with $p=0.4\U{nm}$, $L=10\U{nm}$ and the other with $p=0.4\U{nm}$,
1422 $L=28\U{nm}$. Since the persistence lengths are the same, you would
1423 like to add the contour lengths of all the domains in \emph{both}
1424 states, and plug that total contour length into a single WLC
1425 evaluation (see Section \ref{sec.find_tension}).
1426 <<state definitions>>=
1427 typedef struct state_struct {
1428 char *name; /* name string */
1429 one_dim_fn_t *tension_handler; /* tension handler */
1430 one_dim_fn_t *inverse_tension_handler; /* inverse tension handler */
1431 int tension_group; /* for grouping similar tension states */
1432 void *tension_params; /* pointer to tension parameters */
1433 destroy_data_func_t *destroy_tension;
1434 int num_domains; /* number of domains currently in state */
1435 /* cached lookups generated from the above data */
1436 int num_out_trans; /* length of out_trans array */
1437 struct transition_struct **out_trans; /* transitions leaving this state */
1438 int num_grouped_states; /* length of grouped states array */
1439 struct state_struct **grouped_states; /* states grouped with this one (doesn't include self) */
1442 /* get the state data for the current list node */
1443 #define S(list) ((state_t *)(list)->d)
1445 @ [[tension_params]] is a pointer to the parameters used by the
1446 handler function when determining the tension. [[destroy_tension]]
1447 points to a function for freeing [[tension_params]]. We it in the
1448 state structure so that [[destroy_state]] and will have an easier job
1451 [[create_]] and [[destroy_state]] are simple wrappers around
1452 [[malloc]] and [[free]].
1453 <<create/destroy state>>=
1454 state_t *create_state(char *name,
1455 one_dim_fn_t *tension_handler,
1456 one_dim_fn_t *inverse_tension_handler,
1458 void *tension_params,
1459 destroy_data_func_t *destroy_tension,
1462 state_t *ret = (state_t *)malloc(sizeof(state_t));
1463 assert(ret != NULL);
1464 /* make a copy of the name, so we aren't dependent on the original */
1465 ret->name = (char *)malloc(sizeof(char)*(strlen(name)+1));
1466 assert(ret->name != NULL);
1467 strcpy(ret->name, name); /* we know ret->name is long enough */
1469 ret->tension_handler = tension_handler;
1470 ret->inverse_tension_handler = inverse_tension_handler;
1471 ret->tension_group = tension_group;
1472 ret->tension_params = tension_params;
1473 ret->destroy_tension = destroy_tension;
1474 ret->num_domains = num_domains;
1476 ret->num_out_trans = 0;
1477 ret->out_trans = NULL;
1478 ret->num_grouped_states = 0;
1479 ret->grouped_states = NULL;
1482 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);
1487 void destroy_state(state_t *state)
1491 fprintf(stderr, "freeing state %s (%p)\n", state->name, state);
1495 if (state->destroy_tension)
1496 (*state->destroy_tension)(state->tension_params);
1497 if (state->out_trans)
1498 free(state->out_trans);
1499 if (state->grouped_states)
1500 free(state->grouped_states);
1507 We'll be storing the states in a list (see Appendix \ref{sec.lists}),
1508 so we also define a convenience function for cleaning up.
1509 <<destroy state list>>=
1510 void destroy_state_list(list_t *state_list)
1512 state_list = head(state_list);
1513 while (state_list != NULL)
1514 destroy_state((state_t *) pop(&state_list));
1519 We can also define a convenient way to get a state index from it's name.
1521 int state_index(list_t *state_list, char *name)
1524 state_list = head(state_list);
1525 while (state_list != NULL) {
1526 if (STRMATCH(S(state_list)->name, name)) return i;
1527 state_list = state_list->next;
1535 \subsection{Transition data}
1536 \label{sec.transition_data}
1538 Once you've created states (Sections \ref{sec.main},
1539 \ref{sec.model_selection}, and \ref{sec.state_data}), you need to
1540 create transitions between them.
1541 <<transition definitions>>=
1542 typedef struct transition_struct {
1543 state_t *initial_state;
1544 state_t *final_state;
1545 k_func_t *k; /* transition rate model */
1546 void *k_params; /* pointer to k parameters */
1547 destroy_data_func_t *destroy_k;
1550 /* get the transition data for the current list node */
1551 #define T(list) ((transition_t *)(list)->d)
1552 @ [[k]] is a pointer to the function determining the transition rate
1553 for a given tension. [[k_params]] and [[destroy_k]] have similar
1554 roles with regards to [[k]] as [[tension_params]] and
1555 [[destroy_tension]] have with regards to [[tension_handler]] (see
1556 Section \ref{sec.state_data}).
1558 [[create_]] and [[destroy_transition]] are simple wrappers around
1559 [[malloc]] and [[free]].
1560 <<create/destroy transition>>=
1561 transition_t *create_transition(state_t *initial_state, state_t *final_state,
1564 destroy_data_func_t *destroy_k)
1566 transition_t *ret = (transition_t *)malloc(sizeof(transition_t));
1567 assert(ret != NULL);
1568 assert(initial_state != NULL);
1569 assert(final_state != NULL);
1570 ret->initial_state = initial_state;
1571 ret->final_state = final_state;
1573 ret->k_params = k_params;
1574 ret->destroy_k = destroy_k;
1576 fprintf(stderr, "generated transition %p from %s to %s, params %p\n", ret, ret->initial_state->name, ret->final_state->name, ret->k_params);
1581 void destroy_transition(transition_t *transition)
1585 fprintf(stderr, "freeing transition %p (%s to %s)\n", transition, transition->initial_state->name, transition->final_state->name);
1587 if (transition->destroy_k)
1588 (*transition->destroy_k)(transition->k_params);
1595 We'll be storing the transitions in a list (see Appendix
1596 \ref{sec.lists}), so we also define a convenience function for
1598 <<destroy transition list>>=
1599 void destroy_transition_list(list_t *transition_list)
1601 transition_list = head(transition_list);
1602 while (transition_list != NULL)
1603 destroy_transition((transition_t *) pop(&transition_list));
1608 \subsection{Graphviz output}
1609 \label{sec.graphviz_output}
1611 <<print state graph>>=
1612 void print_state_graph(FILE *file, list_t *state_list, list_t *transition_list)
1614 state_list = head(state_list);
1615 transition_list = head(transition_list);
1616 fprintf(file, "digraph states {\n node [shape=record,color=black,fillcolor=white,style=filled];\n\n");
1618 fprintf(file, " n%d [label=\"%s\"]\n", state_index(state_list, S(state_list)->name), S(state_list)->name);
1619 if (state_list->next == NULL) break;
1620 state_list = state_list->next;
1622 fprintf(file, "\n");
1623 while (transition_list != NULL) {
1624 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));
1625 transition_list = transition_list->next;
1627 fprintf(file, "}\n");
1631 \subsection{Domain model and group handling}
1633 <<group functions>>=
1634 <<crosslink groups>>
1635 <<get active groups>>
1638 Scan through a list of states and construct an array of all the
1640 <<get active groups>>=
1641 void get_active_groups(list_t *state_list,
1642 int *pNum_active_groups,
1643 one_dim_fn_t ***pPer_group_handlers,
1644 one_dim_fn_t ***pPer_group_inverse_handlers,
1645 void ***pPer_group_data /* pointer to array of pointers to tension_handler_data_t */)
1647 void **active_groups=NULL;
1648 one_dim_fn_t *handler, *old_handler, *inverse_handler, *old_inverse_handler;
1649 one_dim_fn_t **per_group_handlers=NULL, **per_group_inverse_handlers=NULL;
1650 void **per_group_data=NULL; /* array of pointers to tension_handler_data_t */
1651 int i, j, num_domains, group, old_group, num_active_groups=0;
1652 list_t *active_states_list=NULL;
1653 tension_handler_data_t *tdata=NULL;
1656 /* get all the active groups in a list */
1657 state_list = head(state_list);
1659 fprintf(stderr, "%s:\t", __FUNCTION__);
1660 list_print(stderr, state_list, "states");
1662 while (state_list != NULL) {
1663 state = S(state_list);
1664 handler = state->tension_handler;
1665 inverse_handler = state->inverse_tension_handler;
1666 group = state->tension_group;
1667 num_domains = state->num_domains;
1668 if (list_index(active_states_list, state) == -1) {
1669 /* we haven't added this state (or it's associated group) yet */
1670 if (num_domains > 0 && handler != NULL) { /* add the state */
1671 num_active_groups++;
1672 push(&active_states_list, state, 1);
1674 fprintf(stderr, "%s:\tactive %s state (%p) with %d domain(s)\n", __FUNCTION__, state->name, state, state->num_domains);
1676 for (i=0; i < state->num_grouped_states; i++) {
1677 if (state->grouped_states[i]->num_domains > 0) {
1678 /* add this related (and active) state now */
1679 assert(state->grouped_states[i]->tension_handler == handler);
1680 assert(state->grouped_states[i]->tension_group == group);
1681 push(&active_states_list, state->grouped_states[i], 1);
1683 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);
1687 } /* else empty state or NULL handler */
1688 } /* else state already added */
1689 state_list = state_list->next;
1692 fprintf(stderr, "%s:\t", __FUNCTION__);
1693 list_print(stderr, active_states_list, "active states");
1696 assert(num_active_groups <= list_length(active_states_list));
1697 assert(num_active_groups > 0);
1699 /* allocate the arrays we'll be returning */
1700 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1701 assert(per_group_handlers != NULL);
1702 per_group_inverse_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1703 assert(per_group_inverse_handlers != NULL);
1704 per_group_data = (void **)malloc(sizeof(void *)*num_active_groups);
1705 assert(per_group_data != NULL);
1707 fprintf(stderr, "%s:\tallocating per group data %p:", __FUNCTION__, per_group_data);
1709 for (i=0; i<num_active_groups; i++) { /* we might want to recycle any group data structs allocated in the previous round */
1710 per_group_data[i] = malloc(sizeof(tension_handler_data_t));
1711 assert(per_group_data[i] != NULL);
1713 fprintf(stderr, " %d,%p", i, per_group_data[i]);
1717 fprintf(stderr, "\n");
1720 /* populate the arrays we'll be returning */
1721 active_states_list = head(active_states_list);
1723 old_inverse_handler = NULL;
1724 j = -1; /* j is the current active _group_ index */
1725 while (active_states_list != NULL) {
1726 state = (state_t *) pop(&active_states_list);
1727 handler = state->tension_handler;
1728 inverse_handler = state->inverse_tension_handler;
1729 group = state->tension_group;
1730 num_domains = state->num_domains;
1731 if (handler != old_handler || inverse_handler != old_inverse_handler || group != old_group) {
1732 j++; /* move to the next group */
1733 tdata = (tension_handler_data_t *) per_group_data[j];
1734 per_group_handlers[j] = handler;
1735 per_group_inverse_handlers[j] = inverse_handler;
1736 tdata->group_tension_params = NULL;
1738 tdata->persist = NULL; /* we might want to grab any appropriate data from the old persist here... */
1741 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);
1743 for (i=0; i<num_domains; i++) { /* add num_domains copies of the tension params */
1744 push(&tdata->group_tension_params, state->tension_params, 1);
1747 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);
1749 old_handler = handler;
1750 old_inverse_handler = inverse_handler;
1754 /* free old groups */
1755 if (*pPer_group_handlers != NULL)
1756 free(*pPer_group_handlers);
1757 if (*pPer_group_inverse_handlers != NULL)
1758 free(*pPer_group_inverse_handlers);
1759 if (*pPer_group_data != NULL) {
1760 for (i=0; i<*pNum_active_groups; i++)
1761 free((*pPer_group_data)[i]);
1762 free(*pPer_group_data);
1765 /* set new groups */
1767 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]);
1769 *pNum_active_groups = num_active_groups;
1770 *pPer_group_handlers = per_group_handlers;
1771 *pPer_group_inverse_handlers = per_group_inverse_handlers;
1772 *pPer_group_data = per_group_data;
1776 <<crosslink groups>>=
1777 void crosslink_groups(list_t *state_groups, list_t *transition_list)
1779 list_t *list, *out_trans=NULL, *associates=NULL;
1780 one_dim_fn_t *handler;
1783 state_groups = head(state_groups);
1784 while (state_groups != NULL) {
1785 /* find transitions leaving this state */
1787 fprintf(stderr, "finding transitions leaving %s\n", S(state_groups)->name);
1789 list = head(transition_list);
1790 while (list != NULL) {
1791 if (T(list)->initial_state == S(state_groups) && list_index(associates, T(list)->initial_state) != -1) {
1792 push(&out_trans, T(list), 1);
1796 n = list_length(out_trans);
1797 S(state_groups)->num_out_trans = n;
1798 S(state_groups)->out_trans = (transition_t **)malloc(sizeof(transition_t *)*n);
1799 assert(S(state_groups)->out_trans != NULL);
1800 for (i=0; i<n; i++) {
1801 S(state_groups)->out_trans[i] = (transition_t *)pop(&out_trans);
1804 /* find states grouped with this state */
1806 fprintf(stderr, "finding states grouped with %s\n", S(state_groups)->name);
1808 handler = S(state_groups)->tension_handler;
1809 group = S(state_groups)->tension_group;
1810 list = head(state_groups);
1811 while (list != NULL) {
1812 if (S(list)->tension_handler == handler && S(list)->tension_group == group && list != state_groups) {
1813 push(&associates, S(list), 1);
1817 n = list_length(associates);
1818 S(state_groups)->num_grouped_states = n;
1819 S(state_groups)->grouped_states = (state_t **)malloc(sizeof(state_t *)*n);
1820 assert(S(state_groups)->grouped_states != NULL);
1821 for (i=0; i<n; i++) {
1822 S(state_groups)->grouped_states[i] = (state_t *)pop(&associates);
1824 state_groups = state_groups->next;
1830 \section{String parsing}
1832 For handling command line arguments, we need parse delimited strings (e.g.\ [["param1,param2,param3..."]]).
1833 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1834 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1835 need to take care of parsing those parameters themselves.
1836 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1842 <<parse definitions>>
1843 <<parse declarations>>
1847 <<parse module makefile lines>>=
1848 NW_SAWSIM_MODS += parse
1849 CHECK_BINS += check_parse
1853 <<parse definitions>>=
1854 #define SEP ',' /* argument separator character */
1857 <<parse declarations>>=
1858 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1859 int *num, char ***string_array);
1860 extern void free_parsed_list(int num, char **string_array);
1863 [[parse_list_string]] allocates memory, don't forget to free it
1864 afterward with [[free_parsed_list]]. It does not alter the original.
1866 The string may start off with a [[deeper]] character
1867 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1868 [[x,y]], where the pointer is one character in on the copied string.
1869 However, when we go to free the memory, we need a pointer to the
1870 beginning of the string. In order to accommodate this for a string
1871 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1872 the first $N$ elements point to the separated arguments, and let the
1873 last element point to the start of the copied string regardless of
1875 <<parse delimited list string functions>>=
1876 /* TODO, split out into parse.hc */
1877 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1880 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1881 if (string[i] == deeper) {depth++;}
1882 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1888 void parse_list_string(char *string, char sep, char deeper, char shallower,
1889 int *num, char ***string_array)
1891 char *str=NULL, **ret=NULL;
1893 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
1895 *string_array = NULL;
1898 /* make a copy of the string, so we don't change the original */
1899 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1900 assert(str != NULL);
1901 strcpy(str, string); /* we know str is long enough */
1902 /* count the number of regions, so we can allocate pointers to them */
1905 n++; i++; /* move on to next argument */
1906 i += next_delim_index(str+i, sep, deeper, shallower);
1907 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1908 } while (str[i] != '\0');
1909 ret = (char **)malloc(sizeof(char *)*(n+1));
1910 assert(ret != NULL);
1911 /* replace the separators with '\0' & assign pointers */
1912 ret[n] = str; /* point to the front of the copied string */
1915 for(i=1; i<n; i++) {
1916 j += next_delim_index(str+j, sep, deeper, shallower);
1918 ret[i] = str+j; /* point to the separated arguments */
1920 /* strip off bounding braces, if any */
1921 for(i=0; i<n; i++) {
1922 if (ret[i][0] == deeper) {
1926 if (ret[i][strlen(ret[i])-1] == shallower)
1927 ret[i][strlen(ret[i])-1] = '\0';
1930 *string_array = ret;
1933 void free_parsed_list(int num, char **string_array)
1936 /* frees all items, since they were allocated as a single string*/
1937 free(string_array[num]);
1938 /* and free the array of pointers */
1944 \subsection{Parsing implementation}
1950 <<parse delimited list string functions>>
1954 #include <assert.h> /* assert() */
1955 #include <stdlib.h> /* NULL */
1956 #include <stdio.h> /* fprintf(), stdout *//*!!*/
1957 #include <string.h> /* strlen() */
1961 \subsection{Parsing unit tests}
1963 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1965 <<parse check includes>>
1966 <<string comparison definition and globals>>
1967 <<check relative error>>
1968 <<parse test suite>>
1969 <<main check program>>
1972 <<parse check includes>>=
1973 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
1974 #include <stdio.h> /* printf() */
1975 #include <assert.h> /* assert() */
1976 #include <string.h> /* strlen() */
1981 <<parse test suite>>=
1982 <<parse list string tests>>
1983 <<parse suite definition>>
1986 <<parse suite definition>>=
1987 Suite *test_suite (void)
1989 Suite *s = suite_create ("k model");
1990 <<parse list string test case defs>>
1992 <<parse list string test case add>>
1997 <<parse list string tests>>=
2000 START_TEST(test_next_delim_index)
2002 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
2003 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
2004 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
2005 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
2006 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
2010 START_TEST(test_parse_list_null)
2014 parse_list_string(NULL, SEP, '{', '}',
2015 &num_param_args, ¶m_args);
2016 fail_unless(num_param_args == 0, NULL);
2017 fail_unless(param_args == NULL, NULL);
2020 START_TEST(test_parse_list_single_simple)
2024 parse_list_string("arg", SEP, '{', '}',
2025 &num_param_args, ¶m_args);
2026 fail_unless(num_param_args == 1, NULL);
2027 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
2030 START_TEST(test_parse_list_single_compound)
2034 parse_list_string("{x,y,z}", SEP, '{', '}',
2035 &num_param_args, ¶m_args);
2036 fail_unless(num_param_args == 1, NULL);
2037 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
2040 START_TEST(test_parse_list_double_simple)
2044 parse_list_string("abc,def", SEP, '{', '}',
2045 &num_param_args, ¶m_args);
2046 fail_unless(num_param_args == 2, NULL);
2047 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2048 fail_unless(STRMATCH(param_args[1],"def"), NULL);
2051 START_TEST(test_parse_list_compound)
2055 parse_list_string("abc,{def,ghi}", SEP, '{', '}',
2056 &num_param_args, ¶m_args);
2057 fail_unless(num_param_args == 2, NULL);
2058 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2059 fail_unless(STRMATCH(param_args[1],"def,ghi"), NULL);
2063 <<parse list string test case defs>>=
2064 TCase *tc_parse_list_string = tcase_create("parse list string");
2066 <<parse list string test case add>>=
2067 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
2068 tcase_add_test(tc_parse_list_string, test_parse_list_null);
2069 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
2070 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
2071 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
2072 tcase_add_test(tc_parse_list_string, test_parse_list_compound);
2073 suite_add_tcase(s, tc_parse_list_string);
2077 \section{Unit tests}
2079 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2086 <<check relative error>>
2089 <<main check program>>
2101 <<determine dt tests>>
2103 <<does domain transition tests>>
2104 <<randomly unfold domains tests>>
2105 <<suite definition>>
2108 <<suite definition>>=
2109 Suite *test_suite (void)
2111 Suite *s = suite_create ("sawsim");
2112 <<F test case defs>>
2113 <<determine dt test case defs>>
2114 <<happens test case defs>>
2115 <<does domain transition test case defs>>
2116 <<randomly unfold domains test case defs>>
2119 <<determine dt test case add>>
2120 <<happens test case add>>
2121 <<does domain transition test case add>>
2122 <<randomly unfold domains test case add>>
2125 tcase_add_checked_fixture(tc_strip_address,
2126 setup_strip_address,
2127 teardown_strip_address);
2133 <<main check program>>=
2137 Suite *s = test_suite();
2138 SRunner *sr = srunner_create(s);
2139 srunner_run_all(sr, CK_ENV);
2140 number_failed = srunner_ntests_failed(sr);
2142 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
2146 \subsection{F tests}
2148 <<worm-like chain tests>>
2150 <<F test case defs>>=
2151 <<worm-like chain test case def>>
2153 <<F test case add>>=
2154 <<worm-like chain test case add>>
2157 <<worm-like chain tests>>=
2158 <<worm-like chain function>>
2160 START_TEST(test_wlc_at_zero)
2162 double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
2163 fail_unless(abs(wlc(x, T, p, L)) < lim, \
2164 "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
2165 x, T, p, L, abs(wlc(x, T, p, L)), lim);
2169 START_TEST(test_wlc_at_half)
2171 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
2172 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
2173 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
2175 * linear term = x/L = 0.5
2176 * nonlinear + linear = 0.75 + 0.5 = 1.25
2177 * wlc = 10e21*1.25 = 12.5e21
2179 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
2180 "wlc(%g, %g, %g, %g) = %g != %g",
2181 x, T, p, L, wlc(x, T, p, L), 12.5e21);
2185 <<worm-like chain test case def>>=
2186 TCase *tc_wlc = tcase_create("WLC");
2189 <<worm-like chain test case add>>=
2190 tcase_add_test(tc_wlc, test_wlc_at_zero);
2191 tcase_add_test(tc_wlc, test_wlc_at_half);
2192 suite_add_tcase(s, tc_wlc);
2195 \subsection{Model tests}
2197 Check the searching with [[linear_k]].
2198 Check overwhelming force treatment with the heavyside-esque step [[?]].
2200 Hmm, I'm not really sure what I was doing here...
2201 <<determine dt tests>>=
2202 double linear_k(double F, environment_t *env, void *params)
2204 double Fnot = *(double *)params;
2208 START_TEST(test_determine_dt_linear_k)
2211 transition_t unfold={0};
2212 environment_t env = {0};
2213 double F[]={0,1,2,3,4,5,6};
2214 double dt_max=3.0, Fnot=3.0, x=1.0, targ_p=0.1, v=1.0;
2215 list_t *state_list=NULL, *transition_list=NULL;
2218 folded.tension_handler = &hooke_handler;
2219 folded.num_domains = 1;
2220 unfold.initial_state = &folded;
2221 unfold.k = linear_k;
2222 unfold.k_params = &Fnot;
2223 push(&state_list, &folded, 1);
2224 push(&transition_list, &unfold, 1);
2225 for( i=0; i < sizeof(F)/sizeof(double); i++) {
2226 //fail_unless(determine_dt(state_list, transition_list, &env, x, targ_p, dt_max, v)==1, NULL);
2231 <<determine dt test case defs>>=
2232 TCase *tc_determine_dt = tcase_create("Determine dt");
2234 <<determine dt test case add>>=
2235 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
2236 suite_add_tcase(s, tc_determine_dt);
2241 <<happens test case defs>>=
2243 <<happens test case add>>=
2246 <<does domain transition tests>>=
2248 <<does domain transition test case defs>>=
2250 <<does domain transition test case add>>=
2253 <<randomly unfold domains tests>>=
2255 <<randomly unfold domains test case defs>>=
2257 <<randomly unfold domains test case add>>=
2261 \section{Balancing group extensions}
2262 \label{sec.tension_balance}
2264 For a given total extension $x$ (of the piezo), the various domain
2265 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
2266 amounts, and we need to tweak the portion that each extends to balance
2267 the tension amongst the active groups. Since the tension balancing is
2268 separable from the bulk of the simulation, we break it out into a
2269 separate module. The interface is defined in [[tension_balance.h]],
2270 the implementation in [[tension_balance.c]], and the unit testing in
2271 [[check_tension_balance.c]]
2273 <<tension-balance.h>>=
2274 #ifndef TENSION_BALANCE_H
2275 #define TENSION_BALANCE_H
2277 <<tension balance function declaration>>
2278 #endif /* TENSION_BALANCE_H */
2281 <<tension balance functions>>=
2282 <<one dimensional bracket>>
2283 <<one dimensional solve>>
2285 <<group stiffness function>>
2286 <<chain stiffness function>>
2287 <<tension balance function>>
2290 <<tension balance module makefile lines>>=
2291 NW_SAWSIM_MODS += tension_balance
2292 CHECK_BINS += check_tension_balance
2293 check_tension_balance_MODS =
2296 The entire force balancing problem reduces to a solving two nested
2297 one-dimensional root-finding problems. First we define the one
2298 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
2299 populated group). $x(x_0)$ is also strictly monotonically increasing,
2300 so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
2301 stores the last successful value of $x$, since we expect to be taking
2302 small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
2303 indicates that the groups have changed.
2304 <<tension balance function declaration>>=
2305 double tension_balance(int num_tension_groups,
2306 one_dim_fn_t **tension_handlers,
2307 one_dim_fn_t **inverse_tension_handlers,
2308 void **params, /* array of pointers to tension_handler_data_t */
2313 <<tension balance function>>=
2314 double tension_balance(int num_tension_groups,
2315 one_dim_fn_t **tension_handlers,
2316 one_dim_fn_t **inverse_tension_handlers,
2317 void **params, /* array of pointers to tension_handler_data_t */
2322 static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL,NULL};
2323 double F, xo_guess, xo, lb, ub;
2324 double min_relx=1e-6, min_rely=1e-6;
2325 int max_steps=200, i;
2327 assert(num_tension_groups > 0);
2328 assert(tension_handlers != NULL);
2329 assert(params != NULL);
2330 assert(num_tension_groups > 0);
2332 if (last_x == -1 || x_xo_data.xi == NULL) { /* new group setup, reset x_xo_data */
2333 /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
2334 if (x_xo_data.xi != NULL)
2336 /* construct new x_xo_data */
2337 x_xo_data.n_groups = num_tension_groups;
2338 x_xo_data.tension_handler = tension_handlers;
2339 x_xo_data.inverse_tension_handler = inverse_tension_handlers;
2340 x_xo_data.handler_data = params;
2342 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);
2344 x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
2345 assert(x_xo_data.xi != NULL);
2346 for (i=0; i<num_tension_groups; i++)
2347 x_xo_data.xi[i] = last_x;
2350 if (num_tension_groups == 1) { /* only one group, no balancing required */
2352 x_xo_data.xi[0] = xo;
2354 fprintf(stderr, "%s:\tonly one active group, no balancing required\n", __FUNCTION__);
2358 fprintf(stderr, "%s:\tbalancing for x = %g with", __FUNCTION__, x);
2359 for(i=0; i<num_tension_groups; i++) fprintf(stderr, "\tx[%d]=%g", i, x_xo_data.xi[i]);
2360 fprintf(stderr, "\n");
2362 if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
2363 if (x == 0) xo_guess = length_scale;
2364 else xo_guess = x/num_tension_groups;
2366 fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
2368 oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
2369 } else { /* work off the last balanced point */
2371 fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
2374 lb = x_xo_data.xi[0];
2375 ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */
2376 } else if (x < last_x) {
2377 lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
2378 ub = x_xo_data.xi[0];
2379 } else { /* x == last_x */
2380 fprintf(stderr, "not moving\n");
2381 lb= x_xo_data.xi[0];
2382 ub= x_xo_data.xi[0];
2386 fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
2387 __FUNCTION__, x, lb, ub);
2389 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
2390 F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
2392 if (num_tension_groups > 1) {
2393 fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
2394 for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
2395 fprintf(stderr, "\n");
2400 F = (*tension_handlers[0])(xo, params[0]);
2402 if (stiffness != NULL)
2403 *stiffness = chain_stiffness(&x_xo_data, min_relx);
2410 <<tension balance internal definitions>>=
2411 <<x of x0 definitions>>
2414 <<x of x0 definitions>>=
2415 typedef struct x_of_xo_data_struct {
2417 one_dim_fn_t **tension_handler; /* array of fn pointers */
2418 one_dim_fn_t **inverse_tension_handler; /* inverse fn pointers */
2419 void **handler_data; /* array of pointers to tension_handler_data_t */
2420 double *xi; /* array of group extensions */
2425 double x_of_xo(double xo, void *pdata)
2427 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
2428 double F, x=0, xi, xi_guess, lb, ub, slope;
2430 double min_relx=1e-6, min_rely=1e-6;
2432 assert(data->n_groups > 0);
2435 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);
2438 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
2440 fprintf(stderr, "%s:\tfinding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
2444 if (data->xi) data->xi[0] = xo;
2445 for (i=1; i < data->n_groups; i++) {
2446 if (data->inverse_tension_handler[i] != NULL) {
2447 xi = (*data->inverse_tension_handler[i])(F, data->handler_data[i]);
2448 } else { /* invert numerically */
2449 if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
2450 else xi_guess = data->xi[i];
2452 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]);
2454 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
2455 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
2460 fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
2464 fprintf(stderr, "%s:\tfound x(x0=%g) = %g\n", __FUNCTION__, xo, x);
2470 <<group stiffness function>>=
2471 double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double x, double relx)
2473 double dx, xl, F, Fl, stiffness;
2475 assert(relx > 0 && relx < 1);
2476 if (x == 0 || relx == 0) {
2477 dx = 1e-12; /* HACK, how to get length scale? */
2487 F = (*tension_handler)(x, handler_data);
2488 Fl = (*tension_handler)(xl, handler_data);
2489 stiffness = (F-Fl)/dx;
2490 assert(stiffness > 0);
2495 <<chain stiffness function>>=
2496 double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx)
2500 /* add group stiffnesses in series */
2501 for (i=0; i < x_xo_data->n_groups; i++) {
2503 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);
2506 stiffness += 1.0/group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx);
2508 stiffness = 1.0 / stiffness;
2514 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
2515 number of steps for monotonically increasing $f(x)$. Simple
2516 bisection, so it's robust and converges fairly quickly. You can set
2517 the maximum number of steps to take, but if you have enough processor
2518 power, it's better to limit your precision with the relative limits
2519 [[min_relx]] and [[y]]. These ensure that the width of the bracket is
2520 small on the length scale of it's larger side. Note that \emph{both}
2521 relative limits must be satisfied for the search to stop.
2522 <<one dimensional function definition>>=
2523 /* equivalent to gsl_func_t */
2524 typedef double one_dim_fn_t(double x, void *params);
2526 <<one dimensional solve declaration>>=
2527 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2528 double min_relx, double min_rely, int max_steps,
2532 <<one dimensional solve>>=
2533 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2534 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2535 double min_relx, double min_rely, int max_steps,
2538 double yx, ylb, yub, x;
2541 ylb = (*f)(lb, params);
2542 yub = (*f)(ub, params);
2544 /* check some simple cases */
2546 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bracketed */
2547 else return lb; /* any x on the range [lb, ub] would work */
2549 if (ylb == y) { x=lb; yx=ylb; goto end; }
2550 if (yub == y) { x=ub; yx=yub; goto end; }
2553 fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2555 assert(ylb < y); assert(yub > y);
2556 for (i=0; i<max_steps; i++) {
2558 yx = (*f)(x, params);
2560 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);
2562 if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2563 else if (yx > y) { ub=x; yub=yx; }
2564 else /* yx < y */ { lb=x; ylb=yx; }
2569 fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2575 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2576 Generate bracketing $x$ values through bisection or exponential growth.
2577 <<one dimensional bracket declaration>>=
2578 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2581 <<one dimensional bracket>>=
2582 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2584 double yx, step, x=xguess;
2587 fprintf(stderr, "%s:\tbracketing %g (params %p)\n", __FUNCTION__, y, params);
2589 yx = (*f)(x, params);
2591 fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2598 fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2602 if (x == 0) x = length_scale/2.0; /* HACK */
2605 yx = (*f)(x, params);
2607 fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2612 fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2615 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2619 \subsection{Balancing implementation}
2621 <<tension-balance.c>>=
2624 <<tension balance includes>>
2625 #include "tension_balance.h"
2627 double length_scale = 1e-10; /* HACK */
2629 <<tension balance internal definitions>>
2630 <<tension balance functions>>
2633 <<tension balance includes>>=
2634 #include <assert.h> /* assert() */
2635 #include <stdlib.h> /* NULL */
2636 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2637 #include <stdio.h> /* fprintf(), stdout */
2641 \subsection{Balancing unit tests}
2643 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2644 <<check-tension-balance.c>>=
2645 <<tension balance check includes>>
2646 <<tension balance test suite>>
2647 <<main check program>>
2650 <<tension balance check includes>>=
2651 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2652 #include <assert.h> /* assert() */
2655 #include "tension_balance.h"
2658 <<tension balance test suite>>=
2659 <<tension balance function tests>>
2660 <<tension balance suite definition>>
2663 <<tension balance suite definition>>=
2664 Suite *test_suite (void)
2666 Suite *s = suite_create ("tension balance");
2667 <<tension balance function test case defs>>
2669 <<tension balance function test case adds>>
2674 <<tension balance function tests>>=
2675 <<check relative error>>
2677 double hooke(double x, void *pK)
2680 return *((double*)pK) * x;
2683 START_TEST(test_single_function)
2685 double k=5, x=3, last_x=2.0, F;
2686 one_dim_fn_t *handlers[] = {&hooke};
2687 one_dim_fn_t *inverse_handlers[] = {NULL};
2688 void *data[] = {&k};
2689 F = tension_balance(1, handlers, inverse_handlers, data, last_x, x, NULL);
2690 fail_unless(F = k*x, NULL);
2695 We can also test balancing two springs with different spring constants.
2696 The analytic solution for a general number of springs is given in Appendix \ref{sec.math_hooke}.
2697 <<tension balance function tests>>=
2698 START_TEST(test_double_hooke)
2700 double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
2701 one_dim_fn_t *handlers[] = {&hooke, &hooke};
2702 one_dim_fn_t *inverse_handlers[] = {NULL, NULL};
2703 void *data[] = {&k1, &k2};
2704 F = tension_balance(2, handlers, inverse_handlers, data, last_x, x, NULL);
2708 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2709 //CHECK_ERR(1e-6, x1e, xi[0]);
2710 //CHECK_ERR(1e-6, x2e, xi[1]);
2711 CHECK_ERR(1e-6, Fe, F);
2715 TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
2717 <<tension balance function test case defs>>=
2718 TCase *tc_tbfunc = tcase_create("tension balance function");
2721 <<tension balance function test case adds>>=
2722 tcase_add_test(tc_tbfunc, test_single_function);
2723 tcase_add_test(tc_tbfunc, test_double_hooke);
2724 suite_add_tcase(s, tc_tbfunc);
2730 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2731 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2732 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2738 <<list definitions>>
2739 <<list declarations>>
2743 <<list declarations>>=
2744 <<head and tail declarations>>
2745 <<list length declaration>>
2746 <<push declaration>>
2748 <<list destroy declaration>>
2749 <<list to array declaration>>
2750 <<move declaration>>
2751 <<sort declaration>>
2752 <<list index declaration>>
2753 <<list element declaration>>
2754 <<unique declaration>>
2755 <<list print declaration>>
2759 <<create and destroy node>>
2774 <<list module makefile lines>>=
2775 NW_SAWSIM_MODS += list
2776 CHECK_BINS += check_list
2780 <<list definitions>>=
2781 typedef struct list_struct {
2782 struct list_struct *next;
2783 struct list_struct *prev;
2787 <<comparison function typedef>>
2790 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2791 <<head and tail declarations>>=
2792 list_t *head(list_t *list);
2793 list_t *tail(list_t *list);
2796 list_t *head(list_t *list)
2800 while (list->prev != NULL) {
2806 list_t *tail(list_t *list)
2810 while (list->next != NULL) {
2817 <<list length declaration>>=
2818 int list_length(list_t *list);
2821 int list_length(list_t *list)
2828 while (list->next != NULL) {
2836 [[push]] inserts a node relative to the current position in the list
2837 without changing the current position.
2838 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2839 so we set it to that of the pushed domain.
2840 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2841 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2842 <<push declaration>>=
2843 void push(list_t **pList, void *data, int below);
2846 void push(list_t **pList, void *data, int below)
2848 list_t *list, *new_node;
2849 assert(pList != NULL);
2851 new_node = create_node(data);
2854 else if (below==0) { /* insert above */
2855 if (list->prev != NULL)
2856 list->prev->next = new_node;
2857 new_node->prev = list->prev;
2858 list->prev = new_node;
2859 new_node->next = list;
2860 } else { /* insert below */
2861 if (list->next != NULL)
2862 list->next->prev = new_node;
2863 new_node->next = list->next;
2864 list->next = new_node;
2865 new_node->prev = list;
2870 [[pop]] removes the current domain node, moving the current position
2871 to the node after it, unless that node is [[NULL]], in which case move
2872 the current position to the node before it.
2873 <<pop declaration>>=
2874 void *pop(list_t **pList);
2877 void *pop(list_t **pList)
2879 list_t *list, *popped;
2881 assert(pList != NULL);
2883 assert(list != NULL); /* not an empty list */
2885 /* bypass the popped node */
2886 if (list->prev != NULL)
2887 list->prev->next = list->next;
2888 if (list->next != NULL) {
2889 list->next->prev = list->prev;
2890 *pList = list->next;
2892 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
2894 destroy_node(popped);
2899 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
2900 If the cleanup function is [[NULL]], no cleanup function is called.
2901 The nodes are not popped in any particular order.
2902 <<list destroy declaration>>=
2903 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
2906 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
2910 assert(pList != NULL);
2913 assert(list != NULL); /* not an empty list */
2914 while (list != NULL) {
2916 if (destroy != NULL)
2922 [[list_to_array]] converts a list to an array of pointers, preserving order.
2923 <<list to array declaration>>=
2924 void list_to_array(list_t *list, int *length, void ***array);
2927 void list_to_array(list_t *list, int *pLength, void ***pArray)
2931 assert(list != NULL);
2932 assert(pLength != NULL);
2933 assert(pArray != NULL);
2935 length = list_length(list);
2936 array = (void **)malloc(sizeof(void *)*length);
2937 assert(array != NULL);
2940 while (list != NULL) {
2941 array[i++] = list->d;
2949 [[move]] moves the current element along the list in either direction.
2950 <<move declaration>>=
2951 void move(list_t **pList, int downstream);
2954 void move(list_t **pList, int downstream)
2956 list_t *list, *flipper;
2958 assert(pList != NULL);
2959 list = *pList; /* a>B>c>d */
2960 assert(list != NULL); /* not an empty list */
2961 if (downstream == 0)
2962 flipper = list->prev; /* flipper is A */
2964 flipper = list->next; /* flipper is C */
2965 /* ensure there is somewhere to go in the direction we're heading */
2966 assert(flipper != NULL);
2967 /* remove the current element */
2968 data = pop(&list); /* data=B, a>C>d */
2969 /* and put it back in in the correct spot */
2971 if (downstream == 0) {
2972 push(&list, data, 0); /* b>A>c>d */
2973 list = list->prev; /* B>a>c>d */
2975 push(&list, data, 1); /* a>C>b>d */
2976 list = list->next; /* a>c>B>d */
2982 [[sort]] sorts a list from smallest to largest according to a user
2983 specified comparison.
2984 <<comparison function typedef>>=
2985 typedef int (comparison_fn_t)(void *A, void *B);
2988 <<int comparison function>>=
2989 int int_comparison(void *A, void *B)
2994 if (a > b) return 1;
2995 else if (a == b) return 0;
3000 <<sort declaration>>=
3001 void sort(list_t **pList, comparison_fn_t *comp);
3003 Here we use the bubble sort algorithm.
3005 void sort(list_t **pList, comparison_fn_t *comp)
3009 assert(pList != NULL);
3011 assert(list != NULL); /* not an empty list */
3012 while (stable == 0) {
3015 while (list->next != NULL) {
3016 if ((*comp)(list->d, list->next->d) > 0) {
3018 move(&list, 1 /* downstream */);
3028 [[list_index]] finds the location of [[data]] in [[list]]. Returns
3029 [[-1]] if [[data]] is not in [[list]].
3030 <<list index declaration>>=
3031 int list_index(list_t *list, void *data);
3035 int list_index(list_t *list, void *data)
3039 while (list != NULL) {
3040 if (list->d == data) return i;
3049 [[list_element]] returns the data bound to the $i^\text{th}$ node of the list.
3050 <<list element declaration>>=
3051 void *list_element(list_t *list, int i);
3054 void *list_element(list_t *list, int i)
3058 while (list != NULL) {
3059 if (i == j) return list->d;
3069 [[unique]] removes duplicates from a sorted list according to a user
3070 specified comparison. Don't do this unless you have other pointers
3071 any dynamically allocated data in your list, because [[unique]] just
3072 drops any non-unique data without freeing it.
3073 <<unique declaration>>=
3074 void unique(list_t **pList, comparison_fn_t *comp);
3077 void unique(list_t **pList, comparison_fn_t *comp)
3080 assert(pList != NULL);
3082 assert(list != NULL); /* not an empty list */
3084 while (list->next != NULL) {
3085 if ((*comp)(list->d, list->next->d) == 0) {
3094 [[list_print]] prints a list to a [[FILE *]] stream.
3095 <<list print declaration>>=
3096 void list_print(FILE *file, list_t *list, char *list_name);
3099 void list_print(FILE *file, list_t *list, char *list_name)
3101 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
3103 while (list != NULL) {
3104 fprintf(file, " %p", list->d);
3107 fprintf(file, "\n");
3111 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
3112 <<create and destroy node>>=
3113 list_t *create_node(void *data)
3115 list_t *ret = (list_t *)malloc(sizeof(list_t));
3116 assert(ret != NULL);
3123 void destroy_node(list_t *node)
3129 The user must free the data pointed to by the node on their own.
3131 \subsection{List implementation}
3141 #include <assert.h> /* assert() */
3142 #include <stdlib.h> /* malloc(), free(), rand() */
3143 #include <stdio.h> /* fprintf(), stdout, FILE */
3144 #include "global.h" /* destroy_data_func_t */
3147 \subsection{List unit tests}
3149 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3151 <<list check includes>>
3153 <<main check program>>
3156 <<list check includes>>=
3157 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3158 #include <stdio.h> /* FILE */
3164 <<list test suite>>=
3167 <<list suite definition>>
3170 <<list suite definition>>=
3171 Suite *test_suite (void)
3173 Suite *s = suite_create ("list");
3174 <<push test case defs>>
3175 <<pop test case defs>>
3177 <<push test case adds>>
3178 <<pop test case adds>>
3184 START_TEST(test_push)
3189 push(&list, (void *)i, 1);
3190 fail_unless(list == head(list), NULL);
3191 fail_unless( (int)list->d == 0 );
3192 for (i=0; i<n; i++) {
3193 /* we expect 0, n-1, n-2, ..., 1 */
3196 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
3202 <<push test case defs>>=
3203 TCase *tc_push = tcase_create("push");
3206 <<push test case adds>>=
3207 tcase_add_test(tc_push, test_push);
3208 suite_add_tcase(s, tc_push);
3213 <<pop test case defs>>=
3215 <<pop test case adds>>=
3218 \section{Function string evaluation}
3220 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).
3221 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
3222 The solution is to run a scripting language as a subprocess accessed via pipes.
3223 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
3225 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
3226 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
3227 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.
3228 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
3230 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.
3231 Then you can either statically or dynamically link to those hardcoded functions.
3232 While much less flexible, this approach would be a fairly simple-to-implement fallback.
3234 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
3235 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
3238 #ifndef STRING_EVAL_H
3239 #define STRING_EVAL_H
3241 <<string eval setup declaration>>
3242 <<string eval function declaration>>
3243 <<string eval teardown declaration>>
3244 #endif /* STRING_EVAL_H */
3247 <<string eval module makefile lines>>=
3248 NW_SAWSIM_MODS += string_eval
3249 CHECK_BINS += check_string_eval
3250 check_string_eval_MODS =
3253 For an introduction to POSIX process control, see\\
3254 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
3255 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
3256 and of course, the relavent [[man]] pages.
3258 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
3259 [[execvp]] replaces the calling process' program with a new program.
3260 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
3261 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
3262 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
3263 The new program needs command line arguments, just like it would if you were running it from a shell.
3264 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
3265 with the final array entry being a [[NULL]] pointer.
3267 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
3268 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
3269 <<string eval subprocess definitions>>=
3270 #define SUBPROCESS "python"
3271 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
3272 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
3273 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
3275 The [[i]] option lets Python know that it should run in interactive mode.
3276 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
3277 In interactive mode, python acts on every instruction as soon as it is recieved.
3278 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.
3279 %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.
3281 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
3282 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
3283 [[fork]] returns two copies of the same program, executing the original code.
3284 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
3285 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
3287 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.
3288 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
3289 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
3290 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
3291 <<string eval pipe definitions>>=
3292 #define PIPE_READ 0 /* the end you read from */
3293 #define PIPE_WRITE 1 /* the end you write to */
3295 #define STDIN 0 /* initial index of stdin pair */
3296 #define STDOUT 2 /* initial index of stdout pair */
3299 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
3301 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.
3303 <<string eval setup declaration>>=
3304 extern void string_eval_setup(FILE **pIn, FILE **pOut);
3306 <<string eval setup definition>>=
3307 void string_eval_setup(FILE **pIn, FILE **pOut)
3310 int pfd[4]; /* file descriptors for stdin and stdout */
3312 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
3313 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
3315 pid = fork(); /* split process into two copies */
3317 if (pid == -1) { /* fork error */
3318 perror("fork error");
3320 } else if (pid == 0) { /* child */
3321 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
3322 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
3323 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
3324 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
3325 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
3326 perror("exec error"); /* exec shouldn't return */
3328 } else { /* parent */
3329 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
3330 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
3331 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
3332 if ( *pIn == NULL ) {
3333 perror("fdopen (in)");
3336 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
3337 if ( *pOut == NULL ) {
3338 perror("fdopen (out)");
3345 To use the evaluating subprocess, we just pipe in our command, and read out the results.
3346 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
3347 <<string eval function declaration>>=
3348 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
3350 <<string eval function definition>>=
3351 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
3354 rc = fprintf(in, "%s", input);
3355 assert(rc == strlen(input));
3358 alarm(1); /* set a one second timeout on the read */
3359 assert( fgets(output, buflen, out) != NULL );
3360 alarm(0); /* cancel the timeout */
3361 //fprintf(stderr, "eval: %s ----> %s", input, output);
3364 The [[alarm]] calls set and clear a timeout on the returned output.
3365 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.
3366 This protects against invalid input for which a line of output is not printed to [[stdout]].
3367 Other invalid inputs (e.g.\ those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
3368 If you are getting strange results, check your python code seperately. TODO, better error handling.
3370 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
3371 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
3372 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
3373 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
3374 <<string eval teardown declaration>>=
3375 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
3378 <<string eval teardown definition>>=
3379 void string_eval_teardown(FILE **pIn, FILE **pOut)
3384 /* redirect Python's stderr */
3385 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
3389 assert( fclose(*pIn) == 0 );
3391 assert( fclose(*pOut) == 0 );
3394 /* wait for python to exit */
3396 pid = wait(&stat_loc);
3403 if (WIFEXITED(stat_loc)) {
3404 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
3405 } else if (WIFSIGNALED(stat_loc)) {
3406 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
3411 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
3413 \subsection{String evaluation implementation}
3417 <<string eval includes>>
3418 #include "string_eval.h"
3419 <<string eval internal definitions>>
3420 <<string eval functions>>
3423 <<string eval includes>>=
3424 #include <assert.h> /* assert() */
3425 #include <stdlib.h> /* NULL */
3426 #include <stdio.h> /* fprintf(), stdout, fdopen() */
3427 #include <string.h> /* strlen() */
3428 #include <sys/types.h> /* pid_t */
3429 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
3430 #include <wait.h> /* wait() */
3433 <<string eval internal definitions>>=
3434 <<string eval subprocess definitions>>
3435 <<string eval pipe definitions>>
3438 <<string eval functions>>=
3439 <<string eval setup definition>>
3440 <<string eval function definition>>
3441 <<string eval teardown definition>>
3444 \subsection{String evaluation unit tests}
3446 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3447 <<check-string-eval.c>>=
3448 <<string eval check includes>>
3449 <<string comparison definition and globals>>
3450 <<string eval test suite>>
3451 <<main check program>>
3454 <<string eval check includes>>=
3455 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3456 #include <stdio.h> /* printf() */
3457 #include <assert.h> /* assert() */
3458 #include <string.h> /* strlen() */
3459 #include <signal.h> /* SIGKILL */
3461 #include "string_eval.h"
3464 <<string eval test suite>>=
3465 <<string eval tests>>
3466 <<string eval suite definition>>
3469 <<string eval suite definition>>=
3470 Suite *test_suite (void)
3472 Suite *s = suite_create ("string eval");
3473 <<string eval test case defs>>
3475 <<string eval test case add>>
3480 <<string eval tests>>=
3481 START_TEST(test_setup_teardown)
3484 string_eval_setup(&in, &out);
3485 string_eval_teardown(&in, &out);
3488 START_TEST(test_invalid_command)
3491 char input[80], output[80]={};
3492 string_eval_setup(&in, &out);
3493 sprintf(input, "print ABCDefg\n");
3494 string_eval(in, out, input, 80, output);
3495 string_eval_teardown(&in, &out);
3498 START_TEST(test_simple_eval)
3501 char input[80], output[80]={};
3502 string_eval_setup(&in, &out);
3503 sprintf(input, "print 3+4*5\n");
3504 string_eval(in, out, input, 80, output);
3505 fail_unless(STRMATCH(output,"23\n"), NULL);
3506 string_eval_teardown(&in, &out);
3509 START_TEST(test_multiple_evals)
3512 char input[80], output[80]={};
3513 string_eval_setup(&in, &out);
3514 sprintf(input, "print 3+4*5\n");
3515 string_eval(in, out, input, 80, output);
3516 fail_unless(STRMATCH(output,"23\n"), NULL);
3517 sprintf(input, "print (3**2 + 4**2)**0.5\n");
3518 string_eval(in, out, input, 80, output);
3519 fail_unless(STRMATCH(output,"5.0\n"), NULL);
3520 string_eval_teardown(&in, &out);
3523 START_TEST(test_eval_with_state)
3526 char input[80], output[80]={};
3527 string_eval_setup(&in, &out);
3528 sprintf(input, "print 3+4*5\n");
3529 fprintf(in, "A = 3\n");
3530 sprintf(input, "print A*3\n");
3531 string_eval(in, out, input, 80, output);
3532 fail_unless(STRMATCH(output,"9\n"), NULL);
3533 string_eval_teardown(&in, &out);
3537 <<string eval test case defs>>=
3538 TCase *tc_string_eval = tcase_create("string_eval");
3540 <<string eval test case add>>=
3541 tcase_add_test(tc_string_eval, test_setup_teardown);
3542 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
3543 tcase_add_test(tc_string_eval, test_simple_eval);
3544 tcase_add_test(tc_string_eval, test_multiple_evals);
3545 tcase_add_test(tc_string_eval, test_eval_with_state);
3546 suite_add_tcase(s, tc_string_eval);
3550 \section{Accelerating function evaluation}
3552 My first version-0.3 code was running very slowly.
3553 With the options suggested in the help
3554 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
3555 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
3556 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
3557 That's only 15 calls per solution, so the search algorithm seems reasonable.
3558 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
3563 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
3565 #endif /* ACCEL_K_H */
3568 <<accel k module makefile lines>>=
3569 NW_SAWSIM_MODS += accel_k
3570 #CHECK_BINS += check_accel_k
3571 check_accel_k_MODS =
3575 #include <assert.h> /* assert() */
3576 #include <stdlib.h> /* realloc(), free(), NULL */
3577 #include "global.h" /* environment_t */
3578 #include "k_model.h" /* k_func_t */
3579 #include "interp.h" /* interp_* */
3580 #include "accel_k.h"
3582 typedef struct accel_k_struct {
3583 interp_table_t *itable;
3589 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3590 static int num_accels = 0;
3591 static accel_k_t *accels=NULL;
3593 /* Wrap k in the standard f(x) acceleration form */
3594 static double k_wrap(double F, void *params)
3596 accel_k_t *p = (accel_k_t *)params;
3597 return (*p->k)(F, p->env, p->params);
3600 static int k_tol(double FA, double kA, double FB, double kB)
3603 if (FB-FA > 1e-12) {
3604 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3605 return 1; /* unnacceptable */
3607 //printf("acceptable tol\n");
3608 return 0; /* acceptable */
3612 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3615 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3616 assert(accels != NULL);
3617 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3619 accels[i].env = env;
3620 accels[i].params = params;
3627 for (i=0; i<num_accels; i++)
3628 interp_table_free(accels[i].itable);
3630 if (accels) free(accels);
3634 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3637 for (i=0; i<num_accels; i++) {
3638 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3639 /* We've already setup this function.
3640 * WARNING: we're assuming param and env are the same. */
3641 return interp_table_eval(accels[i].itable, accels+i, F);
3644 /* set up a new acceleration environment */
3645 i = add_accel_k(k, env, params);
3646 return interp_table_eval(accels[i].itable, accels+i, F);
3650 \section{Tension models}
3651 \label{sec.tension_models}
3653 TODO: tension model intro.
3654 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.
3656 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3657 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]].
3659 <<tension-model.h>>=
3660 #ifndef TENSION_MODEL_H
3661 #define TENSION_MODEL_H
3663 <<tension handler types>>
3664 <<constant tension model declarations>>
3665 <<hooke tension model declarations>>
3666 <<worm-like chain tension model declarations>>
3667 <<freely-jointed chain tension model declarations>>
3668 <<find tension definitions>>
3669 <<tension model global declarations>>
3670 #endif /* TENSION_MODEL_H */
3673 <<tension model module makefile lines>>=
3674 NW_SAWSIM_MODS += tension_model
3675 #CHECK_BINS += check_tension_model
3676 check_tension_model_MODS =
3678 <<tension model utils makefile lines>>=
3679 TENSION_MODEL_MODS = tension_model parse list tension_balance
3680 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3681 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3682 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3683 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
3684 notangle -Rtension-model-utils.c $< > $@
3685 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
3686 gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3687 clean_tension_model_utils :
3688 rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3689 @ The pipe symbol [[|]] marks the prerequisits following it (in this
3690 case, the directories) as ‘order-only’ prerequisites. The timestamp
3691 on these prerequisits does not effect whether the rules are executed.
3692 This is appropriate for directories, where we don't need to recompile
3693 after an unrelated has been added to the directory, but only when the
3694 source prerequisites change. See the [[make]] documentation for more
3696 (\url{http://www.gnu.org/software/make/manual/html_node/Prerequisite-Types.html}).
3699 \label{sec.null_tension}
3701 For unstretchable domains.
3703 <<null tension model getopt>>=
3704 {NULL, NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3707 \subsection{Constant}
3708 \label{sec.const_tension}
3710 An infinitely stretchable domain providing a constant tension.
3711 Obviously this cannot be inverted, so you can't put this domain in
3712 series with any other domains. We include it mostly for testing
3715 <<constant tension functions>>=
3716 <<constant tension function>>
3717 <<constant tension parameter create/destroy functions>>
3720 <<constant tension model declarations>>=
3721 <<constant tension function declaration>>
3722 <<constant tension parameter create/destroy function declarations>>
3723 <<constant tension model global declarations>>
3727 <<constant tension function declaration>>=
3728 extern double const_tension_handler(double x, void *pdata);
3730 <<constant tension function>>=
3731 double const_tension_handler(double x, void *pdata)
3733 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3734 list_t *list = data->group_tension_params;
3739 assert(list != NULL); /* empty active group?! */
3740 F = ((const_tension_param_t *)list->d)->F;
3742 fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3744 while (list != NULL) {
3745 assert(((const_tension_param_t *)list->d)->F == F);
3753 <<constant tension structure>>=
3754 typedef struct const_tension_param_struct {
3755 double F; /* tension (force) in N */
3756 } const_tension_param_t;
3760 <<constant tension parameter create/destroy function declarations>>=
3761 extern void *string_create_const_tension_param_t(char **param_strings);
3762 extern void destroy_const_tension_param_t(void *p);
3765 <<constant tension parameter create/destroy functions>>=
3766 const_tension_param_t *create_const_tension_param_t(double F)
3768 const_tension_param_t *ret
3769 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3770 assert(ret != NULL);
3775 void *string_create_const_tension_param_t(char **param_strings)
3777 return create_const_tension_param_t(
3778 safe_strtod(param_strings[0],__FUNCTION__));
3781 void destroy_const_tension_param_t(void *p)
3789 <<constant tension model global declarations>>=
3790 extern int num_const_tension_params;
3791 extern char *const_tension_param_descriptions[];
3792 extern char const_tension_param_string[];
3795 <<constant tension model globals>>=
3796 int num_const_tension_params = 1;
3797 char *const_tension_param_descriptions[] = {"tension F, N"};
3798 char const_tension_param_string[] = "0";
3801 <<constant tension model getopt>>=
3802 {&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}
3808 <<hooke functions>>=
3809 <<internal hooke functions>>
3811 <<inverse hooke handler>>
3812 <<hooke parameter create/destroy functions>>
3815 <<hooke tension model declarations>>=
3816 <<hooke tension function declaration>>
3817 <<hooke tension parameter create/destroy function declarations>>
3818 <<hooke tension model global declarations>>
3822 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
3823 The behavior of a series of springs $k_i$ in series is given by
3825 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
3827 For a simple proof, see Appendix \ref{sec.math_hooke}.
3829 <<hooke tension function declaration>>=
3830 extern double hooke_handler(double x, void *pdata);
3831 extern double inverse_hooke_handler(double F, void *pdata);
3834 First we define a function that computes the effective spring constant
3835 (stored in a single [[hooke_param_t]]) for the entire group.
3836 <<internal hooke functions>>=
3837 static void hooke_param_grouper(tension_handler_data_t *data,
3838 hooke_param_t *param) {
3839 list_t *list = data->group_tension_params;
3843 assert(list != NULL); /* empty active group?! */
3844 while (list != NULL) {
3845 assert( ((hooke_param_t *)list->d)->k > 0 );
3846 k += 1.0/ ((hooke_param_t *)list->d)->k;
3848 fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
3849 __FUNCTION__, ((hooke_param_t *)list->d)->k);
3855 fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N (from data %p)\n",
3856 __FUNCTION__, k, x, k*x, data);
3863 Everyone knows Hooke's law: $F=kx$.
3865 double hooke_handler(double x, void *pdata)
3867 hooke_param_t param = {0};
3870 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
3876 Inverting Hooke's law gives $x=F/k$.
3877 <<inverse hooke handler>>=
3878 double inverse_hooke_handler(double F, void *pdata)
3880 hooke_param_t param = {0};
3883 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
3889 Not much to keep track of for a single effective spring, just the
3890 spring constant $k$.
3891 <<hooke structure>>=
3892 typedef struct hooke_param_struct {
3893 double k; /* spring constant in N/m */
3898 We wrap up our Hooke implementation with some book-keeping functions.
3899 <<hooke tension parameter create/destroy function declarations>>=
3900 extern void *string_create_hooke_param_t(char **param_strings);
3901 extern void destroy_hooke_param_t(void *p);
3905 <<hooke parameter create/destroy functions>>=
3906 hooke_param_t *create_hooke_param_t(double k)
3908 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
3909 assert(ret != NULL);
3914 void *string_create_hooke_param_t(char **param_strings)
3916 return create_hooke_param_t(safe_strtod(param_strings[0], __FUNCTION__));
3919 void destroy_hooke_param_t(void *p)
3926 <<hooke tension model global declarations>>=
3927 extern int num_hooke_params;
3928 extern char *hooke_param_descriptions[];
3929 extern char hooke_param_string[];
3932 <<hooke tension model globals>>=
3933 int num_hooke_params = 1;
3934 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
3935 char hooke_param_string[]="0.05";
3938 <<hooke tension model getopt>>=
3939 {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}
3942 \subsection{Worm-like chain}
3945 We can model several unfolded domains with a single worm-like chain.
3946 <<worm-like chain functions>>=
3947 <<internal worm-like chain functions>>
3948 <<worm-like chain handler>>
3949 <<inverse worm-like chain handler>>
3950 <<worm-like chain create/destroy functions>>
3953 <<worm-like chain tension model declarations>>=
3955 <<worm-like chain handler declaration>>
3956 <<inverse worm-like chain handler declaration>>
3957 <<worm-like chain create/destroy function declarations>>
3958 <<worm-like chain tension model global declarations>>
3962 <<internal worm-like chain functions>>=
3963 <<worm-like chain function>>
3964 <<inverse worm-like chain function>>
3965 <<worm-like chain parameter grouper>>
3968 The combination of all unfolded domains is modeled as a worm like chain,
3969 with the force $F_{\text{WLC}}$ approximately given by
3971 F_{\text{WLC}} \approx \frac{k_B T}{p}
3972 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
3974 where $x$ is the distance between the fixed ends,
3975 $k_B$ is Boltzmann's constant,
3976 $T$ is the absolute temperature,
3977 $p$ is the persistence length, and
3978 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
3979 <<worm-like chain function>>=
3980 static double wlc(double x, double T, double p, double L)
3982 double xL=0.0; /* uses global k_B */
3984 assert(T > 0); assert(p > 0); assert(L > 0);
3985 if (x >= L) return HUGE_VAL;
3986 xL = x/L; /* unitless */
3987 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
3988 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
3989 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
3994 Inverting the approximate WLC is fairly straightforward, if a bit tedious.
3995 Letting $x_L \equiv x/L$ and $F_T \equiv Fp/k_B T$ we have
3997 F_T &\approx \frac{1}{4} \p({\frac{1}{(1-x_L)^2} - 1}) + x_L \\
3998 F_T (1-x_L)^2 &\approx \frac{1}{4} - \frac{1}{4} (1-x_L)^2 + x_L(1-x_L)^2 \\
3999 0 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) (1-x_L)^2 + x_L(1-x_L)^2 \\
4000 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) \p({1-2x_L+x_L^2})
4001 + x_L - 2x_L^2 + x_L^3 \\
4002 &\approx \p[{\frac{1}{4} - \p({F_T + \frac{1}{4}})}]
4003 + x_L \p[{2F_T + \frac{1}{2} + 1}]
4004 + x_L^2 \p[{-\p({F_T + \frac{1}{4}}) - 2}]
4007 + x_L \p({2F_T + \frac{3}{2}})
4008 - x_L^2 \p({F_T + \frac{9}{4}})
4010 &\approx ax_L^3 + bx_L^2 + cx_L + d
4012 where $a \equiv 1$, $b \equiv -F_T - 9/4$, $c \equiv 2F_T + 3/2$, and
4014 % From \citet{wikipedia_cubic_function} the discriminant
4016 % \Delta &= -4b^3d + b^2c^2 - 4ac^3 + 18abcd - 27a^2d^2 \\
4017 % &= -4F_T\p({F_T + \frac{9}{4}})^3
4018 % + \p({F_T + \frac{9}{4}})^2 \cdot \p({2F_T + \frac{3}{2}})^2
4019 % - 4 \p({2F_T + \frac{3}{2}})^3
4020 % + 18F_T \p({F_T + \frac{9}{4}}) \cdot \p({2F_T + \frac{3}{2}})
4022 % &= -4F_T \p({F_T^3 + \frac{27}{4}F_T^2 +\frac{243}{16}F_T +\frac{729}{64}})
4023 %% a^3 + 3a^2b + 3ab^2 + b^3
4024 % + \p({F_T^2 + \frac{9}{2}F_T + \frac{81}{16}})
4025 % \cdot \p({4F_T^2 + 6F_T + \frac{9}{4}}) \\
4026 % &\qquad - 4 \p({8F_T^3 + 18F_T^2 + \frac{27}{2}F_T + \frac{27}{8}})
4027 %% a^3 + 3a^2b + 3ab^2 + b^3
4028 % + 18F_T \p({2F_T^2 + 6F_T + \frac{27}{8}})
4030 % &= -4F_T^4 - 27F_T^3 - \frac{243}{4}F_T^2 - \frac{729}{16}F_T
4031 % + 4F_T^4 + 24F_T^3 +\frac{99}{2}F_T^2 +\frac{81}{2}F_T +\frac{729}{64}\\
4032 % &\qquad - 32F_T^3 - 72F_T^2 - 54F_T - \frac{27}{2}
4033 % + 36F_T^3 + 108F_T^2 + \frac{243}{4}F_T
4035 % &= F_T^4 \p({-4 + 4}) + F_T^3 \p({-27 + 24 - 32 + 36})
4036 % + F_T^2 \p({-\frac{243}{4} + \frac{99}{2} - 72 + 108 -27}) \\
4037 % &\qquad + F_T \p({-\frac{729}{16} + \frac{81}{2} - 54 + \frac{243}{4}})
4038 % \p({\frac{729}{64} - \frac{27}{2}}) \\
4039 % &= F_T^3 - \frac{9}{4}F_T^2 + \frac{27}{16}F_T - \frac{135}{64}
4041 We can plug these coefficients into the GSL cubic solver to invert the
4042 WLC\citep{galassi05}. The applicable root is always the one which
4043 comes before the singularity, which will be the smallest real root.
4044 <<inverse worm-like chain function>>=
4045 static double inverse_wlc(double F, double T, double p, double L)
4047 double FT = F*p/(k_B*T); /* uses global k_B */
4048 double xL0, xL1, xL2;
4051 assert(T > 0); assert(p > 0); assert(L > 0);
4054 num_roots = gsl_poly_solve_cubic(-(FT+2.25),2*FT+1.5,-FT, &xL0, &xL1, &xL2);
4055 assert(num_roots > 0); assert(xL0 >= -DOUBLE_PRECISION); assert(xL0 < 1);
4056 if (xL0 < 0) xL0 = 0.0;
4062 First we define a function that computes the effective WLC parameters
4063 (stored in a single [[wlc_param_t]]) for the entire group.
4064 <<worm-like chain parameter grouper>>=
4065 static void wlc_param_grouper(tension_handler_data_t *data,
4066 wlc_param_t *param) {
4067 list_t *list = data->group_tension_params;
4068 double p=0.0, L=0.0;
4071 assert(list != NULL); /* empty active group?! */
4072 p = ((wlc_param_t *)list->d)->p;
4073 while (list != NULL) {
4074 /* enforce consistent p */
4075 assert( ((wlc_param_t *)list->d)->p == p);
4076 L += ((wlc_param_t *)list->d)->L;
4078 fprintf(stderr, "%s: WLC section %g m long in series\n",
4079 __FUNCTION__, ((wlc_param_t *)list->d)->L);
4084 fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
4085 __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
4093 <<worm-like chain handler declaration>>=
4094 extern double wlc_handler(double x, void *pdata);
4097 This model requires that each unfolded domain in the group have the
4098 same persistence length.
4099 <<worm-like chain handler>>=
4100 double wlc_handler(double x, void *pdata)
4102 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4103 wlc_param_t param = {0};
4105 wlc_param_grouper(data, ¶m);
4106 return wlc(x, data->env->T, param.p, param.L);
4111 <<inverse worm-like chain handler declaration>>=
4112 extern double inverse_wlc_handler(double F, void *pdata);
4115 <<inverse worm-like chain handler>>=
4116 double inverse_wlc_handler(double F, void *pdata)
4118 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4119 wlc_param_t param = {0};
4121 wlc_param_grouper(data, ¶m);
4122 return inverse_wlc(F, data->env->T, param.p, param.L);
4127 <<worm-like chain structure>>=
4128 typedef struct wlc_param_struct {
4129 double p; /* WLC persistence length (m) */
4130 double L; /* WLC contour length (m) */
4134 <<worm-like chain create/destroy function declarations>>=
4135 extern void *string_create_wlc_param_t(char **param_strings);
4136 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
4139 <<worm-like chain create/destroy functions>>=
4140 wlc_param_t *create_wlc_param_t(double p, double L)
4142 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
4143 assert(ret != NULL);
4149 void *string_create_wlc_param_t(char **param_strings)
4151 return create_wlc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4152 safe_strtod(param_strings[1], __FUNCTION__));
4155 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
4163 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.
4164 TODO (now we copy the string before parsing, but still do this for clarity).
4166 <<valid string write code>>=
4167 char string[] = "abc";
4170 <<valid string write code>>=
4171 char *string = "abc";
4175 <<worm-like chain tension model global declarations>>=
4176 extern int num_wlc_params;
4177 extern char *wlc_param_descriptions[];
4178 extern char wlc_param_string[];
4180 <<worm-like chain tension model globals>>=
4181 int num_wlc_params = 2;
4182 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
4183 char wlc_param_string[]="0.39e-9,28.4e-9";
4186 <<worm-like chain tension model getopt>>=
4187 {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}
4189 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
4191 \subsection{Freely-jointed chain}
4194 <<freely-jointed chain functions>>=
4195 <<freely-jointed chain function>>
4196 <<freely-jointed chain handler>>
4197 <<freely-jointed chain create/destroy functions>>
4200 <<freely-jointed chain tension model declarations>>=
4201 <<freely-jointed chain handler declaration>>
4202 <<freely-jointed chain create/destroy function declarations>>
4203 <<freely-jointed chain tension model global declarations>>
4207 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
4208 The tension of a chain of $N$ such links is given by
4210 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
4212 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}.
4213 <<freely-jointed chain function>>=
4214 double langevin(double x, void *params)
4217 return 1.0/tanh(x) - 1/x;
4219 <<one dimensional bracket declaration>>
4220 <<one dimensional solve declaration>>
4221 double inv_langevin(double x)
4223 double lb=0.0, ub=1.0, ret;
4224 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
4225 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
4229 double fjc(double x, double T, double l, int N)
4231 double L = l*(double)N;
4232 if (x == 0) return 0;
4233 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
4234 return k_B*T/l * inv_langevin(x/L);
4238 <<freely-jointed chain handler declaration>>=
4239 extern double fjc_handler(double x, void *pdata);
4241 <<freely-jointed chain handler>>=
4242 double fjc_handler(double x, void *pdata)
4244 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4245 list_t *list = data->group_tension_params;
4250 assert(list != NULL); /* empty active group?! */
4251 l = ((fjc_param_t *)list->d)->l;
4252 while (list != NULL) {
4253 /* enforce consistent l */
4254 assert( ((fjc_param_t *)list->d)->l == l);
4255 N += ((fjc_param_t *)list->d)->N;
4257 fprintf(stderr, "%s: FJC section %d links long in series\n",
4258 __FUNCTION__, ((fjc_param_t *)list->d)->N);
4263 fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
4264 __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
4266 return fjc(x, data->env->T, l, N);
4270 <<freely-jointed chain structure>>=
4271 typedef struct fjc_param_struct {
4272 double l; /* FJC link length (m) */
4273 int N; /* FJC number of links */
4277 <<freely-jointed chain create/destroy function declarations>>=
4278 extern void *string_create_fjc_param_t(char **param_strings);
4279 extern void destroy_fjc_param_t(void *p);
4282 <<freely-jointed chain create/destroy functions>>=
4283 fjc_param_t *create_fjc_param_t(double l, double N)
4285 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
4286 assert(ret != NULL);
4292 void *string_create_fjc_param_t(char **param_strings)
4294 return create_fjc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4295 safe_strtod(param_strings[1], __FUNCTION__));
4298 void destroy_fjc_param_t(void *p)
4305 <<freely-jointed chain tension model global declarations>>=
4306 extern int num_fjc_params;
4307 extern char *fjc_param_descriptions[];
4308 extern char fjc_param_string[];
4311 <<freely-jointed chain tension model globals>>=
4312 int num_fjc_params = 2;
4313 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
4314 char fjc_param_string[]="0.5e-9,200";
4317 <<freely-jointed chain tension model getopt>>=
4318 {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}
4321 \subsection{Tension model implementation}
4323 <<tension-model.c>>=
4326 <<tension model includes>>
4327 #include "tension_model.h"
4328 <<tension model internal definitions>>
4329 <<tension model globals>>
4330 <<tension model functions>>
4333 <<tension model includes>>=
4334 #include <assert.h> /* assert() */
4335 #include <stdlib.h> /* NULL, strto*() */
4336 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4337 #include <stdio.h> /* fprintf(), stdout */
4338 #include <errno.h> /* errno, ERANGE */
4341 #include "tension_balance.h" /* oneD_*() */
4344 <<tension model internal definitions>>=
4345 <<constant tension structure>>
4347 <<worm-like chain structure>>
4348 <<freely-jointed chain structure>>
4351 <<tension model globals>>=
4352 <<tension handler array global>>
4353 <<constant tension model globals>>
4354 <<hooke tension model globals>>
4355 <<worm-like chain tension model globals>>
4356 <<freely-jointed chain tension model globals>>
4359 <<tension model functions>>=
4361 <<constant tension functions>>
4363 <<worm-like chain functions>>
4364 <<freely-jointed chain functions>>
4368 \subsection{Utilities}
4370 The tension models can get complicated, and users may want to reassure
4371 themselves that this portion of the input physics are functioning properly. The
4372 stand-alone program [[tension_model_utils]] demonstrates the tension model
4373 interface without the overhead of having to understand a full simulation.
4374 [[tension_model_utils]] takes command line model arguments like the full
4375 simulation, and prints $F(x)$ data to the screen for the selected model over a
4378 <<tension-model-utils.c>>=
4380 <<tension model utility includes>>
4381 <<tension model utility definitions>>
4382 <<tension model utility getopt functions>>
4384 <<tension model utility main>>
4387 <<tension model utility main>>=
4388 int main(int argc, char **argv)
4390 <<tension model getopt array definition>>
4391 tension_model_getopt_t *model = NULL;
4395 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
4396 tension_handler_data_t tdata;
4397 int num_param_args; /* for INIT_MODEL() */
4398 char **param_args; /* for INIT_MODEL() */
4400 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model,
4401 &Fmax, &Xmax, &flags);
4403 if (flags & VFLAG) {
4404 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4406 INIT_MODEL("utility", model, model->params, params);
4407 tdata.group_tension_params = NULL;
4408 push(&tdata.group_tension_params, params, 1);
4410 tdata.persist = NULL;
4411 if (model->handler == NULL) {
4412 printf("No tension function for model %s\n", model->name);
4418 printf("#Distance (m)\tForce (N)\n");
4419 for (i=0; i<=N; i++) {
4420 x = Xmax*i/(double)N;
4421 F = (*model->handler)(x, &tdata);
4422 if (F < 0 || F > Fmax) break;
4423 printf("%g\t%g\n", x, F);
4425 if (flags & VFLAG && i <= N) {
4426 /* explain exit condition */
4428 printf("Impossible force %g\n", F);
4430 printf("Reached large force limit %g > %g\n", F, Fmax);
4433 params = pop(&tdata.group_tension_params);
4434 if (model->destructor)
4435 (*model->destructor)(params);
4440 <<tension model utility includes>>=
4443 #include <string.h> /* strlen() */
4444 #include <assert.h> /* assert() */
4445 #include <errno.h> /* errno, ERANGE */
4449 #include "tension_model.h"
4452 <<tension model utility definitions>>=
4453 <<version definition>>
4454 #define VFLAG 1 /* verbose */
4455 <<string comparison definition and globals>>
4456 <<tension model getopt definitions>>
4457 <<initialize model definition>>
4461 <<tension model utility getopt functions>>=
4463 <<index tension model>>
4464 <<tension model utility help>>
4465 <<tension model utility get options>>
4468 <<tension model utility help>>=
4469 void help(char *prog_name,
4471 int n_tension_models, tension_model_getopt_t *tension_models,
4472 int tension_model, double Fmax, double Xmax)
4475 printf("usage: %s [options]\n", prog_name);
4476 printf("Version %s\n\n", VERSION);
4477 printf("A utility for understanding the available tension models\n\n");
4478 printf("Environmental options:\n");
4479 printf("-T T\tTemperature T (currently %g K)\n", env->T);
4480 printf("-C T\tYou can also set the temperature T in Celsius\n");
4481 printf("Model options:\n");
4482 printf("-m model\tFolded tension model (currently %s)\n",
4483 tension_models[tension_model].name);
4484 printf("-a args\tFolded tension model argument string (currently %s)\n",
4485 tension_models[tension_model].params);
4486 printf("F(x) is calculated for a range of x and printed\n");
4487 printf("For example:\n");
4488 printf(" #Distance (m)\tForce (N)\n");
4489 printf(" 123.456\t7.89\n");
4491 printf("-F\tSet the maximum F value for the standard mode F(x) output (currently %g)\n", Fmax);
4492 printf("-X\tSet the maximum x value for the standard mode F(x) output (currently %g)\n", Xmax);
4493 printf("-V\tChange output to verbose mode\n");
4494 printf("-h\tPrint this help and exit\n");
4496 printf("Tension models:\n");
4497 for (i=0; i<n_tension_models; i++) {
4498 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
4499 for (j=0; j < tension_models[i].num_params ; j++)
4500 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
4501 printf("\t\tdefaults: %s\n", tension_models[i].params);
4506 <<tension model utility get options>>=
4507 void get_options(int argc, char **argv, environment_t *env,
4508 int n_tension_models, tension_model_getopt_t *tension_models,
4509 tension_model_getopt_t **model, double *Fmax, double *Xmax,
4510 unsigned int *flags)
4512 char *prog_name = NULL;
4513 char c, options[] = "T:C:m:a:F:X:Vh";
4514 int tension_model=0, num_strings;
4515 extern char *optarg;
4516 extern int optind, optopt, opterr;
4520 /* setup defaults */
4522 prog_name = argv[0];
4523 env->T = 300.0; /* K */
4524 *Fmax = 1e5; /* N */
4525 *Xmax = 1e-6; /* m */
4527 *model = tension_models;
4529 while ((c=getopt(argc, argv, options)) != -1) {
4531 case 'T': env->T = safe_strtod(optarg, "-T"); break;
4532 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
4534 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
4535 *model = tension_models+tension_model;
4538 tension_models[tension_model].params = optarg;
4540 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
4541 case 'X': *Xmax = safe_strtod(optarg, "-X"); break;
4542 case 'V': *flags |= VFLAG; break;
4544 fprintf(stderr, "unrecognized option '%c'\n", optopt);
4545 /* fall through to default case */
4547 help(prog_name, env, n_tension_models, tension_models, tension_model, *Fmax, *Xmax);
4556 \section{Unfolding rate models}
4557 \label{sec.k_models}
4559 We have two main choices for determining $k$: Bell's model, which assumes the
4560 folded domains exist in a single `folded' state and make instantaneous,
4561 stochastic transitions over a free energy barrier; and Kramers' model, which
4562 treats the unfolding as a thermalized diffusion process.
4563 We deal with these options like we dealt with the different tension models: by bundling all model-specific
4564 parameters into structures, with model specific functions for determining $k$.
4566 <<k func definition>>=
4567 typedef double k_func_t(double F, environment_t *env, void *params);
4570 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
4571 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]].
4573 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
4574 because \LaTeX\ doesn't like underscores outside math mode.
4575 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
4576 You could use spaces instead (`\verb+ +'),
4577 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
4578 which I don't like as much.
4584 <<k func definition>>
4585 <<null k declarations>>
4586 <<const k declarations>>
4587 <<bell k declarations>>
4588 <<kbell k declarations>>
4589 <<kramers k declarations>>
4590 <<kramers integ k declarations>>
4591 #endif /* K_MODEL_H */
4594 <<k model module makefile lines>>=
4595 NW_SAWSIM_MODS += k_model
4596 CHECK_BINS += check_k_model
4597 check_k_model_MODS = parse string_eval
4599 [[check_k_model]] also depends on the parse module.
4601 <<k model utils makefile lines>>=
4602 K_MODEL_MODS = k_model parse string_eval
4603 K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
4604 $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h)
4605 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
4606 notangle -Rk-model-utils.c $< > $@
4607 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
4608 gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
4609 clean_k_model_utils :
4610 rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
4616 For domains that are never folded.
4618 <<null k declarations>>=
4620 <<null k functions>>=
4625 <<null k model getopt>>=
4626 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
4629 \subsection{Constant rate model}
4632 This is a pretty boring model, but it is usefull for testing the engine.
4636 where $k_0$ is the constant unfolding rate.
4638 <<const k functions>>=
4639 <<const k function>>
4640 <<const k structure create/destroy functions>>
4643 <<const k declarations>>=
4644 <<const k function declaration>>
4645 <<const k structure create/destroy function declarations>>
4646 <<const k global declarations>>
4650 <<const k structure>>=
4651 typedef struct const_k_param_struct {
4652 double knot; /* transition rate k_0 (frac population per s) */
4656 <<const k function declaration>>=
4657 double const_k(double F, environment_t *env, void *const_k_params);
4659 <<const k function>>=
4660 double const_k(double F, environment_t *env, void *const_k_params)
4661 { /* Returns the rate constant k in frac pop per s. */
4662 const_k_param_t *p = (const_k_param_t *) const_k_params;
4664 assert(p->knot > 0);
4669 <<const k structure create/destroy function declarations>>=
4670 void *string_create_const_k_param_t(char **param_strings);
4671 void destroy_const_k_param_t(void *p);
4674 <<const k structure create/destroy functions>>=
4675 const_k_param_t *create_const_k_param_t(double knot)
4677 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
4678 assert(ret != NULL);
4683 void *string_create_const_k_param_t(char **param_strings)
4685 return create_const_k_param_t(safe_strtod(param_strings[0], __FUNCTION__));
4688 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
4695 <<const k global declarations>>=
4696 extern int num_const_k_params;
4697 extern char *const_k_param_descriptions[];
4698 extern char const_k_param_string[];
4701 <<const k globals>>=
4702 int num_const_k_params = 1;
4703 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
4704 char const_k_param_string[]="1";
4707 <<const k model getopt>>=
4708 {"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}
4711 \subsection{Bell's model}
4714 Quantitatively, Bell's model gives $k$ as
4716 k = k_0 \cdot e^{F dx / k_B T} \;,
4718 where $k_0$ is the unforced unfolding rate,
4719 $F$ is the applied unfolding force,
4720 $dx$ is the distance to the transition state, and
4721 $k_B T$ is the thermal energy\citep{schlierf06}.
4723 <<bell k functions>>=
4725 <<bell k structure create/destroy functions>>
4728 <<bell k declarations>>=
4729 <<bell k function declaration>>
4730 <<bell k structure create/destroy function declarations>>
4731 <<bell k global declarations>>
4735 <<bell k structure>>=
4736 typedef struct bell_param_struct {
4737 double knot; /* transition rate k_0 (frac population per s) */
4738 double dx; /* `distance to transition state' (nm) */
4742 <<bell k function declaration>>=
4743 double bell_k(double F, environment_t *env, void *bell_params);
4745 <<bell k function>>=
4746 double bell_k(double F, environment_t *env, void *bell_params)
4747 { /* Returns the rate constant k in frac pop per s.
4748 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4749 * uses global k_B in J/K */
4750 bell_param_t *p = (bell_param_t *) bell_params;
4751 assert(F >= 0); assert(env->T > 0);
4753 assert(p->knot > 0); //assert(p->dx >= 0); /* p < 0 for refolding reactions, etc. */
4755 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
4756 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
4757 p->knot * exp(F*p->dx / (k_B*env->T) ));
4758 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
4760 return p->knot * exp(F*p->dx / (k_B*env->T) );
4764 <<bell k structure create/destroy function declarations>>=
4765 void *string_create_bell_param_t(char **param_strings);
4766 void destroy_bell_param_t(void *p);
4769 <<bell k structure create/destroy functions>>=
4770 bell_param_t *create_bell_param_t(double knot, double dx)
4772 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
4773 assert(ret != NULL);
4779 void *string_create_bell_param_t(char **param_strings)
4781 return create_bell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4782 safe_strtod(param_strings[1], __FUNCTION__));
4785 void destroy_bell_param_t(void *p /* bell_param_t * */)
4792 <<bell k global declarations>>=
4793 extern int num_bell_params;
4794 extern char *bell_param_descriptions[];
4795 extern char bell_param_string[];
4799 int num_bell_params = 2;
4800 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4801 char bell_param_string[]="3.3e-4,0.25e-9";
4804 <<bell k model getopt>>=
4805 {"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}
4807 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4808 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4811 \subsection{Linker stiffness corrected Bell model}
4814 Walton et. al showed that the Bell model constant force approximation
4815 doesn't hold for biotin-streptavadin binding, and I extended their
4816 results to I27 for the 2009 Biophysical Society Annual
4817 Meeting\cite{walton08,king09}. More details to follow when I get done
4818 with the conference\ldots
4820 We adjust Bell's model to
4822 k = k_0 \cdot e^{\frac{F dx - \frac{1}{2} \kappa dx^2}{k_B T}} \;,
4824 where $k_0$ is the unforced unfolding rate,
4825 $F$ is the applied unfolding force,
4826 $\kappa$ is the effective spring constant,
4827 $dx$ is the distance to the transition state, and
4828 $k_B T$ is the thermal energy\citep{schlierf06}.
4830 <<kbell k functions>>=
4831 <<kbell k function>>
4832 <<kbell k structure create/destroy functions>>
4835 <<kbell k declarations>>=
4836 <<kbell k function declaration>>
4837 <<kbell k structure create/destroy function declarations>>
4838 <<kbell k global declarations>>
4842 <<kbell k structure>>=
4843 typedef struct kbell_param_struct {
4844 double knot; /* transition rate k_0 (frac population per s) */
4845 double dx; /* `distance to transition state' (nm) */
4849 <<kbell k function declaration>>=
4850 double kbell_k(double F, environment_t *env, void *kbell_params);
4852 <<kbell k function>>=
4853 double kbell_k(double F, environment_t *env, void *kbell_params)
4854 { /* Returns the rate constant k in frac pop per s.
4855 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4856 * uses global k_B in J/K */
4857 kbell_param_t *p = (kbell_param_t *) kbell_params;
4858 assert(F >= 0); assert(env->T > 0); assert(env->stiffness >= 0);
4860 assert(p->knot > 0); assert(p->dx >= 0);
4861 return p->knot * exp((F*p->dx - 0.5*env->stiffness*(p->dx)*(p->dx))/(k_B*env->T) );
4865 <<kbell k structure create/destroy function declarations>>=
4866 void *string_create_kbell_param_t(char **param_strings);
4867 void destroy_kbell_param_t(void *p);
4870 <<kbell k structure create/destroy functions>>=
4871 kbell_param_t *create_kbell_param_t(double knot, double dx)
4873 kbell_param_t *ret = (kbell_param_t *) malloc(sizeof(kbell_param_t));
4874 assert(ret != NULL);
4880 void *string_create_kbell_param_t(char **param_strings)
4882 return create_kbell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4883 safe_strtod(param_strings[1], __FUNCTION__));
4886 void destroy_kbell_param_t(void *p /* kbell_param_t * */)
4893 <<kbell k global declarations>>=
4894 extern int num_kbell_params;
4895 extern char *kbell_param_descriptions[];
4896 extern char kbell_param_string[];
4899 <<kbell k globals>>=
4900 int num_kbell_params = 2;
4901 char *kbell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4902 char kbell_param_string[]="3.3e-4,0.25e-9";
4905 <<kbell k model getopt>>=
4906 {"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}
4908 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4909 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4912 \subsection{Kramer's model}
4915 Kramer's model gives $k$ as
4917 % \frac{1}{k} = \frac{1}{D}
4918 % \int_{x_\text{min}}^{x_\text{max}}
4919 % e^{\frac{-U_F(x)}{k_B T}}
4920 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4923 %where $D$ is the diffusion constant,
4924 %$U_F(x)$ is the free energy along the unfolding coordinate, and
4925 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4926 % before the minimum and after the maximum, respectively \citep{schlierf06}.
4928 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
4930 where $D$ is the diffusion constant,
4932 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
4934 are length scales where
4935 $x_c(F)$ is the location of the bound state and
4936 $x_{ts}(F)$ is the location of the transition state,
4937 $E(F,x)$ is the free energy, and
4938 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
4939 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
4940 \citet{evans97} Eqn.~3).
4942 <<kramers k functions>>=
4943 <<kramers k function>>
4944 <<kramers k structure create/destroy functions>>
4947 <<kramers k declarations>>=
4948 <<kramers k function declaration>>
4949 <<kramers k structure create/destroy function declarations>>
4950 <<kramers k global declarations>>
4954 <<kramers k structure>>=
4955 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
4956 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
4958 typedef struct kramers_param_struct {
4959 double D; /* diffusion rate (um/s) */
4966 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
4967 //kramers_x_func_t *xts; /* function returning transition x (nm) */
4968 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
4969 //kramers_E_func_t *E; /* function returning E (J) */
4970 //double *E_params; /* parametrize the energy landscape */
4971 //int n_E_params; /* length of E_params array */
4975 <<kramers k function declaration>>=
4976 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
4977 extern double kramers_k(double F, environment_t *env, void *kramers_params);
4979 <<kramers k function>>=
4980 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
4982 static char input[160]={0};
4983 static char output[80]={0};
4984 /* setup the environment */
4985 fprintf(in, "F = %g; T = %g\n", F, T);
4986 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4987 string_eval(in, out, input, 80, output);
4988 return safe_strtod(output, __FUNCTION__);
4991 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
4993 static char input[160]={0};
4994 static char output[80]={0};
4995 /* setup the environment */
4996 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
4997 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4998 string_eval(in, out, input, 80, output);
4999 return safe_strtod(output, __FUNCTION__);
5002 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
5004 kramers_param_t *p = (kramers_param_t *) kramers_params;
5005 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
5008 double kramers_k(double F, environment_t *env, void *kramers_params)
5009 { /* Returns the rate constant k in frac pop per s.
5010 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5011 * uses global k_B in J/K */
5012 kramers_param_t *p = (kramers_param_t *) kramers_params;
5013 double kbT, xc, xts, lc, lts, Eb;
5014 assert(F >= 0); assert(env->T > 0);
5017 assert(p->xc != NULL); assert(p->xts != NULL);
5018 assert(p->ddEdxx != NULL); assert(p->E != NULL);
5019 assert(p->in != NULL); assert(p->out != NULL);
5021 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
5022 if (xc == -1.0) return -1.0; /* error (Too much force) */
5023 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
5024 if (xc == -1.0) return -1.0; /* error (Too much force) */
5025 lc = sqrt(2.0*M_PI*kbT /
5026 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
5027 lts = sqrt(-2.0*M_PI*kbT/
5028 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
5029 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
5030 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
5031 //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));
5032 return p->D/(lc*lts) * exp(-Eb/kbT);
5036 <<kramers k structure create/destroy function declarations>>=
5037 void *string_create_kramers_param_t(char **param_strings);
5038 void destroy_kramers_param_t(void *p);
5041 <<kramers k structure create/destroy functions>>=
5042 kramers_param_t *create_kramers_param_t(double D,
5043 char *xc_fn, char *xts_fn,
5044 char *E_fn, char *ddEdxx_fn,
5045 char *global_define_string)
5046 // kramers_x_func_t *xc, kramers_x_func_t *xts,
5047 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
5048 // double *E_params, int n_E_params)
5050 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
5051 assert(ret != NULL);
5052 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5053 assert(ret->xc != NULL);
5054 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5055 assert(ret->xts != NULL);
5056 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
5057 assert(ret->E != NULL);
5058 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
5059 assert(ret->ddEdxx != NULL);
5061 strcpy(ret->xc, xc_fn);
5062 strcpy(ret->xts, xts_fn);
5063 strcpy(ret->E, E_fn);
5064 strcpy(ret->ddEdxx, ddEdxx_fn);
5065 //ret->E_params = E_params;
5066 //ret->n_E_params = n_E_params;
5067 string_eval_setup(&ret->in, &ret->out);
5068 fprintf(ret->in, "kB = %g\n", k_B);
5069 fprintf(ret->in, "%s\n", global_define_string);
5073 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
5074 void *string_create_kramers_param_t(char **param_strings)
5076 return create_kramers_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5084 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
5086 kramers_param_t *pk = (kramers_param_t *)p;
5088 string_eval_teardown(&pk->in, &pk->out);
5090 // free(pk->E_params);
5096 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
5097 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.
5098 We follow \citet{shillcock98} and use
5100 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
5101 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
5104 where TODO, variable meanings.%\citep{shillcock98}.
5105 The first and second derivatives of $E(F,E_b,x,x_s)$ are
5107 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\\
5108 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
5110 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
5111 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
5112 The roots are therefor at
5114 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
5115 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
5116 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
5119 <<kramers k global declarations>>=
5120 extern int num_kramers_params;
5121 extern char *kramers_param_descriptions[];
5122 extern char kramers_param_string[];
5125 <<kramers k globals>>=
5126 int num_kramers_params = 6;
5127 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"};
5128 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)}";
5130 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
5131 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
5132 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
5133 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.
5134 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
5135 It works so long as [[val_1]] is not [[false]].
5137 <<kramers k model getopt>>=
5138 {"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}
5141 \subsection{Kramer's model (natural cubic splines)}
5142 \label{sec.kramers_integ}
5144 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
5145 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
5146 \citet{schlierf06} Eqn.~2)
5148 \frac{1}{k} = \frac{1}{D}
5149 \int_{x_\text{min}}^{x_\text{max}}
5150 e^{\frac{U_F(x)}{k_B T}}
5151 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5154 where $D$ is the diffusion constant,
5155 $U_F(x)$ is the free energy along the unfolding coordinate, and
5156 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
5157 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
5159 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
5160 interpolating between them with cubic splines.
5161 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
5163 <<kramers integ k functions>>=
5164 <<spline functions>>
5165 <<kramers integ A k functions>>
5166 <<kramers integ B k functions>>
5167 <<kramers integ k function>>
5168 <<kramers integ k structure create/destroy functions>>
5171 <<kramers integ k declarations>>=
5172 <<kramers integ k function declaration>>
5173 <<kramers integ k structure create/destroy function declarations>>
5174 <<kramers integ k global declarations>>
5178 <<kramers integ k structure>>=
5179 typedef double func_t(double x, void *params);
5180 typedef struct kramers_integ_param_struct {
5181 double D; /* diffusion rate (um/s) */
5182 double x_min; /* integration bounds */
5184 func_t *E; /* function returning E (J) */
5185 void *E_params; /* parametrize the energy landscape */
5186 destroy_data_func_t *destroy_E_params;
5188 double F; /* for passing into GSL evaluations */
5190 } kramers_integ_param_t;
5193 <<spline param structure>>=
5194 typedef struct spline_param_struct {
5195 int n_knots; /* natural cubic spline knots for U(x) */
5198 gsl_spline *spline; /* GSL spline parameters */
5199 gsl_interp_accel *acc; /* GSL spline acceleration data */
5203 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
5205 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5207 <<kramers integ A k functions>>=
5208 double kramers_integ_k_integrandA(double x, void *params)
5210 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5211 double U = (*p->E)(x, p->E_params);
5212 double Fx = p->F * x;
5213 double kbT = k_B * p->env->T;
5214 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
5215 return exp(-(U-Fx)/kbT);
5217 double kramers_integ_k_integralA(double x, void *params)
5219 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5221 double result, abserr;
5222 size_t neval; /* for qng */
5223 static gsl_integration_workspace *w = NULL;
5224 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5225 f.function = &kramers_integ_k_integrandA;
5227 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5228 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5229 //fprintf(stderr, "integralA = %g\n", result);
5235 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
5236 \int_{x_\text{min}}^{x_\text{max}}
5237 e^{\frac{U_F(x)}{k_B T}}
5238 \text{Integral}_A(x)
5241 <<kramers integ B k functions>>=
5242 double kramers_integ_k_integrandB(double x, void *params)
5244 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5245 double U = (*p->E)(x, p->E_params);
5246 double Fx = p->F * x;
5247 double kbT = k_B * p->env->T;
5248 double intA = kramers_integ_k_integralA(x, params);
5249 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
5250 return exp((U-Fx)/kbT)*intA;
5252 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
5254 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5256 double result, abserr;
5257 size_t neval; /* for qng */
5258 static gsl_integration_workspace *w = NULL;
5259 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5260 f.function = &kramers_integ_k_integrandB;
5264 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5265 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5266 //fprintf(stderr, "integralB = %g\n", result);
5271 With the integrals taken care of, evaluating $k$ is simply
5273 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
5275 <<kramers integ k function declaration>>=
5276 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
5278 <<kramers integ k function>>=
5279 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
5280 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
5281 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
5282 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
5286 <<kramers integ k structure create/destroy function declarations>>=
5287 void *string_create_kramers_integ_param_t(char **param_strings);
5288 void destroy_kramers_integ_param_t(void *p);
5291 <<kramers integ k structure create/destroy functions>>=
5292 kramers_integ_param_t *create_kramers_integ_param_t(double D,
5293 double xmin, double xmax, func_t *E,
5295 destroy_data_func_t *destroy_E_params)
5297 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
5298 assert(ret != NULL);
5303 ret->E_params = E_params;
5304 ret->destroy_E_params = destroy_E_params;
5308 void *string_create_kramers_integ_param_t(char **param_strings)
5310 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
5311 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
5312 void *E_params = string_create_spline_param_t(param_strings+1);
5313 return create_kramers_integ_param_t(
5314 safe_strtod(param_strings[0], __FUNCTION__),
5315 safe_strtod(param_strings[2], __FUNCTION__),
5316 safe_strtod(param_strings[3], __FUNCTION__),
5317 &spline_eval, E_params, destroy_spline_param_t);
5320 void destroy_kramers_integ_param_t(void *params)
5322 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5325 (*p->destroy_E_params)(p->E_params);
5331 Finally we have the GSL spline wrappers:
5332 <<spline functions>>=
5334 <<create/destroy spline>>
5337 <<spline function>>=
5338 double spline_eval(double x, void *spline_params)
5340 spline_param_t *p = (spline_param_t *)(spline_params);
5341 return gsl_spline_eval(p->spline, x, p->acc);
5345 <<create/destroy spline>>=
5346 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
5348 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
5349 assert(ret != NULL);
5350 ret->n_knots = n_knots;
5353 ret->acc = gsl_interp_accel_alloc();
5354 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
5355 gsl_spline_init(ret->spline, x, y, n_knots);
5359 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
5360 void *string_create_spline_param_t(char **param_strings)
5364 double *x=NULL, *y=NULL;
5365 /* break into ordered pairs */
5366 parse_list_string(param_strings[0], SEP, '(', ')',
5367 &num_knots, &knot_coords);
5368 x = (double *)malloc(sizeof(double)*num_knots);
5370 y = (double *)malloc(sizeof(double)*num_knots);
5372 for (i=0; i<num_knots; i++) {
5373 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
5374 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
5377 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
5379 free_parsed_list(num_knots, knot_coords);
5380 return create_spline_param_t(num_knots, x, y);
5383 void destroy_spline_param_t(void *params)
5385 spline_param_t *p = (spline_param_t *)params;
5388 gsl_spline_free(p->spline);
5390 gsl_interp_accel_free(p->acc);
5400 <<kramers integ k global declarations>>=
5401 extern int num_kramers_integ_params;
5402 extern char *kramers_integ_param_descriptions[];
5403 extern char kramers_integ_param_string[];
5406 <<kramers integ k globals>>=
5407 int num_kramers_integ_params = 4;
5408 char *kramers_integ_param_descriptions[] = {
5409 "diffusion D, m^2 / s",
5410 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
5411 "starting integration bound x_min, m",
5412 "ending integration bound x_max, m"};
5413 //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";
5414 //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";
5415 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";
5416 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
5417 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
5418 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
5421 <<kramers integ k model getopt>>=
5422 {"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}
5424 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
5425 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5427 \subsection{Unfolding model implementation}
5431 <<k model includes>>
5432 #include "k_model.h"
5433 <<k model internal definitions>>
5435 <<k model functions>>
5438 <<k model includes>>=
5439 #include <assert.h> /* assert() */
5440 #include <stdlib.h> /* NULL, malloc(), strto*() */
5441 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
5442 #include <stdio.h> /* fprintf(), stdout */
5443 #include <string.h> /* strlen(), strcpy() */
5444 #include <errno.h> /* errno, ERANGE */
5445 #include <gsl/gsl_integration.h>
5446 #include <gsl/gsl_spline.h>
5451 <<k model internal definitions>>=
5452 <<const k structure>>
5453 <<bell k structure>>
5454 <<kbell k structure>>
5455 <<kramers k structure>>
5456 <<kramers integ k structure>>
5457 <<spline param structure>>
5460 <<k model globals>>=
5465 <<kramers k globals>>
5466 <<kramers integ k globals>>
5469 <<k model functions>>=
5471 <<null k functions>>
5472 <<const k functions>>
5473 <<bell k functions>>
5474 <<kbell k functions>>
5475 <<kramers k functions>>
5476 <<kramers integ k functions>>
5479 \subsection{Unfolding model unit tests}
5481 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
5482 <<check-k-model.c>>=
5483 <<k model check includes>>
5484 <<check relative error>>
5486 <<k model test suite>>
5487 <<main check program>>
5490 <<k model check includes>>=
5491 #include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
5492 #include <stdio.h> /* sprintf() */
5493 #include <assert.h> /* assert() */
5494 #include <math.h> /* exp() */
5495 #include <errno.h> /* errno, ERANGE */
5498 #include "k_model.h"
5501 <<k model test suite>>=
5505 <<k model suite definition>>
5508 <<k model suite definition>>=
5509 Suite *test_suite (void)
5511 Suite *s = suite_create ("k model");
5512 <<const k test case defs>>
5513 <<bell k test case defs>>
5515 <<const k test case adds>>
5516 <<bell k test case adds>>
5521 \subsubsection{Constant}
5523 <<const k test case defs>>=
5524 TCase *tc_const_k = tcase_create("Constant k");
5527 <<const k test case adds>>=
5528 tcase_add_test(tc_const_k, test_const_k_create_destroy);
5529 tcase_add_test(tc_const_k, test_const_k_over_range);
5530 suite_add_tcase(s, tc_const_k);
5535 START_TEST(test_const_k_create_destroy)
5537 char *knot[] = {"1","2","3","4","5","6"};
5538 char *params[] = {knot[0]};
5541 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5542 params[0] = knot[i];
5543 p = string_create_const_k_param_t(params);
5544 destroy_const_k_param_t(p);
5549 START_TEST(test_const_k_over_range)
5551 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
5552 char *knot[] = {"1","2","3","4","5","6"};
5553 char *params[] = {knot[0]};
5556 char param_string[80];
5559 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5560 params[0] = knot[i];
5561 p = string_create_const_k_param_t(params);
5562 for ( F=Fm; F<FM; F+=dF ) {
5563 fail_unless(const_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
5564 "const_k(%g, %g, &{%s}) = %g != %s",
5565 F, T, knot[i], const_k(F, &env, p), knot[i]);
5567 destroy_const_k_param_t(p);
5573 \subsubsection{Bell}
5575 <<bell k test case defs>>=
5576 TCase *tc_bell = tcase_create("Bell's k");
5579 <<bell k test case adds>>=
5580 tcase_add_test(tc_bell, test_bell_k_create_destroy);
5581 tcase_add_test(tc_bell, test_bell_k_at_zero);
5582 tcase_add_test(tc_bell, test_bell_k_at_one);
5583 suite_add_tcase(s, tc_bell);
5587 START_TEST(test_bell_k_create_destroy)
5589 char *knot[] = {"1","2","3","4","5","6"};
5591 char *params[] = {knot[0], dx};
5594 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5595 params[0] = knot[i];
5596 p = string_create_bell_param_t(params);
5597 destroy_bell_param_t(p);
5602 START_TEST(test_bell_k_at_zero)
5604 double F=0.0, T=300.0;
5605 char *knot[] = {"1","2","3","4","5","6"};
5607 char *params[] = {knot[0], dx};
5610 char param_string[80];
5613 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5614 params[0] = knot[i];
5615 p = string_create_bell_param_t(params);
5616 fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
5617 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
5618 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
5619 destroy_bell_param_t(p);
5624 START_TEST(test_bell_k_at_one)
5627 char *knot[] = {"1","2","3","4","5","6"};
5629 char *params[] = {knot[0], dx};
5630 double F= k_B*T/safe_strtod(dx, __FUNCTION__);
5633 char param_string[80];
5636 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5637 params[0] = knot[i];
5638 p = string_create_bell_param_t(params);
5639 CHECK_ERR(1e-6, safe_strtod(knot[i], __FUNCTION__)*exp(1.0), bell_k(F, &env, p));
5640 //fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__)*exp(1.0),
5641 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
5642 // F, T, knot[i], dx, bell_k(F, &env, p), safe_strtod(knot[i], __FUNCTION__)*exp(1.0));
5643 destroy_bell_param_t(p);
5649 <<kramers k tests>>=
5652 <<kramers k test case def>>=
5655 <<kramers k test case add>>=
5658 <<k model function tests>>=
5659 <<check relative error>>
5661 START_TEST(test_something)
5663 double k=5, x=3, last_x=2.0, F;
5664 one_dim_fn_t *handlers[] = {&hooke};
5665 void *data[] = {&k};
5666 double xi[] = {0.0};
5668 int new_active_groups = 1;
5669 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
5670 fail_unless(F = k*x, NULL);
5675 \subsection{Utilities}
5677 The unfolding models can get complicated, and are sometimes hard to explain to others.
5678 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
5679 the overhead of having to understand a full simulation.
5680 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
5681 or special mode, where the operation depends on the specific model selected.
5683 <<k-model-utils.c>>=
5685 <<k model utility includes>>
5686 <<k model utility definitions>>
5687 <<k model utility getopt functions>>
5688 <<k model utility multi model E>>
5689 <<k model utility main>>
5692 <<k model utility main>>=
5693 int main(int argc, char **argv)
5695 <<k model getopt array definition>>
5696 k_model_getopt_t *model = NULL;
5701 int num_param_args; /* for INIT_MODEL() */
5702 char **param_args; /* for INIT_MODEL() */
5704 double special_xmin;
5705 double special_xmax;
5706 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
5707 &Fmax, &special_xmin, &special_xmax, &flags);
5708 if (flags & VFLAG) {
5709 printf("#initializing model %s with parameters %s\n", model->name, model->params);
5711 INIT_MODEL("utility", model, model->params, params);
5714 if (model->k == NULL) {
5715 printf("No k function for model %s\n", model->name);
5719 printf("Fmax = %g <= 0\n", Fmax);
5725 printf("#F (N)\tk (%% pop. per s)\n");
5726 for (i=0; i<=N; i++) {
5727 F = Fmax*i/(double)N;
5728 k = (*model->k)(F, &env, params);
5730 printf("%g\t%g\n", F, k);
5735 if (model->k == NULL || model->k == &bell_k) {
5736 printf("No special mode for model %s\n", model->name);
5739 if (special_xmax <= special_xmin) {
5740 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
5744 double Xrng=(special_xmax-special_xmin), x, E;
5746 printf("#x (m)\tE (J)\n");
5747 for (i=0; i<=N; i++) {
5748 x = special_xmin + Xrng*i/(double)N;
5749 E = multi_model_E(model->k, params, &env, x);
5750 printf("%g\t%g\n", x, E);
5757 if (model->destructor)
5758 (*model->destructor)(params);
5763 <<k model utility multi model E>>=
5764 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
5766 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
5768 if (k == kramers_integ_k)
5769 E = (*p->E)(x, p->E_params);
5771 E = kramers_E(0, env, model_params, x);
5777 <<k model utility includes>>=
5780 #include <string.h> /* strlen() */
5781 #include <assert.h> /* assert() */
5782 #include <errno.h> /* errno, ERANGE */
5785 #include "string_eval.h"
5786 #include "k_model.h"
5789 <<k model utility definitions>>=
5790 <<version definition>>
5791 #define VFLAG 1 /* verbose */
5792 enum mode_t {M_K_OF_F, M_SPECIAL};
5793 <<string comparison definition and globals>>
5794 <<k model getopt definitions>>
5795 <<initialize model definition>>
5796 <<kramers k structure>>
5797 <<kramers integ k structure>>
5801 <<k model utility getopt functions>>=
5804 <<k model utility help>>
5805 <<k model utility get options>>
5808 <<k model utility help>>=
5809 void help(char *prog_name,
5811 int n_k_models, k_model_getopt_t *k_models,
5812 int k_model, double Fmax, double special_xmin, double special_xmax)
5815 printf("usage: %s [options]\n", prog_name);
5816 printf("Version %s\n\n", VERSION);
5817 printf("A utility for understanding the available unfolding models\n\n");
5818 printf("Environmental options:\n");
5819 printf("-T T\tTemperature T (currently %g K)\n", env->T);
5820 printf("-C T\tYou can also set the temperature T in Celsius\n");
5821 printf("Model options:\n");
5822 printf("-k model\tTransition rate model (currently %s)\n",
5823 k_models[k_model].name);
5824 printf("-K args\tTransition rate model argument string (currently %s)\n",
5825 k_models[k_model].params);
5826 printf("There are two output modes. In standard mode, k(F) is printed\n");
5827 printf("For example:\n");
5828 printf(" #Force (N)\tk (% pop. per s)\n");
5829 printf(" 123.456\t7.89\n");
5831 printf("In special mode, the output depends on the model.\n");
5832 printf("For models defining an energy function E(x), that function is printed\n");
5833 printf(" #Distance (m)\tE (J)\n");
5834 printf(" 0.0012\t0.0034\n");
5836 printf("-m\tChange output to standard mode\n");
5837 printf("-M\tChange output to special mode\n");
5838 printf("-F\tSet the maximum F value for the standard mode k(F) output (currently %g)\n", Fmax);
5839 printf("-x\tSet the minimum x value for the special mode E(x) output (currently %g)\n", special_xmin);
5840 printf("-X\tSet the maximum x value for the special mode E(x) output (currently %g)\n", special_xmax);
5841 printf("-V\tChange output to verbose mode\n");
5842 printf("-h\tPrint this help and exit\n");
5844 printf("Unfolding rate models:\n");
5845 for (i=0; i<n_k_models; i++) {
5846 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
5847 for (j=0; j < k_models[i].num_params ; j++)
5848 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
5849 printf("\t\tdefaults: %s\n", k_models[i].params);
5854 <<k model utility get options>>=
5855 void get_options(int argc, char **argv, environment_t *env,
5856 int n_k_models, k_model_getopt_t *k_models,
5857 enum mode_t *mode, k_model_getopt_t **model,
5858 double *Fmax, double *special_xmin, double *special_xmax,
5859 unsigned int *flags)
5861 char *prog_name = NULL;
5862 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
5864 extern char *optarg;
5865 extern int optind, optopt, opterr;
5869 /* setup defaults */
5871 prog_name = argv[0];
5872 env->T = 300.0; /* K */
5876 *Fmax = 1e-9; /* N */
5877 *special_xmax = 1e-8;
5878 *special_xmin = 0.0;
5880 while ((c=getopt(argc, argv, options)) != -1) {
5882 case 'T': env->T = safe_strtod(optarg, "-T"); break;
5883 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
5885 k_model = index_k_model(n_k_models, k_models, optarg);
5886 *model = k_models+k_model;
5889 k_models[k_model].params = optarg;
5891 case 'm': *mode = M_K_OF_F; break;
5892 case 'M': *mode = M_SPECIAL; break;
5893 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
5894 case 'x': *special_xmin = safe_strtod(optarg, "-x"); break;
5895 case 'X': *special_xmax = safe_strtod(optarg, "-X"); break;
5896 case 'V': *flags |= VFLAG; break;
5898 fprintf(stderr, "unrecognized option '%c'\n", optopt);
5899 /* fall through to default case */
5901 help(prog_name, env, n_k_models, k_models, k_model, *Fmax, *special_xmin, *special_xmax);
5910 \section{\LaTeX\ documentation}
5912 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
5913 The comment blocks are extracted (with nicely formatted code blocks), using
5914 <<latex makefile lines>>=
5915 $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5916 noweave -latex -index -delay $< > $@
5917 $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR)
5921 <<latex makefile lines>>=
5922 $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
5924 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5925 $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
5926 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5927 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5928 mv $(BUILD_DIR)/sawsim.pdf $@
5930 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
5931 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
5932 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
5934 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
5935 <<latex makefile lines>>=
5937 rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
5938 $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
5939 $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.toc \
5940 $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib
5941 clean_latex : semi-clean_latex
5942 rm -f $(DOC_DIR)/sawsim.pdf
5947 [[make]] is a common utility on *nix systems for managing dependencies while building software.
5948 In this case, we have several \emph{targets} that we'd like to build:
5949 the main simulation program \prog;
5950 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
5951 the documentation [[sawsim.pdf]];
5952 and the [[Makefile]] itself.
5953 Besides the generated files, there is the utility target
5954 [[clean]] for removing all generated files except [[Makefile]].
5956 % [[dist]] for generating a distributable tar file.
5958 Extract the makefile with
5959 `[[notangle -Rmakefile src/sawsim.nw | sed 's/ /\t/' > Makefile]]'.
5960 The sed is needed because notangle replaces tabs with eight spaces,
5961 but [[make]] needs tabs.
5963 # Decide what directories to put things in
5968 # Note: Cannot use, for example, `./build', because make eats the `./'
5969 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5971 # Modules (X.c, X.h) defined in the noweb file
5974 # Modules defined outside the noweb file
5975 FREE_SAWSIM_MODS = interp tavl
5977 <<list module makefile lines>>
5978 <<tension balance module makefile lines>>
5979 <<tension model module makefile lines>>
5980 <<k model module makefile lines>>
5981 <<parse module makefile lines>>
5982 <<string eval module makefile lines>>
5983 <<accel k module makefile lines>>
5985 SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS)
5987 # Everything needed for sawsim under one roof. sawsim.c must come first
5988 SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \
5989 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h)
5990 # Libraries needed to compile sawsim
5991 LIBS = gsl gslcblas m
5992 CHECK_LIBS = gsl gslcblas m check
5993 # The non-check binaries generated
5994 BINS = sawsim tension_model_utils k_model_utils sawsim_profile
5997 # Define the major targets
5998 all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
6000 view : $(DOC_DIR)/sawsim.pdf
6002 profile : $(BIN_DIR)/sawsim_profile | $(BIN_DIR)
6003 $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \
6004 -mnull -Mwlc -A0.39e-9,28e-9 -F8
6005 gprof $(BIN_DIR)/sawsim_profile gmon.out > $@
6006 check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim
6007 $(SHELL) -e -c 'for B in $^; do ./$$B; done'
6008 clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \
6009 clean_tension_model_utils \
6010 clean_k_model_utils clean_latex clean_check_sawsim
6011 rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \
6012 $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \
6013 $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
6014 $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
6015 $(BUILD_DIR)/global.h ./gmon.out
6016 $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
6018 # Various builds of sawsim
6019 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
6020 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6021 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
6022 gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6023 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
6024 gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6026 # Create the directories
6027 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
6030 # Copy over the source external to sawsim
6031 # Note: Cannot use, for example, `./build', because make eats the `./'
6032 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
6034 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6038 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6042 ## Basic source generated with noweb
6043 # The central files sawsim.c and global.h...
6044 $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6046 $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6047 notangle -Rglobal.h $< > $@
6048 # ... and the modules
6049 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6050 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6051 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6052 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6053 $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6054 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6055 # Note: I use `_' as a space in my file names, but noweb doesn't like
6056 # that so we use `-' to name the noweb roots and substitute here.
6058 ## Building the unit-test programs
6060 $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
6061 notangle -Rchecks $< > $@
6062 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
6063 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
6064 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR)
6065 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
6066 clean_check_sawsim :
6067 rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
6068 # ... and the modules
6070 $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
6071 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \
6072 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \
6073 $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6074 $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6075 $$(BUILD_DIR)/global.h | $(BIN_DIR)
6076 gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
6077 $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
6079 # todo: clean up the required modules to
6080 $(CHECK_BINS:%=clean_%) :
6081 rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
6083 # Cleaning up the modules
6085 $(SAWSIM_MODS:%=clean_%) :
6086 rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h)
6088 <<tension model utils makefile lines>>
6089 <<k model utils makefile lines>>
6090 <<latex makefile lines>>
6092 Makefile : $(SRC_DIR)/sawsim.nw
6093 notangle -Rmakefile $< | sed 's/ /\t/' > $@
6095 This is fairly self-explanatory. We compile a dynamically linked
6096 version ([[sawsim]]) and a statically linked version
6097 ([[sawsim_static]]). The static version is about 9 times as big, but
6098 it runs without needing dynamic GSL libraries to link to, so it's a
6099 better format for distributing to a cluster for parallel evaluation.
6103 \subsection{Hookean springs in series}
6104 \label{sec.math_hooke}
6107 The effective spring constant for $n$ springs $k_i$ in series is given by
6109 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6115 F &= k_i x_i &&\forall i \le n \\
6116 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
6117 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
6118 F &= k_1 x_1 = k_\text{eff} x \\
6119 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
6120 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6125 \addcontentsline{toc}{section}{References}
6126 \bibliography{sawsim}