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}
130 Sawim has been published\citep{king10}! It is, as far as I know, the
131 only open-source Monte Carlo force spectroscopy simulator, but similar
132 closed source simulators have been around since the seminal paper by
133 \citet{rief97a}. In chronological order, major contributions have
137 \item \citet{rief97a} \citeyear{rief97a}:
138 \item \citet{rief97b} \citeyear{rief97b}:
139 \item \citet{kellermayer97} \citeyear{kellermayer97}:
140 \item \citet{rief98} \citeyear{rief98}:
141 first paper to focus mainly on the simulation
142 \item \citet{oberhauser98} \citeyear{oberhauser98}:
143 \item \citet{carrion-vazquez99a} \citeyear{carrion-vazquez99a}:
144 \item \citet{best02} \citeyear{best02}:
145 pointed out large fit valley
146 \item \citet{zinober02} \citeyear{zinober02}:
147 introduced the scaffold effect
148 \item \citet{jollymore09} \citeyear{jollymore09}:
149 \item \citet{king10} \citeyear{king10}:
154 \subsection{About this document}
156 This document was written in \citetalias{sw:noweb} with code blocks in
157 \lang\ and comment blocks in \LaTeX.
159 \subsection{System Requirements and Dependencies}
162 \subsection{Change Log}
164 Version 0.0 used only the unfolded domain WLC to determine the tension
165 and had a fixed timestep $dt$.
167 Version 0.1 added adaptive timesteps, adjusting so that the unfolding
168 probability for a given domain was constant.
170 Version 0.2 added Kramers' model unfolding rate calculations, and a
171 switch to select Bell's or Kramers' model. This induced a major
172 rewrite, introducing the tension group abstraction, and a split to
173 multiple \lang\ source files. Also added a unit-test suites for
174 sanity checking, and switched to SI units throughout.
176 Version 0.3 added integrated Kramers' model unfolding rate
177 calculations. Also replaced some of my hand-coded routines with GNU
178 Scientific Library \citepalias{galassi05} calls.
180 Version 0.4 added Python evaluation for the saddle-point Kramers'
181 model unfolding rate calculations, so that the model functions could
182 be controlled from the command line. Also added interpolating lookup
183 tables to accelerate slow unfolding rate model evaluation, and command
184 line control of tension groups.
186 Version 0.5 added the stiffness environmental parameter and it's
187 associated unfolding models.
189 Version 0.6 generalized from two state unfolding to arbitrary
190 $N$-state rate reactions. This allows simulations of
191 unfolding-refolding, metastable intermediates, etc., but required
192 reorganizing the command line interface.
194 Version 0.7 added tension model inverses, which seems to reduce
195 computation time by about a factor of 2. I was expecting more, but
196 I'll take what I can get.
198 Version 0.8 fixed a minor bug in unfolding probability for
199 multi-domain groups. The probability of at least one domain unfolding
200 had been calculated as $P_N=NP_1$, rather than $P_N=1-(1-P_1)^N$.
201 However, for small $P$ the two are equivalent.
203 Version 0.9 added the [[-P]] option to sawsim, setting the target
204 probability for the per-state transition $P_N$, not the per-domain
207 Version 0.10 added the [[-F]] option to sawsim, setting a limit on the
208 force change $dF$ in a single timestep for continuous pulling
209 simulations. It also added the piston tension model (Section
210 \ref{sec.piston_tension}), and adjusted the stiffness calculations to
211 handle domains with undefined stiffness.
213 <<version definition>>=
214 #define VERSION "0.10"
220 sawsim - simulating single molecule mechanical unfolding.
221 Copyright (C) 2008-2010, William Trevor King
223 This program is free software; you can redistribute it and/or
224 modify it under the terms of the GNU General Public License as
225 published by the Free Software Foundation; either version 3 of the
226 License, or (at your option) any later version.
228 This program is distributed in the hope that it will be useful, but
229 WITHOUT ANY WARRANTY; without even the implied warranty of
230 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
231 See the GNU General Public License for more details.
233 You should have received a copy of the GNU General Public License
234 along with this program. If not, see <http://www.gnu.org/licenses/>.
236 The author may be contacted at <wking@drexel.edu> on the Internet, or
237 write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
238 Philadelphia PA 19104, USA.
251 The unfolding system is stored as a chain of \emph{domains} (Figure
252 \ref{fig.domain_chain}). Domains can be folded, globular protein
253 domains, unfolded protein linkers, AFM cantilevers, or other
254 stretchable system components. The system is modeled as graph of
255 possible domain states with transitions between them (Figure
256 \ref{fig.domain_states}).
260 \subfloat[][]{\label{fig.domain_chain}
261 \begin{tikzpicture}[->,node distance=2cm,font=\footnotesize]
262 \tikzstyle{every state}=[fill,draw=red!50,very thick,fill=red!20]
263 \node[state] (A) {domain 1};
264 \node[state] (B) [below of=A] {domain 2};
265 \node[state] (C) [below of=B] {.~.~.};
266 \node[state] (D) [below of=C] {domain $N$};
267 \node (S) [below of=D] {Surface};
268 \node (E) [above of=A] {};
270 \path[-] (A) edge (B)
271 (B) edge node [right] {Tension} (C)
274 \path[->,green] (A) edge node [right,black] {Extension} (E);
278 \subfloat[][]{\label{fig.domain_states}
279 \begin{tikzpicture}[->,node distance=2.5cm,shorten <=1pt,shorten >=1pt,font=\footnotesize]
280 \tikzstyle{every state}=[fill,draw=blue!50,very thick,fill=blue!20]
281 \node[state] (A) {cantilever};
282 \node[state] (C) [below of=A] {transition};
283 \node[state] (B) [left of=C] {folded};
284 \node[state] (D) [right of=C] {unfolded};
286 \path (B) edge [bend left] node [above] {$k_1$} (C)
287 (C) edge [bend left] node [below] {$k_1'$} (B)
288 edge [bend left] node [above] {$k_2$} (D)
289 (D) edge [bend left] node [below] {$k_2'$} (C);
292 \caption{\subref{fig.domain_chain} Extending a chain of domains. One end
293 of the chain is fixed, while the other is extended at a constant
294 speed. The domains are coupled with rigid linkers, so the domains
295 themselves must stretch to accomodate the extension.
296 \subref{fig.domain_states} Each domain exists in a discrete state. At
297 each timestep, it may transition into another state following a
298 user-defined state matrix such as this one, showing a metastable
299 transition state and an explicit ``cantilever'' domain.}
303 Each domain contributes to the total tension, which depends on the
304 chain extension. The particular model for the tension as a function
305 of extension is handled generally with function pointers. So far the
306 following models have been implemented:
308 \item Null (Appendix \ref{sec.null_tension}),
309 \item Constant (Appendix \ref{sec.const_tension}),
310 \item Hookean spring (Appendix \ref{sec.hooke}),
311 \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
312 \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
313 \item Piston (Appendix \ref{sec.piston_tension}),
316 The tension model and parameters used for a particular domain depend
317 on the domain's current state. The transition rates between states
318 are also handled generally with function pointers, with
321 \item Null (Appendix \ref{sec.null_k}),
322 \item Bell's model (Appendix \ref{sec.bell}),
323 \item Spring constant corrected Bell model (Appendix \ref{sec.kbell}),
324 \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
325 \item Kramers' saddle point model (Appendix \ref{sec.kramers}).
328 The unfolding of the chain domains is modeled in two stages. First
329 the tension of the chain is calculated. Then the tension (assumed to
330 be constant over the length of the chain, see Section
331 \ref{sec.timescales}), is applied to each domain, and used to
332 calculate the probability of that domain changing state in the next
333 timestep $dt$. Then the time is stepped forward, and the process is
334 repeated. The simulation is complete when a pre-set time limit
335 [[t_max]] is reached or a pre-selected state [[stop_state]] is empty.
337 int main(int argc, char **argv)
339 <<tension model getopt array definition>>
340 <<k model getopt array definition>>
342 /* variables initialized in get_options() */
343 list_t *state_list=NULL, *transition_list=NULL;
344 environment_t env={0};
345 double P, dF, t_max, dt_max, v, x_step;
346 state_t *stop_state=NULL;
348 /* variables used in the simulation loop */
349 double t=0, x=0, dt, F; /* time, distance, timestep, force */
350 int transition=1; /* boolean marking a transition for tension and timestep optimization */
352 get_options(argc, argv, &P, &dF, &t_max, &dt_max, &v, &x_step,
353 &stop_state, &env, NUM_TENSION_MODELS, tension_models,
354 NUM_K_MODELS, k_models, &state_list, &transition_list);
357 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Distance (m)\tForce (N)\tStiffness (N/m)\n");
358 if (flags & FLAG_OUTPUT_TRANSITION_FORCES) printf("#Force (N)\tInitial state\tFinal state\n");
359 while ((t_max > 0 && t <= t_max) || (stop_state != NULL && stop_state->num_domains > 0)) {
360 F = find_tension(state_list, &env, x, transition, 0);
362 dt = determine_dt(state_list, transition_list, &env, F, x, dF, P, dt_max, v, transition);
364 dt = determine_dt(state_list, transition_list, &env, F, x, dF, P, dt_max, 0, transition);
365 transition=random_transitions(transition_list, F, dt, &env);
366 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness);
371 if (dt == dt_max) { /* step completed */
374 } else { /* still working on this step */
379 destroy_state_list(state_list);
380 destroy_transition_list(transition_list);
384 @ The meat of the simulation is bundled into the three functions
385 [[find_tension]], [[determine_dt]], and [[random_transitions]].
386 [[find_tension]] is discussed in Section \ref{sec.find_tension},
387 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
388 [[random_transitions]] in Section \ref{sec.transition_rate}.
390 The stretched distance is extended in one of two modes depending on
391 [[xstep]]. If [[xstep == 0]], then the pulling is continuous, at
392 least within the limits of the inherent discretization of a time
393 stepped approach. At any rate, the timestep is limited by the maximum
394 allowed timestep [[dt_max]] and the maximum allowed unfolding
395 probability in a given step [[P]]. In the second mode [[xstep > 0]],
396 and the end to end distance increases in discrete steps of that
397 length. The time between steps is chosen to maintain an average
398 unfolding velocity [[v]]. This approach is less work to simulate than
399 smooth pulling and also closer to the reality of many experiments, but
400 it is practically impossible to treat analytically. With both pulling
401 styles implemented, the effects of the discretization can be easily
402 calculated, bridging the gap between theory and experiment.
404 Environmental parameters important in determining reaction rates and
405 tensions (e.g.\ temperature, pH) are stored in a single structure to
406 facilitate extension to more complicated models in the future.
407 <<environment definition>>=
408 typedef struct environment_struct {
418 #define DOUBLE_PRECISION 1e-12
420 <<environment definition>>
421 <<one dimensional function definition>>
422 <<create/destroy definitions>>
424 #endif /* GLOBAL_H */
426 Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
427 included multiple times.
429 \section{Simulation functions}
431 <<simulation functions>>=
435 <<does domain transition>>
436 <<randomly transition domains>>
440 \label{sec.find_tension}
442 Because the stretched system may be made up of several parts (folded
443 domains, unfolded domains, spring-like cantilever, \ldots), we will
444 assign the domains to states with tension models and groups. The
445 models are used to determine the tension of all the domain in a given
446 state following some model (Hookean spring, WLC, \ldots). The domains
447 are assumed to be commutative, so ordering is ignored. The
448 interactions between the groups are assumed to be linear, but the
449 interactions between domains of the same group need not be. Each
450 model has a tension handler function, which gives the tension $F$
451 needed to stretch a given group of domains a distance $x$.
453 Groups represent collections of states which the model should treat as
454 a single entity. To understand the purpose of groups, consider a
455 system with two unfolded domain states (e.g.\ for two different
456 proteins) that were both modeled as WLCs. If both unfolded states had
457 the same persistence length, it would be wise to assign them to the
458 same group. This would reduce both the computational cost of
459 balancing the tension and the error associated with the inter-group
460 interaction linearity. Note that group numbers only matter
461 \emph{within} model classes, since grouping domains with seperate
462 models doesn't make sense. Within designated groups, the tension
463 parameters for each domain are still checked for compatibility, so if
464 you accidentally assign incompatible domains to the same group you
465 should get an appropriate error message.
467 <<find tension definitions>>=
468 #define NUM_TENSION_MODELS 6
472 <<tension handler array global declaration>>=
473 extern one_dim_fn_t *tension_handlers[];
476 <<tension handler array global>>=
477 one_dim_fn_t *tension_handlers[] = {
479 &const_tension_handler,
487 <<tension model global declarations>>=
488 <<tension handler array global declaration>>
491 <<tension handler types>>=
492 typedef struct tension_handler_data_struct {
493 list_t *group_tension_params; /* one item for each domain in the group */
496 } tension_handler_data_t;
500 After sorting the chain into separate groups $G_i$, with tension
501 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
502 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
503 \forall i,j$. Note that there may be several states within each
504 group. [[find_tension]] is basically a wrapper to feed proper input
505 into the [[tension_balance]] function. [[transition]] should be set
506 to 0 if there were no transitions since the previous call to
507 [[find_tension]] to support some optimizations that assume a small
508 increase in tension between steps (only true if no transition has
509 occured). While we're messing with the tension handlers and if
510 [[const_env==0]], we also grab a measure of the current linker
511 stiffness for the environmental variable, which is needed by some of
512 the unfolding rate models (e.g.\ the linker-stiffness corrected Bell
513 model, Appendix \ref{sec.kbell}).
515 double find_tension(list_t *state_list, environment_t *env, double x,
516 int transition, int const_env)
517 { // TODO: !cont_env -> step_call, only alter env or statics if step_call==1
518 static int num_active_groups=0;
519 static one_dim_fn_t **per_group_handlers = NULL;
520 static one_dim_fn_t **per_group_inverse_handlers = NULL;
521 static void **per_group_data = NULL; /* array of pointers to tension_handler_data_t */
522 static double last_x;
523 tension_handler_data_t *tdata;
524 double F, *pStiffness;
528 fprintf(stderr, "%s:\tfinding tension at distance %g%s\n", __FUNCTION__, x, transition ? " after a transition" : "");
531 if (transition != 0 || num_active_groups == 0) { /* setup statics */
532 /* get new statics, freeing the old ones if they aren't NULL */
533 get_active_groups(state_list, &num_active_groups, &per_group_handlers, &per_group_inverse_handlers, &per_group_data);
535 /* fill in the group handlers and params */
537 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]);
539 for (i=0; i<num_active_groups; i++) {
540 tdata = (tension_handler_data_t *) per_group_data[i];
542 fprintf(stderr, "%s:\tactive group %d of %d (data %p, param list %p)", __FUNCTION__, i, num_active_groups, tdata, tdata->group_tension_params);
543 list_print(stderr, tdata->group_tension_params, " parameters");
546 /* tdata->persist continues unchanged */
549 } /* else continue with the current statics */
554 pStiffness = &env->stiffness;
556 F = tension_balance(num_active_groups, per_group_handlers, per_group_inverse_handlers, (void **)per_group_data, last_x, x, pStiffness);
562 @ For the moment, we will restrict our group boundaries to a single
563 dimension, so $\sum_i x_i = x$, our total extension (it is this
564 restriction that causes the groups to interact linearly). We'll also
565 restrict our extensions to all be positive. With these restrictions,
566 the problem of balancing the tensions reduces to solving a system of
567 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
568 the number of active groups. In general this can be fairly
569 complicated, but without much loss of practicality we can restrict
570 ourselves to strictly monotonically increasing, non-negative tension
571 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
572 simpler. For example, it guarantees the existence of a unique, real
573 solution for finite forces. See Appendix \ref{sec.tension_balance}
574 for the tension balancing implementation.
576 Each group must define a parameter structure if the tension model is
578 a function to create the parameter structure,
579 a function to destroy the parameter structure, and
580 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as its data argument.
581 The tension-specific parameter structures are contained in the domain
582 group field. For optimization, a group may also define an inverse
583 tension function as an optimization. See the Section
584 \ref{sec.model_selection} for a discussion on the command line
585 framework and inverse function discussion. See the worm-like chain
586 implementation for an example (Section \ref{sec.wlc}).
588 The major limitation of our approach is that all of our tension models
589 are Markovian (memory-less), which is not really a problem for the
590 chains (see \citet{evans99} for WLC thermalization rate calculations).
592 \subsection{Transition rate}
593 \label{sec.transition_rate}
595 Each state transition is modeled as unimolecular, first order reaction
597 \text{State 1} \xrightarrow{k} \text{State 2}
599 With the rate constant $k$ defined by
601 \frac{d[\text{State 2}]}{dt} = -\frac{d[\text{State 1}]}{dt} = k [\text{State 1}].
603 So the probability of a given domain transitioning in a short time
609 For a group of $N$ identical domains, and therefore identical
610 unfolding probabilities $P_1$, the probability of $n$ domains
611 unfolding in a given timestep is given by the binomial distribution
613 P(n) = \frac{N!}{n!(N-n)!}P_1^n(1-P_1)^{(N-n)}.
615 The probability that \emph{none} of these domains unfold is then
619 so the probability that \emph{at least one} domain unfolds is
623 Note that for small $P$,
625 p(n>0) = 1-(1-NP_1+\mathcal{O}(P_1^2)) = NP_1 - \mathcal{O}(P^2)
628 This conversion from $P_1$ to $P_N$ is handled by the [[pN]] macro
630 /* find multi-domain transition rate from the single domain rate */
631 #define PN(P,N) (1.0-pow(1.0-(P), (N)))
635 We can also define a macro to find the approprate $dt$ to achieve a
636 target $P_N$ with a given $k$ and $N$.
638 P_N &= 1-(1-P_1)^N \\
639 1-P_1 &= (1-P_N)^{1/N} \\
640 P_1 &= 1-(1-P_N)^{1/N}
643 #define P1(PN,N) (1.0-pow(1.0-(PN), 1.0/(N)))
646 We take some time to discuss the meaning of $p(n>1)$
647 (i.e. multi-unfolding timesteps) in Section \ref{sec.timescales}.
649 <<does domain transition>>=
650 int domain_transitions(double F, double dt, environment_t *env, transition_t *transition)
651 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to transition structure */
653 k = accel_k(transition->k, F, env, transition->k_params);
654 //(*transition->k)(F, env, domain->k_params);
655 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
656 return happens(PN(k*dt, N_INIT(transition)));
658 @ [[happens]] is a random decision making function defined in Appendix \ref{sec.utils}.
660 <<randomly transition domains>>=
661 int random_transitions(list_t *transition_list, double tension,
662 double dt, environment_t *env)
663 { /* returns 1 if there was a transition and 0 otherwise */
664 while (transition_list != NULL) {
665 if (T(transition_list)->initial_state->num_domains > 0
666 && domain_transitions(tension, dt, env, T(transition_list))) {
667 if (flags & FLAG_OUTPUT_TRANSITION_FORCES)
668 fprintf(stdout, "%g\t%s\t%s\n", tension, T(transition_list)->initial_state->name, T(transition_list)->final_state->name);
669 T(transition_list)->initial_state->num_domains--;
670 T(transition_list)->final_state->num_domains++;
671 return 1; /* our one domain has transitioned, stop looking for others */
673 transition_list = transition_list->next;
679 \subsection{Timescales}
680 \label{sec.timescales}
682 The simulation assumes that chain equilibration is instantaneous
683 relative to the simulation timestep $dt$, so that tension is uniform
684 along the chain. The quality of this assumption depends on your
685 particular chain. For example, a damped spring thermalizes on a
686 timescale of order $\tau = 1/\gamma$ where the damping ratio $\gamma
687 \equiv \omega_0(-\zeta \pm \sqrt{\zeta^2-1}$, the natural angular
688 frequency $omega_0 \equiv \sqrt{k/m}$, $\zeta \equiv b/2\sqrt{km}$,
689 and $b$ sets the drag $F_b = -bv$. For our cantilevers $\tau$ is
690 on the order of milliseconds, which is longer than a timestep.
691 However, the approximation is still reasonable, because there is
692 rarely unfolding during the cantilever return. You could, of course,
693 take cantilever, WLC, etc.\ response times into effect, but that
694 would significantly complicate a simulation that seems to work
695 well enough without bothering :p. For WLC and FJC relaxation times,
696 see the Appendix of \citet{evans99} and \citet{kroy07}.
698 Besides assuming our timestep is much greater than equilibration
699 times, we also force it to be short enough so that only one domain may
700 unfold in a given timestep. For an unfolding timescale of $dt_u =
701 1/k$, we adapt our timesteps to keep $dt \ll dt_u$, so the probability
702 of two domains unfolding in the same timestep is small (see
703 [[determine_dt]] in Section \ref{sec.adaptive_dt} and [[P]] in
704 [[main()]] in Section \ref{sec.main} set by [[-P]] in
705 [[get_options()]] in Section \ref{sec.get_options}). This approach
706 breaks down as the adaptive timestep scheme approaches the transition
707 timestep $dt_t$, but $dt_t \sim 1\U{$\mu$s}$ for Ig27-like proteins
708 \citep{klimov00}, so this shouldn't be a problem. To reassure
709 yourself, you can ask the simulator to print the smallest timestep
710 that was used with TODO.
712 Even if the likelyhood of two domains unfolding in the same timestep
713 is small, if you run long enough simulations it will eventually occur.
714 In this case, we assume that $dt_t \ll dt$, so even if two domains
715 would unfold if we stuck to our unfolding rate $k$ for an entire
716 timestep $dt$, in ``reality'' only one of those domains unfolds.
717 Because we do not know when in the timestep the unfolding took place,
718 we move on to the next timestep without checking for further
719 unfoldings. This ``unchecked timestep portion'' should not contribute
720 significant errors to the calculation, because if $dt$ was small
721 enough that unfolding was unlikely at the pre-unfolding force, it
722 should be even more unlikely at the post-unfolding force, especially
723 over only a fraction of the timestep. The only dangerous cases for
724 the current implementation are those where unfolding intermediates are
725 much less stable than their precursor states, in which case an
726 unfolding event during the remaining timestep may actually be likely.
727 A simple workaround for such cases is to pick the value for [[P]] to
728 be small enough that the $dt \ll dt_u$ assumption survives even under
729 a transition populating the unstable intermediate.
731 \subsection{Adaptive timesteps}
732 \label{sec.adaptive_dt}
734 We'd like to pick $dt$ so the probability of changing state
735 (e.g.\ unfolding) in the next timestep is small. If we don't adapt our
736 timestep, we also risk breaking our approximation that $F(x' \in [x,
737 x+v \cdot dt]) = F(x)$ or that only one domain changes state in a
738 given timestep. The simulation would have been accurate for
739 sufficiently small fixed timesteps, but adaptive timesteps will allow
740 us to move through low tension regions in fewer steps, leading to a
741 more efficient simulation.
743 Consider the two state folded $\rightarrow$ unfolded transition.
744 Because $F(x)$ is monotonically increasing between unfoldings,
745 excessively large timesteps will lead to erroneously small unfolding
746 rates (an thus, higher average unfolding force).
748 The actual adaptive timestep implementation is not particularly
749 interesting, since we are only required to reduce $dt$ to somewhere
750 below a set threshold, so I've removed it to Appendix
751 \ref{sec.adaptive_dt_implementation}.
757 The models provide the physics of the simulation, but the simulation
758 allows interchangeable models, and we are currently focusing on the
759 simulation itself, so we remove the actual model implementation to the
762 The tension models are defined in Appendix \ref{sec.tension_models},
763 and unfolding models are defined in Appendix \ref{sec.k_models}.
766 #define k_B 1.3806503e-23 /* J/K */
771 \section{Command line interface}
773 <<get option functions>>=
775 <<index tension model>>
780 \subsection{Model selection}
781 \label{sec.model_selection}
783 The main difficulty with the command line interface in \prog\ is
784 developing an intuitive method for accessing the possibly large number
785 of available models. We'll treat this generally by defining an array
786 of available models, containing their important parameters.
788 <<get options definitions>>=
789 <<model getopt definitions>>
792 <<create/destroy definitions>>=
793 typedef void *create_data_func_t(char **param_strings);
794 typedef void destroy_data_func_t(void *param_struct);
797 <<model getopt definitions>>=
798 <<tension model getopt definitions>>
799 <<k model getopt definitions>>
802 In order to choose models, we need a method of assembling a reaction
803 graph such as that in Figure \ref{fig.domain_states} and population
804 the graph with a set of domains. First we declare the domain states
807 [[-s <state-name>,<model>[,<group>] -S <model-params> [-N <initial-population]]] \\
811 [[-s unfolded,wlc,1 -S 1e-8,4e-10 -N 5]]
813 This creates the state named [[unfolded]], using the WLC model and one
814 as a group number (defaults to 0, see Section \ref{sec.find_tension}).
815 The tension parameters are then set to [[1e-8,4e-10]], which in the
816 case of the WLC are the contour and persistence lengths respectively.
817 Finally we populate the state by adding five domains to the state.
818 The name of the state is importent for identifying states later.
820 Once the domains have been created and populated, you can added
821 transition rates following
823 [[-k <initial-state>,<final-state>,<model> -K <model-params>]] \\
827 [[-k folded,unfolded,bell -K 3.3e-4,0.25e-9]]
829 This creates a reaction going from the [[folded]] state to the
830 [[unfolded]] state, with the rate constant given by the Bell model
831 (Appendix \ref{sec.bell}). The unfolding parameters are then set to
832 [[3.3e-4,0.25e-9]], which in the case of the Bell model are the
833 unforced rate and transition state distance respectively.
835 \subsubsection{Tension}
837 To access the various tension models from the command line, we define
838 a structure that contains all the neccessary information about the
840 <<tension model getopt definitions>>=
841 typedef struct tension_model_getopt_struct {
842 one_dim_fn_t *handler; /* fn implementing the model on a group */
843 one_dim_fn_t *inverse_handler; /* fn implementing inverse on group */
844 char *name; /* name string identifying the model */
845 char *description; /* a brief description of the model */
846 int num_params; /* number of extra params for the model */
847 char **param_descriptions; /* descriptions of extra params */
848 char *params; /* default values for the extra params */
849 create_data_func_t *creator; /* param-string -> model-data-struct fn */
850 destroy_data_func_t *destructor; /* fn to free the model data structure */
851 } tension_model_getopt_t;
852 @ Set [[inverse_handler]] to [[NULL]] if you wish to invert using
853 bisection. Obviously, this will be slower than computing an
854 analytical inverse and more prone to rounding errors, but it may be
855 easier if a closed-form inverse does not exist.
857 <<tension model getopt array definition>>=
858 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
859 <<null tension model getopt>>,
860 <<constant tension model getopt>>,
861 <<hooke tension model getopt>>,
862 <<worm-like chain tension model getopt>>,
863 <<freely-jointed chain tension model getopt>>,
864 <<piston tension model getopt>>
868 \subsubsection{Unfolding rate}
870 <<k model getopt definitions>>=
871 #define NUM_K_MODELS 6
873 typedef struct k_model_getopt_struct {
878 char **param_descriptions;
880 create_data_func_t *creator;
881 destroy_data_func_t *destructor;
885 <<k model getopt array definition>>=
886 k_model_getopt_t k_models[NUM_K_MODELS] = {
887 <<null k model getopt>>,
888 <<const k model getopt>>,
889 <<bell k model getopt>>,
890 <<kbell k model getopt>>,
891 <<kramers k model getopt>>,
892 <<kramers integ k model getopt>>
899 void help(char *prog_name, double P, double dF,
900 double t_max, double dt_max, double v, double x_step,
901 state_t *stop_state, environment_t *env,
902 int n_tension_models, tension_model_getopt_t *tension_models,
903 int n_k_models, k_model_getopt_t *k_models,
904 int k_model, list_t *state_list)
909 if (state_list != NULL)
910 state = S(tail(state_list));
912 printf("usage: %s [options]\n", prog_name);
913 printf("Version %s\n\n", VERSION);
914 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
916 printf("Simulation options:\n");
918 printf("-q state\tStop simulation when this state empties (- for no limit, currently %s)\n", stop_state ? stop_state->name : "-" );
919 printf("-t tmax\tMaximum simulated time (-1 for no limit, currently %g s)\n", t_max);
920 printf("-d dt\tMaximum allowed timestep dt (currently %g s)\n", dt_max);
921 printf("-P P\tMaximum transition probability for dt (currently %g)\n", P);
922 printf("-F dF\tMaximum change in force over dt. Only relevant if dx > 0. (currently %g N)\n", dF);
923 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
924 printf("-x dx\tPulling step size (currently %g m)\n", x_step);
925 printf("\tWhen dx = 0, the pulling is continuous.\n");
926 printf("\tWhen dx > 0, the pulling occurs in discrete steps.\n");
928 printf("Environmental options:\n");
929 printf("-T T\tTemperature T (currently %g K)\n", env->T);
930 printf("-C T\tYou can also set the temperature T in Celsius\n");
932 printf("State creation options:\n");
933 printf("-s name,model[[,group],params]\tCreate a new state.\n");
934 printf("Once you've set up the state, you may populate it with domains\n");
935 printf("-N n\tAdd n domains in the current state (%s)\n", state ? state->name : "-");
937 printf("Transition creation options:\n");
938 printf("-k initial_state,final_state,model[,params]\tCreate a new transition\n");
940 printf("To double check your reaction graph, you can produce graphviz dot output\n");
941 printf("describing the current states and transitions.\n");
942 printf("-D\tPrint dot description of states and transitions and exit.\n");
944 printf("Simulation output mode options:\n");
945 printf("There are two output modes. In standard mode, only the transition\n");
946 printf("events are printed. For example:\n");
947 printf(" #Force (N)\tInitial state\tFinal state\n");
948 printf(" 123.456e-12\tfolded\tunfolded\n");
950 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
951 printf(" #Distance (m)\tForce (N)\tStiffness (N/m)\n");
952 printf(" 0.001\t0.0023\t0.002\n");
954 printf("-V\tChange output to verbose mode.\n");
955 printf("-h\tPrint this help and exit.\n");
958 printf("Tension models:\n");
959 for (i=0; i<n_tension_models; i++) {
960 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
961 for (j=0; j < tension_models[i].num_params ; j++)
962 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
963 printf("\t\tdefaults: %s\n", tension_models[i].params);
965 printf("Unfolding rate models:\n");
966 for (i=0; i<n_k_models; i++) {
967 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
968 for (j=0; j < k_models[i].num_params ; j++)
969 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
970 printf("\t\tdefaults: %s\n", k_models[i].params);
973 printf("Examples:\n");
974 printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded. Stop once there are no more folded domains.\n");
975 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);
976 printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded. Stop after 0.25 seconds\n");
977 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);
981 \subsection{Get options}
982 \label{sec.get_options}
986 void get_options(int argc, char **argv, double *pP, double *pDF,
987 double *pT_max, double *pDt_max, double *pV,
988 double *pX_step, state_t **pStop_state, environment_t *env,
989 int n_tension_models, tension_model_getopt_t *tension_models,
990 int n_k_models, k_model_getopt_t *k_models,
991 list_t **pState_list, list_t **pTransition_list)
993 char *prog_name = NULL;
994 char c, options[] = "q:t:d:P:F:v:x:T:C:s:N:k:DVh";
995 int tension_model=0, k_model=0, initial_state=-1, final_state=-1;
996 char *state_name=NULL;
997 int i, n, tension_group=0;
998 list_t *temp_list=NULL;
1001 void *model_params; /* for INIT_MODEL() */
1002 int num_param_args; /* for INIT_MODEL() */
1003 char **param_args; /* for INIT_MODEL() */
1004 extern char *optarg;
1005 extern int optind, optopt, opterr;
1010 for (i=0; i<n_tension_models; i++) {
1011 fprintf(stderr, "tension model %i %s:\t model's handler %p\t inverse handler %p\n",
1012 i, tension_models[i].name, tension_models[i].handler, tension_models[i].inverse_handler);
1013 assert(tension_models[i].handler == tension_handlers[i]);
1017 /* setup defaults */
1018 flags = FLAG_OUTPUT_TRANSITION_FORCES;
1019 prog_name = argv[0];
1020 *pStop_state = NULL;
1021 *pT_max = -1; /* s, -1 removes this limit */
1022 *pDt_max = 0.001; /* s */
1023 *pP = 1e-3; /* % pop per s */
1024 *pDF = 1e-12; /* N */
1025 *pV = 1e-6; /* m/s */
1026 *pX_step = 0.0; /* m */
1027 env->T = 300.0; /* K */
1028 *pState_list = NULL;
1029 *pTransition_list = NULL;
1031 while ((c=getopt(argc, argv, options)) != -1) {
1034 if (STRMATCH(optarg, "-")) {
1035 *pStop_state = NULL;
1037 temp_list = head(*pState_list);
1038 while (temp_list != NULL) {
1039 if (STRMATCH(S(temp_list)->name, optarg)) {
1040 *pStop_state = S(temp_list);
1043 temp_list = temp_list->next;
1047 case 't': *pT_max = safe_strtod(optarg, "-t"); break;
1048 case 'd': *pDt_max = safe_strtod(optarg, "-d"); break;
1049 case 'P': *pP = safe_strtod(optarg, "-P"); break;
1050 case 'F': *pDF = safe_strtod(optarg, "-F"); break;
1051 case 'v': *pV = safe_strtod(optarg, "-v"); break;
1052 case 'x': *pX_step = safe_strtod(optarg, "-x"); break;
1053 case 'T': env->T = safe_strtod(optarg, "-T"); break;
1054 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
1056 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
1057 assert(num_strings >= 2 && num_strings <= 4);
1059 state_name = strings[0];
1060 if (state_index(*pState_list, state_name) != -1) {
1061 fprintf(stderr, "Invalid option '%s', repeated state name '%s'\n", optarg, state_name);
1064 tension_model = index_tension_model(n_tension_models, tension_models, strings[1]);
1065 if (num_strings == 4)
1066 tension_group = safe_strtoi(strings[2], optarg);
1069 if (num_strings >= 3)
1070 tension_models[tension_model].params = strings[num_strings-1];
1072 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);
1074 INIT_MODEL(state_name, &tension_models[tension_model], tension_models[tension_model].params, model_params);
1075 push(pState_list, create_state(state_name,
1076 tension_models[tension_model].handler,
1077 tension_models[tension_model].inverse_handler,
1080 tension_models[tension_model].destructor,
1082 *pState_list = tail(*pState_list);
1084 free_parsed_list(num_strings, strings);
1087 n = safe_strtoi(optarg, "-N");
1089 fprintf(stderr, "%s:\tadding %d domains to %s\n", __FUNCTION__, n, S(*pState_list)->name);
1091 S(*pState_list)->num_domains += n;
1094 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
1095 assert(num_strings >= 3 && num_strings <= 4);
1096 initial_state = state_index(*pState_list, strings[0]);
1097 final_state = state_index(*pState_list, strings[1]);
1098 k_model = index_k_model(n_k_models, k_models, strings[2]);
1100 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);
1102 assert(initial_state != final_state);
1103 if (num_strings == 4)
1104 k_models[k_model].params = strings[3];
1105 INIT_MODEL(state_name, &k_models[k_model], k_models[k_model].params, model_params);
1106 push(pTransition_list, create_transition(list_element(*pState_list, initial_state),
1107 list_element(*pState_list, final_state),
1108 k_models[k_model].k,
1110 k_models[k_model].destructor), 1);
1111 free_parsed_list(num_strings, strings);
1114 print_state_graph(stdout, *pState_list, *pTransition_list);
1118 flags = FLAG_OUTPUT_FULL_CURVE;
1121 fprintf(stderr, "unrecognized option '%c'\n", optopt);
1122 /* fall through to default case */
1124 help(prog_name, *pP, *pDF, *pT_max, *pDt_max, *pV, *pX_step,
1125 *pStop_state, env, n_tension_models, tension_models,
1126 n_k_models, k_models, k_model, *pState_list);
1130 /* check the arguments */
1131 assert(*pState_list != NULL); /* no domains! */
1132 assert(*pP > 0.0); assert(*pP < 1.0);
1133 assert(*pDt_max > 0.0);
1135 assert(env->T > 0.0);
1137 crosslink_groups(*pState_list, *pTransition_list);
1143 <<index tension model>>=
1144 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
1147 for (i=0; i<n_models; i++)
1148 if (STRMATCH(models[i].name, name))
1150 if (i == n_models) {
1151 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
1159 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
1162 for (i=0; i<n_models; i++)
1163 if (STRMATCH(models[i].name, name))
1165 if (i == n_models) {
1166 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
1173 <<initialize model definition>>=
1174 /* requires int num_param_args and char **param_args in the current scope
1176 * INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
1177 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
1179 #define INIT_MODEL(role, model, param_string, param_pointer) \
1181 parse_list_string((param_string), SEP, '{', '}', \
1182 &num_param_args, ¶m_args); \
1183 if (num_param_args != (model)->num_params) { \
1185 "%s model %s expected %d params,\n", \
1186 (role), (model)->name, (model)->num_params); \
1188 "not the %d params in '%s'\n", \
1189 num_param_args, (param_string)); \
1190 assert(num_param_args == (model)->num_params); \
1192 if ((model)->creator) \
1193 param_pointer = (*(model)->creator)(param_args); \
1195 param_pointer = NULL; \
1196 free_parsed_list(num_param_args, param_args); \
1200 First we define some safe string parsing functions. Currently these
1201 abort on any irregularity, but in the future they could print error
1202 messages, etc. [[static]] because the functions are currently
1203 defined in each [[*.c]] file that uses them, but perhaps they should
1204 be moved to [[utils.h/c]] or some such instead\ldots
1207 static int safe_strtoi(const char *s, const char *description) {
1211 ret = strtol(s, &endp, 10);
1212 if (endp[0] != '\0') { //strlen(endp) > 0) {
1213 fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1214 endp, s, description);
1215 assert(1==0); //strlen(endp) == 0);
1217 fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1220 } else if (errno == ERANGE) {
1221 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1222 assert(errno != ERANGE);
1227 static double safe_strtod(const char *s, const char *description) {
1231 ret = strtod(s, &endp);
1232 if (endp[0] != '\0') { //strlen(endp) > 0) {
1233 fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1234 endp, s, description);
1235 assert(1==0); //strlen(endp) == 0);
1237 fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1240 } else if (errno == ERANGE) {
1241 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1242 assert(errno != ERANGE);
1251 \addcontentsline{toc}{section}{Appendicies}
1253 \section{sawsim.c details}
1257 The general layout of our simulation code is:
1268 We include [[math.h]], so don't forget to link to the libm with `-lm'.
1270 #include <assert.h> /* assert() */
1271 #include <stdlib.h> /* malloc(), free(), rand(), strto*() */
1272 #include <stdio.h> /* fprintf(), stdout */
1273 #include <string.h> /* strlen, strtok() */
1274 #include <math.h> /* exp(), M_PI, sqrt() */
1275 #include <sys/types.h> /* pid_t (returned by getpid()) */
1276 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
1277 #include <errno.h> /* errno, ERANGE (for safe_strto*()) */
1280 #include "tension_balance.h"
1281 #include "k_model.h"
1282 #include "tension_model.h"
1284 #include "accel_k.h"
1289 <<version definition>>
1290 <<flag definitions>>
1291 <<max/min definitions>>
1292 <<string comparison definition and globals>>
1293 <<initialize model definition>>
1294 <<get options definitions>>
1295 <<state definitions>>
1296 <<transition definitions>>
1305 <<create/destroy state>>
1306 <<destroy state list>>
1308 <<create/destroy transition>>
1309 <<destroy transition list>>
1310 <<print state graph>>
1312 <<simulation functions>>
1313 <<boilerplate functions>>
1316 <<boilerplate functions>>=
1318 <<get option functions>>
1324 srand(getpid()*time(NULL)); /* seed rand() */
1328 <<flag definitions>>=
1329 /* in octal b/c of prefixed '0' */
1330 #define FLAG_OUTPUT_FULL_CURVE 01
1331 #define FLAG_OUTPUT_TRANSITION_FORCES 02
1335 static unsigned long int flags = 0;
1338 \subsection{Utilities}
1341 <<max/min definitions>>=
1342 #define MAX(a,b) ((a)>(b) ? (a) : (b))
1343 #define MIN(a,b) ((a)<(b) ? (a) : (b))
1346 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
1347 <<string comparison definition and globals>>=
1348 // Check if two strings match, return 1 if they do
1349 static char *temp_string_A;
1350 static char *temp_string_B;
1351 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
1352 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
1353 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
1354 /* +1 to also compare the '\0' */
1357 We also define a macro for our [[check]] unit testing
1358 <<check relative error>>=
1359 #define CHECK_ERR(max_err, expected, received) \
1361 fail_unless( (received-expected)/expected < max_err, \
1362 "relative error %g >= %g in %s (Expected %g, received %g)", \
1363 (received-expected)/expected, max_err, #received, \
1364 expected, received); \
1365 fail_unless(-(received-expected)/expected < max_err, \
1366 "relative error %g >= %g in %s (Expected %g, received %g)", \
1367 -(received-expected)/expected, max_err, #received, \
1368 expected, received); \
1374 int happens(double probability)
1376 assert(probability >= 0.0); assert(probability <= 1.0);
1377 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*/
1382 \subsection{Adaptive timesteps}
1383 \label{sec.adaptive_dt_implementation}
1385 $F(x)$ increases with $x$, possibly exploding (e.g.\ in the worm-like
1386 chain model), so basing the timestep on the unfolding probability at
1387 the current tension is dangerous, and we need to search for a $dt$ for
1388 which $P_N(F(x+v*dt))<P_\text{max}$ (See Section
1389 \ref{sec.transition_rate} for a discussion of $P_N$). There are two
1390 cases to consider. In the most common, no domains have unfolded since
1391 the last step, and we expect the next step to be slightly shorter than
1392 the previous one. In the less common, domains did unfold in the last
1393 step, and we expect the next step to be considerably longer than the
1396 double search_dt(transition_t *transition,
1398 environment_t *env, double x,
1399 double max_prob, double max_dt, double v)
1400 { /* Returns a valid timestep dt in seconds for the given transition.
1401 * Takes a list of all states, a pointer env to the environmental data,
1402 * a starting separation x in m, a max_prob between 0 and 1,
1403 * max_dt in s, stretching velocity v in m/s.
1405 double F, k, P, dtCur, dtU, dtUCur, dtL, dt;
1407 /* get upper bound using the current position */
1408 F = find_tension(state_list, env, x, 0, 1); /* BUG. repeated calculation */
1409 //printf("Start with x = %g (F = %g)\n", x, F);
1410 k = accel_k(transition->k, F, env, transition->k_params);
1411 //printf("x %g\tF %g\tk %g\n", x, F, k);
1412 dtU = P1(max_prob, N_INIT(transition)) / k; /* upper bound on dt */
1414 //printf("overshot max_dt\n");
1417 if (v == 0) /* discrete stepping, so force is temporarily constant */
1420 /* set a lower bound on dt too */
1423 /* The dt determined above may produce illegitimate forces or ks.
1424 * Reduce the upper bound until we have valid ks. */
1426 F = find_tension(state_list, env, x+v*dt, 0, 1);
1427 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1430 F = find_tension(state_list, env, x+v*dt, 0, 1);
1432 //printf("Try for dt = %g (F = %g)\n", dt, F);
1433 k = accel_k(transition->k, F, env, transition->k_params);
1434 /* returned k may be -1.0 */
1435 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1436 while (k == -1.0) { /* reduce step until we hit a valid k */
1438 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1439 F = find_tension(state_list, env, x+v*dt, 0, 1);
1440 //printf("Try for dt = %g (F = %g)\n", dt, F);
1441 k = accel_k(transition->k, F, env, transition->k_params);
1442 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1444 assert(dtU > 1e-14); /* timestep to valid k too small */
1445 /* safe timestep back from x+dtU */
1446 dtUCur = P1(max_prob, N_INIT(transition)) / k;
1448 return dt; /* dtU is safe. */
1450 /* dtU wasn't safe, lets see what would be. */
1451 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1452 dt = (dtU + dtL) / 2.0;
1453 F = find_tension(state_list, env, x+v*dt, 0, 1);
1454 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1455 k = accel_k(transition->k, F, env, transition->k_params);
1456 dtCur = P1(max_prob, N_INIT(transition)) / k;
1457 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1458 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1460 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1461 dtU = dt; /* ... stepping out only dtCur would be safe */
1464 break; /* dtCur = dt */
1466 return MAX(dtUCur, dtL);
1471 Determine $dt$ for an array of potentially different folded domains.
1472 We apply [[search_dt()]] to each populated domain to find the maximum
1473 $dt$ the domain can handle. The final $dt$ is less than the
1474 individual domain $dt$s (to satisfy [[search_dt()]], see above). If
1475 $v>0$ (continuous pulling), $dt$ is also reduced to satisfy
1476 $|F(x+v*dt)-F(x)|<dF_\text{max}$, which limits errors due to our
1477 transition rate assumption that the force does not change appreciably
1478 over a single timestep.
1481 double determine_dt(list_t *state_list, list_t *transition_list,
1482 environment_t *env, double F, double x, double max_dF,
1483 double max_prob, double dt_max, double v, int transition)
1484 { /* Returns the timestep dt in seconds.
1485 * Takes the state and transition lists, a pointer to the environment,
1486 * the force F(x) in N, an extension x in m, max_dF in N,
1487 * max_prob between 0 and 1, dt_max in s, v in m/s, and transition in {0,1}*/
1488 double dt=dt_max, new_dt, F_new;
1489 assert(max_dF > 0.0);
1490 assert(max_prob > 0.0); assert(max_prob < 1.0);
1491 assert(dt_max > 0.0);
1492 assert(state_list != NULL);
1493 assert(transition_list != NULL);
1494 transition_list = head(transition_list);
1497 /* Ensure |dF| < max_dF */
1498 F_new = find_tension(state_list, env, x+v*dt, transition, 1);
1499 while (fabs(F_new - F) > max_dF) {
1501 F_new = find_tension(state_list, env, x+v*dt, transition, 1);
1505 /* Ensure all individual domains pass search_dt() */
1506 while (transition_list != NULL) {
1507 if (T(transition_list)->initial_state->num_domains > 0) {
1508 new_dt = search_dt(T(transition_list),state_list,env,x,max_prob,dt,v);
1509 dt = MIN(dt, new_dt);
1511 transition_list = transition_list->next;
1518 \subsection{State data}
1519 \label{sec.state_data}
1521 The domains exist in one of an array of possible states. Each state
1522 has a tension model which determines it's force/extension curve
1523 (possibly [[null]] for rigid states, see Appendix
1524 \ref{sec.null_tension}).
1526 Groups are, for example, for two domain states with WLC tensions; one
1527 with $p=0.4\U{nm}$, $L=10\U{nm}$ and the other with $p=0.4\U{nm}$,
1528 $L=28\U{nm}$. Since the persistence lengths are the same, you would
1529 like to add the contour lengths of all the domains in \emph{both}
1530 states, and plug that total contour length into a single WLC
1531 evaluation (see Section \ref{sec.find_tension}).
1532 <<state definitions>>=
1533 typedef struct state_struct {
1534 char *name; /* name string */
1535 one_dim_fn_t *tension_handler; /* tension handler */
1536 one_dim_fn_t *inverse_tension_handler; /* inverse tension handler */
1537 int tension_group; /* for grouping similar tension states */
1538 void *tension_params; /* pointer to tension parameters */
1539 destroy_data_func_t *destroy_tension;
1540 int num_domains; /* number of domains currently in state */
1541 /* cached lookups generated from the above data */
1542 int num_out_trans; /* length of out_trans array */
1543 struct transition_struct **out_trans; /* transitions leaving this state */
1544 int num_grouped_states; /* length of grouped states array */
1545 struct state_struct **grouped_states; /* states grouped with this one (doesn't include self) */
1548 /* get the state data for the current list node */
1549 #define S(list) ((state_t *)(list)->d)
1551 @ [[tension_params]] is a pointer to the parameters used by the
1552 handler function when determining the tension. [[destroy_tension]]
1553 points to a function for freeing [[tension_params]]. We it in the
1554 state structure so that [[destroy_state]] and will have an easier job
1557 [[create_]] and [[destroy_state]] are simple wrappers around
1558 [[malloc]] and [[free]].
1559 <<create/destroy state>>=
1560 state_t *create_state(char *name,
1561 one_dim_fn_t *tension_handler,
1562 one_dim_fn_t *inverse_tension_handler,
1564 void *tension_params,
1565 destroy_data_func_t *destroy_tension,
1568 state_t *ret = (state_t *)malloc(sizeof(state_t));
1569 assert(ret != NULL);
1570 /* make a copy of the name, so we aren't dependent on the original */
1571 ret->name = (char *)malloc(sizeof(char)*(strlen(name)+1));
1572 assert(ret->name != NULL);
1573 strcpy(ret->name, name); /* we know ret->name is long enough */
1575 ret->tension_handler = tension_handler;
1576 ret->inverse_tension_handler = inverse_tension_handler;
1577 ret->tension_group = tension_group;
1578 ret->tension_params = tension_params;
1579 ret->destroy_tension = destroy_tension;
1580 ret->num_domains = num_domains;
1582 ret->num_out_trans = 0;
1583 ret->out_trans = NULL;
1584 ret->num_grouped_states = 0;
1585 ret->grouped_states = NULL;
1588 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);
1593 void destroy_state(state_t *state)
1597 fprintf(stderr, "freeing state %s (%p)\n", state->name, state);
1601 if (state->destroy_tension)
1602 (*state->destroy_tension)(state->tension_params);
1603 if (state->out_trans)
1604 free(state->out_trans);
1605 if (state->grouped_states)
1606 free(state->grouped_states);
1613 We'll be storing the states in a list (see Appendix \ref{sec.lists}),
1614 so we also define a convenience function for cleaning up.
1615 <<destroy state list>>=
1616 void destroy_state_list(list_t *state_list)
1618 state_list = head(state_list);
1619 while (state_list != NULL)
1620 destroy_state((state_t *) pop(&state_list));
1625 We can also define a convenient way to get a state index from it's name.
1627 int state_index(list_t *state_list, char *name)
1630 state_list = head(state_list);
1631 while (state_list != NULL) {
1632 if (STRMATCH(S(state_list)->name, name)) return i;
1633 state_list = state_list->next;
1641 \subsection{Transition data}
1642 \label{sec.transition_data}
1644 Once you've created states (Sections \ref{sec.main},
1645 \ref{sec.model_selection}, and \ref{sec.state_data}), you need to
1646 create transitions between them.
1647 <<transition definitions>>=
1648 typedef struct transition_struct {
1649 state_t *initial_state;
1650 state_t *final_state;
1651 k_func_t *k; /* transition rate model */
1652 void *k_params; /* pointer to k parameters */
1653 destroy_data_func_t *destroy_k;
1656 /* get the transition data for the current list node */
1657 #define T(list) ((transition_t *)(list)->d)
1659 /* get the number of initial state population for a transition state */
1660 #define N_INIT(transition) ((transition)->initial_state->num_domains)
1664 @ [[k]] is a pointer to the function determining the transition rate
1665 for a given tension. [[k_params]] and [[destroy_k]] have similar
1666 roles with regards to [[k]] as [[tension_params]] and
1667 [[destroy_tension]] have with regards to [[tension_handler]] (see
1668 Section \ref{sec.state_data}).
1670 [[create_]] and [[destroy_transition]] are simple wrappers around
1671 [[malloc]] and [[free]].
1672 <<create/destroy transition>>=
1673 transition_t *create_transition(state_t *initial_state, state_t *final_state,
1676 destroy_data_func_t *destroy_k)
1678 transition_t *ret = (transition_t *)malloc(sizeof(transition_t));
1679 assert(ret != NULL);
1680 assert(initial_state != NULL);
1681 assert(final_state != NULL);
1682 ret->initial_state = initial_state;
1683 ret->final_state = final_state;
1685 ret->k_params = k_params;
1686 ret->destroy_k = destroy_k;
1688 fprintf(stderr, "generated transition %p from %s to %s, params %p\n", ret, ret->initial_state->name, ret->final_state->name, ret->k_params);
1693 void destroy_transition(transition_t *transition)
1697 fprintf(stderr, "freeing transition %p (%s to %s)\n", transition, transition->initial_state->name, transition->final_state->name);
1699 if (transition->destroy_k)
1700 (*transition->destroy_k)(transition->k_params);
1707 We'll be storing the transitions in a list (see Appendix
1708 \ref{sec.lists}), so we also define a convenience function for
1710 <<destroy transition list>>=
1711 void destroy_transition_list(list_t *transition_list)
1713 transition_list = head(transition_list);
1714 while (transition_list != NULL)
1715 destroy_transition((transition_t *) pop(&transition_list));
1720 \subsection{Graphviz output}
1721 \label{sec.graphviz_output}
1723 <<print state graph>>=
1724 void print_state_graph(FILE *file, list_t *state_list, list_t *transition_list)
1726 state_list = head(state_list);
1727 transition_list = head(transition_list);
1728 fprintf(file, "digraph states {\n node [shape=record,color=black,fillcolor=white,style=filled];\n\n");
1730 fprintf(file, " n%d [label=\"%s\"]\n", state_index(state_list, S(state_list)->name), S(state_list)->name);
1731 if (state_list->next == NULL) break;
1732 state_list = state_list->next;
1734 fprintf(file, "\n");
1735 while (transition_list != NULL) {
1736 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));
1737 transition_list = transition_list->next;
1739 fprintf(file, "}\n");
1743 \subsection{Domain model and group handling}
1745 <<group functions>>=
1746 <<crosslink groups>>
1747 <<get active groups>>
1750 Scan through a list of states and construct an array of all the
1752 <<get active groups>>=
1753 void get_active_groups(list_t *state_list,
1754 int *pNum_active_groups,
1755 one_dim_fn_t ***pPer_group_handlers,
1756 one_dim_fn_t ***pPer_group_inverse_handlers,
1757 void ***pPer_group_data /* pointer to array of pointers to tension_handler_data_t */)
1759 void **active_groups=NULL;
1760 one_dim_fn_t *handler, *old_handler, *inverse_handler, *old_inverse_handler;
1761 one_dim_fn_t **per_group_handlers=NULL, **per_group_inverse_handlers=NULL;
1762 void **per_group_data=NULL; /* array of pointers to tension_handler_data_t */
1763 int i, j, num_domains, group, old_group, num_active_groups=0;
1764 list_t *active_states_list=NULL;
1765 tension_handler_data_t *tdata=NULL;
1768 /* get all the active groups in a list */
1769 state_list = head(state_list);
1771 fprintf(stderr, "%s:\t", __FUNCTION__);
1772 list_print(stderr, state_list, "states");
1774 while (state_list != NULL) {
1775 state = S(state_list);
1776 handler = state->tension_handler;
1777 inverse_handler = state->inverse_tension_handler;
1778 group = state->tension_group;
1779 num_domains = state->num_domains;
1780 if (list_index(active_states_list, state) == -1) {
1781 /* we haven't added this state (or it's associated group) yet */
1782 if (num_domains > 0 && handler != NULL) { /* add the state */
1783 num_active_groups++;
1784 push(&active_states_list, state, 1);
1786 fprintf(stderr, "%s:\tactive %s state (%p) with %d domain(s)\n", __FUNCTION__, state->name, state, state->num_domains);
1788 for (i=0; i < state->num_grouped_states; i++) {
1789 if (state->grouped_states[i]->num_domains > 0) {
1790 /* add this related (and active) state now */
1791 assert(state->grouped_states[i]->tension_handler == handler);
1792 assert(state->grouped_states[i]->tension_group == group);
1793 push(&active_states_list, state->grouped_states[i], 1);
1795 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);
1799 } /* else empty state or NULL handler */
1800 } /* else state already added */
1801 state_list = state_list->next;
1804 fprintf(stderr, "%s:\t", __FUNCTION__);
1805 list_print(stderr, active_states_list, "active states");
1808 assert(num_active_groups <= list_length(active_states_list));
1809 assert(num_active_groups > 0);
1811 /* allocate the arrays we'll be returning */
1812 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1813 assert(per_group_handlers != NULL);
1814 per_group_inverse_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1815 assert(per_group_inverse_handlers != NULL);
1816 per_group_data = (void **)malloc(sizeof(void *)*num_active_groups);
1817 assert(per_group_data != NULL);
1819 fprintf(stderr, "%s:\tallocating per group data %p:", __FUNCTION__, per_group_data);
1821 for (i=0; i<num_active_groups; i++) { /* we might want to recycle any group data structs allocated in the previous round */
1822 per_group_data[i] = malloc(sizeof(tension_handler_data_t));
1823 assert(per_group_data[i] != NULL);
1825 fprintf(stderr, " %d,%p", i, per_group_data[i]);
1829 fprintf(stderr, "\n");
1832 /* populate the arrays we'll be returning */
1833 active_states_list = head(active_states_list);
1835 old_inverse_handler = NULL;
1836 j = -1; /* j is the current active _group_ index */
1837 while (active_states_list != NULL) {
1838 state = (state_t *) pop(&active_states_list);
1839 handler = state->tension_handler;
1840 inverse_handler = state->inverse_tension_handler;
1841 group = state->tension_group;
1842 num_domains = state->num_domains;
1843 if (handler != old_handler || inverse_handler != old_inverse_handler || group != old_group) {
1844 j++; /* move to the next group */
1845 tdata = (tension_handler_data_t *) per_group_data[j];
1846 per_group_handlers[j] = handler;
1847 per_group_inverse_handlers[j] = inverse_handler;
1848 tdata->group_tension_params = NULL;
1850 tdata->persist = NULL; /* we might want to grab any appropriate data from the old persist here... */
1853 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);
1855 for (i=0; i<num_domains; i++) { /* add num_domains copies of the tension params */
1856 push(&tdata->group_tension_params, state->tension_params, 1);
1859 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);
1861 old_handler = handler;
1862 old_inverse_handler = inverse_handler;
1866 /* free old groups */
1867 if (*pPer_group_handlers != NULL)
1868 free(*pPer_group_handlers);
1869 if (*pPer_group_inverse_handlers != NULL)
1870 free(*pPer_group_inverse_handlers);
1871 if (*pPer_group_data != NULL) {
1872 for (i=0; i<*pNum_active_groups; i++)
1873 free((*pPer_group_data)[i]);
1874 free(*pPer_group_data);
1877 /* set new groups */
1879 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]);
1881 *pNum_active_groups = num_active_groups;
1882 *pPer_group_handlers = per_group_handlers;
1883 *pPer_group_inverse_handlers = per_group_inverse_handlers;
1884 *pPer_group_data = per_group_data;
1888 <<crosslink groups>>=
1889 void crosslink_groups(list_t *state_groups, list_t *transition_list)
1891 list_t *list, *out_trans=NULL, *associates=NULL;
1892 one_dim_fn_t *handler;
1895 state_groups = head(state_groups);
1896 while (state_groups != NULL) {
1897 /* find transitions leaving this state */
1899 fprintf(stderr, "finding transitions leaving %s\n", S(state_groups)->name);
1901 list = head(transition_list);
1902 while (list != NULL) {
1903 if (T(list)->initial_state == S(state_groups) && list_index(associates, T(list)->initial_state) != -1) {
1904 push(&out_trans, T(list), 1);
1908 n = list_length(out_trans);
1909 S(state_groups)->num_out_trans = n;
1910 S(state_groups)->out_trans = (transition_t **)malloc(sizeof(transition_t *)*n);
1911 assert(S(state_groups)->out_trans != NULL);
1912 for (i=0; i<n; i++) {
1913 S(state_groups)->out_trans[i] = (transition_t *)pop(&out_trans);
1916 /* find states grouped with this state */
1918 fprintf(stderr, "finding states grouped with %s\n", S(state_groups)->name);
1920 handler = S(state_groups)->tension_handler;
1921 group = S(state_groups)->tension_group;
1922 list = head(state_groups);
1923 while (list != NULL) {
1924 if (S(list)->tension_handler == handler && S(list)->tension_group == group && list != state_groups) {
1925 push(&associates, S(list), 1);
1929 n = list_length(associates);
1930 S(state_groups)->num_grouped_states = n;
1931 S(state_groups)->grouped_states = (state_t **)malloc(sizeof(state_t *)*n);
1932 assert(S(state_groups)->grouped_states != NULL);
1933 for (i=0; i<n; i++) {
1934 S(state_groups)->grouped_states[i] = (state_t *)pop(&associates);
1936 state_groups = state_groups->next;
1942 \section{String parsing}
1944 For handling command line arguments, we need parse delimited strings (e.g.\ [["param1,param2,param3..."]]).
1945 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1946 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1947 need to take care of parsing those parameters themselves.
1948 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1954 <<parse definitions>>
1955 <<parse declarations>>
1959 <<parse module makefile lines>>=
1960 NW_SAWSIM_MODS += parse
1961 CHECK_BINS += check_parse
1965 <<parse definitions>>=
1966 #define SEP ',' /* argument separator character */
1969 <<parse declarations>>=
1970 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1971 int *num, char ***string_array);
1972 extern void free_parsed_list(int num, char **string_array);
1975 [[parse_list_string]] allocates memory, don't forget to free it
1976 afterward with [[free_parsed_list]]. It does not alter the original.
1978 The string may start off with a [[deeper]] character
1979 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1980 [[x,y]], where the pointer is one character in on the copied string.
1981 However, when we go to free the memory, we need a pointer to the
1982 beginning of the string. In order to accommodate this for a string
1983 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1984 the first $N$ elements point to the separated arguments, and let the
1985 last element point to the start of the copied string regardless of
1987 <<parse delimited list string functions>>=
1988 /* TODO, split out into parse.hc */
1989 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1992 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1993 if (string[i] == deeper) {depth++;}
1994 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
2000 void parse_list_string(char *string, char sep, char deeper, char shallower,
2001 int *num, char ***string_array)
2003 char *str=NULL, **ret=NULL;
2005 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
2007 *string_array = NULL;
2010 /* make a copy of the string, so we don't change the original */
2011 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
2012 assert(str != NULL);
2013 strcpy(str, string); /* we know str is long enough */
2014 /* count the number of regions, so we can allocate pointers to them */
2017 n++; i++; /* move on to next argument */
2018 i += next_delim_index(str+i, sep, deeper, shallower);
2019 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
2020 } while (str[i] != '\0');
2021 ret = (char **)malloc(sizeof(char *)*(n+1));
2022 assert(ret != NULL);
2023 /* replace the separators with '\0' & assign pointers */
2024 ret[n] = str; /* point to the front of the copied string */
2027 for(i=1; i<n; i++) {
2028 j += next_delim_index(str+j, sep, deeper, shallower);
2030 ret[i] = str+j; /* point to the separated arguments */
2032 /* strip off bounding braces, if any */
2033 for(i=0; i<n; i++) {
2034 if (ret[i][0] == deeper) {
2038 if (ret[i][strlen(ret[i])-1] == shallower)
2039 ret[i][strlen(ret[i])-1] = '\0';
2042 *string_array = ret;
2045 void free_parsed_list(int num, char **string_array)
2048 /* frees all items, since they were allocated as a single string*/
2049 free(string_array[num]);
2050 /* and free the array of pointers */
2056 \subsection{Parsing implementation}
2062 <<parse delimited list string functions>>
2066 #include <assert.h> /* assert() */
2067 #include <stdlib.h> /* NULL */
2068 #include <stdio.h> /* fprintf(), stdout *//*!!*/
2069 #include <string.h> /* strlen() */
2073 \subsection{Parsing unit tests}
2075 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2078 <<parse check includes>>
2079 <<string comparison definition and globals>>
2080 <<check relative error>>
2081 <<parse test suite>>
2082 <<main check program>>
2085 <<parse check includes>>=
2086 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2087 #include <stdio.h> /* printf() */
2088 #include <assert.h> /* assert() */
2089 #include <string.h> /* strlen() */
2094 <<parse test suite>>=
2095 <<parse list string tests>>
2096 <<parse suite definition>>
2099 <<parse suite definition>>=
2100 Suite *test_suite (void)
2102 Suite *s = suite_create ("parse");
2103 <<parse list string test case defs>>
2105 <<parse list string test case adds>>
2110 <<parse list string tests>>=
2113 START_TEST(test_next_delim_index)
2115 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
2116 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
2117 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
2118 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
2119 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
2123 START_TEST(test_parse_list_null)
2127 parse_list_string(NULL, SEP, '{', '}',
2128 &num_param_args, ¶m_args);
2129 fail_unless(num_param_args == 0, NULL);
2130 fail_unless(param_args == NULL, NULL);
2133 START_TEST(test_parse_list_single_simple)
2137 parse_list_string("arg", SEP, '{', '}',
2138 &num_param_args, ¶m_args);
2139 fail_unless(num_param_args == 1, NULL);
2140 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
2143 START_TEST(test_parse_list_single_compound)
2147 parse_list_string("{x,y,z}", SEP, '{', '}',
2148 &num_param_args, ¶m_args);
2149 fail_unless(num_param_args == 1, NULL);
2150 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
2153 START_TEST(test_parse_list_double_simple)
2157 parse_list_string("abc,def", SEP, '{', '}',
2158 &num_param_args, ¶m_args);
2159 fail_unless(num_param_args == 2, NULL);
2160 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2161 fail_unless(STRMATCH(param_args[1],"def"), NULL);
2164 START_TEST(test_parse_list_compound)
2168 parse_list_string("abc,{def,ghi}", SEP, '{', '}',
2169 &num_param_args, ¶m_args);
2170 fail_unless(num_param_args == 2, NULL);
2171 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2172 fail_unless(STRMATCH(param_args[1],"def,ghi"), NULL);
2176 <<parse list string test case defs>>=
2177 TCase *tc_parse_list_string = tcase_create("parse list string");
2179 <<parse list string test case adds>>=
2180 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
2181 tcase_add_test(tc_parse_list_string, test_parse_list_null);
2182 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
2183 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
2184 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
2185 tcase_add_test(tc_parse_list_string, test_parse_list_compound);
2186 suite_add_tcase(s, tc_parse_list_string);
2190 \section{Unit tests}
2192 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2199 <<check relative error>>
2202 <<main check program>>
2213 <<determine dt tests>>
2215 <<does domain transition tests>>
2216 <<randomly unfold domains tests>>
2217 <<suite definition>>
2220 <<suite definition>>=
2221 Suite *test_suite (void)
2223 Suite *s = suite_create ("sawsim");
2224 <<determine dt test case defs>>
2225 <<happens test case defs>>
2226 <<does domain transition test case defs>>
2227 <<randomly unfold domains test case defs>>
2229 <<determine dt test case adds>>
2230 <<happens test case adds>>
2231 <<does domain transition test case adds>>
2232 <<randomly unfold domains test case adds>>
2235 tcase_add_checked_fixture(tc_strip_address,
2236 setup_strip_address,
2237 teardown_strip_address);
2243 <<main check program>>=
2247 Suite *s = test_suite();
2248 SRunner *sr = srunner_create(s);
2249 srunner_run_all(sr, CK_ENV);
2250 number_failed = srunner_ntests_failed(sr);
2252 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
2257 \subsection{Model tests}
2259 Check the searching with [[linear_k]].
2260 Check overwhelming force treatment with the heavyside-esque step [[?]].
2262 Hmm, I'm not really sure what I was doing here...
2263 <<determine dt tests>>=
2264 double linear_k(double F, environment_t *env, void *params)
2266 double Fnot = *(double *)params;
2270 START_TEST(test_determine_dt_linear_k)
2273 transition_t unfold={0};
2274 environment_t env = {0};
2275 double F[]={0,1,2,3,4,5,6};
2276 double dt_max=3.0, Fnot=3.0, x=1.0, max_p=0.1, max_dF=0.1, v=1.0;
2277 list_t *state_list=NULL, *transition_list=NULL;
2280 folded.tension_handler = &hooke_handler;
2281 folded.num_domains = 1;
2282 unfold.initial_state = &folded;
2283 unfold.k = linear_k;
2284 unfold.k_params = &Fnot;
2285 push(&state_list, &folded, 1);
2286 push(&transition_list, &unfold, 1);
2287 for( i=0; i < sizeof(F)/sizeof(double); i++) {
2288 //fail_unless(determine_dt(state_list, transition_list, &env, x, max_p, max_dF, dt_max, v)==1, NULL);
2293 <<determine dt test case defs>>=
2294 TCase *tc_determine_dt = tcase_create("Determine dt");
2296 <<determine dt test case adds>>=
2297 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
2298 suite_add_tcase(s, tc_determine_dt);
2303 <<happens test case defs>>=
2305 <<happens test case adds>>=
2308 <<does domain transition tests>>=
2310 <<does domain transition test case defs>>=
2312 <<does domain transition test case adds>>=
2315 <<randomly unfold domains tests>>=
2317 <<randomly unfold domains test case defs>>=
2319 <<randomly unfold domains test case adds>>=
2323 \section{Balancing group extensions}
2324 \label{sec.tension_balance}
2326 For a given total extension $x$ (of the piezo), the various domain
2327 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
2328 amounts, and we need to tweak the portion that each extends to balance
2329 the tension amongst the active groups. Since the tension balancing is
2330 separable from the bulk of the simulation, we break it out into a
2331 separate module. The interface is defined in [[tension_balance.h]],
2332 the implementation in [[tension_balance.c]], and the unit testing in
2333 [[check_tension_balance.c]]
2335 <<tension-balance.h>>=
2336 #ifndef TENSION_BALANCE_H
2337 #define TENSION_BALANCE_H
2339 <<tension balance function declaration>>
2340 #endif /* TENSION_BALANCE_H */
2343 <<tension balance functions>>=
2344 <<one dimensional bracket>>
2345 <<one dimensional solve>>
2347 <<group stiffness function>>
2348 <<chain stiffness function>>
2349 <<full chain stiffness function>>
2350 <<tension balance function>>
2353 <<tension balance module makefile lines>>=
2354 NW_SAWSIM_MODS += tension_balance
2355 CHECK_BINS += check_tension_balance
2356 check_tension_balance_MODS =
2359 The entire force balancing problem reduces to a solving two nested
2360 one-dimensional root-finding problems. First we define the one
2361 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
2362 populated group). $x(x_0)$ is also strictly monotonically increasing,
2363 so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
2364 stores the last successful value of $x$, since we expect to be taking
2365 small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
2366 indicates that the groups have changed.
2367 <<tension balance function declaration>>=
2368 double tension_balance(int num_tension_groups,
2369 one_dim_fn_t **tension_handlers,
2370 one_dim_fn_t **inverse_tension_handlers,
2371 void **params, /* array of pointers to tension_handler_data_t */
2376 <<tension balance function>>=
2377 double tension_balance(int num_tension_groups,
2378 one_dim_fn_t **tension_handlers,
2379 one_dim_fn_t **inverse_tension_handlers,
2380 void **params, /* array of pointers to tension_handler_data_t */
2385 static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL,NULL};
2386 double F, xo_guess, xo, lb, ub;
2387 double min_relx=1e-6, min_rely=1e-6;
2388 int max_steps=200, i;
2390 assert(num_tension_groups > 0);
2391 assert(tension_handlers != NULL);
2392 assert(params != NULL);
2393 assert(num_tension_groups > 0);
2395 if (last_x == -1 || x_xo_data.xi == NULL) { /* new group setup, reset x_xo_data */
2396 /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
2397 if (x_xo_data.xi != NULL)
2399 /* construct new x_xo_data */
2400 x_xo_data.n_groups = num_tension_groups;
2401 x_xo_data.tension_handler = tension_handlers;
2402 x_xo_data.inverse_tension_handler = inverse_tension_handlers;
2403 x_xo_data.handler_data = params;
2405 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);
2407 x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
2408 assert(x_xo_data.xi != NULL);
2409 for (i=0; i<num_tension_groups; i++)
2410 x_xo_data.xi[i] = last_x;
2413 if (num_tension_groups == 1) { /* only one group, no balancing required */
2415 x_xo_data.xi[0] = xo;
2417 fprintf(stderr, "%s:\tonly one active group, no balancing required\n", __FUNCTION__);
2421 fprintf(stderr, "%s:\tbalancing for x = %g with", __FUNCTION__, x);
2422 for(i=0; i<num_tension_groups; i++) fprintf(stderr, "\tx[%d]=%g", i, x_xo_data.xi[i]);
2423 fprintf(stderr, "\n");
2425 if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
2426 if (x == 0) xo_guess = length_scale;
2427 else xo_guess = x/num_tension_groups;
2429 fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
2431 oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
2432 } else { /* work off the last balanced point */
2434 fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
2437 lb = x_xo_data.xi[0];
2438 ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */
2439 } else if (x < last_x) {
2440 lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
2441 ub = x_xo_data.xi[0];
2442 } else { /* x == last_x */
2444 fprintf(stderr, "not moving\n");
2446 lb = x_xo_data.xi[0];
2447 ub = x_xo_data.xi[0];
2451 fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
2452 __FUNCTION__, x, lb, ub);
2454 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
2455 F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
2457 if (num_tension_groups > 1) {
2458 fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
2459 for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
2460 fprintf(stderr, "\n");
2465 F = (*tension_handlers[0])(xo, params[0]);
2467 if (stiffness != NULL) {
2468 *stiffness = chain_stiffness(&x_xo_data, min_relx);
2469 if (*stiffness == 0.0) { /* re-calculate based on full chain */
2470 *stiffness = full_chain_stiffness(num_tension_groups, tension_handlers,
2471 inverse_tension_handlers, params,
2472 last_x, x, min_relx, F);
2480 <<tension balance internal definitions>>=
2481 <<x of x0 definitions>>
2484 <<x of x0 definitions>>=
2485 typedef struct x_of_xo_data_struct {
2487 one_dim_fn_t **tension_handler; /* array of fn pointers */
2488 one_dim_fn_t **inverse_tension_handler; /* inverse fn pointers */
2489 void **handler_data; /* array of pointers to tension_handler_data_t */
2490 double *xi; /* array of group extensions */
2495 double x_of_xo(double xo, void *pdata)
2497 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
2498 double F, x=0, xi, xi_guess, lb, ub, slope;
2500 double min_relx=1e-6, min_rely=1e-6;
2502 assert(data->n_groups > 0);
2505 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);
2508 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
2510 fprintf(stderr, "%s:\tfinding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
2514 if (data->xi) data->xi[0] = xo;
2515 for (i=1; i < data->n_groups; i++) {
2516 if (data->inverse_tension_handler[i] != NULL) {
2517 xi = (*data->inverse_tension_handler[i])(F, data->handler_data[i]);
2518 } else { /* invert numerically */
2519 if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
2520 else xi_guess = data->xi[i];
2522 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]);
2524 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
2525 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
2530 fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
2534 fprintf(stderr, "%s:\tfound x(x0=%g) = %g\n", __FUNCTION__, xo, x);
2540 Determine the stiffness of a single tension group by taking a small
2541 step [[dx]] back from the position [[x]] for which we want the
2542 stiffness. The touchy part is determining a good [[dx]] and ensuring
2543 the whole interval is on [[x>=0]].
2544 <<group stiffness function>>=
2545 double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double x, double relx)
2547 double dx, xl, F, Fl, stiffness;
2549 assert(relx > 0 && relx < 1);
2550 if (x == 0 || relx == 0) {
2551 dx = 1e-12; /* HACK, how to get length scale? */
2561 F = (*tension_handler)(x, handler_data);
2562 Fl = (*tension_handler)(xl, handler_data);
2563 stiffness = (F-Fl)/dx;
2564 assert(stiffness >= 0);
2569 Attempt to calculate the chain stiffness by adding group stiffnesses
2570 as springs in series. In the case where a group has undefined
2571 stiffness (e.g. piston tension, Section \ref{sec.piston_tension}),
2572 this algorithm breaks down. In that case, [[tension_balance()]]
2573 (\ref{sec.tension_balance}) should use [[full_chain_stiffness()]]
2574 which uses the full chain's $F(x)$ rather than that of the individual
2575 domains'. We attempt the series approach first because, lacking the
2576 numerical inversion inside [[tension_balance()]], it is both faster
2577 and more accurate than the full-chain derivative.
2578 <<chain stiffness function>>=
2579 double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx)
2581 double stiffness, gstiff;
2583 /* add group stiffnesses in series */
2584 for (i=0; i < x_xo_data->n_groups; i++) {
2586 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);
2589 gstiff = group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx);
2592 stiffness += 1.0/gstiff;
2594 stiffness = 1.0 / stiffness;
2600 Determine the chain stiffness using only [[tension_balance()]]. This
2601 works around difficulties with tension models that have undefined
2602 stiffness (see discussion for [[chain_stiffness()]] above).
2603 <<full chain stiffness function>>=
2604 double full_chain_stiffness(int num_tension_groups,
2605 one_dim_fn_t **tension_handlers,
2606 one_dim_fn_t **inverse_tension_handlers,
2607 void **params, /* array of pointers to tension_handler_data_t */
2611 double F /* already evaluated F(x) */)
2613 double dx, xl, Fl, stiffness;
2616 assert(relx > 0 && relx < 1);
2618 /* Other option for dx that we ignore because of our poor tension_balance()
2619 * resolution (i.e. bad slopes if you zoom in too close):
2620 * if (last_x != -1.0 && last_x != x) // last_x limited
2621 * dx fabs(x-last_x);
2624 * if (1==1) { // constant max-value limited
2626 * dx = MIN(dx, dx_p);
2628 * if (x != 0 && relx != 0) { // relx limited
2630 * dx = MIN(dx, dx_p);
2632 * TODO, determine which of (min_relx, min_rely, max_steps) is actually
2633 * limiting tension_balance.
2635 dx = 1e-12; /* HACK, how to get length scale? */
2639 Fl = tension_balance(num_tension_groups, tension_handlers,
2640 inverse_tension_handlers, params,
2646 F = tension_balance(num_tension_groups, tension_handlers,
2647 inverse_tension_handlers, params,
2650 stiffness = (F-Fl)/dx;
2651 assert(stiffness >= 0);
2657 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
2658 number of steps for monotonically increasing $f(x)$. Simple
2659 bisection, so it's robust and converges fairly quickly. You can set
2660 the maximum number of steps to take, but if you have enough processor
2661 power, it's better to limit your precision with the relative limits
2662 [[min_relx]] and [[y]]. These ensure that the width of the bracket is
2663 small on the length scale of it's larger side. Note that \emph{both}
2664 relative limits must be satisfied for the search to stop.
2665 <<one dimensional function definition>>=
2666 /* equivalent to gsl_func_t */
2667 typedef double one_dim_fn_t(double x, void *params);
2669 <<one dimensional solve declaration>>=
2670 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2671 double min_relx, double min_rely, int max_steps,
2675 <<one dimensional solve>>=
2676 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2677 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2678 double min_relx, double min_rely, int max_steps,
2681 double yx, ylb, yub, x;
2684 ylb = (*f)(lb, params);
2685 yub = (*f)(ub, params);
2687 /* check some simple cases */
2689 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bracketed */
2690 else return lb; /* any x on the range [lb, ub] would work */
2692 if (ylb == y) { x=lb; yx=ylb; goto end; }
2693 if (yub == y) { x=ub; yx=yub; goto end; }
2696 fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2698 assert(ylb < y); assert(yub > y);
2699 for (i=0; i<max_steps; i++) {
2701 yx = (*f)(x, params);
2703 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);
2705 if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2706 else if (yx > y) { ub=x; yub=yx; }
2707 else /* yx < y */ { lb=x; ylb=yx; }
2712 fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2718 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2719 Generate bracketing $x$ values through bisection or exponential growth.
2720 <<one dimensional bracket declaration>>=
2721 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2724 <<one dimensional bracket>>=
2725 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2727 double yx, step, x=xguess;
2730 fprintf(stderr, "%s:\tbracketing %g (params %p)\n", __FUNCTION__, y, params);
2732 yx = (*f)(x, params);
2734 fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2741 fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2745 if (x == 0) x = length_scale/2.0; /* HACK */
2748 yx = (*f)(x, params);
2750 fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2755 fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2758 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2762 \subsection{Balancing implementation}
2764 <<tension-balance.c>>=
2767 <<tension balance includes>>
2768 #include "tension_balance.h"
2770 double length_scale = 1e-10; /* HACK */
2772 <<tension balance internal definitions>>
2773 <<tension balance functions>>
2776 <<tension balance includes>>=
2777 #include <assert.h> /* assert() */
2778 #include <stdlib.h> /* NULL */
2779 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2780 #include <stdio.h> /* fprintf(), stdout */
2784 \subsection{Balancing unit tests}
2786 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2787 <<check-tension-balance.c>>=
2789 <<tension balance check includes>>
2790 <<tension balance test suite>>
2791 <<main check program>>
2794 <<tension balance check includes>>=
2795 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2796 #include <assert.h> /* assert() */
2799 #include "tension_balance.h"
2802 <<tension balance test suite>>=
2803 <<tension balance function tests>>
2804 <<tension balance suite definition>>
2807 <<tension balance suite definition>>=
2808 Suite *test_suite (void)
2810 Suite *s = suite_create ("tension balance");
2811 <<tension balance function test case defs>>
2813 <<tension balance function test case adds>>
2818 <<tension balance function tests>>=
2819 <<check relative error>>
2821 double hooke(double x, void *pK)
2824 return *((double*)pK) * x;
2827 START_TEST(test_single_function)
2829 double k=5, x=3, last_x=2.0, F;
2830 one_dim_fn_t *handlers[] = {&hooke};
2831 one_dim_fn_t *inverse_handlers[] = {NULL};
2832 void *data[] = {&k};
2833 F = tension_balance(1, handlers, inverse_handlers, data, last_x, x, NULL);
2834 fail_unless(F = k*x, NULL);
2839 We can also test balancing two springs with different spring constants.
2840 The analytic solution for a general number of springs is given in Appendix \ref{sec.math_hooke}.
2841 <<tension balance function tests>>=
2842 START_TEST(test_double_hooke)
2844 double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
2845 one_dim_fn_t *handlers[] = {&hooke, &hooke};
2846 one_dim_fn_t *inverse_handlers[] = {NULL, NULL};
2847 void *data[] = {&k1, &k2};
2848 F = tension_balance(2, handlers, inverse_handlers, data, last_x, x, NULL);
2852 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2853 //CHECK_ERR(1e-6, x1e, xi[0]);
2854 //CHECK_ERR(1e-6, x2e, xi[1]);
2855 CHECK_ERR(1e-6, Fe, F);
2859 TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
2861 <<tension balance function test case defs>>=
2862 TCase *tc_tbfunc = tcase_create("tension balance function");
2865 <<tension balance function test case adds>>=
2866 tcase_add_test(tc_tbfunc, test_single_function);
2867 tcase_add_test(tc_tbfunc, test_double_hooke);
2868 suite_add_tcase(s, tc_tbfunc);
2874 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2875 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2876 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2882 <<list definitions>>
2883 <<list declarations>>
2887 <<list declarations>>=
2888 <<head and tail declarations>>
2889 <<list length declaration>>
2890 <<push declaration>>
2892 <<list destroy declaration>>
2893 <<list to array declaration>>
2894 <<move declaration>>
2895 <<sort declaration>>
2896 <<list index declaration>>
2897 <<list element declaration>>
2898 <<unique declaration>>
2899 <<list print declaration>>
2903 <<create and destroy node>>
2918 <<list module makefile lines>>=
2919 NW_SAWSIM_MODS += list
2920 CHECK_BINS += check_list
2924 <<list definitions>>=
2925 typedef struct list_struct {
2926 struct list_struct *next;
2927 struct list_struct *prev;
2931 <<comparison function typedef>>
2934 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2935 <<head and tail declarations>>=
2936 list_t *head(list_t *list);
2937 list_t *tail(list_t *list);
2940 list_t *head(list_t *list)
2944 while (list->prev != NULL) {
2950 list_t *tail(list_t *list)
2954 while (list->next != NULL) {
2961 <<list length declaration>>=
2962 int list_length(list_t *list);
2965 int list_length(list_t *list)
2972 while (list->next != NULL) {
2980 [[push]] inserts a node relative to the current position in the list
2981 without changing the current position.
2982 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2983 so we set it to that of the pushed domain.
2984 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2985 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2986 <<push declaration>>=
2987 void push(list_t **pList, void *data, int below);
2990 void push(list_t **pList, void *data, int below)
2992 list_t *list, *new_node;
2993 assert(pList != NULL);
2995 new_node = create_node(data);
2998 else if (below==0) { /* insert above */
2999 if (list->prev != NULL)
3000 list->prev->next = new_node;
3001 new_node->prev = list->prev;
3002 list->prev = new_node;
3003 new_node->next = list;
3004 } else { /* insert below */
3005 if (list->next != NULL)
3006 list->next->prev = new_node;
3007 new_node->next = list->next;
3008 list->next = new_node;
3009 new_node->prev = list;
3014 [[pop]] removes the current domain node, moving the current position
3015 to the node after it, unless that node is [[NULL]], in which case move
3016 the current position to the node before it.
3017 <<pop declaration>>=
3018 void *pop(list_t **pList);
3021 void *pop(list_t **pList)
3023 list_t *list, *popped;
3025 assert(pList != NULL);
3027 assert(list != NULL); /* not an empty list */
3029 /* bypass the popped node */
3030 if (list->prev != NULL)
3031 list->prev->next = list->next;
3032 if (list->next != NULL) {
3033 list->next->prev = list->prev;
3034 *pList = list->next;
3036 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
3038 destroy_node(popped);
3043 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
3044 If the cleanup function is [[NULL]], no cleanup function is called.
3045 The nodes are not popped in any particular order.
3046 <<list destroy declaration>>=
3047 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
3050 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
3054 assert(pList != NULL);
3057 assert(list != NULL); /* not an empty list */
3058 while (list != NULL) {
3060 if (destroy != NULL)
3066 [[list_to_array]] converts a list to an array of pointers, preserving order.
3067 <<list to array declaration>>=
3068 void list_to_array(list_t *list, int *length, void ***array);
3071 void list_to_array(list_t *list, int *pLength, void ***pArray)
3075 assert(list != NULL);
3076 assert(pLength != NULL);
3077 assert(pArray != NULL);
3079 length = list_length(list);
3080 array = (void **)malloc(sizeof(void *)*length);
3081 assert(array != NULL);
3084 while (list != NULL) {
3085 array[i++] = list->d;
3093 [[move]] moves the current element along the list in either direction.
3094 <<move declaration>>=
3095 void move(list_t **pList, int downstream);
3098 void move(list_t **pList, int downstream)
3100 list_t *list, *flipper;
3102 assert(pList != NULL);
3103 list = *pList; /* a>B>c>d */
3104 assert(list != NULL); /* not an empty list */
3105 if (downstream == 0)
3106 flipper = list->prev; /* flipper is A */
3108 flipper = list->next; /* flipper is C */
3109 /* ensure there is somewhere to go in the direction we're heading */
3110 assert(flipper != NULL);
3111 /* remove the current element */
3112 data = pop(&list); /* data=B, a>C>d */
3113 /* and put it back in in the correct spot */
3115 if (downstream == 0) {
3116 push(&list, data, 0); /* b>A>c>d */
3117 list = list->prev; /* B>a>c>d */
3119 push(&list, data, 1); /* a>C>b>d */
3120 list = list->next; /* a>c>B>d */
3126 [[sort]] sorts a list from smallest to largest according to a user
3127 specified comparison.
3128 <<comparison function typedef>>=
3129 typedef int (comparison_fn_t)(void *A, void *B);
3132 <<int comparison function>>=
3133 int int_comparison(void *A, void *B)
3138 if (a > b) return 1;
3139 else if (a == b) return 0;
3144 <<sort declaration>>=
3145 void sort(list_t **pList, comparison_fn_t *comp);
3147 Here we use the bubble sort algorithm.
3149 void sort(list_t **pList, comparison_fn_t *comp)
3153 assert(pList != NULL);
3155 assert(list != NULL); /* not an empty list */
3156 while (stable == 0) {
3159 while (list->next != NULL) {
3160 if ((*comp)(list->d, list->next->d) > 0) {
3162 move(&list, 1 /* downstream */);
3172 [[list_index]] finds the location of [[data]] in [[list]]. Returns
3173 [[-1]] if [[data]] is not in [[list]].
3174 <<list index declaration>>=
3175 int list_index(list_t *list, void *data);
3179 int list_index(list_t *list, void *data)
3183 while (list != NULL) {
3184 if (list->d == data) return i;
3193 [[list_element]] returns the data bound to the $i^\text{th}$ node of the list.
3194 <<list element declaration>>=
3195 void *list_element(list_t *list, int i);
3198 void *list_element(list_t *list, int i)
3202 while (list != NULL) {
3203 if (i == j) return list->d;
3213 [[unique]] removes duplicates from a sorted list according to a user
3214 specified comparison. Don't do this unless you have other pointers
3215 any dynamically allocated data in your list, because [[unique]] just
3216 drops any non-unique data without freeing it.
3217 <<unique declaration>>=
3218 void unique(list_t **pList, comparison_fn_t *comp);
3221 void unique(list_t **pList, comparison_fn_t *comp)
3224 assert(pList != NULL);
3226 assert(list != NULL); /* not an empty list */
3228 while (list->next != NULL) {
3229 if ((*comp)(list->d, list->next->d) == 0) {
3238 [[list_print]] prints a list to a [[FILE *]] stream.
3239 <<list print declaration>>=
3240 void list_print(FILE *file, list_t *list, char *list_name);
3243 void list_print(FILE *file, list_t *list, char *list_name)
3245 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
3247 while (list != NULL) {
3248 fprintf(file, " %p", list->d);
3251 fprintf(file, "\n");
3255 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
3256 <<create and destroy node>>=
3257 list_t *create_node(void *data)
3259 list_t *ret = (list_t *)malloc(sizeof(list_t));
3260 assert(ret != NULL);
3267 void destroy_node(list_t *node)
3273 The user must free the data pointed to by the node on their own.
3275 \subsection{List implementation}
3285 #include <assert.h> /* assert() */
3286 #include <stdlib.h> /* malloc(), free(), rand() */
3287 #include <stdio.h> /* fprintf(), stdout, FILE */
3288 #include "global.h" /* destroy_data_func_t */
3291 \subsection{List unit tests}
3293 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3296 <<list check includes>>
3298 <<main check program>>
3301 <<list check includes>>=
3302 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3303 #include <stdio.h> /* FILE */
3309 <<list test suite>>=
3312 <<list suite definition>>
3315 <<list suite definition>>=
3316 Suite *test_suite (void)
3318 Suite *s = suite_create ("list");
3319 <<push test case defs>>
3320 <<pop test case defs>>
3322 <<push test case adds>>
3323 <<pop test case adds>>
3329 START_TEST(test_push)
3334 push(&list, (void *)i, 1);
3335 fail_unless(list == head(list), NULL);
3336 fail_unless( (int)list->d == 0 );
3337 for (i=0; i<n; i++) {
3338 /* we expect 0, n-1, n-2, ..., 1 */
3341 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
3347 <<push test case defs>>=
3348 TCase *tc_push = tcase_create("push");
3351 <<push test case adds>>=
3352 tcase_add_test(tc_push, test_push);
3353 suite_add_tcase(s, tc_push);
3358 <<pop test case defs>>=
3360 <<pop test case adds>>=
3363 \section{Function string evaluation}
3365 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).
3366 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
3367 The solution is to run a scripting language as a subprocess accessed via pipes.
3368 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
3370 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
3371 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
3372 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.
3373 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
3375 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.
3376 Then you can either statically or dynamically link to those hardcoded functions.
3377 While much less flexible, this approach would be a fairly simple-to-implement fallback.
3379 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
3380 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
3383 #ifndef STRING_EVAL_H
3384 #define STRING_EVAL_H
3386 <<string eval setup declaration>>
3387 <<string eval function declaration>>
3388 <<string eval teardown declaration>>
3389 #endif /* STRING_EVAL_H */
3392 <<string eval module makefile lines>>=
3393 NW_SAWSIM_MODS += string_eval
3394 CHECK_BINS += check_string_eval
3395 check_string_eval_MODS =
3398 For an introduction to POSIX process control, see\\
3399 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
3400 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
3401 and of course, the relavent [[man]] pages.
3403 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
3404 [[execvp]] replaces the calling process' program with a new program.
3405 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
3406 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
3407 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
3408 The new program needs command line arguments, just like it would if you were running it from a shell.
3409 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
3410 with the final array entry being a [[NULL]] pointer.
3412 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
3413 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
3414 <<string eval subprocess definitions>>=
3415 #define SUBPROCESS "python"
3416 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
3417 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
3418 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
3420 The [[i]] option lets Python know that it should run in interactive mode.
3421 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
3422 In interactive mode, python acts on every instruction as soon as it is recieved.
3423 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.
3424 %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.
3426 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
3427 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
3428 [[fork]] returns two copies of the same program, executing the original code.
3429 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
3430 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
3432 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.
3433 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
3434 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
3435 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
3436 <<string eval pipe definitions>>=
3437 #define PIPE_READ 0 /* the end you read from */
3438 #define PIPE_WRITE 1 /* the end you write to */
3440 #define STDIN 0 /* initial index of stdin pair */
3441 #define STDOUT 2 /* initial index of stdout pair */
3444 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
3446 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.
3448 <<string eval setup declaration>>=
3449 extern void string_eval_setup(FILE **pIn, FILE **pOut);
3451 <<string eval setup definition>>=
3452 void string_eval_setup(FILE **pIn, FILE **pOut)
3455 int pfd[4]; /* file descriptors for stdin and stdout */
3457 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
3458 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
3460 pid = fork(); /* split process into two copies */
3462 if (pid == -1) { /* fork error */
3463 perror("fork error");
3465 } else if (pid == 0) { /* child */
3466 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
3467 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
3468 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
3469 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
3470 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
3471 perror("exec error"); /* exec shouldn't return */
3473 } else { /* parent */
3474 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
3475 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
3476 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
3477 if ( *pIn == NULL ) {
3478 perror("fdopen (in)");
3481 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
3482 if ( *pOut == NULL ) {
3483 perror("fdopen (out)");
3490 To use the evaluating subprocess, we just pipe in our command, and read out the results.
3491 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
3492 <<string eval function declaration>>=
3493 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
3495 <<string eval function definition>>=
3496 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
3499 rc = fprintf(in, "%s", input);
3500 assert(rc == strlen(input));
3503 alarm(1); /* set a one second timeout on the read */
3504 assert( fgets(output, buflen, out) != NULL );
3505 alarm(0); /* cancel the timeout */
3506 //fprintf(stderr, "eval: %s ----> %s", input, output);
3509 The [[alarm]] calls set and clear a timeout on the returned output.
3510 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.
3511 This protects against invalid input for which a line of output is not printed to [[stdout]].
3512 Other invalid inputs (e.g.\ those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
3513 If you are getting strange results, check your python code seperately. TODO, better error handling.
3515 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
3516 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
3517 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
3518 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
3519 <<string eval teardown declaration>>=
3520 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
3523 <<string eval teardown definition>>=
3524 void string_eval_teardown(FILE **pIn, FILE **pOut)
3529 /* redirect Python's stderr */
3530 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
3534 assert( fclose(*pIn) == 0 );
3536 assert( fclose(*pOut) == 0 );
3539 /* wait for python to exit */
3541 pid = wait(&stat_loc);
3548 if (WIFEXITED(stat_loc)) {
3549 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
3550 } else if (WIFSIGNALED(stat_loc)) {
3551 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
3556 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
3558 \subsection{String evaluation implementation}
3562 <<string eval includes>>
3563 #include "string_eval.h"
3564 <<string eval internal definitions>>
3565 <<string eval functions>>
3568 <<string eval includes>>=
3569 #include <assert.h> /* assert() */
3570 #include <stdlib.h> /* NULL */
3571 #include <stdio.h> /* fprintf(), stdout, fdopen() */
3572 #include <string.h> /* strlen() */
3573 #include <sys/types.h> /* pid_t */
3574 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
3575 #include <wait.h> /* wait() */
3578 <<string eval internal definitions>>=
3579 <<string eval subprocess definitions>>
3580 <<string eval pipe definitions>>
3583 <<string eval functions>>=
3584 <<string eval setup definition>>
3585 <<string eval function definition>>
3586 <<string eval teardown definition>>
3589 \subsection{String evaluation unit tests}
3591 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3592 <<check-string-eval.c>>=
3594 <<string eval check includes>>
3595 <<string comparison definition and globals>>
3596 <<string eval test suite>>
3597 <<main check program>>
3600 <<string eval check includes>>=
3601 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3602 #include <stdio.h> /* printf() */
3603 #include <assert.h> /* assert() */
3604 #include <string.h> /* strlen() */
3605 #include <signal.h> /* SIGKILL */
3607 #include "string_eval.h"
3610 <<string eval test suite>>=
3611 <<string eval tests>>
3612 <<string eval suite definition>>
3615 <<string eval suite definition>>=
3616 Suite *test_suite (void)
3618 Suite *s = suite_create ("string eval");
3619 <<string eval test case defs>>
3621 <<string eval test case adds>>
3626 <<string eval tests>>=
3627 START_TEST(test_setup_teardown)
3630 string_eval_setup(&in, &out);
3631 string_eval_teardown(&in, &out);
3634 START_TEST(test_invalid_command)
3637 char input[80], output[80]={};
3638 string_eval_setup(&in, &out);
3639 sprintf(input, "print undefined_name__traceback_should_be_printed_to_stderr\n");
3640 string_eval(in, out, input, 80, output);
3641 string_eval_teardown(&in, &out);
3644 START_TEST(test_simple_eval)
3647 char input[80], output[80]={};
3648 string_eval_setup(&in, &out);
3649 sprintf(input, "print 3+4*5\n");
3650 string_eval(in, out, input, 80, output);
3651 fail_unless(STRMATCH(output,"23\n"), NULL);
3652 string_eval_teardown(&in, &out);
3655 START_TEST(test_multiple_evals)
3658 char input[80], output[80]={};
3659 string_eval_setup(&in, &out);
3660 sprintf(input, "print 3+4*5\n");
3661 string_eval(in, out, input, 80, output);
3662 fail_unless(STRMATCH(output,"23\n"), NULL);
3663 sprintf(input, "print (3**2 + 4**2)**0.5\n");
3664 string_eval(in, out, input, 80, output);
3665 fail_unless(STRMATCH(output,"5.0\n"), NULL);
3666 string_eval_teardown(&in, &out);
3669 START_TEST(test_eval_with_state)
3672 char input[80], output[80]={};
3673 string_eval_setup(&in, &out);
3674 sprintf(input, "print 3+4*5\n");
3675 fprintf(in, "A = 3\n");
3676 sprintf(input, "print A*3\n");
3677 string_eval(in, out, input, 80, output);
3678 fail_unless(STRMATCH(output,"9\n"), NULL);
3679 string_eval_teardown(&in, &out);
3683 <<string eval test case defs>>=
3684 TCase *tc_string_eval = tcase_create("string_eval");
3686 <<string eval test case adds>>=
3687 tcase_add_test(tc_string_eval, test_setup_teardown);
3688 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
3689 tcase_add_test(tc_string_eval, test_simple_eval);
3690 tcase_add_test(tc_string_eval, test_multiple_evals);
3691 tcase_add_test(tc_string_eval, test_eval_with_state);
3692 suite_add_tcase(s, tc_string_eval);
3696 \section{Accelerating function evaluation}
3698 My first version-0.3 code was running very slowly.
3699 With the options suggested in the help
3700 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
3701 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
3702 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
3703 That's only 15 calls per solution, so the search algorithm seems reasonable.
3704 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
3710 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
3712 #endif /* ACCEL_K_H */
3715 <<accel k module makefile lines>>=
3716 NW_SAWSIM_MODS += accel_k
3717 CHECK_BINS += check_accel_k
3718 check_accel_k_MODS = interp k_model parse string_eval tavl
3723 #include <assert.h> /* assert() */
3724 #include <stdlib.h> /* realloc(), free(), NULL */
3725 #include "global.h" /* environment_t */
3726 #include "k_model.h" /* k_func_t */
3727 #include "interp.h" /* interp_* */
3728 #include "accel_k.h"
3730 typedef struct accel_k_struct {
3731 interp_table_t *itable;
3737 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3738 static int num_accels = 0;
3739 static accel_k_t *accels=NULL;
3741 /* Wrap k in the standard f(x) acceleration form */
3742 static double k_wrap(double F, void *params)
3744 accel_k_t *p = (accel_k_t *)params;
3745 return (*p->k)(F, p->env, p->params);
3748 static int k_tol(double FA, double kA, double FB, double kB)
3751 if (FB-FA > 1e-12) {
3752 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3753 return 1; /* unnacceptable */
3755 //printf("acceptable tol\n");
3756 return 0; /* acceptable */
3760 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3763 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3764 assert(accels != NULL);
3765 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3767 accels[i].env = env;
3768 accels[i].params = params;
3775 for (i=0; i<num_accels; i++)
3776 interp_table_free(accels[i].itable);
3778 if (accels) free(accels);
3782 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3785 for (i=0; i<num_accels; i++) {
3786 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3787 /* We've already setup this function.
3788 * WARNING: we're assuming param and env are the same. */
3789 return interp_table_eval(accels[i].itable, accels+i, F);
3792 /* set up a new acceleration environment */
3793 i = add_accel_k(k, env, params);
3794 return interp_table_eval(accels[i].itable, accels+i, F);
3798 \subsection{Unit tests}
3800 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3801 <<check-accel-k.c>>=
3803 <<accel k check includes>>
3804 <<check relative error>>
3806 <<accel k test suite>>
3807 <<main check program>>
3810 <<accel k check includes>>=
3811 #include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
3812 #include <stdio.h> /* sprintf() */
3813 #include <assert.h> /* assert() */
3814 #include <math.h> /* exp() */
3815 #include <errno.h> /* errno, ERANGE */
3818 #include "k_model.h"
3819 #include "accel_k.h"
3823 <<accel k test suite>>=
3824 <<accel k suite tests>>
3825 <<accel k suite definition>>
3828 <<accel k suite definition>>=
3829 Suite *test_suite (void)
3831 Suite *s = suite_create ("accelerated k model");
3832 <<accel k test case defs>>
3834 <<accel k test case adds>>
3839 <<accel k test case defs>>=
3840 TCase *tc_accel_k = tcase_create("accelerated k model");
3843 <<accel k test case adds>>=
3844 tcase_add_test(tc_accel_k, test_accel_k_accuracy);
3845 suite_add_tcase(s, tc_accel_k);
3848 <<accel k suite tests>>=
3849 START_TEST(test_accel_k_accuracy)
3851 double k, k_accel, F, Fm=0.0, FM=300e-12, dF=0.5e-12;
3852 k_func_t *k_func = bell_k;
3853 char *param_strings[] = {"3.3e-4", "0.25e-9"};
3854 void *params = string_create_bell_param_t(param_strings);
3855 environment_t env = {};
3858 for(F=Fm; F <= FM; F+=dF) {
3859 k = k_func(F, &env, params);
3860 k_accel = accel_k(k_func, F, &env, params);
3861 CHECK_ERR(1e-8, k, k_accel);
3863 destroy_bell_param_t(params);
3869 \section{Tension models}
3870 \label{sec.tension_models}
3872 TODO: tension model intro.
3873 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.
3875 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3876 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]].
3878 <<tension-model.h>>=
3879 #ifndef TENSION_MODEL_H
3880 #define TENSION_MODEL_H
3882 <<tension handler types>>
3883 <<constant tension model declarations>>
3884 <<hooke tension model declarations>>
3885 <<worm-like chain tension model declarations>>
3886 <<freely-jointed chain tension model declarations>>
3887 <<piston tension model declarations>>
3888 <<find tension definitions>>
3889 <<tension model global declarations>>
3890 #endif /* TENSION_MODEL_H */
3893 <<tension model module makefile lines>>=
3894 NW_SAWSIM_MODS += tension_model
3895 CHECK_BINS += check_tension_model
3896 check_tension_model_MODS = list tension_balance
3898 <<tension model utils makefile lines>>=
3899 TENSION_MODEL_MODS = tension_model parse list tension_balance
3900 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3901 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3902 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3903 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
3904 notangle -Rtension-model-utils.c $< > $@
3905 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
3906 $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3907 clean_tension_model_utils :
3908 rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3909 @ The pipe symbol [[|]] marks the prerequisits following it (in this
3910 case, the directories) as ‘order-only’ prerequisites. The timestamp
3911 on these prerequisits does not effect whether the rules are executed.
3912 This is appropriate for directories, where we don't need to recompile
3913 after an unrelated has been added to the directory, but only when the
3914 source prerequisites change. See the [[make]] documentation for more
3916 (\url{http://www.gnu.org/software/make/manual/html_node/Prerequisite-Types.html}).
3919 \label{sec.null_tension}
3921 For unstretchable domains.
3923 <<null tension model getopt>>=
3924 {NULL, NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3927 \subsection{Constant}
3928 \label{sec.const_tension}
3930 An infinitely stretchable domain providing a constant tension.
3931 Obviously this cannot be inverted, so you can't put this domain in
3932 series with any other domains. We include it mostly for testing
3935 <<constant tension functions>>=
3936 <<constant tension function>>
3937 <<constant tension parameter create/destroy functions>>
3940 <<constant tension model declarations>>=
3941 <<constant tension function declaration>>
3942 <<constant tension parameter create/destroy function declarations>>
3943 <<constant tension model global declarations>>
3947 <<constant tension function declaration>>=
3948 extern double const_tension_handler(double x, void *pdata);
3950 <<constant tension function>>=
3951 double const_tension_handler(double x, void *pdata)
3953 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3954 list_t *list = data->group_tension_params;
3959 assert(list != NULL); /* empty active group?! */
3960 F = ((const_tension_param_t *)list->d)->F;
3962 fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3964 while (list != NULL) {
3965 assert(((const_tension_param_t *)list->d)->F == F);
3973 <<constant tension structure>>=
3974 typedef struct const_tension_param_struct {
3975 double F; /* tension (force) in N */
3976 } const_tension_param_t;
3980 <<constant tension parameter create/destroy function declarations>>=
3981 extern void *string_create_const_tension_param_t(char **param_strings);
3982 extern void destroy_const_tension_param_t(void *p);
3985 <<constant tension parameter create/destroy functions>>=
3986 const_tension_param_t *create_const_tension_param_t(double F)
3988 const_tension_param_t *ret
3989 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3990 assert(ret != NULL);
3995 void *string_create_const_tension_param_t(char **param_strings)
3997 return create_const_tension_param_t(
3998 safe_strtod(param_strings[0],__FUNCTION__));
4001 void destroy_const_tension_param_t(void *p)
4009 <<constant tension model global declarations>>=
4010 extern int num_const_tension_params;
4011 extern char *const_tension_param_descriptions[];
4012 extern char const_tension_param_string[];
4015 <<constant tension model globals>>=
4016 int num_const_tension_params = 1;
4017 char *const_tension_param_descriptions[] = {"tension F, N"};
4018 char const_tension_param_string[] = "0";
4021 <<constant tension model getopt>>=
4022 {&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}
4028 <<hooke functions>>=
4029 <<internal hooke functions>>
4031 <<inverse hooke handler>>
4032 <<hooke parameter create/destroy functions>>
4035 <<hooke tension model declarations>>=
4036 <<hooke tension function declaration>>
4037 <<hooke tension parameter create/destroy function declarations>>
4038 <<hooke tension model global declarations>>
4042 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
4043 The behavior of a series of springs $k_i$ in series is given by
4045 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
4047 For a simple proof, see Appendix \ref{sec.math_hooke}.
4049 <<hooke tension function declaration>>=
4050 extern double hooke_handler(double x, void *pdata);
4051 extern double inverse_hooke_handler(double F, void *pdata);
4054 First we define a function that computes the effective spring constant
4055 (stored in a single [[hooke_param_t]]) for the entire group.
4056 <<internal hooke functions>>=
4057 static void hooke_param_grouper(tension_handler_data_t *data,
4058 hooke_param_t *param) {
4059 list_t *list = data->group_tension_params;
4063 assert(list != NULL); /* empty active group?! */
4064 while (list != NULL) {
4065 assert( ((hooke_param_t *)list->d)->k > 0 );
4066 k += 1.0/ ((hooke_param_t *)list->d)->k;
4068 fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
4069 __FUNCTION__, ((hooke_param_t *)list->d)->k);
4075 fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N (from data %p)\n",
4076 __FUNCTION__, k, x, k*x, data);
4083 Everyone knows Hooke's law: $F=kx$.
4085 double hooke_handler(double x, void *pdata)
4087 hooke_param_t param = {0};
4090 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
4096 Inverting Hooke's law gives $x=F/k$.
4097 <<inverse hooke handler>>=
4098 double inverse_hooke_handler(double F, void *pdata)
4100 hooke_param_t param = {0};
4103 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
4109 Not much to keep track of for a single effective spring, just the
4110 spring constant $k$.
4111 <<hooke structure>>=
4112 typedef struct hooke_param_struct {
4113 double k; /* spring constant in N/m */
4118 We wrap up our Hooke implementation with some book-keeping functions.
4119 <<hooke tension parameter create/destroy function declarations>>=
4120 extern void *string_create_hooke_param_t(char **param_strings);
4121 extern void destroy_hooke_param_t(void *p);
4125 <<hooke parameter create/destroy functions>>=
4126 hooke_param_t *create_hooke_param_t(double k)
4128 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
4129 assert(ret != NULL);
4134 void *string_create_hooke_param_t(char **param_strings)
4136 return create_hooke_param_t(safe_strtod(param_strings[0], __FUNCTION__));
4139 void destroy_hooke_param_t(void *p)
4146 <<hooke tension model global declarations>>=
4147 extern int num_hooke_params;
4148 extern char *hooke_param_descriptions[];
4149 extern char hooke_param_string[];
4152 <<hooke tension model globals>>=
4153 int num_hooke_params = 1;
4154 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
4155 char hooke_param_string[]="0.05";
4158 <<hooke tension model getopt>>=
4159 {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}
4162 \subsection{Worm-like chain}
4165 We can model several unfolded domains with a single worm-like chain.
4166 <<worm-like chain functions>>=
4167 <<internal worm-like chain functions>>
4168 <<worm-like chain handler>>
4169 <<inverse worm-like chain handler>>
4170 <<worm-like chain create/destroy functions>>
4173 <<worm-like chain tension model declarations>>=
4175 <<worm-like chain handler declaration>>
4176 <<inverse worm-like chain handler declaration>>
4177 <<worm-like chain create/destroy function declarations>>
4178 <<worm-like chain tension model global declarations>>
4182 <<internal worm-like chain functions>>=
4183 <<worm-like chain function>>
4184 <<inverse worm-like chain function>>
4185 <<worm-like chain parameter grouper>>
4188 The combination of all unfolded domains is modeled as a worm like chain,
4189 with the force $F_{\text{WLC}}$ approximately given by
4191 F_{\text{WLC}} \approx \frac{k_B T}{p}
4192 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
4194 where $x$ is the distance between the fixed ends,
4195 $k_B$ is Boltzmann's constant,
4196 $T$ is the absolute temperature,
4197 $p$ is the persistence length, and
4198 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
4199 <<worm-like chain function>>=
4200 static double wlc(double x, double T, double p, double L)
4202 double xL=0.0; /* uses global k_B */
4204 assert(T > 0); assert(p > 0); assert(L > 0);
4205 if (x >= L) return HUGE_VAL;
4206 xL = x/L; /* unitless */
4207 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
4208 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
4209 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
4214 Inverting the approximate WLC is fairly straightforward, if a bit tedious.
4215 Letting $x_L \equiv x/L$ and $F_T \equiv Fp/k_B T$ we have
4217 F_T &\approx \frac{1}{4} \p({\frac{1}{(1-x_L)^2} - 1}) + x_L \\
4218 F_T (1-x_L)^2 &\approx \frac{1}{4} - \frac{1}{4} (1-x_L)^2 + x_L(1-x_L)^2 \\
4219 0 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) (1-x_L)^2 + x_L(1-x_L)^2 \\
4220 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) \p({1-2x_L+x_L^2})
4221 + x_L - 2x_L^2 + x_L^3 \\
4222 &\approx \p[{\frac{1}{4} - \p({F_T + \frac{1}{4}})}]
4223 + x_L \p[{2F_T + \frac{1}{2} + 1}]
4224 + x_L^2 \p[{-\p({F_T + \frac{1}{4}}) - 2}]
4227 + x_L \p({2F_T + \frac{3}{2}})
4228 - x_L^2 \p({F_T + \frac{9}{4}})
4230 &\approx ax_L^3 + bx_L^2 + cx_L + d
4232 where $a \equiv 1$, $b \equiv -F_T - 9/4$, $c \equiv 2F_T + 3/2$, and
4234 % From \citet{wikipedia_cubic_function} the discriminant
4236 % \Delta &= -4b^3d + b^2c^2 - 4ac^3 + 18abcd - 27a^2d^2 \\
4237 % &= -4F_T\p({F_T + \frac{9}{4}})^3
4238 % + \p({F_T + \frac{9}{4}})^2 \cdot \p({2F_T + \frac{3}{2}})^2
4239 % - 4 \p({2F_T + \frac{3}{2}})^3
4240 % + 18F_T \p({F_T + \frac{9}{4}}) \cdot \p({2F_T + \frac{3}{2}})
4242 % &= -4F_T \p({F_T^3 + \frac{27}{4}F_T^2 +\frac{243}{16}F_T +\frac{729}{64}})
4243 %% a^3 + 3a^2b + 3ab^2 + b^3
4244 % + \p({F_T^2 + \frac{9}{2}F_T + \frac{81}{16}})
4245 % \cdot \p({4F_T^2 + 6F_T + \frac{9}{4}}) \\
4246 % &\qquad - 4 \p({8F_T^3 + 18F_T^2 + \frac{27}{2}F_T + \frac{27}{8}})
4247 %% a^3 + 3a^2b + 3ab^2 + b^3
4248 % + 18F_T \p({2F_T^2 + 6F_T + \frac{27}{8}})
4250 % &= -4F_T^4 - 27F_T^3 - \frac{243}{4}F_T^2 - \frac{729}{16}F_T
4251 % + 4F_T^4 + 24F_T^3 +\frac{99}{2}F_T^2 +\frac{81}{2}F_T +\frac{729}{64}\\
4252 % &\qquad - 32F_T^3 - 72F_T^2 - 54F_T - \frac{27}{2}
4253 % + 36F_T^3 + 108F_T^2 + \frac{243}{4}F_T
4255 % &= F_T^4 \p({-4 + 4}) + F_T^3 \p({-27 + 24 - 32 + 36})
4256 % + F_T^2 \p({-\frac{243}{4} + \frac{99}{2} - 72 + 108 -27}) \\
4257 % &\qquad + F_T \p({-\frac{729}{16} + \frac{81}{2} - 54 + \frac{243}{4}})
4258 % \p({\frac{729}{64} - \frac{27}{2}}) \\
4259 % &= F_T^3 - \frac{9}{4}F_T^2 + \frac{27}{16}F_T - \frac{135}{64}
4261 We can plug these coefficients into the GSL cubic solver to invert the
4262 WLC\citep{galassi05}. The applicable root is always the one which
4263 comes before the singularity, which will be the smallest real root.
4264 <<inverse worm-like chain function>>=
4265 static double inverse_wlc(double F, double T, double p, double L)
4267 double FT = F*p/(k_B*T); /* uses global k_B */
4268 double xL0, xL1, xL2;
4271 assert(T > 0); assert(p > 0); assert(L > 0);
4274 num_roots = gsl_poly_solve_cubic(-(FT+2.25),2*FT+1.5,-FT, &xL0, &xL1, &xL2);
4275 assert(num_roots > 0); assert(xL0 >= -DOUBLE_PRECISION); assert(xL0 < 1);
4276 if (xL0 < 0) xL0 = 0.0;
4282 First we define a function that computes the effective WLC parameters
4283 (stored in a single [[wlc_param_t]]) for the entire group.
4284 <<worm-like chain parameter grouper>>=
4285 static void wlc_param_grouper(tension_handler_data_t *data,
4286 wlc_param_t *param) {
4287 list_t *list = data->group_tension_params;
4288 double p=0.0, L=0.0;
4291 assert(list != NULL); /* empty active group?! */
4292 p = ((wlc_param_t *)list->d)->p;
4293 while (list != NULL) {
4294 /* enforce consistent p */
4295 assert( ((wlc_param_t *)list->d)->p == p);
4296 L += ((wlc_param_t *)list->d)->L;
4298 fprintf(stderr, "%s: WLC section %g m long in series\n",
4299 __FUNCTION__, ((wlc_param_t *)list->d)->L);
4304 fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
4305 __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
4313 <<worm-like chain handler declaration>>=
4314 extern double wlc_handler(double x, void *pdata);
4317 This model requires that each unfolded domain in the group have the
4318 same persistence length.
4319 <<worm-like chain handler>>=
4320 double wlc_handler(double x, void *pdata)
4322 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4323 wlc_param_t param = {0};
4326 wlc_param_grouper(data, ¶m);
4327 return wlc(x, data->env->T, param.p, param.L);
4332 <<inverse worm-like chain handler declaration>>=
4333 extern double inverse_wlc_handler(double F, void *pdata);
4336 <<inverse worm-like chain handler>>=
4337 double inverse_wlc_handler(double F, void *pdata)
4339 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4340 wlc_param_t param = {0};
4342 wlc_param_grouper(data, ¶m);
4343 return inverse_wlc(F, data->env->T, param.p, param.L);
4348 <<worm-like chain structure>>=
4349 typedef struct wlc_param_struct {
4350 double p; /* WLC persistence length (m) */
4351 double L; /* WLC contour length (m) */
4355 <<worm-like chain create/destroy function declarations>>=
4356 extern void *string_create_wlc_param_t(char **param_strings);
4357 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
4360 <<worm-like chain create/destroy functions>>=
4361 wlc_param_t *create_wlc_param_t(double p, double L)
4363 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
4364 assert(ret != NULL);
4370 void *string_create_wlc_param_t(char **param_strings)
4372 return create_wlc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4373 safe_strtod(param_strings[1], __FUNCTION__));
4376 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
4384 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.
4385 TODO (now we copy the string before parsing, but still do this for clarity).
4387 <<valid string write code>>=
4388 char string[] = "abc";
4391 <<valid string write code>>=
4392 char *string = "abc";
4396 <<worm-like chain tension model global declarations>>=
4397 extern int num_wlc_params;
4398 extern char *wlc_param_descriptions[];
4399 extern char wlc_param_string[];
4401 <<worm-like chain tension model globals>>=
4402 int num_wlc_params = 2;
4403 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
4404 char wlc_param_string[]="0.39e-9,28.4e-9";
4407 <<worm-like chain tension model getopt>>=
4408 {&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}
4410 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
4412 \subsection{Freely-jointed chain}
4415 <<freely-jointed chain functions>>=
4416 <<freely-jointed chain function>>
4417 <<freely-jointed chain handler>>
4418 <<freely-jointed chain create/destroy functions>>
4421 <<freely-jointed chain tension model declarations>>=
4422 <<freely-jointed chain handler declaration>>
4423 <<freely-jointed chain create/destroy function declarations>>
4424 <<freely-jointed chain tension model global declarations>>
4428 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
4429 The tension of a chain of $N$ such links is given by
4431 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
4433 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}.
4434 <<freely-jointed chain function>>=
4435 static double langevin(double x, void *params)
4438 return 1.0/tanh(x) - 1/x;
4441 static double inv_langevin(double x)
4443 double lb=0.0, ub=1.0, ret;
4444 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
4445 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
4449 static double fjc(double x, double T, double l, int N)
4451 double L = l*(double)N;
4452 if (x == 0) return 0;
4453 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
4454 return k_B*T/l * inv_langevin(x/L);
4459 <<freely-jointed chain handler declaration>>=
4460 extern double fjc_handler(double x, void *pdata);
4462 <<freely-jointed chain handler>>=
4463 double fjc_handler(double x, void *pdata)
4465 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4466 list_t *list = data->group_tension_params;
4471 assert(list != NULL); /* empty active group?! */
4472 l = ((fjc_param_t *)list->d)->l;
4473 while (list != NULL) {
4474 /* enforce consistent l */
4475 assert( ((fjc_param_t *)list->d)->l == l);
4476 N += ((fjc_param_t *)list->d)->N;
4478 fprintf(stderr, "%s: FJC section %d links long in series\n",
4479 __FUNCTION__, ((fjc_param_t *)list->d)->N);
4484 fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
4485 __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
4487 return fjc(x, data->env->T, l, N);
4491 <<freely-jointed chain structure>>=
4492 typedef struct fjc_param_struct {
4493 double l; /* FJC link length (m) */
4494 int N; /* FJC number of links */
4498 <<freely-jointed chain create/destroy function declarations>>=
4499 extern void *string_create_fjc_param_t(char **param_strings);
4500 extern void destroy_fjc_param_t(void *p);
4503 <<freely-jointed chain create/destroy functions>>=
4504 fjc_param_t *create_fjc_param_t(double l, double N)
4506 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
4507 assert(ret != NULL);
4513 void *string_create_fjc_param_t(char **param_strings)
4515 return create_fjc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4516 safe_strtod(param_strings[1], __FUNCTION__));
4519 void destroy_fjc_param_t(void *p)
4526 <<freely-jointed chain tension model global declarations>>=
4527 extern int num_fjc_params;
4528 extern char *fjc_param_descriptions[];
4529 extern char fjc_param_string[];
4532 <<freely-jointed chain tension model globals>>=
4533 int num_fjc_params = 2;
4534 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
4535 char fjc_param_string[]="0.5e-9,200";
4538 <<freely-jointed chain tension model getopt>>=
4539 {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}
4543 \label{sec.piston_tension}
4545 An infinitely stretchable domain with no tension for extensions less
4546 than a particular threshold $L$, and infinite tension for $x>L$. The
4547 tension at $x=L$ is undefined, but may be any positive value. The
4548 model is called the ``piston'' model because it resembles a
4549 frictionless piston sliding in a rigid cylinder of length $L$.
4552 0 & \text{if $x<L$}, \\
4553 \infty & \text{if $x>L$}.
4557 Note that because of it's infinte stiffness at $x=L$, fully extended
4558 domains with this tension model will be identical to completely rigid
4559 domains. However, there is a distinction between this model and rigid
4560 domains for $x<L$, because any reactions that occur at $F=0$
4561 (e.g. refolding) will have more time to occur.
4563 Because the tension at $x=L$ is undefined, the user must make sure
4564 domains with this tension model are never the initial active group,
4565 because it would break [[tension_balance()]] in [[find_tension()]]
4566 (see Section \ref{sec.tension_balance}).
4568 <<piston tension functions>>=
4569 <<piston tension parameter grouper>>
4570 <<piston tension handler>>
4571 <<inverse piston tension handler>>
4572 <<piston tension parameter create/destroy functions>>
4575 <<piston tension model declarations>>=
4576 <<piston tension handler declaration>>
4577 <<inverse piston tension handler declaration>>
4578 <<piston tension parameter create/destroy function declarations>>
4579 <<piston tension model global declarations>>
4583 First we define a function that computes the effective piston parameters
4584 (stored in a single [[piston_tension_param_t]]) for the entire group.
4585 <<piston tension parameter grouper>>=
4586 static void piston_tension_param_grouper(tension_handler_data_t *data,
4587 piston_tension_param_t *param) {
4588 list_t *list = data->group_tension_params;
4592 assert(list != NULL); /* empty active group?! */
4593 while (list != NULL) {
4594 L += ((piston_tension_param_t *)list->d)->L;
4600 <<piston tension handler declaration>>=
4601 extern double piston_tension_handler(double x, void *pdata);
4603 <<piston tension handler>>=
4604 double piston_tension_handler(double x, void *pdata)
4606 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4607 piston_tension_param_t param = {0};
4610 piston_tension_param_grouper(data, ¶m);
4611 if (x <= param.L) F = 0;
4614 fprintf(stderr, "%s: x %g, L %g, tension %g N\n", __FUNCTION__, x, param.L, F);
4621 We abuse [[x_of_xo()]] and [[oneD_solve()]] with this inverse
4622 definition (see Section \ref{sec.tension_balance}), because the
4623 $x(F=0)<=L$, but our function returns $x(F=0)=0$. This is fine when
4624 the total extension \emph{is} zero, but cheats a bit when there is
4625 some total extension. For any non-zero extension, the initial active
4626 group produces some tension (we hope. True for all our other tension
4627 models). This tension fully extends the piston group (of which there
4628 should be only one, since all piston states can get lumped into the
4629 same tension group.). If the total extension is $\le$ the target
4630 extension, the full extension is accurate, so the inverse was valid
4631 after all. If not, [[oneD_solve()]] will halve the extension of the
4632 first group, reducing the overall tension. For total extension less
4633 than the piston extension, this bisection continues for [[max_steps]],
4634 leaving $x_0$ reduced by a factor of $2^\text{max\_steps}$, a very
4635 small number. Because of this, the returned tension $F_0(x_0)$ is
4636 very small, even though the total extension found by [[x_of_xo()]]
4637 is still $>$ the piston length.
4639 While this works, it is clearly not the most elegant or robust
4640 solution. TODO: think of (and implement) a better procedure.
4642 <<inverse piston tension handler declaration>>=
4643 extern double inverse_piston_tension_handler(double x, void *pdata);
4645 <<inverse piston tension handler>>=
4646 double inverse_piston_tension_handler(double F, void *pdata)
4648 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4649 piston_tension_param_t param = {0};
4652 piston_tension_param_grouper(data, ¶m);
4653 if (F == 0.0) return 0;
4659 <<piston tension structure>>=
4660 typedef struct piston_tension_param_struct {
4661 double L; /* length of piston in m */
4662 } piston_tension_param_t;
4666 <<piston tension parameter create/destroy function declarations>>=
4667 extern void *string_create_piston_tension_param_t(char **param_strings);
4668 extern void destroy_piston_tension_param_t(void *p);
4672 <<piston tension parameter create/destroy functions>>=
4673 piston_tension_param_t *create_piston_tension_param_t(double L)
4675 piston_tension_param_t *ret
4676 = (piston_tension_param_t *) malloc(sizeof(piston_tension_param_t));
4677 assert(ret != NULL);
4682 void *string_create_piston_tension_param_t(char **param_strings)
4684 return create_piston_tension_param_t(
4685 safe_strtod(param_strings[0],__FUNCTION__));
4688 void destroy_piston_tension_param_t(void *p)
4696 <<piston tension model global declarations>>=
4697 extern int num_piston_tension_params;
4698 extern char *piston_tension_param_descriptions[];
4699 extern char piston_tension_param_string[];
4703 <<piston tension model globals>>=
4704 int num_piston_tension_params = 1;
4705 char *piston_tension_param_descriptions[] = {"length L, m"};
4706 char piston_tension_param_string[] = "0";
4710 <<piston tension model getopt>>=
4711 {&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}
4714 \subsection{Tension model implementation}
4716 <<tension-model.c>>=
4719 <<tension model includes>>
4720 #include "tension_model.h"
4721 <<tension model internal definitions>>
4722 <<tension model globals>>
4723 <<tension model functions>>
4726 <<tension model includes>>=
4727 #include <assert.h> /* assert() */
4728 #include <stdlib.h> /* NULL, strto*() */
4729 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4730 #include <stdio.h> /* fprintf(), stdout */
4731 #include <errno.h> /* errno, ERANGE */
4734 #include "tension_balance.h" /* oneD_*() */
4737 <<tension model internal definitions>>=
4738 <<constant tension structure>>
4740 <<worm-like chain structure>>
4741 <<freely-jointed chain structure>>
4742 <<piston tension structure>>
4745 <<tension model globals>>=
4746 <<tension handler array global>>
4747 <<constant tension model globals>>
4748 <<hooke tension model globals>>
4749 <<worm-like chain tension model globals>>
4750 <<freely-jointed chain tension model globals>>
4751 <<piston tension model globals>>
4754 <<tension model functions>>=
4756 <<constant tension functions>>
4758 <<worm-like chain functions>>
4759 <<freely-jointed chain functions>>
4760 <<piston tension functions>>
4763 \subsection{Utilities}
4765 The tension models can get complicated, and users may want to reassure
4766 themselves that this portion of the input physics are functioning properly. The
4767 stand-alone program [[tension_model_utils]] demonstrates the tension model
4768 interface without the overhead of having to understand a full simulation.
4769 [[tension_model_utils]] takes command line model arguments like the full
4770 simulation, and prints $F(x)$ data to the screen for the selected model over a
4773 <<tension-model-utils.c>>=
4775 <<tension model utility includes>>
4776 <<tension model utility definitions>>
4777 <<tension model utility getopt functions>>
4779 <<tension model utility main>>
4782 <<tension model utility main>>=
4783 int main(int argc, char **argv)
4785 <<tension model getopt array definition>>
4786 tension_model_getopt_t *model = NULL;
4790 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
4791 tension_handler_data_t tdata;
4792 int num_param_args; /* for INIT_MODEL() */
4793 char **param_args; /* for INIT_MODEL() */
4795 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model,
4796 &Fmax, &Xmax, &flags);
4798 if (flags & VFLAG) {
4799 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4801 INIT_MODEL("utility", model, model->params, params);
4802 tdata.group_tension_params = NULL;
4803 push(&tdata.group_tension_params, params, 1);
4805 tdata.persist = NULL;
4806 if (model->handler == NULL) {
4807 printf("No tension function for model %s\n", model->name);
4813 printf("#Distance (m)\tForce (N)\n");
4814 for (i=0; i<=N; i++) {
4815 x = Xmax*i/(double)N;
4816 F = (*model->handler)(x, &tdata);
4817 if (F < 0 || F > Fmax) break;
4818 printf("%g\t%g\n", x, F);
4820 if (flags & VFLAG && i <= N) {
4821 /* explain exit condition */
4823 printf("Impossible force %g\n", F);
4825 printf("Reached large force limit %g > %g\n", F, Fmax);
4828 params = pop(&tdata.group_tension_params);
4829 if (model->destructor)
4830 (*model->destructor)(params);
4835 <<tension model utility includes>>=
4838 #include <string.h> /* strlen() */
4839 #include <assert.h> /* assert() */
4840 #include <errno.h> /* errno, ERANGE */
4844 #include "tension_model.h"
4847 <<tension model utility definitions>>=
4848 <<version definition>>
4849 #define VFLAG 1 /* verbose */
4850 <<string comparison definition and globals>>
4851 <<tension model getopt definitions>>
4852 <<initialize model definition>>
4856 <<tension model utility getopt functions>>=
4858 <<index tension model>>
4859 <<tension model utility help>>
4860 <<tension model utility get options>>
4863 <<tension model utility help>>=
4864 void help(char *prog_name,
4866 int n_tension_models, tension_model_getopt_t *tension_models,
4867 int tension_model, double Fmax, double Xmax)
4870 printf("usage: %s [options]\n", prog_name);
4871 printf("Version %s\n\n", VERSION);
4872 printf("A utility for understanding the available tension models\n\n");
4873 printf("Environmental options:\n");
4874 printf("-T T\tTemperature T (currently %g K)\n", env->T);
4875 printf("-C T\tYou can also set the temperature T in Celsius\n");
4876 printf("Model options:\n");
4877 printf("-m model\tFolded tension model (currently %s)\n",
4878 tension_models[tension_model].name);
4879 printf("-a args\tFolded tension model argument string (currently %s)\n",
4880 tension_models[tension_model].params);
4881 printf("F(x) is calculated for a range of x and printed\n");
4882 printf("For example:\n");
4883 printf(" #Distance (m)\tForce (N)\n");
4884 printf(" 123.456\t7.89\n");
4886 printf("-F\tSet the maximum F value for the standard mode F(x) output (currently %g)\n", Fmax);
4887 printf("-X\tSet the maximum x value for the standard mode F(x) output (currently %g)\n", Xmax);
4888 printf("-V\tChange output to verbose mode\n");
4889 printf("-h\tPrint this help and exit\n");
4891 printf("Tension models:\n");
4892 for (i=0; i<n_tension_models; i++) {
4893 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
4894 for (j=0; j < tension_models[i].num_params ; j++)
4895 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
4896 printf("\t\tdefaults: %s\n", tension_models[i].params);
4901 <<tension model utility get options>>=
4902 void get_options(int argc, char **argv, environment_t *env,
4903 int n_tension_models, tension_model_getopt_t *tension_models,
4904 tension_model_getopt_t **model, double *Fmax, double *Xmax,
4905 unsigned int *flags)
4907 char *prog_name = NULL;
4908 char c, options[] = "T:C:m:a:F:X:Vh";
4909 int tension_model=0, num_strings;
4910 extern char *optarg;
4911 extern int optind, optopt, opterr;
4915 /* setup defaults */
4917 prog_name = argv[0];
4918 env->T = 300.0; /* K */
4919 *Fmax = 1e5; /* N */
4920 *Xmax = 1e-6; /* m */
4922 *model = tension_models;
4924 while ((c=getopt(argc, argv, options)) != -1) {
4926 case 'T': env->T = safe_strtod(optarg, "-T"); break;
4927 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
4929 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
4930 *model = tension_models+tension_model;
4933 tension_models[tension_model].params = optarg;
4935 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
4936 case 'X': *Xmax = safe_strtod(optarg, "-X"); break;
4937 case 'V': *flags |= VFLAG; break;
4939 fprintf(stderr, "unrecognized option '%c'\n", optopt);
4940 /* fall through to default case */
4942 help(prog_name, env, n_tension_models, tension_models, tension_model, *Fmax, *Xmax);
4950 \subsection{Tension model unit tests}
4952 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
4953 <<check-tension-model.c>>=
4955 <<tension model check includes>>
4956 <<check relative error>>
4958 <<tension model test suite>>
4959 <<main check program>>
4962 <<tension model check includes>>=
4963 #include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
4964 #include <stdio.h> /* sprintf() */
4965 #include <assert.h> /* assert() */
4966 #include <math.h> /* exp() */
4967 #include <errno.h> /* errno, ERANGE */
4971 #include "tension_balance.h" /* oneD_*() */
4972 #include "tension_model.h"
4975 <<tension model test suite>>=
4977 <<const tension tests>>
4979 <<worm-like chain tests>>
4980 <<freely-jointed chain tests>>
4982 <<tension model suite definition>>
4985 <<tension model suite definition>>=
4986 Suite *test_suite (void)
4988 Suite *s = suite_create ("tension model");
4989 <<const tension test case defs>>
4990 <<hooke test case defs>>
4991 <<worm-like chain test case defs>>
4992 <<freely-jointed chain test case defs>>
4993 <<piston test case defs>>
4995 <<const tension test case adds>>
4996 <<hooke test case adds>>
4997 <<worm-like chain test case adds>>
4998 <<freely-jointed chain test case adds>>
4999 <<piston test case adds>>
5004 \subsubsection{Constant}
5006 <<const tension test case defs>>=
5007 TCase *tc_const_tension = tcase_create("Constant tension");
5010 <<const tension test case adds>>=
5011 tcase_add_test(tc_const_tension, test_const_tension_create_destroy);
5012 suite_add_tcase(s, tc_const_tension);
5015 <<const tension tests>>=
5016 START_TEST(test_const_tension_create_destroy)
5018 char *tension[] = {"1","2.2", "3"};
5019 char *params[] = {tension[0]};
5022 for( i=0; i < sizeof(tension)/sizeof(char *); i++) {
5023 params[0] = tension[i];
5024 p = string_create_const_tension_param_t(params);
5025 destroy_const_tension_param_t(p);
5030 @ TODO: In order to test the constant tension handler itself, we'd
5031 have to construct a group.
5034 \subsubsection{Hooke}
5036 <<hooke test case defs>>=
5037 TCase *tc_hooke = tcase_create("Hooke");
5040 <<hooke test case adds>>=
5041 tcase_add_test(tc_hooke, test_hooke_create_destroy);
5042 suite_add_tcase(s, tc_hooke);
5047 START_TEST(test_hooke_create_destroy)
5049 char *k[] = {"1","2.2", "3"};
5050 char *params[] = {k[0]};
5053 for( i=0; i < sizeof(k)/sizeof(char *); i++) {
5055 p = string_create_hooke_param_t(params);
5056 destroy_hooke_param_t(p);
5061 @ TODO: In order to test the Hooke tension handler itself, we'd
5062 have to construct a group.
5065 \subsubsection{Worm-like chain}
5067 <<worm-like chain test case defs>>=
5068 TCase *tc_wlc = tcase_create("WLC");
5071 <<worm-like chain test case adds>>=
5072 tcase_add_test(tc_wlc, test_wlc_at_zero);
5073 tcase_add_test(tc_wlc, test_wlc_at_half);
5074 suite_add_tcase(s, tc_wlc);
5078 <<worm-like chain tests>>=
5079 <<worm-like chain function>>
5080 START_TEST(test_wlc_at_zero)
5082 double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
5083 fail_unless(abs(wlc(x, T, p, L)) < lim, \
5084 "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
5085 x, T, p, L, abs(wlc(x, T, p, L)), lim);
5089 START_TEST(test_wlc_at_half)
5091 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
5092 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
5093 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
5095 * linear term = x/L = 0.5
5096 * nonlinear + linear = 0.75 + 0.5 = 1.25
5097 * wlc = 10e21*1.25 = 12.5e21
5099 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
5100 "wlc(%g, %g, %g, %g) = %g != %g",
5101 x, T, p, L, wlc(x, T, p, L), 12.5e21);
5108 \subsubsection{Freely-jointed chain}
5110 <<freely-jointed chain test case defs>>=
5111 TCase *tc_fjc = tcase_create("FJC");
5114 <<freely-jointed chain test case adds>>=
5115 tcase_add_test(tc_fjc, test_fjc_at_zero);
5116 tcase_add_test(tc_fjc, test_fjc_at_half);
5117 suite_add_tcase(s, tc_fjc);
5121 <<freely-jointed chain tests>>=
5122 <<freely-jointed chain function>>
5123 START_TEST(test_fjc_at_zero)
5126 double T=1.0, l=1.0, p=0.1, x=0.0, lim=1e-30;
5127 fail_unless(abs(fjc(x, T, l, N)) < lim, \
5128 "abs(fjc(x=%g,T=%g,l=%g,N=%d)) = %g >= %g",
5129 x, T, l, N, abs(fjc(x, T, l, N)), lim);
5133 START_TEST(test_fjc_at_half)
5136 double T=1.0/k_B, l=1.0, x=5.0;
5137 /* prefactor = k_B T / l = k_B (1.0/k_B) / 1.0 = 1.0 J/nm = 1.0e21 pN
5138 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
5140 * linear term = x/L = 0.5
5141 * nonlinear + linear = 0.75 + 0.5 = 1.25
5142 * fjc = 10e21*1.25 = 12.5e21
5144 fail_unless(fjc(x, T, l, N)-12.5e21 < 1e16,
5145 "fjc(%g, %g, %g, %d) = %g != %g",
5146 x, T, l, N, fjc(x, T, l, N), 12.5e21);
5153 \subsubsection{Piston}
5155 <<piston test case defs>>=
5156 TCase *tc_piston = tcase_create("Piston");
5159 <<piston test case adds>>=
5160 tcase_add_test(tc_piston, test_piston_create_destroy);
5161 suite_add_tcase(s, tc_piston);
5166 START_TEST(test_piston_create_destroy)
5168 char *L[] = {"1","2.2", "3"};
5169 char *params[] = {L[0]};
5172 for( i=0; i < sizeof(L)/sizeof(char *); i++) {
5174 p = string_create_piston_tension_param_t(params);
5175 destroy_piston_tension_param_t(p);
5180 @ TODO: In order to test the piston tension handler itself, we'd
5181 have to construct a group.
5184 \section{Unfolding rate models}
5185 \label{sec.k_models}
5187 We have two main choices for determining $k$: Bell's model, which assumes the
5188 folded domains exist in a single `folded' state and make instantaneous,
5189 stochastic transitions over a free energy barrier; and Kramers' model, which
5190 treats the unfolding as a thermalized diffusion process.
5191 We deal with these options like we dealt with the different tension models: by bundling all model-specific
5192 parameters into structures, with model specific functions for determining $k$.
5194 <<k func definition>>=
5195 typedef double k_func_t(double F, environment_t *env, void *params);
5198 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
5199 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]].
5201 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
5202 because \LaTeX\ doesn't like underscores outside math mode.
5203 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
5204 You could use spaces instead (`\verb+ +'),
5205 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
5206 which I don't like as much.
5212 <<k func definition>>
5213 <<null k declarations>>
5214 <<const k declarations>>
5215 <<bell k declarations>>
5216 <<kbell k declarations>>
5217 <<kramers k declarations>>
5218 <<kramers integ k declarations>>
5219 #endif /* K_MODEL_H */
5222 <<k model module makefile lines>>=
5223 NW_SAWSIM_MODS += k_model
5224 CHECK_BINS += check_k_model
5225 check_k_model_MODS = parse string_eval
5227 [[check_k_model]] also depends on the parse module.
5229 <<k model utils makefile lines>>=
5230 K_MODEL_MODS = k_model parse string_eval
5231 K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
5232 $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h)
5233 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
5234 notangle -Rk-model-utils.c $< > $@
5235 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
5236 $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5237 clean_k_model_utils :
5238 rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
5244 For domains that are never folded.
5246 <<null k declarations>>=
5248 <<null k functions>>=
5253 <<null k model getopt>>=
5254 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
5257 \subsection{Constant rate model}
5260 This is a pretty boring model, but it is usefull for testing the engine.
5264 where $k_0$ is the constant unfolding rate.
5266 <<const k functions>>=
5267 <<const k function>>
5268 <<const k structure create/destroy functions>>
5271 <<const k declarations>>=
5272 <<const k function declaration>>
5273 <<const k structure create/destroy function declarations>>
5274 <<const k global declarations>>
5278 <<const k structure>>=
5279 typedef struct const_k_param_struct {
5280 double knot; /* transition rate k_0 (frac population per s) */
5284 <<const k function declaration>>=
5285 double const_k(double F, environment_t *env, void *const_k_params);
5287 <<const k function>>=
5288 double const_k(double F, environment_t *env, void *const_k_params)
5289 { /* Returns the rate constant k in frac pop per s. */
5290 const_k_param_t *p = (const_k_param_t *) const_k_params;
5292 assert(p->knot > 0);
5297 <<const k structure create/destroy function declarations>>=
5298 void *string_create_const_k_param_t(char **param_strings);
5299 void destroy_const_k_param_t(void *p);
5302 <<const k structure create/destroy functions>>=
5303 const_k_param_t *create_const_k_param_t(double knot)
5305 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
5306 assert(ret != NULL);
5311 void *string_create_const_k_param_t(char **param_strings)
5313 return create_const_k_param_t(safe_strtod(param_strings[0], __FUNCTION__));
5316 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
5323 <<const k global declarations>>=
5324 extern int num_const_k_params;
5325 extern char *const_k_param_descriptions[];
5326 extern char const_k_param_string[];
5329 <<const k globals>>=
5330 int num_const_k_params = 1;
5331 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
5332 char const_k_param_string[]="1";
5335 <<const k model getopt>>=
5336 {"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}
5339 \subsection{Bell's model}
5342 Quantitatively, Bell's model gives $k$ as
5344 k = k_0 \cdot e^{F dx / k_B T} \;,
5346 where $k_0$ is the unforced unfolding rate,
5347 $F$ is the applied unfolding force,
5348 $dx$ is the distance to the transition state, and
5349 $k_B T$ is the thermal energy\citep{schlierf06}.
5351 <<bell k functions>>=
5353 <<bell k structure create/destroy functions>>
5356 <<bell k declarations>>=
5357 <<bell k function declaration>>
5358 <<bell k structure create/destroy function declarations>>
5359 <<bell k global declarations>>
5363 <<bell k structure>>=
5364 typedef struct bell_param_struct {
5365 double knot; /* transition rate k_0 (frac population per s) */
5366 double dx; /* `distance to transition state' (nm) */
5370 <<bell k function declaration>>=
5371 double bell_k(double F, environment_t *env, void *bell_params);
5373 <<bell k function>>=
5374 double bell_k(double F, environment_t *env, void *bell_params)
5375 { /* Returns the rate constant k in frac pop per s.
5376 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5377 * uses global k_B in J/K */
5378 bell_param_t *p = (bell_param_t *) bell_params;
5379 assert(F >= 0); assert(env->T > 0);
5381 assert(p->knot > 0); //assert(p->dx >= 0); /* p < 0 for refolding reactions, etc. */
5383 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
5384 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
5385 p->knot * exp(F*p->dx / (k_B*env->T) ));
5386 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
5388 return p->knot * exp(F*p->dx / (k_B*env->T) );
5392 <<bell k structure create/destroy function declarations>>=
5393 void *string_create_bell_param_t(char **param_strings);
5394 void destroy_bell_param_t(void *p);
5397 <<bell k structure create/destroy functions>>=
5398 bell_param_t *create_bell_param_t(double knot, double dx)
5400 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
5401 assert(ret != NULL);
5407 void *string_create_bell_param_t(char **param_strings)
5409 return create_bell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5410 safe_strtod(param_strings[1], __FUNCTION__));
5413 void destroy_bell_param_t(void *p /* bell_param_t * */)
5420 <<bell k global declarations>>=
5421 extern int num_bell_params;
5422 extern char *bell_param_descriptions[];
5423 extern char bell_param_string[];
5427 int num_bell_params = 2;
5428 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
5429 char bell_param_string[]="3.3e-4,0.25e-9";
5432 <<bell k model getopt>>=
5433 {"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}
5435 Initialized with Titin I27 parameters\citep{carrion-vazquez99a}.
5436 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5439 \subsection{Linker stiffness corrected Bell model}
5442 Walton et. al showed that the Bell model constant force approximation
5443 doesn't hold for biotin-streptavadin binding, and I extended their
5444 results to I27 for the 2009 Biophysical Society Annual
5445 Meeting\cite{walton08,king09}. More details to follow when I get done
5446 with the conference\ldots
5448 We adjust Bell's model to
5450 k = k_0 \cdot e^{\frac{F dx - \frac{1}{2} \kappa dx^2}{k_B T}} \;,
5452 where $k_0$ is the unforced unfolding rate,
5453 $F$ is the applied unfolding force,
5454 $\kappa$ is the effective spring constant,
5455 $dx$ is the distance to the transition state, and
5456 $k_B T$ is the thermal energy\citep{schlierf06}.
5458 <<kbell k functions>>=
5459 <<kbell k function>>
5460 <<kbell k structure create/destroy functions>>
5463 <<kbell k declarations>>=
5464 <<kbell k function declaration>>
5465 <<kbell k structure create/destroy function declarations>>
5466 <<kbell k global declarations>>
5470 <<kbell k structure>>=
5471 typedef struct kbell_param_struct {
5472 double knot; /* transition rate k_0 (frac population per s) */
5473 double dx; /* `distance to transition state' (nm) */
5477 <<kbell k function declaration>>=
5478 double kbell_k(double F, environment_t *env, void *kbell_params);
5480 <<kbell k function>>=
5481 double kbell_k(double F, environment_t *env, void *kbell_params)
5482 { /* Returns the rate constant k in frac pop per s.
5483 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5484 * uses global k_B in J/K */
5485 kbell_param_t *p = (kbell_param_t *) kbell_params;
5486 assert(F >= 0); assert(env->T > 0); assert(env->stiffness >= 0);
5488 assert(p->knot > 0); assert(p->dx >= 0);
5489 return p->knot * exp((F*p->dx - 0.5*env->stiffness*(p->dx)*(p->dx))/(k_B*env->T) );
5493 <<kbell k structure create/destroy function declarations>>=
5494 void *string_create_kbell_param_t(char **param_strings);
5495 void destroy_kbell_param_t(void *p);
5498 <<kbell k structure create/destroy functions>>=
5499 kbell_param_t *create_kbell_param_t(double knot, double dx)
5501 kbell_param_t *ret = (kbell_param_t *) malloc(sizeof(kbell_param_t));
5502 assert(ret != NULL);
5508 void *string_create_kbell_param_t(char **param_strings)
5510 return create_kbell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5511 safe_strtod(param_strings[1], __FUNCTION__));
5514 void destroy_kbell_param_t(void *p /* kbell_param_t * */)
5521 <<kbell k global declarations>>=
5522 extern int num_kbell_params;
5523 extern char *kbell_param_descriptions[];
5524 extern char kbell_param_string[];
5527 <<kbell k globals>>=
5528 int num_kbell_params = 2;
5529 char *kbell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
5530 char kbell_param_string[]="3.3e-4,0.25e-9";
5533 <<kbell k model getopt>>=
5534 {"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}
5536 Initialized with Titin I27 parameters\citep{carrion-vazquez99a}.
5537 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5540 \subsection{Kramer's model}
5543 Kramer's model gives $k$ as
5545 % \frac{1}{k} = \frac{1}{D}
5546 % \int_{x_\text{min}}^{x_\text{max}}
5547 % e^{\frac{-U_F(x)}{k_B T}}
5548 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5551 %where $D$ is the diffusion constant,
5552 %$U_F(x)$ is the free energy along the unfolding coordinate, and
5553 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
5554 % before the minimum and after the maximum, respectively \citep{schlierf06}.
5556 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
5558 where $D$ is the diffusion constant,
5560 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
5562 are length scales where
5563 $x_c(F)$ is the location of the bound state and
5564 $x_{ts}(F)$ is the location of the transition state,
5565 $E(F,x)$ is the free energy, and
5566 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
5567 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
5568 \citet{evans97} Eqn.~3).
5570 <<kramers k functions>>=
5571 <<kramers k function>>
5572 <<kramers k structure create/destroy functions>>
5575 <<kramers k declarations>>=
5576 <<kramers k function declaration>>
5577 <<kramers k structure create/destroy function declarations>>
5578 <<kramers k global declarations>>
5582 <<kramers k structure>>=
5583 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
5584 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
5586 typedef struct kramers_param_struct {
5587 double D; /* diffusion rate (um/s) */
5594 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
5595 //kramers_x_func_t *xts; /* function returning transition x (nm) */
5596 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
5597 //kramers_E_func_t *E; /* function returning E (J) */
5598 //double *E_params; /* parametrize the energy landscape */
5599 //int n_E_params; /* length of E_params array */
5603 <<kramers k function declaration>>=
5604 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
5605 extern double kramers_k(double F, environment_t *env, void *kramers_params);
5607 <<kramers k function>>=
5608 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
5610 static char input[160]={0};
5611 static char output[80]={0};
5612 /* setup the environment */
5613 fprintf(in, "F = %g; T = %g\n", F, T);
5614 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
5615 string_eval(in, out, input, 80, output);
5616 return safe_strtod(output, __FUNCTION__);
5619 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
5621 static char input[160]={0};
5622 static char output[80]={0};
5623 /* setup the environment */
5624 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
5625 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
5626 string_eval(in, out, input, 80, output);
5627 return safe_strtod(output, __FUNCTION__);
5630 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
5632 kramers_param_t *p = (kramers_param_t *) kramers_params;
5633 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
5636 double kramers_k(double F, environment_t *env, void *kramers_params)
5637 { /* Returns the rate constant k in frac pop per s.
5638 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5639 * uses global k_B in J/K */
5640 kramers_param_t *p = (kramers_param_t *) kramers_params;
5641 double kbT, xc, xts, lc, lts, Eb;
5642 assert(F >= 0); assert(env->T > 0);
5645 assert(p->xc != NULL); assert(p->xts != NULL);
5646 assert(p->ddEdxx != NULL); assert(p->E != NULL);
5647 assert(p->in != NULL); assert(p->out != NULL);
5649 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
5650 if (xc == -1.0) return -1.0; /* error (Too much force) */
5651 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
5652 if (xc == -1.0) return -1.0; /* error (Too much force) */
5653 lc = sqrt(2.0*M_PI*kbT /
5654 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
5655 lts = sqrt(-2.0*M_PI*kbT/
5656 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
5657 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
5658 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
5659 //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));
5660 return p->D/(lc*lts) * exp(-Eb/kbT);
5664 <<kramers k structure create/destroy function declarations>>=
5665 void *string_create_kramers_param_t(char **param_strings);
5666 void destroy_kramers_param_t(void *p);
5669 <<kramers k structure create/destroy functions>>=
5670 kramers_param_t *create_kramers_param_t(double D,
5671 char *xc_fn, char *xts_fn,
5672 char *E_fn, char *ddEdxx_fn,
5673 char *global_define_string)
5674 // kramers_x_func_t *xc, kramers_x_func_t *xts,
5675 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
5676 // double *E_params, int n_E_params)
5678 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
5679 assert(ret != NULL);
5680 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5681 assert(ret->xc != NULL);
5682 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5683 assert(ret->xts != NULL);
5684 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
5685 assert(ret->E != NULL);
5686 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
5687 assert(ret->ddEdxx != NULL);
5689 strcpy(ret->xc, xc_fn);
5690 strcpy(ret->xts, xts_fn);
5691 strcpy(ret->E, E_fn);
5692 strcpy(ret->ddEdxx, ddEdxx_fn);
5693 //ret->E_params = E_params;
5694 //ret->n_E_params = n_E_params;
5695 string_eval_setup(&ret->in, &ret->out);
5696 fprintf(ret->in, "kB = %g\n", k_B);
5697 fprintf(ret->in, "%s\n", global_define_string);
5701 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
5702 void *string_create_kramers_param_t(char **param_strings)
5704 return create_kramers_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5712 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
5714 kramers_param_t *pk = (kramers_param_t *)p;
5716 string_eval_teardown(&pk->in, &pk->out);
5718 // free(pk->E_params);
5724 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
5725 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.
5726 We follow \citet{shillcock98} and use
5728 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
5729 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
5732 where TODO, variable meanings.%\citep{shillcock98}.
5733 The first and second derivatives of $E(F,E_b,x,x_s)$ are
5735 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\\
5736 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
5738 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
5739 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
5740 The roots are therefor at
5742 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
5743 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
5744 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
5747 <<kramers k global declarations>>=
5748 extern int num_kramers_params;
5749 extern char *kramers_param_descriptions[];
5750 extern char kramers_param_string[];
5753 <<kramers k globals>>=
5754 int num_kramers_params = 6;
5755 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"};
5756 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)}";
5758 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
5759 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
5760 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
5761 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.
5762 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
5763 It works so long as [[val_1]] is not [[false]].
5765 <<kramers k model getopt>>=
5766 {"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}
5769 \subsection{Kramer's model (natural cubic splines)}
5770 \label{sec.kramers_integ}
5772 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
5773 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
5774 \citet{schlierf06} Eqn.~2)
5776 \frac{1}{k} = \frac{1}{D}
5777 \int_{x_\text{min}}^{x_\text{max}}
5778 e^{\frac{U_F(x)}{k_B T}}
5779 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5782 where $D$ is the diffusion constant,
5783 $U_F(x)$ is the free energy along the unfolding coordinate, and
5784 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
5785 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
5787 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
5788 interpolating between them with cubic splines.
5789 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
5791 <<kramers integ k functions>>=
5792 <<spline functions>>
5793 <<kramers integ A k functions>>
5794 <<kramers integ B k functions>>
5795 <<kramers integ k function>>
5796 <<kramers integ k structure create/destroy functions>>
5799 <<kramers integ k declarations>>=
5800 <<kramers integ k function declaration>>
5801 <<kramers integ k structure create/destroy function declarations>>
5802 <<kramers integ k global declarations>>
5806 <<kramers integ k structure>>=
5807 typedef double func_t(double x, void *params);
5808 typedef struct kramers_integ_param_struct {
5809 double D; /* diffusion rate (um/s) */
5810 double x_min; /* integration bounds */
5812 func_t *E; /* function returning E (J) */
5813 void *E_params; /* parametrize the energy landscape */
5814 destroy_data_func_t *destroy_E_params;
5816 double F; /* for passing into GSL evaluations */
5818 } kramers_integ_param_t;
5821 <<spline param structure>>=
5822 typedef struct spline_param_struct {
5823 int n_knots; /* natural cubic spline knots for U(x) */
5826 gsl_spline *spline; /* GSL spline parameters */
5827 gsl_interp_accel *acc; /* GSL spline acceleration data */
5831 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
5833 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5835 <<kramers integ A k functions>>=
5836 double kramers_integ_k_integrandA(double x, void *params)
5838 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5839 double U = (*p->E)(x, p->E_params);
5840 double Fx = p->F * x;
5841 double kbT = k_B * p->env->T;
5842 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
5843 return exp(-(U-Fx)/kbT);
5845 double kramers_integ_k_integralA(double x, void *params)
5847 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5849 double result, abserr;
5850 size_t neval; /* for qng */
5851 static gsl_integration_workspace *w = NULL;
5852 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5853 f.function = &kramers_integ_k_integrandA;
5855 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5856 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5857 //fprintf(stderr, "integralA = %g\n", result);
5863 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
5864 \int_{x_\text{min}}^{x_\text{max}}
5865 e^{\frac{U_F(x)}{k_B T}}
5866 \text{Integral}_A(x)
5869 <<kramers integ B k functions>>=
5870 double kramers_integ_k_integrandB(double x, void *params)
5872 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5873 double U = (*p->E)(x, p->E_params);
5874 double Fx = p->F * x;
5875 double kbT = k_B * p->env->T;
5876 double intA = kramers_integ_k_integralA(x, params);
5877 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
5878 return exp((U-Fx)/kbT)*intA;
5880 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
5882 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5884 double result, abserr;
5885 size_t neval; /* for qng */
5886 static gsl_integration_workspace *w = NULL;
5887 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5888 f.function = &kramers_integ_k_integrandB;
5892 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5893 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5894 //fprintf(stderr, "integralB = %g\n", result);
5899 With the integrals taken care of, evaluating $k$ is simply
5901 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
5903 <<kramers integ k function declaration>>=
5904 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
5906 <<kramers integ k function>>=
5907 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
5908 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
5909 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
5910 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
5914 <<kramers integ k structure create/destroy function declarations>>=
5915 void *string_create_kramers_integ_param_t(char **param_strings);
5916 void destroy_kramers_integ_param_t(void *p);
5919 <<kramers integ k structure create/destroy functions>>=
5920 kramers_integ_param_t *create_kramers_integ_param_t(double D,
5921 double xmin, double xmax, func_t *E,
5923 destroy_data_func_t *destroy_E_params)
5925 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
5926 assert(ret != NULL);
5931 ret->E_params = E_params;
5932 ret->destroy_E_params = destroy_E_params;
5936 void *string_create_kramers_integ_param_t(char **param_strings)
5938 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
5939 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
5940 void *E_params = string_create_spline_param_t(param_strings+1);
5941 return create_kramers_integ_param_t(
5942 safe_strtod(param_strings[0], __FUNCTION__),
5943 safe_strtod(param_strings[2], __FUNCTION__),
5944 safe_strtod(param_strings[3], __FUNCTION__),
5945 &spline_eval, E_params, destroy_spline_param_t);
5948 void destroy_kramers_integ_param_t(void *params)
5950 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5953 (*p->destroy_E_params)(p->E_params);
5959 Finally we have the GSL spline wrappers:
5960 <<spline functions>>=
5962 <<create/destroy spline>>
5965 <<spline function>>=
5966 double spline_eval(double x, void *spline_params)
5968 spline_param_t *p = (spline_param_t *)(spline_params);
5969 return gsl_spline_eval(p->spline, x, p->acc);
5973 <<create/destroy spline>>=
5974 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
5976 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
5977 assert(ret != NULL);
5978 ret->n_knots = n_knots;
5981 ret->acc = gsl_interp_accel_alloc();
5982 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
5983 gsl_spline_init(ret->spline, x, y, n_knots);
5987 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
5988 void *string_create_spline_param_t(char **param_strings)
5992 double *x=NULL, *y=NULL;
5993 /* break into ordered pairs */
5994 parse_list_string(param_strings[0], SEP, '(', ')',
5995 &num_knots, &knot_coords);
5996 x = (double *)malloc(sizeof(double)*num_knots);
5998 y = (double *)malloc(sizeof(double)*num_knots);
6000 for (i=0; i<num_knots; i++) {
6001 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
6002 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
6005 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
6007 free_parsed_list(num_knots, knot_coords);
6008 return create_spline_param_t(num_knots, x, y);
6011 void destroy_spline_param_t(void *params)
6013 spline_param_t *p = (spline_param_t *)params;
6016 gsl_spline_free(p->spline);
6018 gsl_interp_accel_free(p->acc);
6028 <<kramers integ k global declarations>>=
6029 extern int num_kramers_integ_params;
6030 extern char *kramers_integ_param_descriptions[];
6031 extern char kramers_integ_param_string[];
6034 <<kramers integ k globals>>=
6035 int num_kramers_integ_params = 4;
6036 char *kramers_integ_param_descriptions[] = {
6037 "diffusion D, m^2 / s",
6038 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
6039 "starting integration bound x_min, m",
6040 "ending integration bound x_max, m"};
6041 //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";
6042 //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";
6043 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";
6044 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
6045 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
6046 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
6049 <<kramers integ k model getopt>>=
6050 {"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}
6052 Initialized with Titin I27 parameters\citep{carrion-vazquez99a}.
6053 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
6055 \subsection{Unfolding model implementation}
6059 <<k model includes>>
6060 #include "k_model.h"
6061 <<k model internal definitions>>
6063 <<k model functions>>
6066 <<k model includes>>=
6067 #include <assert.h> /* assert() */
6068 #include <stdlib.h> /* NULL, malloc(), strto*() */
6069 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
6070 #include <stdio.h> /* fprintf(), stdout */
6071 #include <string.h> /* strlen(), strcpy() */
6072 #include <errno.h> /* errno, ERANGE */
6073 #include <gsl/gsl_integration.h>
6074 #include <gsl/gsl_spline.h>
6079 <<k model internal definitions>>=
6080 <<const k structure>>
6081 <<bell k structure>>
6082 <<kbell k structure>>
6083 <<kramers k structure>>
6084 <<kramers integ k structure>>
6085 <<spline param structure>>
6088 <<k model globals>>=
6093 <<kramers k globals>>
6094 <<kramers integ k globals>>
6097 <<k model functions>>=
6099 <<null k functions>>
6100 <<const k functions>>
6101 <<bell k functions>>
6102 <<kbell k functions>>
6103 <<kramers k functions>>
6104 <<kramers integ k functions>>
6107 \subsection{Unfolding model unit tests}
6109 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
6110 <<check-k-model.c>>=
6112 <<k model check includes>>
6113 <<check relative error>>
6115 <<k model test suite>>
6116 <<main check program>>
6119 <<k model check includes>>=
6120 #include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
6121 #include <stdio.h> /* sprintf() */
6122 #include <assert.h> /* assert() */
6123 #include <math.h> /* exp() */
6124 #include <errno.h> /* errno, ERANGE */
6127 #include "k_model.h"
6130 <<k model test suite>>=
6134 <<k model suite definition>>
6137 <<k model suite definition>>=
6138 Suite *test_suite (void)
6140 Suite *s = suite_create ("k model");
6141 <<const k test case defs>>
6142 <<bell k test case defs>>
6144 <<const k test case adds>>
6145 <<bell k test case adds>>
6150 \subsubsection{Constant}
6152 <<const k test case defs>>=
6153 TCase *tc_const_k = tcase_create("Constant k");
6156 <<const k test case adds>>=
6157 tcase_add_test(tc_const_k, test_const_k_create_destroy);
6158 tcase_add_test(tc_const_k, test_const_k_over_range);
6159 suite_add_tcase(s, tc_const_k);
6164 START_TEST(test_const_k_create_destroy)
6166 char *knot[] = {"1","2","3","4","5","6"};
6167 char *params[] = {knot[0]};
6170 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
6171 params[0] = knot[i];
6172 p = string_create_const_k_param_t(params);
6173 destroy_const_k_param_t(p);
6178 START_TEST(test_const_k_over_range)
6180 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
6181 char *knot[] = {"1","2","3","4","5","6"};
6182 char *params[] = {knot[0]};
6185 char param_string[80];
6188 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
6189 params[0] = knot[i];
6190 p = string_create_const_k_param_t(params);
6191 for ( F=Fm; F<FM; F+=dF ) {
6192 fail_unless(const_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
6193 "const_k(%g, %g, &{%s}) = %g != %s",
6194 F, T, knot[i], const_k(F, &env, p), knot[i]);
6196 destroy_const_k_param_t(p);
6202 \subsubsection{Bell}
6204 <<bell k test case defs>>=
6205 TCase *tc_bell = tcase_create("Bell's k");
6208 <<bell k test case adds>>=
6209 tcase_add_test(tc_bell, test_bell_k_create_destroy);
6210 tcase_add_test(tc_bell, test_bell_k_at_zero);
6211 tcase_add_test(tc_bell, test_bell_k_at_one);
6212 suite_add_tcase(s, tc_bell);
6216 START_TEST(test_bell_k_create_destroy)
6218 char *knot[] = {"1","2","3","4","5","6"};
6220 char *params[] = {knot[0], dx};
6223 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
6224 params[0] = knot[i];
6225 p = string_create_bell_param_t(params);
6226 destroy_bell_param_t(p);
6231 START_TEST(test_bell_k_at_zero)
6233 double F=0.0, T=300.0;
6234 char *knot[] = {"1","2","3","4","5","6"};
6236 char *params[] = {knot[0], dx};
6239 char param_string[80];
6242 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
6243 params[0] = knot[i];
6244 p = string_create_bell_param_t(params);
6245 fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
6246 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
6247 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
6248 destroy_bell_param_t(p);
6253 START_TEST(test_bell_k_at_one)
6256 char *knot[] = {"1","2","3","4","5","6"};
6258 char *params[] = {knot[0], dx};
6259 double F= k_B*T/safe_strtod(dx, __FUNCTION__);
6262 char param_string[80];
6265 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
6266 params[0] = knot[i];
6267 p = string_create_bell_param_t(params);
6268 CHECK_ERR(1e-6, safe_strtod(knot[i], __FUNCTION__)*exp(1.0), bell_k(F, &env, p));
6269 //fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__)*exp(1.0),
6270 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
6271 // F, T, knot[i], dx, bell_k(F, &env, p), safe_strtod(knot[i], __FUNCTION__)*exp(1.0));
6272 destroy_bell_param_t(p);
6278 <<kramers k tests>>=
6281 <<kramers k test case defs>>=
6284 <<kramers k test case adds>>=
6287 <<k model function tests>>=
6288 <<check relative error>>
6290 START_TEST(test_something)
6292 double k=5, x=3, last_x=2.0, F;
6293 one_dim_fn_t *handlers[] = {&hooke};
6294 void *data[] = {&k};
6295 double xi[] = {0.0};
6297 int new_active_groups = 1;
6298 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
6299 fail_unless(F = k*x, NULL);
6304 \subsection{Utilities}
6306 The unfolding models can get complicated, and are sometimes hard to explain to others.
6307 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
6308 the overhead of having to understand a full simulation.
6309 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
6310 or special mode, where the operation depends on the specific model selected.
6312 <<k-model-utils.c>>=
6314 <<k model utility includes>>
6315 <<k model utility definitions>>
6316 <<k model utility getopt functions>>
6317 <<k model utility multi model E>>
6318 <<k model utility main>>
6321 <<k model utility main>>=
6322 int main(int argc, char **argv)
6324 <<k model getopt array definition>>
6325 k_model_getopt_t *model = NULL;
6330 int num_param_args; /* for INIT_MODEL() */
6331 char **param_args; /* for INIT_MODEL() */
6333 double special_xmin;
6334 double special_xmax;
6335 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
6336 &Fmax, &special_xmin, &special_xmax, &flags);
6337 if (flags & VFLAG) {
6338 printf("#initializing model %s with parameters %s\n", model->name, model->params);
6340 INIT_MODEL("utility", model, model->params, params);
6343 if (model->k == NULL) {
6344 printf("No k function for model %s\n", model->name);
6348 printf("Fmax = %g <= 0\n", Fmax);
6354 printf("#F (N)\tk (%% pop. per s)\n");
6355 for (i=0; i<=N; i++) {
6356 F = Fmax*i/(double)N;
6357 k = (*model->k)(F, &env, params);
6359 printf("%g\t%g\n", F, k);
6364 if (model->k == NULL || model->k == &bell_k) {
6365 printf("No special mode for model %s\n", model->name);
6368 if (special_xmax <= special_xmin) {
6369 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
6373 double Xrng=(special_xmax-special_xmin), x, E;
6375 printf("#x (m)\tE (J)\n");
6376 for (i=0; i<=N; i++) {
6377 x = special_xmin + Xrng*i/(double)N;
6378 E = multi_model_E(model->k, params, &env, x);
6379 printf("%g\t%g\n", x, E);
6386 if (model->destructor)
6387 (*model->destructor)(params);
6392 <<k model utility multi model E>>=
6393 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
6395 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
6397 if (k == kramers_integ_k)
6398 E = (*p->E)(x, p->E_params);
6400 E = kramers_E(0, env, model_params, x);
6406 <<k model utility includes>>=
6409 #include <string.h> /* strlen() */
6410 #include <assert.h> /* assert() */
6411 #include <errno.h> /* errno, ERANGE */
6414 #include "string_eval.h"
6415 #include "k_model.h"
6418 <<k model utility definitions>>=
6419 <<version definition>>
6420 #define VFLAG 1 /* verbose */
6421 enum mode_t {M_K_OF_F, M_SPECIAL};
6422 <<string comparison definition and globals>>
6423 <<k model getopt definitions>>
6424 <<initialize model definition>>
6425 <<kramers k structure>>
6426 <<kramers integ k structure>>
6430 <<k model utility getopt functions>>=
6433 <<k model utility help>>
6434 <<k model utility get options>>
6437 <<k model utility help>>=
6438 void help(char *prog_name,
6440 int n_k_models, k_model_getopt_t *k_models,
6441 int k_model, double Fmax, double special_xmin, double special_xmax)
6444 printf("usage: %s [options]\n", prog_name);
6445 printf("Version %s\n\n", VERSION);
6446 printf("A utility for understanding the available unfolding models\n\n");
6447 printf("Environmental options:\n");
6448 printf("-T T\tTemperature T (currently %g K)\n", env->T);
6449 printf("-C T\tYou can also set the temperature T in Celsius\n");
6450 printf("Model options:\n");
6451 printf("-k model\tTransition rate model (currently %s)\n",
6452 k_models[k_model].name);
6453 printf("-K args\tTransition rate model argument string (currently %s)\n",
6454 k_models[k_model].params);
6455 printf("There are two output modes. In standard mode, k(F) is printed\n");
6456 printf("For example:\n");
6457 printf(" #Force (N)\tk (%% pop. per s)\n");
6458 printf(" 123.456\t7.89\n");
6460 printf("In special mode, the output depends on the model.\n");
6461 printf("For models defining an energy function E(x), that function is printed\n");
6462 printf(" #Distance (m)\tE (J)\n");
6463 printf(" 0.0012\t0.0034\n");
6465 printf("-m\tChange output to standard mode\n");
6466 printf("-M\tChange output to special mode\n");
6467 printf("-F\tSet the maximum F value for the standard mode k(F) output (currently %g)\n", Fmax);
6468 printf("-x\tSet the minimum x value for the special mode E(x) output (currently %g)\n", special_xmin);
6469 printf("-X\tSet the maximum x value for the special mode E(x) output (currently %g)\n", special_xmax);
6470 printf("-V\tChange output to verbose mode\n");
6471 printf("-h\tPrint this help and exit\n");
6473 printf("Unfolding rate models:\n");
6474 for (i=0; i<n_k_models; i++) {
6475 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
6476 for (j=0; j < k_models[i].num_params ; j++)
6477 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
6478 printf("\t\tdefaults: %s\n", k_models[i].params);
6483 <<k model utility get options>>=
6484 void get_options(int argc, char **argv, environment_t *env,
6485 int n_k_models, k_model_getopt_t *k_models,
6486 enum mode_t *mode, k_model_getopt_t **model,
6487 double *Fmax, double *special_xmin, double *special_xmax,
6488 unsigned int *flags)
6490 char *prog_name = NULL;
6491 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
6493 extern char *optarg;
6494 extern int optind, optopt, opterr;
6498 /* setup defaults */
6500 prog_name = argv[0];
6501 env->T = 300.0; /* K */
6505 *Fmax = 1e-9; /* N */
6506 *special_xmax = 1e-8;
6507 *special_xmin = 0.0;
6509 while ((c=getopt(argc, argv, options)) != -1) {
6511 case 'T': env->T = safe_strtod(optarg, "-T"); break;
6512 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
6514 k_model = index_k_model(n_k_models, k_models, optarg);
6515 *model = k_models+k_model;
6518 k_models[k_model].params = optarg;
6520 case 'm': *mode = M_K_OF_F; break;
6521 case 'M': *mode = M_SPECIAL; break;
6522 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
6523 case 'x': *special_xmin = safe_strtod(optarg, "-x"); break;
6524 case 'X': *special_xmax = safe_strtod(optarg, "-X"); break;
6525 case 'V': *flags |= VFLAG; break;
6527 fprintf(stderr, "unrecognized option '%c'\n", optopt);
6528 /* fall through to default case */
6530 help(prog_name, env, n_k_models, k_models, k_model, *Fmax, *special_xmin, *special_xmax);
6539 \section{\LaTeX\ documentation}
6541 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
6542 The comment blocks are extracted (with nicely formatted code blocks), using
6543 <<latex makefile lines>>=
6544 $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6545 noweave -latex -index -delay $< > $@
6546 $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR)
6550 <<latex makefile lines>>=
6551 $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
6553 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
6554 $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
6555 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
6556 $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
6557 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
6558 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
6559 mv $(BUILD_DIR)/sawsim.pdf $@
6561 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
6562 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
6563 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
6565 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
6566 <<latex makefile lines>>=
6568 rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
6569 $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
6570 $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.toc \
6571 $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib
6572 clean_latex : semi-clean_latex
6573 rm -f $(DOC_DIR)/sawsim.pdf
6578 [[make]] is a common utility on *nix systems for managing dependencies while building software.
6579 In this case, we have several \emph{targets} that we'd like to build:
6580 the main simulation program \prog;
6581 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
6582 the documentation [[sawsim.pdf]];
6583 and the [[Makefile]] itself.
6584 Besides the generated files, there is the utility target
6585 [[clean]] for removing all generated files except [[Makefile]].
6587 % [[dist]] for generating a distributable tar file.
6589 Extract the makefile with
6590 `[[notangle -Rmakefile src/sawsim.nw | sed 's/ /\t/' > Makefile]]'.
6591 The sed is needed because notangle replaces tabs with eight spaces,
6592 but [[make]] needs tabs.
6594 # Customize compilation
6599 # Decide what directories to put things in
6604 # Note: Cannot use, for example, `./build', because make eats the `./'
6605 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
6607 # Modules (X.c, X.h) defined in the noweb file
6609 # Check binaries (check_*) defined in the noweb file
6611 # Modules defined outside the noweb file
6612 FREE_SAWSIM_MODS = interp tavl
6614 <<list module makefile lines>>
6615 <<tension balance module makefile lines>>
6616 <<tension model module makefile lines>>
6617 <<k model module makefile lines>>
6618 <<parse module makefile lines>>
6619 <<string eval module makefile lines>>
6620 <<accel k module makefile lines>>
6622 SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS)
6624 # Everything needed for sawsim under one roof. sawsim.c must come first
6625 SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \
6626 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h)
6627 # Libraries needed to compile sawsim
6628 LIBS = gsl gslcblas m
6629 CHECK_LIBS = gsl gslcblas m check
6630 # The non-check binaries generated
6631 BINS = sawsim tension_model_utils k_model_utils sawsim_profile
6634 # Define the major targets
6635 all : ./Makefile all_bin all_doc
6637 .PHONY: all_bin all_doc
6638 all_bin : $(BINS:%=$(BIN_DIR)/%)
6639 all_doc : $(DOCS:%=$(DOC_DIR)/%)
6641 view : $(DOC_DIR)/sawsim.pdf
6643 profile : $(BIN_DIR)/sawsim_profile | $(BIN_DIR)
6644 $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \
6645 -mnull -Mwlc -A0.39e-9,28e-9 -F8
6646 gprof $(BIN_DIR)/sawsim_profile gmon.out > $@
6647 check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim
6648 $(SHELL) -e -c 'for B in $^; do ./$$B; done'
6649 clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \
6650 clean_tension_model_utils \
6651 clean_k_model_utils clean_latex clean_check_sawsim
6652 rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \
6653 $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \
6654 $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
6655 $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
6656 $(BUILD_DIR)/global.h ./gmon.out
6657 $(SHELL) -e -c "rmdir $(BUILD_DIR) $(DOC_DIR)"
6659 # Various builds of sawsim
6660 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
6661 $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6662 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
6663 $(CC) $(CFLAGS) $(LDFLAGS) -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6664 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
6665 $(CC) $(CFLAGS) $(LDFLAGS) -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6667 # Create the directories
6668 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
6671 # Copy over the source external to sawsim
6672 # Note: Cannot use, for example, `./build', because make eats the `./'
6673 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
6675 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6679 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6683 ## Basic source generated with noweb
6684 # The central files sawsim.c and global.h...
6685 $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6687 $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6688 notangle -Rglobal.h $< > $@
6689 # ... and the modules
6690 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6691 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6692 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6693 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6694 $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6695 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6696 # Note: I use `_' as a space in my file names, but noweb doesn't like
6697 # that so we use `-' to name the noweb roots and substitute here.
6699 ## Building the unit-test programs
6701 $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
6702 notangle -Rchecks $< > $@
6703 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
6704 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
6705 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR)
6706 $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
6707 clean_check_sawsim :
6708 rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
6709 # ... and the modules
6711 $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
6712 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \
6713 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \
6714 $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6715 $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6716 $$(BUILD_DIR)/global.h | $(BIN_DIR)
6717 $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
6718 $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
6720 # TODO: clean up the required modules too
6721 $(CHECK_BINS:%=clean_%) :
6722 rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
6724 # Cleaning up the modules
6726 $(SAWSIM_MODS:%=clean_%) :
6727 rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h)
6729 <<tension model utils makefile lines>>
6730 <<k model utils makefile lines>>
6731 <<latex makefile lines>>
6733 Makefile : $(SRC_DIR)/sawsim.nw
6734 notangle -Rmakefile $< | sed 's/ /\t/' > $@
6736 This is fairly self-explanatory. We compile a dynamically linked
6737 version ([[sawsim]]) and a statically linked version
6738 ([[sawsim_static]]). The static version is about 9 times as big, but
6739 it runs without needing dynamic GSL libraries to link to, so it's a
6740 better format for distributing to a cluster for parallel evaluation.
6744 \subsection{Hookean springs in series}
6745 \label{sec.math_hooke}
6748 The effective spring constant for $n$ springs $k_i$ in series is given by
6750 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6756 F &= k_i x_i &&\forall i \le n \\
6757 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
6758 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
6759 F &= k_1 x_1 = k_\text{eff} x \\
6760 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
6761 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6766 \addcontentsline{toc}{section}{References}
6767 \bibliography{sawsim}