1 %%%%%%%%%%%%%%%%%%%%%%%%% Start LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%
2 % we have our own preamble,
3 % use `noweave -delay` to not print noweb's automatic one
4 % -*- mode: Noweb; noweb-code-mode: c-mode -*-
5 \documentclass[letterpaper, 11pt]{article}
8 \noweboptions{smallcode}
10 \usepackage{url} % needed for natbib for underscores in urls?
11 \usepackage[breaklinks,colorlinks]{hyperref} % for internal links, \phantomsection
12 % breaklinks breaks long links
13 % colorlinks colors link text instead of boxing it
14 \usepackage{amsmath} % for \text{}, \xrightarrow[sub]{super}
15 \usepackage{amsthm} % for theorem and proof environments
16 \newtheorem{theorem}{Theorem}
18 \usepackage{subfig} % Support side-by-side figures
19 \usepackage{pgf} % Fancy graphics
20 \usepackage{tikz} % A nice, inline PGF frontend
21 \usetikzlibrary{automata} % Graph-theory library
23 \usepackage{doc} % BibTeX
24 \usepackage[super,sort&compress]{natbib} % fancy citation extensions
25 % super selects citations in superscript mode
26 % sort&compress automatically sorts and compresses compound citations (\cite{a,b,...})
28 \usepackage[mathletters]{ucs} % use math fonts to allow Greek UTF-8 in the text
29 \usepackage[utf8x]{inputenc} % allow UTF-8 characters
31 %\bibliographystyle{ieeetr} % pick the bibliography style, short and sweet
32 %\bibliographystyle{plain} % pick the bibliography style, includes dates
33 \bibliographystyle{plainnat}
34 \defcitealias{sw:noweb}{{\tt noweb}}
35 \defcitealias{galassi05}{GSL}
36 \defcitealias{sw:check}{{\tt check}}
37 \defcitealias{sw:python}{Python}
42 \textheight 9.5in % leave a bit of extra room for page numbers (at top of page)
45 \setlength{\parindent}{0pt}
46 \setlength{\parskip}{5pt}
48 % For some reason, the \maketitle wants to go on it's own page
49 % so we'll just hardcode the following in our first page.
50 %\title{Sawsim: a sawtooth protein unfolding simulator}
51 %\author{W.~Trevor King}
54 \newcommand{\prog}{[[sawsim]]}
55 \newcommand{\lang}{[[C]]}
56 \newcommand{\p}[3]{\left#1 #2 \right#3} % e.g. \p({complicated expr})
57 \newcommand{\Th}{\ensuremath{\sup{\text{th}}}} % ordinal th
58 \newcommand{\U}[1]{\ensuremath{\text{ #1}}} % units: $1\U{m/s}$
60 % single spaced lists, from Alvin Alexander
61 % http://www.devdaily.com/blog/post/latex/control-line-spacing-in-itemize-enumerate-tags/
62 \newenvironment{packed_item}{
64 \setlength{\itemsep}{1pt}
65 \setlength{\parskip}{0pt}
66 \setlength{\parsep}{0pt}
70 \nwfilename{sawsim.nw}
75 Sawsim: a sawtooth protein unfolding simulator \\
83 %%%%%%%%%%%%%%%%%%%%%%%%%%%% End LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%%
85 \section{Introduction}
87 The unfolding of multi-globular protein strings is a stochastic
88 process. In the \prog\ program, we use Monte Carlo methods to
89 simulate this unfolding through an explicit time-stepping approach.
90 We develop a program in \lang\ to simulate probability of unfolding
91 under a constant extension velocity of a chain of globular domains.
93 In order to extract physical meaning from many mechanical unfolding
94 simulations, it is necessary to compare the experimental results
95 against the results of simulations using different models and
96 parameters. The model and parameters that best match the experimental
97 data provide the ‘best’ model for that experiment.
99 Previous mechanical unfolding studies have rolled their own
100 simulations for modelling their experiment, and there has been little
101 exchange and standardization of these simulations between groups.
102 The \prog\ program attempts to fill this gap, providing a flexible,
103 extensible, \emph{community} program, to make it easier for groups
104 to share the burden and increase the transparency of data analysis.
106 ‘Sawsim’ is blend of ‘sawtooth simulation’.
108 \subsection{Related work}
112 \subsection{About this document}
114 This document was written in \citetalias{sw:noweb} with code blocks in
115 \lang\ and comment blocks in \LaTeX.
117 \subsection{System Requirements and Dependencies}
120 \subsection{Change Log}
122 Version 0.0 used only the unfolded domain WLC to determine the tension
123 and had a fixed timestep $dt$.
125 Version 0.1 added adaptive timesteps, adjusting so that the unfolding
126 probability for a given domain was constant.
128 Version 0.2 added Kramers' model unfolding rate calculations, and a
129 switch to select Bell's or Kramers' model. This induced a major
130 rewrite, introducing the tension group abstraction, and a split to
131 multiple \lang\ source files. Also added a unit-test suites for
132 sanity checking, and switched to SI units throughout.
134 Version 0.3 added integrated Kramers' model unfolding rate
135 calculations. Also replaced some of my hand-coded routines with GNU
136 Scientific Library \citepalias{galassi05} calls.
138 Version 0.4 added Python evaluation for the saddle-point Kramers'
139 model unfolding rate calculations, so that the model functions could
140 be controlled from the command line. Also added interpolating lookup
141 tables to accelerate slow unfolding rate model evaluation, and command
142 line control of tension groups.
144 Version 0.5 added the stiffness environmental parameter and it's
145 associated unfolding models.
147 Version 0.6 generalized from two state unfolding to arbitrary
148 $N$-state rate reactions. This allows simulations of
149 unfolding-refolding, metastable intermediates, etc., but required
150 reorganizing the command line interface.
152 Version 0.7 added tension model inverses, which seems to reduce
153 computation time by about a factor of 2. I was expecting more, but
154 I'll take what I can get.
156 Version 0.8 fixed a minor bug in unfolding probability for
157 multi-domain groups. The probability of at least one domain unfolding
158 had been calculated as $P_N=NP_1$, rather than $P_N=1-(1-P_1)^N$.
159 However, for small $P$ the two are equivalent.
161 Version 0.9 the -P option to sawsim now sets the target probability
162 for the per-state transition $P_N$, not the per-domain transisiton
165 <<version definition>>=
166 #define VERSION "0.9"
172 sawsim - program for simulating single molecule mechanical unfolding.
173 Copyright (C) 2008-2009, William Trevor King
175 This program is free software; you can redistribute it and/or
176 modify it under the terms of the GNU General Public License as
177 published by the Free Software Foundation; either version 3 of the
178 License, or (at your option) any later version.
180 This program is distributed in the hope that it will be useful, but
181 WITHOUT ANY WARRANTY; without even the implied warranty of
182 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
183 See the GNU General Public License for more details.
185 You should have received a copy of the GNU General Public License
186 along with this program; if not, write to the Free Software
187 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
190 The author may be contacted at <wking@drexel.edu> on the Internet, or
191 write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
192 Philadelphia PA 19104, USA.
205 The unfolding system is stored as a chain of \emph{domains} (Figure
206 \ref{fig.domain_chain}). Domains can be folded, globular protein
207 domains, unfolded protein linkers, AFM cantilevers, or other
208 stretchable system components. The system is modeled as graph of
209 possible domain states with transitions between them (Figure
210 \ref{fig.domain_states}).
214 \subfloat[][]{\label{fig.domain_chain}
215 \begin{tikzpicture}[->,node distance=2cm,font=\footnotesize]
216 \tikzstyle{every state}=[fill,draw=red!50,very thick,fill=red!20]
217 \node[state] (A) {domain 1};
218 \node[state] (B) [below of=A] {domain 2};
219 \node[state] (C) [below of=B] {.~.~.};
220 \node[state] (D) [below of=C] {domain $N$};
221 \node (S) [below of=D] {Surface};
222 \node (E) [above of=A] {};
224 \path[-] (A) edge (B)
225 (B) edge node [right] {Tension} (C)
228 \path[->,green] (A) edge node [right,black] {Extension} (E);
232 \subfloat[][]{\label{fig.domain_states}
233 \begin{tikzpicture}[->,node distance=2.5cm,shorten <=1pt,shorten >=1pt,font=\footnotesize]
234 \tikzstyle{every state}=[fill,draw=blue!50,very thick,fill=blue!20]
235 \node[state] (A) {cantilever};
236 \node[state] (C) [below of=A] {transition};
237 \node[state] (B) [left of=C] {folded};
238 \node[state] (D) [right of=C] {unfolded};
240 \path (B) edge [bend left] node [above] {$k_1$} (C)
241 (C) edge [bend left] node [below] {$k_1'$} (B)
242 edge [bend left] node [above] {$k_2$} (D)
243 (D) edge [bend left] node [below] {$k_2'$} (C);
246 \caption{\subref{fig.domain_chain} Extending a chain of domains. One end
247 of the chain is fixed, while the other is extended at a constant
248 speed. The domains are coupled with rigid linkers, so the domains
249 themselves must stretch to accomodate the extension.
250 \subref{fig.domain_states} Each domain exists in a discrete state. At
251 each timestep, it may transition into another state following a
252 user-defined state matrix such as this one, showing a metastable
253 transition state and an explicit ``cantilever'' domain.}
257 Each domain contributes to the total tension, which depends on the
258 chain extension. The particular model for the tension as a function
259 of extension is handled generally with function pointers. So far the
260 following models have been implemented:
262 \item Null (Appendix \ref{sec.null_tension}),
263 \item Constant (Appendix \ref{sec.const_tension}),
264 \item Hookean spring (Appendix \ref{sec.hooke}),
265 \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
266 \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
267 \item Piston (Appendix \ref{sec.piston_tension}),
270 The tension model and parameters used for a particular domain depend
271 on the domain's current state. The transition rates between states
272 are also handled generally with function pointers, with
275 \item Null (Appendix \ref{sec.null_k}),
276 \item Bell's model (Appendix \ref{sec.bell}),
277 \item Spring constant corrected Bell model (Appendix \ref{sec.kbell}),
278 \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
279 \item Kramers' saddle point model (Appendix \ref{sec.kramers}).
282 The unfolding of the chain domains is modeled in two stages. First
283 the tension of the chain is calculated. Then the tension (assumed to
284 be constant over the length of the chain, see Section
285 \ref{sec.timescales}), is applied to each domain, and used to
286 calculate the probability of that domain changing state in the next
287 timestep $dt$. Then the time is stepped forward, and the process is
288 repeated. The simulation is complete when a pre-set time limit
289 [[t_max]] is reached or a pre-selected state [[stop_state]] is empty.
291 int main(int argc, char **argv)
293 <<tension model getopt array definition>>
294 <<k model getopt array definition>>
296 /* variables initialized in get_options() */
297 list_t *state_list=NULL, *transition_list=NULL;
298 environment_t env={0};
299 double P, t_max, dt_max, v, x_step;
300 state_t *stop_state=NULL;
302 /* variables used in the simulation loop */
303 double t=0, x=0, dt, F; /* time, distance, timestep, force */
304 int transition=1; /* boolean marking a transition for tension optimization */
306 get_options(argc, argv, &P, &t_max, &dt_max, &v, &x_step, &stop_state, &env,
307 NUM_TENSION_MODELS, tension_models, NUM_K_MODELS, k_models,
308 &state_list, &transition_list);
311 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Distance (m)\tForce (N)\tStiffness (N/m)\n");
312 if (flags & FLAG_OUTPUT_TRANSITION_FORCES) printf("#Force (N)\tInitial state\tFinal state\n");
313 while ((t_max > 0 && t <= t_max) || (stop_state != NULL && stop_state->num_domains > 0)) {
314 F = find_tension(state_list, &env, x, transition);
316 dt = determine_dt(state_list, transition_list, &env, x, P, dt_max, v);
318 dt = determine_dt(state_list, transition_list, &env, x, P, dt_max, 0);
319 transition=random_transitions(transition_list, F, dt, &env);
320 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness);
325 if (dt == dt_max) { /* step completed */
328 } else { /* still working on this step */
333 destroy_state_list(state_list);
334 destroy_transition_list(transition_list);
338 @ The meat of the simulation is bundled into the three functions
339 [[find_tension]], [[determine_dt]], and [[random_transitions]].
340 [[find_tension]] is discussed in Section \ref{sec.find_tension},
341 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
342 [[random_transitions]] in Section \ref{sec.transition_rate}.
344 The stretched distance is extended in one of two modes depending on
345 [[xstep]]. If [[xstep == 0]], then the pulling is continuous, at
346 least within the limits of the inherent discretization of a time
347 stepped approach. At any rate, the timestep is limited by the maximum
348 allowed timestep [[dt_max]] and the maximum allowed unfolding
349 probability in a given step [[P]]. In the second mode [[xstep > 0]],
350 and the end to end distance increases in discrete steps of that
351 length. The time between steps is chosen to maintain an average
352 unfolding velocity [[v]]. This approach is less work to simulate than
353 smooth pulling and also closer to the reality of many experiments, but
354 it is practically impossible to treat analytically. With both pulling
355 styles implemented, the effects of the discretization can be easily
356 calculated, bridging the gap between theory and experiment.
358 Environmental parameters important in determining reaction rates and
359 tensions (e.g.\ temperature, pH) are stored in a single structure to
360 facilitate extension to more complicated models in the future.
361 <<environment definition>>=
362 typedef struct environment_struct {
372 #define DOUBLE_PRECISION 1e-12
374 <<environment definition>>
375 <<one dimensional function definition>>
376 <<create/destroy definitions>>
378 #endif /* GLOBAL_H */
380 Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
381 included multiple times.
383 \section{Simulation functions}
385 <<simulation functions>>=
389 <<does domain transition>>
390 <<randomly transition domains>>
394 \label{sec.find_tension}
396 Because the stretched system may be made up of several parts (folded
397 domains, unfolded domains, spring-like cantilever, \ldots), we will
398 assign the domains to states with tension models and groups. The
399 models are used to determine the tension of all the domain in a given
400 state following some model (Hookean spring, WLC, \ldots). The domains
401 are assumed to be commutative, so ordering is ignored. The
402 interactions between the groups are assumed to be linear, but the
403 interactions between domains of the same group need not be. Each
404 model has a tension handler function, which gives the tension $F$
405 needed to stretch a given group of domains a distance $x$.
407 Groups represent collections of states which the model should treat as
408 a single entity. To understand the purpose of groups, consider a
409 system with two unfolded domain states (e.g.\ for two different
410 proteins) that were both modeled as WLCs. If both unfolded states had
411 the same persistence length, it would be wise to assign them to the
412 same group. This would reduce both the computational cost of
413 balancing the tension and the error associated with the inter-group
414 interaction linearity. Note that group numbers only matter
415 \emph{within} model classes, since grouping domains with seperate
416 models doesn't make sense. Within designated groups, the tension
417 parameters for each domain are still checked for compatibility, so if
418 you accidentally assign incompatible domains to the same group you
419 should get an appropriate error message.
421 <<find tension definitions>>=
422 #define NUM_TENSION_MODELS 6
426 <<tension handler array global declaration>>=
427 extern one_dim_fn_t *tension_handlers[];
430 <<tension handler array global>>=
431 one_dim_fn_t *tension_handlers[] = {
433 &const_tension_handler,
441 <<tension model global declarations>>=
442 <<tension handler array global declaration>>
445 <<tension handler types>>=
446 typedef struct tension_handler_data_struct {
447 list_t *group_tension_params; /* one item for each domain in the group */
450 } tension_handler_data_t;
454 After sorting the chain into separate groups $G_i$, with tension
455 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
456 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
457 \forall i,j$. Note that there may be several states within each
458 group. [[find_tension]] is basically a wrapper to feed proper input
459 into the [[tension_balance]] function. [[transition]] should be set
460 to 0 if there were no transitions since the previous call to
461 [[find_tension]] to support some optimizations that assume a small
462 increase in tension between steps (only true if no transition has
463 occured). While we're messing with the tension handlers, we also grab
464 a measure of the current linker stiffness for the environmental
465 variable, which is needed by some of the unfolding rate models
466 (e.g.\ the linker-stiffness corrected Bell model, Appendix
469 double find_tension(list_t *state_list, environment_t *env, double x, int transition)
471 static int num_active_groups=0;
472 static one_dim_fn_t **per_group_handlers = NULL;
473 static one_dim_fn_t **per_group_inverse_handlers = NULL;
474 static void **per_group_data = NULL; /* array of pointers to tension_handler_data_t */
475 static double last_x;
476 tension_handler_data_t *tdata;
481 fprintf(stderr, "%s:\tfinding tension at distance %g%s\n", __FUNCTION__, x, transition ? " after a transition" : "");
484 if (transition != 0 || num_active_groups == 0) { /* setup statics */
485 /* get new statics, freeing the old ones if they aren't NULL */
486 get_active_groups(state_list, &num_active_groups, &per_group_handlers, &per_group_inverse_handlers, &per_group_data);
488 /* fill in the group handlers and params */
490 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]);
492 for (i=0; i<num_active_groups; i++) {
493 tdata = (tension_handler_data_t *) per_group_data[i];
495 fprintf(stderr, "%s:\tactive group %d of %d (data %p, param list %p)", __FUNCTION__, i, num_active_groups, tdata, tdata->group_tension_params);
496 list_print(stderr, tdata->group_tension_params, " parameters");
499 /* tdata->persist continues unchanged */
502 } /* else continue with the current statics */
504 F = tension_balance(num_active_groups, per_group_handlers, per_group_inverse_handlers, (void **)per_group_data, last_x, x, &env->stiffness);
510 @ For the moment, we will restrict our group boundaries to a single
511 dimension, so $\sum_i x_i = x$, our total extension (it is this
512 restriction that causes the groups to interact linearly). We'll also
513 restrict our extensions to all be positive. With these restrictions,
514 the problem of balancing the tensions reduces to solving a system of
515 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
516 the number of active groups. In general this can be fairly
517 complicated, but without much loss of practicality we can restrict
518 ourselves to strictly monotonically increasing, non-negative tension
519 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
520 simpler. For example, it guarantees the existence of a unique, real
521 solution for finite forces. See Appendix \ref{sec.tension_balance}
522 for the tension balancing implementation.
524 Each group must define a parameter structure if the tension model is
526 a function to create the parameter structure,
527 a function to destroy the parameter structure, and
528 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as its data argument.
529 The tension-specific parameter structures are contained in the domain
530 group field. For optimization, a group may also define an inverse
531 tension function as an optimization. See the Section
532 \ref{sec.model_selection} for a discussion on the command line
533 framework and inverse function discussion. See the worm-like chain
534 implementation for an example (Section \ref{sec.wlc}).
536 The major limitation of our approach is that all of our tension models
537 are Markovian (memory-less), which is not really a problem for the
538 chains (see \citet{evans99} for WLC thermalization rate calculations).
540 \subsection{Transition rate}
541 \label{sec.transition_rate}
543 Each state transition is modeled as unimolecular, first order reaction
545 \text{State 1} \xrightarrow{k} \text{State 2}
547 With the rate constant $k$ defined by
549 \frac{d[\text{State 2}]}{dt} = -\frac{d[\text{State 1}]}{dt} = k [\text{State 1}].
551 So the probability of a given domain transitioning in a short time
557 For a group of $N$ identical domains, and therefore identical
558 unfolding probabilities $P_1$, the probability of $n$ domains
559 unfolding in a given timestep is given by the binomial distribution
561 P(n) = \frac{N!}{n!(N-n)!}P_1^n(1-P_1)^(N-n).
563 The probability that \emph{none} of these domains unfold is then
567 so the probability that \emph{at least one} domain unfolds is
571 Note that for small $P$,
573 p(n>0) = 1-(1-NP_1+\mathcal{O}(P_1^2)) = NP_1 - \mathcal{O}(P^2)
576 This conversion from $P_1$ to $P_N$ is handled by the [[pN]] macro
578 /* find multi-domain transition rate from the single domain rate */
579 #define PN(P,N) (1.0-pow(1.0-(P), (N)))
583 We can also define a macro to find the approprate $dt$ to achieve a
584 target $P_N$ with a given $k$ and $N$.
586 P_N &= 1-(1-P_1)^N \\
587 1-P_1 &= (1-P_N)^{1/N} \\
588 P_1 &= 1-(1-P_N)^{1/N}
591 #define P1(PN,N) (1.0-pow(1.0-(PN), 1.0/(N)))
594 We take some time to discuss the meaning of $p(n>1)$
595 (i.e. multi-unfolding timesteps) in Section \ref{sec.timescales}.
597 <<does domain transition>>=
598 int domain_transitions(double F, double dt, environment_t *env, transition_t *transition)
599 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to transition structure */
601 k = accel_k(transition->k, F, env, transition->k_params);
602 //(*transition->k)(F, env, domain->k_params);
603 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
604 return happens(PN(k*dt, N_INIT(transition)));
606 @ [[happens]] is a random decision making function defined in Appendix \ref{sec.utils}.
608 <<randomly transition domains>>=
609 int random_transitions(list_t *transition_list, double tension,
610 double dt, environment_t *env)
611 { /* returns 1 if there was a transition and 0 otherwise */
612 while (transition_list != NULL) {
613 if (T(transition_list)->initial_state->num_domains > 0
614 && domain_transitions(tension, dt, env, T(transition_list))) {
615 if (flags & FLAG_OUTPUT_TRANSITION_FORCES)
616 fprintf(stdout, "%g\t%s\t%s\n", tension, T(transition_list)->initial_state->name, T(transition_list)->final_state->name);
617 T(transition_list)->initial_state->num_domains--;
618 T(transition_list)->final_state->num_domains++;
619 return 1; /* our one domain has transitioned, stop looking for others */
621 transition_list = transition_list->next;
627 \subsection{Timescales}
628 \label{sec.timescales}
630 The simulation assumes that chain equilibration is instantaneous
631 relative to the simulation timestep $dt$, so that tension is uniform
632 along the chain. The quality of this assumption depends on your
633 particular chain. For example, a damped spring thermalizes on a
634 timescale of order $\tau = 1/\gamma$ where the damping ratio $\gamma
635 \equiv \omega_0(-\zeta \pm \sqrt{\zeta^2-1}$, the natural angular
636 frequency $omega_0 \equiv \sqrt{k/m}$, $\zeta \equiv b/2\sqrt{km}$,
637 and $b$ sets the drag $F_b = -bv$. For our cantilevers $\tau$ is
638 on the order of milliseconds, which is longer than a timestep.
639 However, the approximation is still reasonable, because there is
640 rarely unfolding during the cantilever return. You could, of course,
641 take cantilever, WLC, etc.\ response times into effect, but that
642 would significantly complicate a simulation that seems to work
643 well enough without bothering :p. For WLC and FJC relaxation times,
644 see the Appendix of \citet{evans99} and \citet{kroy07}.
646 Besides assuming our timestep is much greater than equilibration
647 times, we also force it to be short enough so that only one domain may
648 unfold in a given timestep. For an unfolding timescale of $dt_u =
649 1/k$, we adapt our timesteps to keep $dt \ll dt_u$, so the probability
650 of two domains unfolding in the same timestep is small (see
651 [[determine_dt]] in Section \ref{sec.adaptive_dt} and [[P]] in
652 [[main()]] in Section \ref{sec.main} set by [[-P]] in
653 [[get_options()]] in Section \ref{sec.get_options}). This approach
654 breaks down as the adaptive timestep scheme approaches the transition
655 timestep $dt_t$, but $dt_t \sim 1\U{$\mu$s}$ for Ig27-like proteins
656 \citep{klimov00}, so this shouldn't be a problem. To reassure
657 yourself, you can ask the simulator to print the smallest timestep
658 that was used with TODO.
660 Even if the likelyhood of two domains unfolding in the same timestep
661 is small, if you run long enough simulations it will eventually occur.
662 In this case, we assume that $dt_t \ll dt$, so even if two domains
663 would unfold if we stuck to our unfolding rate $k$ for an entire
664 timestep $dt$, in ``reality'' only one of those domains unfolds.
665 Because we do not know when in the timestep the unfolding took place,
666 we move on to the next timestep without checking for further
667 unfoldings. This ``unchecked timestep portion'' should not contribute
668 significant errors to the calculation, because if $dt$ was small
669 enough that unfolding was unlikely at the pre-unfolding force, it
670 should be even more unlikely at the post-unfolding force, especially
671 over only a fraction of the timestep. The only dangerous cases for
672 the current implementation are those where unfolding intermediates are
673 much less stable than their precursor states, in which case an
674 unfolding event during the remaining timestep may actually be likely.
675 A simple workaround for such cases is to pick the value for [[P]] to
676 be small enough that the $dt \ll dt_u$ assumption survives even under
677 a transition populating the unstable intermediate.
679 \subsection{Adaptive timesteps}
680 \label{sec.adaptive_dt}
682 We'd like to pick $dt$ so the probability of changing state
683 (e.g.\ unfolding) in the next timestep is small. If we don't adapt our
684 timestep, we also risk breaking our approximation that $F(x' \in [x,
685 x+v \cdot dt]) = F(x)$ or that only one domain changes state in a
686 given timestep. The simulation would have been accurate for
687 sufficiently small fixed timesteps, but adaptive timesteps will allow
688 us to move through low tension regions in fewer steps, leading to a
689 more efficient simulation.
691 Consider the two state folded $\rightarrow$ unfolded transition.
692 Because $F(x)$ is monotonically increasing between unfoldings,
693 excessively large timesteps will lead to erroneously small unfolding
694 rates (an thus, higher average unfolding force).
696 The actual adaptive timestep implementation is not particularly
697 interesting, since we are only required to reduce $dt$ to somewhere
698 below a set threshold, so I've removed it to Appendix
699 \ref{sec.adaptive_dt_implementation}.
705 The models provide the physics of the simulation, but the simulation
706 allows interchangeable models, and we are currently focusing on the
707 simulation itself, so we remove the actual model implementation to the
710 The tension models are defined in Appendix \ref{sec.tension_models}
711 and unfolding models are defined in Appendix \ref{sec.k_models}.
714 #define k_B 1.3806503e-23 /* J/K */
718 \section{Command line interface}
720 <<get option functions>>=
722 <<index tension model>>
727 \subsection{Model selection}
728 \label{sec.model_selection}
730 The main difficulty with the command line interface in \prog\ is
731 developing an intuitive method for accessing the possibly large number
732 of available models. We'll treat this generally by defining an array
733 of available models, containing their important parameters.
735 <<get options definitions>>=
736 <<model getopt definitions>>
739 <<create/destroy definitions>>=
740 typedef void *create_data_func_t(char **param_strings);
741 typedef void destroy_data_func_t(void *param_struct);
744 <<model getopt definitions>>=
745 <<tension model getopt definitions>>
746 <<k model getopt definitions>>
749 In order to choose models, we need a method of assembling a reaction
750 graph such as that in Figure \ref{fig.domain_states} and population
751 the graph with a set of domains. First we declare the domain states
754 [[-s <state-name>,<model>[,<group>] -S <model-params> [-N <initial-population]]] \\
758 [[-s unfolded,wlc,1 -S 1e-8,4e-10 -N 5]]
760 This creates the state named [[unfolded]], using the WLC model and one
761 as a group number (defaults to 0, see Section \ref{sec.find_tension}).
762 The tension parameters are then set to [[1e-8,4e-10]], which in the
763 case of the WLC are the contour and persistence lengths respectively.
764 Finally we populate the state by adding five domains to the state.
765 The name of the state is importent for identifying states later.
767 Once the domains have been created and populated, you can added
768 transition rates following
770 [[-k <initial-state>,<final-state>,<model> -K <model-params>]] \\
774 [[-k folded,unfolded,bell -K 3.3e-4,0.25e-9]]
776 This creates a reaction going from the [[folded]] state to the
777 [[unfolded]] state, with the rate constant given by the Bell model
778 (Appendix \ref{sec.bell}). The unfolding parameters are then set to
779 [[3.3e-4,0.25e-9]], which in the case of the Bell model are the
780 unforced rate and transition state distance respectively.
782 \subsubsection{Tension}
784 To access the various tension models from the command line, we define
785 a structure that contains all the neccessary information about the
787 <<tension model getopt definitions>>=
788 typedef struct tension_model_getopt_struct {
789 one_dim_fn_t *handler; /* fn implementing the model on a group */
790 one_dim_fn_t *inverse_handler; /* fn implementing inverse on group */
791 char *name; /* name string identifying the model */
792 char *description; /* a brief description of the model */
793 int num_params; /* number of extra params for the model */
794 char **param_descriptions; /* descriptions of extra params */
795 char *params; /* default values for the extra params */
796 create_data_func_t *creator; /* param-string -> model-data-struct fn */
797 destroy_data_func_t *destructor; /* fn to free the model data structure */
798 } tension_model_getopt_t;
799 @ Set [[inverse_handler]] to [[NULL]] if you wish to invert using
800 bisection. Obviously, this will be slower than computing an
801 analytical inverse and more prone to rounding errors, but it may be
802 easier if a closed-form inverse does not exist.
804 <<tension model getopt array definition>>=
805 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
806 <<null tension model getopt>>,
807 <<constant tension model getopt>>,
808 <<hooke tension model getopt>>,
809 <<worm-like chain tension model getopt>>,
810 <<freely-jointed chain tension model getopt>>,
811 <<piston tension model getopt>>
815 \subsubsection{Unfolding rate}
817 <<k model getopt definitions>>=
818 #define NUM_K_MODELS 6
820 typedef struct k_model_getopt_struct {
825 char **param_descriptions;
827 create_data_func_t *creator;
828 destroy_data_func_t *destructor;
832 <<k model getopt array definition>>=
833 k_model_getopt_t k_models[NUM_K_MODELS] = {
834 <<null k model getopt>>,
835 <<const k model getopt>>,
836 <<bell k model getopt>>,
837 <<kbell k model getopt>>,
838 <<kramers k model getopt>>,
839 <<kramers integ k model getopt>>
846 void help(char *prog_name, double P, double t_max, double dt_max, double v, double x_step,
847 state_t *stop_state, environment_t *env,
848 int n_tension_models, tension_model_getopt_t *tension_models,
849 int n_k_models, k_model_getopt_t *k_models,
850 int k_model, list_t *state_list)
855 if (state_list != NULL)
856 state = S(tail(state_list));
858 printf("usage: %s [options]\n", prog_name);
859 printf("Version %s\n\n", VERSION);
860 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
862 printf("Simulation options:\n");
864 printf("-q state\tStop simulation when this state empties (- for no limit, currently %s)\n", stop_state ? stop_state->name : "-" );
865 printf("-t tmax\tMaximum simulated time (-1 for no limit, currently %g s)\n", t_max);
866 printf("-d dt\tMaximum allowed timestep dt (currently %g s)\n", dt_max);
867 printf("-P P\tTarget probability for dt (currently %g)\n", P);
868 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
869 printf("-x dx\tPulling step size (currently %g m)\n", x_step);
870 printf("\tWhen dx = 0, the pulling is continuous.\n");
871 printf("\tWhen dx > 0, the pulling occurs in discrete steps.\n");
873 printf("Environmental options:\n");
874 printf("-T T\tTemperature T (currently %g K)\n", env->T);
875 printf("-C T\tYou can also set the temperature T in Celsius\n");
877 printf("State creation options:\n");
878 printf("-s name,model[[,group],params]\tCreate a new state.\n");
879 printf("Once you've set up the state, you may populate it with domains\n");
880 printf("-N n\tAdd n domains in the current state (%s)\n", state ? state->name : "-");
882 printf("Transition creation options:\n");
883 printf("-k initial_state,final_state,model[,params]\tCreate a new transition\n");
885 printf("To double check your reaction graph, you can produce graphviz dot output\n");
886 printf("describing the current states and transitions.\n");
887 printf("-D\tPrint dot description of states and transitions and exit.\n");
889 printf("Simulation output mode options:\n");
890 printf("There are two output modes. In standard mode, only the transition\n");
891 printf("events are printed. For example:\n");
892 printf(" #Force (N)\tInitial state\tFinal state\n");
893 printf(" 123.456e-12\tfolded\tunfolded\n");
895 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
896 printf(" #Distance (m)\tForce (N)\tStiffness (N/m)\n");
897 printf(" 0.001\t0.0023\t0.002\n");
899 printf("-V\tChange output to verbose mode.\n");
900 printf("-h\tPrint this help and exit.\n");
903 printf("Tension models:\n");
904 for (i=0; i<n_tension_models; i++) {
905 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
906 for (j=0; j < tension_models[i].num_params ; j++)
907 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
908 printf("\t\tdefaults: %s\n", tension_models[i].params);
910 printf("Unfolding rate models:\n");
911 for (i=0; i<n_k_models; i++) {
912 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
913 for (j=0; j < k_models[i].num_params ; j++)
914 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
915 printf("\t\tdefaults: %s\n", k_models[i].params);
918 printf("Examples:\n");
919 printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded. Stop once there are no more folded domains.\n");
920 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);
921 printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded. Stop after 0.25 seconds\n");
922 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);
926 \subsection{Get options}
927 \label{sec.get_options}
931 void get_options(int argc, char **argv,
932 double *pP, double *pT_max, double *pDt_max, double *pV,
933 double *pX_step, state_t **pStop_state, environment_t *env,
934 int n_tension_models, tension_model_getopt_t *tension_models,
935 int n_k_models, k_model_getopt_t *k_models,
936 list_t **pState_list, list_t **pTransition_list)
938 char *prog_name = NULL;
939 char c, options[] = "q:t:d:P:v:x:T:C:s:N:k:DVh";
940 int tension_model=0, k_model=0, initial_state=-1, final_state=-1;
941 char *state_name=NULL;
942 int i, n, tension_group=0;
943 list_t *temp_list=NULL;
946 void *model_params; /* for INIT_MODEL() */
947 int num_param_args; /* for INIT_MODEL() */
948 char **param_args; /* for INIT_MODEL() */
950 extern int optind, optopt, opterr;
955 for (i=0; i<n_tension_models; i++) {
956 fprintf(stderr, "tension model %i %s:\t model's handler %p\t inverse handler %p\n",
957 i, tension_models[i].name, tension_models[i].handler, tension_models[i].inverse_handler);
958 assert(tension_models[i].handler == tension_handlers[i]);
963 flags = FLAG_OUTPUT_TRANSITION_FORCES;
966 *pT_max = -1; /* s, -1 removes this limit */
967 *pDt_max = 0.001; /* s */
968 *pP = 1e-3; /* % pop per s */
969 *pV = 1e-6; /* m/s */
970 *pX_step = 0.0; /* m */
971 env->T = 300.0; /* K */
973 *pTransition_list = NULL;
975 while ((c=getopt(argc, argv, options)) != -1) {
978 if (STRMATCH(optarg, "-")) {
981 temp_list = head(*pState_list);
982 while (temp_list != NULL) {
983 if (STRMATCH(S(temp_list)->name, optarg)) {
984 *pStop_state = S(temp_list);
987 temp_list = temp_list->next;
991 case 't': *pT_max = safe_strtod(optarg, "-t"); break;
992 case 'd': *pDt_max = safe_strtod(optarg, "-d"); break;
993 case 'P': *pP = safe_strtod(optarg, "-P"); break;
994 case 'v': *pV = safe_strtod(optarg, "-v"); break;
995 case 'x': *pX_step = safe_strtod(optarg, "-x"); break;
996 case 'T': env->T = safe_strtod(optarg, "-T"); break;
997 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
999 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
1000 assert(num_strings >= 2 && num_strings <= 4);
1002 state_name = strings[0];
1003 if (state_index(*pState_list, state_name) != -1) {
1004 fprintf(stderr, "Invalid option '%s', repeated state name '%s'\n", optarg, state_name);
1007 tension_model = index_tension_model(n_tension_models, tension_models, strings[1]);
1008 if (num_strings == 4)
1009 tension_group = safe_strtoi(strings[2], optarg);
1012 if (num_strings >= 3)
1013 tension_models[tension_model].params = strings[num_strings-1];
1015 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);
1017 INIT_MODEL(state_name, &tension_models[tension_model], tension_models[tension_model].params, model_params);
1018 push(pState_list, create_state(state_name,
1019 tension_models[tension_model].handler,
1020 tension_models[tension_model].inverse_handler,
1023 tension_models[tension_model].destructor,
1025 *pState_list = tail(*pState_list);
1027 free_parsed_list(num_strings, strings);
1030 n = safe_strtoi(optarg, "-N");
1032 fprintf(stderr, "%s:\tadding %d domains to %s\n", __FUNCTION__, n, S(*pState_list)->name);
1034 S(*pState_list)->num_domains += n;
1037 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
1038 assert(num_strings >= 3 && num_strings <= 4);
1039 initial_state = state_index(*pState_list, strings[0]);
1040 final_state = state_index(*pState_list, strings[1]);
1041 k_model = index_k_model(n_k_models, k_models, strings[2]);
1043 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);
1045 assert(initial_state != final_state);
1046 if (num_strings == 4)
1047 k_models[k_model].params = strings[3];
1048 INIT_MODEL(state_name, &k_models[k_model], k_models[k_model].params, model_params);
1049 push(pTransition_list, create_transition(list_element(*pState_list, initial_state),
1050 list_element(*pState_list, final_state),
1051 k_models[k_model].k,
1053 k_models[k_model].destructor), 1);
1054 free_parsed_list(num_strings, strings);
1057 print_state_graph(stdout, *pState_list, *pTransition_list);
1061 flags = FLAG_OUTPUT_FULL_CURVE;
1064 fprintf(stderr, "unrecognized option '%c'\n", optopt);
1065 /* fall through to default case */
1067 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);
1071 /* check the arguments */
1072 assert(*pState_list != NULL); /* no domains! */
1073 assert(*pP > 0.0); assert(*pP < 1.0);
1074 assert(*pDt_max > 0.0);
1076 assert(env->T > 0.0);
1078 crosslink_groups(*pState_list, *pTransition_list);
1084 <<index tension model>>=
1085 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
1088 for (i=0; i<n_models; i++)
1089 if (STRMATCH(models[i].name, name))
1091 if (i == n_models) {
1092 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
1100 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
1103 for (i=0; i<n_models; i++)
1104 if (STRMATCH(models[i].name, name))
1106 if (i == n_models) {
1107 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
1114 <<initialize model definition>>=
1115 /* requires int num_param_args and char **param_args in the current scope
1117 * INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
1118 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
1120 #define INIT_MODEL(role, model, param_string, param_pointer) \
1122 parse_list_string((param_string), SEP, '{', '}', \
1123 &num_param_args, ¶m_args); \
1124 if (num_param_args != (model)->num_params) { \
1126 "%s model %s expected %d params,\n", \
1127 (role), (model)->name, (model)->num_params); \
1129 "not the %d params in '%s'\n", \
1130 num_param_args, (param_string)); \
1131 assert(num_param_args == (model)->num_params); \
1133 if ((model)->creator) \
1134 param_pointer = (*(model)->creator)(param_args); \
1136 param_pointer = NULL; \
1137 free_parsed_list(num_param_args, param_args); \
1141 First we define some safe string parsing functions. Currently these
1142 abort on any irregularity, but in the future they could print error
1143 messages, etc. [[static]] because the functions are currently
1144 defined in each [[*.c]] file that uses them, but perhaps they should
1145 be moved to [[utils.h/c]] or some such instead\ldots
1148 static int safe_strtoi(const char *s, const char *description) {
1152 ret = strtol(s, &endp, 10);
1153 if (endp[0] != '\0') { //strlen(endp) > 0) {
1154 fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1155 endp, s, description);
1156 assert(1==0); //strlen(endp) == 0);
1158 fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1161 } else if (errno == ERANGE) {
1162 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1163 assert(errno != ERANGE);
1168 static double safe_strtod(const char *s, const char *description) {
1172 ret = strtod(s, &endp);
1173 if (endp[0] != '\0') { //strlen(endp) > 0) {
1174 fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1175 endp, s, description);
1176 assert(1==0); //strlen(endp) == 0);
1178 fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1181 } else if (errno == ERANGE) {
1182 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1183 assert(errno != ERANGE);
1192 \addcontentsline{toc}{section}{Appendicies}
1194 \section{sawsim.c details}
1198 The general layout of our simulation code is:
1209 We include [[math.h]], so don't forget to link to the libm with `-lm'.
1211 #include <assert.h> /* assert() */
1212 #include <stdlib.h> /* malloc(), free(), rand(), strto*() */
1213 #include <stdio.h> /* fprintf(), stdout */
1214 #include <string.h> /* strlen, strtok() */
1215 #include <math.h> /* exp(), M_PI, sqrt() */
1216 #include <sys/types.h> /* pid_t (returned by getpid()) */
1217 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
1218 #include <errno.h> /* errno, ERANGE (for safe_strto*()) */
1221 #include "tension_balance.h"
1222 #include "k_model.h"
1223 #include "tension_model.h"
1225 #include "accel_k.h"
1229 <<version definition>>
1230 <<flag definitions>>
1231 <<max/min definitions>>
1232 <<string comparison definition and globals>>
1233 <<initialize model definition>>
1234 <<get options definitions>>
1235 <<state definitions>>
1236 <<transition definitions>>
1245 <<create/destroy state>>
1246 <<destroy state list>>
1248 <<create/destroy transition>>
1249 <<destroy transition list>>
1250 <<print state graph>>
1252 <<simulation functions>>
1253 <<boilerplate functions>>
1256 <<boilerplate functions>>=
1258 <<get option functions>>
1264 srand(getpid()*time(NULL)); /* seed rand() */
1268 <<flag definitions>>=
1269 /* in octal b/c of prefixed '0' */
1270 #define FLAG_OUTPUT_FULL_CURVE 01
1271 #define FLAG_OUTPUT_TRANSITION_FORCES 02
1275 static unsigned long int flags = 0;
1278 \subsection{Utilities}
1281 <<max/min definitions>>=
1282 #define MAX(a,b) ((a)>(b) ? (a) : (b))
1283 #define MIN(a,b) ((a)<(b) ? (a) : (b))
1286 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
1287 <<string comparison definition and globals>>=
1288 // Check if two strings match, return 1 if they do
1289 static char *temp_string_A;
1290 static char *temp_string_B;
1291 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
1292 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
1293 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
1294 /* +1 to also compare the '\0' */
1297 We also define a macro for our [[check]] unit testing
1298 <<check relative error>>=
1299 #define CHECK_ERR(max_err, expected, received) \
1301 fail_unless( (received-expected)/expected < max_err, \
1302 "relative error %g >= %g in %s (Expected %g, received %g)", \
1303 (received-expected)/expected, max_err, #received, \
1304 expected, received); \
1305 fail_unless(-(received-expected)/expected < max_err, \
1306 "relative error %g >= %g in %s (Expected %g, received %g)", \
1307 -(received-expected)/expected, max_err, #received, \
1308 expected, received); \
1313 int happens(double probability)
1315 assert(probability >= 0.0); assert(probability <= 1.0);
1316 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*/
1321 \subsection{Adaptive timesteps}
1322 \label{sec.adaptive_dt_implementation}
1324 $F(x)$ increases with $x$, possibly exploding (e.g.\ in the worm-like
1325 chain model), so basing the timestep on the the unfolding probability
1326 at the current tension is dangerous, and we need to search for a $dt$
1327 for which $P_N(F(x+v*dt)) < P_\text{target}$. There are two cases to
1328 consider. In the most common, no domains have unfolded since the last
1329 step, and we expect the next step to be slightly shorter than the
1330 previous one. In the less common, domains did unfold in the last
1331 step, and we expect the next step to be considerably longer than the
1334 double search_dt(transition_t *transition,
1336 environment_t *env, double x,
1337 double target_prob, double max_dt, double v)
1338 { /* Returns a valid timestep dt in seconds for the given transition.
1339 * Takes a list of all states, a pointer env to the environmental data,
1340 * a starting separation x in m, a target_prob between 0 and 1,
1341 * max_dt in s, stretching velocity v in m/s.
1343 double F, k, P, dtCur, dtU, dtUCur, dtL, dt;
1345 /* get upper bound using the current position */
1346 F = find_tension(state_list, env, x, 0); /* BUG. repeated calculation */
1347 //printf("Start with x = %g (F = %g)\n", x, F);
1348 k = accel_k(transition->k, F, env, transition->k_params);
1349 //printf("x %g\tF %g\tk %g\n", x, F, k);
1350 dtU = P1(target_prob, N_INIT(transition)) / k; /* upper bound on dt */
1352 //printf("overshot max_dt\n");
1355 if (v == 0) /* discrete stepping, so force is temporarily constant */
1358 /* set a lower bound on dt too */
1361 /* The dt determined above may produce illegitimate forces or ks.
1362 * Reduce the upper bound until we have valid ks. */
1364 F = find_tension(state_list, env, x+v*dt, 0);
1365 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1368 F = find_tension(state_list, env, x+v*dt, 0);
1370 //printf("Try for dt = %g (F = %g)\n", dt, F);
1371 k = accel_k(transition->k, F, env, transition->k_params);
1372 /* returned k may be -1.0 */
1373 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1374 while (k == -1.0) { /* reduce step until we hit a valid k */
1376 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1377 F = find_tension(state_list, env, x+v*dt, 0);
1378 //printf("Try for dt = %g (F = %g)\n", dt, F);
1379 k = accel_k(transition->k, F, env, transition->k_params);
1380 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1382 assert(dtU > 1e-14); /* timestep to valid k too small */
1383 /* safe timestep back from x+dtU */
1384 dtUCur = P1(target_prob, N_INIT(transition)) / k;
1386 return dt; /* dtU is safe. */
1388 /* dtU wasn't safe, lets see what would be. */
1389 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1390 dt = (dtU + dtL) / 2.0;
1391 F = find_tension(state_list, env, x+v*dt, 0);
1392 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1393 k = accel_k(transition->k, F, env, transition->k_params);
1394 dtCur = P1(target_prob, N_INIT(transition)) / k;
1395 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1396 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1398 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1399 dtU = dt; /* ... stepping out only dtCur would be safe */
1402 break; /* dtCur = dt */
1404 return MAX(dtUCur, dtL);
1409 To determine $dt$ for an array of potentially different folded domains,
1410 we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
1413 double determine_dt(list_t *state_list, list_t *transition_list,
1414 environment_t *env, double x,
1415 double target_prob, double dt_max, double v)
1416 { /* Returns the timestep dt in seconds.
1417 * Takes the state and transition lists, a pointer to the environment,
1418 * an extension x in m, target_prob between 0 and 1, dt_max in s,
1420 double dt=dt_max, new_dt;
1421 assert(target_prob > 0.0); assert(target_prob < 1.0);
1422 assert(dt_max > 0.0);
1423 assert(state_list != NULL);
1424 assert(transition_list != NULL);
1425 transition_list = head(transition_list);
1426 while (transition_list != NULL) {
1427 if (T(transition_list)->initial_state->num_domains > 0) {
1428 new_dt = search_dt(T(transition_list), state_list, env, x, target_prob, dt, v);
1429 dt = MIN(dt, new_dt);
1431 transition_list = transition_list->next;
1438 \subsection{State data}
1439 \label{sec.state_data}
1441 The domains exist in one of an array of possible states. Each state
1442 has a tension model which determines it's force/extension curve
1443 (possibly [[null]] for rigid states, see Appendix
1444 \ref{sec.null_tension}).
1446 Groups are, for example, for two domain states with WLC tensions; one
1447 with $p=0.4\U{nm}$, $L=10\U{nm}$ and the other with $p=0.4\U{nm}$,
1448 $L=28\U{nm}$. Since the persistence lengths are the same, you would
1449 like to add the contour lengths of all the domains in \emph{both}
1450 states, and plug that total contour length into a single WLC
1451 evaluation (see Section \ref{sec.find_tension}).
1452 <<state definitions>>=
1453 typedef struct state_struct {
1454 char *name; /* name string */
1455 one_dim_fn_t *tension_handler; /* tension handler */
1456 one_dim_fn_t *inverse_tension_handler; /* inverse tension handler */
1457 int tension_group; /* for grouping similar tension states */
1458 void *tension_params; /* pointer to tension parameters */
1459 destroy_data_func_t *destroy_tension;
1460 int num_domains; /* number of domains currently in state */
1461 /* cached lookups generated from the above data */
1462 int num_out_trans; /* length of out_trans array */
1463 struct transition_struct **out_trans; /* transitions leaving this state */
1464 int num_grouped_states; /* length of grouped states array */
1465 struct state_struct **grouped_states; /* states grouped with this one (doesn't include self) */
1468 /* get the state data for the current list node */
1469 #define S(list) ((state_t *)(list)->d)
1471 @ [[tension_params]] is a pointer to the parameters used by the
1472 handler function when determining the tension. [[destroy_tension]]
1473 points to a function for freeing [[tension_params]]. We it in the
1474 state structure so that [[destroy_state]] and will have an easier job
1477 [[create_]] and [[destroy_state]] are simple wrappers around
1478 [[malloc]] and [[free]].
1479 <<create/destroy state>>=
1480 state_t *create_state(char *name,
1481 one_dim_fn_t *tension_handler,
1482 one_dim_fn_t *inverse_tension_handler,
1484 void *tension_params,
1485 destroy_data_func_t *destroy_tension,
1488 state_t *ret = (state_t *)malloc(sizeof(state_t));
1489 assert(ret != NULL);
1490 /* make a copy of the name, so we aren't dependent on the original */
1491 ret->name = (char *)malloc(sizeof(char)*(strlen(name)+1));
1492 assert(ret->name != NULL);
1493 strcpy(ret->name, name); /* we know ret->name is long enough */
1495 ret->tension_handler = tension_handler;
1496 ret->inverse_tension_handler = inverse_tension_handler;
1497 ret->tension_group = tension_group;
1498 ret->tension_params = tension_params;
1499 ret->destroy_tension = destroy_tension;
1500 ret->num_domains = num_domains;
1502 ret->num_out_trans = 0;
1503 ret->out_trans = NULL;
1504 ret->num_grouped_states = 0;
1505 ret->grouped_states = NULL;
1508 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);
1513 void destroy_state(state_t *state)
1517 fprintf(stderr, "freeing state %s (%p)\n", state->name, state);
1521 if (state->destroy_tension)
1522 (*state->destroy_tension)(state->tension_params);
1523 if (state->out_trans)
1524 free(state->out_trans);
1525 if (state->grouped_states)
1526 free(state->grouped_states);
1533 We'll be storing the states in a list (see Appendix \ref{sec.lists}),
1534 so we also define a convenience function for cleaning up.
1535 <<destroy state list>>=
1536 void destroy_state_list(list_t *state_list)
1538 state_list = head(state_list);
1539 while (state_list != NULL)
1540 destroy_state((state_t *) pop(&state_list));
1545 We can also define a convenient way to get a state index from it's name.
1547 int state_index(list_t *state_list, char *name)
1550 state_list = head(state_list);
1551 while (state_list != NULL) {
1552 if (STRMATCH(S(state_list)->name, name)) return i;
1553 state_list = state_list->next;
1561 \subsection{Transition data}
1562 \label{sec.transition_data}
1564 Once you've created states (Sections \ref{sec.main},
1565 \ref{sec.model_selection}, and \ref{sec.state_data}), you need to
1566 create transitions between them.
1567 <<transition definitions>>=
1568 typedef struct transition_struct {
1569 state_t *initial_state;
1570 state_t *final_state;
1571 k_func_t *k; /* transition rate model */
1572 void *k_params; /* pointer to k parameters */
1573 destroy_data_func_t *destroy_k;
1576 /* get the transition data for the current list node */
1577 #define T(list) ((transition_t *)(list)->d)
1579 /* get the number of initial state population for a transition state */
1580 #define N_INIT(transition) ((transition)->initial_state->num_domains)
1584 @ [[k]] is a pointer to the function determining the transition rate
1585 for a given tension. [[k_params]] and [[destroy_k]] have similar
1586 roles with regards to [[k]] as [[tension_params]] and
1587 [[destroy_tension]] have with regards to [[tension_handler]] (see
1588 Section \ref{sec.state_data}).
1590 [[create_]] and [[destroy_transition]] are simple wrappers around
1591 [[malloc]] and [[free]].
1592 <<create/destroy transition>>=
1593 transition_t *create_transition(state_t *initial_state, state_t *final_state,
1596 destroy_data_func_t *destroy_k)
1598 transition_t *ret = (transition_t *)malloc(sizeof(transition_t));
1599 assert(ret != NULL);
1600 assert(initial_state != NULL);
1601 assert(final_state != NULL);
1602 ret->initial_state = initial_state;
1603 ret->final_state = final_state;
1605 ret->k_params = k_params;
1606 ret->destroy_k = destroy_k;
1608 fprintf(stderr, "generated transition %p from %s to %s, params %p\n", ret, ret->initial_state->name, ret->final_state->name, ret->k_params);
1613 void destroy_transition(transition_t *transition)
1617 fprintf(stderr, "freeing transition %p (%s to %s)\n", transition, transition->initial_state->name, transition->final_state->name);
1619 if (transition->destroy_k)
1620 (*transition->destroy_k)(transition->k_params);
1627 We'll be storing the transitions in a list (see Appendix
1628 \ref{sec.lists}), so we also define a convenience function for
1630 <<destroy transition list>>=
1631 void destroy_transition_list(list_t *transition_list)
1633 transition_list = head(transition_list);
1634 while (transition_list != NULL)
1635 destroy_transition((transition_t *) pop(&transition_list));
1640 \subsection{Graphviz output}
1641 \label{sec.graphviz_output}
1643 <<print state graph>>=
1644 void print_state_graph(FILE *file, list_t *state_list, list_t *transition_list)
1646 state_list = head(state_list);
1647 transition_list = head(transition_list);
1648 fprintf(file, "digraph states {\n node [shape=record,color=black,fillcolor=white,style=filled];\n\n");
1650 fprintf(file, " n%d [label=\"%s\"]\n", state_index(state_list, S(state_list)->name), S(state_list)->name);
1651 if (state_list->next == NULL) break;
1652 state_list = state_list->next;
1654 fprintf(file, "\n");
1655 while (transition_list != NULL) {
1656 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));
1657 transition_list = transition_list->next;
1659 fprintf(file, "}\n");
1663 \subsection{Domain model and group handling}
1665 <<group functions>>=
1666 <<crosslink groups>>
1667 <<get active groups>>
1670 Scan through a list of states and construct an array of all the
1672 <<get active groups>>=
1673 void get_active_groups(list_t *state_list,
1674 int *pNum_active_groups,
1675 one_dim_fn_t ***pPer_group_handlers,
1676 one_dim_fn_t ***pPer_group_inverse_handlers,
1677 void ***pPer_group_data /* pointer to array of pointers to tension_handler_data_t */)
1679 void **active_groups=NULL;
1680 one_dim_fn_t *handler, *old_handler, *inverse_handler, *old_inverse_handler;
1681 one_dim_fn_t **per_group_handlers=NULL, **per_group_inverse_handlers=NULL;
1682 void **per_group_data=NULL; /* array of pointers to tension_handler_data_t */
1683 int i, j, num_domains, group, old_group, num_active_groups=0;
1684 list_t *active_states_list=NULL;
1685 tension_handler_data_t *tdata=NULL;
1688 /* get all the active groups in a list */
1689 state_list = head(state_list);
1691 fprintf(stderr, "%s:\t", __FUNCTION__);
1692 list_print(stderr, state_list, "states");
1694 while (state_list != NULL) {
1695 state = S(state_list);
1696 handler = state->tension_handler;
1697 inverse_handler = state->inverse_tension_handler;
1698 group = state->tension_group;
1699 num_domains = state->num_domains;
1700 if (list_index(active_states_list, state) == -1) {
1701 /* we haven't added this state (or it's associated group) yet */
1702 if (num_domains > 0 && handler != NULL) { /* add the state */
1703 num_active_groups++;
1704 push(&active_states_list, state, 1);
1706 fprintf(stderr, "%s:\tactive %s state (%p) with %d domain(s)\n", __FUNCTION__, state->name, state, state->num_domains);
1708 for (i=0; i < state->num_grouped_states; i++) {
1709 if (state->grouped_states[i]->num_domains > 0) {
1710 /* add this related (and active) state now */
1711 assert(state->grouped_states[i]->tension_handler == handler);
1712 assert(state->grouped_states[i]->tension_group == group);
1713 push(&active_states_list, state->grouped_states[i], 1);
1715 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);
1719 } /* else empty state or NULL handler */
1720 } /* else state already added */
1721 state_list = state_list->next;
1724 fprintf(stderr, "%s:\t", __FUNCTION__);
1725 list_print(stderr, active_states_list, "active states");
1728 assert(num_active_groups <= list_length(active_states_list));
1729 assert(num_active_groups > 0);
1731 /* allocate the arrays we'll be returning */
1732 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1733 assert(per_group_handlers != NULL);
1734 per_group_inverse_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1735 assert(per_group_inverse_handlers != NULL);
1736 per_group_data = (void **)malloc(sizeof(void *)*num_active_groups);
1737 assert(per_group_data != NULL);
1739 fprintf(stderr, "%s:\tallocating per group data %p:", __FUNCTION__, per_group_data);
1741 for (i=0; i<num_active_groups; i++) { /* we might want to recycle any group data structs allocated in the previous round */
1742 per_group_data[i] = malloc(sizeof(tension_handler_data_t));
1743 assert(per_group_data[i] != NULL);
1745 fprintf(stderr, " %d,%p", i, per_group_data[i]);
1749 fprintf(stderr, "\n");
1752 /* populate the arrays we'll be returning */
1753 active_states_list = head(active_states_list);
1755 old_inverse_handler = NULL;
1756 j = -1; /* j is the current active _group_ index */
1757 while (active_states_list != NULL) {
1758 state = (state_t *) pop(&active_states_list);
1759 handler = state->tension_handler;
1760 inverse_handler = state->inverse_tension_handler;
1761 group = state->tension_group;
1762 num_domains = state->num_domains;
1763 if (handler != old_handler || inverse_handler != old_inverse_handler || group != old_group) {
1764 j++; /* move to the next group */
1765 tdata = (tension_handler_data_t *) per_group_data[j];
1766 per_group_handlers[j] = handler;
1767 per_group_inverse_handlers[j] = inverse_handler;
1768 tdata->group_tension_params = NULL;
1770 tdata->persist = NULL; /* we might want to grab any appropriate data from the old persist here... */
1773 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);
1775 for (i=0; i<num_domains; i++) { /* add num_domains copies of the tension params */
1776 push(&tdata->group_tension_params, state->tension_params, 1);
1779 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);
1781 old_handler = handler;
1782 old_inverse_handler = inverse_handler;
1786 /* free old groups */
1787 if (*pPer_group_handlers != NULL)
1788 free(*pPer_group_handlers);
1789 if (*pPer_group_inverse_handlers != NULL)
1790 free(*pPer_group_inverse_handlers);
1791 if (*pPer_group_data != NULL) {
1792 for (i=0; i<*pNum_active_groups; i++)
1793 free((*pPer_group_data)[i]);
1794 free(*pPer_group_data);
1797 /* set new groups */
1799 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]);
1801 *pNum_active_groups = num_active_groups;
1802 *pPer_group_handlers = per_group_handlers;
1803 *pPer_group_inverse_handlers = per_group_inverse_handlers;
1804 *pPer_group_data = per_group_data;
1808 <<crosslink groups>>=
1809 void crosslink_groups(list_t *state_groups, list_t *transition_list)
1811 list_t *list, *out_trans=NULL, *associates=NULL;
1812 one_dim_fn_t *handler;
1815 state_groups = head(state_groups);
1816 while (state_groups != NULL) {
1817 /* find transitions leaving this state */
1819 fprintf(stderr, "finding transitions leaving %s\n", S(state_groups)->name);
1821 list = head(transition_list);
1822 while (list != NULL) {
1823 if (T(list)->initial_state == S(state_groups) && list_index(associates, T(list)->initial_state) != -1) {
1824 push(&out_trans, T(list), 1);
1828 n = list_length(out_trans);
1829 S(state_groups)->num_out_trans = n;
1830 S(state_groups)->out_trans = (transition_t **)malloc(sizeof(transition_t *)*n);
1831 assert(S(state_groups)->out_trans != NULL);
1832 for (i=0; i<n; i++) {
1833 S(state_groups)->out_trans[i] = (transition_t *)pop(&out_trans);
1836 /* find states grouped with this state */
1838 fprintf(stderr, "finding states grouped with %s\n", S(state_groups)->name);
1840 handler = S(state_groups)->tension_handler;
1841 group = S(state_groups)->tension_group;
1842 list = head(state_groups);
1843 while (list != NULL) {
1844 if (S(list)->tension_handler == handler && S(list)->tension_group == group && list != state_groups) {
1845 push(&associates, S(list), 1);
1849 n = list_length(associates);
1850 S(state_groups)->num_grouped_states = n;
1851 S(state_groups)->grouped_states = (state_t **)malloc(sizeof(state_t *)*n);
1852 assert(S(state_groups)->grouped_states != NULL);
1853 for (i=0; i<n; i++) {
1854 S(state_groups)->grouped_states[i] = (state_t *)pop(&associates);
1856 state_groups = state_groups->next;
1862 \section{String parsing}
1864 For handling command line arguments, we need parse delimited strings (e.g.\ [["param1,param2,param3..."]]).
1865 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1866 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1867 need to take care of parsing those parameters themselves.
1868 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1874 <<parse definitions>>
1875 <<parse declarations>>
1879 <<parse module makefile lines>>=
1880 NW_SAWSIM_MODS += parse
1881 CHECK_BINS += check_parse
1885 <<parse definitions>>=
1886 #define SEP ',' /* argument separator character */
1889 <<parse declarations>>=
1890 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1891 int *num, char ***string_array);
1892 extern void free_parsed_list(int num, char **string_array);
1895 [[parse_list_string]] allocates memory, don't forget to free it
1896 afterward with [[free_parsed_list]]. It does not alter the original.
1898 The string may start off with a [[deeper]] character
1899 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1900 [[x,y]], where the pointer is one character in on the copied string.
1901 However, when we go to free the memory, we need a pointer to the
1902 beginning of the string. In order to accommodate this for a string
1903 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1904 the first $N$ elements point to the separated arguments, and let the
1905 last element point to the start of the copied string regardless of
1907 <<parse delimited list string functions>>=
1908 /* TODO, split out into parse.hc */
1909 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1912 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1913 if (string[i] == deeper) {depth++;}
1914 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1920 void parse_list_string(char *string, char sep, char deeper, char shallower,
1921 int *num, char ***string_array)
1923 char *str=NULL, **ret=NULL;
1925 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
1927 *string_array = NULL;
1930 /* make a copy of the string, so we don't change the original */
1931 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1932 assert(str != NULL);
1933 strcpy(str, string); /* we know str is long enough */
1934 /* count the number of regions, so we can allocate pointers to them */
1937 n++; i++; /* move on to next argument */
1938 i += next_delim_index(str+i, sep, deeper, shallower);
1939 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1940 } while (str[i] != '\0');
1941 ret = (char **)malloc(sizeof(char *)*(n+1));
1942 assert(ret != NULL);
1943 /* replace the separators with '\0' & assign pointers */
1944 ret[n] = str; /* point to the front of the copied string */
1947 for(i=1; i<n; i++) {
1948 j += next_delim_index(str+j, sep, deeper, shallower);
1950 ret[i] = str+j; /* point to the separated arguments */
1952 /* strip off bounding braces, if any */
1953 for(i=0; i<n; i++) {
1954 if (ret[i][0] == deeper) {
1958 if (ret[i][strlen(ret[i])-1] == shallower)
1959 ret[i][strlen(ret[i])-1] = '\0';
1962 *string_array = ret;
1965 void free_parsed_list(int num, char **string_array)
1968 /* frees all items, since they were allocated as a single string*/
1969 free(string_array[num]);
1970 /* and free the array of pointers */
1976 \subsection{Parsing implementation}
1982 <<parse delimited list string functions>>
1986 #include <assert.h> /* assert() */
1987 #include <stdlib.h> /* NULL */
1988 #include <stdio.h> /* fprintf(), stdout *//*!!*/
1989 #include <string.h> /* strlen() */
1993 \subsection{Parsing unit tests}
1995 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1997 <<parse check includes>>
1998 <<string comparison definition and globals>>
1999 <<check relative error>>
2000 <<parse test suite>>
2001 <<main check program>>
2004 <<parse check includes>>=
2005 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2006 #include <stdio.h> /* printf() */
2007 #include <assert.h> /* assert() */
2008 #include <string.h> /* strlen() */
2013 <<parse test suite>>=
2014 <<parse list string tests>>
2015 <<parse suite definition>>
2018 <<parse suite definition>>=
2019 Suite *test_suite (void)
2021 Suite *s = suite_create ("k model");
2022 <<parse list string test case defs>>
2024 <<parse list string test case add>>
2029 <<parse list string tests>>=
2032 START_TEST(test_next_delim_index)
2034 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
2035 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
2036 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
2037 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
2038 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
2042 START_TEST(test_parse_list_null)
2046 parse_list_string(NULL, SEP, '{', '}',
2047 &num_param_args, ¶m_args);
2048 fail_unless(num_param_args == 0, NULL);
2049 fail_unless(param_args == NULL, NULL);
2052 START_TEST(test_parse_list_single_simple)
2056 parse_list_string("arg", SEP, '{', '}',
2057 &num_param_args, ¶m_args);
2058 fail_unless(num_param_args == 1, NULL);
2059 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
2062 START_TEST(test_parse_list_single_compound)
2066 parse_list_string("{x,y,z}", SEP, '{', '}',
2067 &num_param_args, ¶m_args);
2068 fail_unless(num_param_args == 1, NULL);
2069 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
2072 START_TEST(test_parse_list_double_simple)
2076 parse_list_string("abc,def", SEP, '{', '}',
2077 &num_param_args, ¶m_args);
2078 fail_unless(num_param_args == 2, NULL);
2079 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2080 fail_unless(STRMATCH(param_args[1],"def"), NULL);
2083 START_TEST(test_parse_list_compound)
2087 parse_list_string("abc,{def,ghi}", SEP, '{', '}',
2088 &num_param_args, ¶m_args);
2089 fail_unless(num_param_args == 2, NULL);
2090 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2091 fail_unless(STRMATCH(param_args[1],"def,ghi"), NULL);
2095 <<parse list string test case defs>>=
2096 TCase *tc_parse_list_string = tcase_create("parse list string");
2098 <<parse list string test case add>>=
2099 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
2100 tcase_add_test(tc_parse_list_string, test_parse_list_null);
2101 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
2102 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
2103 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
2104 tcase_add_test(tc_parse_list_string, test_parse_list_compound);
2105 suite_add_tcase(s, tc_parse_list_string);
2109 \section{Unit tests}
2111 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2118 <<check relative error>>
2121 <<main check program>>
2133 <<determine dt tests>>
2135 <<does domain transition tests>>
2136 <<randomly unfold domains tests>>
2137 <<suite definition>>
2140 <<suite definition>>=
2141 Suite *test_suite (void)
2143 Suite *s = suite_create ("sawsim");
2144 <<F test case defs>>
2145 <<determine dt test case defs>>
2146 <<happens test case defs>>
2147 <<does domain transition test case defs>>
2148 <<randomly unfold domains test case defs>>
2151 <<determine dt test case add>>
2152 <<happens test case add>>
2153 <<does domain transition test case add>>
2154 <<randomly unfold domains test case add>>
2157 tcase_add_checked_fixture(tc_strip_address,
2158 setup_strip_address,
2159 teardown_strip_address);
2165 <<main check program>>=
2169 Suite *s = test_suite();
2170 SRunner *sr = srunner_create(s);
2171 srunner_run_all(sr, CK_ENV);
2172 number_failed = srunner_ntests_failed(sr);
2174 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
2178 \subsection{F tests}
2180 <<worm-like chain tests>>
2182 <<F test case defs>>=
2183 <<worm-like chain test case def>>
2185 <<F test case add>>=
2186 <<worm-like chain test case add>>
2189 <<worm-like chain tests>>=
2190 <<worm-like chain function>>
2192 START_TEST(test_wlc_at_zero)
2194 double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
2195 fail_unless(abs(wlc(x, T, p, L)) < lim, \
2196 "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
2197 x, T, p, L, abs(wlc(x, T, p, L)), lim);
2201 START_TEST(test_wlc_at_half)
2203 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
2204 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
2205 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
2207 * linear term = x/L = 0.5
2208 * nonlinear + linear = 0.75 + 0.5 = 1.25
2209 * wlc = 10e21*1.25 = 12.5e21
2211 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
2212 "wlc(%g, %g, %g, %g) = %g != %g",
2213 x, T, p, L, wlc(x, T, p, L), 12.5e21);
2217 <<worm-like chain test case def>>=
2218 TCase *tc_wlc = tcase_create("WLC");
2221 <<worm-like chain test case add>>=
2222 tcase_add_test(tc_wlc, test_wlc_at_zero);
2223 tcase_add_test(tc_wlc, test_wlc_at_half);
2224 suite_add_tcase(s, tc_wlc);
2227 \subsection{Model tests}
2229 Check the searching with [[linear_k]].
2230 Check overwhelming force treatment with the heavyside-esque step [[?]].
2232 Hmm, I'm not really sure what I was doing here...
2233 <<determine dt tests>>=
2234 double linear_k(double F, environment_t *env, void *params)
2236 double Fnot = *(double *)params;
2240 START_TEST(test_determine_dt_linear_k)
2243 transition_t unfold={0};
2244 environment_t env = {0};
2245 double F[]={0,1,2,3,4,5,6};
2246 double dt_max=3.0, Fnot=3.0, x=1.0, targ_p=0.1, v=1.0;
2247 list_t *state_list=NULL, *transition_list=NULL;
2250 folded.tension_handler = &hooke_handler;
2251 folded.num_domains = 1;
2252 unfold.initial_state = &folded;
2253 unfold.k = linear_k;
2254 unfold.k_params = &Fnot;
2255 push(&state_list, &folded, 1);
2256 push(&transition_list, &unfold, 1);
2257 for( i=0; i < sizeof(F)/sizeof(double); i++) {
2258 //fail_unless(determine_dt(state_list, transition_list, &env, x, targ_p, dt_max, v)==1, NULL);
2263 <<determine dt test case defs>>=
2264 TCase *tc_determine_dt = tcase_create("Determine dt");
2266 <<determine dt test case add>>=
2267 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
2268 suite_add_tcase(s, tc_determine_dt);
2273 <<happens test case defs>>=
2275 <<happens test case add>>=
2278 <<does domain transition tests>>=
2280 <<does domain transition test case defs>>=
2282 <<does domain transition test case add>>=
2285 <<randomly unfold domains tests>>=
2287 <<randomly unfold domains test case defs>>=
2289 <<randomly unfold domains test case add>>=
2293 \section{Balancing group extensions}
2294 \label{sec.tension_balance}
2296 For a given total extension $x$ (of the piezo), the various domain
2297 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
2298 amounts, and we need to tweak the portion that each extends to balance
2299 the tension amongst the active groups. Since the tension balancing is
2300 separable from the bulk of the simulation, we break it out into a
2301 separate module. The interface is defined in [[tension_balance.h]],
2302 the implementation in [[tension_balance.c]], and the unit testing in
2303 [[check_tension_balance.c]]
2305 <<tension-balance.h>>=
2306 #ifndef TENSION_BALANCE_H
2307 #define TENSION_BALANCE_H
2309 <<tension balance function declaration>>
2310 #endif /* TENSION_BALANCE_H */
2313 <<tension balance functions>>=
2314 <<one dimensional bracket>>
2315 <<one dimensional solve>>
2317 <<group stiffness function>>
2318 <<chain stiffness function>>
2319 <<tension balance function>>
2322 <<tension balance module makefile lines>>=
2323 NW_SAWSIM_MODS += tension_balance
2324 CHECK_BINS += check_tension_balance
2325 check_tension_balance_MODS =
2328 The entire force balancing problem reduces to a solving two nested
2329 one-dimensional root-finding problems. First we define the one
2330 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
2331 populated group). $x(x_0)$ is also strictly monotonically increasing,
2332 so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
2333 stores the last successful value of $x$, since we expect to be taking
2334 small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
2335 indicates that the groups have changed.
2336 <<tension balance function declaration>>=
2337 double tension_balance(int num_tension_groups,
2338 one_dim_fn_t **tension_handlers,
2339 one_dim_fn_t **inverse_tension_handlers,
2340 void **params, /* array of pointers to tension_handler_data_t */
2345 <<tension balance function>>=
2346 double tension_balance(int num_tension_groups,
2347 one_dim_fn_t **tension_handlers,
2348 one_dim_fn_t **inverse_tension_handlers,
2349 void **params, /* array of pointers to tension_handler_data_t */
2354 static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL,NULL};
2355 double F, xo_guess, xo, lb, ub;
2356 double min_relx=1e-6, min_rely=1e-6;
2357 int max_steps=200, i;
2359 assert(num_tension_groups > 0);
2360 assert(tension_handlers != NULL);
2361 assert(params != NULL);
2362 assert(num_tension_groups > 0);
2364 if (last_x == -1 || x_xo_data.xi == NULL) { /* new group setup, reset x_xo_data */
2365 /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
2366 if (x_xo_data.xi != NULL)
2368 /* construct new x_xo_data */
2369 x_xo_data.n_groups = num_tension_groups;
2370 x_xo_data.tension_handler = tension_handlers;
2371 x_xo_data.inverse_tension_handler = inverse_tension_handlers;
2372 x_xo_data.handler_data = params;
2374 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);
2376 x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
2377 assert(x_xo_data.xi != NULL);
2378 for (i=0; i<num_tension_groups; i++)
2379 x_xo_data.xi[i] = last_x;
2382 if (num_tension_groups == 1) { /* only one group, no balancing required */
2384 x_xo_data.xi[0] = xo;
2386 fprintf(stderr, "%s:\tonly one active group, no balancing required\n", __FUNCTION__);
2390 fprintf(stderr, "%s:\tbalancing for x = %g with", __FUNCTION__, x);
2391 for(i=0; i<num_tension_groups; i++) fprintf(stderr, "\tx[%d]=%g", i, x_xo_data.xi[i]);
2392 fprintf(stderr, "\n");
2394 if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
2395 if (x == 0) xo_guess = length_scale;
2396 else xo_guess = x/num_tension_groups;
2398 fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
2400 oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
2401 } else { /* work off the last balanced point */
2403 fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
2406 lb = x_xo_data.xi[0];
2407 ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */
2408 } else if (x < last_x) {
2409 lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
2410 ub = x_xo_data.xi[0];
2411 } else { /* x == last_x */
2412 fprintf(stderr, "not moving\n");
2413 lb= x_xo_data.xi[0];
2414 ub= x_xo_data.xi[0];
2418 fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
2419 __FUNCTION__, x, lb, ub);
2421 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
2422 F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
2424 if (num_tension_groups > 1) {
2425 fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
2426 for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
2427 fprintf(stderr, "\n");
2432 F = (*tension_handlers[0])(xo, params[0]);
2434 if (stiffness != NULL)
2435 *stiffness = chain_stiffness(&x_xo_data, min_relx);
2442 <<tension balance internal definitions>>=
2443 <<x of x0 definitions>>
2446 <<x of x0 definitions>>=
2447 typedef struct x_of_xo_data_struct {
2449 one_dim_fn_t **tension_handler; /* array of fn pointers */
2450 one_dim_fn_t **inverse_tension_handler; /* inverse fn pointers */
2451 void **handler_data; /* array of pointers to tension_handler_data_t */
2452 double *xi; /* array of group extensions */
2457 double x_of_xo(double xo, void *pdata)
2459 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
2460 double F, x=0, xi, xi_guess, lb, ub, slope;
2462 double min_relx=1e-6, min_rely=1e-6;
2464 assert(data->n_groups > 0);
2467 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);
2470 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
2472 fprintf(stderr, "%s:\tfinding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
2476 if (data->xi) data->xi[0] = xo;
2477 for (i=1; i < data->n_groups; i++) {
2478 if (data->inverse_tension_handler[i] != NULL) {
2479 xi = (*data->inverse_tension_handler[i])(F, data->handler_data[i]);
2480 } else { /* invert numerically */
2481 if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
2482 else xi_guess = data->xi[i];
2484 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]);
2486 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
2487 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
2492 fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
2496 fprintf(stderr, "%s:\tfound x(x0=%g) = %g\n", __FUNCTION__, xo, x);
2502 <<group stiffness function>>=
2503 double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double x, double relx)
2505 double dx, xl, F, Fl, stiffness;
2507 assert(relx > 0 && relx < 1);
2508 if (x == 0 || relx == 0) {
2509 dx = 1e-12; /* HACK, how to get length scale? */
2519 F = (*tension_handler)(x, handler_data);
2520 Fl = (*tension_handler)(xl, handler_data);
2521 stiffness = (F-Fl)/dx;
2522 assert(stiffness >= 0);
2527 <<chain stiffness function>>=
2528 double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx)
2532 /* add group stiffnesses in series */
2533 for (i=0; i < x_xo_data->n_groups; i++) {
2535 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);
2538 stiffness += 1.0/group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx);
2540 stiffness = 1.0 / stiffness;
2546 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
2547 number of steps for monotonically increasing $f(x)$. Simple
2548 bisection, so it's robust and converges fairly quickly. You can set
2549 the maximum number of steps to take, but if you have enough processor
2550 power, it's better to limit your precision with the relative limits
2551 [[min_relx]] and [[y]]. These ensure that the width of the bracket is
2552 small on the length scale of it's larger side. Note that \emph{both}
2553 relative limits must be satisfied for the search to stop.
2554 <<one dimensional function definition>>=
2555 /* equivalent to gsl_func_t */
2556 typedef double one_dim_fn_t(double x, void *params);
2558 <<one dimensional solve declaration>>=
2559 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2560 double min_relx, double min_rely, int max_steps,
2564 <<one dimensional solve>>=
2565 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2566 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2567 double min_relx, double min_rely, int max_steps,
2570 double yx, ylb, yub, x;
2573 ylb = (*f)(lb, params);
2574 yub = (*f)(ub, params);
2576 /* check some simple cases */
2578 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bracketed */
2579 else return lb; /* any x on the range [lb, ub] would work */
2581 if (ylb == y) { x=lb; yx=ylb; goto end; }
2582 if (yub == y) { x=ub; yx=yub; goto end; }
2585 fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2587 assert(ylb < y); assert(yub > y);
2588 for (i=0; i<max_steps; i++) {
2590 yx = (*f)(x, params);
2592 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);
2594 if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2595 else if (yx > y) { ub=x; yub=yx; }
2596 else /* yx < y */ { lb=x; ylb=yx; }
2601 fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2607 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2608 Generate bracketing $x$ values through bisection or exponential growth.
2609 <<one dimensional bracket declaration>>=
2610 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2613 <<one dimensional bracket>>=
2614 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2616 double yx, step, x=xguess;
2619 fprintf(stderr, "%s:\tbracketing %g (params %p)\n", __FUNCTION__, y, params);
2621 yx = (*f)(x, params);
2623 fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2630 fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2634 if (x == 0) x = length_scale/2.0; /* HACK */
2637 yx = (*f)(x, params);
2639 fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2644 fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2647 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2651 \subsection{Balancing implementation}
2653 <<tension-balance.c>>=
2656 <<tension balance includes>>
2657 #include "tension_balance.h"
2659 double length_scale = 1e-10; /* HACK */
2661 <<tension balance internal definitions>>
2662 <<tension balance functions>>
2665 <<tension balance includes>>=
2666 #include <assert.h> /* assert() */
2667 #include <stdlib.h> /* NULL */
2668 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2669 #include <stdio.h> /* fprintf(), stdout */
2673 \subsection{Balancing unit tests}
2675 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2676 <<check-tension-balance.c>>=
2677 <<tension balance check includes>>
2678 <<tension balance test suite>>
2679 <<main check program>>
2682 <<tension balance check includes>>=
2683 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2684 #include <assert.h> /* assert() */
2687 #include "tension_balance.h"
2690 <<tension balance test suite>>=
2691 <<tension balance function tests>>
2692 <<tension balance suite definition>>
2695 <<tension balance suite definition>>=
2696 Suite *test_suite (void)
2698 Suite *s = suite_create ("tension balance");
2699 <<tension balance function test case defs>>
2701 <<tension balance function test case adds>>
2706 <<tension balance function tests>>=
2707 <<check relative error>>
2709 double hooke(double x, void *pK)
2712 return *((double*)pK) * x;
2715 START_TEST(test_single_function)
2717 double k=5, x=3, last_x=2.0, F;
2718 one_dim_fn_t *handlers[] = {&hooke};
2719 one_dim_fn_t *inverse_handlers[] = {NULL};
2720 void *data[] = {&k};
2721 F = tension_balance(1, handlers, inverse_handlers, data, last_x, x, NULL);
2722 fail_unless(F = k*x, NULL);
2727 We can also test balancing two springs with different spring constants.
2728 The analytic solution for a general number of springs is given in Appendix \ref{sec.math_hooke}.
2729 <<tension balance function tests>>=
2730 START_TEST(test_double_hooke)
2732 double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
2733 one_dim_fn_t *handlers[] = {&hooke, &hooke};
2734 one_dim_fn_t *inverse_handlers[] = {NULL, NULL};
2735 void *data[] = {&k1, &k2};
2736 F = tension_balance(2, handlers, inverse_handlers, data, last_x, x, NULL);
2740 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2741 //CHECK_ERR(1e-6, x1e, xi[0]);
2742 //CHECK_ERR(1e-6, x2e, xi[1]);
2743 CHECK_ERR(1e-6, Fe, F);
2747 TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
2749 <<tension balance function test case defs>>=
2750 TCase *tc_tbfunc = tcase_create("tension balance function");
2753 <<tension balance function test case adds>>=
2754 tcase_add_test(tc_tbfunc, test_single_function);
2755 tcase_add_test(tc_tbfunc, test_double_hooke);
2756 suite_add_tcase(s, tc_tbfunc);
2762 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2763 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2764 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2770 <<list definitions>>
2771 <<list declarations>>
2775 <<list declarations>>=
2776 <<head and tail declarations>>
2777 <<list length declaration>>
2778 <<push declaration>>
2780 <<list destroy declaration>>
2781 <<list to array declaration>>
2782 <<move declaration>>
2783 <<sort declaration>>
2784 <<list index declaration>>
2785 <<list element declaration>>
2786 <<unique declaration>>
2787 <<list print declaration>>
2791 <<create and destroy node>>
2806 <<list module makefile lines>>=
2807 NW_SAWSIM_MODS += list
2808 CHECK_BINS += check_list
2812 <<list definitions>>=
2813 typedef struct list_struct {
2814 struct list_struct *next;
2815 struct list_struct *prev;
2819 <<comparison function typedef>>
2822 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2823 <<head and tail declarations>>=
2824 list_t *head(list_t *list);
2825 list_t *tail(list_t *list);
2828 list_t *head(list_t *list)
2832 while (list->prev != NULL) {
2838 list_t *tail(list_t *list)
2842 while (list->next != NULL) {
2849 <<list length declaration>>=
2850 int list_length(list_t *list);
2853 int list_length(list_t *list)
2860 while (list->next != NULL) {
2868 [[push]] inserts a node relative to the current position in the list
2869 without changing the current position.
2870 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2871 so we set it to that of the pushed domain.
2872 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2873 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2874 <<push declaration>>=
2875 void push(list_t **pList, void *data, int below);
2878 void push(list_t **pList, void *data, int below)
2880 list_t *list, *new_node;
2881 assert(pList != NULL);
2883 new_node = create_node(data);
2886 else if (below==0) { /* insert above */
2887 if (list->prev != NULL)
2888 list->prev->next = new_node;
2889 new_node->prev = list->prev;
2890 list->prev = new_node;
2891 new_node->next = list;
2892 } else { /* insert below */
2893 if (list->next != NULL)
2894 list->next->prev = new_node;
2895 new_node->next = list->next;
2896 list->next = new_node;
2897 new_node->prev = list;
2902 [[pop]] removes the current domain node, moving the current position
2903 to the node after it, unless that node is [[NULL]], in which case move
2904 the current position to the node before it.
2905 <<pop declaration>>=
2906 void *pop(list_t **pList);
2909 void *pop(list_t **pList)
2911 list_t *list, *popped;
2913 assert(pList != NULL);
2915 assert(list != NULL); /* not an empty list */
2917 /* bypass the popped node */
2918 if (list->prev != NULL)
2919 list->prev->next = list->next;
2920 if (list->next != NULL) {
2921 list->next->prev = list->prev;
2922 *pList = list->next;
2924 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
2926 destroy_node(popped);
2931 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
2932 If the cleanup function is [[NULL]], no cleanup function is called.
2933 The nodes are not popped in any particular order.
2934 <<list destroy declaration>>=
2935 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
2938 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
2942 assert(pList != NULL);
2945 assert(list != NULL); /* not an empty list */
2946 while (list != NULL) {
2948 if (destroy != NULL)
2954 [[list_to_array]] converts a list to an array of pointers, preserving order.
2955 <<list to array declaration>>=
2956 void list_to_array(list_t *list, int *length, void ***array);
2959 void list_to_array(list_t *list, int *pLength, void ***pArray)
2963 assert(list != NULL);
2964 assert(pLength != NULL);
2965 assert(pArray != NULL);
2967 length = list_length(list);
2968 array = (void **)malloc(sizeof(void *)*length);
2969 assert(array != NULL);
2972 while (list != NULL) {
2973 array[i++] = list->d;
2981 [[move]] moves the current element along the list in either direction.
2982 <<move declaration>>=
2983 void move(list_t **pList, int downstream);
2986 void move(list_t **pList, int downstream)
2988 list_t *list, *flipper;
2990 assert(pList != NULL);
2991 list = *pList; /* a>B>c>d */
2992 assert(list != NULL); /* not an empty list */
2993 if (downstream == 0)
2994 flipper = list->prev; /* flipper is A */
2996 flipper = list->next; /* flipper is C */
2997 /* ensure there is somewhere to go in the direction we're heading */
2998 assert(flipper != NULL);
2999 /* remove the current element */
3000 data = pop(&list); /* data=B, a>C>d */
3001 /* and put it back in in the correct spot */
3003 if (downstream == 0) {
3004 push(&list, data, 0); /* b>A>c>d */
3005 list = list->prev; /* B>a>c>d */
3007 push(&list, data, 1); /* a>C>b>d */
3008 list = list->next; /* a>c>B>d */
3014 [[sort]] sorts a list from smallest to largest according to a user
3015 specified comparison.
3016 <<comparison function typedef>>=
3017 typedef int (comparison_fn_t)(void *A, void *B);
3020 <<int comparison function>>=
3021 int int_comparison(void *A, void *B)
3026 if (a > b) return 1;
3027 else if (a == b) return 0;
3032 <<sort declaration>>=
3033 void sort(list_t **pList, comparison_fn_t *comp);
3035 Here we use the bubble sort algorithm.
3037 void sort(list_t **pList, comparison_fn_t *comp)
3041 assert(pList != NULL);
3043 assert(list != NULL); /* not an empty list */
3044 while (stable == 0) {
3047 while (list->next != NULL) {
3048 if ((*comp)(list->d, list->next->d) > 0) {
3050 move(&list, 1 /* downstream */);
3060 [[list_index]] finds the location of [[data]] in [[list]]. Returns
3061 [[-1]] if [[data]] is not in [[list]].
3062 <<list index declaration>>=
3063 int list_index(list_t *list, void *data);
3067 int list_index(list_t *list, void *data)
3071 while (list != NULL) {
3072 if (list->d == data) return i;
3081 [[list_element]] returns the data bound to the $i^\text{th}$ node of the list.
3082 <<list element declaration>>=
3083 void *list_element(list_t *list, int i);
3086 void *list_element(list_t *list, int i)
3090 while (list != NULL) {
3091 if (i == j) return list->d;
3101 [[unique]] removes duplicates from a sorted list according to a user
3102 specified comparison. Don't do this unless you have other pointers
3103 any dynamically allocated data in your list, because [[unique]] just
3104 drops any non-unique data without freeing it.
3105 <<unique declaration>>=
3106 void unique(list_t **pList, comparison_fn_t *comp);
3109 void unique(list_t **pList, comparison_fn_t *comp)
3112 assert(pList != NULL);
3114 assert(list != NULL); /* not an empty list */
3116 while (list->next != NULL) {
3117 if ((*comp)(list->d, list->next->d) == 0) {
3126 [[list_print]] prints a list to a [[FILE *]] stream.
3127 <<list print declaration>>=
3128 void list_print(FILE *file, list_t *list, char *list_name);
3131 void list_print(FILE *file, list_t *list, char *list_name)
3133 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
3135 while (list != NULL) {
3136 fprintf(file, " %p", list->d);
3139 fprintf(file, "\n");
3143 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
3144 <<create and destroy node>>=
3145 list_t *create_node(void *data)
3147 list_t *ret = (list_t *)malloc(sizeof(list_t));
3148 assert(ret != NULL);
3155 void destroy_node(list_t *node)
3161 The user must free the data pointed to by the node on their own.
3163 \subsection{List implementation}
3173 #include <assert.h> /* assert() */
3174 #include <stdlib.h> /* malloc(), free(), rand() */
3175 #include <stdio.h> /* fprintf(), stdout, FILE */
3176 #include "global.h" /* destroy_data_func_t */
3179 \subsection{List unit tests}
3181 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3183 <<list check includes>>
3185 <<main check program>>
3188 <<list check includes>>=
3189 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3190 #include <stdio.h> /* FILE */
3196 <<list test suite>>=
3199 <<list suite definition>>
3202 <<list suite definition>>=
3203 Suite *test_suite (void)
3205 Suite *s = suite_create ("list");
3206 <<push test case defs>>
3207 <<pop test case defs>>
3209 <<push test case adds>>
3210 <<pop test case adds>>
3216 START_TEST(test_push)
3221 push(&list, (void *)i, 1);
3222 fail_unless(list == head(list), NULL);
3223 fail_unless( (int)list->d == 0 );
3224 for (i=0; i<n; i++) {
3225 /* we expect 0, n-1, n-2, ..., 1 */
3228 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
3234 <<push test case defs>>=
3235 TCase *tc_push = tcase_create("push");
3238 <<push test case adds>>=
3239 tcase_add_test(tc_push, test_push);
3240 suite_add_tcase(s, tc_push);
3245 <<pop test case defs>>=
3247 <<pop test case adds>>=
3250 \section{Function string evaluation}
3252 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).
3253 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
3254 The solution is to run a scripting language as a subprocess accessed via pipes.
3255 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
3257 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
3258 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
3259 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.
3260 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
3262 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.
3263 Then you can either statically or dynamically link to those hardcoded functions.
3264 While much less flexible, this approach would be a fairly simple-to-implement fallback.
3266 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
3267 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
3270 #ifndef STRING_EVAL_H
3271 #define STRING_EVAL_H
3273 <<string eval setup declaration>>
3274 <<string eval function declaration>>
3275 <<string eval teardown declaration>>
3276 #endif /* STRING_EVAL_H */
3279 <<string eval module makefile lines>>=
3280 NW_SAWSIM_MODS += string_eval
3281 CHECK_BINS += check_string_eval
3282 check_string_eval_MODS =
3285 For an introduction to POSIX process control, see\\
3286 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
3287 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
3288 and of course, the relavent [[man]] pages.
3290 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
3291 [[execvp]] replaces the calling process' program with a new program.
3292 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
3293 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
3294 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
3295 The new program needs command line arguments, just like it would if you were running it from a shell.
3296 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
3297 with the final array entry being a [[NULL]] pointer.
3299 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
3300 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
3301 <<string eval subprocess definitions>>=
3302 #define SUBPROCESS "python"
3303 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
3304 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
3305 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
3307 The [[i]] option lets Python know that it should run in interactive mode.
3308 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
3309 In interactive mode, python acts on every instruction as soon as it is recieved.
3310 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.
3311 %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.
3313 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
3314 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
3315 [[fork]] returns two copies of the same program, executing the original code.
3316 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
3317 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
3319 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.
3320 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
3321 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
3322 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
3323 <<string eval pipe definitions>>=
3324 #define PIPE_READ 0 /* the end you read from */
3325 #define PIPE_WRITE 1 /* the end you write to */
3327 #define STDIN 0 /* initial index of stdin pair */
3328 #define STDOUT 2 /* initial index of stdout pair */
3331 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
3333 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.
3335 <<string eval setup declaration>>=
3336 extern void string_eval_setup(FILE **pIn, FILE **pOut);
3338 <<string eval setup definition>>=
3339 void string_eval_setup(FILE **pIn, FILE **pOut)
3342 int pfd[4]; /* file descriptors for stdin and stdout */
3344 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
3345 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
3347 pid = fork(); /* split process into two copies */
3349 if (pid == -1) { /* fork error */
3350 perror("fork error");
3352 } else if (pid == 0) { /* child */
3353 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
3354 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
3355 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
3356 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
3357 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
3358 perror("exec error"); /* exec shouldn't return */
3360 } else { /* parent */
3361 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
3362 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
3363 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
3364 if ( *pIn == NULL ) {
3365 perror("fdopen (in)");
3368 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
3369 if ( *pOut == NULL ) {
3370 perror("fdopen (out)");
3377 To use the evaluating subprocess, we just pipe in our command, and read out the results.
3378 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
3379 <<string eval function declaration>>=
3380 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
3382 <<string eval function definition>>=
3383 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
3386 rc = fprintf(in, "%s", input);
3387 assert(rc == strlen(input));
3390 alarm(1); /* set a one second timeout on the read */
3391 assert( fgets(output, buflen, out) != NULL );
3392 alarm(0); /* cancel the timeout */
3393 //fprintf(stderr, "eval: %s ----> %s", input, output);
3396 The [[alarm]] calls set and clear a timeout on the returned output.
3397 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.
3398 This protects against invalid input for which a line of output is not printed to [[stdout]].
3399 Other invalid inputs (e.g.\ those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
3400 If you are getting strange results, check your python code seperately. TODO, better error handling.
3402 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
3403 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
3404 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
3405 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
3406 <<string eval teardown declaration>>=
3407 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
3410 <<string eval teardown definition>>=
3411 void string_eval_teardown(FILE **pIn, FILE **pOut)
3416 /* redirect Python's stderr */
3417 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
3421 assert( fclose(*pIn) == 0 );
3423 assert( fclose(*pOut) == 0 );
3426 /* wait for python to exit */
3428 pid = wait(&stat_loc);
3435 if (WIFEXITED(stat_loc)) {
3436 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
3437 } else if (WIFSIGNALED(stat_loc)) {
3438 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
3443 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
3445 \subsection{String evaluation implementation}
3449 <<string eval includes>>
3450 #include "string_eval.h"
3451 <<string eval internal definitions>>
3452 <<string eval functions>>
3455 <<string eval includes>>=
3456 #include <assert.h> /* assert() */
3457 #include <stdlib.h> /* NULL */
3458 #include <stdio.h> /* fprintf(), stdout, fdopen() */
3459 #include <string.h> /* strlen() */
3460 #include <sys/types.h> /* pid_t */
3461 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
3462 #include <wait.h> /* wait() */
3465 <<string eval internal definitions>>=
3466 <<string eval subprocess definitions>>
3467 <<string eval pipe definitions>>
3470 <<string eval functions>>=
3471 <<string eval setup definition>>
3472 <<string eval function definition>>
3473 <<string eval teardown definition>>
3476 \subsection{String evaluation unit tests}
3478 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3479 <<check-string-eval.c>>=
3480 <<string eval check includes>>
3481 <<string comparison definition and globals>>
3482 <<string eval test suite>>
3483 <<main check program>>
3486 <<string eval check includes>>=
3487 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3488 #include <stdio.h> /* printf() */
3489 #include <assert.h> /* assert() */
3490 #include <string.h> /* strlen() */
3491 #include <signal.h> /* SIGKILL */
3493 #include "string_eval.h"
3496 <<string eval test suite>>=
3497 <<string eval tests>>
3498 <<string eval suite definition>>
3501 <<string eval suite definition>>=
3502 Suite *test_suite (void)
3504 Suite *s = suite_create ("string eval");
3505 <<string eval test case defs>>
3507 <<string eval test case add>>
3512 <<string eval tests>>=
3513 START_TEST(test_setup_teardown)
3516 string_eval_setup(&in, &out);
3517 string_eval_teardown(&in, &out);
3520 START_TEST(test_invalid_command)
3523 char input[80], output[80]={};
3524 string_eval_setup(&in, &out);
3525 sprintf(input, "print ABCDefg\n");
3526 string_eval(in, out, input, 80, output);
3527 string_eval_teardown(&in, &out);
3530 START_TEST(test_simple_eval)
3533 char input[80], output[80]={};
3534 string_eval_setup(&in, &out);
3535 sprintf(input, "print 3+4*5\n");
3536 string_eval(in, out, input, 80, output);
3537 fail_unless(STRMATCH(output,"23\n"), NULL);
3538 string_eval_teardown(&in, &out);
3541 START_TEST(test_multiple_evals)
3544 char input[80], output[80]={};
3545 string_eval_setup(&in, &out);
3546 sprintf(input, "print 3+4*5\n");
3547 string_eval(in, out, input, 80, output);
3548 fail_unless(STRMATCH(output,"23\n"), NULL);
3549 sprintf(input, "print (3**2 + 4**2)**0.5\n");
3550 string_eval(in, out, input, 80, output);
3551 fail_unless(STRMATCH(output,"5.0\n"), NULL);
3552 string_eval_teardown(&in, &out);
3555 START_TEST(test_eval_with_state)
3558 char input[80], output[80]={};
3559 string_eval_setup(&in, &out);
3560 sprintf(input, "print 3+4*5\n");
3561 fprintf(in, "A = 3\n");
3562 sprintf(input, "print A*3\n");
3563 string_eval(in, out, input, 80, output);
3564 fail_unless(STRMATCH(output,"9\n"), NULL);
3565 string_eval_teardown(&in, &out);
3569 <<string eval test case defs>>=
3570 TCase *tc_string_eval = tcase_create("string_eval");
3572 <<string eval test case add>>=
3573 tcase_add_test(tc_string_eval, test_setup_teardown);
3574 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
3575 tcase_add_test(tc_string_eval, test_simple_eval);
3576 tcase_add_test(tc_string_eval, test_multiple_evals);
3577 tcase_add_test(tc_string_eval, test_eval_with_state);
3578 suite_add_tcase(s, tc_string_eval);
3582 \section{Accelerating function evaluation}
3584 My first version-0.3 code was running very slowly.
3585 With the options suggested in the help
3586 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
3587 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
3588 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
3589 That's only 15 calls per solution, so the search algorithm seems reasonable.
3590 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
3595 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
3597 #endif /* ACCEL_K_H */
3600 <<accel k module makefile lines>>=
3601 NW_SAWSIM_MODS += accel_k
3602 #CHECK_BINS += check_accel_k
3603 check_accel_k_MODS =
3607 #include <assert.h> /* assert() */
3608 #include <stdlib.h> /* realloc(), free(), NULL */
3609 #include "global.h" /* environment_t */
3610 #include "k_model.h" /* k_func_t */
3611 #include "interp.h" /* interp_* */
3612 #include "accel_k.h"
3614 typedef struct accel_k_struct {
3615 interp_table_t *itable;
3621 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3622 static int num_accels = 0;
3623 static accel_k_t *accels=NULL;
3625 /* Wrap k in the standard f(x) acceleration form */
3626 static double k_wrap(double F, void *params)
3628 accel_k_t *p = (accel_k_t *)params;
3629 return (*p->k)(F, p->env, p->params);
3632 static int k_tol(double FA, double kA, double FB, double kB)
3635 if (FB-FA > 1e-12) {
3636 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3637 return 1; /* unnacceptable */
3639 //printf("acceptable tol\n");
3640 return 0; /* acceptable */
3644 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3647 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3648 assert(accels != NULL);
3649 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3651 accels[i].env = env;
3652 accels[i].params = params;
3659 for (i=0; i<num_accels; i++)
3660 interp_table_free(accels[i].itable);
3662 if (accels) free(accels);
3666 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3669 for (i=0; i<num_accels; i++) {
3670 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3671 /* We've already setup this function.
3672 * WARNING: we're assuming param and env are the same. */
3673 return interp_table_eval(accels[i].itable, accels+i, F);
3676 /* set up a new acceleration environment */
3677 i = add_accel_k(k, env, params);
3678 return interp_table_eval(accels[i].itable, accels+i, F);
3682 \section{Tension models}
3683 \label{sec.tension_models}
3685 TODO: tension model intro.
3686 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.
3688 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3689 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]].
3691 <<tension-model.h>>=
3692 #ifndef TENSION_MODEL_H
3693 #define TENSION_MODEL_H
3695 <<tension handler types>>
3696 <<constant tension model declarations>>
3697 <<hooke tension model declarations>>
3698 <<worm-like chain tension model declarations>>
3699 <<freely-jointed chain tension model declarations>>
3700 <<piston tension model declarations>>
3701 <<find tension definitions>>
3702 <<tension model global declarations>>
3703 #endif /* TENSION_MODEL_H */
3706 <<tension model module makefile lines>>=
3707 NW_SAWSIM_MODS += tension_model
3708 #CHECK_BINS += check_tension_model
3709 check_tension_model_MODS =
3711 <<tension model utils makefile lines>>=
3712 TENSION_MODEL_MODS = tension_model parse list tension_balance
3713 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3714 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3715 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3716 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
3717 notangle -Rtension-model-utils.c $< > $@
3718 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
3719 gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3720 clean_tension_model_utils :
3721 rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3722 @ The pipe symbol [[|]] marks the prerequisits following it (in this
3723 case, the directories) as ‘order-only’ prerequisites. The timestamp
3724 on these prerequisits does not effect whether the rules are executed.
3725 This is appropriate for directories, where we don't need to recompile
3726 after an unrelated has been added to the directory, but only when the
3727 source prerequisites change. See the [[make]] documentation for more
3729 (\url{http://www.gnu.org/software/make/manual/html_node/Prerequisite-Types.html}).
3732 \label{sec.null_tension}
3734 For unstretchable domains.
3736 <<null tension model getopt>>=
3737 {NULL, NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3740 \subsection{Constant}
3741 \label{sec.const_tension}
3743 An infinitely stretchable domain providing a constant tension.
3744 Obviously this cannot be inverted, so you can't put this domain in
3745 series with any other domains. We include it mostly for testing
3748 <<constant tension functions>>=
3749 <<constant tension function>>
3750 <<constant tension parameter create/destroy functions>>
3753 <<constant tension model declarations>>=
3754 <<constant tension function declaration>>
3755 <<constant tension parameter create/destroy function declarations>>
3756 <<constant tension model global declarations>>
3760 <<constant tension function declaration>>=
3761 extern double const_tension_handler(double x, void *pdata);
3763 <<constant tension function>>=
3764 double const_tension_handler(double x, void *pdata)
3766 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3767 list_t *list = data->group_tension_params;
3772 assert(list != NULL); /* empty active group?! */
3773 F = ((const_tension_param_t *)list->d)->F;
3775 fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3777 while (list != NULL) {
3778 assert(((const_tension_param_t *)list->d)->F == F);
3786 <<constant tension structure>>=
3787 typedef struct const_tension_param_struct {
3788 double F; /* tension (force) in N */
3789 } const_tension_param_t;
3793 <<constant tension parameter create/destroy function declarations>>=
3794 extern void *string_create_const_tension_param_t(char **param_strings);
3795 extern void destroy_const_tension_param_t(void *p);
3798 <<constant tension parameter create/destroy functions>>=
3799 const_tension_param_t *create_const_tension_param_t(double F)
3801 const_tension_param_t *ret
3802 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3803 assert(ret != NULL);
3808 void *string_create_const_tension_param_t(char **param_strings)
3810 return create_const_tension_param_t(
3811 safe_strtod(param_strings[0],__FUNCTION__));
3814 void destroy_const_tension_param_t(void *p)
3822 <<constant tension model global declarations>>=
3823 extern int num_const_tension_params;
3824 extern char *const_tension_param_descriptions[];
3825 extern char const_tension_param_string[];
3828 <<constant tension model globals>>=
3829 int num_const_tension_params = 1;
3830 char *const_tension_param_descriptions[] = {"tension F, N"};
3831 char const_tension_param_string[] = "0";
3834 <<constant tension model getopt>>=
3835 {&const_tension_handler, NULL, "const", "an infinitely stretchable domain with constant tension", num_const_tension_params, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t}
3841 <<hooke functions>>=
3842 <<internal hooke functions>>
3844 <<inverse hooke handler>>
3845 <<hooke parameter create/destroy functions>>
3848 <<hooke tension model declarations>>=
3849 <<hooke tension function declaration>>
3850 <<hooke tension parameter create/destroy function declarations>>
3851 <<hooke tension model global declarations>>
3855 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
3856 The behavior of a series of springs $k_i$ in series is given by
3858 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
3860 For a simple proof, see Appendix \ref{sec.math_hooke}.
3862 <<hooke tension function declaration>>=
3863 extern double hooke_handler(double x, void *pdata);
3864 extern double inverse_hooke_handler(double F, void *pdata);
3867 First we define a function that computes the effective spring constant
3868 (stored in a single [[hooke_param_t]]) for the entire group.
3869 <<internal hooke functions>>=
3870 static void hooke_param_grouper(tension_handler_data_t *data,
3871 hooke_param_t *param) {
3872 list_t *list = data->group_tension_params;
3876 assert(list != NULL); /* empty active group?! */
3877 while (list != NULL) {
3878 assert( ((hooke_param_t *)list->d)->k > 0 );
3879 k += 1.0/ ((hooke_param_t *)list->d)->k;
3881 fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
3882 __FUNCTION__, ((hooke_param_t *)list->d)->k);
3888 fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N (from data %p)\n",
3889 __FUNCTION__, k, x, k*x, data);
3896 Everyone knows Hooke's law: $F=kx$.
3898 double hooke_handler(double x, void *pdata)
3900 hooke_param_t param = {0};
3903 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
3909 Inverting Hooke's law gives $x=F/k$.
3910 <<inverse hooke handler>>=
3911 double inverse_hooke_handler(double F, void *pdata)
3913 hooke_param_t param = {0};
3916 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
3922 Not much to keep track of for a single effective spring, just the
3923 spring constant $k$.
3924 <<hooke structure>>=
3925 typedef struct hooke_param_struct {
3926 double k; /* spring constant in N/m */
3931 We wrap up our Hooke implementation with some book-keeping functions.
3932 <<hooke tension parameter create/destroy function declarations>>=
3933 extern void *string_create_hooke_param_t(char **param_strings);
3934 extern void destroy_hooke_param_t(void *p);
3938 <<hooke parameter create/destroy functions>>=
3939 hooke_param_t *create_hooke_param_t(double k)
3941 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
3942 assert(ret != NULL);
3947 void *string_create_hooke_param_t(char **param_strings)
3949 return create_hooke_param_t(safe_strtod(param_strings[0], __FUNCTION__));
3952 void destroy_hooke_param_t(void *p)
3959 <<hooke tension model global declarations>>=
3960 extern int num_hooke_params;
3961 extern char *hooke_param_descriptions[];
3962 extern char hooke_param_string[];
3965 <<hooke tension model globals>>=
3966 int num_hooke_params = 1;
3967 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
3968 char hooke_param_string[]="0.05";
3971 <<hooke tension model getopt>>=
3972 {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}
3975 \subsection{Worm-like chain}
3978 We can model several unfolded domains with a single worm-like chain.
3979 <<worm-like chain functions>>=
3980 <<internal worm-like chain functions>>
3981 <<worm-like chain handler>>
3982 <<inverse worm-like chain handler>>
3983 <<worm-like chain create/destroy functions>>
3986 <<worm-like chain tension model declarations>>=
3988 <<worm-like chain handler declaration>>
3989 <<inverse worm-like chain handler declaration>>
3990 <<worm-like chain create/destroy function declarations>>
3991 <<worm-like chain tension model global declarations>>
3995 <<internal worm-like chain functions>>=
3996 <<worm-like chain function>>
3997 <<inverse worm-like chain function>>
3998 <<worm-like chain parameter grouper>>
4001 The combination of all unfolded domains is modeled as a worm like chain,
4002 with the force $F_{\text{WLC}}$ approximately given by
4004 F_{\text{WLC}} \approx \frac{k_B T}{p}
4005 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
4007 where $x$ is the distance between the fixed ends,
4008 $k_B$ is Boltzmann's constant,
4009 $T$ is the absolute temperature,
4010 $p$ is the persistence length, and
4011 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
4012 <<worm-like chain function>>=
4013 static double wlc(double x, double T, double p, double L)
4015 double xL=0.0; /* uses global k_B */
4017 assert(T > 0); assert(p > 0); assert(L > 0);
4018 if (x >= L) return HUGE_VAL;
4019 xL = x/L; /* unitless */
4020 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
4021 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
4022 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
4027 Inverting the approximate WLC is fairly straightforward, if a bit tedious.
4028 Letting $x_L \equiv x/L$ and $F_T \equiv Fp/k_B T$ we have
4030 F_T &\approx \frac{1}{4} \p({\frac{1}{(1-x_L)^2} - 1}) + x_L \\
4031 F_T (1-x_L)^2 &\approx \frac{1}{4} - \frac{1}{4} (1-x_L)^2 + x_L(1-x_L)^2 \\
4032 0 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) (1-x_L)^2 + x_L(1-x_L)^2 \\
4033 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) \p({1-2x_L+x_L^2})
4034 + x_L - 2x_L^2 + x_L^3 \\
4035 &\approx \p[{\frac{1}{4} - \p({F_T + \frac{1}{4}})}]
4036 + x_L \p[{2F_T + \frac{1}{2} + 1}]
4037 + x_L^2 \p[{-\p({F_T + \frac{1}{4}}) - 2}]
4040 + x_L \p({2F_T + \frac{3}{2}})
4041 - x_L^2 \p({F_T + \frac{9}{4}})
4043 &\approx ax_L^3 + bx_L^2 + cx_L + d
4045 where $a \equiv 1$, $b \equiv -F_T - 9/4$, $c \equiv 2F_T + 3/2$, and
4047 % From \citet{wikipedia_cubic_function} the discriminant
4049 % \Delta &= -4b^3d + b^2c^2 - 4ac^3 + 18abcd - 27a^2d^2 \\
4050 % &= -4F_T\p({F_T + \frac{9}{4}})^3
4051 % + \p({F_T + \frac{9}{4}})^2 \cdot \p({2F_T + \frac{3}{2}})^2
4052 % - 4 \p({2F_T + \frac{3}{2}})^3
4053 % + 18F_T \p({F_T + \frac{9}{4}}) \cdot \p({2F_T + \frac{3}{2}})
4055 % &= -4F_T \p({F_T^3 + \frac{27}{4}F_T^2 +\frac{243}{16}F_T +\frac{729}{64}})
4056 %% a^3 + 3a^2b + 3ab^2 + b^3
4057 % + \p({F_T^2 + \frac{9}{2}F_T + \frac{81}{16}})
4058 % \cdot \p({4F_T^2 + 6F_T + \frac{9}{4}}) \\
4059 % &\qquad - 4 \p({8F_T^3 + 18F_T^2 + \frac{27}{2}F_T + \frac{27}{8}})
4060 %% a^3 + 3a^2b + 3ab^2 + b^3
4061 % + 18F_T \p({2F_T^2 + 6F_T + \frac{27}{8}})
4063 % &= -4F_T^4 - 27F_T^3 - \frac{243}{4}F_T^2 - \frac{729}{16}F_T
4064 % + 4F_T^4 + 24F_T^3 +\frac{99}{2}F_T^2 +\frac{81}{2}F_T +\frac{729}{64}\\
4065 % &\qquad - 32F_T^3 - 72F_T^2 - 54F_T - \frac{27}{2}
4066 % + 36F_T^3 + 108F_T^2 + \frac{243}{4}F_T
4068 % &= F_T^4 \p({-4 + 4}) + F_T^3 \p({-27 + 24 - 32 + 36})
4069 % + F_T^2 \p({-\frac{243}{4} + \frac{99}{2} - 72 + 108 -27}) \\
4070 % &\qquad + F_T \p({-\frac{729}{16} + \frac{81}{2} - 54 + \frac{243}{4}})
4071 % \p({\frac{729}{64} - \frac{27}{2}}) \\
4072 % &= F_T^3 - \frac{9}{4}F_T^2 + \frac{27}{16}F_T - \frac{135}{64}
4074 We can plug these coefficients into the GSL cubic solver to invert the
4075 WLC\citep{galassi05}. The applicable root is always the one which
4076 comes before the singularity, which will be the smallest real root.
4077 <<inverse worm-like chain function>>=
4078 static double inverse_wlc(double F, double T, double p, double L)
4080 double FT = F*p/(k_B*T); /* uses global k_B */
4081 double xL0, xL1, xL2;
4084 assert(T > 0); assert(p > 0); assert(L > 0);
4087 num_roots = gsl_poly_solve_cubic(-(FT+2.25),2*FT+1.5,-FT, &xL0, &xL1, &xL2);
4088 assert(num_roots > 0); assert(xL0 >= -DOUBLE_PRECISION); assert(xL0 < 1);
4089 if (xL0 < 0) xL0 = 0.0;
4095 First we define a function that computes the effective WLC parameters
4096 (stored in a single [[wlc_param_t]]) for the entire group.
4097 <<worm-like chain parameter grouper>>=
4098 static void wlc_param_grouper(tension_handler_data_t *data,
4099 wlc_param_t *param) {
4100 list_t *list = data->group_tension_params;
4101 double p=0.0, L=0.0;
4104 assert(list != NULL); /* empty active group?! */
4105 p = ((wlc_param_t *)list->d)->p;
4106 while (list != NULL) {
4107 /* enforce consistent p */
4108 assert( ((wlc_param_t *)list->d)->p == p);
4109 L += ((wlc_param_t *)list->d)->L;
4111 fprintf(stderr, "%s: WLC section %g m long in series\n",
4112 __FUNCTION__, ((wlc_param_t *)list->d)->L);
4117 fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
4118 __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
4126 <<worm-like chain handler declaration>>=
4127 extern double wlc_handler(double x, void *pdata);
4130 This model requires that each unfolded domain in the group have the
4131 same persistence length.
4132 <<worm-like chain handler>>=
4133 double wlc_handler(double x, void *pdata)
4135 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4136 wlc_param_t param = {0};
4139 wlc_param_grouper(data, ¶m);
4140 return wlc(x, data->env->T, param.p, param.L);
4145 <<inverse worm-like chain handler declaration>>=
4146 extern double inverse_wlc_handler(double F, void *pdata);
4149 <<inverse worm-like chain handler>>=
4150 double inverse_wlc_handler(double F, void *pdata)
4152 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4153 wlc_param_t param = {0};
4155 wlc_param_grouper(data, ¶m);
4156 return inverse_wlc(F, data->env->T, param.p, param.L);
4161 <<worm-like chain structure>>=
4162 typedef struct wlc_param_struct {
4163 double p; /* WLC persistence length (m) */
4164 double L; /* WLC contour length (m) */
4168 <<worm-like chain create/destroy function declarations>>=
4169 extern void *string_create_wlc_param_t(char **param_strings);
4170 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
4173 <<worm-like chain create/destroy functions>>=
4174 wlc_param_t *create_wlc_param_t(double p, double L)
4176 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
4177 assert(ret != NULL);
4183 void *string_create_wlc_param_t(char **param_strings)
4185 return create_wlc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4186 safe_strtod(param_strings[1], __FUNCTION__));
4189 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
4197 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.
4198 TODO (now we copy the string before parsing, but still do this for clarity).
4200 <<valid string write code>>=
4201 char string[] = "abc";
4204 <<valid string write code>>=
4205 char *string = "abc";
4209 <<worm-like chain tension model global declarations>>=
4210 extern int num_wlc_params;
4211 extern char *wlc_param_descriptions[];
4212 extern char wlc_param_string[];
4214 <<worm-like chain tension model globals>>=
4215 int num_wlc_params = 2;
4216 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
4217 char wlc_param_string[]="0.39e-9,28.4e-9";
4220 <<worm-like chain tension model getopt>>=
4221 {&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}
4223 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
4225 \subsection{Freely-jointed chain}
4228 <<freely-jointed chain functions>>=
4229 <<freely-jointed chain function>>
4230 <<freely-jointed chain handler>>
4231 <<freely-jointed chain create/destroy functions>>
4234 <<freely-jointed chain tension model declarations>>=
4235 <<freely-jointed chain handler declaration>>
4236 <<freely-jointed chain create/destroy function declarations>>
4237 <<freely-jointed chain tension model global declarations>>
4241 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
4242 The tension of a chain of $N$ such links is given by
4244 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
4246 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}.
4247 <<freely-jointed chain function>>=
4248 double langevin(double x, void *params)
4251 return 1.0/tanh(x) - 1/x;
4253 <<one dimensional bracket declaration>>
4254 <<one dimensional solve declaration>>
4255 double inv_langevin(double x)
4257 double lb=0.0, ub=1.0, ret;
4258 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
4259 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
4263 double fjc(double x, double T, double l, int N)
4265 double L = l*(double)N;
4266 if (x == 0) return 0;
4267 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
4268 return k_B*T/l * inv_langevin(x/L);
4272 <<freely-jointed chain handler declaration>>=
4273 extern double fjc_handler(double x, void *pdata);
4275 <<freely-jointed chain handler>>=
4276 double fjc_handler(double x, void *pdata)
4278 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4279 list_t *list = data->group_tension_params;
4284 assert(list != NULL); /* empty active group?! */
4285 l = ((fjc_param_t *)list->d)->l;
4286 while (list != NULL) {
4287 /* enforce consistent l */
4288 assert( ((fjc_param_t *)list->d)->l == l);
4289 N += ((fjc_param_t *)list->d)->N;
4291 fprintf(stderr, "%s: FJC section %d links long in series\n",
4292 __FUNCTION__, ((fjc_param_t *)list->d)->N);
4297 fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
4298 __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
4300 return fjc(x, data->env->T, l, N);
4304 <<freely-jointed chain structure>>=
4305 typedef struct fjc_param_struct {
4306 double l; /* FJC link length (m) */
4307 int N; /* FJC number of links */
4311 <<freely-jointed chain create/destroy function declarations>>=
4312 extern void *string_create_fjc_param_t(char **param_strings);
4313 extern void destroy_fjc_param_t(void *p);
4316 <<freely-jointed chain create/destroy functions>>=
4317 fjc_param_t *create_fjc_param_t(double l, double N)
4319 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
4320 assert(ret != NULL);
4326 void *string_create_fjc_param_t(char **param_strings)
4328 return create_fjc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4329 safe_strtod(param_strings[1], __FUNCTION__));
4332 void destroy_fjc_param_t(void *p)
4339 <<freely-jointed chain tension model global declarations>>=
4340 extern int num_fjc_params;
4341 extern char *fjc_param_descriptions[];
4342 extern char fjc_param_string[];
4345 <<freely-jointed chain tension model globals>>=
4346 int num_fjc_params = 2;
4347 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
4348 char fjc_param_string[]="0.5e-9,200";
4351 <<freely-jointed chain tension model getopt>>=
4352 {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}
4356 \label{sec.piston_tension}
4358 An infinitely stretchable domain with no tension for extensions less
4359 than a particular threshold $L$, and infinite tension for $x>L$. The
4360 tension at $x=L$ is undefined, but may be any positive value. The
4361 model is called the ``piston'' model because it resembles a
4362 frictionless piston sliding in a rigid cylinder of length $L$.
4364 Because the tension at $x=L$ is undefined, the user must make sure
4365 domains with this tension model are never the initial active group,
4366 because it would break [[tension_balance()]] in [[find_tension()]]
4367 (see Section \ref{sec.tension_balance}).
4369 <<piston tension functions>>=
4370 <<piston tension parameter grouper>>
4371 <<piston tension handler>>
4372 <<inverse piston tension handler>>
4373 <<piston tension parameter create/destroy functions>>
4376 <<piston tension model declarations>>=
4377 <<piston tension handler declaration>>
4378 <<inverse piston tension handler declaration>>
4379 <<piston tension parameter create/destroy function declarations>>
4380 <<piston tension model global declarations>>
4384 First we define a function that computes the effective piston parameters
4385 (stored in a single [[piston_tension_param_t]]) for the entire group.
4386 <<piston tension parameter grouper>>=
4387 static void piston_tension_param_grouper(tension_handler_data_t *data,
4388 piston_tension_param_t *param) {
4389 list_t *list = data->group_tension_params;
4393 assert(list != NULL); /* empty active group?! */
4394 while (list != NULL) {
4395 L += ((piston_tension_param_t *)list->d)->L;
4401 <<piston tension handler declaration>>=
4402 extern double piston_tension_handler(double x, void *pdata);
4404 <<piston tension handler>>=
4405 double piston_tension_handler(double x, void *pdata)
4407 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4408 piston_tension_param_t param = {0};
4411 piston_tension_param_grouper(data, ¶m);
4412 if (x <= param.L) F = 0;
4415 fprintf(stderr, "%s: x %g, L %g, tension %g N\n", __FUNCTION__, x, param.L, F);
4422 <<inverse piston tension handler declaration>>=
4423 extern double inverse_piston_tension_handler(double x, void *pdata);
4425 <<inverse piston tension handler>>=
4426 double inverse_piston_tension_handler(double F, void *pdata)
4428 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4429 piston_tension_param_t param = {0};
4432 piston_tension_param_grouper(data, ¶m);
4433 if (F == 0.0) return 0;
4439 <<piston tension structure>>=
4440 typedef struct piston_tension_param_struct {
4441 double L; /* length of piston in m */
4442 } piston_tension_param_t;
4446 <<piston tension parameter create/destroy function declarations>>=
4447 extern void *string_create_piston_tension_param_t(char **param_strings);
4448 extern void destroy_piston_tension_param_t(void *p);
4452 <<piston tension parameter create/destroy functions>>=
4453 piston_tension_param_t *create_piston_tension_param_t(double L)
4455 piston_tension_param_t *ret
4456 = (piston_tension_param_t *) malloc(sizeof(piston_tension_param_t));
4457 assert(ret != NULL);
4462 void *string_create_piston_tension_param_t(char **param_strings)
4464 return create_piston_tension_param_t(
4465 safe_strtod(param_strings[0],__FUNCTION__));
4468 void destroy_piston_tension_param_t(void *p)
4476 <<piston tension model global declarations>>=
4477 extern int num_piston_tension_params;
4478 extern char *piston_tension_param_descriptions[];
4479 extern char piston_tension_param_string[];
4483 <<piston tension model globals>>=
4484 int num_piston_tension_params = 1;
4485 char *piston_tension_param_descriptions[] = {"length L, m"};
4486 char piston_tension_param_string[] = "0";
4490 <<piston tension model getopt>>=
4491 {&piston_tension_handler, &inverse_piston_tension_handler, "piston", "a domain that slides without tension for x<L, but has infinite tension for x>L.", num_piston_tension_params, piston_tension_param_descriptions, piston_tension_param_string, &string_create_piston_tension_param_t, &destroy_piston_tension_param_t}
4494 \subsection{Tension model implementation}
4496 <<tension-model.c>>=
4499 <<tension model includes>>
4500 #include "tension_model.h"
4501 <<tension model internal definitions>>
4502 <<tension model globals>>
4503 <<tension model functions>>
4506 <<tension model includes>>=
4507 #include <assert.h> /* assert() */
4508 #include <stdlib.h> /* NULL, strto*() */
4509 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4510 #include <stdio.h> /* fprintf(), stdout */
4511 #include <errno.h> /* errno, ERANGE */
4514 #include "tension_balance.h" /* oneD_*() */
4517 <<tension model internal definitions>>=
4518 <<constant tension structure>>
4520 <<worm-like chain structure>>
4521 <<freely-jointed chain structure>>
4522 <<piston tension structure>>
4525 <<tension model globals>>=
4526 <<tension handler array global>>
4527 <<constant tension model globals>>
4528 <<hooke tension model globals>>
4529 <<worm-like chain tension model globals>>
4530 <<freely-jointed chain tension model globals>>
4531 <<piston tension model globals>>
4534 <<tension model functions>>=
4536 <<constant tension functions>>
4538 <<worm-like chain functions>>
4539 <<freely-jointed chain functions>>
4540 <<piston tension functions>>
4543 \subsection{Utilities}
4545 The tension models can get complicated, and users may want to reassure
4546 themselves that this portion of the input physics are functioning properly. The
4547 stand-alone program [[tension_model_utils]] demonstrates the tension model
4548 interface without the overhead of having to understand a full simulation.
4549 [[tension_model_utils]] takes command line model arguments like the full
4550 simulation, and prints $F(x)$ data to the screen for the selected model over a
4553 <<tension-model-utils.c>>=
4555 <<tension model utility includes>>
4556 <<tension model utility definitions>>
4557 <<tension model utility getopt functions>>
4559 <<tension model utility main>>
4562 <<tension model utility main>>=
4563 int main(int argc, char **argv)
4565 <<tension model getopt array definition>>
4566 tension_model_getopt_t *model = NULL;
4570 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
4571 tension_handler_data_t tdata;
4572 int num_param_args; /* for INIT_MODEL() */
4573 char **param_args; /* for INIT_MODEL() */
4575 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model,
4576 &Fmax, &Xmax, &flags);
4578 if (flags & VFLAG) {
4579 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4581 INIT_MODEL("utility", model, model->params, params);
4582 tdata.group_tension_params = NULL;
4583 push(&tdata.group_tension_params, params, 1);
4585 tdata.persist = NULL;
4586 if (model->handler == NULL) {
4587 printf("No tension function for model %s\n", model->name);
4593 printf("#Distance (m)\tForce (N)\n");
4594 for (i=0; i<=N; i++) {
4595 x = Xmax*i/(double)N;
4596 F = (*model->handler)(x, &tdata);
4597 if (F < 0 || F > Fmax) break;
4598 printf("%g\t%g\n", x, F);
4600 if (flags & VFLAG && i <= N) {
4601 /* explain exit condition */
4603 printf("Impossible force %g\n", F);
4605 printf("Reached large force limit %g > %g\n", F, Fmax);
4608 params = pop(&tdata.group_tension_params);
4609 if (model->destructor)
4610 (*model->destructor)(params);
4615 <<tension model utility includes>>=
4618 #include <string.h> /* strlen() */
4619 #include <assert.h> /* assert() */
4620 #include <errno.h> /* errno, ERANGE */
4624 #include "tension_model.h"
4627 <<tension model utility definitions>>=
4628 <<version definition>>
4629 #define VFLAG 1 /* verbose */
4630 <<string comparison definition and globals>>
4631 <<tension model getopt definitions>>
4632 <<initialize model definition>>
4636 <<tension model utility getopt functions>>=
4638 <<index tension model>>
4639 <<tension model utility help>>
4640 <<tension model utility get options>>
4643 <<tension model utility help>>=
4644 void help(char *prog_name,
4646 int n_tension_models, tension_model_getopt_t *tension_models,
4647 int tension_model, double Fmax, double Xmax)
4650 printf("usage: %s [options]\n", prog_name);
4651 printf("Version %s\n\n", VERSION);
4652 printf("A utility for understanding the available tension models\n\n");
4653 printf("Environmental options:\n");
4654 printf("-T T\tTemperature T (currently %g K)\n", env->T);
4655 printf("-C T\tYou can also set the temperature T in Celsius\n");
4656 printf("Model options:\n");
4657 printf("-m model\tFolded tension model (currently %s)\n",
4658 tension_models[tension_model].name);
4659 printf("-a args\tFolded tension model argument string (currently %s)\n",
4660 tension_models[tension_model].params);
4661 printf("F(x) is calculated for a range of x and printed\n");
4662 printf("For example:\n");
4663 printf(" #Distance (m)\tForce (N)\n");
4664 printf(" 123.456\t7.89\n");
4666 printf("-F\tSet the maximum F value for the standard mode F(x) output (currently %g)\n", Fmax);
4667 printf("-X\tSet the maximum x value for the standard mode F(x) output (currently %g)\n", Xmax);
4668 printf("-V\tChange output to verbose mode\n");
4669 printf("-h\tPrint this help and exit\n");
4671 printf("Tension models:\n");
4672 for (i=0; i<n_tension_models; i++) {
4673 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
4674 for (j=0; j < tension_models[i].num_params ; j++)
4675 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
4676 printf("\t\tdefaults: %s\n", tension_models[i].params);
4681 <<tension model utility get options>>=
4682 void get_options(int argc, char **argv, environment_t *env,
4683 int n_tension_models, tension_model_getopt_t *tension_models,
4684 tension_model_getopt_t **model, double *Fmax, double *Xmax,
4685 unsigned int *flags)
4687 char *prog_name = NULL;
4688 char c, options[] = "T:C:m:a:F:X:Vh";
4689 int tension_model=0, num_strings;
4690 extern char *optarg;
4691 extern int optind, optopt, opterr;
4695 /* setup defaults */
4697 prog_name = argv[0];
4698 env->T = 300.0; /* K */
4699 *Fmax = 1e5; /* N */
4700 *Xmax = 1e-6; /* m */
4702 *model = tension_models;
4704 while ((c=getopt(argc, argv, options)) != -1) {
4706 case 'T': env->T = safe_strtod(optarg, "-T"); break;
4707 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
4709 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
4710 *model = tension_models+tension_model;
4713 tension_models[tension_model].params = optarg;
4715 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
4716 case 'X': *Xmax = safe_strtod(optarg, "-X"); break;
4717 case 'V': *flags |= VFLAG; break;
4719 fprintf(stderr, "unrecognized option '%c'\n", optopt);
4720 /* fall through to default case */
4722 help(prog_name, env, n_tension_models, tension_models, tension_model, *Fmax, *Xmax);
4731 \section{Unfolding rate models}
4732 \label{sec.k_models}
4734 We have two main choices for determining $k$: Bell's model, which assumes the
4735 folded domains exist in a single `folded' state and make instantaneous,
4736 stochastic transitions over a free energy barrier; and Kramers' model, which
4737 treats the unfolding as a thermalized diffusion process.
4738 We deal with these options like we dealt with the different tension models: by bundling all model-specific
4739 parameters into structures, with model specific functions for determining $k$.
4741 <<k func definition>>=
4742 typedef double k_func_t(double F, environment_t *env, void *params);
4745 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
4746 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]].
4748 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
4749 because \LaTeX\ doesn't like underscores outside math mode.
4750 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
4751 You could use spaces instead (`\verb+ +'),
4752 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
4753 which I don't like as much.
4759 <<k func definition>>
4760 <<null k declarations>>
4761 <<const k declarations>>
4762 <<bell k declarations>>
4763 <<kbell k declarations>>
4764 <<kramers k declarations>>
4765 <<kramers integ k declarations>>
4766 #endif /* K_MODEL_H */
4769 <<k model module makefile lines>>=
4770 NW_SAWSIM_MODS += k_model
4771 CHECK_BINS += check_k_model
4772 check_k_model_MODS = parse string_eval
4774 [[check_k_model]] also depends on the parse module.
4776 <<k model utils makefile lines>>=
4777 K_MODEL_MODS = k_model parse string_eval
4778 K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
4779 $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h)
4780 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
4781 notangle -Rk-model-utils.c $< > $@
4782 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
4783 gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
4784 clean_k_model_utils :
4785 rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
4791 For domains that are never folded.
4793 <<null k declarations>>=
4795 <<null k functions>>=
4800 <<null k model getopt>>=
4801 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
4804 \subsection{Constant rate model}
4807 This is a pretty boring model, but it is usefull for testing the engine.
4811 where $k_0$ is the constant unfolding rate.
4813 <<const k functions>>=
4814 <<const k function>>
4815 <<const k structure create/destroy functions>>
4818 <<const k declarations>>=
4819 <<const k function declaration>>
4820 <<const k structure create/destroy function declarations>>
4821 <<const k global declarations>>
4825 <<const k structure>>=
4826 typedef struct const_k_param_struct {
4827 double knot; /* transition rate k_0 (frac population per s) */
4831 <<const k function declaration>>=
4832 double const_k(double F, environment_t *env, void *const_k_params);
4834 <<const k function>>=
4835 double const_k(double F, environment_t *env, void *const_k_params)
4836 { /* Returns the rate constant k in frac pop per s. */
4837 const_k_param_t *p = (const_k_param_t *) const_k_params;
4839 assert(p->knot > 0);
4844 <<const k structure create/destroy function declarations>>=
4845 void *string_create_const_k_param_t(char **param_strings);
4846 void destroy_const_k_param_t(void *p);
4849 <<const k structure create/destroy functions>>=
4850 const_k_param_t *create_const_k_param_t(double knot)
4852 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
4853 assert(ret != NULL);
4858 void *string_create_const_k_param_t(char **param_strings)
4860 return create_const_k_param_t(safe_strtod(param_strings[0], __FUNCTION__));
4863 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
4870 <<const k global declarations>>=
4871 extern int num_const_k_params;
4872 extern char *const_k_param_descriptions[];
4873 extern char const_k_param_string[];
4876 <<const k globals>>=
4877 int num_const_k_params = 1;
4878 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
4879 char const_k_param_string[]="1";
4882 <<const k model getopt>>=
4883 {"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}
4886 \subsection{Bell's model}
4889 Quantitatively, Bell's model gives $k$ as
4891 k = k_0 \cdot e^{F dx / k_B T} \;,
4893 where $k_0$ is the unforced unfolding rate,
4894 $F$ is the applied unfolding force,
4895 $dx$ is the distance to the transition state, and
4896 $k_B T$ is the thermal energy\citep{schlierf06}.
4898 <<bell k functions>>=
4900 <<bell k structure create/destroy functions>>
4903 <<bell k declarations>>=
4904 <<bell k function declaration>>
4905 <<bell k structure create/destroy function declarations>>
4906 <<bell k global declarations>>
4910 <<bell k structure>>=
4911 typedef struct bell_param_struct {
4912 double knot; /* transition rate k_0 (frac population per s) */
4913 double dx; /* `distance to transition state' (nm) */
4917 <<bell k function declaration>>=
4918 double bell_k(double F, environment_t *env, void *bell_params);
4920 <<bell k function>>=
4921 double bell_k(double F, environment_t *env, void *bell_params)
4922 { /* Returns the rate constant k in frac pop per s.
4923 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4924 * uses global k_B in J/K */
4925 bell_param_t *p = (bell_param_t *) bell_params;
4926 assert(F >= 0); assert(env->T > 0);
4928 assert(p->knot > 0); //assert(p->dx >= 0); /* p < 0 for refolding reactions, etc. */
4930 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
4931 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
4932 p->knot * exp(F*p->dx / (k_B*env->T) ));
4933 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
4935 return p->knot * exp(F*p->dx / (k_B*env->T) );
4939 <<bell k structure create/destroy function declarations>>=
4940 void *string_create_bell_param_t(char **param_strings);
4941 void destroy_bell_param_t(void *p);
4944 <<bell k structure create/destroy functions>>=
4945 bell_param_t *create_bell_param_t(double knot, double dx)
4947 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
4948 assert(ret != NULL);
4954 void *string_create_bell_param_t(char **param_strings)
4956 return create_bell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4957 safe_strtod(param_strings[1], __FUNCTION__));
4960 void destroy_bell_param_t(void *p /* bell_param_t * */)
4967 <<bell k global declarations>>=
4968 extern int num_bell_params;
4969 extern char *bell_param_descriptions[];
4970 extern char bell_param_string[];
4974 int num_bell_params = 2;
4975 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4976 char bell_param_string[]="3.3e-4,0.25e-9";
4979 <<bell k model getopt>>=
4980 {"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}
4982 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4983 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4986 \subsection{Linker stiffness corrected Bell model}
4989 Walton et. al showed that the Bell model constant force approximation
4990 doesn't hold for biotin-streptavadin binding, and I extended their
4991 results to I27 for the 2009 Biophysical Society Annual
4992 Meeting\cite{walton08,king09}. More details to follow when I get done
4993 with the conference\ldots
4995 We adjust Bell's model to
4997 k = k_0 \cdot e^{\frac{F dx - \frac{1}{2} \kappa dx^2}{k_B T}} \;,
4999 where $k_0$ is the unforced unfolding rate,
5000 $F$ is the applied unfolding force,
5001 $\kappa$ is the effective spring constant,
5002 $dx$ is the distance to the transition state, and
5003 $k_B T$ is the thermal energy\citep{schlierf06}.
5005 <<kbell k functions>>=
5006 <<kbell k function>>
5007 <<kbell k structure create/destroy functions>>
5010 <<kbell k declarations>>=
5011 <<kbell k function declaration>>
5012 <<kbell k structure create/destroy function declarations>>
5013 <<kbell k global declarations>>
5017 <<kbell k structure>>=
5018 typedef struct kbell_param_struct {
5019 double knot; /* transition rate k_0 (frac population per s) */
5020 double dx; /* `distance to transition state' (nm) */
5024 <<kbell k function declaration>>=
5025 double kbell_k(double F, environment_t *env, void *kbell_params);
5027 <<kbell k function>>=
5028 double kbell_k(double F, environment_t *env, void *kbell_params)
5029 { /* Returns the rate constant k in frac pop per s.
5030 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5031 * uses global k_B in J/K */
5032 kbell_param_t *p = (kbell_param_t *) kbell_params;
5033 assert(F >= 0); assert(env->T > 0); assert(env->stiffness >= 0);
5035 assert(p->knot > 0); assert(p->dx >= 0);
5036 return p->knot * exp((F*p->dx - 0.5*env->stiffness*(p->dx)*(p->dx))/(k_B*env->T) );
5040 <<kbell k structure create/destroy function declarations>>=
5041 void *string_create_kbell_param_t(char **param_strings);
5042 void destroy_kbell_param_t(void *p);
5045 <<kbell k structure create/destroy functions>>=
5046 kbell_param_t *create_kbell_param_t(double knot, double dx)
5048 kbell_param_t *ret = (kbell_param_t *) malloc(sizeof(kbell_param_t));
5049 assert(ret != NULL);
5055 void *string_create_kbell_param_t(char **param_strings)
5057 return create_kbell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5058 safe_strtod(param_strings[1], __FUNCTION__));
5061 void destroy_kbell_param_t(void *p /* kbell_param_t * */)
5068 <<kbell k global declarations>>=
5069 extern int num_kbell_params;
5070 extern char *kbell_param_descriptions[];
5071 extern char kbell_param_string[];
5074 <<kbell k globals>>=
5075 int num_kbell_params = 2;
5076 char *kbell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
5077 char kbell_param_string[]="3.3e-4,0.25e-9";
5080 <<kbell k model getopt>>=
5081 {"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}
5083 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
5084 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5087 \subsection{Kramer's model}
5090 Kramer's model gives $k$ as
5092 % \frac{1}{k} = \frac{1}{D}
5093 % \int_{x_\text{min}}^{x_\text{max}}
5094 % e^{\frac{-U_F(x)}{k_B T}}
5095 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5098 %where $D$ is the diffusion constant,
5099 %$U_F(x)$ is the free energy along the unfolding coordinate, and
5100 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
5101 % before the minimum and after the maximum, respectively \citep{schlierf06}.
5103 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
5105 where $D$ is the diffusion constant,
5107 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
5109 are length scales where
5110 $x_c(F)$ is the location of the bound state and
5111 $x_{ts}(F)$ is the location of the transition state,
5112 $E(F,x)$ is the free energy, and
5113 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
5114 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
5115 \citet{evans97} Eqn.~3).
5117 <<kramers k functions>>=
5118 <<kramers k function>>
5119 <<kramers k structure create/destroy functions>>
5122 <<kramers k declarations>>=
5123 <<kramers k function declaration>>
5124 <<kramers k structure create/destroy function declarations>>
5125 <<kramers k global declarations>>
5129 <<kramers k structure>>=
5130 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
5131 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
5133 typedef struct kramers_param_struct {
5134 double D; /* diffusion rate (um/s) */
5141 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
5142 //kramers_x_func_t *xts; /* function returning transition x (nm) */
5143 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
5144 //kramers_E_func_t *E; /* function returning E (J) */
5145 //double *E_params; /* parametrize the energy landscape */
5146 //int n_E_params; /* length of E_params array */
5150 <<kramers k function declaration>>=
5151 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
5152 extern double kramers_k(double F, environment_t *env, void *kramers_params);
5154 <<kramers k function>>=
5155 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
5157 static char input[160]={0};
5158 static char output[80]={0};
5159 /* setup the environment */
5160 fprintf(in, "F = %g; T = %g\n", F, T);
5161 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
5162 string_eval(in, out, input, 80, output);
5163 return safe_strtod(output, __FUNCTION__);
5166 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
5168 static char input[160]={0};
5169 static char output[80]={0};
5170 /* setup the environment */
5171 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
5172 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
5173 string_eval(in, out, input, 80, output);
5174 return safe_strtod(output, __FUNCTION__);
5177 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
5179 kramers_param_t *p = (kramers_param_t *) kramers_params;
5180 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
5183 double kramers_k(double F, environment_t *env, void *kramers_params)
5184 { /* Returns the rate constant k in frac pop per s.
5185 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5186 * uses global k_B in J/K */
5187 kramers_param_t *p = (kramers_param_t *) kramers_params;
5188 double kbT, xc, xts, lc, lts, Eb;
5189 assert(F >= 0); assert(env->T > 0);
5192 assert(p->xc != NULL); assert(p->xts != NULL);
5193 assert(p->ddEdxx != NULL); assert(p->E != NULL);
5194 assert(p->in != NULL); assert(p->out != NULL);
5196 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
5197 if (xc == -1.0) return -1.0; /* error (Too much force) */
5198 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
5199 if (xc == -1.0) return -1.0; /* error (Too much force) */
5200 lc = sqrt(2.0*M_PI*kbT /
5201 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
5202 lts = sqrt(-2.0*M_PI*kbT/
5203 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
5204 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
5205 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
5206 //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));
5207 return p->D/(lc*lts) * exp(-Eb/kbT);
5211 <<kramers k structure create/destroy function declarations>>=
5212 void *string_create_kramers_param_t(char **param_strings);
5213 void destroy_kramers_param_t(void *p);
5216 <<kramers k structure create/destroy functions>>=
5217 kramers_param_t *create_kramers_param_t(double D,
5218 char *xc_fn, char *xts_fn,
5219 char *E_fn, char *ddEdxx_fn,
5220 char *global_define_string)
5221 // kramers_x_func_t *xc, kramers_x_func_t *xts,
5222 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
5223 // double *E_params, int n_E_params)
5225 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
5226 assert(ret != NULL);
5227 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5228 assert(ret->xc != NULL);
5229 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5230 assert(ret->xts != NULL);
5231 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
5232 assert(ret->E != NULL);
5233 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
5234 assert(ret->ddEdxx != NULL);
5236 strcpy(ret->xc, xc_fn);
5237 strcpy(ret->xts, xts_fn);
5238 strcpy(ret->E, E_fn);
5239 strcpy(ret->ddEdxx, ddEdxx_fn);
5240 //ret->E_params = E_params;
5241 //ret->n_E_params = n_E_params;
5242 string_eval_setup(&ret->in, &ret->out);
5243 fprintf(ret->in, "kB = %g\n", k_B);
5244 fprintf(ret->in, "%s\n", global_define_string);
5248 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
5249 void *string_create_kramers_param_t(char **param_strings)
5251 return create_kramers_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5259 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
5261 kramers_param_t *pk = (kramers_param_t *)p;
5263 string_eval_teardown(&pk->in, &pk->out);
5265 // free(pk->E_params);
5271 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
5272 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.
5273 We follow \citet{shillcock98} and use
5275 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
5276 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
5279 where TODO, variable meanings.%\citep{shillcock98}.
5280 The first and second derivatives of $E(F,E_b,x,x_s)$ are
5282 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\\
5283 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
5285 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
5286 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
5287 The roots are therefor at
5289 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
5290 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
5291 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
5294 <<kramers k global declarations>>=
5295 extern int num_kramers_params;
5296 extern char *kramers_param_descriptions[];
5297 extern char kramers_param_string[];
5300 <<kramers k globals>>=
5301 int num_kramers_params = 6;
5302 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"};
5303 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)}";
5305 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
5306 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
5307 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
5308 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.
5309 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
5310 It works so long as [[val_1]] is not [[false]].
5312 <<kramers k model getopt>>=
5313 {"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}
5316 \subsection{Kramer's model (natural cubic splines)}
5317 \label{sec.kramers_integ}
5319 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
5320 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
5321 \citet{schlierf06} Eqn.~2)
5323 \frac{1}{k} = \frac{1}{D}
5324 \int_{x_\text{min}}^{x_\text{max}}
5325 e^{\frac{U_F(x)}{k_B T}}
5326 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5329 where $D$ is the diffusion constant,
5330 $U_F(x)$ is the free energy along the unfolding coordinate, and
5331 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
5332 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
5334 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
5335 interpolating between them with cubic splines.
5336 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
5338 <<kramers integ k functions>>=
5339 <<spline functions>>
5340 <<kramers integ A k functions>>
5341 <<kramers integ B k functions>>
5342 <<kramers integ k function>>
5343 <<kramers integ k structure create/destroy functions>>
5346 <<kramers integ k declarations>>=
5347 <<kramers integ k function declaration>>
5348 <<kramers integ k structure create/destroy function declarations>>
5349 <<kramers integ k global declarations>>
5353 <<kramers integ k structure>>=
5354 typedef double func_t(double x, void *params);
5355 typedef struct kramers_integ_param_struct {
5356 double D; /* diffusion rate (um/s) */
5357 double x_min; /* integration bounds */
5359 func_t *E; /* function returning E (J) */
5360 void *E_params; /* parametrize the energy landscape */
5361 destroy_data_func_t *destroy_E_params;
5363 double F; /* for passing into GSL evaluations */
5365 } kramers_integ_param_t;
5368 <<spline param structure>>=
5369 typedef struct spline_param_struct {
5370 int n_knots; /* natural cubic spline knots for U(x) */
5373 gsl_spline *spline; /* GSL spline parameters */
5374 gsl_interp_accel *acc; /* GSL spline acceleration data */
5378 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
5380 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5382 <<kramers integ A k functions>>=
5383 double kramers_integ_k_integrandA(double x, void *params)
5385 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5386 double U = (*p->E)(x, p->E_params);
5387 double Fx = p->F * x;
5388 double kbT = k_B * p->env->T;
5389 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
5390 return exp(-(U-Fx)/kbT);
5392 double kramers_integ_k_integralA(double x, void *params)
5394 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5396 double result, abserr;
5397 size_t neval; /* for qng */
5398 static gsl_integration_workspace *w = NULL;
5399 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5400 f.function = &kramers_integ_k_integrandA;
5402 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5403 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5404 //fprintf(stderr, "integralA = %g\n", result);
5410 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
5411 \int_{x_\text{min}}^{x_\text{max}}
5412 e^{\frac{U_F(x)}{k_B T}}
5413 \text{Integral}_A(x)
5416 <<kramers integ B k functions>>=
5417 double kramers_integ_k_integrandB(double x, void *params)
5419 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5420 double U = (*p->E)(x, p->E_params);
5421 double Fx = p->F * x;
5422 double kbT = k_B * p->env->T;
5423 double intA = kramers_integ_k_integralA(x, params);
5424 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
5425 return exp((U-Fx)/kbT)*intA;
5427 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
5429 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5431 double result, abserr;
5432 size_t neval; /* for qng */
5433 static gsl_integration_workspace *w = NULL;
5434 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5435 f.function = &kramers_integ_k_integrandB;
5439 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5440 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5441 //fprintf(stderr, "integralB = %g\n", result);
5446 With the integrals taken care of, evaluating $k$ is simply
5448 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
5450 <<kramers integ k function declaration>>=
5451 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
5453 <<kramers integ k function>>=
5454 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
5455 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
5456 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
5457 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
5461 <<kramers integ k structure create/destroy function declarations>>=
5462 void *string_create_kramers_integ_param_t(char **param_strings);
5463 void destroy_kramers_integ_param_t(void *p);
5466 <<kramers integ k structure create/destroy functions>>=
5467 kramers_integ_param_t *create_kramers_integ_param_t(double D,
5468 double xmin, double xmax, func_t *E,
5470 destroy_data_func_t *destroy_E_params)
5472 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
5473 assert(ret != NULL);
5478 ret->E_params = E_params;
5479 ret->destroy_E_params = destroy_E_params;
5483 void *string_create_kramers_integ_param_t(char **param_strings)
5485 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
5486 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
5487 void *E_params = string_create_spline_param_t(param_strings+1);
5488 return create_kramers_integ_param_t(
5489 safe_strtod(param_strings[0], __FUNCTION__),
5490 safe_strtod(param_strings[2], __FUNCTION__),
5491 safe_strtod(param_strings[3], __FUNCTION__),
5492 &spline_eval, E_params, destroy_spline_param_t);
5495 void destroy_kramers_integ_param_t(void *params)
5497 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5500 (*p->destroy_E_params)(p->E_params);
5506 Finally we have the GSL spline wrappers:
5507 <<spline functions>>=
5509 <<create/destroy spline>>
5512 <<spline function>>=
5513 double spline_eval(double x, void *spline_params)
5515 spline_param_t *p = (spline_param_t *)(spline_params);
5516 return gsl_spline_eval(p->spline, x, p->acc);
5520 <<create/destroy spline>>=
5521 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
5523 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
5524 assert(ret != NULL);
5525 ret->n_knots = n_knots;
5528 ret->acc = gsl_interp_accel_alloc();
5529 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
5530 gsl_spline_init(ret->spline, x, y, n_knots);
5534 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
5535 void *string_create_spline_param_t(char **param_strings)
5539 double *x=NULL, *y=NULL;
5540 /* break into ordered pairs */
5541 parse_list_string(param_strings[0], SEP, '(', ')',
5542 &num_knots, &knot_coords);
5543 x = (double *)malloc(sizeof(double)*num_knots);
5545 y = (double *)malloc(sizeof(double)*num_knots);
5547 for (i=0; i<num_knots; i++) {
5548 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
5549 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
5552 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
5554 free_parsed_list(num_knots, knot_coords);
5555 return create_spline_param_t(num_knots, x, y);
5558 void destroy_spline_param_t(void *params)
5560 spline_param_t *p = (spline_param_t *)params;
5563 gsl_spline_free(p->spline);
5565 gsl_interp_accel_free(p->acc);
5575 <<kramers integ k global declarations>>=
5576 extern int num_kramers_integ_params;
5577 extern char *kramers_integ_param_descriptions[];
5578 extern char kramers_integ_param_string[];
5581 <<kramers integ k globals>>=
5582 int num_kramers_integ_params = 4;
5583 char *kramers_integ_param_descriptions[] = {
5584 "diffusion D, m^2 / s",
5585 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
5586 "starting integration bound x_min, m",
5587 "ending integration bound x_max, m"};
5588 //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";
5589 //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";
5590 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";
5591 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
5592 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
5593 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
5596 <<kramers integ k model getopt>>=
5597 {"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}
5599 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
5600 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5602 \subsection{Unfolding model implementation}
5606 <<k model includes>>
5607 #include "k_model.h"
5608 <<k model internal definitions>>
5610 <<k model functions>>
5613 <<k model includes>>=
5614 #include <assert.h> /* assert() */
5615 #include <stdlib.h> /* NULL, malloc(), strto*() */
5616 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
5617 #include <stdio.h> /* fprintf(), stdout */
5618 #include <string.h> /* strlen(), strcpy() */
5619 #include <errno.h> /* errno, ERANGE */
5620 #include <gsl/gsl_integration.h>
5621 #include <gsl/gsl_spline.h>
5626 <<k model internal definitions>>=
5627 <<const k structure>>
5628 <<bell k structure>>
5629 <<kbell k structure>>
5630 <<kramers k structure>>
5631 <<kramers integ k structure>>
5632 <<spline param structure>>
5635 <<k model globals>>=
5640 <<kramers k globals>>
5641 <<kramers integ k globals>>
5644 <<k model functions>>=
5646 <<null k functions>>
5647 <<const k functions>>
5648 <<bell k functions>>
5649 <<kbell k functions>>
5650 <<kramers k functions>>
5651 <<kramers integ k functions>>
5654 \subsection{Unfolding model unit tests}
5656 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
5657 <<check-k-model.c>>=
5658 <<k model check includes>>
5659 <<check relative error>>
5661 <<k model test suite>>
5662 <<main check program>>
5665 <<k model check includes>>=
5666 #include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
5667 #include <stdio.h> /* sprintf() */
5668 #include <assert.h> /* assert() */
5669 #include <math.h> /* exp() */
5670 #include <errno.h> /* errno, ERANGE */
5673 #include "k_model.h"
5676 <<k model test suite>>=
5680 <<k model suite definition>>
5683 <<k model suite definition>>=
5684 Suite *test_suite (void)
5686 Suite *s = suite_create ("k model");
5687 <<const k test case defs>>
5688 <<bell k test case defs>>
5690 <<const k test case adds>>
5691 <<bell k test case adds>>
5696 \subsubsection{Constant}
5698 <<const k test case defs>>=
5699 TCase *tc_const_k = tcase_create("Constant k");
5702 <<const k test case adds>>=
5703 tcase_add_test(tc_const_k, test_const_k_create_destroy);
5704 tcase_add_test(tc_const_k, test_const_k_over_range);
5705 suite_add_tcase(s, tc_const_k);
5710 START_TEST(test_const_k_create_destroy)
5712 char *knot[] = {"1","2","3","4","5","6"};
5713 char *params[] = {knot[0]};
5716 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5717 params[0] = knot[i];
5718 p = string_create_const_k_param_t(params);
5719 destroy_const_k_param_t(p);
5724 START_TEST(test_const_k_over_range)
5726 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
5727 char *knot[] = {"1","2","3","4","5","6"};
5728 char *params[] = {knot[0]};
5731 char param_string[80];
5734 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5735 params[0] = knot[i];
5736 p = string_create_const_k_param_t(params);
5737 for ( F=Fm; F<FM; F+=dF ) {
5738 fail_unless(const_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
5739 "const_k(%g, %g, &{%s}) = %g != %s",
5740 F, T, knot[i], const_k(F, &env, p), knot[i]);
5742 destroy_const_k_param_t(p);
5748 \subsubsection{Bell}
5750 <<bell k test case defs>>=
5751 TCase *tc_bell = tcase_create("Bell's k");
5754 <<bell k test case adds>>=
5755 tcase_add_test(tc_bell, test_bell_k_create_destroy);
5756 tcase_add_test(tc_bell, test_bell_k_at_zero);
5757 tcase_add_test(tc_bell, test_bell_k_at_one);
5758 suite_add_tcase(s, tc_bell);
5762 START_TEST(test_bell_k_create_destroy)
5764 char *knot[] = {"1","2","3","4","5","6"};
5766 char *params[] = {knot[0], dx};
5769 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5770 params[0] = knot[i];
5771 p = string_create_bell_param_t(params);
5772 destroy_bell_param_t(p);
5777 START_TEST(test_bell_k_at_zero)
5779 double F=0.0, T=300.0;
5780 char *knot[] = {"1","2","3","4","5","6"};
5782 char *params[] = {knot[0], dx};
5785 char param_string[80];
5788 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5789 params[0] = knot[i];
5790 p = string_create_bell_param_t(params);
5791 fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
5792 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
5793 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
5794 destroy_bell_param_t(p);
5799 START_TEST(test_bell_k_at_one)
5802 char *knot[] = {"1","2","3","4","5","6"};
5804 char *params[] = {knot[0], dx};
5805 double F= k_B*T/safe_strtod(dx, __FUNCTION__);
5808 char param_string[80];
5811 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5812 params[0] = knot[i];
5813 p = string_create_bell_param_t(params);
5814 CHECK_ERR(1e-6, safe_strtod(knot[i], __FUNCTION__)*exp(1.0), bell_k(F, &env, p));
5815 //fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__)*exp(1.0),
5816 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
5817 // F, T, knot[i], dx, bell_k(F, &env, p), safe_strtod(knot[i], __FUNCTION__)*exp(1.0));
5818 destroy_bell_param_t(p);
5824 <<kramers k tests>>=
5827 <<kramers k test case def>>=
5830 <<kramers k test case add>>=
5833 <<k model function tests>>=
5834 <<check relative error>>
5836 START_TEST(test_something)
5838 double k=5, x=3, last_x=2.0, F;
5839 one_dim_fn_t *handlers[] = {&hooke};
5840 void *data[] = {&k};
5841 double xi[] = {0.0};
5843 int new_active_groups = 1;
5844 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
5845 fail_unless(F = k*x, NULL);
5850 \subsection{Utilities}
5852 The unfolding models can get complicated, and are sometimes hard to explain to others.
5853 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
5854 the overhead of having to understand a full simulation.
5855 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
5856 or special mode, where the operation depends on the specific model selected.
5858 <<k-model-utils.c>>=
5860 <<k model utility includes>>
5861 <<k model utility definitions>>
5862 <<k model utility getopt functions>>
5863 <<k model utility multi model E>>
5864 <<k model utility main>>
5867 <<k model utility main>>=
5868 int main(int argc, char **argv)
5870 <<k model getopt array definition>>
5871 k_model_getopt_t *model = NULL;
5876 int num_param_args; /* for INIT_MODEL() */
5877 char **param_args; /* for INIT_MODEL() */
5879 double special_xmin;
5880 double special_xmax;
5881 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
5882 &Fmax, &special_xmin, &special_xmax, &flags);
5883 if (flags & VFLAG) {
5884 printf("#initializing model %s with parameters %s\n", model->name, model->params);
5886 INIT_MODEL("utility", model, model->params, params);
5889 if (model->k == NULL) {
5890 printf("No k function for model %s\n", model->name);
5894 printf("Fmax = %g <= 0\n", Fmax);
5900 printf("#F (N)\tk (%% pop. per s)\n");
5901 for (i=0; i<=N; i++) {
5902 F = Fmax*i/(double)N;
5903 k = (*model->k)(F, &env, params);
5905 printf("%g\t%g\n", F, k);
5910 if (model->k == NULL || model->k == &bell_k) {
5911 printf("No special mode for model %s\n", model->name);
5914 if (special_xmax <= special_xmin) {
5915 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
5919 double Xrng=(special_xmax-special_xmin), x, E;
5921 printf("#x (m)\tE (J)\n");
5922 for (i=0; i<=N; i++) {
5923 x = special_xmin + Xrng*i/(double)N;
5924 E = multi_model_E(model->k, params, &env, x);
5925 printf("%g\t%g\n", x, E);
5932 if (model->destructor)
5933 (*model->destructor)(params);
5938 <<k model utility multi model E>>=
5939 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
5941 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
5943 if (k == kramers_integ_k)
5944 E = (*p->E)(x, p->E_params);
5946 E = kramers_E(0, env, model_params, x);
5952 <<k model utility includes>>=
5955 #include <string.h> /* strlen() */
5956 #include <assert.h> /* assert() */
5957 #include <errno.h> /* errno, ERANGE */
5960 #include "string_eval.h"
5961 #include "k_model.h"
5964 <<k model utility definitions>>=
5965 <<version definition>>
5966 #define VFLAG 1 /* verbose */
5967 enum mode_t {M_K_OF_F, M_SPECIAL};
5968 <<string comparison definition and globals>>
5969 <<k model getopt definitions>>
5970 <<initialize model definition>>
5971 <<kramers k structure>>
5972 <<kramers integ k structure>>
5976 <<k model utility getopt functions>>=
5979 <<k model utility help>>
5980 <<k model utility get options>>
5983 <<k model utility help>>=
5984 void help(char *prog_name,
5986 int n_k_models, k_model_getopt_t *k_models,
5987 int k_model, double Fmax, double special_xmin, double special_xmax)
5990 printf("usage: %s [options]\n", prog_name);
5991 printf("Version %s\n\n", VERSION);
5992 printf("A utility for understanding the available unfolding models\n\n");
5993 printf("Environmental options:\n");
5994 printf("-T T\tTemperature T (currently %g K)\n", env->T);
5995 printf("-C T\tYou can also set the temperature T in Celsius\n");
5996 printf("Model options:\n");
5997 printf("-k model\tTransition rate model (currently %s)\n",
5998 k_models[k_model].name);
5999 printf("-K args\tTransition rate model argument string (currently %s)\n",
6000 k_models[k_model].params);
6001 printf("There are two output modes. In standard mode, k(F) is printed\n");
6002 printf("For example:\n");
6003 printf(" #Force (N)\tk (% pop. per s)\n");
6004 printf(" 123.456\t7.89\n");
6006 printf("In special mode, the output depends on the model.\n");
6007 printf("For models defining an energy function E(x), that function is printed\n");
6008 printf(" #Distance (m)\tE (J)\n");
6009 printf(" 0.0012\t0.0034\n");
6011 printf("-m\tChange output to standard mode\n");
6012 printf("-M\tChange output to special mode\n");
6013 printf("-F\tSet the maximum F value for the standard mode k(F) output (currently %g)\n", Fmax);
6014 printf("-x\tSet the minimum x value for the special mode E(x) output (currently %g)\n", special_xmin);
6015 printf("-X\tSet the maximum x value for the special mode E(x) output (currently %g)\n", special_xmax);
6016 printf("-V\tChange output to verbose mode\n");
6017 printf("-h\tPrint this help and exit\n");
6019 printf("Unfolding rate models:\n");
6020 for (i=0; i<n_k_models; i++) {
6021 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
6022 for (j=0; j < k_models[i].num_params ; j++)
6023 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
6024 printf("\t\tdefaults: %s\n", k_models[i].params);
6029 <<k model utility get options>>=
6030 void get_options(int argc, char **argv, environment_t *env,
6031 int n_k_models, k_model_getopt_t *k_models,
6032 enum mode_t *mode, k_model_getopt_t **model,
6033 double *Fmax, double *special_xmin, double *special_xmax,
6034 unsigned int *flags)
6036 char *prog_name = NULL;
6037 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
6039 extern char *optarg;
6040 extern int optind, optopt, opterr;
6044 /* setup defaults */
6046 prog_name = argv[0];
6047 env->T = 300.0; /* K */
6051 *Fmax = 1e-9; /* N */
6052 *special_xmax = 1e-8;
6053 *special_xmin = 0.0;
6055 while ((c=getopt(argc, argv, options)) != -1) {
6057 case 'T': env->T = safe_strtod(optarg, "-T"); break;
6058 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
6060 k_model = index_k_model(n_k_models, k_models, optarg);
6061 *model = k_models+k_model;
6064 k_models[k_model].params = optarg;
6066 case 'm': *mode = M_K_OF_F; break;
6067 case 'M': *mode = M_SPECIAL; break;
6068 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
6069 case 'x': *special_xmin = safe_strtod(optarg, "-x"); break;
6070 case 'X': *special_xmax = safe_strtod(optarg, "-X"); break;
6071 case 'V': *flags |= VFLAG; break;
6073 fprintf(stderr, "unrecognized option '%c'\n", optopt);
6074 /* fall through to default case */
6076 help(prog_name, env, n_k_models, k_models, k_model, *Fmax, *special_xmin, *special_xmax);
6085 \section{\LaTeX\ documentation}
6087 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
6088 The comment blocks are extracted (with nicely formatted code blocks), using
6089 <<latex makefile lines>>=
6090 $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6091 noweave -latex -index -delay $< > $@
6092 $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR)
6096 <<latex makefile lines>>=
6097 $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
6099 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
6100 $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
6101 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
6102 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
6103 mv $(BUILD_DIR)/sawsim.pdf $@
6105 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
6106 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
6107 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
6109 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
6110 <<latex makefile lines>>=
6112 rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
6113 $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
6114 $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.toc \
6115 $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib
6116 clean_latex : semi-clean_latex
6117 rm -f $(DOC_DIR)/sawsim.pdf
6122 [[make]] is a common utility on *nix systems for managing dependencies while building software.
6123 In this case, we have several \emph{targets} that we'd like to build:
6124 the main simulation program \prog;
6125 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
6126 the documentation [[sawsim.pdf]];
6127 and the [[Makefile]] itself.
6128 Besides the generated files, there is the utility target
6129 [[clean]] for removing all generated files except [[Makefile]].
6131 % [[dist]] for generating a distributable tar file.
6133 Extract the makefile with
6134 `[[notangle -Rmakefile src/sawsim.nw | sed 's/ /\t/' > Makefile]]'.
6135 The sed is needed because notangle replaces tabs with eight spaces,
6136 but [[make]] needs tabs.
6138 # Decide what directories to put things in
6143 # Note: Cannot use, for example, `./build', because make eats the `./'
6144 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
6146 # Modules (X.c, X.h) defined in the noweb file
6149 # Modules defined outside the noweb file
6150 FREE_SAWSIM_MODS = interp tavl
6152 <<list module makefile lines>>
6153 <<tension balance module makefile lines>>
6154 <<tension model module makefile lines>>
6155 <<k model module makefile lines>>
6156 <<parse module makefile lines>>
6157 <<string eval module makefile lines>>
6158 <<accel k module makefile lines>>
6160 SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS)
6162 # Everything needed for sawsim under one roof. sawsim.c must come first
6163 SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \
6164 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h)
6165 # Libraries needed to compile sawsim
6166 LIBS = gsl gslcblas m
6167 CHECK_LIBS = gsl gslcblas m check
6168 # The non-check binaries generated
6169 BINS = sawsim tension_model_utils k_model_utils sawsim_profile
6172 # Define the major targets
6173 all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
6175 view : $(DOC_DIR)/sawsim.pdf
6177 profile : $(BIN_DIR)/sawsim_profile | $(BIN_DIR)
6178 $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \
6179 -mnull -Mwlc -A0.39e-9,28e-9 -F8
6180 gprof $(BIN_DIR)/sawsim_profile gmon.out > $@
6181 check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim
6182 $(SHELL) -e -c 'for B in $^; do ./$$B; done'
6183 clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \
6184 clean_tension_model_utils \
6185 clean_k_model_utils clean_latex clean_check_sawsim
6186 rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \
6187 $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \
6188 $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
6189 $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
6190 $(BUILD_DIR)/global.h ./gmon.out
6191 $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
6193 # Various builds of sawsim
6194 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
6195 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6196 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
6197 gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6198 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
6199 gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6201 # Create the directories
6202 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
6205 # Copy over the source external to sawsim
6206 # Note: Cannot use, for example, `./build', because make eats the `./'
6207 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
6209 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6213 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6217 ## Basic source generated with noweb
6218 # The central files sawsim.c and global.h...
6219 $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6221 $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6222 notangle -Rglobal.h $< > $@
6223 # ... and the modules
6224 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6225 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6226 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6227 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6228 $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6229 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6230 # Note: I use `_' as a space in my file names, but noweb doesn't like
6231 # that so we use `-' to name the noweb roots and substitute here.
6233 ## Building the unit-test programs
6235 $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
6236 notangle -Rchecks $< > $@
6237 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
6238 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
6239 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR)
6240 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
6241 clean_check_sawsim :
6242 rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
6243 # ... and the modules
6245 $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
6246 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \
6247 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \
6248 $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6249 $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6250 $$(BUILD_DIR)/global.h | $(BIN_DIR)
6251 gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
6252 $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
6254 # todo: clean up the required modules to
6255 $(CHECK_BINS:%=clean_%) :
6256 rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
6258 # Cleaning up the modules
6260 $(SAWSIM_MODS:%=clean_%) :
6261 rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h)
6263 <<tension model utils makefile lines>>
6264 <<k model utils makefile lines>>
6265 <<latex makefile lines>>
6267 Makefile : $(SRC_DIR)/sawsim.nw
6268 notangle -Rmakefile $< | sed 's/ /\t/' > $@
6270 This is fairly self-explanatory. We compile a dynamically linked
6271 version ([[sawsim]]) and a statically linked version
6272 ([[sawsim_static]]). The static version is about 9 times as big, but
6273 it runs without needing dynamic GSL libraries to link to, so it's a
6274 better format for distributing to a cluster for parallel evaluation.
6278 \subsection{Hookean springs in series}
6279 \label{sec.math_hooke}
6282 The effective spring constant for $n$ springs $k_i$ in series is given by
6284 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6290 F &= k_i x_i &&\forall i \le n \\
6291 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
6292 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
6293 F &= k_1 x_1 = k_\text{eff} x \\
6294 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
6295 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6300 \addcontentsline{toc}{section}{References}
6301 \bibliography{sawsim}