1 % sawsim - simulating single molecule mechanical unfolding.
2 % Copyright (C) 2008-2010 William Trevor King
4 % This program is free software: you can redistribute it and/or modify
5 % it under the terms of the GNU General Public License as published by
6 % the Free Software Foundation, either version 3 of the License, or
7 % (at your option) any later version.
9 % This program is distributed in the hope that it will be useful,
10 % but WITHOUT ANY WARRANTY; without even the implied warranty of
11 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 % GNU General Public License for more details.
14 % You should have received a copy of the GNU General Public License
15 % along with this program. If not, see <http://www.gnu.org/licenses/>.
17 % The author may be contacted at <wking@drexel.edu> on the Internet, or
18 % write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
19 % Philadelphia PA 19104, USA.
21 %%%%%%%%%%%%%%%%%%%%%%%%% Start LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%
22 % we have our own preamble,
23 % use `noweave -delay` to not print noweb's automatic one
24 % -*- mode: Noweb; noweb-code-mode: c-mode -*-
25 \documentclass[letterpaper, 11pt]{article}
28 \noweboptions{smallcode}
30 \usepackage{url} % needed for natbib for underscores in urls?
31 \usepackage[breaklinks,colorlinks]{hyperref} % for internal links, \phantomsection
32 % breaklinks breaks long links
33 % colorlinks colors link text instead of boxing it
34 \usepackage{amsmath} % for \text{}, \xrightarrow[sub]{super}
35 \usepackage{amsthm} % for theorem and proof environments
36 \newtheorem{theorem}{Theorem}
38 \usepackage{subfig} % Support side-by-side figures
39 \usepackage{pgf} % Fancy graphics
40 \usepackage{tikz} % A nice, inline PGF frontend
41 \usetikzlibrary{automata} % Graph-theory library
43 \usepackage{doc} % BibTeX
44 \usepackage[super,sort&compress]{natbib} % fancy citation extensions
45 % super selects citations in superscript mode
46 % sort&compress automatically sorts and compresses compound citations (\cite{a,b,...})
48 \usepackage[mathletters]{ucs} % use math fonts to allow Greek UTF-8 in the text
49 \usepackage[utf8x]{inputenc} % allow UTF-8 characters
51 %\bibliographystyle{ieeetr} % pick the bibliography style, short and sweet
52 %\bibliographystyle{plain} % pick the bibliography style, includes dates
53 \bibliographystyle{plainnat}
54 \defcitealias{sw:noweb}{{\tt noweb}}
55 \defcitealias{galassi05}{GSL}
56 \defcitealias{sw:check}{{\tt check}}
57 \defcitealias{sw:python}{Python}
62 \textheight 9.5in % leave a bit of extra room for page numbers (at top of page)
65 \setlength{\parindent}{0pt}
66 \setlength{\parskip}{5pt}
68 % For some reason, the \maketitle wants to go on it's own page
69 % so we'll just hardcode the following in our first page.
70 %\title{Sawsim: a sawtooth protein unfolding simulator}
71 %\author{W.~Trevor King}
74 \newcommand{\prog}{[[sawsim]]}
75 \newcommand{\lang}{[[C]]}
76 \newcommand{\p}[3]{\left#1 #2 \right#3} % e.g. \p({complicated expr})
77 \newcommand{\Th}{\ensuremath{\sup{\text{th}}}} % ordinal th
78 \newcommand{\U}[1]{\ensuremath{\text{ #1}}} % units: $1\U{m/s}$
80 % single spaced lists, from Alvin Alexander
81 % http://www.devdaily.com/blog/post/latex/control-line-spacing-in-itemize-enumerate-tags/
82 \newenvironment{packed_item}{
84 \setlength{\itemsep}{1pt}
85 \setlength{\parskip}{0pt}
86 \setlength{\parsep}{0pt}
90 \nwfilename{sawsim.nw}
95 Sawsim: a sawtooth protein unfolding simulator \\
103 %%%%%%%%%%%%%%%%%%%%%%%%%%%% End LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%%
105 \section{Introduction}
107 The unfolding of multi-globular protein strings is a stochastic
108 process. In the \prog\ program, we use Monte Carlo methods to
109 simulate this unfolding through an explicit time-stepping approach.
110 We develop a program in \lang\ to simulate probability of unfolding
111 under a constant extension velocity of a chain of globular domains.
113 In order to extract physical meaning from many mechanical unfolding
114 simulations, it is necessary to compare the experimental results
115 against the results of simulations using different models and
116 parameters. The model and parameters that best match the experimental
117 data provide the ‘best’ model for that experiment.
119 Previous mechanical unfolding studies have rolled their own
120 simulations for modelling their experiment, and there has been little
121 exchange and standardization of these simulations between groups.
122 The \prog\ program attempts to fill this gap, providing a flexible,
123 extensible, \emph{community} program, to make it easier for groups
124 to share the burden and increase the transparency of data analysis.
126 ‘Sawsim’ is blend of ‘sawtooth simulation’.
128 \subsection{Related work}
132 \subsection{About this document}
134 This document was written in \citetalias{sw:noweb} with code blocks in
135 \lang\ and comment blocks in \LaTeX.
137 \subsection{System Requirements and Dependencies}
140 \subsection{Change Log}
142 Version 0.0 used only the unfolded domain WLC to determine the tension
143 and had a fixed timestep $dt$.
145 Version 0.1 added adaptive timesteps, adjusting so that the unfolding
146 probability for a given domain was constant.
148 Version 0.2 added Kramers' model unfolding rate calculations, and a
149 switch to select Bell's or Kramers' model. This induced a major
150 rewrite, introducing the tension group abstraction, and a split to
151 multiple \lang\ source files. Also added a unit-test suites for
152 sanity checking, and switched to SI units throughout.
154 Version 0.3 added integrated Kramers' model unfolding rate
155 calculations. Also replaced some of my hand-coded routines with GNU
156 Scientific Library \citepalias{galassi05} calls.
158 Version 0.4 added Python evaluation for the saddle-point Kramers'
159 model unfolding rate calculations, so that the model functions could
160 be controlled from the command line. Also added interpolating lookup
161 tables to accelerate slow unfolding rate model evaluation, and command
162 line control of tension groups.
164 Version 0.5 added the stiffness environmental parameter and it's
165 associated unfolding models.
167 Version 0.6 generalized from two state unfolding to arbitrary
168 $N$-state rate reactions. This allows simulations of
169 unfolding-refolding, metastable intermediates, etc., but required
170 reorganizing the command line interface.
172 Version 0.7 added tension model inverses, which seems to reduce
173 computation time by about a factor of 2. I was expecting more, but
174 I'll take what I can get.
176 Version 0.8 fixed a minor bug in unfolding probability for
177 multi-domain groups. The probability of at least one domain unfolding
178 had been calculated as $P_N=NP_1$, rather than $P_N=1-(1-P_1)^N$.
179 However, for small $P$ the two are equivalent.
181 Version 0.9 added the [[-P]] option to sawsim, setting the target
182 probability for the per-state transition $P_N$, not the per-domain
185 Version 0.10 added the [[-F]] option to sawsim, setting a limit on the
186 force change $dF$ in a single timestep for continuous pulling
187 simulations. It also added the piston tension model (Section
188 \ref{sec.piston_tension}), and adjusted the stiffness calculations to
189 handle domains with undefined stiffness.
191 <<version definition>>=
192 #define VERSION "0.10"
198 sawsim - simulating single molecule mechanical unfolding.
199 Copyright (C) 2008-2010, William Trevor King
201 This program is free software; you can redistribute it and/or
202 modify it under the terms of the GNU General Public License as
203 published by the Free Software Foundation; either version 3 of the
204 License, or (at your option) any later version.
206 This program is distributed in the hope that it will be useful, but
207 WITHOUT ANY WARRANTY; without even the implied warranty of
208 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
209 See the GNU General Public License for more details.
211 You should have received a copy of the GNU General Public License
212 along with this program. If not, see <http://www.gnu.org/licenses/>.
214 The author may be contacted at <wking@drexel.edu> on the Internet, or
215 write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
216 Philadelphia PA 19104, USA.
229 The unfolding system is stored as a chain of \emph{domains} (Figure
230 \ref{fig.domain_chain}). Domains can be folded, globular protein
231 domains, unfolded protein linkers, AFM cantilevers, or other
232 stretchable system components. The system is modeled as graph of
233 possible domain states with transitions between them (Figure
234 \ref{fig.domain_states}).
238 \subfloat[][]{\label{fig.domain_chain}
239 \begin{tikzpicture}[->,node distance=2cm,font=\footnotesize]
240 \tikzstyle{every state}=[fill,draw=red!50,very thick,fill=red!20]
241 \node[state] (A) {domain 1};
242 \node[state] (B) [below of=A] {domain 2};
243 \node[state] (C) [below of=B] {.~.~.};
244 \node[state] (D) [below of=C] {domain $N$};
245 \node (S) [below of=D] {Surface};
246 \node (E) [above of=A] {};
248 \path[-] (A) edge (B)
249 (B) edge node [right] {Tension} (C)
252 \path[->,green] (A) edge node [right,black] {Extension} (E);
256 \subfloat[][]{\label{fig.domain_states}
257 \begin{tikzpicture}[->,node distance=2.5cm,shorten <=1pt,shorten >=1pt,font=\footnotesize]
258 \tikzstyle{every state}=[fill,draw=blue!50,very thick,fill=blue!20]
259 \node[state] (A) {cantilever};
260 \node[state] (C) [below of=A] {transition};
261 \node[state] (B) [left of=C] {folded};
262 \node[state] (D) [right of=C] {unfolded};
264 \path (B) edge [bend left] node [above] {$k_1$} (C)
265 (C) edge [bend left] node [below] {$k_1'$} (B)
266 edge [bend left] node [above] {$k_2$} (D)
267 (D) edge [bend left] node [below] {$k_2'$} (C);
270 \caption{\subref{fig.domain_chain} Extending a chain of domains. One end
271 of the chain is fixed, while the other is extended at a constant
272 speed. The domains are coupled with rigid linkers, so the domains
273 themselves must stretch to accomodate the extension.
274 \subref{fig.domain_states} Each domain exists in a discrete state. At
275 each timestep, it may transition into another state following a
276 user-defined state matrix such as this one, showing a metastable
277 transition state and an explicit ``cantilever'' domain.}
281 Each domain contributes to the total tension, which depends on the
282 chain extension. The particular model for the tension as a function
283 of extension is handled generally with function pointers. So far the
284 following models have been implemented:
286 \item Null (Appendix \ref{sec.null_tension}),
287 \item Constant (Appendix \ref{sec.const_tension}),
288 \item Hookean spring (Appendix \ref{sec.hooke}),
289 \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
290 \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
291 \item Piston (Appendix \ref{sec.piston_tension}),
294 The tension model and parameters used for a particular domain depend
295 on the domain's current state. The transition rates between states
296 are also handled generally with function pointers, with
299 \item Null (Appendix \ref{sec.null_k}),
300 \item Bell's model (Appendix \ref{sec.bell}),
301 \item Spring constant corrected Bell model (Appendix \ref{sec.kbell}),
302 \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
303 \item Kramers' saddle point model (Appendix \ref{sec.kramers}).
306 The unfolding of the chain domains is modeled in two stages. First
307 the tension of the chain is calculated. Then the tension (assumed to
308 be constant over the length of the chain, see Section
309 \ref{sec.timescales}), is applied to each domain, and used to
310 calculate the probability of that domain changing state in the next
311 timestep $dt$. Then the time is stepped forward, and the process is
312 repeated. The simulation is complete when a pre-set time limit
313 [[t_max]] is reached or a pre-selected state [[stop_state]] is empty.
315 int main(int argc, char **argv)
317 <<tension model getopt array definition>>
318 <<k model getopt array definition>>
320 /* variables initialized in get_options() */
321 list_t *state_list=NULL, *transition_list=NULL;
322 environment_t env={0};
323 double P, dF, t_max, dt_max, v, x_step;
324 state_t *stop_state=NULL;
326 /* variables used in the simulation loop */
327 double t=0, x=0, dt, F; /* time, distance, timestep, force */
328 int transition=1; /* boolean marking a transition for tension and timestep optimization */
330 get_options(argc, argv, &P, &dF, &t_max, &dt_max, &v, &x_step,
331 &stop_state, &env, NUM_TENSION_MODELS, tension_models,
332 NUM_K_MODELS, k_models, &state_list, &transition_list);
335 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Distance (m)\tForce (N)\tStiffness (N/m)\n");
336 if (flags & FLAG_OUTPUT_TRANSITION_FORCES) printf("#Force (N)\tInitial state\tFinal state\n");
337 while ((t_max > 0 && t <= t_max) || (stop_state != NULL && stop_state->num_domains > 0)) {
338 F = find_tension(state_list, &env, x, transition, 0);
340 dt = determine_dt(state_list, transition_list, &env, F, x, dF, P, dt_max, v, transition);
342 dt = determine_dt(state_list, transition_list, &env, F, x, dF, P, dt_max, 0, transition);
343 transition=random_transitions(transition_list, F, dt, &env);
344 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness);
349 if (dt == dt_max) { /* step completed */
352 } else { /* still working on this step */
357 destroy_state_list(state_list);
358 destroy_transition_list(transition_list);
362 @ The meat of the simulation is bundled into the three functions
363 [[find_tension]], [[determine_dt]], and [[random_transitions]].
364 [[find_tension]] is discussed in Section \ref{sec.find_tension},
365 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
366 [[random_transitions]] in Section \ref{sec.transition_rate}.
368 The stretched distance is extended in one of two modes depending on
369 [[xstep]]. If [[xstep == 0]], then the pulling is continuous, at
370 least within the limits of the inherent discretization of a time
371 stepped approach. At any rate, the timestep is limited by the maximum
372 allowed timestep [[dt_max]] and the maximum allowed unfolding
373 probability in a given step [[P]]. In the second mode [[xstep > 0]],
374 and the end to end distance increases in discrete steps of that
375 length. The time between steps is chosen to maintain an average
376 unfolding velocity [[v]]. This approach is less work to simulate than
377 smooth pulling and also closer to the reality of many experiments, but
378 it is practically impossible to treat analytically. With both pulling
379 styles implemented, the effects of the discretization can be easily
380 calculated, bridging the gap between theory and experiment.
382 Environmental parameters important in determining reaction rates and
383 tensions (e.g.\ temperature, pH) are stored in a single structure to
384 facilitate extension to more complicated models in the future.
385 <<environment definition>>=
386 typedef struct environment_struct {
396 #define DOUBLE_PRECISION 1e-12
398 <<environment definition>>
399 <<one dimensional function definition>>
400 <<create/destroy definitions>>
402 #endif /* GLOBAL_H */
404 Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
405 included multiple times.
407 \section{Simulation functions}
409 <<simulation functions>>=
413 <<does domain transition>>
414 <<randomly transition domains>>
418 \label{sec.find_tension}
420 Because the stretched system may be made up of several parts (folded
421 domains, unfolded domains, spring-like cantilever, \ldots), we will
422 assign the domains to states with tension models and groups. The
423 models are used to determine the tension of all the domain in a given
424 state following some model (Hookean spring, WLC, \ldots). The domains
425 are assumed to be commutative, so ordering is ignored. The
426 interactions between the groups are assumed to be linear, but the
427 interactions between domains of the same group need not be. Each
428 model has a tension handler function, which gives the tension $F$
429 needed to stretch a given group of domains a distance $x$.
431 Groups represent collections of states which the model should treat as
432 a single entity. To understand the purpose of groups, consider a
433 system with two unfolded domain states (e.g.\ for two different
434 proteins) that were both modeled as WLCs. If both unfolded states had
435 the same persistence length, it would be wise to assign them to the
436 same group. This would reduce both the computational cost of
437 balancing the tension and the error associated with the inter-group
438 interaction linearity. Note that group numbers only matter
439 \emph{within} model classes, since grouping domains with seperate
440 models doesn't make sense. Within designated groups, the tension
441 parameters for each domain are still checked for compatibility, so if
442 you accidentally assign incompatible domains to the same group you
443 should get an appropriate error message.
445 <<find tension definitions>>=
446 #define NUM_TENSION_MODELS 6
450 <<tension handler array global declaration>>=
451 extern one_dim_fn_t *tension_handlers[];
454 <<tension handler array global>>=
455 one_dim_fn_t *tension_handlers[] = {
457 &const_tension_handler,
465 <<tension model global declarations>>=
466 <<tension handler array global declaration>>
469 <<tension handler types>>=
470 typedef struct tension_handler_data_struct {
471 list_t *group_tension_params; /* one item for each domain in the group */
474 } tension_handler_data_t;
478 After sorting the chain into separate groups $G_i$, with tension
479 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
480 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
481 \forall i,j$. Note that there may be several states within each
482 group. [[find_tension]] is basically a wrapper to feed proper input
483 into the [[tension_balance]] function. [[transition]] should be set
484 to 0 if there were no transitions since the previous call to
485 [[find_tension]] to support some optimizations that assume a small
486 increase in tension between steps (only true if no transition has
487 occured). While we're messing with the tension handlers and if
488 [[const_env==0]], we also grab a measure of the current linker
489 stiffness for the environmental variable, which is needed by some of
490 the unfolding rate models (e.g.\ the linker-stiffness corrected Bell
491 model, Appendix \ref{sec.kbell}).
493 double find_tension(list_t *state_list, environment_t *env, double x,
494 int transition, int const_env)
495 { // TODO: !cont_env -> step_call, only alter env or statics if step_call==1
496 static int num_active_groups=0;
497 static one_dim_fn_t **per_group_handlers = NULL;
498 static one_dim_fn_t **per_group_inverse_handlers = NULL;
499 static void **per_group_data = NULL; /* array of pointers to tension_handler_data_t */
500 static double last_x;
501 tension_handler_data_t *tdata;
502 double F, *pStiffness;
506 fprintf(stderr, "%s:\tfinding tension at distance %g%s\n", __FUNCTION__, x, transition ? " after a transition" : "");
509 if (transition != 0 || num_active_groups == 0) { /* setup statics */
510 /* get new statics, freeing the old ones if they aren't NULL */
511 get_active_groups(state_list, &num_active_groups, &per_group_handlers, &per_group_inverse_handlers, &per_group_data);
513 /* fill in the group handlers and params */
515 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]);
517 for (i=0; i<num_active_groups; i++) {
518 tdata = (tension_handler_data_t *) per_group_data[i];
520 fprintf(stderr, "%s:\tactive group %d of %d (data %p, param list %p)", __FUNCTION__, i, num_active_groups, tdata, tdata->group_tension_params);
521 list_print(stderr, tdata->group_tension_params, " parameters");
524 /* tdata->persist continues unchanged */
527 } /* else continue with the current statics */
532 pStiffness = &env->stiffness;
534 F = tension_balance(num_active_groups, per_group_handlers, per_group_inverse_handlers, (void **)per_group_data, last_x, x, pStiffness);
540 @ For the moment, we will restrict our group boundaries to a single
541 dimension, so $\sum_i x_i = x$, our total extension (it is this
542 restriction that causes the groups to interact linearly). We'll also
543 restrict our extensions to all be positive. With these restrictions,
544 the problem of balancing the tensions reduces to solving a system of
545 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
546 the number of active groups. In general this can be fairly
547 complicated, but without much loss of practicality we can restrict
548 ourselves to strictly monotonically increasing, non-negative tension
549 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
550 simpler. For example, it guarantees the existence of a unique, real
551 solution for finite forces. See Appendix \ref{sec.tension_balance}
552 for the tension balancing implementation.
554 Each group must define a parameter structure if the tension model is
556 a function to create the parameter structure,
557 a function to destroy the parameter structure, and
558 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as its data argument.
559 The tension-specific parameter structures are contained in the domain
560 group field. For optimization, a group may also define an inverse
561 tension function as an optimization. See the Section
562 \ref{sec.model_selection} for a discussion on the command line
563 framework and inverse function discussion. See the worm-like chain
564 implementation for an example (Section \ref{sec.wlc}).
566 The major limitation of our approach is that all of our tension models
567 are Markovian (memory-less), which is not really a problem for the
568 chains (see \citet{evans99} for WLC thermalization rate calculations).
570 \subsection{Transition rate}
571 \label{sec.transition_rate}
573 Each state transition is modeled as unimolecular, first order reaction
575 \text{State 1} \xrightarrow{k} \text{State 2}
577 With the rate constant $k$ defined by
579 \frac{d[\text{State 2}]}{dt} = -\frac{d[\text{State 1}]}{dt} = k [\text{State 1}].
581 So the probability of a given domain transitioning in a short time
587 For a group of $N$ identical domains, and therefore identical
588 unfolding probabilities $P_1$, the probability of $n$ domains
589 unfolding in a given timestep is given by the binomial distribution
591 P(n) = \frac{N!}{n!(N-n)!}P_1^n(1-P_1)^{(N-n)}.
593 The probability that \emph{none} of these domains unfold is then
597 so the probability that \emph{at least one} domain unfolds is
601 Note that for small $P$,
603 p(n>0) = 1-(1-NP_1+\mathcal{O}(P_1^2)) = NP_1 - \mathcal{O}(P^2)
606 This conversion from $P_1$ to $P_N$ is handled by the [[pN]] macro
608 /* find multi-domain transition rate from the single domain rate */
609 #define PN(P,N) (1.0-pow(1.0-(P), (N)))
613 We can also define a macro to find the approprate $dt$ to achieve a
614 target $P_N$ with a given $k$ and $N$.
616 P_N &= 1-(1-P_1)^N \\
617 1-P_1 &= (1-P_N)^{1/N} \\
618 P_1 &= 1-(1-P_N)^{1/N}
621 #define P1(PN,N) (1.0-pow(1.0-(PN), 1.0/(N)))
624 We take some time to discuss the meaning of $p(n>1)$
625 (i.e. multi-unfolding timesteps) in Section \ref{sec.timescales}.
627 <<does domain transition>>=
628 int domain_transitions(double F, double dt, environment_t *env, transition_t *transition)
629 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to transition structure */
631 k = accel_k(transition->k, F, env, transition->k_params);
632 //(*transition->k)(F, env, domain->k_params);
633 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
634 return happens(PN(k*dt, N_INIT(transition)));
636 @ [[happens]] is a random decision making function defined in Appendix \ref{sec.utils}.
638 <<randomly transition domains>>=
639 int random_transitions(list_t *transition_list, double tension,
640 double dt, environment_t *env)
641 { /* returns 1 if there was a transition and 0 otherwise */
642 while (transition_list != NULL) {
643 if (T(transition_list)->initial_state->num_domains > 0
644 && domain_transitions(tension, dt, env, T(transition_list))) {
645 if (flags & FLAG_OUTPUT_TRANSITION_FORCES)
646 fprintf(stdout, "%g\t%s\t%s\n", tension, T(transition_list)->initial_state->name, T(transition_list)->final_state->name);
647 T(transition_list)->initial_state->num_domains--;
648 T(transition_list)->final_state->num_domains++;
649 return 1; /* our one domain has transitioned, stop looking for others */
651 transition_list = transition_list->next;
657 \subsection{Timescales}
658 \label{sec.timescales}
660 The simulation assumes that chain equilibration is instantaneous
661 relative to the simulation timestep $dt$, so that tension is uniform
662 along the chain. The quality of this assumption depends on your
663 particular chain. For example, a damped spring thermalizes on a
664 timescale of order $\tau = 1/\gamma$ where the damping ratio $\gamma
665 \equiv \omega_0(-\zeta \pm \sqrt{\zeta^2-1}$, the natural angular
666 frequency $omega_0 \equiv \sqrt{k/m}$, $\zeta \equiv b/2\sqrt{km}$,
667 and $b$ sets the drag $F_b = -bv$. For our cantilevers $\tau$ is
668 on the order of milliseconds, which is longer than a timestep.
669 However, the approximation is still reasonable, because there is
670 rarely unfolding during the cantilever return. You could, of course,
671 take cantilever, WLC, etc.\ response times into effect, but that
672 would significantly complicate a simulation that seems to work
673 well enough without bothering :p. For WLC and FJC relaxation times,
674 see the Appendix of \citet{evans99} and \citet{kroy07}.
676 Besides assuming our timestep is much greater than equilibration
677 times, we also force it to be short enough so that only one domain may
678 unfold in a given timestep. For an unfolding timescale of $dt_u =
679 1/k$, we adapt our timesteps to keep $dt \ll dt_u$, so the probability
680 of two domains unfolding in the same timestep is small (see
681 [[determine_dt]] in Section \ref{sec.adaptive_dt} and [[P]] in
682 [[main()]] in Section \ref{sec.main} set by [[-P]] in
683 [[get_options()]] in Section \ref{sec.get_options}). This approach
684 breaks down as the adaptive timestep scheme approaches the transition
685 timestep $dt_t$, but $dt_t \sim 1\U{$\mu$s}$ for Ig27-like proteins
686 \citep{klimov00}, so this shouldn't be a problem. To reassure
687 yourself, you can ask the simulator to print the smallest timestep
688 that was used with TODO.
690 Even if the likelyhood of two domains unfolding in the same timestep
691 is small, if you run long enough simulations it will eventually occur.
692 In this case, we assume that $dt_t \ll dt$, so even if two domains
693 would unfold if we stuck to our unfolding rate $k$ for an entire
694 timestep $dt$, in ``reality'' only one of those domains unfolds.
695 Because we do not know when in the timestep the unfolding took place,
696 we move on to the next timestep without checking for further
697 unfoldings. This ``unchecked timestep portion'' should not contribute
698 significant errors to the calculation, because if $dt$ was small
699 enough that unfolding was unlikely at the pre-unfolding force, it
700 should be even more unlikely at the post-unfolding force, especially
701 over only a fraction of the timestep. The only dangerous cases for
702 the current implementation are those where unfolding intermediates are
703 much less stable than their precursor states, in which case an
704 unfolding event during the remaining timestep may actually be likely.
705 A simple workaround for such cases is to pick the value for [[P]] to
706 be small enough that the $dt \ll dt_u$ assumption survives even under
707 a transition populating the unstable intermediate.
709 \subsection{Adaptive timesteps}
710 \label{sec.adaptive_dt}
712 We'd like to pick $dt$ so the probability of changing state
713 (e.g.\ unfolding) in the next timestep is small. If we don't adapt our
714 timestep, we also risk breaking our approximation that $F(x' \in [x,
715 x+v \cdot dt]) = F(x)$ or that only one domain changes state in a
716 given timestep. The simulation would have been accurate for
717 sufficiently small fixed timesteps, but adaptive timesteps will allow
718 us to move through low tension regions in fewer steps, leading to a
719 more efficient simulation.
721 Consider the two state folded $\rightarrow$ unfolded transition.
722 Because $F(x)$ is monotonically increasing between unfoldings,
723 excessively large timesteps will lead to erroneously small unfolding
724 rates (an thus, higher average unfolding force).
726 The actual adaptive timestep implementation is not particularly
727 interesting, since we are only required to reduce $dt$ to somewhere
728 below a set threshold, so I've removed it to Appendix
729 \ref{sec.adaptive_dt_implementation}.
735 The models provide the physics of the simulation, but the simulation
736 allows interchangeable models, and we are currently focusing on the
737 simulation itself, so we remove the actual model implementation to the
740 The tension models are defined in Appendix \ref{sec.tension_models},
741 and unfolding models are defined in Appendix \ref{sec.k_models}.
744 #define k_B 1.3806503e-23 /* J/K */
749 \section{Command line interface}
751 <<get option functions>>=
753 <<index tension model>>
758 \subsection{Model selection}
759 \label{sec.model_selection}
761 The main difficulty with the command line interface in \prog\ is
762 developing an intuitive method for accessing the possibly large number
763 of available models. We'll treat this generally by defining an array
764 of available models, containing their important parameters.
766 <<get options definitions>>=
767 <<model getopt definitions>>
770 <<create/destroy definitions>>=
771 typedef void *create_data_func_t(char **param_strings);
772 typedef void destroy_data_func_t(void *param_struct);
775 <<model getopt definitions>>=
776 <<tension model getopt definitions>>
777 <<k model getopt definitions>>
780 In order to choose models, we need a method of assembling a reaction
781 graph such as that in Figure \ref{fig.domain_states} and population
782 the graph with a set of domains. First we declare the domain states
785 [[-s <state-name>,<model>[,<group>] -S <model-params> [-N <initial-population]]] \\
789 [[-s unfolded,wlc,1 -S 1e-8,4e-10 -N 5]]
791 This creates the state named [[unfolded]], using the WLC model and one
792 as a group number (defaults to 0, see Section \ref{sec.find_tension}).
793 The tension parameters are then set to [[1e-8,4e-10]], which in the
794 case of the WLC are the contour and persistence lengths respectively.
795 Finally we populate the state by adding five domains to the state.
796 The name of the state is importent for identifying states later.
798 Once the domains have been created and populated, you can added
799 transition rates following
801 [[-k <initial-state>,<final-state>,<model> -K <model-params>]] \\
805 [[-k folded,unfolded,bell -K 3.3e-4,0.25e-9]]
807 This creates a reaction going from the [[folded]] state to the
808 [[unfolded]] state, with the rate constant given by the Bell model
809 (Appendix \ref{sec.bell}). The unfolding parameters are then set to
810 [[3.3e-4,0.25e-9]], which in the case of the Bell model are the
811 unforced rate and transition state distance respectively.
813 \subsubsection{Tension}
815 To access the various tension models from the command line, we define
816 a structure that contains all the neccessary information about the
818 <<tension model getopt definitions>>=
819 typedef struct tension_model_getopt_struct {
820 one_dim_fn_t *handler; /* fn implementing the model on a group */
821 one_dim_fn_t *inverse_handler; /* fn implementing inverse on group */
822 char *name; /* name string identifying the model */
823 char *description; /* a brief description of the model */
824 int num_params; /* number of extra params for the model */
825 char **param_descriptions; /* descriptions of extra params */
826 char *params; /* default values for the extra params */
827 create_data_func_t *creator; /* param-string -> model-data-struct fn */
828 destroy_data_func_t *destructor; /* fn to free the model data structure */
829 } tension_model_getopt_t;
830 @ Set [[inverse_handler]] to [[NULL]] if you wish to invert using
831 bisection. Obviously, this will be slower than computing an
832 analytical inverse and more prone to rounding errors, but it may be
833 easier if a closed-form inverse does not exist.
835 <<tension model getopt array definition>>=
836 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
837 <<null tension model getopt>>,
838 <<constant tension model getopt>>,
839 <<hooke tension model getopt>>,
840 <<worm-like chain tension model getopt>>,
841 <<freely-jointed chain tension model getopt>>,
842 <<piston tension model getopt>>
846 \subsubsection{Unfolding rate}
848 <<k model getopt definitions>>=
849 #define NUM_K_MODELS 6
851 typedef struct k_model_getopt_struct {
856 char **param_descriptions;
858 create_data_func_t *creator;
859 destroy_data_func_t *destructor;
863 <<k model getopt array definition>>=
864 k_model_getopt_t k_models[NUM_K_MODELS] = {
865 <<null k model getopt>>,
866 <<const k model getopt>>,
867 <<bell k model getopt>>,
868 <<kbell k model getopt>>,
869 <<kramers k model getopt>>,
870 <<kramers integ k model getopt>>
877 void help(char *prog_name, double P, double dF,
878 double t_max, double dt_max, double v, double x_step,
879 state_t *stop_state, environment_t *env,
880 int n_tension_models, tension_model_getopt_t *tension_models,
881 int n_k_models, k_model_getopt_t *k_models,
882 int k_model, list_t *state_list)
887 if (state_list != NULL)
888 state = S(tail(state_list));
890 printf("usage: %s [options]\n", prog_name);
891 printf("Version %s\n\n", VERSION);
892 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
894 printf("Simulation options:\n");
896 printf("-q state\tStop simulation when this state empties (- for no limit, currently %s)\n", stop_state ? stop_state->name : "-" );
897 printf("-t tmax\tMaximum simulated time (-1 for no limit, currently %g s)\n", t_max);
898 printf("-d dt\tMaximum allowed timestep dt (currently %g s)\n", dt_max);
899 printf("-P P\tMaximum transition probability for dt (currently %g)\n", P);
900 printf("-F dF\tMaximum change in force over dt. Only relevant if dx > 0. (currently %g N)\n", dF);
901 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
902 printf("-x dx\tPulling step size (currently %g m)\n", x_step);
903 printf("\tWhen dx = 0, the pulling is continuous.\n");
904 printf("\tWhen dx > 0, the pulling occurs in discrete steps.\n");
906 printf("Environmental options:\n");
907 printf("-T T\tTemperature T (currently %g K)\n", env->T);
908 printf("-C T\tYou can also set the temperature T in Celsius\n");
910 printf("State creation options:\n");
911 printf("-s name,model[[,group],params]\tCreate a new state.\n");
912 printf("Once you've set up the state, you may populate it with domains\n");
913 printf("-N n\tAdd n domains in the current state (%s)\n", state ? state->name : "-");
915 printf("Transition creation options:\n");
916 printf("-k initial_state,final_state,model[,params]\tCreate a new transition\n");
918 printf("To double check your reaction graph, you can produce graphviz dot output\n");
919 printf("describing the current states and transitions.\n");
920 printf("-D\tPrint dot description of states and transitions and exit.\n");
922 printf("Simulation output mode options:\n");
923 printf("There are two output modes. In standard mode, only the transition\n");
924 printf("events are printed. For example:\n");
925 printf(" #Force (N)\tInitial state\tFinal state\n");
926 printf(" 123.456e-12\tfolded\tunfolded\n");
928 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
929 printf(" #Distance (m)\tForce (N)\tStiffness (N/m)\n");
930 printf(" 0.001\t0.0023\t0.002\n");
932 printf("-V\tChange output to verbose mode.\n");
933 printf("-h\tPrint this help and exit.\n");
936 printf("Tension models:\n");
937 for (i=0; i<n_tension_models; i++) {
938 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
939 for (j=0; j < tension_models[i].num_params ; j++)
940 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
941 printf("\t\tdefaults: %s\n", tension_models[i].params);
943 printf("Unfolding rate models:\n");
944 for (i=0; i<n_k_models; i++) {
945 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
946 for (j=0; j < k_models[i].num_params ; j++)
947 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
948 printf("\t\tdefaults: %s\n", k_models[i].params);
951 printf("Examples:\n");
952 printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded. Stop once there are no more folded domains.\n");
953 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);
954 printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded. Stop after 0.25 seconds\n");
955 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);
959 \subsection{Get options}
960 \label{sec.get_options}
964 void get_options(int argc, char **argv, double *pP, double *pDF,
965 double *pT_max, double *pDt_max, double *pV,
966 double *pX_step, state_t **pStop_state, environment_t *env,
967 int n_tension_models, tension_model_getopt_t *tension_models,
968 int n_k_models, k_model_getopt_t *k_models,
969 list_t **pState_list, list_t **pTransition_list)
971 char *prog_name = NULL;
972 char c, options[] = "q:t:d:P:F:v:x:T:C:s:N:k:DVh";
973 int tension_model=0, k_model=0, initial_state=-1, final_state=-1;
974 char *state_name=NULL;
975 int i, n, tension_group=0;
976 list_t *temp_list=NULL;
979 void *model_params; /* for INIT_MODEL() */
980 int num_param_args; /* for INIT_MODEL() */
981 char **param_args; /* for INIT_MODEL() */
983 extern int optind, optopt, opterr;
988 for (i=0; i<n_tension_models; i++) {
989 fprintf(stderr, "tension model %i %s:\t model's handler %p\t inverse handler %p\n",
990 i, tension_models[i].name, tension_models[i].handler, tension_models[i].inverse_handler);
991 assert(tension_models[i].handler == tension_handlers[i]);
996 flags = FLAG_OUTPUT_TRANSITION_FORCES;
999 *pT_max = -1; /* s, -1 removes this limit */
1000 *pDt_max = 0.001; /* s */
1001 *pP = 1e-3; /* % pop per s */
1002 *pDF = 1e-12; /* N */
1003 *pV = 1e-6; /* m/s */
1004 *pX_step = 0.0; /* m */
1005 env->T = 300.0; /* K */
1006 *pState_list = NULL;
1007 *pTransition_list = NULL;
1009 while ((c=getopt(argc, argv, options)) != -1) {
1012 if (STRMATCH(optarg, "-")) {
1013 *pStop_state = NULL;
1015 temp_list = head(*pState_list);
1016 while (temp_list != NULL) {
1017 if (STRMATCH(S(temp_list)->name, optarg)) {
1018 *pStop_state = S(temp_list);
1021 temp_list = temp_list->next;
1025 case 't': *pT_max = safe_strtod(optarg, "-t"); break;
1026 case 'd': *pDt_max = safe_strtod(optarg, "-d"); break;
1027 case 'P': *pP = safe_strtod(optarg, "-P"); break;
1028 case 'F': *pDF = safe_strtod(optarg, "-F"); break;
1029 case 'v': *pV = safe_strtod(optarg, "-v"); break;
1030 case 'x': *pX_step = safe_strtod(optarg, "-x"); break;
1031 case 'T': env->T = safe_strtod(optarg, "-T"); break;
1032 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
1034 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
1035 assert(num_strings >= 2 && num_strings <= 4);
1037 state_name = strings[0];
1038 if (state_index(*pState_list, state_name) != -1) {
1039 fprintf(stderr, "Invalid option '%s', repeated state name '%s'\n", optarg, state_name);
1042 tension_model = index_tension_model(n_tension_models, tension_models, strings[1]);
1043 if (num_strings == 4)
1044 tension_group = safe_strtoi(strings[2], optarg);
1047 if (num_strings >= 3)
1048 tension_models[tension_model].params = strings[num_strings-1];
1050 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);
1052 INIT_MODEL(state_name, &tension_models[tension_model], tension_models[tension_model].params, model_params);
1053 push(pState_list, create_state(state_name,
1054 tension_models[tension_model].handler,
1055 tension_models[tension_model].inverse_handler,
1058 tension_models[tension_model].destructor,
1060 *pState_list = tail(*pState_list);
1062 free_parsed_list(num_strings, strings);
1065 n = safe_strtoi(optarg, "-N");
1067 fprintf(stderr, "%s:\tadding %d domains to %s\n", __FUNCTION__, n, S(*pState_list)->name);
1069 S(*pState_list)->num_domains += n;
1072 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
1073 assert(num_strings >= 3 && num_strings <= 4);
1074 initial_state = state_index(*pState_list, strings[0]);
1075 final_state = state_index(*pState_list, strings[1]);
1076 k_model = index_k_model(n_k_models, k_models, strings[2]);
1078 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);
1080 assert(initial_state != final_state);
1081 if (num_strings == 4)
1082 k_models[k_model].params = strings[3];
1083 INIT_MODEL(state_name, &k_models[k_model], k_models[k_model].params, model_params);
1084 push(pTransition_list, create_transition(list_element(*pState_list, initial_state),
1085 list_element(*pState_list, final_state),
1086 k_models[k_model].k,
1088 k_models[k_model].destructor), 1);
1089 free_parsed_list(num_strings, strings);
1092 print_state_graph(stdout, *pState_list, *pTransition_list);
1096 flags = FLAG_OUTPUT_FULL_CURVE;
1099 fprintf(stderr, "unrecognized option '%c'\n", optopt);
1100 /* fall through to default case */
1102 help(prog_name, *pP, *pDF, *pT_max, *pDt_max, *pV, *pX_step,
1103 *pStop_state, env, n_tension_models, tension_models,
1104 n_k_models, k_models, k_model, *pState_list);
1108 /* check the arguments */
1109 assert(*pState_list != NULL); /* no domains! */
1110 assert(*pP > 0.0); assert(*pP < 1.0);
1111 assert(*pDt_max > 0.0);
1113 assert(env->T > 0.0);
1115 crosslink_groups(*pState_list, *pTransition_list);
1121 <<index tension model>>=
1122 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
1125 for (i=0; i<n_models; i++)
1126 if (STRMATCH(models[i].name, name))
1128 if (i == n_models) {
1129 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
1137 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
1140 for (i=0; i<n_models; i++)
1141 if (STRMATCH(models[i].name, name))
1143 if (i == n_models) {
1144 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
1151 <<initialize model definition>>=
1152 /* requires int num_param_args and char **param_args in the current scope
1154 * INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
1155 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
1157 #define INIT_MODEL(role, model, param_string, param_pointer) \
1159 parse_list_string((param_string), SEP, '{', '}', \
1160 &num_param_args, ¶m_args); \
1161 if (num_param_args != (model)->num_params) { \
1163 "%s model %s expected %d params,\n", \
1164 (role), (model)->name, (model)->num_params); \
1166 "not the %d params in '%s'\n", \
1167 num_param_args, (param_string)); \
1168 assert(num_param_args == (model)->num_params); \
1170 if ((model)->creator) \
1171 param_pointer = (*(model)->creator)(param_args); \
1173 param_pointer = NULL; \
1174 free_parsed_list(num_param_args, param_args); \
1178 First we define some safe string parsing functions. Currently these
1179 abort on any irregularity, but in the future they could print error
1180 messages, etc. [[static]] because the functions are currently
1181 defined in each [[*.c]] file that uses them, but perhaps they should
1182 be moved to [[utils.h/c]] or some such instead\ldots
1185 static int safe_strtoi(const char *s, const char *description) {
1189 ret = strtol(s, &endp, 10);
1190 if (endp[0] != '\0') { //strlen(endp) > 0) {
1191 fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1192 endp, s, description);
1193 assert(1==0); //strlen(endp) == 0);
1195 fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1198 } else if (errno == ERANGE) {
1199 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1200 assert(errno != ERANGE);
1205 static double safe_strtod(const char *s, const char *description) {
1209 ret = strtod(s, &endp);
1210 if (endp[0] != '\0') { //strlen(endp) > 0) {
1211 fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1212 endp, s, description);
1213 assert(1==0); //strlen(endp) == 0);
1215 fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1218 } else if (errno == ERANGE) {
1219 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1220 assert(errno != ERANGE);
1229 \addcontentsline{toc}{section}{Appendicies}
1231 \section{sawsim.c details}
1235 The general layout of our simulation code is:
1246 We include [[math.h]], so don't forget to link to the libm with `-lm'.
1248 #include <assert.h> /* assert() */
1249 #include <stdlib.h> /* malloc(), free(), rand(), strto*() */
1250 #include <stdio.h> /* fprintf(), stdout */
1251 #include <string.h> /* strlen, strtok() */
1252 #include <math.h> /* exp(), M_PI, sqrt() */
1253 #include <sys/types.h> /* pid_t (returned by getpid()) */
1254 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
1255 #include <errno.h> /* errno, ERANGE (for safe_strto*()) */
1258 #include "tension_balance.h"
1259 #include "k_model.h"
1260 #include "tension_model.h"
1262 #include "accel_k.h"
1267 <<version definition>>
1268 <<flag definitions>>
1269 <<max/min definitions>>
1270 <<string comparison definition and globals>>
1271 <<initialize model definition>>
1272 <<get options definitions>>
1273 <<state definitions>>
1274 <<transition definitions>>
1283 <<create/destroy state>>
1284 <<destroy state list>>
1286 <<create/destroy transition>>
1287 <<destroy transition list>>
1288 <<print state graph>>
1290 <<simulation functions>>
1291 <<boilerplate functions>>
1294 <<boilerplate functions>>=
1296 <<get option functions>>
1302 srand(getpid()*time(NULL)); /* seed rand() */
1306 <<flag definitions>>=
1307 /* in octal b/c of prefixed '0' */
1308 #define FLAG_OUTPUT_FULL_CURVE 01
1309 #define FLAG_OUTPUT_TRANSITION_FORCES 02
1313 static unsigned long int flags = 0;
1316 \subsection{Utilities}
1319 <<max/min definitions>>=
1320 #define MAX(a,b) ((a)>(b) ? (a) : (b))
1321 #define MIN(a,b) ((a)<(b) ? (a) : (b))
1324 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
1325 <<string comparison definition and globals>>=
1326 // Check if two strings match, return 1 if they do
1327 static char *temp_string_A;
1328 static char *temp_string_B;
1329 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
1330 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
1331 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
1332 /* +1 to also compare the '\0' */
1335 We also define a macro for our [[check]] unit testing
1336 <<check relative error>>=
1337 #define CHECK_ERR(max_err, expected, received) \
1339 fail_unless( (received-expected)/expected < max_err, \
1340 "relative error %g >= %g in %s (Expected %g, received %g)", \
1341 (received-expected)/expected, max_err, #received, \
1342 expected, received); \
1343 fail_unless(-(received-expected)/expected < max_err, \
1344 "relative error %g >= %g in %s (Expected %g, received %g)", \
1345 -(received-expected)/expected, max_err, #received, \
1346 expected, received); \
1352 int happens(double probability)
1354 assert(probability >= 0.0); assert(probability <= 1.0);
1355 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*/
1360 \subsection{Adaptive timesteps}
1361 \label{sec.adaptive_dt_implementation}
1363 $F(x)$ increases with $x$, possibly exploding (e.g.\ in the worm-like
1364 chain model), so basing the timestep on the unfolding probability at
1365 the current tension is dangerous, and we need to search for a $dt$ for
1366 which $P_N(F(x+v*dt))<P_\text{max}$ (See Section
1367 \ref{sec.transition_rate} for a discussion of $P_N$). There are two
1368 cases to consider. In the most common, no domains have unfolded since
1369 the last step, and we expect the next step to be slightly shorter than
1370 the previous one. In the less common, domains did unfold in the last
1371 step, and we expect the next step to be considerably longer than the
1374 double search_dt(transition_t *transition,
1376 environment_t *env, double x,
1377 double max_prob, double max_dt, double v)
1378 { /* Returns a valid timestep dt in seconds for the given transition.
1379 * Takes a list of all states, a pointer env to the environmental data,
1380 * a starting separation x in m, a max_prob between 0 and 1,
1381 * max_dt in s, stretching velocity v in m/s.
1383 double F, k, P, dtCur, dtU, dtUCur, dtL, dt;
1385 /* get upper bound using the current position */
1386 F = find_tension(state_list, env, x, 0, 1); /* BUG. repeated calculation */
1387 //printf("Start with x = %g (F = %g)\n", x, F);
1388 k = accel_k(transition->k, F, env, transition->k_params);
1389 //printf("x %g\tF %g\tk %g\n", x, F, k);
1390 dtU = P1(max_prob, N_INIT(transition)) / k; /* upper bound on dt */
1392 //printf("overshot max_dt\n");
1395 if (v == 0) /* discrete stepping, so force is temporarily constant */
1398 /* set a lower bound on dt too */
1401 /* The dt determined above may produce illegitimate forces or ks.
1402 * Reduce the upper bound until we have valid ks. */
1404 F = find_tension(state_list, env, x+v*dt, 0, 1);
1405 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1408 F = find_tension(state_list, env, x+v*dt, 0, 1);
1410 //printf("Try for dt = %g (F = %g)\n", dt, F);
1411 k = accel_k(transition->k, F, env, transition->k_params);
1412 /* returned k may be -1.0 */
1413 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1414 while (k == -1.0) { /* reduce step until we hit a valid k */
1416 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1417 F = find_tension(state_list, env, x+v*dt, 0, 1);
1418 //printf("Try for dt = %g (F = %g)\n", dt, F);
1419 k = accel_k(transition->k, F, env, transition->k_params);
1420 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1422 assert(dtU > 1e-14); /* timestep to valid k too small */
1423 /* safe timestep back from x+dtU */
1424 dtUCur = P1(max_prob, N_INIT(transition)) / k;
1426 return dt; /* dtU is safe. */
1428 /* dtU wasn't safe, lets see what would be. */
1429 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1430 dt = (dtU + dtL) / 2.0;
1431 F = find_tension(state_list, env, x+v*dt, 0, 1);
1432 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1433 k = accel_k(transition->k, F, env, transition->k_params);
1434 dtCur = P1(max_prob, N_INIT(transition)) / k;
1435 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1436 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1438 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1439 dtU = dt; /* ... stepping out only dtCur would be safe */
1442 break; /* dtCur = dt */
1444 return MAX(dtUCur, dtL);
1449 Determine $dt$ for an array of potentially different folded domains.
1450 We apply [[search_dt()]] to each populated domain to find the maximum
1451 $dt$ the domain can handle. The final $dt$ is less than the
1452 individual domain $dt$s (to satisfy [[search_dt()]], see above). If
1453 $v>0$ (continuous pulling), $dt$ is also reduced to satisfy
1454 $|F(x+v*dt)-F(x)|<dF_\text{max}$, which limits errors due to our
1455 transition rate assumption that the force does not change appreciably
1456 over a single timestep.
1459 double determine_dt(list_t *state_list, list_t *transition_list,
1460 environment_t *env, double F, double x, double max_dF,
1461 double max_prob, double dt_max, double v, int transition)
1462 { /* Returns the timestep dt in seconds.
1463 * Takes the state and transition lists, a pointer to the environment,
1464 * the force F(x) in N, an extension x in m, max_dF in N,
1465 * max_prob between 0 and 1, dt_max in s, v in m/s, and transition in {0,1}*/
1466 double dt=dt_max, new_dt, F_new;
1467 assert(max_dF > 0.0);
1468 assert(max_prob > 0.0); assert(max_prob < 1.0);
1469 assert(dt_max > 0.0);
1470 assert(state_list != NULL);
1471 assert(transition_list != NULL);
1472 transition_list = head(transition_list);
1475 /* Ensure |dF| < max_dF */
1476 F_new = find_tension(state_list, env, x+v*dt, transition, 1);
1477 while (fabs(F_new - F) > max_dF) {
1479 F_new = find_tension(state_list, env, x+v*dt, transition, 1);
1483 /* Ensure all individual domains pass search_dt() */
1484 while (transition_list != NULL) {
1485 if (T(transition_list)->initial_state->num_domains > 0) {
1486 new_dt = search_dt(T(transition_list),state_list,env,x,max_prob,dt,v);
1487 dt = MIN(dt, new_dt);
1489 transition_list = transition_list->next;
1496 \subsection{State data}
1497 \label{sec.state_data}
1499 The domains exist in one of an array of possible states. Each state
1500 has a tension model which determines it's force/extension curve
1501 (possibly [[null]] for rigid states, see Appendix
1502 \ref{sec.null_tension}).
1504 Groups are, for example, for two domain states with WLC tensions; one
1505 with $p=0.4\U{nm}$, $L=10\U{nm}$ and the other with $p=0.4\U{nm}$,
1506 $L=28\U{nm}$. Since the persistence lengths are the same, you would
1507 like to add the contour lengths of all the domains in \emph{both}
1508 states, and plug that total contour length into a single WLC
1509 evaluation (see Section \ref{sec.find_tension}).
1510 <<state definitions>>=
1511 typedef struct state_struct {
1512 char *name; /* name string */
1513 one_dim_fn_t *tension_handler; /* tension handler */
1514 one_dim_fn_t *inverse_tension_handler; /* inverse tension handler */
1515 int tension_group; /* for grouping similar tension states */
1516 void *tension_params; /* pointer to tension parameters */
1517 destroy_data_func_t *destroy_tension;
1518 int num_domains; /* number of domains currently in state */
1519 /* cached lookups generated from the above data */
1520 int num_out_trans; /* length of out_trans array */
1521 struct transition_struct **out_trans; /* transitions leaving this state */
1522 int num_grouped_states; /* length of grouped states array */
1523 struct state_struct **grouped_states; /* states grouped with this one (doesn't include self) */
1526 /* get the state data for the current list node */
1527 #define S(list) ((state_t *)(list)->d)
1529 @ [[tension_params]] is a pointer to the parameters used by the
1530 handler function when determining the tension. [[destroy_tension]]
1531 points to a function for freeing [[tension_params]]. We it in the
1532 state structure so that [[destroy_state]] and will have an easier job
1535 [[create_]] and [[destroy_state]] are simple wrappers around
1536 [[malloc]] and [[free]].
1537 <<create/destroy state>>=
1538 state_t *create_state(char *name,
1539 one_dim_fn_t *tension_handler,
1540 one_dim_fn_t *inverse_tension_handler,
1542 void *tension_params,
1543 destroy_data_func_t *destroy_tension,
1546 state_t *ret = (state_t *)malloc(sizeof(state_t));
1547 assert(ret != NULL);
1548 /* make a copy of the name, so we aren't dependent on the original */
1549 ret->name = (char *)malloc(sizeof(char)*(strlen(name)+1));
1550 assert(ret->name != NULL);
1551 strcpy(ret->name, name); /* we know ret->name is long enough */
1553 ret->tension_handler = tension_handler;
1554 ret->inverse_tension_handler = inverse_tension_handler;
1555 ret->tension_group = tension_group;
1556 ret->tension_params = tension_params;
1557 ret->destroy_tension = destroy_tension;
1558 ret->num_domains = num_domains;
1560 ret->num_out_trans = 0;
1561 ret->out_trans = NULL;
1562 ret->num_grouped_states = 0;
1563 ret->grouped_states = NULL;
1566 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);
1571 void destroy_state(state_t *state)
1575 fprintf(stderr, "freeing state %s (%p)\n", state->name, state);
1579 if (state->destroy_tension)
1580 (*state->destroy_tension)(state->tension_params);
1581 if (state->out_trans)
1582 free(state->out_trans);
1583 if (state->grouped_states)
1584 free(state->grouped_states);
1591 We'll be storing the states in a list (see Appendix \ref{sec.lists}),
1592 so we also define a convenience function for cleaning up.
1593 <<destroy state list>>=
1594 void destroy_state_list(list_t *state_list)
1596 state_list = head(state_list);
1597 while (state_list != NULL)
1598 destroy_state((state_t *) pop(&state_list));
1603 We can also define a convenient way to get a state index from it's name.
1605 int state_index(list_t *state_list, char *name)
1608 state_list = head(state_list);
1609 while (state_list != NULL) {
1610 if (STRMATCH(S(state_list)->name, name)) return i;
1611 state_list = state_list->next;
1619 \subsection{Transition data}
1620 \label{sec.transition_data}
1622 Once you've created states (Sections \ref{sec.main},
1623 \ref{sec.model_selection}, and \ref{sec.state_data}), you need to
1624 create transitions between them.
1625 <<transition definitions>>=
1626 typedef struct transition_struct {
1627 state_t *initial_state;
1628 state_t *final_state;
1629 k_func_t *k; /* transition rate model */
1630 void *k_params; /* pointer to k parameters */
1631 destroy_data_func_t *destroy_k;
1634 /* get the transition data for the current list node */
1635 #define T(list) ((transition_t *)(list)->d)
1637 /* get the number of initial state population for a transition state */
1638 #define N_INIT(transition) ((transition)->initial_state->num_domains)
1642 @ [[k]] is a pointer to the function determining the transition rate
1643 for a given tension. [[k_params]] and [[destroy_k]] have similar
1644 roles with regards to [[k]] as [[tension_params]] and
1645 [[destroy_tension]] have with regards to [[tension_handler]] (see
1646 Section \ref{sec.state_data}).
1648 [[create_]] and [[destroy_transition]] are simple wrappers around
1649 [[malloc]] and [[free]].
1650 <<create/destroy transition>>=
1651 transition_t *create_transition(state_t *initial_state, state_t *final_state,
1654 destroy_data_func_t *destroy_k)
1656 transition_t *ret = (transition_t *)malloc(sizeof(transition_t));
1657 assert(ret != NULL);
1658 assert(initial_state != NULL);
1659 assert(final_state != NULL);
1660 ret->initial_state = initial_state;
1661 ret->final_state = final_state;
1663 ret->k_params = k_params;
1664 ret->destroy_k = destroy_k;
1666 fprintf(stderr, "generated transition %p from %s to %s, params %p\n", ret, ret->initial_state->name, ret->final_state->name, ret->k_params);
1671 void destroy_transition(transition_t *transition)
1675 fprintf(stderr, "freeing transition %p (%s to %s)\n", transition, transition->initial_state->name, transition->final_state->name);
1677 if (transition->destroy_k)
1678 (*transition->destroy_k)(transition->k_params);
1685 We'll be storing the transitions in a list (see Appendix
1686 \ref{sec.lists}), so we also define a convenience function for
1688 <<destroy transition list>>=
1689 void destroy_transition_list(list_t *transition_list)
1691 transition_list = head(transition_list);
1692 while (transition_list != NULL)
1693 destroy_transition((transition_t *) pop(&transition_list));
1698 \subsection{Graphviz output}
1699 \label{sec.graphviz_output}
1701 <<print state graph>>=
1702 void print_state_graph(FILE *file, list_t *state_list, list_t *transition_list)
1704 state_list = head(state_list);
1705 transition_list = head(transition_list);
1706 fprintf(file, "digraph states {\n node [shape=record,color=black,fillcolor=white,style=filled];\n\n");
1708 fprintf(file, " n%d [label=\"%s\"]\n", state_index(state_list, S(state_list)->name), S(state_list)->name);
1709 if (state_list->next == NULL) break;
1710 state_list = state_list->next;
1712 fprintf(file, "\n");
1713 while (transition_list != NULL) {
1714 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));
1715 transition_list = transition_list->next;
1717 fprintf(file, "}\n");
1721 \subsection{Domain model and group handling}
1723 <<group functions>>=
1724 <<crosslink groups>>
1725 <<get active groups>>
1728 Scan through a list of states and construct an array of all the
1730 <<get active groups>>=
1731 void get_active_groups(list_t *state_list,
1732 int *pNum_active_groups,
1733 one_dim_fn_t ***pPer_group_handlers,
1734 one_dim_fn_t ***pPer_group_inverse_handlers,
1735 void ***pPer_group_data /* pointer to array of pointers to tension_handler_data_t */)
1737 void **active_groups=NULL;
1738 one_dim_fn_t *handler, *old_handler, *inverse_handler, *old_inverse_handler;
1739 one_dim_fn_t **per_group_handlers=NULL, **per_group_inverse_handlers=NULL;
1740 void **per_group_data=NULL; /* array of pointers to tension_handler_data_t */
1741 int i, j, num_domains, group, old_group, num_active_groups=0;
1742 list_t *active_states_list=NULL;
1743 tension_handler_data_t *tdata=NULL;
1746 /* get all the active groups in a list */
1747 state_list = head(state_list);
1749 fprintf(stderr, "%s:\t", __FUNCTION__);
1750 list_print(stderr, state_list, "states");
1752 while (state_list != NULL) {
1753 state = S(state_list);
1754 handler = state->tension_handler;
1755 inverse_handler = state->inverse_tension_handler;
1756 group = state->tension_group;
1757 num_domains = state->num_domains;
1758 if (list_index(active_states_list, state) == -1) {
1759 /* we haven't added this state (or it's associated group) yet */
1760 if (num_domains > 0 && handler != NULL) { /* add the state */
1761 num_active_groups++;
1762 push(&active_states_list, state, 1);
1764 fprintf(stderr, "%s:\tactive %s state (%p) with %d domain(s)\n", __FUNCTION__, state->name, state, state->num_domains);
1766 for (i=0; i < state->num_grouped_states; i++) {
1767 if (state->grouped_states[i]->num_domains > 0) {
1768 /* add this related (and active) state now */
1769 assert(state->grouped_states[i]->tension_handler == handler);
1770 assert(state->grouped_states[i]->tension_group == group);
1771 push(&active_states_list, state->grouped_states[i], 1);
1773 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);
1777 } /* else empty state or NULL handler */
1778 } /* else state already added */
1779 state_list = state_list->next;
1782 fprintf(stderr, "%s:\t", __FUNCTION__);
1783 list_print(stderr, active_states_list, "active states");
1786 assert(num_active_groups <= list_length(active_states_list));
1787 assert(num_active_groups > 0);
1789 /* allocate the arrays we'll be returning */
1790 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1791 assert(per_group_handlers != NULL);
1792 per_group_inverse_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1793 assert(per_group_inverse_handlers != NULL);
1794 per_group_data = (void **)malloc(sizeof(void *)*num_active_groups);
1795 assert(per_group_data != NULL);
1797 fprintf(stderr, "%s:\tallocating per group data %p:", __FUNCTION__, per_group_data);
1799 for (i=0; i<num_active_groups; i++) { /* we might want to recycle any group data structs allocated in the previous round */
1800 per_group_data[i] = malloc(sizeof(tension_handler_data_t));
1801 assert(per_group_data[i] != NULL);
1803 fprintf(stderr, " %d,%p", i, per_group_data[i]);
1807 fprintf(stderr, "\n");
1810 /* populate the arrays we'll be returning */
1811 active_states_list = head(active_states_list);
1813 old_inverse_handler = NULL;
1814 j = -1; /* j is the current active _group_ index */
1815 while (active_states_list != NULL) {
1816 state = (state_t *) pop(&active_states_list);
1817 handler = state->tension_handler;
1818 inverse_handler = state->inverse_tension_handler;
1819 group = state->tension_group;
1820 num_domains = state->num_domains;
1821 if (handler != old_handler || inverse_handler != old_inverse_handler || group != old_group) {
1822 j++; /* move to the next group */
1823 tdata = (tension_handler_data_t *) per_group_data[j];
1824 per_group_handlers[j] = handler;
1825 per_group_inverse_handlers[j] = inverse_handler;
1826 tdata->group_tension_params = NULL;
1828 tdata->persist = NULL; /* we might want to grab any appropriate data from the old persist here... */
1831 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);
1833 for (i=0; i<num_domains; i++) { /* add num_domains copies of the tension params */
1834 push(&tdata->group_tension_params, state->tension_params, 1);
1837 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);
1839 old_handler = handler;
1840 old_inverse_handler = inverse_handler;
1844 /* free old groups */
1845 if (*pPer_group_handlers != NULL)
1846 free(*pPer_group_handlers);
1847 if (*pPer_group_inverse_handlers != NULL)
1848 free(*pPer_group_inverse_handlers);
1849 if (*pPer_group_data != NULL) {
1850 for (i=0; i<*pNum_active_groups; i++)
1851 free((*pPer_group_data)[i]);
1852 free(*pPer_group_data);
1855 /* set new groups */
1857 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]);
1859 *pNum_active_groups = num_active_groups;
1860 *pPer_group_handlers = per_group_handlers;
1861 *pPer_group_inverse_handlers = per_group_inverse_handlers;
1862 *pPer_group_data = per_group_data;
1866 <<crosslink groups>>=
1867 void crosslink_groups(list_t *state_groups, list_t *transition_list)
1869 list_t *list, *out_trans=NULL, *associates=NULL;
1870 one_dim_fn_t *handler;
1873 state_groups = head(state_groups);
1874 while (state_groups != NULL) {
1875 /* find transitions leaving this state */
1877 fprintf(stderr, "finding transitions leaving %s\n", S(state_groups)->name);
1879 list = head(transition_list);
1880 while (list != NULL) {
1881 if (T(list)->initial_state == S(state_groups) && list_index(associates, T(list)->initial_state) != -1) {
1882 push(&out_trans, T(list), 1);
1886 n = list_length(out_trans);
1887 S(state_groups)->num_out_trans = n;
1888 S(state_groups)->out_trans = (transition_t **)malloc(sizeof(transition_t *)*n);
1889 assert(S(state_groups)->out_trans != NULL);
1890 for (i=0; i<n; i++) {
1891 S(state_groups)->out_trans[i] = (transition_t *)pop(&out_trans);
1894 /* find states grouped with this state */
1896 fprintf(stderr, "finding states grouped with %s\n", S(state_groups)->name);
1898 handler = S(state_groups)->tension_handler;
1899 group = S(state_groups)->tension_group;
1900 list = head(state_groups);
1901 while (list != NULL) {
1902 if (S(list)->tension_handler == handler && S(list)->tension_group == group && list != state_groups) {
1903 push(&associates, S(list), 1);
1907 n = list_length(associates);
1908 S(state_groups)->num_grouped_states = n;
1909 S(state_groups)->grouped_states = (state_t **)malloc(sizeof(state_t *)*n);
1910 assert(S(state_groups)->grouped_states != NULL);
1911 for (i=0; i<n; i++) {
1912 S(state_groups)->grouped_states[i] = (state_t *)pop(&associates);
1914 state_groups = state_groups->next;
1920 \section{String parsing}
1922 For handling command line arguments, we need parse delimited strings (e.g.\ [["param1,param2,param3..."]]).
1923 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1924 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1925 need to take care of parsing those parameters themselves.
1926 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1932 <<parse definitions>>
1933 <<parse declarations>>
1937 <<parse module makefile lines>>=
1938 NW_SAWSIM_MODS += parse
1939 CHECK_BINS += check_parse
1943 <<parse definitions>>=
1944 #define SEP ',' /* argument separator character */
1947 <<parse declarations>>=
1948 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1949 int *num, char ***string_array);
1950 extern void free_parsed_list(int num, char **string_array);
1953 [[parse_list_string]] allocates memory, don't forget to free it
1954 afterward with [[free_parsed_list]]. It does not alter the original.
1956 The string may start off with a [[deeper]] character
1957 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1958 [[x,y]], where the pointer is one character in on the copied string.
1959 However, when we go to free the memory, we need a pointer to the
1960 beginning of the string. In order to accommodate this for a string
1961 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1962 the first $N$ elements point to the separated arguments, and let the
1963 last element point to the start of the copied string regardless of
1965 <<parse delimited list string functions>>=
1966 /* TODO, split out into parse.hc */
1967 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1970 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1971 if (string[i] == deeper) {depth++;}
1972 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1978 void parse_list_string(char *string, char sep, char deeper, char shallower,
1979 int *num, char ***string_array)
1981 char *str=NULL, **ret=NULL;
1983 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
1985 *string_array = NULL;
1988 /* make a copy of the string, so we don't change the original */
1989 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1990 assert(str != NULL);
1991 strcpy(str, string); /* we know str is long enough */
1992 /* count the number of regions, so we can allocate pointers to them */
1995 n++; i++; /* move on to next argument */
1996 i += next_delim_index(str+i, sep, deeper, shallower);
1997 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1998 } while (str[i] != '\0');
1999 ret = (char **)malloc(sizeof(char *)*(n+1));
2000 assert(ret != NULL);
2001 /* replace the separators with '\0' & assign pointers */
2002 ret[n] = str; /* point to the front of the copied string */
2005 for(i=1; i<n; i++) {
2006 j += next_delim_index(str+j, sep, deeper, shallower);
2008 ret[i] = str+j; /* point to the separated arguments */
2010 /* strip off bounding braces, if any */
2011 for(i=0; i<n; i++) {
2012 if (ret[i][0] == deeper) {
2016 if (ret[i][strlen(ret[i])-1] == shallower)
2017 ret[i][strlen(ret[i])-1] = '\0';
2020 *string_array = ret;
2023 void free_parsed_list(int num, char **string_array)
2026 /* frees all items, since they were allocated as a single string*/
2027 free(string_array[num]);
2028 /* and free the array of pointers */
2034 \subsection{Parsing implementation}
2040 <<parse delimited list string functions>>
2044 #include <assert.h> /* assert() */
2045 #include <stdlib.h> /* NULL */
2046 #include <stdio.h> /* fprintf(), stdout *//*!!*/
2047 #include <string.h> /* strlen() */
2051 \subsection{Parsing unit tests}
2053 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2056 <<parse check includes>>
2057 <<string comparison definition and globals>>
2058 <<check relative error>>
2059 <<parse test suite>>
2060 <<main check program>>
2063 <<parse check includes>>=
2064 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2065 #include <stdio.h> /* printf() */
2066 #include <assert.h> /* assert() */
2067 #include <string.h> /* strlen() */
2072 <<parse test suite>>=
2073 <<parse list string tests>>
2074 <<parse suite definition>>
2077 <<parse suite definition>>=
2078 Suite *test_suite (void)
2080 Suite *s = suite_create ("parse");
2081 <<parse list string test case defs>>
2083 <<parse list string test case adds>>
2088 <<parse list string tests>>=
2091 START_TEST(test_next_delim_index)
2093 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
2094 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
2095 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
2096 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
2097 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
2101 START_TEST(test_parse_list_null)
2105 parse_list_string(NULL, SEP, '{', '}',
2106 &num_param_args, ¶m_args);
2107 fail_unless(num_param_args == 0, NULL);
2108 fail_unless(param_args == NULL, NULL);
2111 START_TEST(test_parse_list_single_simple)
2115 parse_list_string("arg", SEP, '{', '}',
2116 &num_param_args, ¶m_args);
2117 fail_unless(num_param_args == 1, NULL);
2118 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
2121 START_TEST(test_parse_list_single_compound)
2125 parse_list_string("{x,y,z}", SEP, '{', '}',
2126 &num_param_args, ¶m_args);
2127 fail_unless(num_param_args == 1, NULL);
2128 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
2131 START_TEST(test_parse_list_double_simple)
2135 parse_list_string("abc,def", SEP, '{', '}',
2136 &num_param_args, ¶m_args);
2137 fail_unless(num_param_args == 2, NULL);
2138 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2139 fail_unless(STRMATCH(param_args[1],"def"), NULL);
2142 START_TEST(test_parse_list_compound)
2146 parse_list_string("abc,{def,ghi}", SEP, '{', '}',
2147 &num_param_args, ¶m_args);
2148 fail_unless(num_param_args == 2, NULL);
2149 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2150 fail_unless(STRMATCH(param_args[1],"def,ghi"), NULL);
2154 <<parse list string test case defs>>=
2155 TCase *tc_parse_list_string = tcase_create("parse list string");
2157 <<parse list string test case adds>>=
2158 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
2159 tcase_add_test(tc_parse_list_string, test_parse_list_null);
2160 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
2161 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
2162 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
2163 tcase_add_test(tc_parse_list_string, test_parse_list_compound);
2164 suite_add_tcase(s, tc_parse_list_string);
2168 \section{Unit tests}
2170 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2177 <<check relative error>>
2180 <<main check program>>
2191 <<determine dt tests>>
2193 <<does domain transition tests>>
2194 <<randomly unfold domains tests>>
2195 <<suite definition>>
2198 <<suite definition>>=
2199 Suite *test_suite (void)
2201 Suite *s = suite_create ("sawsim");
2202 <<determine dt test case defs>>
2203 <<happens test case defs>>
2204 <<does domain transition test case defs>>
2205 <<randomly unfold domains test case defs>>
2207 <<determine dt test case adds>>
2208 <<happens test case adds>>
2209 <<does domain transition test case adds>>
2210 <<randomly unfold domains test case adds>>
2213 tcase_add_checked_fixture(tc_strip_address,
2214 setup_strip_address,
2215 teardown_strip_address);
2221 <<main check program>>=
2225 Suite *s = test_suite();
2226 SRunner *sr = srunner_create(s);
2227 srunner_run_all(sr, CK_ENV);
2228 number_failed = srunner_ntests_failed(sr);
2230 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
2235 \subsection{Model tests}
2237 Check the searching with [[linear_k]].
2238 Check overwhelming force treatment with the heavyside-esque step [[?]].
2240 Hmm, I'm not really sure what I was doing here...
2241 <<determine dt tests>>=
2242 double linear_k(double F, environment_t *env, void *params)
2244 double Fnot = *(double *)params;
2248 START_TEST(test_determine_dt_linear_k)
2251 transition_t unfold={0};
2252 environment_t env = {0};
2253 double F[]={0,1,2,3,4,5,6};
2254 double dt_max=3.0, Fnot=3.0, x=1.0, max_p=0.1, max_dF=0.1, v=1.0;
2255 list_t *state_list=NULL, *transition_list=NULL;
2258 folded.tension_handler = &hooke_handler;
2259 folded.num_domains = 1;
2260 unfold.initial_state = &folded;
2261 unfold.k = linear_k;
2262 unfold.k_params = &Fnot;
2263 push(&state_list, &folded, 1);
2264 push(&transition_list, &unfold, 1);
2265 for( i=0; i < sizeof(F)/sizeof(double); i++) {
2266 //fail_unless(determine_dt(state_list, transition_list, &env, x, max_p, max_dF, dt_max, v)==1, NULL);
2271 <<determine dt test case defs>>=
2272 TCase *tc_determine_dt = tcase_create("Determine dt");
2274 <<determine dt test case adds>>=
2275 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
2276 suite_add_tcase(s, tc_determine_dt);
2281 <<happens test case defs>>=
2283 <<happens test case adds>>=
2286 <<does domain transition tests>>=
2288 <<does domain transition test case defs>>=
2290 <<does domain transition test case adds>>=
2293 <<randomly unfold domains tests>>=
2295 <<randomly unfold domains test case defs>>=
2297 <<randomly unfold domains test case adds>>=
2301 \section{Balancing group extensions}
2302 \label{sec.tension_balance}
2304 For a given total extension $x$ (of the piezo), the various domain
2305 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
2306 amounts, and we need to tweak the portion that each extends to balance
2307 the tension amongst the active groups. Since the tension balancing is
2308 separable from the bulk of the simulation, we break it out into a
2309 separate module. The interface is defined in [[tension_balance.h]],
2310 the implementation in [[tension_balance.c]], and the unit testing in
2311 [[check_tension_balance.c]]
2313 <<tension-balance.h>>=
2314 #ifndef TENSION_BALANCE_H
2315 #define TENSION_BALANCE_H
2317 <<tension balance function declaration>>
2318 #endif /* TENSION_BALANCE_H */
2321 <<tension balance functions>>=
2322 <<one dimensional bracket>>
2323 <<one dimensional solve>>
2325 <<group stiffness function>>
2326 <<chain stiffness function>>
2327 <<full chain stiffness function>>
2328 <<tension balance function>>
2331 <<tension balance module makefile lines>>=
2332 NW_SAWSIM_MODS += tension_balance
2333 CHECK_BINS += check_tension_balance
2334 check_tension_balance_MODS =
2337 The entire force balancing problem reduces to a solving two nested
2338 one-dimensional root-finding problems. First we define the one
2339 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
2340 populated group). $x(x_0)$ is also strictly monotonically increasing,
2341 so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
2342 stores the last successful value of $x$, since we expect to be taking
2343 small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
2344 indicates that the groups have changed.
2345 <<tension balance function declaration>>=
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 <<tension balance function>>=
2355 double tension_balance(int num_tension_groups,
2356 one_dim_fn_t **tension_handlers,
2357 one_dim_fn_t **inverse_tension_handlers,
2358 void **params, /* array of pointers to tension_handler_data_t */
2363 static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL,NULL};
2364 double F, xo_guess, xo, lb, ub;
2365 double min_relx=1e-6, min_rely=1e-6;
2366 int max_steps=200, i;
2368 assert(num_tension_groups > 0);
2369 assert(tension_handlers != NULL);
2370 assert(params != NULL);
2371 assert(num_tension_groups > 0);
2373 if (last_x == -1 || x_xo_data.xi == NULL) { /* new group setup, reset x_xo_data */
2374 /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
2375 if (x_xo_data.xi != NULL)
2377 /* construct new x_xo_data */
2378 x_xo_data.n_groups = num_tension_groups;
2379 x_xo_data.tension_handler = tension_handlers;
2380 x_xo_data.inverse_tension_handler = inverse_tension_handlers;
2381 x_xo_data.handler_data = params;
2383 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);
2385 x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
2386 assert(x_xo_data.xi != NULL);
2387 for (i=0; i<num_tension_groups; i++)
2388 x_xo_data.xi[i] = last_x;
2391 if (num_tension_groups == 1) { /* only one group, no balancing required */
2393 x_xo_data.xi[0] = xo;
2395 fprintf(stderr, "%s:\tonly one active group, no balancing required\n", __FUNCTION__);
2399 fprintf(stderr, "%s:\tbalancing for x = %g with", __FUNCTION__, x);
2400 for(i=0; i<num_tension_groups; i++) fprintf(stderr, "\tx[%d]=%g", i, x_xo_data.xi[i]);
2401 fprintf(stderr, "\n");
2403 if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
2404 if (x == 0) xo_guess = length_scale;
2405 else xo_guess = x/num_tension_groups;
2407 fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
2409 oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
2410 } else { /* work off the last balanced point */
2412 fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
2415 lb = x_xo_data.xi[0];
2416 ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */
2417 } else if (x < last_x) {
2418 lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
2419 ub = x_xo_data.xi[0];
2420 } else { /* x == last_x */
2422 fprintf(stderr, "not moving\n");
2424 lb = x_xo_data.xi[0];
2425 ub = x_xo_data.xi[0];
2429 fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
2430 __FUNCTION__, x, lb, ub);
2432 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
2433 F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
2435 if (num_tension_groups > 1) {
2436 fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
2437 for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
2438 fprintf(stderr, "\n");
2443 F = (*tension_handlers[0])(xo, params[0]);
2445 if (stiffness != NULL) {
2446 *stiffness = chain_stiffness(&x_xo_data, min_relx);
2447 if (*stiffness == 0.0) { /* re-calculate based on full chain */
2448 *stiffness = full_chain_stiffness(num_tension_groups, tension_handlers,
2449 inverse_tension_handlers, params,
2450 last_x, x, min_relx, F);
2458 <<tension balance internal definitions>>=
2459 <<x of x0 definitions>>
2462 <<x of x0 definitions>>=
2463 typedef struct x_of_xo_data_struct {
2465 one_dim_fn_t **tension_handler; /* array of fn pointers */
2466 one_dim_fn_t **inverse_tension_handler; /* inverse fn pointers */
2467 void **handler_data; /* array of pointers to tension_handler_data_t */
2468 double *xi; /* array of group extensions */
2473 double x_of_xo(double xo, void *pdata)
2475 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
2476 double F, x=0, xi, xi_guess, lb, ub, slope;
2478 double min_relx=1e-6, min_rely=1e-6;
2480 assert(data->n_groups > 0);
2483 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);
2486 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
2488 fprintf(stderr, "%s:\tfinding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
2492 if (data->xi) data->xi[0] = xo;
2493 for (i=1; i < data->n_groups; i++) {
2494 if (data->inverse_tension_handler[i] != NULL) {
2495 xi = (*data->inverse_tension_handler[i])(F, data->handler_data[i]);
2496 } else { /* invert numerically */
2497 if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
2498 else xi_guess = data->xi[i];
2500 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]);
2502 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
2503 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
2508 fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
2512 fprintf(stderr, "%s:\tfound x(x0=%g) = %g\n", __FUNCTION__, xo, x);
2518 Determine the stiffness of a single tension group by taking a small
2519 step [[dx]] back from the position [[x]] for which we want the
2520 stiffness. The touchy part is determining a good [[dx]] and ensuring
2521 the whole interval is on [[x>=0]].
2522 <<group stiffness function>>=
2523 double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double x, double relx)
2525 double dx, xl, F, Fl, stiffness;
2527 assert(relx > 0 && relx < 1);
2528 if (x == 0 || relx == 0) {
2529 dx = 1e-12; /* HACK, how to get length scale? */
2539 F = (*tension_handler)(x, handler_data);
2540 Fl = (*tension_handler)(xl, handler_data);
2541 stiffness = (F-Fl)/dx;
2542 assert(stiffness >= 0);
2547 Attempt to calculate the chain stiffness by adding group stiffnesses
2548 as springs in series. In the case where a group has undefined
2549 stiffness (e.g. piston tension, Section \ref{sec.piston_tension}),
2550 this algorithm breaks down. In that case, [[tension_balance()]]
2551 (\ref{sec.tension_balance}) should use [[full_chain_stiffness()]]
2552 which uses the full chain's $F(x)$ rather than that of the individual
2553 domains'. We attempt the series approach first because, lacking the
2554 numerical inversion inside [[tension_balance()]], it is both faster
2555 and more accurate than the full-chain derivative.
2556 <<chain stiffness function>>=
2557 double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx)
2559 double stiffness, gstiff;
2561 /* add group stiffnesses in series */
2562 for (i=0; i < x_xo_data->n_groups; i++) {
2564 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);
2567 gstiff = group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx);
2570 stiffness += 1.0/gstiff;
2572 stiffness = 1.0 / stiffness;
2578 Determine the chain stiffness using only [[tension_balance()]]. This
2579 works around difficulties with tension models that have undefined
2580 stiffness (see discussion for [[chain_stiffness()]] above).
2581 <<full chain stiffness function>>=
2582 double full_chain_stiffness(int num_tension_groups,
2583 one_dim_fn_t **tension_handlers,
2584 one_dim_fn_t **inverse_tension_handlers,
2585 void **params, /* array of pointers to tension_handler_data_t */
2589 double F /* already evaluated F(x) */)
2591 double dx, xl, Fl, stiffness;
2594 assert(relx > 0 && relx < 1);
2596 /* Other option for dx that we ignore because of our poor tension_balance()
2597 * resolution (i.e. bad slopes if you zoom in too close):
2598 * if (last_x != -1.0 && last_x != x) // last_x limited
2599 * dx fabs(x-last_x);
2602 * if (1==1) { // constant max-value limited
2604 * dx = MIN(dx, dx_p);
2606 * if (x != 0 && relx != 0) { // relx limited
2608 * dx = MIN(dx, dx_p);
2610 * TODO, determine which of (min_relx, min_rely, max_steps) is actually
2611 * limiting tension_balance.
2613 dx = 1e-12; /* HACK, how to get length scale? */
2617 Fl = tension_balance(num_tension_groups, tension_handlers,
2618 inverse_tension_handlers, params,
2624 F = tension_balance(num_tension_groups, tension_handlers,
2625 inverse_tension_handlers, params,
2628 stiffness = (F-Fl)/dx;
2629 assert(stiffness >= 0);
2635 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
2636 number of steps for monotonically increasing $f(x)$. Simple
2637 bisection, so it's robust and converges fairly quickly. You can set
2638 the maximum number of steps to take, but if you have enough processor
2639 power, it's better to limit your precision with the relative limits
2640 [[min_relx]] and [[y]]. These ensure that the width of the bracket is
2641 small on the length scale of it's larger side. Note that \emph{both}
2642 relative limits must be satisfied for the search to stop.
2643 <<one dimensional function definition>>=
2644 /* equivalent to gsl_func_t */
2645 typedef double one_dim_fn_t(double x, void *params);
2647 <<one dimensional solve declaration>>=
2648 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2649 double min_relx, double min_rely, int max_steps,
2653 <<one dimensional solve>>=
2654 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2655 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2656 double min_relx, double min_rely, int max_steps,
2659 double yx, ylb, yub, x;
2662 ylb = (*f)(lb, params);
2663 yub = (*f)(ub, params);
2665 /* check some simple cases */
2667 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bracketed */
2668 else return lb; /* any x on the range [lb, ub] would work */
2670 if (ylb == y) { x=lb; yx=ylb; goto end; }
2671 if (yub == y) { x=ub; yx=yub; goto end; }
2674 fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2676 assert(ylb < y); assert(yub > y);
2677 for (i=0; i<max_steps; i++) {
2679 yx = (*f)(x, params);
2681 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);
2683 if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2684 else if (yx > y) { ub=x; yub=yx; }
2685 else /* yx < y */ { lb=x; ylb=yx; }
2690 fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2696 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2697 Generate bracketing $x$ values through bisection or exponential growth.
2698 <<one dimensional bracket declaration>>=
2699 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2702 <<one dimensional bracket>>=
2703 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2705 double yx, step, x=xguess;
2708 fprintf(stderr, "%s:\tbracketing %g (params %p)\n", __FUNCTION__, y, params);
2710 yx = (*f)(x, params);
2712 fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2719 fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2723 if (x == 0) x = length_scale/2.0; /* HACK */
2726 yx = (*f)(x, params);
2728 fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2733 fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2736 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2740 \subsection{Balancing implementation}
2742 <<tension-balance.c>>=
2745 <<tension balance includes>>
2746 #include "tension_balance.h"
2748 double length_scale = 1e-10; /* HACK */
2750 <<tension balance internal definitions>>
2751 <<tension balance functions>>
2754 <<tension balance includes>>=
2755 #include <assert.h> /* assert() */
2756 #include <stdlib.h> /* NULL */
2757 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2758 #include <stdio.h> /* fprintf(), stdout */
2762 \subsection{Balancing unit tests}
2764 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2765 <<check-tension-balance.c>>=
2767 <<tension balance check includes>>
2768 <<tension balance test suite>>
2769 <<main check program>>
2772 <<tension balance check includes>>=
2773 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2774 #include <assert.h> /* assert() */
2777 #include "tension_balance.h"
2780 <<tension balance test suite>>=
2781 <<tension balance function tests>>
2782 <<tension balance suite definition>>
2785 <<tension balance suite definition>>=
2786 Suite *test_suite (void)
2788 Suite *s = suite_create ("tension balance");
2789 <<tension balance function test case defs>>
2791 <<tension balance function test case adds>>
2796 <<tension balance function tests>>=
2797 <<check relative error>>
2799 double hooke(double x, void *pK)
2802 return *((double*)pK) * x;
2805 START_TEST(test_single_function)
2807 double k=5, x=3, last_x=2.0, F;
2808 one_dim_fn_t *handlers[] = {&hooke};
2809 one_dim_fn_t *inverse_handlers[] = {NULL};
2810 void *data[] = {&k};
2811 F = tension_balance(1, handlers, inverse_handlers, data, last_x, x, NULL);
2812 fail_unless(F = k*x, NULL);
2817 We can also test balancing two springs with different spring constants.
2818 The analytic solution for a general number of springs is given in Appendix \ref{sec.math_hooke}.
2819 <<tension balance function tests>>=
2820 START_TEST(test_double_hooke)
2822 double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
2823 one_dim_fn_t *handlers[] = {&hooke, &hooke};
2824 one_dim_fn_t *inverse_handlers[] = {NULL, NULL};
2825 void *data[] = {&k1, &k2};
2826 F = tension_balance(2, handlers, inverse_handlers, data, last_x, x, NULL);
2830 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2831 //CHECK_ERR(1e-6, x1e, xi[0]);
2832 //CHECK_ERR(1e-6, x2e, xi[1]);
2833 CHECK_ERR(1e-6, Fe, F);
2837 TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
2839 <<tension balance function test case defs>>=
2840 TCase *tc_tbfunc = tcase_create("tension balance function");
2843 <<tension balance function test case adds>>=
2844 tcase_add_test(tc_tbfunc, test_single_function);
2845 tcase_add_test(tc_tbfunc, test_double_hooke);
2846 suite_add_tcase(s, tc_tbfunc);
2852 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2853 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2854 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2860 <<list definitions>>
2861 <<list declarations>>
2865 <<list declarations>>=
2866 <<head and tail declarations>>
2867 <<list length declaration>>
2868 <<push declaration>>
2870 <<list destroy declaration>>
2871 <<list to array declaration>>
2872 <<move declaration>>
2873 <<sort declaration>>
2874 <<list index declaration>>
2875 <<list element declaration>>
2876 <<unique declaration>>
2877 <<list print declaration>>
2881 <<create and destroy node>>
2896 <<list module makefile lines>>=
2897 NW_SAWSIM_MODS += list
2898 CHECK_BINS += check_list
2902 <<list definitions>>=
2903 typedef struct list_struct {
2904 struct list_struct *next;
2905 struct list_struct *prev;
2909 <<comparison function typedef>>
2912 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2913 <<head and tail declarations>>=
2914 list_t *head(list_t *list);
2915 list_t *tail(list_t *list);
2918 list_t *head(list_t *list)
2922 while (list->prev != NULL) {
2928 list_t *tail(list_t *list)
2932 while (list->next != NULL) {
2939 <<list length declaration>>=
2940 int list_length(list_t *list);
2943 int list_length(list_t *list)
2950 while (list->next != NULL) {
2958 [[push]] inserts a node relative to the current position in the list
2959 without changing the current position.
2960 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2961 so we set it to that of the pushed domain.
2962 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2963 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2964 <<push declaration>>=
2965 void push(list_t **pList, void *data, int below);
2968 void push(list_t **pList, void *data, int below)
2970 list_t *list, *new_node;
2971 assert(pList != NULL);
2973 new_node = create_node(data);
2976 else if (below==0) { /* insert above */
2977 if (list->prev != NULL)
2978 list->prev->next = new_node;
2979 new_node->prev = list->prev;
2980 list->prev = new_node;
2981 new_node->next = list;
2982 } else { /* insert below */
2983 if (list->next != NULL)
2984 list->next->prev = new_node;
2985 new_node->next = list->next;
2986 list->next = new_node;
2987 new_node->prev = list;
2992 [[pop]] removes the current domain node, moving the current position
2993 to the node after it, unless that node is [[NULL]], in which case move
2994 the current position to the node before it.
2995 <<pop declaration>>=
2996 void *pop(list_t **pList);
2999 void *pop(list_t **pList)
3001 list_t *list, *popped;
3003 assert(pList != NULL);
3005 assert(list != NULL); /* not an empty list */
3007 /* bypass the popped node */
3008 if (list->prev != NULL)
3009 list->prev->next = list->next;
3010 if (list->next != NULL) {
3011 list->next->prev = list->prev;
3012 *pList = list->next;
3014 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
3016 destroy_node(popped);
3021 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
3022 If the cleanup function is [[NULL]], no cleanup function is called.
3023 The nodes are not popped in any particular order.
3024 <<list destroy declaration>>=
3025 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
3028 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
3032 assert(pList != NULL);
3035 assert(list != NULL); /* not an empty list */
3036 while (list != NULL) {
3038 if (destroy != NULL)
3044 [[list_to_array]] converts a list to an array of pointers, preserving order.
3045 <<list to array declaration>>=
3046 void list_to_array(list_t *list, int *length, void ***array);
3049 void list_to_array(list_t *list, int *pLength, void ***pArray)
3053 assert(list != NULL);
3054 assert(pLength != NULL);
3055 assert(pArray != NULL);
3057 length = list_length(list);
3058 array = (void **)malloc(sizeof(void *)*length);
3059 assert(array != NULL);
3062 while (list != NULL) {
3063 array[i++] = list->d;
3071 [[move]] moves the current element along the list in either direction.
3072 <<move declaration>>=
3073 void move(list_t **pList, int downstream);
3076 void move(list_t **pList, int downstream)
3078 list_t *list, *flipper;
3080 assert(pList != NULL);
3081 list = *pList; /* a>B>c>d */
3082 assert(list != NULL); /* not an empty list */
3083 if (downstream == 0)
3084 flipper = list->prev; /* flipper is A */
3086 flipper = list->next; /* flipper is C */
3087 /* ensure there is somewhere to go in the direction we're heading */
3088 assert(flipper != NULL);
3089 /* remove the current element */
3090 data = pop(&list); /* data=B, a>C>d */
3091 /* and put it back in in the correct spot */
3093 if (downstream == 0) {
3094 push(&list, data, 0); /* b>A>c>d */
3095 list = list->prev; /* B>a>c>d */
3097 push(&list, data, 1); /* a>C>b>d */
3098 list = list->next; /* a>c>B>d */
3104 [[sort]] sorts a list from smallest to largest according to a user
3105 specified comparison.
3106 <<comparison function typedef>>=
3107 typedef int (comparison_fn_t)(void *A, void *B);
3110 <<int comparison function>>=
3111 int int_comparison(void *A, void *B)
3116 if (a > b) return 1;
3117 else if (a == b) return 0;
3122 <<sort declaration>>=
3123 void sort(list_t **pList, comparison_fn_t *comp);
3125 Here we use the bubble sort algorithm.
3127 void sort(list_t **pList, comparison_fn_t *comp)
3131 assert(pList != NULL);
3133 assert(list != NULL); /* not an empty list */
3134 while (stable == 0) {
3137 while (list->next != NULL) {
3138 if ((*comp)(list->d, list->next->d) > 0) {
3140 move(&list, 1 /* downstream */);
3150 [[list_index]] finds the location of [[data]] in [[list]]. Returns
3151 [[-1]] if [[data]] is not in [[list]].
3152 <<list index declaration>>=
3153 int list_index(list_t *list, void *data);
3157 int list_index(list_t *list, void *data)
3161 while (list != NULL) {
3162 if (list->d == data) return i;
3171 [[list_element]] returns the data bound to the $i^\text{th}$ node of the list.
3172 <<list element declaration>>=
3173 void *list_element(list_t *list, int i);
3176 void *list_element(list_t *list, int i)
3180 while (list != NULL) {
3181 if (i == j) return list->d;
3191 [[unique]] removes duplicates from a sorted list according to a user
3192 specified comparison. Don't do this unless you have other pointers
3193 any dynamically allocated data in your list, because [[unique]] just
3194 drops any non-unique data without freeing it.
3195 <<unique declaration>>=
3196 void unique(list_t **pList, comparison_fn_t *comp);
3199 void unique(list_t **pList, comparison_fn_t *comp)
3202 assert(pList != NULL);
3204 assert(list != NULL); /* not an empty list */
3206 while (list->next != NULL) {
3207 if ((*comp)(list->d, list->next->d) == 0) {
3216 [[list_print]] prints a list to a [[FILE *]] stream.
3217 <<list print declaration>>=
3218 void list_print(FILE *file, list_t *list, char *list_name);
3221 void list_print(FILE *file, list_t *list, char *list_name)
3223 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
3225 while (list != NULL) {
3226 fprintf(file, " %p", list->d);
3229 fprintf(file, "\n");
3233 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
3234 <<create and destroy node>>=
3235 list_t *create_node(void *data)
3237 list_t *ret = (list_t *)malloc(sizeof(list_t));
3238 assert(ret != NULL);
3245 void destroy_node(list_t *node)
3251 The user must free the data pointed to by the node on their own.
3253 \subsection{List implementation}
3263 #include <assert.h> /* assert() */
3264 #include <stdlib.h> /* malloc(), free(), rand() */
3265 #include <stdio.h> /* fprintf(), stdout, FILE */
3266 #include "global.h" /* destroy_data_func_t */
3269 \subsection{List unit tests}
3271 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3274 <<list check includes>>
3276 <<main check program>>
3279 <<list check includes>>=
3280 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3281 #include <stdio.h> /* FILE */
3287 <<list test suite>>=
3290 <<list suite definition>>
3293 <<list suite definition>>=
3294 Suite *test_suite (void)
3296 Suite *s = suite_create ("list");
3297 <<push test case defs>>
3298 <<pop test case defs>>
3300 <<push test case adds>>
3301 <<pop test case adds>>
3307 START_TEST(test_push)
3312 push(&list, (void *)i, 1);
3313 fail_unless(list == head(list), NULL);
3314 fail_unless( (int)list->d == 0 );
3315 for (i=0; i<n; i++) {
3316 /* we expect 0, n-1, n-2, ..., 1 */
3319 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
3325 <<push test case defs>>=
3326 TCase *tc_push = tcase_create("push");
3329 <<push test case adds>>=
3330 tcase_add_test(tc_push, test_push);
3331 suite_add_tcase(s, tc_push);
3336 <<pop test case defs>>=
3338 <<pop test case adds>>=
3341 \section{Function string evaluation}
3343 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).
3344 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
3345 The solution is to run a scripting language as a subprocess accessed via pipes.
3346 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
3348 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
3349 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
3350 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.
3351 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
3353 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.
3354 Then you can either statically or dynamically link to those hardcoded functions.
3355 While much less flexible, this approach would be a fairly simple-to-implement fallback.
3357 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
3358 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
3361 #ifndef STRING_EVAL_H
3362 #define STRING_EVAL_H
3364 <<string eval setup declaration>>
3365 <<string eval function declaration>>
3366 <<string eval teardown declaration>>
3367 #endif /* STRING_EVAL_H */
3370 <<string eval module makefile lines>>=
3371 NW_SAWSIM_MODS += string_eval
3372 CHECK_BINS += check_string_eval
3373 check_string_eval_MODS =
3376 For an introduction to POSIX process control, see\\
3377 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
3378 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
3379 and of course, the relavent [[man]] pages.
3381 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
3382 [[execvp]] replaces the calling process' program with a new program.
3383 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
3384 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
3385 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
3386 The new program needs command line arguments, just like it would if you were running it from a shell.
3387 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
3388 with the final array entry being a [[NULL]] pointer.
3390 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
3391 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
3392 <<string eval subprocess definitions>>=
3393 #define SUBPROCESS "python"
3394 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
3395 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
3396 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
3398 The [[i]] option lets Python know that it should run in interactive mode.
3399 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
3400 In interactive mode, python acts on every instruction as soon as it is recieved.
3401 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.
3402 %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.
3404 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
3405 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
3406 [[fork]] returns two copies of the same program, executing the original code.
3407 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
3408 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
3410 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.
3411 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
3412 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
3413 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
3414 <<string eval pipe definitions>>=
3415 #define PIPE_READ 0 /* the end you read from */
3416 #define PIPE_WRITE 1 /* the end you write to */
3418 #define STDIN 0 /* initial index of stdin pair */
3419 #define STDOUT 2 /* initial index of stdout pair */
3422 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
3424 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.
3426 <<string eval setup declaration>>=
3427 extern void string_eval_setup(FILE **pIn, FILE **pOut);
3429 <<string eval setup definition>>=
3430 void string_eval_setup(FILE **pIn, FILE **pOut)
3433 int pfd[4]; /* file descriptors for stdin and stdout */
3435 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
3436 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
3438 pid = fork(); /* split process into two copies */
3440 if (pid == -1) { /* fork error */
3441 perror("fork error");
3443 } else if (pid == 0) { /* child */
3444 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
3445 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
3446 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
3447 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
3448 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
3449 perror("exec error"); /* exec shouldn't return */
3451 } else { /* parent */
3452 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
3453 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
3454 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
3455 if ( *pIn == NULL ) {
3456 perror("fdopen (in)");
3459 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
3460 if ( *pOut == NULL ) {
3461 perror("fdopen (out)");
3468 To use the evaluating subprocess, we just pipe in our command, and read out the results.
3469 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
3470 <<string eval function declaration>>=
3471 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
3473 <<string eval function definition>>=
3474 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
3477 rc = fprintf(in, "%s", input);
3478 assert(rc == strlen(input));
3481 alarm(1); /* set a one second timeout on the read */
3482 assert( fgets(output, buflen, out) != NULL );
3483 alarm(0); /* cancel the timeout */
3484 //fprintf(stderr, "eval: %s ----> %s", input, output);
3487 The [[alarm]] calls set and clear a timeout on the returned output.
3488 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.
3489 This protects against invalid input for which a line of output is not printed to [[stdout]].
3490 Other invalid inputs (e.g.\ those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
3491 If you are getting strange results, check your python code seperately. TODO, better error handling.
3493 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
3494 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
3495 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
3496 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
3497 <<string eval teardown declaration>>=
3498 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
3501 <<string eval teardown definition>>=
3502 void string_eval_teardown(FILE **pIn, FILE **pOut)
3507 /* redirect Python's stderr */
3508 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
3512 assert( fclose(*pIn) == 0 );
3514 assert( fclose(*pOut) == 0 );
3517 /* wait for python to exit */
3519 pid = wait(&stat_loc);
3526 if (WIFEXITED(stat_loc)) {
3527 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
3528 } else if (WIFSIGNALED(stat_loc)) {
3529 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
3534 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
3536 \subsection{String evaluation implementation}
3540 <<string eval includes>>
3541 #include "string_eval.h"
3542 <<string eval internal definitions>>
3543 <<string eval functions>>
3546 <<string eval includes>>=
3547 #include <assert.h> /* assert() */
3548 #include <stdlib.h> /* NULL */
3549 #include <stdio.h> /* fprintf(), stdout, fdopen() */
3550 #include <string.h> /* strlen() */
3551 #include <sys/types.h> /* pid_t */
3552 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
3553 #include <wait.h> /* wait() */
3556 <<string eval internal definitions>>=
3557 <<string eval subprocess definitions>>
3558 <<string eval pipe definitions>>
3561 <<string eval functions>>=
3562 <<string eval setup definition>>
3563 <<string eval function definition>>
3564 <<string eval teardown definition>>
3567 \subsection{String evaluation unit tests}
3569 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3570 <<check-string-eval.c>>=
3572 <<string eval check includes>>
3573 <<string comparison definition and globals>>
3574 <<string eval test suite>>
3575 <<main check program>>
3578 <<string eval check includes>>=
3579 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3580 #include <stdio.h> /* printf() */
3581 #include <assert.h> /* assert() */
3582 #include <string.h> /* strlen() */
3583 #include <signal.h> /* SIGKILL */
3585 #include "string_eval.h"
3588 <<string eval test suite>>=
3589 <<string eval tests>>
3590 <<string eval suite definition>>
3593 <<string eval suite definition>>=
3594 Suite *test_suite (void)
3596 Suite *s = suite_create ("string eval");
3597 <<string eval test case defs>>
3599 <<string eval test case adds>>
3604 <<string eval tests>>=
3605 START_TEST(test_setup_teardown)
3608 string_eval_setup(&in, &out);
3609 string_eval_teardown(&in, &out);
3612 START_TEST(test_invalid_command)
3615 char input[80], output[80]={};
3616 string_eval_setup(&in, &out);
3617 sprintf(input, "print undefined_name__traceback_should_be_printed_to_stderr\n");
3618 string_eval(in, out, input, 80, output);
3619 string_eval_teardown(&in, &out);
3622 START_TEST(test_simple_eval)
3625 char input[80], output[80]={};
3626 string_eval_setup(&in, &out);
3627 sprintf(input, "print 3+4*5\n");
3628 string_eval(in, out, input, 80, output);
3629 fail_unless(STRMATCH(output,"23\n"), NULL);
3630 string_eval_teardown(&in, &out);
3633 START_TEST(test_multiple_evals)
3636 char input[80], output[80]={};
3637 string_eval_setup(&in, &out);
3638 sprintf(input, "print 3+4*5\n");
3639 string_eval(in, out, input, 80, output);
3640 fail_unless(STRMATCH(output,"23\n"), NULL);
3641 sprintf(input, "print (3**2 + 4**2)**0.5\n");
3642 string_eval(in, out, input, 80, output);
3643 fail_unless(STRMATCH(output,"5.0\n"), NULL);
3644 string_eval_teardown(&in, &out);
3647 START_TEST(test_eval_with_state)
3650 char input[80], output[80]={};
3651 string_eval_setup(&in, &out);
3652 sprintf(input, "print 3+4*5\n");
3653 fprintf(in, "A = 3\n");
3654 sprintf(input, "print A*3\n");
3655 string_eval(in, out, input, 80, output);
3656 fail_unless(STRMATCH(output,"9\n"), NULL);
3657 string_eval_teardown(&in, &out);
3661 <<string eval test case defs>>=
3662 TCase *tc_string_eval = tcase_create("string_eval");
3664 <<string eval test case adds>>=
3665 tcase_add_test(tc_string_eval, test_setup_teardown);
3666 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
3667 tcase_add_test(tc_string_eval, test_simple_eval);
3668 tcase_add_test(tc_string_eval, test_multiple_evals);
3669 tcase_add_test(tc_string_eval, test_eval_with_state);
3670 suite_add_tcase(s, tc_string_eval);
3674 \section{Accelerating function evaluation}
3676 My first version-0.3 code was running very slowly.
3677 With the options suggested in the help
3678 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
3679 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
3680 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
3681 That's only 15 calls per solution, so the search algorithm seems reasonable.
3682 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
3688 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
3690 #endif /* ACCEL_K_H */
3693 <<accel k module makefile lines>>=
3694 NW_SAWSIM_MODS += accel_k
3695 CHECK_BINS += check_accel_k
3696 check_accel_k_MODS = interp k_model parse string_eval tavl
3701 #include <assert.h> /* assert() */
3702 #include <stdlib.h> /* realloc(), free(), NULL */
3703 #include "global.h" /* environment_t */
3704 #include "k_model.h" /* k_func_t */
3705 #include "interp.h" /* interp_* */
3706 #include "accel_k.h"
3708 typedef struct accel_k_struct {
3709 interp_table_t *itable;
3715 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3716 static int num_accels = 0;
3717 static accel_k_t *accels=NULL;
3719 /* Wrap k in the standard f(x) acceleration form */
3720 static double k_wrap(double F, void *params)
3722 accel_k_t *p = (accel_k_t *)params;
3723 return (*p->k)(F, p->env, p->params);
3726 static int k_tol(double FA, double kA, double FB, double kB)
3729 if (FB-FA > 1e-12) {
3730 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3731 return 1; /* unnacceptable */
3733 //printf("acceptable tol\n");
3734 return 0; /* acceptable */
3738 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3741 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3742 assert(accels != NULL);
3743 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3745 accels[i].env = env;
3746 accels[i].params = params;
3753 for (i=0; i<num_accels; i++)
3754 interp_table_free(accels[i].itable);
3756 if (accels) free(accels);
3760 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3763 for (i=0; i<num_accels; i++) {
3764 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3765 /* We've already setup this function.
3766 * WARNING: we're assuming param and env are the same. */
3767 return interp_table_eval(accels[i].itable, accels+i, F);
3770 /* set up a new acceleration environment */
3771 i = add_accel_k(k, env, params);
3772 return interp_table_eval(accels[i].itable, accels+i, F);
3776 \subsection{Unit tests}
3778 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3779 <<check-accel-k.c>>=
3781 <<accel k check includes>>
3782 <<check relative error>>
3784 <<accel k test suite>>
3785 <<main check program>>
3788 <<accel k check includes>>=
3789 #include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
3790 #include <stdio.h> /* sprintf() */
3791 #include <assert.h> /* assert() */
3792 #include <math.h> /* exp() */
3793 #include <errno.h> /* errno, ERANGE */
3796 #include "k_model.h"
3797 #include "accel_k.h"
3801 <<accel k test suite>>=
3802 <<accel k suite tests>>
3803 <<accel k suite definition>>
3806 <<accel k suite definition>>=
3807 Suite *test_suite (void)
3809 Suite *s = suite_create ("accelerated k model");
3810 <<accel k test case defs>>
3812 <<accel k test case adds>>
3817 <<accel k test case defs>>=
3818 TCase *tc_accel_k = tcase_create("accelerated k model");
3821 <<accel k test case adds>>=
3822 tcase_add_test(tc_accel_k, test_accel_k_accuracy);
3823 suite_add_tcase(s, tc_accel_k);
3826 <<accel k suite tests>>=
3827 START_TEST(test_accel_k_accuracy)
3829 double k, k_accel, F, Fm=0.0, FM=300e-12, dF=0.5e-12;
3830 k_func_t *k_func = bell_k;
3831 char *param_strings[] = {"3.3e-4", "0.25e-9"};
3832 void *params = string_create_bell_param_t(param_strings);
3833 environment_t env = {};
3836 for(F=Fm; F <= FM; F+=dF) {
3837 k = k_func(F, &env, params);
3838 k_accel = accel_k(k_func, F, &env, params);
3839 CHECK_ERR(1e-8, k, k_accel);
3841 destroy_bell_param_t(params);
3847 \section{Tension models}
3848 \label{sec.tension_models}
3850 TODO: tension model intro.
3851 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.
3853 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3854 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]].
3856 <<tension-model.h>>=
3857 #ifndef TENSION_MODEL_H
3858 #define TENSION_MODEL_H
3860 <<tension handler types>>
3861 <<constant tension model declarations>>
3862 <<hooke tension model declarations>>
3863 <<worm-like chain tension model declarations>>
3864 <<freely-jointed chain tension model declarations>>
3865 <<piston tension model declarations>>
3866 <<find tension definitions>>
3867 <<tension model global declarations>>
3868 #endif /* TENSION_MODEL_H */
3871 <<tension model module makefile lines>>=
3872 NW_SAWSIM_MODS += tension_model
3873 CHECK_BINS += check_tension_model
3874 check_tension_model_MODS = list tension_balance
3876 <<tension model utils makefile lines>>=
3877 TENSION_MODEL_MODS = tension_model parse list tension_balance
3878 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3879 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3880 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3881 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
3882 notangle -Rtension-model-utils.c $< > $@
3883 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
3884 $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3885 clean_tension_model_utils :
3886 rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3887 @ The pipe symbol [[|]] marks the prerequisits following it (in this
3888 case, the directories) as ‘order-only’ prerequisites. The timestamp
3889 on these prerequisits does not effect whether the rules are executed.
3890 This is appropriate for directories, where we don't need to recompile
3891 after an unrelated has been added to the directory, but only when the
3892 source prerequisites change. See the [[make]] documentation for more
3894 (\url{http://www.gnu.org/software/make/manual/html_node/Prerequisite-Types.html}).
3897 \label{sec.null_tension}
3899 For unstretchable domains.
3901 <<null tension model getopt>>=
3902 {NULL, NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3905 \subsection{Constant}
3906 \label{sec.const_tension}
3908 An infinitely stretchable domain providing a constant tension.
3909 Obviously this cannot be inverted, so you can't put this domain in
3910 series with any other domains. We include it mostly for testing
3913 <<constant tension functions>>=
3914 <<constant tension function>>
3915 <<constant tension parameter create/destroy functions>>
3918 <<constant tension model declarations>>=
3919 <<constant tension function declaration>>
3920 <<constant tension parameter create/destroy function declarations>>
3921 <<constant tension model global declarations>>
3925 <<constant tension function declaration>>=
3926 extern double const_tension_handler(double x, void *pdata);
3928 <<constant tension function>>=
3929 double const_tension_handler(double x, void *pdata)
3931 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3932 list_t *list = data->group_tension_params;
3937 assert(list != NULL); /* empty active group?! */
3938 F = ((const_tension_param_t *)list->d)->F;
3940 fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3942 while (list != NULL) {
3943 assert(((const_tension_param_t *)list->d)->F == F);
3951 <<constant tension structure>>=
3952 typedef struct const_tension_param_struct {
3953 double F; /* tension (force) in N */
3954 } const_tension_param_t;
3958 <<constant tension parameter create/destroy function declarations>>=
3959 extern void *string_create_const_tension_param_t(char **param_strings);
3960 extern void destroy_const_tension_param_t(void *p);
3963 <<constant tension parameter create/destroy functions>>=
3964 const_tension_param_t *create_const_tension_param_t(double F)
3966 const_tension_param_t *ret
3967 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3968 assert(ret != NULL);
3973 void *string_create_const_tension_param_t(char **param_strings)
3975 return create_const_tension_param_t(
3976 safe_strtod(param_strings[0],__FUNCTION__));
3979 void destroy_const_tension_param_t(void *p)
3987 <<constant tension model global declarations>>=
3988 extern int num_const_tension_params;
3989 extern char *const_tension_param_descriptions[];
3990 extern char const_tension_param_string[];
3993 <<constant tension model globals>>=
3994 int num_const_tension_params = 1;
3995 char *const_tension_param_descriptions[] = {"tension F, N"};
3996 char const_tension_param_string[] = "0";
3999 <<constant tension model getopt>>=
4000 {&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}
4006 <<hooke functions>>=
4007 <<internal hooke functions>>
4009 <<inverse hooke handler>>
4010 <<hooke parameter create/destroy functions>>
4013 <<hooke tension model declarations>>=
4014 <<hooke tension function declaration>>
4015 <<hooke tension parameter create/destroy function declarations>>
4016 <<hooke tension model global declarations>>
4020 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
4021 The behavior of a series of springs $k_i$ in series is given by
4023 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
4025 For a simple proof, see Appendix \ref{sec.math_hooke}.
4027 <<hooke tension function declaration>>=
4028 extern double hooke_handler(double x, void *pdata);
4029 extern double inverse_hooke_handler(double F, void *pdata);
4032 First we define a function that computes the effective spring constant
4033 (stored in a single [[hooke_param_t]]) for the entire group.
4034 <<internal hooke functions>>=
4035 static void hooke_param_grouper(tension_handler_data_t *data,
4036 hooke_param_t *param) {
4037 list_t *list = data->group_tension_params;
4041 assert(list != NULL); /* empty active group?! */
4042 while (list != NULL) {
4043 assert( ((hooke_param_t *)list->d)->k > 0 );
4044 k += 1.0/ ((hooke_param_t *)list->d)->k;
4046 fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
4047 __FUNCTION__, ((hooke_param_t *)list->d)->k);
4053 fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N (from data %p)\n",
4054 __FUNCTION__, k, x, k*x, data);
4061 Everyone knows Hooke's law: $F=kx$.
4063 double hooke_handler(double x, void *pdata)
4065 hooke_param_t param = {0};
4068 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
4074 Inverting Hooke's law gives $x=F/k$.
4075 <<inverse hooke handler>>=
4076 double inverse_hooke_handler(double F, void *pdata)
4078 hooke_param_t param = {0};
4081 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
4087 Not much to keep track of for a single effective spring, just the
4088 spring constant $k$.
4089 <<hooke structure>>=
4090 typedef struct hooke_param_struct {
4091 double k; /* spring constant in N/m */
4096 We wrap up our Hooke implementation with some book-keeping functions.
4097 <<hooke tension parameter create/destroy function declarations>>=
4098 extern void *string_create_hooke_param_t(char **param_strings);
4099 extern void destroy_hooke_param_t(void *p);
4103 <<hooke parameter create/destroy functions>>=
4104 hooke_param_t *create_hooke_param_t(double k)
4106 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
4107 assert(ret != NULL);
4112 void *string_create_hooke_param_t(char **param_strings)
4114 return create_hooke_param_t(safe_strtod(param_strings[0], __FUNCTION__));
4117 void destroy_hooke_param_t(void *p)
4124 <<hooke tension model global declarations>>=
4125 extern int num_hooke_params;
4126 extern char *hooke_param_descriptions[];
4127 extern char hooke_param_string[];
4130 <<hooke tension model globals>>=
4131 int num_hooke_params = 1;
4132 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
4133 char hooke_param_string[]="0.05";
4136 <<hooke tension model getopt>>=
4137 {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}
4140 \subsection{Worm-like chain}
4143 We can model several unfolded domains with a single worm-like chain.
4144 <<worm-like chain functions>>=
4145 <<internal worm-like chain functions>>
4146 <<worm-like chain handler>>
4147 <<inverse worm-like chain handler>>
4148 <<worm-like chain create/destroy functions>>
4151 <<worm-like chain tension model declarations>>=
4153 <<worm-like chain handler declaration>>
4154 <<inverse worm-like chain handler declaration>>
4155 <<worm-like chain create/destroy function declarations>>
4156 <<worm-like chain tension model global declarations>>
4160 <<internal worm-like chain functions>>=
4161 <<worm-like chain function>>
4162 <<inverse worm-like chain function>>
4163 <<worm-like chain parameter grouper>>
4166 The combination of all unfolded domains is modeled as a worm like chain,
4167 with the force $F_{\text{WLC}}$ approximately given by
4169 F_{\text{WLC}} \approx \frac{k_B T}{p}
4170 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
4172 where $x$ is the distance between the fixed ends,
4173 $k_B$ is Boltzmann's constant,
4174 $T$ is the absolute temperature,
4175 $p$ is the persistence length, and
4176 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
4177 <<worm-like chain function>>=
4178 static double wlc(double x, double T, double p, double L)
4180 double xL=0.0; /* uses global k_B */
4182 assert(T > 0); assert(p > 0); assert(L > 0);
4183 if (x >= L) return HUGE_VAL;
4184 xL = x/L; /* unitless */
4185 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
4186 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
4187 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
4192 Inverting the approximate WLC is fairly straightforward, if a bit tedious.
4193 Letting $x_L \equiv x/L$ and $F_T \equiv Fp/k_B T$ we have
4195 F_T &\approx \frac{1}{4} \p({\frac{1}{(1-x_L)^2} - 1}) + x_L \\
4196 F_T (1-x_L)^2 &\approx \frac{1}{4} - \frac{1}{4} (1-x_L)^2 + x_L(1-x_L)^2 \\
4197 0 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) (1-x_L)^2 + x_L(1-x_L)^2 \\
4198 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) \p({1-2x_L+x_L^2})
4199 + x_L - 2x_L^2 + x_L^3 \\
4200 &\approx \p[{\frac{1}{4} - \p({F_T + \frac{1}{4}})}]
4201 + x_L \p[{2F_T + \frac{1}{2} + 1}]
4202 + x_L^2 \p[{-\p({F_T + \frac{1}{4}}) - 2}]
4205 + x_L \p({2F_T + \frac{3}{2}})
4206 - x_L^2 \p({F_T + \frac{9}{4}})
4208 &\approx ax_L^3 + bx_L^2 + cx_L + d
4210 where $a \equiv 1$, $b \equiv -F_T - 9/4$, $c \equiv 2F_T + 3/2$, and
4212 % From \citet{wikipedia_cubic_function} the discriminant
4214 % \Delta &= -4b^3d + b^2c^2 - 4ac^3 + 18abcd - 27a^2d^2 \\
4215 % &= -4F_T\p({F_T + \frac{9}{4}})^3
4216 % + \p({F_T + \frac{9}{4}})^2 \cdot \p({2F_T + \frac{3}{2}})^2
4217 % - 4 \p({2F_T + \frac{3}{2}})^3
4218 % + 18F_T \p({F_T + \frac{9}{4}}) \cdot \p({2F_T + \frac{3}{2}})
4220 % &= -4F_T \p({F_T^3 + \frac{27}{4}F_T^2 +\frac{243}{16}F_T +\frac{729}{64}})
4221 %% a^3 + 3a^2b + 3ab^2 + b^3
4222 % + \p({F_T^2 + \frac{9}{2}F_T + \frac{81}{16}})
4223 % \cdot \p({4F_T^2 + 6F_T + \frac{9}{4}}) \\
4224 % &\qquad - 4 \p({8F_T^3 + 18F_T^2 + \frac{27}{2}F_T + \frac{27}{8}})
4225 %% a^3 + 3a^2b + 3ab^2 + b^3
4226 % + 18F_T \p({2F_T^2 + 6F_T + \frac{27}{8}})
4228 % &= -4F_T^4 - 27F_T^3 - \frac{243}{4}F_T^2 - \frac{729}{16}F_T
4229 % + 4F_T^4 + 24F_T^3 +\frac{99}{2}F_T^2 +\frac{81}{2}F_T +\frac{729}{64}\\
4230 % &\qquad - 32F_T^3 - 72F_T^2 - 54F_T - \frac{27}{2}
4231 % + 36F_T^3 + 108F_T^2 + \frac{243}{4}F_T
4233 % &= F_T^4 \p({-4 + 4}) + F_T^3 \p({-27 + 24 - 32 + 36})
4234 % + F_T^2 \p({-\frac{243}{4} + \frac{99}{2} - 72 + 108 -27}) \\
4235 % &\qquad + F_T \p({-\frac{729}{16} + \frac{81}{2} - 54 + \frac{243}{4}})
4236 % \p({\frac{729}{64} - \frac{27}{2}}) \\
4237 % &= F_T^3 - \frac{9}{4}F_T^2 + \frac{27}{16}F_T - \frac{135}{64}
4239 We can plug these coefficients into the GSL cubic solver to invert the
4240 WLC\citep{galassi05}. The applicable root is always the one which
4241 comes before the singularity, which will be the smallest real root.
4242 <<inverse worm-like chain function>>=
4243 static double inverse_wlc(double F, double T, double p, double L)
4245 double FT = F*p/(k_B*T); /* uses global k_B */
4246 double xL0, xL1, xL2;
4249 assert(T > 0); assert(p > 0); assert(L > 0);
4252 num_roots = gsl_poly_solve_cubic(-(FT+2.25),2*FT+1.5,-FT, &xL0, &xL1, &xL2);
4253 assert(num_roots > 0); assert(xL0 >= -DOUBLE_PRECISION); assert(xL0 < 1);
4254 if (xL0 < 0) xL0 = 0.0;
4260 First we define a function that computes the effective WLC parameters
4261 (stored in a single [[wlc_param_t]]) for the entire group.
4262 <<worm-like chain parameter grouper>>=
4263 static void wlc_param_grouper(tension_handler_data_t *data,
4264 wlc_param_t *param) {
4265 list_t *list = data->group_tension_params;
4266 double p=0.0, L=0.0;
4269 assert(list != NULL); /* empty active group?! */
4270 p = ((wlc_param_t *)list->d)->p;
4271 while (list != NULL) {
4272 /* enforce consistent p */
4273 assert( ((wlc_param_t *)list->d)->p == p);
4274 L += ((wlc_param_t *)list->d)->L;
4276 fprintf(stderr, "%s: WLC section %g m long in series\n",
4277 __FUNCTION__, ((wlc_param_t *)list->d)->L);
4282 fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
4283 __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
4291 <<worm-like chain handler declaration>>=
4292 extern double wlc_handler(double x, void *pdata);
4295 This model requires that each unfolded domain in the group have the
4296 same persistence length.
4297 <<worm-like chain handler>>=
4298 double wlc_handler(double x, void *pdata)
4300 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4301 wlc_param_t param = {0};
4304 wlc_param_grouper(data, ¶m);
4305 return wlc(x, data->env->T, param.p, param.L);
4310 <<inverse worm-like chain handler declaration>>=
4311 extern double inverse_wlc_handler(double F, void *pdata);
4314 <<inverse worm-like chain handler>>=
4315 double inverse_wlc_handler(double F, void *pdata)
4317 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4318 wlc_param_t param = {0};
4320 wlc_param_grouper(data, ¶m);
4321 return inverse_wlc(F, data->env->T, param.p, param.L);
4326 <<worm-like chain structure>>=
4327 typedef struct wlc_param_struct {
4328 double p; /* WLC persistence length (m) */
4329 double L; /* WLC contour length (m) */
4333 <<worm-like chain create/destroy function declarations>>=
4334 extern void *string_create_wlc_param_t(char **param_strings);
4335 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
4338 <<worm-like chain create/destroy functions>>=
4339 wlc_param_t *create_wlc_param_t(double p, double L)
4341 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
4342 assert(ret != NULL);
4348 void *string_create_wlc_param_t(char **param_strings)
4350 return create_wlc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4351 safe_strtod(param_strings[1], __FUNCTION__));
4354 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
4362 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.
4363 TODO (now we copy the string before parsing, but still do this for clarity).
4365 <<valid string write code>>=
4366 char string[] = "abc";
4369 <<valid string write code>>=
4370 char *string = "abc";
4374 <<worm-like chain tension model global declarations>>=
4375 extern int num_wlc_params;
4376 extern char *wlc_param_descriptions[];
4377 extern char wlc_param_string[];
4379 <<worm-like chain tension model globals>>=
4380 int num_wlc_params = 2;
4381 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
4382 char wlc_param_string[]="0.39e-9,28.4e-9";
4385 <<worm-like chain tension model getopt>>=
4386 {&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}
4388 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
4390 \subsection{Freely-jointed chain}
4393 <<freely-jointed chain functions>>=
4394 <<freely-jointed chain function>>
4395 <<freely-jointed chain handler>>
4396 <<freely-jointed chain create/destroy functions>>
4399 <<freely-jointed chain tension model declarations>>=
4400 <<freely-jointed chain handler declaration>>
4401 <<freely-jointed chain create/destroy function declarations>>
4402 <<freely-jointed chain tension model global declarations>>
4406 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
4407 The tension of a chain of $N$ such links is given by
4409 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
4411 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}.
4412 <<freely-jointed chain function>>=
4413 static double langevin(double x, void *params)
4416 return 1.0/tanh(x) - 1/x;
4419 static double inv_langevin(double x)
4421 double lb=0.0, ub=1.0, ret;
4422 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
4423 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
4427 static double fjc(double x, double T, double l, int N)
4429 double L = l*(double)N;
4430 if (x == 0) return 0;
4431 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
4432 return k_B*T/l * inv_langevin(x/L);
4437 <<freely-jointed chain handler declaration>>=
4438 extern double fjc_handler(double x, void *pdata);
4440 <<freely-jointed chain handler>>=
4441 double fjc_handler(double x, void *pdata)
4443 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4444 list_t *list = data->group_tension_params;
4449 assert(list != NULL); /* empty active group?! */
4450 l = ((fjc_param_t *)list->d)->l;
4451 while (list != NULL) {
4452 /* enforce consistent l */
4453 assert( ((fjc_param_t *)list->d)->l == l);
4454 N += ((fjc_param_t *)list->d)->N;
4456 fprintf(stderr, "%s: FJC section %d links long in series\n",
4457 __FUNCTION__, ((fjc_param_t *)list->d)->N);
4462 fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
4463 __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
4465 return fjc(x, data->env->T, l, N);
4469 <<freely-jointed chain structure>>=
4470 typedef struct fjc_param_struct {
4471 double l; /* FJC link length (m) */
4472 int N; /* FJC number of links */
4476 <<freely-jointed chain create/destroy function declarations>>=
4477 extern void *string_create_fjc_param_t(char **param_strings);
4478 extern void destroy_fjc_param_t(void *p);
4481 <<freely-jointed chain create/destroy functions>>=
4482 fjc_param_t *create_fjc_param_t(double l, double N)
4484 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
4485 assert(ret != NULL);
4491 void *string_create_fjc_param_t(char **param_strings)
4493 return create_fjc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4494 safe_strtod(param_strings[1], __FUNCTION__));
4497 void destroy_fjc_param_t(void *p)
4504 <<freely-jointed chain tension model global declarations>>=
4505 extern int num_fjc_params;
4506 extern char *fjc_param_descriptions[];
4507 extern char fjc_param_string[];
4510 <<freely-jointed chain tension model globals>>=
4511 int num_fjc_params = 2;
4512 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
4513 char fjc_param_string[]="0.5e-9,200";
4516 <<freely-jointed chain tension model getopt>>=
4517 {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}
4521 \label{sec.piston_tension}
4523 An infinitely stretchable domain with no tension for extensions less
4524 than a particular threshold $L$, and infinite tension for $x>L$. The
4525 tension at $x=L$ is undefined, but may be any positive value. The
4526 model is called the ``piston'' model because it resembles a
4527 frictionless piston sliding in a rigid cylinder of length $L$.
4530 0 & \text{if $x<L$}, \\
4531 \infty & \text{if $x>L$}.
4535 Note that because of it's infinte stiffness at $x=L$, fully extended
4536 domains with this tension model will be identical to completely rigid
4537 domains. However, there is a distinction between this model and rigid
4538 domains for $x<L$, because any reactions that occur at $F=0$
4539 (e.g. refolding) will have more time to occur.
4541 Because the tension at $x=L$ is undefined, the user must make sure
4542 domains with this tension model are never the initial active group,
4543 because it would break [[tension_balance()]] in [[find_tension()]]
4544 (see Section \ref{sec.tension_balance}).
4546 <<piston tension functions>>=
4547 <<piston tension parameter grouper>>
4548 <<piston tension handler>>
4549 <<inverse piston tension handler>>
4550 <<piston tension parameter create/destroy functions>>
4553 <<piston tension model declarations>>=
4554 <<piston tension handler declaration>>
4555 <<inverse piston tension handler declaration>>
4556 <<piston tension parameter create/destroy function declarations>>
4557 <<piston tension model global declarations>>
4561 First we define a function that computes the effective piston parameters
4562 (stored in a single [[piston_tension_param_t]]) for the entire group.
4563 <<piston tension parameter grouper>>=
4564 static void piston_tension_param_grouper(tension_handler_data_t *data,
4565 piston_tension_param_t *param) {
4566 list_t *list = data->group_tension_params;
4570 assert(list != NULL); /* empty active group?! */
4571 while (list != NULL) {
4572 L += ((piston_tension_param_t *)list->d)->L;
4578 <<piston tension handler declaration>>=
4579 extern double piston_tension_handler(double x, void *pdata);
4581 <<piston tension handler>>=
4582 double piston_tension_handler(double x, void *pdata)
4584 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4585 piston_tension_param_t param = {0};
4588 piston_tension_param_grouper(data, ¶m);
4589 if (x <= param.L) F = 0;
4592 fprintf(stderr, "%s: x %g, L %g, tension %g N\n", __FUNCTION__, x, param.L, F);
4599 We abuse [[x_of_xo()]] and [[oneD_solve()]] with this inverse
4600 definition (see Section \ref{sec.tension_balance}), because the
4601 $x(F=0)<=L$, but our function returns $x(F=0)=0$. This is fine when
4602 the total extension \emph{is} zero, but cheats a bit when there is
4603 some total extension. For any non-zero extension, the initial active
4604 group produces some tension (we hope. True for all our other tension
4605 models). This tension fully extends the piston group (of which there
4606 should be only one, since all piston states can get lumped into the
4607 same tension group.). If the total extension is $\le$ the target
4608 extension, the full extension is accurate, so the inverse was valid
4609 after all. If not, [[oneD_solve()]] will halve the extension of the
4610 first group, reducing the overall tension. For total extension less
4611 than the piston extension, this bisection continues for [[max_steps]],
4612 leaving $x_0$ reduced by a factor of $2^\text{max\_steps}$, a very
4613 small number. Because of this, the returned tension $F_0(x_0)$ is
4614 very small, even though the total extension found by [[x_of_xo()]]
4615 is still $>$ the piston length.
4617 While this works, it is clearly not the most elegant or robust
4618 solution. TODO: think of (and implement) a better procedure.
4620 <<inverse piston tension handler declaration>>=
4621 extern double inverse_piston_tension_handler(double x, void *pdata);
4623 <<inverse piston tension handler>>=
4624 double inverse_piston_tension_handler(double F, void *pdata)
4626 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4627 piston_tension_param_t param = {0};
4630 piston_tension_param_grouper(data, ¶m);
4631 if (F == 0.0) return 0;
4637 <<piston tension structure>>=
4638 typedef struct piston_tension_param_struct {
4639 double L; /* length of piston in m */
4640 } piston_tension_param_t;
4644 <<piston tension parameter create/destroy function declarations>>=
4645 extern void *string_create_piston_tension_param_t(char **param_strings);
4646 extern void destroy_piston_tension_param_t(void *p);
4650 <<piston tension parameter create/destroy functions>>=
4651 piston_tension_param_t *create_piston_tension_param_t(double L)
4653 piston_tension_param_t *ret
4654 = (piston_tension_param_t *) malloc(sizeof(piston_tension_param_t));
4655 assert(ret != NULL);
4660 void *string_create_piston_tension_param_t(char **param_strings)
4662 return create_piston_tension_param_t(
4663 safe_strtod(param_strings[0],__FUNCTION__));
4666 void destroy_piston_tension_param_t(void *p)
4674 <<piston tension model global declarations>>=
4675 extern int num_piston_tension_params;
4676 extern char *piston_tension_param_descriptions[];
4677 extern char piston_tension_param_string[];
4681 <<piston tension model globals>>=
4682 int num_piston_tension_params = 1;
4683 char *piston_tension_param_descriptions[] = {"length L, m"};
4684 char piston_tension_param_string[] = "0";
4688 <<piston tension model getopt>>=
4689 {&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}
4692 \subsection{Tension model implementation}
4694 <<tension-model.c>>=
4697 <<tension model includes>>
4698 #include "tension_model.h"
4699 <<tension model internal definitions>>
4700 <<tension model globals>>
4701 <<tension model functions>>
4704 <<tension model includes>>=
4705 #include <assert.h> /* assert() */
4706 #include <stdlib.h> /* NULL, strto*() */
4707 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4708 #include <stdio.h> /* fprintf(), stdout */
4709 #include <errno.h> /* errno, ERANGE */
4712 #include "tension_balance.h" /* oneD_*() */
4715 <<tension model internal definitions>>=
4716 <<constant tension structure>>
4718 <<worm-like chain structure>>
4719 <<freely-jointed chain structure>>
4720 <<piston tension structure>>
4723 <<tension model globals>>=
4724 <<tension handler array global>>
4725 <<constant tension model globals>>
4726 <<hooke tension model globals>>
4727 <<worm-like chain tension model globals>>
4728 <<freely-jointed chain tension model globals>>
4729 <<piston tension model globals>>
4732 <<tension model functions>>=
4734 <<constant tension functions>>
4736 <<worm-like chain functions>>
4737 <<freely-jointed chain functions>>
4738 <<piston tension functions>>
4741 \subsection{Utilities}
4743 The tension models can get complicated, and users may want to reassure
4744 themselves that this portion of the input physics are functioning properly. The
4745 stand-alone program [[tension_model_utils]] demonstrates the tension model
4746 interface without the overhead of having to understand a full simulation.
4747 [[tension_model_utils]] takes command line model arguments like the full
4748 simulation, and prints $F(x)$ data to the screen for the selected model over a
4751 <<tension-model-utils.c>>=
4753 <<tension model utility includes>>
4754 <<tension model utility definitions>>
4755 <<tension model utility getopt functions>>
4757 <<tension model utility main>>
4760 <<tension model utility main>>=
4761 int main(int argc, char **argv)
4763 <<tension model getopt array definition>>
4764 tension_model_getopt_t *model = NULL;
4768 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
4769 tension_handler_data_t tdata;
4770 int num_param_args; /* for INIT_MODEL() */
4771 char **param_args; /* for INIT_MODEL() */
4773 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model,
4774 &Fmax, &Xmax, &flags);
4776 if (flags & VFLAG) {
4777 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4779 INIT_MODEL("utility", model, model->params, params);
4780 tdata.group_tension_params = NULL;
4781 push(&tdata.group_tension_params, params, 1);
4783 tdata.persist = NULL;
4784 if (model->handler == NULL) {
4785 printf("No tension function for model %s\n", model->name);
4791 printf("#Distance (m)\tForce (N)\n");
4792 for (i=0; i<=N; i++) {
4793 x = Xmax*i/(double)N;
4794 F = (*model->handler)(x, &tdata);
4795 if (F < 0 || F > Fmax) break;
4796 printf("%g\t%g\n", x, F);
4798 if (flags & VFLAG && i <= N) {
4799 /* explain exit condition */
4801 printf("Impossible force %g\n", F);
4803 printf("Reached large force limit %g > %g\n", F, Fmax);
4806 params = pop(&tdata.group_tension_params);
4807 if (model->destructor)
4808 (*model->destructor)(params);
4813 <<tension model utility includes>>=
4816 #include <string.h> /* strlen() */
4817 #include <assert.h> /* assert() */
4818 #include <errno.h> /* errno, ERANGE */
4822 #include "tension_model.h"
4825 <<tension model utility definitions>>=
4826 <<version definition>>
4827 #define VFLAG 1 /* verbose */
4828 <<string comparison definition and globals>>
4829 <<tension model getopt definitions>>
4830 <<initialize model definition>>
4834 <<tension model utility getopt functions>>=
4836 <<index tension model>>
4837 <<tension model utility help>>
4838 <<tension model utility get options>>
4841 <<tension model utility help>>=
4842 void help(char *prog_name,
4844 int n_tension_models, tension_model_getopt_t *tension_models,
4845 int tension_model, double Fmax, double Xmax)
4848 printf("usage: %s [options]\n", prog_name);
4849 printf("Version %s\n\n", VERSION);
4850 printf("A utility for understanding the available tension models\n\n");
4851 printf("Environmental options:\n");
4852 printf("-T T\tTemperature T (currently %g K)\n", env->T);
4853 printf("-C T\tYou can also set the temperature T in Celsius\n");
4854 printf("Model options:\n");
4855 printf("-m model\tFolded tension model (currently %s)\n",
4856 tension_models[tension_model].name);
4857 printf("-a args\tFolded tension model argument string (currently %s)\n",
4858 tension_models[tension_model].params);
4859 printf("F(x) is calculated for a range of x and printed\n");
4860 printf("For example:\n");
4861 printf(" #Distance (m)\tForce (N)\n");
4862 printf(" 123.456\t7.89\n");
4864 printf("-F\tSet the maximum F value for the standard mode F(x) output (currently %g)\n", Fmax);
4865 printf("-X\tSet the maximum x value for the standard mode F(x) output (currently %g)\n", Xmax);
4866 printf("-V\tChange output to verbose mode\n");
4867 printf("-h\tPrint this help and exit\n");
4869 printf("Tension models:\n");
4870 for (i=0; i<n_tension_models; i++) {
4871 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
4872 for (j=0; j < tension_models[i].num_params ; j++)
4873 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
4874 printf("\t\tdefaults: %s\n", tension_models[i].params);
4879 <<tension model utility get options>>=
4880 void get_options(int argc, char **argv, environment_t *env,
4881 int n_tension_models, tension_model_getopt_t *tension_models,
4882 tension_model_getopt_t **model, double *Fmax, double *Xmax,
4883 unsigned int *flags)
4885 char *prog_name = NULL;
4886 char c, options[] = "T:C:m:a:F:X:Vh";
4887 int tension_model=0, num_strings;
4888 extern char *optarg;
4889 extern int optind, optopt, opterr;
4893 /* setup defaults */
4895 prog_name = argv[0];
4896 env->T = 300.0; /* K */
4897 *Fmax = 1e5; /* N */
4898 *Xmax = 1e-6; /* m */
4900 *model = tension_models;
4902 while ((c=getopt(argc, argv, options)) != -1) {
4904 case 'T': env->T = safe_strtod(optarg, "-T"); break;
4905 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
4907 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
4908 *model = tension_models+tension_model;
4911 tension_models[tension_model].params = optarg;
4913 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
4914 case 'X': *Xmax = safe_strtod(optarg, "-X"); break;
4915 case 'V': *flags |= VFLAG; break;
4917 fprintf(stderr, "unrecognized option '%c'\n", optopt);
4918 /* fall through to default case */
4920 help(prog_name, env, n_tension_models, tension_models, tension_model, *Fmax, *Xmax);
4928 \subsection{Tension model unit tests}
4930 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
4931 <<check-tension-model.c>>=
4933 <<tension model check includes>>
4934 <<check relative error>>
4936 <<tension model test suite>>
4937 <<main check program>>
4940 <<tension model check includes>>=
4941 #include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
4942 #include <stdio.h> /* sprintf() */
4943 #include <assert.h> /* assert() */
4944 #include <math.h> /* exp() */
4945 #include <errno.h> /* errno, ERANGE */
4949 #include "tension_balance.h" /* oneD_*() */
4950 #include "tension_model.h"
4953 <<tension model test suite>>=
4955 <<const tension tests>>
4957 <<worm-like chain tests>>
4958 <<freely-jointed chain tests>>
4960 <<tension model suite definition>>
4963 <<tension model suite definition>>=
4964 Suite *test_suite (void)
4966 Suite *s = suite_create ("tension model");
4967 <<const tension test case defs>>
4968 <<hooke test case defs>>
4969 <<worm-like chain test case defs>>
4970 <<freely-jointed chain test case defs>>
4971 <<piston test case defs>>
4973 <<const tension test case adds>>
4974 <<hooke test case adds>>
4975 <<worm-like chain test case adds>>
4976 <<freely-jointed chain test case adds>>
4977 <<piston test case adds>>
4982 \subsubsection{Constant}
4984 <<const tension test case defs>>=
4985 TCase *tc_const_tension = tcase_create("Constant tension");
4988 <<const tension test case adds>>=
4989 tcase_add_test(tc_const_tension, test_const_tension_create_destroy);
4990 suite_add_tcase(s, tc_const_tension);
4993 <<const tension tests>>=
4994 START_TEST(test_const_tension_create_destroy)
4996 char *tension[] = {"1","2.2", "3"};
4997 char *params[] = {tension[0]};
5000 for( i=0; i < sizeof(tension)/sizeof(char *); i++) {
5001 params[0] = tension[i];
5002 p = string_create_const_tension_param_t(params);
5003 destroy_const_tension_param_t(p);
5008 @ TODO: In order to test the constant tension handler itself, we'd
5009 have to construct a group.
5012 \subsubsection{Hooke}
5014 <<hooke test case defs>>=
5015 TCase *tc_hooke = tcase_create("Hooke");
5018 <<hooke test case adds>>=
5019 tcase_add_test(tc_hooke, test_hooke_create_destroy);
5020 suite_add_tcase(s, tc_hooke);
5025 START_TEST(test_hooke_create_destroy)
5027 char *k[] = {"1","2.2", "3"};
5028 char *params[] = {k[0]};
5031 for( i=0; i < sizeof(k)/sizeof(char *); i++) {
5033 p = string_create_hooke_param_t(params);
5034 destroy_hooke_param_t(p);
5039 @ TODO: In order to test the Hooke tension handler itself, we'd
5040 have to construct a group.
5043 \subsubsection{Worm-like chain}
5045 <<worm-like chain test case defs>>=
5046 TCase *tc_wlc = tcase_create("WLC");
5049 <<worm-like chain test case adds>>=
5050 tcase_add_test(tc_wlc, test_wlc_at_zero);
5051 tcase_add_test(tc_wlc, test_wlc_at_half);
5052 suite_add_tcase(s, tc_wlc);
5056 <<worm-like chain tests>>=
5057 <<worm-like chain function>>
5058 START_TEST(test_wlc_at_zero)
5060 double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
5061 fail_unless(abs(wlc(x, T, p, L)) < lim, \
5062 "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
5063 x, T, p, L, abs(wlc(x, T, p, L)), lim);
5067 START_TEST(test_wlc_at_half)
5069 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
5070 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
5071 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
5073 * linear term = x/L = 0.5
5074 * nonlinear + linear = 0.75 + 0.5 = 1.25
5075 * wlc = 10e21*1.25 = 12.5e21
5077 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
5078 "wlc(%g, %g, %g, %g) = %g != %g",
5079 x, T, p, L, wlc(x, T, p, L), 12.5e21);
5086 \subsubsection{Freely-jointed chain}
5088 <<freely-jointed chain test case defs>>=
5089 TCase *tc_fjc = tcase_create("FJC");
5092 <<freely-jointed chain test case adds>>=
5093 tcase_add_test(tc_fjc, test_fjc_at_zero);
5094 tcase_add_test(tc_fjc, test_fjc_at_half);
5095 suite_add_tcase(s, tc_fjc);
5099 <<freely-jointed chain tests>>=
5100 <<freely-jointed chain function>>
5101 START_TEST(test_fjc_at_zero)
5104 double T=1.0, l=1.0, p=0.1, x=0.0, lim=1e-30;
5105 fail_unless(abs(fjc(x, T, l, N)) < lim, \
5106 "abs(fjc(x=%g,T=%g,l=%g,N=%d)) = %g >= %g",
5107 x, T, l, N, abs(fjc(x, T, l, N)), lim);
5111 START_TEST(test_fjc_at_half)
5114 double T=1.0/k_B, l=1.0, x=5.0;
5115 /* prefactor = k_B T / l = k_B (1.0/k_B) / 1.0 = 1.0 J/nm = 1.0e21 pN
5116 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
5118 * linear term = x/L = 0.5
5119 * nonlinear + linear = 0.75 + 0.5 = 1.25
5120 * fjc = 10e21*1.25 = 12.5e21
5122 fail_unless(fjc(x, T, l, N)-12.5e21 < 1e16,
5123 "fjc(%g, %g, %g, %d) = %g != %g",
5124 x, T, l, N, fjc(x, T, l, N), 12.5e21);
5131 \subsubsection{Piston}
5133 <<piston test case defs>>=
5134 TCase *tc_piston = tcase_create("Piston");
5137 <<piston test case adds>>=
5138 tcase_add_test(tc_piston, test_piston_create_destroy);
5139 suite_add_tcase(s, tc_piston);
5144 START_TEST(test_piston_create_destroy)
5146 char *L[] = {"1","2.2", "3"};
5147 char *params[] = {L[0]};
5150 for( i=0; i < sizeof(L)/sizeof(char *); i++) {
5152 p = string_create_piston_tension_param_t(params);
5153 destroy_piston_tension_param_t(p);
5158 @ TODO: In order to test the piston tension handler itself, we'd
5159 have to construct a group.
5162 \section{Unfolding rate models}
5163 \label{sec.k_models}
5165 We have two main choices for determining $k$: Bell's model, which assumes the
5166 folded domains exist in a single `folded' state and make instantaneous,
5167 stochastic transitions over a free energy barrier; and Kramers' model, which
5168 treats the unfolding as a thermalized diffusion process.
5169 We deal with these options like we dealt with the different tension models: by bundling all model-specific
5170 parameters into structures, with model specific functions for determining $k$.
5172 <<k func definition>>=
5173 typedef double k_func_t(double F, environment_t *env, void *params);
5176 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
5177 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]].
5179 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
5180 because \LaTeX\ doesn't like underscores outside math mode.
5181 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
5182 You could use spaces instead (`\verb+ +'),
5183 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
5184 which I don't like as much.
5190 <<k func definition>>
5191 <<null k declarations>>
5192 <<const k declarations>>
5193 <<bell k declarations>>
5194 <<kbell k declarations>>
5195 <<kramers k declarations>>
5196 <<kramers integ k declarations>>
5197 #endif /* K_MODEL_H */
5200 <<k model module makefile lines>>=
5201 NW_SAWSIM_MODS += k_model
5202 CHECK_BINS += check_k_model
5203 check_k_model_MODS = parse string_eval
5205 [[check_k_model]] also depends on the parse module.
5207 <<k model utils makefile lines>>=
5208 K_MODEL_MODS = k_model parse string_eval
5209 K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
5210 $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h)
5211 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
5212 notangle -Rk-model-utils.c $< > $@
5213 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
5214 $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5215 clean_k_model_utils :
5216 rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
5222 For domains that are never folded.
5224 <<null k declarations>>=
5226 <<null k functions>>=
5231 <<null k model getopt>>=
5232 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
5235 \subsection{Constant rate model}
5238 This is a pretty boring model, but it is usefull for testing the engine.
5242 where $k_0$ is the constant unfolding rate.
5244 <<const k functions>>=
5245 <<const k function>>
5246 <<const k structure create/destroy functions>>
5249 <<const k declarations>>=
5250 <<const k function declaration>>
5251 <<const k structure create/destroy function declarations>>
5252 <<const k global declarations>>
5256 <<const k structure>>=
5257 typedef struct const_k_param_struct {
5258 double knot; /* transition rate k_0 (frac population per s) */
5262 <<const k function declaration>>=
5263 double const_k(double F, environment_t *env, void *const_k_params);
5265 <<const k function>>=
5266 double const_k(double F, environment_t *env, void *const_k_params)
5267 { /* Returns the rate constant k in frac pop per s. */
5268 const_k_param_t *p = (const_k_param_t *) const_k_params;
5270 assert(p->knot > 0);
5275 <<const k structure create/destroy function declarations>>=
5276 void *string_create_const_k_param_t(char **param_strings);
5277 void destroy_const_k_param_t(void *p);
5280 <<const k structure create/destroy functions>>=
5281 const_k_param_t *create_const_k_param_t(double knot)
5283 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
5284 assert(ret != NULL);
5289 void *string_create_const_k_param_t(char **param_strings)
5291 return create_const_k_param_t(safe_strtod(param_strings[0], __FUNCTION__));
5294 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
5301 <<const k global declarations>>=
5302 extern int num_const_k_params;
5303 extern char *const_k_param_descriptions[];
5304 extern char const_k_param_string[];
5307 <<const k globals>>=
5308 int num_const_k_params = 1;
5309 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
5310 char const_k_param_string[]="1";
5313 <<const k model getopt>>=
5314 {"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}
5317 \subsection{Bell's model}
5320 Quantitatively, Bell's model gives $k$ as
5322 k = k_0 \cdot e^{F dx / k_B T} \;,
5324 where $k_0$ is the unforced unfolding rate,
5325 $F$ is the applied unfolding force,
5326 $dx$ is the distance to the transition state, and
5327 $k_B T$ is the thermal energy\citep{schlierf06}.
5329 <<bell k functions>>=
5331 <<bell k structure create/destroy functions>>
5334 <<bell k declarations>>=
5335 <<bell k function declaration>>
5336 <<bell k structure create/destroy function declarations>>
5337 <<bell k global declarations>>
5341 <<bell k structure>>=
5342 typedef struct bell_param_struct {
5343 double knot; /* transition rate k_0 (frac population per s) */
5344 double dx; /* `distance to transition state' (nm) */
5348 <<bell k function declaration>>=
5349 double bell_k(double F, environment_t *env, void *bell_params);
5351 <<bell k function>>=
5352 double bell_k(double F, environment_t *env, void *bell_params)
5353 { /* Returns the rate constant k in frac pop per s.
5354 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5355 * uses global k_B in J/K */
5356 bell_param_t *p = (bell_param_t *) bell_params;
5357 assert(F >= 0); assert(env->T > 0);
5359 assert(p->knot > 0); //assert(p->dx >= 0); /* p < 0 for refolding reactions, etc. */
5361 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
5362 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
5363 p->knot * exp(F*p->dx / (k_B*env->T) ));
5364 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
5366 return p->knot * exp(F*p->dx / (k_B*env->T) );
5370 <<bell k structure create/destroy function declarations>>=
5371 void *string_create_bell_param_t(char **param_strings);
5372 void destroy_bell_param_t(void *p);
5375 <<bell k structure create/destroy functions>>=
5376 bell_param_t *create_bell_param_t(double knot, double dx)
5378 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
5379 assert(ret != NULL);
5385 void *string_create_bell_param_t(char **param_strings)
5387 return create_bell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5388 safe_strtod(param_strings[1], __FUNCTION__));
5391 void destroy_bell_param_t(void *p /* bell_param_t * */)
5398 <<bell k global declarations>>=
5399 extern int num_bell_params;
5400 extern char *bell_param_descriptions[];
5401 extern char bell_param_string[];
5405 int num_bell_params = 2;
5406 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
5407 char bell_param_string[]="3.3e-4,0.25e-9";
5410 <<bell k model getopt>>=
5411 {"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}
5413 Initialized with Titin I27 parameters\citep{carrion-vazquez99a}.
5414 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5417 \subsection{Linker stiffness corrected Bell model}
5420 Walton et. al showed that the Bell model constant force approximation
5421 doesn't hold for biotin-streptavadin binding, and I extended their
5422 results to I27 for the 2009 Biophysical Society Annual
5423 Meeting\cite{walton08,king09}. More details to follow when I get done
5424 with the conference\ldots
5426 We adjust Bell's model to
5428 k = k_0 \cdot e^{\frac{F dx - \frac{1}{2} \kappa dx^2}{k_B T}} \;,
5430 where $k_0$ is the unforced unfolding rate,
5431 $F$ is the applied unfolding force,
5432 $\kappa$ is the effective spring constant,
5433 $dx$ is the distance to the transition state, and
5434 $k_B T$ is the thermal energy\citep{schlierf06}.
5436 <<kbell k functions>>=
5437 <<kbell k function>>
5438 <<kbell k structure create/destroy functions>>
5441 <<kbell k declarations>>=
5442 <<kbell k function declaration>>
5443 <<kbell k structure create/destroy function declarations>>
5444 <<kbell k global declarations>>
5448 <<kbell k structure>>=
5449 typedef struct kbell_param_struct {
5450 double knot; /* transition rate k_0 (frac population per s) */
5451 double dx; /* `distance to transition state' (nm) */
5455 <<kbell k function declaration>>=
5456 double kbell_k(double F, environment_t *env, void *kbell_params);
5458 <<kbell k function>>=
5459 double kbell_k(double F, environment_t *env, void *kbell_params)
5460 { /* Returns the rate constant k in frac pop per s.
5461 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5462 * uses global k_B in J/K */
5463 kbell_param_t *p = (kbell_param_t *) kbell_params;
5464 assert(F >= 0); assert(env->T > 0); assert(env->stiffness >= 0);
5466 assert(p->knot > 0); assert(p->dx >= 0);
5467 return p->knot * exp((F*p->dx - 0.5*env->stiffness*(p->dx)*(p->dx))/(k_B*env->T) );
5471 <<kbell k structure create/destroy function declarations>>=
5472 void *string_create_kbell_param_t(char **param_strings);
5473 void destroy_kbell_param_t(void *p);
5476 <<kbell k structure create/destroy functions>>=
5477 kbell_param_t *create_kbell_param_t(double knot, double dx)
5479 kbell_param_t *ret = (kbell_param_t *) malloc(sizeof(kbell_param_t));
5480 assert(ret != NULL);
5486 void *string_create_kbell_param_t(char **param_strings)
5488 return create_kbell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5489 safe_strtod(param_strings[1], __FUNCTION__));
5492 void destroy_kbell_param_t(void *p /* kbell_param_t * */)
5499 <<kbell k global declarations>>=
5500 extern int num_kbell_params;
5501 extern char *kbell_param_descriptions[];
5502 extern char kbell_param_string[];
5505 <<kbell k globals>>=
5506 int num_kbell_params = 2;
5507 char *kbell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
5508 char kbell_param_string[]="3.3e-4,0.25e-9";
5511 <<kbell k model getopt>>=
5512 {"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}
5514 Initialized with Titin I27 parameters\citep{carrion-vazquez99a}.
5515 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5518 \subsection{Kramer's model}
5521 Kramer's model gives $k$ as
5523 % \frac{1}{k} = \frac{1}{D}
5524 % \int_{x_\text{min}}^{x_\text{max}}
5525 % e^{\frac{-U_F(x)}{k_B T}}
5526 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5529 %where $D$ is the diffusion constant,
5530 %$U_F(x)$ is the free energy along the unfolding coordinate, and
5531 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
5532 % before the minimum and after the maximum, respectively \citep{schlierf06}.
5534 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
5536 where $D$ is the diffusion constant,
5538 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
5540 are length scales where
5541 $x_c(F)$ is the location of the bound state and
5542 $x_{ts}(F)$ is the location of the transition state,
5543 $E(F,x)$ is the free energy, and
5544 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
5545 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
5546 \citet{evans97} Eqn.~3).
5548 <<kramers k functions>>=
5549 <<kramers k function>>
5550 <<kramers k structure create/destroy functions>>
5553 <<kramers k declarations>>=
5554 <<kramers k function declaration>>
5555 <<kramers k structure create/destroy function declarations>>
5556 <<kramers k global declarations>>
5560 <<kramers k structure>>=
5561 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
5562 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
5564 typedef struct kramers_param_struct {
5565 double D; /* diffusion rate (um/s) */
5572 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
5573 //kramers_x_func_t *xts; /* function returning transition x (nm) */
5574 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
5575 //kramers_E_func_t *E; /* function returning E (J) */
5576 //double *E_params; /* parametrize the energy landscape */
5577 //int n_E_params; /* length of E_params array */
5581 <<kramers k function declaration>>=
5582 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
5583 extern double kramers_k(double F, environment_t *env, void *kramers_params);
5585 <<kramers k function>>=
5586 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
5588 static char input[160]={0};
5589 static char output[80]={0};
5590 /* setup the environment */
5591 fprintf(in, "F = %g; T = %g\n", F, T);
5592 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
5593 string_eval(in, out, input, 80, output);
5594 return safe_strtod(output, __FUNCTION__);
5597 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
5599 static char input[160]={0};
5600 static char output[80]={0};
5601 /* setup the environment */
5602 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
5603 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
5604 string_eval(in, out, input, 80, output);
5605 return safe_strtod(output, __FUNCTION__);
5608 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
5610 kramers_param_t *p = (kramers_param_t *) kramers_params;
5611 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
5614 double kramers_k(double F, environment_t *env, void *kramers_params)
5615 { /* Returns the rate constant k in frac pop per s.
5616 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5617 * uses global k_B in J/K */
5618 kramers_param_t *p = (kramers_param_t *) kramers_params;
5619 double kbT, xc, xts, lc, lts, Eb;
5620 assert(F >= 0); assert(env->T > 0);
5623 assert(p->xc != NULL); assert(p->xts != NULL);
5624 assert(p->ddEdxx != NULL); assert(p->E != NULL);
5625 assert(p->in != NULL); assert(p->out != NULL);
5627 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
5628 if (xc == -1.0) return -1.0; /* error (Too much force) */
5629 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
5630 if (xc == -1.0) return -1.0; /* error (Too much force) */
5631 lc = sqrt(2.0*M_PI*kbT /
5632 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
5633 lts = sqrt(-2.0*M_PI*kbT/
5634 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
5635 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
5636 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
5637 //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));
5638 return p->D/(lc*lts) * exp(-Eb/kbT);
5642 <<kramers k structure create/destroy function declarations>>=
5643 void *string_create_kramers_param_t(char **param_strings);
5644 void destroy_kramers_param_t(void *p);
5647 <<kramers k structure create/destroy functions>>=
5648 kramers_param_t *create_kramers_param_t(double D,
5649 char *xc_fn, char *xts_fn,
5650 char *E_fn, char *ddEdxx_fn,
5651 char *global_define_string)
5652 // kramers_x_func_t *xc, kramers_x_func_t *xts,
5653 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
5654 // double *E_params, int n_E_params)
5656 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
5657 assert(ret != NULL);
5658 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5659 assert(ret->xc != NULL);
5660 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5661 assert(ret->xts != NULL);
5662 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
5663 assert(ret->E != NULL);
5664 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
5665 assert(ret->ddEdxx != NULL);
5667 strcpy(ret->xc, xc_fn);
5668 strcpy(ret->xts, xts_fn);
5669 strcpy(ret->E, E_fn);
5670 strcpy(ret->ddEdxx, ddEdxx_fn);
5671 //ret->E_params = E_params;
5672 //ret->n_E_params = n_E_params;
5673 string_eval_setup(&ret->in, &ret->out);
5674 fprintf(ret->in, "kB = %g\n", k_B);
5675 fprintf(ret->in, "%s\n", global_define_string);
5679 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
5680 void *string_create_kramers_param_t(char **param_strings)
5682 return create_kramers_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5690 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
5692 kramers_param_t *pk = (kramers_param_t *)p;
5694 string_eval_teardown(&pk->in, &pk->out);
5696 // free(pk->E_params);
5702 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
5703 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.
5704 We follow \citet{shillcock98} and use
5706 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
5707 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
5710 where TODO, variable meanings.%\citep{shillcock98}.
5711 The first and second derivatives of $E(F,E_b,x,x_s)$ are
5713 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\\
5714 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
5716 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
5717 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
5718 The roots are therefor at
5720 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
5721 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
5722 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
5725 <<kramers k global declarations>>=
5726 extern int num_kramers_params;
5727 extern char *kramers_param_descriptions[];
5728 extern char kramers_param_string[];
5731 <<kramers k globals>>=
5732 int num_kramers_params = 6;
5733 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"};
5734 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)}";
5736 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
5737 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
5738 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
5739 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.
5740 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
5741 It works so long as [[val_1]] is not [[false]].
5743 <<kramers k model getopt>>=
5744 {"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}
5747 \subsection{Kramer's model (natural cubic splines)}
5748 \label{sec.kramers_integ}
5750 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
5751 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
5752 \citet{schlierf06} Eqn.~2)
5754 \frac{1}{k} = \frac{1}{D}
5755 \int_{x_\text{min}}^{x_\text{max}}
5756 e^{\frac{U_F(x)}{k_B T}}
5757 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5760 where $D$ is the diffusion constant,
5761 $U_F(x)$ is the free energy along the unfolding coordinate, and
5762 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
5763 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
5765 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
5766 interpolating between them with cubic splines.
5767 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
5769 <<kramers integ k functions>>=
5770 <<spline functions>>
5771 <<kramers integ A k functions>>
5772 <<kramers integ B k functions>>
5773 <<kramers integ k function>>
5774 <<kramers integ k structure create/destroy functions>>
5777 <<kramers integ k declarations>>=
5778 <<kramers integ k function declaration>>
5779 <<kramers integ k structure create/destroy function declarations>>
5780 <<kramers integ k global declarations>>
5784 <<kramers integ k structure>>=
5785 typedef double func_t(double x, void *params);
5786 typedef struct kramers_integ_param_struct {
5787 double D; /* diffusion rate (um/s) */
5788 double x_min; /* integration bounds */
5790 func_t *E; /* function returning E (J) */
5791 void *E_params; /* parametrize the energy landscape */
5792 destroy_data_func_t *destroy_E_params;
5794 double F; /* for passing into GSL evaluations */
5796 } kramers_integ_param_t;
5799 <<spline param structure>>=
5800 typedef struct spline_param_struct {
5801 int n_knots; /* natural cubic spline knots for U(x) */
5804 gsl_spline *spline; /* GSL spline parameters */
5805 gsl_interp_accel *acc; /* GSL spline acceleration data */
5809 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
5811 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5813 <<kramers integ A k functions>>=
5814 double kramers_integ_k_integrandA(double x, void *params)
5816 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5817 double U = (*p->E)(x, p->E_params);
5818 double Fx = p->F * x;
5819 double kbT = k_B * p->env->T;
5820 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
5821 return exp(-(U-Fx)/kbT);
5823 double kramers_integ_k_integralA(double x, void *params)
5825 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5827 double result, abserr;
5828 size_t neval; /* for qng */
5829 static gsl_integration_workspace *w = NULL;
5830 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5831 f.function = &kramers_integ_k_integrandA;
5833 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5834 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5835 //fprintf(stderr, "integralA = %g\n", result);
5841 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
5842 \int_{x_\text{min}}^{x_\text{max}}
5843 e^{\frac{U_F(x)}{k_B T}}
5844 \text{Integral}_A(x)
5847 <<kramers integ B k functions>>=
5848 double kramers_integ_k_integrandB(double x, void *params)
5850 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5851 double U = (*p->E)(x, p->E_params);
5852 double Fx = p->F * x;
5853 double kbT = k_B * p->env->T;
5854 double intA = kramers_integ_k_integralA(x, params);
5855 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
5856 return exp((U-Fx)/kbT)*intA;
5858 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
5860 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5862 double result, abserr;
5863 size_t neval; /* for qng */
5864 static gsl_integration_workspace *w = NULL;
5865 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5866 f.function = &kramers_integ_k_integrandB;
5870 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5871 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5872 //fprintf(stderr, "integralB = %g\n", result);
5877 With the integrals taken care of, evaluating $k$ is simply
5879 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
5881 <<kramers integ k function declaration>>=
5882 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
5884 <<kramers integ k function>>=
5885 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
5886 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
5887 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
5888 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
5892 <<kramers integ k structure create/destroy function declarations>>=
5893 void *string_create_kramers_integ_param_t(char **param_strings);
5894 void destroy_kramers_integ_param_t(void *p);
5897 <<kramers integ k structure create/destroy functions>>=
5898 kramers_integ_param_t *create_kramers_integ_param_t(double D,
5899 double xmin, double xmax, func_t *E,
5901 destroy_data_func_t *destroy_E_params)
5903 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
5904 assert(ret != NULL);
5909 ret->E_params = E_params;
5910 ret->destroy_E_params = destroy_E_params;
5914 void *string_create_kramers_integ_param_t(char **param_strings)
5916 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
5917 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
5918 void *E_params = string_create_spline_param_t(param_strings+1);
5919 return create_kramers_integ_param_t(
5920 safe_strtod(param_strings[0], __FUNCTION__),
5921 safe_strtod(param_strings[2], __FUNCTION__),
5922 safe_strtod(param_strings[3], __FUNCTION__),
5923 &spline_eval, E_params, destroy_spline_param_t);
5926 void destroy_kramers_integ_param_t(void *params)
5928 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5931 (*p->destroy_E_params)(p->E_params);
5937 Finally we have the GSL spline wrappers:
5938 <<spline functions>>=
5940 <<create/destroy spline>>
5943 <<spline function>>=
5944 double spline_eval(double x, void *spline_params)
5946 spline_param_t *p = (spline_param_t *)(spline_params);
5947 return gsl_spline_eval(p->spline, x, p->acc);
5951 <<create/destroy spline>>=
5952 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
5954 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
5955 assert(ret != NULL);
5956 ret->n_knots = n_knots;
5959 ret->acc = gsl_interp_accel_alloc();
5960 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
5961 gsl_spline_init(ret->spline, x, y, n_knots);
5965 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
5966 void *string_create_spline_param_t(char **param_strings)
5970 double *x=NULL, *y=NULL;
5971 /* break into ordered pairs */
5972 parse_list_string(param_strings[0], SEP, '(', ')',
5973 &num_knots, &knot_coords);
5974 x = (double *)malloc(sizeof(double)*num_knots);
5976 y = (double *)malloc(sizeof(double)*num_knots);
5978 for (i=0; i<num_knots; i++) {
5979 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
5980 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
5983 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
5985 free_parsed_list(num_knots, knot_coords);
5986 return create_spline_param_t(num_knots, x, y);
5989 void destroy_spline_param_t(void *params)
5991 spline_param_t *p = (spline_param_t *)params;
5994 gsl_spline_free(p->spline);
5996 gsl_interp_accel_free(p->acc);
6006 <<kramers integ k global declarations>>=
6007 extern int num_kramers_integ_params;
6008 extern char *kramers_integ_param_descriptions[];
6009 extern char kramers_integ_param_string[];
6012 <<kramers integ k globals>>=
6013 int num_kramers_integ_params = 4;
6014 char *kramers_integ_param_descriptions[] = {
6015 "diffusion D, m^2 / s",
6016 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
6017 "starting integration bound x_min, m",
6018 "ending integration bound x_max, m"};
6019 //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";
6020 //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";
6021 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";
6022 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
6023 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
6024 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
6027 <<kramers integ k model getopt>>=
6028 {"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}
6030 Initialized with Titin I27 parameters\citep{carrion-vazquez99a}.
6031 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
6033 \subsection{Unfolding model implementation}
6037 <<k model includes>>
6038 #include "k_model.h"
6039 <<k model internal definitions>>
6041 <<k model functions>>
6044 <<k model includes>>=
6045 #include <assert.h> /* assert() */
6046 #include <stdlib.h> /* NULL, malloc(), strto*() */
6047 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
6048 #include <stdio.h> /* fprintf(), stdout */
6049 #include <string.h> /* strlen(), strcpy() */
6050 #include <errno.h> /* errno, ERANGE */
6051 #include <gsl/gsl_integration.h>
6052 #include <gsl/gsl_spline.h>
6057 <<k model internal definitions>>=
6058 <<const k structure>>
6059 <<bell k structure>>
6060 <<kbell k structure>>
6061 <<kramers k structure>>
6062 <<kramers integ k structure>>
6063 <<spline param structure>>
6066 <<k model globals>>=
6071 <<kramers k globals>>
6072 <<kramers integ k globals>>
6075 <<k model functions>>=
6077 <<null k functions>>
6078 <<const k functions>>
6079 <<bell k functions>>
6080 <<kbell k functions>>
6081 <<kramers k functions>>
6082 <<kramers integ k functions>>
6085 \subsection{Unfolding model unit tests}
6087 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
6088 <<check-k-model.c>>=
6090 <<k model check includes>>
6091 <<check relative error>>
6093 <<k model test suite>>
6094 <<main check program>>
6097 <<k model check includes>>=
6098 #include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
6099 #include <stdio.h> /* sprintf() */
6100 #include <assert.h> /* assert() */
6101 #include <math.h> /* exp() */
6102 #include <errno.h> /* errno, ERANGE */
6105 #include "k_model.h"
6108 <<k model test suite>>=
6112 <<k model suite definition>>
6115 <<k model suite definition>>=
6116 Suite *test_suite (void)
6118 Suite *s = suite_create ("k model");
6119 <<const k test case defs>>
6120 <<bell k test case defs>>
6122 <<const k test case adds>>
6123 <<bell k test case adds>>
6128 \subsubsection{Constant}
6130 <<const k test case defs>>=
6131 TCase *tc_const_k = tcase_create("Constant k");
6134 <<const k test case adds>>=
6135 tcase_add_test(tc_const_k, test_const_k_create_destroy);
6136 tcase_add_test(tc_const_k, test_const_k_over_range);
6137 suite_add_tcase(s, tc_const_k);
6142 START_TEST(test_const_k_create_destroy)
6144 char *knot[] = {"1","2","3","4","5","6"};
6145 char *params[] = {knot[0]};
6148 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
6149 params[0] = knot[i];
6150 p = string_create_const_k_param_t(params);
6151 destroy_const_k_param_t(p);
6156 START_TEST(test_const_k_over_range)
6158 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
6159 char *knot[] = {"1","2","3","4","5","6"};
6160 char *params[] = {knot[0]};
6163 char param_string[80];
6166 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
6167 params[0] = knot[i];
6168 p = string_create_const_k_param_t(params);
6169 for ( F=Fm; F<FM; F+=dF ) {
6170 fail_unless(const_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
6171 "const_k(%g, %g, &{%s}) = %g != %s",
6172 F, T, knot[i], const_k(F, &env, p), knot[i]);
6174 destroy_const_k_param_t(p);
6180 \subsubsection{Bell}
6182 <<bell k test case defs>>=
6183 TCase *tc_bell = tcase_create("Bell's k");
6186 <<bell k test case adds>>=
6187 tcase_add_test(tc_bell, test_bell_k_create_destroy);
6188 tcase_add_test(tc_bell, test_bell_k_at_zero);
6189 tcase_add_test(tc_bell, test_bell_k_at_one);
6190 suite_add_tcase(s, tc_bell);
6194 START_TEST(test_bell_k_create_destroy)
6196 char *knot[] = {"1","2","3","4","5","6"};
6198 char *params[] = {knot[0], dx};
6201 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
6202 params[0] = knot[i];
6203 p = string_create_bell_param_t(params);
6204 destroy_bell_param_t(p);
6209 START_TEST(test_bell_k_at_zero)
6211 double F=0.0, T=300.0;
6212 char *knot[] = {"1","2","3","4","5","6"};
6214 char *params[] = {knot[0], dx};
6217 char param_string[80];
6220 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
6221 params[0] = knot[i];
6222 p = string_create_bell_param_t(params);
6223 fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
6224 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
6225 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
6226 destroy_bell_param_t(p);
6231 START_TEST(test_bell_k_at_one)
6234 char *knot[] = {"1","2","3","4","5","6"};
6236 char *params[] = {knot[0], dx};
6237 double F= k_B*T/safe_strtod(dx, __FUNCTION__);
6240 char param_string[80];
6243 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
6244 params[0] = knot[i];
6245 p = string_create_bell_param_t(params);
6246 CHECK_ERR(1e-6, safe_strtod(knot[i], __FUNCTION__)*exp(1.0), bell_k(F, &env, p));
6247 //fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__)*exp(1.0),
6248 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
6249 // F, T, knot[i], dx, bell_k(F, &env, p), safe_strtod(knot[i], __FUNCTION__)*exp(1.0));
6250 destroy_bell_param_t(p);
6256 <<kramers k tests>>=
6259 <<kramers k test case defs>>=
6262 <<kramers k test case adds>>=
6265 <<k model function tests>>=
6266 <<check relative error>>
6268 START_TEST(test_something)
6270 double k=5, x=3, last_x=2.0, F;
6271 one_dim_fn_t *handlers[] = {&hooke};
6272 void *data[] = {&k};
6273 double xi[] = {0.0};
6275 int new_active_groups = 1;
6276 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
6277 fail_unless(F = k*x, NULL);
6282 \subsection{Utilities}
6284 The unfolding models can get complicated, and are sometimes hard to explain to others.
6285 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
6286 the overhead of having to understand a full simulation.
6287 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
6288 or special mode, where the operation depends on the specific model selected.
6290 <<k-model-utils.c>>=
6292 <<k model utility includes>>
6293 <<k model utility definitions>>
6294 <<k model utility getopt functions>>
6295 <<k model utility multi model E>>
6296 <<k model utility main>>
6299 <<k model utility main>>=
6300 int main(int argc, char **argv)
6302 <<k model getopt array definition>>
6303 k_model_getopt_t *model = NULL;
6308 int num_param_args; /* for INIT_MODEL() */
6309 char **param_args; /* for INIT_MODEL() */
6311 double special_xmin;
6312 double special_xmax;
6313 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
6314 &Fmax, &special_xmin, &special_xmax, &flags);
6315 if (flags & VFLAG) {
6316 printf("#initializing model %s with parameters %s\n", model->name, model->params);
6318 INIT_MODEL("utility", model, model->params, params);
6321 if (model->k == NULL) {
6322 printf("No k function for model %s\n", model->name);
6326 printf("Fmax = %g <= 0\n", Fmax);
6332 printf("#F (N)\tk (%% pop. per s)\n");
6333 for (i=0; i<=N; i++) {
6334 F = Fmax*i/(double)N;
6335 k = (*model->k)(F, &env, params);
6337 printf("%g\t%g\n", F, k);
6342 if (model->k == NULL || model->k == &bell_k) {
6343 printf("No special mode for model %s\n", model->name);
6346 if (special_xmax <= special_xmin) {
6347 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
6351 double Xrng=(special_xmax-special_xmin), x, E;
6353 printf("#x (m)\tE (J)\n");
6354 for (i=0; i<=N; i++) {
6355 x = special_xmin + Xrng*i/(double)N;
6356 E = multi_model_E(model->k, params, &env, x);
6357 printf("%g\t%g\n", x, E);
6364 if (model->destructor)
6365 (*model->destructor)(params);
6370 <<k model utility multi model E>>=
6371 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
6373 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
6375 if (k == kramers_integ_k)
6376 E = (*p->E)(x, p->E_params);
6378 E = kramers_E(0, env, model_params, x);
6384 <<k model utility includes>>=
6387 #include <string.h> /* strlen() */
6388 #include <assert.h> /* assert() */
6389 #include <errno.h> /* errno, ERANGE */
6392 #include "string_eval.h"
6393 #include "k_model.h"
6396 <<k model utility definitions>>=
6397 <<version definition>>
6398 #define VFLAG 1 /* verbose */
6399 enum mode_t {M_K_OF_F, M_SPECIAL};
6400 <<string comparison definition and globals>>
6401 <<k model getopt definitions>>
6402 <<initialize model definition>>
6403 <<kramers k structure>>
6404 <<kramers integ k structure>>
6408 <<k model utility getopt functions>>=
6411 <<k model utility help>>
6412 <<k model utility get options>>
6415 <<k model utility help>>=
6416 void help(char *prog_name,
6418 int n_k_models, k_model_getopt_t *k_models,
6419 int k_model, double Fmax, double special_xmin, double special_xmax)
6422 printf("usage: %s [options]\n", prog_name);
6423 printf("Version %s\n\n", VERSION);
6424 printf("A utility for understanding the available unfolding models\n\n");
6425 printf("Environmental options:\n");
6426 printf("-T T\tTemperature T (currently %g K)\n", env->T);
6427 printf("-C T\tYou can also set the temperature T in Celsius\n");
6428 printf("Model options:\n");
6429 printf("-k model\tTransition rate model (currently %s)\n",
6430 k_models[k_model].name);
6431 printf("-K args\tTransition rate model argument string (currently %s)\n",
6432 k_models[k_model].params);
6433 printf("There are two output modes. In standard mode, k(F) is printed\n");
6434 printf("For example:\n");
6435 printf(" #Force (N)\tk (%% pop. per s)\n");
6436 printf(" 123.456\t7.89\n");
6438 printf("In special mode, the output depends on the model.\n");
6439 printf("For models defining an energy function E(x), that function is printed\n");
6440 printf(" #Distance (m)\tE (J)\n");
6441 printf(" 0.0012\t0.0034\n");
6443 printf("-m\tChange output to standard mode\n");
6444 printf("-M\tChange output to special mode\n");
6445 printf("-F\tSet the maximum F value for the standard mode k(F) output (currently %g)\n", Fmax);
6446 printf("-x\tSet the minimum x value for the special mode E(x) output (currently %g)\n", special_xmin);
6447 printf("-X\tSet the maximum x value for the special mode E(x) output (currently %g)\n", special_xmax);
6448 printf("-V\tChange output to verbose mode\n");
6449 printf("-h\tPrint this help and exit\n");
6451 printf("Unfolding rate models:\n");
6452 for (i=0; i<n_k_models; i++) {
6453 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
6454 for (j=0; j < k_models[i].num_params ; j++)
6455 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
6456 printf("\t\tdefaults: %s\n", k_models[i].params);
6461 <<k model utility get options>>=
6462 void get_options(int argc, char **argv, environment_t *env,
6463 int n_k_models, k_model_getopt_t *k_models,
6464 enum mode_t *mode, k_model_getopt_t **model,
6465 double *Fmax, double *special_xmin, double *special_xmax,
6466 unsigned int *flags)
6468 char *prog_name = NULL;
6469 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
6471 extern char *optarg;
6472 extern int optind, optopt, opterr;
6476 /* setup defaults */
6478 prog_name = argv[0];
6479 env->T = 300.0; /* K */
6483 *Fmax = 1e-9; /* N */
6484 *special_xmax = 1e-8;
6485 *special_xmin = 0.0;
6487 while ((c=getopt(argc, argv, options)) != -1) {
6489 case 'T': env->T = safe_strtod(optarg, "-T"); break;
6490 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
6492 k_model = index_k_model(n_k_models, k_models, optarg);
6493 *model = k_models+k_model;
6496 k_models[k_model].params = optarg;
6498 case 'm': *mode = M_K_OF_F; break;
6499 case 'M': *mode = M_SPECIAL; break;
6500 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
6501 case 'x': *special_xmin = safe_strtod(optarg, "-x"); break;
6502 case 'X': *special_xmax = safe_strtod(optarg, "-X"); break;
6503 case 'V': *flags |= VFLAG; break;
6505 fprintf(stderr, "unrecognized option '%c'\n", optopt);
6506 /* fall through to default case */
6508 help(prog_name, env, n_k_models, k_models, k_model, *Fmax, *special_xmin, *special_xmax);
6517 \section{\LaTeX\ documentation}
6519 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
6520 The comment blocks are extracted (with nicely formatted code blocks), using
6521 <<latex makefile lines>>=
6522 $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6523 noweave -latex -index -delay $< > $@
6524 $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR)
6528 <<latex makefile lines>>=
6529 $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
6531 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
6532 $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
6533 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
6534 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
6535 mv $(BUILD_DIR)/sawsim.pdf $@
6537 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
6538 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
6539 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
6541 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
6542 <<latex makefile lines>>=
6544 rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
6545 $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
6546 $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.toc \
6547 $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib
6548 clean_latex : semi-clean_latex
6549 rm -f $(DOC_DIR)/sawsim.pdf
6554 [[make]] is a common utility on *nix systems for managing dependencies while building software.
6555 In this case, we have several \emph{targets} that we'd like to build:
6556 the main simulation program \prog;
6557 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
6558 the documentation [[sawsim.pdf]];
6559 and the [[Makefile]] itself.
6560 Besides the generated files, there is the utility target
6561 [[clean]] for removing all generated files except [[Makefile]].
6563 % [[dist]] for generating a distributable tar file.
6565 Extract the makefile with
6566 `[[notangle -Rmakefile src/sawsim.nw | sed 's/ /\t/' > Makefile]]'.
6567 The sed is needed because notangle replaces tabs with eight spaces,
6568 but [[make]] needs tabs.
6570 # Customize compilation
6575 # Decide what directories to put things in
6580 # Note: Cannot use, for example, `./build', because make eats the `./'
6581 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
6583 # Modules (X.c, X.h) defined in the noweb file
6585 # Check binaries (check_*) defined in the noweb file
6587 # Modules defined outside the noweb file
6588 FREE_SAWSIM_MODS = interp tavl
6590 <<list module makefile lines>>
6591 <<tension balance module makefile lines>>
6592 <<tension model module makefile lines>>
6593 <<k model module makefile lines>>
6594 <<parse module makefile lines>>
6595 <<string eval module makefile lines>>
6596 <<accel k module makefile lines>>
6598 SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS)
6600 # Everything needed for sawsim under one roof. sawsim.c must come first
6601 SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \
6602 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h)
6603 # Libraries needed to compile sawsim
6604 LIBS = gsl gslcblas m
6605 CHECK_LIBS = gsl gslcblas m check
6606 # The non-check binaries generated
6607 BINS = sawsim tension_model_utils k_model_utils sawsim_profile
6610 # Define the major targets
6611 all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
6613 view : $(DOC_DIR)/sawsim.pdf
6615 profile : $(BIN_DIR)/sawsim_profile | $(BIN_DIR)
6616 $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \
6617 -mnull -Mwlc -A0.39e-9,28e-9 -F8
6618 gprof $(BIN_DIR)/sawsim_profile gmon.out > $@
6619 check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim
6620 $(SHELL) -e -c 'for B in $^; do ./$$B; done'
6621 clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \
6622 clean_tension_model_utils \
6623 clean_k_model_utils clean_latex clean_check_sawsim
6624 rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \
6625 $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \
6626 $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
6627 $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
6628 $(BUILD_DIR)/global.h ./gmon.out
6629 $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
6631 # Various builds of sawsim
6632 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
6633 $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6634 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
6635 $(CC) $(CFLAGS) $(LDFLAGS) -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6636 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
6637 $(CC) $(CFLAGS) $(LDFLAGS) -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6639 # Create the directories
6640 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
6643 # Copy over the source external to sawsim
6644 # Note: Cannot use, for example, `./build', because make eats the `./'
6645 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
6647 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6651 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6655 ## Basic source generated with noweb
6656 # The central files sawsim.c and global.h...
6657 $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6659 $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6660 notangle -Rglobal.h $< > $@
6661 # ... and the modules
6662 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6663 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6664 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6665 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6666 $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6667 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6668 # Note: I use `_' as a space in my file names, but noweb doesn't like
6669 # that so we use `-' to name the noweb roots and substitute here.
6671 ## Building the unit-test programs
6673 $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
6674 notangle -Rchecks $< > $@
6675 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
6676 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
6677 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR)
6678 $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
6679 clean_check_sawsim :
6680 rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
6681 # ... and the modules
6683 $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
6684 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \
6685 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \
6686 $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6687 $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6688 $$(BUILD_DIR)/global.h | $(BIN_DIR)
6689 $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
6690 $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
6692 # TODO: clean up the required modules too
6693 $(CHECK_BINS:%=clean_%) :
6694 rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
6696 # Cleaning up the modules
6698 $(SAWSIM_MODS:%=clean_%) :
6699 rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h)
6701 <<tension model utils makefile lines>>
6702 <<k model utils makefile lines>>
6703 <<latex makefile lines>>
6705 Makefile : $(SRC_DIR)/sawsim.nw
6706 notangle -Rmakefile $< | sed 's/ /\t/' > $@
6708 This is fairly self-explanatory. We compile a dynamically linked
6709 version ([[sawsim]]) and a statically linked version
6710 ([[sawsim_static]]). The static version is about 9 times as big, but
6711 it runs without needing dynamic GSL libraries to link to, so it's a
6712 better format for distributing to a cluster for parallel evaluation.
6716 \subsection{Hookean springs in series}
6717 \label{sec.math_hooke}
6720 The effective spring constant for $n$ springs $k_i$ in series is given by
6722 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6728 F &= k_i x_i &&\forall i \le n \\
6729 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
6730 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
6731 F &= k_1 x_1 = k_\text{eff} x \\
6732 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
6733 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6738 \addcontentsline{toc}{section}{References}
6739 \bibliography{sawsim}