1 % sawsim - simulating single molecule mechanical unfolding.
2 % Copyright (C) 2008-2010 William Trevor King
4 % This program is free software: you can redistribute it and/or modify
5 % it under the terms of the GNU General Public License as published by
6 % the Free Software Foundation, either version 3 of the License, or
7 % (at your option) any later version.
9 % This program is distributed in the hope that it will be useful,
10 % but WITHOUT ANY WARRANTY; without even the implied warranty of
11 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 % GNU General Public License for more details.
14 % You should have received a copy of the GNU General Public License
15 % along with this program. If not, see <http://www.gnu.org/licenses/>.
17 % The author may be contacted at <wking@drexel.edu> on the Internet, or
18 % write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
19 % Philadelphia PA 19104, USA.
21 %%%%%%%%%%%%%%%%%%%%%%%%% Start LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%
22 % we have our own preamble,
23 % use `noweave -delay` to not print noweb's automatic one
24 % -*- mode: Noweb; noweb-code-mode: c-mode -*-
25 \documentclass[letterpaper, 11pt]{article}
28 \noweboptions{smallcode}
30 \usepackage{url} % needed for natbib for underscores in urls?
31 \usepackage[breaklinks,colorlinks]{hyperref} % for internal links, \phantomsection
32 % breaklinks breaks long links
33 % colorlinks colors link text instead of boxing it
34 \usepackage{amsmath} % for \text{}, \xrightarrow[sub]{super}
35 \usepackage{amsthm} % for theorem and proof environments
36 \newtheorem{theorem}{Theorem}
38 \usepackage{subfig} % Support side-by-side figures
39 \usepackage{pgf} % Fancy graphics
40 \usepackage{tikz} % A nice, inline PGF frontend
41 \usetikzlibrary{automata} % Graph-theory library
43 \usepackage{doc} % BibTeX
44 \usepackage[super,sort&compress]{natbib} % fancy citation extensions
45 % super selects citations in superscript mode
46 % sort&compress automatically sorts and compresses compound citations (\cite{a,b,...})
48 \usepackage[mathletters]{ucs} % use math fonts to allow Greek UTF-8 in the text
49 \usepackage[utf8x]{inputenc} % allow UTF-8 characters
51 %\bibliographystyle{ieeetr} % pick the bibliography style, short and sweet
52 %\bibliographystyle{plain} % pick the bibliography style, includes dates
53 \bibliographystyle{plainnat}
54 \defcitealias{sw:noweb}{{\tt noweb}}
55 \defcitealias{galassi05}{GSL}
56 \defcitealias{sw:check}{{\tt check}}
57 \defcitealias{sw:python}{Python}
62 \textheight 9.5in % leave a bit of extra room for page numbers (at top of page)
65 \setlength{\parindent}{0pt}
66 \setlength{\parskip}{5pt}
68 % For some reason, the \maketitle wants to go on it's own page
69 % so we'll just hardcode the following in our first page.
70 %\title{Sawsim: a sawtooth protein unfolding simulator}
71 %\author{W.~Trevor King}
74 \newcommand{\prog}{[[sawsim]]}
75 \newcommand{\lang}{[[C]]}
76 \newcommand{\p}[3]{\left#1 #2 \right#3} % e.g. \p({complicated expr})
77 \newcommand{\Th}{\ensuremath{\sup{\text{th}}}} % ordinal th
78 \newcommand{\U}[1]{\ensuremath{\text{ #1}}} % units: $1\U{m/s}$
80 % single spaced lists, from Alvin Alexander
81 % http://www.devdaily.com/blog/post/latex/control-line-spacing-in-itemize-enumerate-tags/
82 \newenvironment{packed_item}{
84 \setlength{\itemsep}{1pt}
85 \setlength{\parskip}{0pt}
86 \setlength{\parsep}{0pt}
90 \nwfilename{sawsim.nw}
95 Sawsim: a sawtooth protein unfolding simulator \\
103 %%%%%%%%%%%%%%%%%%%%%%%%%%%% End LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%%
105 \section{Introduction}
107 The unfolding of multi-globular protein strings is a stochastic
108 process. In the \prog\ program, we use Monte Carlo methods to
109 simulate this unfolding through an explicit time-stepping approach.
110 We develop a program in \lang\ to simulate probability of unfolding
111 under a constant extension velocity of a chain of globular domains.
113 In order to extract physical meaning from many mechanical unfolding
114 simulations, it is necessary to compare the experimental results
115 against the results of simulations using different models and
116 parameters. The model and parameters that best match the experimental
117 data provide the ‘best’ model for that experiment.
119 Previous mechanical unfolding studies have rolled their own
120 simulations for modelling their experiment, and there has been little
121 exchange and standardization of these simulations between groups.
122 The \prog\ program attempts to fill this gap, providing a flexible,
123 extensible, \emph{community} program, to make it easier for groups
124 to share the burden and increase the transparency of data analysis.
126 ‘Sawsim’ is blend of ‘sawtooth simulation’.
128 \subsection{Related work}
132 \subsection{About this document}
134 This document was written in \citetalias{sw:noweb} with code blocks in
135 \lang\ and comment blocks in \LaTeX.
137 \subsection{System Requirements and Dependencies}
140 \subsection{Change Log}
142 Version 0.0 used only the unfolded domain WLC to determine the tension
143 and had a fixed timestep $dt$.
145 Version 0.1 added adaptive timesteps, adjusting so that the unfolding
146 probability for a given domain was constant.
148 Version 0.2 added Kramers' model unfolding rate calculations, and a
149 switch to select Bell's or Kramers' model. This induced a major
150 rewrite, introducing the tension group abstraction, and a split to
151 multiple \lang\ source files. Also added a unit-test suites for
152 sanity checking, and switched to SI units throughout.
154 Version 0.3 added integrated Kramers' model unfolding rate
155 calculations. Also replaced some of my hand-coded routines with GNU
156 Scientific Library \citepalias{galassi05} calls.
158 Version 0.4 added Python evaluation for the saddle-point Kramers'
159 model unfolding rate calculations, so that the model functions could
160 be controlled from the command line. Also added interpolating lookup
161 tables to accelerate slow unfolding rate model evaluation, and command
162 line control of tension groups.
164 Version 0.5 added the stiffness environmental parameter and it's
165 associated unfolding models.
167 Version 0.6 generalized from two state unfolding to arbitrary
168 $N$-state rate reactions. This allows simulations of
169 unfolding-refolding, metastable intermediates, etc., but required
170 reorganizing the command line interface.
172 Version 0.7 added tension model inverses, which seems to reduce
173 computation time by about a factor of 2. I was expecting more, but
174 I'll take what I can get.
176 Version 0.8 fixed a minor bug in unfolding probability for
177 multi-domain groups. The probability of at least one domain unfolding
178 had been calculated as $P_N=NP_1$, rather than $P_N=1-(1-P_1)^N$.
179 However, for small $P$ the two are equivalent.
181 Version 0.9 added the [[-P]] option to sawsim, setting the target
182 probability for the per-state transition $P_N$, not the per-domain
185 Version 0.10 added the [[-F]] option to sawsim, setting a limit on the
186 force change $dF$ in a single timestep for continuous pulling
187 simulations. It also added the piston tension model (Section
188 \ref{sec.piston_tension}), and adjusted the stiffness calculations to
189 handle domains with undefined stiffness.
191 <<version definition>>=
192 #define VERSION "0.10"
198 sawsim - simulating single molecule mechanical unfolding.
199 Copyright (C) 2008-2010, William Trevor King
201 This program is free software; you can redistribute it and/or
202 modify it under the terms of the GNU General Public License as
203 published by the Free Software Foundation; either version 3 of the
204 License, or (at your option) any later version.
206 This program is distributed in the hope that it will be useful, but
207 WITHOUT ANY WARRANTY; without even the implied warranty of
208 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
209 See the GNU General Public License for more details.
211 You should have received a copy of the GNU General Public License
212 along with this program. If not, see <http://www.gnu.org/licenses/>.
214 The author may be contacted at <wking@drexel.edu> on the Internet, or
215 write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
216 Philadelphia PA 19104, USA.
229 The unfolding system is stored as a chain of \emph{domains} (Figure
230 \ref{fig.domain_chain}). Domains can be folded, globular protein
231 domains, unfolded protein linkers, AFM cantilevers, or other
232 stretchable system components. The system is modeled as graph of
233 possible domain states with transitions between them (Figure
234 \ref{fig.domain_states}).
238 \subfloat[][]{\label{fig.domain_chain}
239 \begin{tikzpicture}[->,node distance=2cm,font=\footnotesize]
240 \tikzstyle{every state}=[fill,draw=red!50,very thick,fill=red!20]
241 \node[state] (A) {domain 1};
242 \node[state] (B) [below of=A] {domain 2};
243 \node[state] (C) [below of=B] {.~.~.};
244 \node[state] (D) [below of=C] {domain $N$};
245 \node (S) [below of=D] {Surface};
246 \node (E) [above of=A] {};
248 \path[-] (A) edge (B)
249 (B) edge node [right] {Tension} (C)
252 \path[->,green] (A) edge node [right,black] {Extension} (E);
256 \subfloat[][]{\label{fig.domain_states}
257 \begin{tikzpicture}[->,node distance=2.5cm,shorten <=1pt,shorten >=1pt,font=\footnotesize]
258 \tikzstyle{every state}=[fill,draw=blue!50,very thick,fill=blue!20]
259 \node[state] (A) {cantilever};
260 \node[state] (C) [below of=A] {transition};
261 \node[state] (B) [left of=C] {folded};
262 \node[state] (D) [right of=C] {unfolded};
264 \path (B) edge [bend left] node [above] {$k_1$} (C)
265 (C) edge [bend left] node [below] {$k_1'$} (B)
266 edge [bend left] node [above] {$k_2$} (D)
267 (D) edge [bend left] node [below] {$k_2'$} (C);
270 \caption{\subref{fig.domain_chain} Extending a chain of domains. One end
271 of the chain is fixed, while the other is extended at a constant
272 speed. The domains are coupled with rigid linkers, so the domains
273 themselves must stretch to accomodate the extension.
274 \subref{fig.domain_states} Each domain exists in a discrete state. At
275 each timestep, it may transition into another state following a
276 user-defined state matrix such as this one, showing a metastable
277 transition state and an explicit ``cantilever'' domain.}
281 Each domain contributes to the total tension, which depends on the
282 chain extension. The particular model for the tension as a function
283 of extension is handled generally with function pointers. So far the
284 following models have been implemented:
286 \item Null (Appendix \ref{sec.null_tension}),
287 \item Constant (Appendix \ref{sec.const_tension}),
288 \item Hookean spring (Appendix \ref{sec.hooke}),
289 \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
290 \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
291 \item Piston (Appendix \ref{sec.piston_tension}),
294 The tension model and parameters used for a particular domain depend
295 on the domain's current state. The transition rates between states
296 are also handled generally with function pointers, with
299 \item Null (Appendix \ref{sec.null_k}),
300 \item Bell's model (Appendix \ref{sec.bell}),
301 \item Spring constant corrected Bell model (Appendix \ref{sec.kbell}),
302 \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
303 \item Kramers' saddle point model (Appendix \ref{sec.kramers}).
306 The unfolding of the chain domains is modeled in two stages. First
307 the tension of the chain is calculated. Then the tension (assumed to
308 be constant over the length of the chain, see Section
309 \ref{sec.timescales}), is applied to each domain, and used to
310 calculate the probability of that domain changing state in the next
311 timestep $dt$. Then the time is stepped forward, and the process is
312 repeated. The simulation is complete when a pre-set time limit
313 [[t_max]] is reached or a pre-selected state [[stop_state]] is empty.
315 int main(int argc, char **argv)
317 <<tension model getopt array definition>>
318 <<k model getopt array definition>>
320 /* variables initialized in get_options() */
321 list_t *state_list=NULL, *transition_list=NULL;
322 environment_t env={0};
323 double P, dF, t_max, dt_max, v, x_step;
324 state_t *stop_state=NULL;
326 /* variables used in the simulation loop */
327 double t=0, x=0, dt, F; /* time, distance, timestep, force */
328 int transition=1; /* boolean marking a transition for tension and timestep optimization */
330 get_options(argc, argv, &P, &dF, &t_max, &dt_max, &v, &x_step,
331 &stop_state, &env, NUM_TENSION_MODELS, tension_models,
332 NUM_K_MODELS, k_models, &state_list, &transition_list);
335 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Distance (m)\tForce (N)\tStiffness (N/m)\n");
336 if (flags & FLAG_OUTPUT_TRANSITION_FORCES) printf("#Force (N)\tInitial state\tFinal state\n");
337 while ((t_max > 0 && t <= t_max) || (stop_state != NULL && stop_state->num_domains > 0)) {
338 F = find_tension(state_list, &env, x, transition, 0);
340 dt = determine_dt(state_list, transition_list, &env, F, x, dF, P, dt_max, v, transition);
342 dt = determine_dt(state_list, transition_list, &env, F, x, dF, P, dt_max, 0, transition);
343 transition=random_transitions(transition_list, F, dt, &env);
344 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness);
349 if (dt == dt_max) { /* step completed */
352 } else { /* still working on this step */
357 destroy_state_list(state_list);
358 destroy_transition_list(transition_list);
362 @ The meat of the simulation is bundled into the three functions
363 [[find_tension]], [[determine_dt]], and [[random_transitions]].
364 [[find_tension]] is discussed in Section \ref{sec.find_tension},
365 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
366 [[random_transitions]] in Section \ref{sec.transition_rate}.
368 The stretched distance is extended in one of two modes depending on
369 [[xstep]]. If [[xstep == 0]], then the pulling is continuous, at
370 least within the limits of the inherent discretization of a time
371 stepped approach. At any rate, the timestep is limited by the maximum
372 allowed timestep [[dt_max]] and the maximum allowed unfolding
373 probability in a given step [[P]]. In the second mode [[xstep > 0]],
374 and the end to end distance increases in discrete steps of that
375 length. The time between steps is chosen to maintain an average
376 unfolding velocity [[v]]. This approach is less work to simulate than
377 smooth pulling and also closer to the reality of many experiments, but
378 it is practically impossible to treat analytically. With both pulling
379 styles implemented, the effects of the discretization can be easily
380 calculated, bridging the gap between theory and experiment.
382 Environmental parameters important in determining reaction rates and
383 tensions (e.g.\ temperature, pH) are stored in a single structure to
384 facilitate extension to more complicated models in the future.
385 <<environment definition>>=
386 typedef struct environment_struct {
396 #define DOUBLE_PRECISION 1e-12
398 <<environment definition>>
399 <<one dimensional function definition>>
400 <<create/destroy definitions>>
402 #endif /* GLOBAL_H */
404 Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
405 included multiple times.
407 \section{Simulation functions}
409 <<simulation functions>>=
413 <<does domain transition>>
414 <<randomly transition domains>>
418 \label{sec.find_tension}
420 Because the stretched system may be made up of several parts (folded
421 domains, unfolded domains, spring-like cantilever, \ldots), we will
422 assign the domains to states with tension models and groups. The
423 models are used to determine the tension of all the domain in a given
424 state following some model (Hookean spring, WLC, \ldots). The domains
425 are assumed to be commutative, so ordering is ignored. The
426 interactions between the groups are assumed to be linear, but the
427 interactions between domains of the same group need not be. Each
428 model has a tension handler function, which gives the tension $F$
429 needed to stretch a given group of domains a distance $x$.
431 Groups represent collections of states which the model should treat as
432 a single entity. To understand the purpose of groups, consider a
433 system with two unfolded domain states (e.g.\ for two different
434 proteins) that were both modeled as WLCs. If both unfolded states had
435 the same persistence length, it would be wise to assign them to the
436 same group. This would reduce both the computational cost of
437 balancing the tension and the error associated with the inter-group
438 interaction linearity. Note that group numbers only matter
439 \emph{within} model classes, since grouping domains with seperate
440 models doesn't make sense. Within designated groups, the tension
441 parameters for each domain are still checked for compatibility, so if
442 you accidentally assign incompatible domains to the same group you
443 should get an appropriate error message.
445 <<find tension definitions>>=
446 #define NUM_TENSION_MODELS 6
450 <<tension handler array global declaration>>=
451 extern one_dim_fn_t *tension_handlers[];
454 <<tension handler array global>>=
455 one_dim_fn_t *tension_handlers[] = {
457 &const_tension_handler,
465 <<tension model global declarations>>=
466 <<tension handler array global declaration>>
469 <<tension handler types>>=
470 typedef struct tension_handler_data_struct {
471 list_t *group_tension_params; /* one item for each domain in the group */
474 } tension_handler_data_t;
478 After sorting the chain into separate groups $G_i$, with tension
479 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
480 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
481 \forall i,j$. Note that there may be several states within each
482 group. [[find_tension]] is basically a wrapper to feed proper input
483 into the [[tension_balance]] function. [[transition]] should be set
484 to 0 if there were no transitions since the previous call to
485 [[find_tension]] to support some optimizations that assume a small
486 increase in tension between steps (only true if no transition has
487 occured). While we're messing with the tension handlers and if
488 [[const_env==0]], we also grab a measure of the current linker
489 stiffness for the environmental variable, which is needed by some of
490 the unfolding rate models (e.g.\ the linker-stiffness corrected Bell
491 model, Appendix \ref{sec.kbell}).
493 double find_tension(list_t *state_list, environment_t *env, double x,
494 int transition, int const_env)
495 { // TODO: !cont_env -> step_call, only alter env or statics if step_call==1
496 static int num_active_groups=0;
497 static one_dim_fn_t **per_group_handlers = NULL;
498 static one_dim_fn_t **per_group_inverse_handlers = NULL;
499 static void **per_group_data = NULL; /* array of pointers to tension_handler_data_t */
500 static double last_x;
501 tension_handler_data_t *tdata;
502 double F, *pStiffness;
506 fprintf(stderr, "%s:\tfinding tension at distance %g%s\n", __FUNCTION__, x, transition ? " after a transition" : "");
509 if (transition != 0 || num_active_groups == 0) { /* setup statics */
510 /* get new statics, freeing the old ones if they aren't NULL */
511 get_active_groups(state_list, &num_active_groups, &per_group_handlers, &per_group_inverse_handlers, &per_group_data);
513 /* fill in the group handlers and params */
515 fprintf(stderr, "%s:\t%d groups, tension handler array %p, inverse tension handler array %p, handler data array %p, (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_inverse_handlers, per_group_data, per_group_data[0]);
517 for (i=0; i<num_active_groups; i++) {
518 tdata = (tension_handler_data_t *) per_group_data[i];
520 fprintf(stderr, "%s:\tactive group %d of %d (data %p, param list %p)", __FUNCTION__, i, num_active_groups, tdata, tdata->group_tension_params);
521 list_print(stderr, tdata->group_tension_params, " parameters");
524 /* tdata->persist continues unchanged */
527 } /* else continue with the current statics */
532 pStiffness = &env->stiffness;
534 F = tension_balance(num_active_groups, per_group_handlers, per_group_inverse_handlers, (void **)per_group_data, last_x, x, pStiffness);
540 @ For the moment, we will restrict our group boundaries to a single
541 dimension, so $\sum_i x_i = x$, our total extension (it is this
542 restriction that causes the groups to interact linearly). We'll also
543 restrict our extensions to all be positive. With these restrictions,
544 the problem of balancing the tensions reduces to solving a system of
545 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
546 the number of active groups. In general this can be fairly
547 complicated, but without much loss of practicality we can restrict
548 ourselves to strictly monotonically increasing, non-negative tension
549 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
550 simpler. For example, it guarantees the existence of a unique, real
551 solution for finite forces. See Appendix \ref{sec.tension_balance}
552 for the tension balancing implementation.
554 Each group must define a parameter structure if the tension model is
556 a function to create the parameter structure,
557 a function to destroy the parameter structure, and
558 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as its data argument.
559 The tension-specific parameter structures are contained in the domain
560 group field. For optimization, a group may also define an inverse
561 tension function as an optimization. See the Section
562 \ref{sec.model_selection} for a discussion on the command line
563 framework and inverse function discussion. See the worm-like chain
564 implementation for an example (Section \ref{sec.wlc}).
566 The major limitation of our approach is that all of our tension models
567 are Markovian (memory-less), which is not really a problem for the
568 chains (see \citet{evans99} for WLC thermalization rate calculations).
570 \subsection{Transition rate}
571 \label{sec.transition_rate}
573 Each state transition is modeled as unimolecular, first order reaction
575 \text{State 1} \xrightarrow{k} \text{State 2}
577 With the rate constant $k$ defined by
579 \frac{d[\text{State 2}]}{dt} = -\frac{d[\text{State 1}]}{dt} = k [\text{State 1}].
581 So the probability of a given domain transitioning in a short time
587 For a group of $N$ identical domains, and therefore identical
588 unfolding probabilities $P_1$, the probability of $n$ domains
589 unfolding in a given timestep is given by the binomial distribution
591 P(n) = \frac{N!}{n!(N-n)!}P_1^n(1-P_1)^{(N-n)}.
593 The probability that \emph{none} of these domains unfold is then
597 so the probability that \emph{at least one} domain unfolds is
601 Note that for small $P$,
603 p(n>0) = 1-(1-NP_1+\mathcal{O}(P_1^2)) = NP_1 - \mathcal{O}(P^2)
606 This conversion from $P_1$ to $P_N$ is handled by the [[pN]] macro
608 /* find multi-domain transition rate from the single domain rate */
609 #define PN(P,N) (1.0-pow(1.0-(P), (N)))
613 We can also define a macro to find the approprate $dt$ to achieve a
614 target $P_N$ with a given $k$ and $N$.
616 P_N &= 1-(1-P_1)^N \\
617 1-P_1 &= (1-P_N)^{1/N} \\
618 P_1 &= 1-(1-P_N)^{1/N}
621 #define P1(PN,N) (1.0-pow(1.0-(PN), 1.0/(N)))
624 We take some time to discuss the meaning of $p(n>1)$
625 (i.e. multi-unfolding timesteps) in Section \ref{sec.timescales}.
627 <<does domain transition>>=
628 int domain_transitions(double F, double dt, environment_t *env, transition_t *transition)
629 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to transition structure */
631 k = accel_k(transition->k, F, env, transition->k_params);
632 //(*transition->k)(F, env, domain->k_params);
633 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
634 return happens(PN(k*dt, N_INIT(transition)));
636 @ [[happens]] is a random decision making function defined in Appendix \ref{sec.utils}.
638 <<randomly transition domains>>=
639 int random_transitions(list_t *transition_list, double tension,
640 double dt, environment_t *env)
641 { /* returns 1 if there was a transition and 0 otherwise */
642 while (transition_list != NULL) {
643 if (T(transition_list)->initial_state->num_domains > 0
644 && domain_transitions(tension, dt, env, T(transition_list))) {
645 if (flags & FLAG_OUTPUT_TRANSITION_FORCES)
646 fprintf(stdout, "%g\t%s\t%s\n", tension, T(transition_list)->initial_state->name, T(transition_list)->final_state->name);
647 T(transition_list)->initial_state->num_domains--;
648 T(transition_list)->final_state->num_domains++;
649 return 1; /* our one domain has transitioned, stop looking for others */
651 transition_list = transition_list->next;
657 \subsection{Timescales}
658 \label{sec.timescales}
660 The simulation assumes that chain equilibration is instantaneous
661 relative to the simulation timestep $dt$, so that tension is uniform
662 along the chain. The quality of this assumption depends on your
663 particular chain. For example, a damped spring thermalizes on a
664 timescale of order $\tau = 1/\gamma$ where the damping ratio $\gamma
665 \equiv \omega_0(-\zeta \pm \sqrt{\zeta^2-1}$, the natural angular
666 frequency $omega_0 \equiv \sqrt{k/m}$, $\zeta \equiv b/2\sqrt{km}$,
667 and $b$ sets the drag $F_b = -bv$. For our cantilevers $\tau$ is
668 on the order of milliseconds, which is longer than a timestep.
669 However, the approximation is still reasonable, because there is
670 rarely unfolding during the cantilever return. You could, of course,
671 take cantilever, WLC, etc.\ response times into effect, but that
672 would significantly complicate a simulation that seems to work
673 well enough without bothering :p. For WLC and FJC relaxation times,
674 see the Appendix of \citet{evans99} and \citet{kroy07}.
676 Besides assuming our timestep is much greater than equilibration
677 times, we also force it to be short enough so that only one domain may
678 unfold in a given timestep. For an unfolding timescale of $dt_u =
679 1/k$, we adapt our timesteps to keep $dt \ll dt_u$, so the probability
680 of two domains unfolding in the same timestep is small (see
681 [[determine_dt]] in Section \ref{sec.adaptive_dt} and [[P]] in
682 [[main()]] in Section \ref{sec.main} set by [[-P]] in
683 [[get_options()]] in Section \ref{sec.get_options}). This approach
684 breaks down as the adaptive timestep scheme approaches the transition
685 timestep $dt_t$, but $dt_t \sim 1\U{$\mu$s}$ for Ig27-like proteins
686 \citep{klimov00}, so this shouldn't be a problem. To reassure
687 yourself, you can ask the simulator to print the smallest timestep
688 that was used with TODO.
690 Even if the likelyhood of two domains unfolding in the same timestep
691 is small, if you run long enough simulations it will eventually occur.
692 In this case, we assume that $dt_t \ll dt$, so even if two domains
693 would unfold if we stuck to our unfolding rate $k$ for an entire
694 timestep $dt$, in ``reality'' only one of those domains unfolds.
695 Because we do not know when in the timestep the unfolding took place,
696 we move on to the next timestep without checking for further
697 unfoldings. This ``unchecked timestep portion'' should not contribute
698 significant errors to the calculation, because if $dt$ was small
699 enough that unfolding was unlikely at the pre-unfolding force, it
700 should be even more unlikely at the post-unfolding force, especially
701 over only a fraction of the timestep. The only dangerous cases for
702 the current implementation are those where unfolding intermediates are
703 much less stable than their precursor states, in which case an
704 unfolding event during the remaining timestep may actually be likely.
705 A simple workaround for such cases is to pick the value for [[P]] to
706 be small enough that the $dt \ll dt_u$ assumption survives even under
707 a transition populating the unstable intermediate.
709 \subsection{Adaptive timesteps}
710 \label{sec.adaptive_dt}
712 We'd like to pick $dt$ so the probability of changing state
713 (e.g.\ unfolding) in the next timestep is small. If we don't adapt our
714 timestep, we also risk breaking our approximation that $F(x' \in [x,
715 x+v \cdot dt]) = F(x)$ or that only one domain changes state in a
716 given timestep. The simulation would have been accurate for
717 sufficiently small fixed timesteps, but adaptive timesteps will allow
718 us to move through low tension regions in fewer steps, leading to a
719 more efficient simulation.
721 Consider the two state folded $\rightarrow$ unfolded transition.
722 Because $F(x)$ is monotonically increasing between unfoldings,
723 excessively large timesteps will lead to erroneously small unfolding
724 rates (an thus, higher average unfolding force).
726 The actual adaptive timestep implementation is not particularly
727 interesting, since we are only required to reduce $dt$ to somewhere
728 below a set threshold, so I've removed it to Appendix
729 \ref{sec.adaptive_dt_implementation}.
735 The models provide the physics of the simulation, but the simulation
736 allows interchangeable models, and we are currently focusing on the
737 simulation itself, so we remove the actual model implementation to the
740 The tension models are defined in Appendix \ref{sec.tension_models}
741 and unfolding models are defined in Appendix \ref{sec.k_models}.
744 #define k_B 1.3806503e-23 /* J/K */
748 \section{Command line interface}
750 <<get option functions>>=
752 <<index tension model>>
757 \subsection{Model selection}
758 \label{sec.model_selection}
760 The main difficulty with the command line interface in \prog\ is
761 developing an intuitive method for accessing the possibly large number
762 of available models. We'll treat this generally by defining an array
763 of available models, containing their important parameters.
765 <<get options definitions>>=
766 <<model getopt definitions>>
769 <<create/destroy definitions>>=
770 typedef void *create_data_func_t(char **param_strings);
771 typedef void destroy_data_func_t(void *param_struct);
774 <<model getopt definitions>>=
775 <<tension model getopt definitions>>
776 <<k model getopt definitions>>
779 In order to choose models, we need a method of assembling a reaction
780 graph such as that in Figure \ref{fig.domain_states} and population
781 the graph with a set of domains. First we declare the domain states
784 [[-s <state-name>,<model>[,<group>] -S <model-params> [-N <initial-population]]] \\
788 [[-s unfolded,wlc,1 -S 1e-8,4e-10 -N 5]]
790 This creates the state named [[unfolded]], using the WLC model and one
791 as a group number (defaults to 0, see Section \ref{sec.find_tension}).
792 The tension parameters are then set to [[1e-8,4e-10]], which in the
793 case of the WLC are the contour and persistence lengths respectively.
794 Finally we populate the state by adding five domains to the state.
795 The name of the state is importent for identifying states later.
797 Once the domains have been created and populated, you can added
798 transition rates following
800 [[-k <initial-state>,<final-state>,<model> -K <model-params>]] \\
804 [[-k folded,unfolded,bell -K 3.3e-4,0.25e-9]]
806 This creates a reaction going from the [[folded]] state to the
807 [[unfolded]] state, with the rate constant given by the Bell model
808 (Appendix \ref{sec.bell}). The unfolding parameters are then set to
809 [[3.3e-4,0.25e-9]], which in the case of the Bell model are the
810 unforced rate and transition state distance respectively.
812 \subsubsection{Tension}
814 To access the various tension models from the command line, we define
815 a structure that contains all the neccessary information about the
817 <<tension model getopt definitions>>=
818 typedef struct tension_model_getopt_struct {
819 one_dim_fn_t *handler; /* fn implementing the model on a group */
820 one_dim_fn_t *inverse_handler; /* fn implementing inverse on group */
821 char *name; /* name string identifying the model */
822 char *description; /* a brief description of the model */
823 int num_params; /* number of extra params for the model */
824 char **param_descriptions; /* descriptions of extra params */
825 char *params; /* default values for the extra params */
826 create_data_func_t *creator; /* param-string -> model-data-struct fn */
827 destroy_data_func_t *destructor; /* fn to free the model data structure */
828 } tension_model_getopt_t;
829 @ Set [[inverse_handler]] to [[NULL]] if you wish to invert using
830 bisection. Obviously, this will be slower than computing an
831 analytical inverse and more prone to rounding errors, but it may be
832 easier if a closed-form inverse does not exist.
834 <<tension model getopt array definition>>=
835 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
836 <<null tension model getopt>>,
837 <<constant tension model getopt>>,
838 <<hooke tension model getopt>>,
839 <<worm-like chain tension model getopt>>,
840 <<freely-jointed chain tension model getopt>>,
841 <<piston tension model getopt>>
845 \subsubsection{Unfolding rate}
847 <<k model getopt definitions>>=
848 #define NUM_K_MODELS 6
850 typedef struct k_model_getopt_struct {
855 char **param_descriptions;
857 create_data_func_t *creator;
858 destroy_data_func_t *destructor;
862 <<k model getopt array definition>>=
863 k_model_getopt_t k_models[NUM_K_MODELS] = {
864 <<null k model getopt>>,
865 <<const k model getopt>>,
866 <<bell k model getopt>>,
867 <<kbell k model getopt>>,
868 <<kramers k model getopt>>,
869 <<kramers integ k model getopt>>
876 void help(char *prog_name, double P, double dF,
877 double t_max, double dt_max, double v, double x_step,
878 state_t *stop_state, environment_t *env,
879 int n_tension_models, tension_model_getopt_t *tension_models,
880 int n_k_models, k_model_getopt_t *k_models,
881 int k_model, list_t *state_list)
886 if (state_list != NULL)
887 state = S(tail(state_list));
889 printf("usage: %s [options]\n", prog_name);
890 printf("Version %s\n\n", VERSION);
891 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
893 printf("Simulation options:\n");
895 printf("-q state\tStop simulation when this state empties (- for no limit, currently %s)\n", stop_state ? stop_state->name : "-" );
896 printf("-t tmax\tMaximum simulated time (-1 for no limit, currently %g s)\n", t_max);
897 printf("-d dt\tMaximum allowed timestep dt (currently %g s)\n", dt_max);
898 printf("-P P\tMaximum transition probability for dt (currently %g)\n", P);
899 printf("-F dF\tMaximum change in force over dt. Only relevant if dx > 0. (currently %g N)\n", dF);
900 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
901 printf("-x dx\tPulling step size (currently %g m)\n", x_step);
902 printf("\tWhen dx = 0, the pulling is continuous.\n");
903 printf("\tWhen dx > 0, the pulling occurs in discrete steps.\n");
905 printf("Environmental options:\n");
906 printf("-T T\tTemperature T (currently %g K)\n", env->T);
907 printf("-C T\tYou can also set the temperature T in Celsius\n");
909 printf("State creation options:\n");
910 printf("-s name,model[[,group],params]\tCreate a new state.\n");
911 printf("Once you've set up the state, you may populate it with domains\n");
912 printf("-N n\tAdd n domains in the current state (%s)\n", state ? state->name : "-");
914 printf("Transition creation options:\n");
915 printf("-k initial_state,final_state,model[,params]\tCreate a new transition\n");
917 printf("To double check your reaction graph, you can produce graphviz dot output\n");
918 printf("describing the current states and transitions.\n");
919 printf("-D\tPrint dot description of states and transitions and exit.\n");
921 printf("Simulation output mode options:\n");
922 printf("There are two output modes. In standard mode, only the transition\n");
923 printf("events are printed. For example:\n");
924 printf(" #Force (N)\tInitial state\tFinal state\n");
925 printf(" 123.456e-12\tfolded\tunfolded\n");
927 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
928 printf(" #Distance (m)\tForce (N)\tStiffness (N/m)\n");
929 printf(" 0.001\t0.0023\t0.002\n");
931 printf("-V\tChange output to verbose mode.\n");
932 printf("-h\tPrint this help and exit.\n");
935 printf("Tension models:\n");
936 for (i=0; i<n_tension_models; i++) {
937 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
938 for (j=0; j < tension_models[i].num_params ; j++)
939 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
940 printf("\t\tdefaults: %s\n", tension_models[i].params);
942 printf("Unfolding rate models:\n");
943 for (i=0; i<n_k_models; i++) {
944 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
945 for (j=0; j < k_models[i].num_params ; j++)
946 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
947 printf("\t\tdefaults: %s\n", k_models[i].params);
950 printf("Examples:\n");
951 printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded. Stop once there are no more folded domains.\n");
952 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);
953 printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded. Stop after 0.25 seconds\n");
954 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);
958 \subsection{Get options}
959 \label{sec.get_options}
963 void get_options(int argc, char **argv, double *pP, double *pDF,
964 double *pT_max, double *pDt_max, double *pV,
965 double *pX_step, state_t **pStop_state, environment_t *env,
966 int n_tension_models, tension_model_getopt_t *tension_models,
967 int n_k_models, k_model_getopt_t *k_models,
968 list_t **pState_list, list_t **pTransition_list)
970 char *prog_name = NULL;
971 char c, options[] = "q:t:d:P:F:v:x:T:C:s:N:k:DVh";
972 int tension_model=0, k_model=0, initial_state=-1, final_state=-1;
973 char *state_name=NULL;
974 int i, n, tension_group=0;
975 list_t *temp_list=NULL;
978 void *model_params; /* for INIT_MODEL() */
979 int num_param_args; /* for INIT_MODEL() */
980 char **param_args; /* for INIT_MODEL() */
982 extern int optind, optopt, opterr;
987 for (i=0; i<n_tension_models; i++) {
988 fprintf(stderr, "tension model %i %s:\t model's handler %p\t inverse handler %p\n",
989 i, tension_models[i].name, tension_models[i].handler, tension_models[i].inverse_handler);
990 assert(tension_models[i].handler == tension_handlers[i]);
995 flags = FLAG_OUTPUT_TRANSITION_FORCES;
998 *pT_max = -1; /* s, -1 removes this limit */
999 *pDt_max = 0.001; /* s */
1000 *pP = 1e-3; /* % pop per s */
1001 *pDF = 1e-12; /* N */
1002 *pV = 1e-6; /* m/s */
1003 *pX_step = 0.0; /* m */
1004 env->T = 300.0; /* K */
1005 *pState_list = NULL;
1006 *pTransition_list = NULL;
1008 while ((c=getopt(argc, argv, options)) != -1) {
1011 if (STRMATCH(optarg, "-")) {
1012 *pStop_state = NULL;
1014 temp_list = head(*pState_list);
1015 while (temp_list != NULL) {
1016 if (STRMATCH(S(temp_list)->name, optarg)) {
1017 *pStop_state = S(temp_list);
1020 temp_list = temp_list->next;
1024 case 't': *pT_max = safe_strtod(optarg, "-t"); break;
1025 case 'd': *pDt_max = safe_strtod(optarg, "-d"); break;
1026 case 'P': *pP = safe_strtod(optarg, "-P"); break;
1027 case 'F': *pDF = safe_strtod(optarg, "-F"); break;
1028 case 'v': *pV = safe_strtod(optarg, "-v"); break;
1029 case 'x': *pX_step = safe_strtod(optarg, "-x"); break;
1030 case 'T': env->T = safe_strtod(optarg, "-T"); break;
1031 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
1033 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
1034 assert(num_strings >= 2 && num_strings <= 4);
1036 state_name = strings[0];
1037 if (state_index(*pState_list, state_name) != -1) {
1038 fprintf(stderr, "Invalid option '%s', repeated state name '%s'\n", optarg, state_name);
1041 tension_model = index_tension_model(n_tension_models, tension_models, strings[1]);
1042 if (num_strings == 4)
1043 tension_group = safe_strtoi(strings[2], optarg);
1046 if (num_strings >= 3)
1047 tension_models[tension_model].params = strings[num_strings-1];
1049 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);
1051 INIT_MODEL(state_name, &tension_models[tension_model], tension_models[tension_model].params, model_params);
1052 push(pState_list, create_state(state_name,
1053 tension_models[tension_model].handler,
1054 tension_models[tension_model].inverse_handler,
1057 tension_models[tension_model].destructor,
1059 *pState_list = tail(*pState_list);
1061 free_parsed_list(num_strings, strings);
1064 n = safe_strtoi(optarg, "-N");
1066 fprintf(stderr, "%s:\tadding %d domains to %s\n", __FUNCTION__, n, S(*pState_list)->name);
1068 S(*pState_list)->num_domains += n;
1071 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
1072 assert(num_strings >= 3 && num_strings <= 4);
1073 initial_state = state_index(*pState_list, strings[0]);
1074 final_state = state_index(*pState_list, strings[1]);
1075 k_model = index_k_model(n_k_models, k_models, strings[2]);
1077 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);
1079 assert(initial_state != final_state);
1080 if (num_strings == 4)
1081 k_models[k_model].params = strings[3];
1082 INIT_MODEL(state_name, &k_models[k_model], k_models[k_model].params, model_params);
1083 push(pTransition_list, create_transition(list_element(*pState_list, initial_state),
1084 list_element(*pState_list, final_state),
1085 k_models[k_model].k,
1087 k_models[k_model].destructor), 1);
1088 free_parsed_list(num_strings, strings);
1091 print_state_graph(stdout, *pState_list, *pTransition_list);
1095 flags = FLAG_OUTPUT_FULL_CURVE;
1098 fprintf(stderr, "unrecognized option '%c'\n", optopt);
1099 /* fall through to default case */
1101 help(prog_name, *pP, *pDF, *pT_max, *pDt_max, *pV, *pX_step,
1102 *pStop_state, env, n_tension_models, tension_models,
1103 n_k_models, k_models, k_model, *pState_list);
1107 /* check the arguments */
1108 assert(*pState_list != NULL); /* no domains! */
1109 assert(*pP > 0.0); assert(*pP < 1.0);
1110 assert(*pDt_max > 0.0);
1112 assert(env->T > 0.0);
1114 crosslink_groups(*pState_list, *pTransition_list);
1120 <<index tension model>>=
1121 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
1124 for (i=0; i<n_models; i++)
1125 if (STRMATCH(models[i].name, name))
1127 if (i == n_models) {
1128 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
1136 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
1139 for (i=0; i<n_models; i++)
1140 if (STRMATCH(models[i].name, name))
1142 if (i == n_models) {
1143 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
1150 <<initialize model definition>>=
1151 /* requires int num_param_args and char **param_args in the current scope
1153 * INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
1154 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
1156 #define INIT_MODEL(role, model, param_string, param_pointer) \
1158 parse_list_string((param_string), SEP, '{', '}', \
1159 &num_param_args, ¶m_args); \
1160 if (num_param_args != (model)->num_params) { \
1162 "%s model %s expected %d params,\n", \
1163 (role), (model)->name, (model)->num_params); \
1165 "not the %d params in '%s'\n", \
1166 num_param_args, (param_string)); \
1167 assert(num_param_args == (model)->num_params); \
1169 if ((model)->creator) \
1170 param_pointer = (*(model)->creator)(param_args); \
1172 param_pointer = NULL; \
1173 free_parsed_list(num_param_args, param_args); \
1177 First we define some safe string parsing functions. Currently these
1178 abort on any irregularity, but in the future they could print error
1179 messages, etc. [[static]] because the functions are currently
1180 defined in each [[*.c]] file that uses them, but perhaps they should
1181 be moved to [[utils.h/c]] or some such instead\ldots
1184 static int safe_strtoi(const char *s, const char *description) {
1188 ret = strtol(s, &endp, 10);
1189 if (endp[0] != '\0') { //strlen(endp) > 0) {
1190 fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1191 endp, s, description);
1192 assert(1==0); //strlen(endp) == 0);
1194 fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1197 } else if (errno == ERANGE) {
1198 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1199 assert(errno != ERANGE);
1204 static double safe_strtod(const char *s, const char *description) {
1208 ret = strtod(s, &endp);
1209 if (endp[0] != '\0') { //strlen(endp) > 0) {
1210 fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1211 endp, s, description);
1212 assert(1==0); //strlen(endp) == 0);
1214 fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1217 } else if (errno == ERANGE) {
1218 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1219 assert(errno != ERANGE);
1228 \addcontentsline{toc}{section}{Appendicies}
1230 \section{sawsim.c details}
1234 The general layout of our simulation code is:
1245 We include [[math.h]], so don't forget to link to the libm with `-lm'.
1247 #include <assert.h> /* assert() */
1248 #include <stdlib.h> /* malloc(), free(), rand(), strto*() */
1249 #include <stdio.h> /* fprintf(), stdout */
1250 #include <string.h> /* strlen, strtok() */
1251 #include <math.h> /* exp(), M_PI, sqrt() */
1252 #include <sys/types.h> /* pid_t (returned by getpid()) */
1253 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
1254 #include <errno.h> /* errno, ERANGE (for safe_strto*()) */
1257 #include "tension_balance.h"
1258 #include "k_model.h"
1259 #include "tension_model.h"
1261 #include "accel_k.h"
1265 <<version definition>>
1266 <<flag definitions>>
1267 <<max/min definitions>>
1268 <<string comparison definition and globals>>
1269 <<initialize model definition>>
1270 <<get options definitions>>
1271 <<state definitions>>
1272 <<transition definitions>>
1281 <<create/destroy state>>
1282 <<destroy state list>>
1284 <<create/destroy transition>>
1285 <<destroy transition list>>
1286 <<print state graph>>
1288 <<simulation functions>>
1289 <<boilerplate functions>>
1292 <<boilerplate functions>>=
1294 <<get option functions>>
1300 srand(getpid()*time(NULL)); /* seed rand() */
1304 <<flag definitions>>=
1305 /* in octal b/c of prefixed '0' */
1306 #define FLAG_OUTPUT_FULL_CURVE 01
1307 #define FLAG_OUTPUT_TRANSITION_FORCES 02
1311 static unsigned long int flags = 0;
1314 \subsection{Utilities}
1317 <<max/min definitions>>=
1318 #define MAX(a,b) ((a)>(b) ? (a) : (b))
1319 #define MIN(a,b) ((a)<(b) ? (a) : (b))
1322 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
1323 <<string comparison definition and globals>>=
1324 // Check if two strings match, return 1 if they do
1325 static char *temp_string_A;
1326 static char *temp_string_B;
1327 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
1328 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
1329 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
1330 /* +1 to also compare the '\0' */
1333 We also define a macro for our [[check]] unit testing
1334 <<check relative error>>=
1335 #define CHECK_ERR(max_err, expected, received) \
1337 fail_unless( (received-expected)/expected < max_err, \
1338 "relative error %g >= %g in %s (Expected %g, received %g)", \
1339 (received-expected)/expected, max_err, #received, \
1340 expected, received); \
1341 fail_unless(-(received-expected)/expected < max_err, \
1342 "relative error %g >= %g in %s (Expected %g, received %g)", \
1343 -(received-expected)/expected, max_err, #received, \
1344 expected, received); \
1349 int happens(double probability)
1351 assert(probability >= 0.0); assert(probability <= 1.0);
1352 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*/
1357 \subsection{Adaptive timesteps}
1358 \label{sec.adaptive_dt_implementation}
1360 $F(x)$ increases with $x$, possibly exploding (e.g.\ in the worm-like
1361 chain model), so basing the timestep on the unfolding probability at
1362 the current tension is dangerous, and we need to search for a $dt$ for
1363 which $P_N(F(x+v*dt))<P_\text{max}$ (See Section
1364 \ref{sec.transition_rate} for a discussion of $P_N$). There are two
1365 cases to consider. In the most common, no domains have unfolded since
1366 the last step, and we expect the next step to be slightly shorter than
1367 the previous one. In the less common, domains did unfold in the last
1368 step, and we expect the next step to be considerably longer than the
1371 double search_dt(transition_t *transition,
1373 environment_t *env, double x,
1374 double max_prob, double max_dt, double v)
1375 { /* Returns a valid timestep dt in seconds for the given transition.
1376 * Takes a list of all states, a pointer env to the environmental data,
1377 * a starting separation x in m, a max_prob between 0 and 1,
1378 * max_dt in s, stretching velocity v in m/s.
1380 double F, k, P, dtCur, dtU, dtUCur, dtL, dt;
1382 /* get upper bound using the current position */
1383 F = find_tension(state_list, env, x, 0, 1); /* BUG. repeated calculation */
1384 //printf("Start with x = %g (F = %g)\n", x, F);
1385 k = accel_k(transition->k, F, env, transition->k_params);
1386 //printf("x %g\tF %g\tk %g\n", x, F, k);
1387 dtU = P1(max_prob, N_INIT(transition)) / k; /* upper bound on dt */
1389 //printf("overshot max_dt\n");
1392 if (v == 0) /* discrete stepping, so force is temporarily constant */
1395 /* set a lower bound on dt too */
1398 /* The dt determined above may produce illegitimate forces or ks.
1399 * Reduce the upper bound until we have valid ks. */
1401 F = find_tension(state_list, env, x+v*dt, 0, 1);
1402 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1405 F = find_tension(state_list, env, x+v*dt, 0, 1);
1407 //printf("Try for dt = %g (F = %g)\n", dt, F);
1408 k = accel_k(transition->k, F, env, transition->k_params);
1409 /* returned k may be -1.0 */
1410 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1411 while (k == -1.0) { /* reduce step until we hit a valid k */
1413 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1414 F = find_tension(state_list, env, x+v*dt, 0, 1);
1415 //printf("Try for dt = %g (F = %g)\n", dt, F);
1416 k = accel_k(transition->k, F, env, transition->k_params);
1417 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1419 assert(dtU > 1e-14); /* timestep to valid k too small */
1420 /* safe timestep back from x+dtU */
1421 dtUCur = P1(max_prob, N_INIT(transition)) / k;
1423 return dt; /* dtU is safe. */
1425 /* dtU wasn't safe, lets see what would be. */
1426 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1427 dt = (dtU + dtL) / 2.0;
1428 F = find_tension(state_list, env, x+v*dt, 0, 1);
1429 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1430 k = accel_k(transition->k, F, env, transition->k_params);
1431 dtCur = P1(max_prob, N_INIT(transition)) / k;
1432 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1433 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1435 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1436 dtU = dt; /* ... stepping out only dtCur would be safe */
1439 break; /* dtCur = dt */
1441 return MAX(dtUCur, dtL);
1446 Determine $dt$ for an array of potentially different folded domains.
1447 We apply [[search_dt()]] to each populated domain to find the maximum
1448 $dt$ the domain can handle. The final $dt$ is less than the
1449 individual domain $dt$s (to satisfy [[search_dt()]], see above). If
1450 $v>0$ (continuous pulling), $dt$ is also reduced to satisfy
1451 $|F(x+v*dt)-F(x)|<dF_\text{max}$, which limits errors due to our
1452 transition rate assumption that the force does not change appreciably
1453 over a single timestep.
1456 double determine_dt(list_t *state_list, list_t *transition_list,
1457 environment_t *env, double F, double x, double max_dF,
1458 double max_prob, double dt_max, double v, int transition)
1459 { /* Returns the timestep dt in seconds.
1460 * Takes the state and transition lists, a pointer to the environment,
1461 * the force F(x) in N, an extension x in m, max_dF in N,
1462 * max_prob between 0 and 1, dt_max in s, v in m/s, and transition in {0,1}*/
1463 double dt=dt_max, new_dt, F_new;
1464 assert(max_dF > 0.0);
1465 assert(max_prob > 0.0); assert(max_prob < 1.0);
1466 assert(dt_max > 0.0);
1467 assert(state_list != NULL);
1468 assert(transition_list != NULL);
1469 transition_list = head(transition_list);
1472 /* Ensure |dF| < max_dF */
1473 F_new = find_tension(state_list, env, x+v*dt, transition, 1);
1474 while (fabs(F_new - F) > max_dF) {
1476 F_new = find_tension(state_list, env, x+v*dt, transition, 1);
1480 /* Ensure all individual domains pass search_dt() */
1481 while (transition_list != NULL) {
1482 if (T(transition_list)->initial_state->num_domains > 0) {
1483 new_dt = search_dt(T(transition_list),state_list,env,x,max_prob,dt,v);
1484 dt = MIN(dt, new_dt);
1486 transition_list = transition_list->next;
1493 \subsection{State data}
1494 \label{sec.state_data}
1496 The domains exist in one of an array of possible states. Each state
1497 has a tension model which determines it's force/extension curve
1498 (possibly [[null]] for rigid states, see Appendix
1499 \ref{sec.null_tension}).
1501 Groups are, for example, for two domain states with WLC tensions; one
1502 with $p=0.4\U{nm}$, $L=10\U{nm}$ and the other with $p=0.4\U{nm}$,
1503 $L=28\U{nm}$. Since the persistence lengths are the same, you would
1504 like to add the contour lengths of all the domains in \emph{both}
1505 states, and plug that total contour length into a single WLC
1506 evaluation (see Section \ref{sec.find_tension}).
1507 <<state definitions>>=
1508 typedef struct state_struct {
1509 char *name; /* name string */
1510 one_dim_fn_t *tension_handler; /* tension handler */
1511 one_dim_fn_t *inverse_tension_handler; /* inverse tension handler */
1512 int tension_group; /* for grouping similar tension states */
1513 void *tension_params; /* pointer to tension parameters */
1514 destroy_data_func_t *destroy_tension;
1515 int num_domains; /* number of domains currently in state */
1516 /* cached lookups generated from the above data */
1517 int num_out_trans; /* length of out_trans array */
1518 struct transition_struct **out_trans; /* transitions leaving this state */
1519 int num_grouped_states; /* length of grouped states array */
1520 struct state_struct **grouped_states; /* states grouped with this one (doesn't include self) */
1523 /* get the state data for the current list node */
1524 #define S(list) ((state_t *)(list)->d)
1526 @ [[tension_params]] is a pointer to the parameters used by the
1527 handler function when determining the tension. [[destroy_tension]]
1528 points to a function for freeing [[tension_params]]. We it in the
1529 state structure so that [[destroy_state]] and will have an easier job
1532 [[create_]] and [[destroy_state]] are simple wrappers around
1533 [[malloc]] and [[free]].
1534 <<create/destroy state>>=
1535 state_t *create_state(char *name,
1536 one_dim_fn_t *tension_handler,
1537 one_dim_fn_t *inverse_tension_handler,
1539 void *tension_params,
1540 destroy_data_func_t *destroy_tension,
1543 state_t *ret = (state_t *)malloc(sizeof(state_t));
1544 assert(ret != NULL);
1545 /* make a copy of the name, so we aren't dependent on the original */
1546 ret->name = (char *)malloc(sizeof(char)*(strlen(name)+1));
1547 assert(ret->name != NULL);
1548 strcpy(ret->name, name); /* we know ret->name is long enough */
1550 ret->tension_handler = tension_handler;
1551 ret->inverse_tension_handler = inverse_tension_handler;
1552 ret->tension_group = tension_group;
1553 ret->tension_params = tension_params;
1554 ret->destroy_tension = destroy_tension;
1555 ret->num_domains = num_domains;
1557 ret->num_out_trans = 0;
1558 ret->out_trans = NULL;
1559 ret->num_grouped_states = 0;
1560 ret->grouped_states = NULL;
1563 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);
1568 void destroy_state(state_t *state)
1572 fprintf(stderr, "freeing state %s (%p)\n", state->name, state);
1576 if (state->destroy_tension)
1577 (*state->destroy_tension)(state->tension_params);
1578 if (state->out_trans)
1579 free(state->out_trans);
1580 if (state->grouped_states)
1581 free(state->grouped_states);
1588 We'll be storing the states in a list (see Appendix \ref{sec.lists}),
1589 so we also define a convenience function for cleaning up.
1590 <<destroy state list>>=
1591 void destroy_state_list(list_t *state_list)
1593 state_list = head(state_list);
1594 while (state_list != NULL)
1595 destroy_state((state_t *) pop(&state_list));
1600 We can also define a convenient way to get a state index from it's name.
1602 int state_index(list_t *state_list, char *name)
1605 state_list = head(state_list);
1606 while (state_list != NULL) {
1607 if (STRMATCH(S(state_list)->name, name)) return i;
1608 state_list = state_list->next;
1616 \subsection{Transition data}
1617 \label{sec.transition_data}
1619 Once you've created states (Sections \ref{sec.main},
1620 \ref{sec.model_selection}, and \ref{sec.state_data}), you need to
1621 create transitions between them.
1622 <<transition definitions>>=
1623 typedef struct transition_struct {
1624 state_t *initial_state;
1625 state_t *final_state;
1626 k_func_t *k; /* transition rate model */
1627 void *k_params; /* pointer to k parameters */
1628 destroy_data_func_t *destroy_k;
1631 /* get the transition data for the current list node */
1632 #define T(list) ((transition_t *)(list)->d)
1634 /* get the number of initial state population for a transition state */
1635 #define N_INIT(transition) ((transition)->initial_state->num_domains)
1639 @ [[k]] is a pointer to the function determining the transition rate
1640 for a given tension. [[k_params]] and [[destroy_k]] have similar
1641 roles with regards to [[k]] as [[tension_params]] and
1642 [[destroy_tension]] have with regards to [[tension_handler]] (see
1643 Section \ref{sec.state_data}).
1645 [[create_]] and [[destroy_transition]] are simple wrappers around
1646 [[malloc]] and [[free]].
1647 <<create/destroy transition>>=
1648 transition_t *create_transition(state_t *initial_state, state_t *final_state,
1651 destroy_data_func_t *destroy_k)
1653 transition_t *ret = (transition_t *)malloc(sizeof(transition_t));
1654 assert(ret != NULL);
1655 assert(initial_state != NULL);
1656 assert(final_state != NULL);
1657 ret->initial_state = initial_state;
1658 ret->final_state = final_state;
1660 ret->k_params = k_params;
1661 ret->destroy_k = destroy_k;
1663 fprintf(stderr, "generated transition %p from %s to %s, params %p\n", ret, ret->initial_state->name, ret->final_state->name, ret->k_params);
1668 void destroy_transition(transition_t *transition)
1672 fprintf(stderr, "freeing transition %p (%s to %s)\n", transition, transition->initial_state->name, transition->final_state->name);
1674 if (transition->destroy_k)
1675 (*transition->destroy_k)(transition->k_params);
1682 We'll be storing the transitions in a list (see Appendix
1683 \ref{sec.lists}), so we also define a convenience function for
1685 <<destroy transition list>>=
1686 void destroy_transition_list(list_t *transition_list)
1688 transition_list = head(transition_list);
1689 while (transition_list != NULL)
1690 destroy_transition((transition_t *) pop(&transition_list));
1695 \subsection{Graphviz output}
1696 \label{sec.graphviz_output}
1698 <<print state graph>>=
1699 void print_state_graph(FILE *file, list_t *state_list, list_t *transition_list)
1701 state_list = head(state_list);
1702 transition_list = head(transition_list);
1703 fprintf(file, "digraph states {\n node [shape=record,color=black,fillcolor=white,style=filled];\n\n");
1705 fprintf(file, " n%d [label=\"%s\"]\n", state_index(state_list, S(state_list)->name), S(state_list)->name);
1706 if (state_list->next == NULL) break;
1707 state_list = state_list->next;
1709 fprintf(file, "\n");
1710 while (transition_list != NULL) {
1711 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));
1712 transition_list = transition_list->next;
1714 fprintf(file, "}\n");
1718 \subsection{Domain model and group handling}
1720 <<group functions>>=
1721 <<crosslink groups>>
1722 <<get active groups>>
1725 Scan through a list of states and construct an array of all the
1727 <<get active groups>>=
1728 void get_active_groups(list_t *state_list,
1729 int *pNum_active_groups,
1730 one_dim_fn_t ***pPer_group_handlers,
1731 one_dim_fn_t ***pPer_group_inverse_handlers,
1732 void ***pPer_group_data /* pointer to array of pointers to tension_handler_data_t */)
1734 void **active_groups=NULL;
1735 one_dim_fn_t *handler, *old_handler, *inverse_handler, *old_inverse_handler;
1736 one_dim_fn_t **per_group_handlers=NULL, **per_group_inverse_handlers=NULL;
1737 void **per_group_data=NULL; /* array of pointers to tension_handler_data_t */
1738 int i, j, num_domains, group, old_group, num_active_groups=0;
1739 list_t *active_states_list=NULL;
1740 tension_handler_data_t *tdata=NULL;
1743 /* get all the active groups in a list */
1744 state_list = head(state_list);
1746 fprintf(stderr, "%s:\t", __FUNCTION__);
1747 list_print(stderr, state_list, "states");
1749 while (state_list != NULL) {
1750 state = S(state_list);
1751 handler = state->tension_handler;
1752 inverse_handler = state->inverse_tension_handler;
1753 group = state->tension_group;
1754 num_domains = state->num_domains;
1755 if (list_index(active_states_list, state) == -1) {
1756 /* we haven't added this state (or it's associated group) yet */
1757 if (num_domains > 0 && handler != NULL) { /* add the state */
1758 num_active_groups++;
1759 push(&active_states_list, state, 1);
1761 fprintf(stderr, "%s:\tactive %s state (%p) with %d domain(s)\n", __FUNCTION__, state->name, state, state->num_domains);
1763 for (i=0; i < state->num_grouped_states; i++) {
1764 if (state->grouped_states[i]->num_domains > 0) {
1765 /* add this related (and active) state now */
1766 assert(state->grouped_states[i]->tension_handler == handler);
1767 assert(state->grouped_states[i]->tension_group == group);
1768 push(&active_states_list, state->grouped_states[i], 1);
1770 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);
1774 } /* else empty state or NULL handler */
1775 } /* else state already added */
1776 state_list = state_list->next;
1779 fprintf(stderr, "%s:\t", __FUNCTION__);
1780 list_print(stderr, active_states_list, "active states");
1783 assert(num_active_groups <= list_length(active_states_list));
1784 assert(num_active_groups > 0);
1786 /* allocate the arrays we'll be returning */
1787 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1788 assert(per_group_handlers != NULL);
1789 per_group_inverse_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1790 assert(per_group_inverse_handlers != NULL);
1791 per_group_data = (void **)malloc(sizeof(void *)*num_active_groups);
1792 assert(per_group_data != NULL);
1794 fprintf(stderr, "%s:\tallocating per group data %p:", __FUNCTION__, per_group_data);
1796 for (i=0; i<num_active_groups; i++) { /* we might want to recycle any group data structs allocated in the previous round */
1797 per_group_data[i] = malloc(sizeof(tension_handler_data_t));
1798 assert(per_group_data[i] != NULL);
1800 fprintf(stderr, " %d,%p", i, per_group_data[i]);
1804 fprintf(stderr, "\n");
1807 /* populate the arrays we'll be returning */
1808 active_states_list = head(active_states_list);
1810 old_inverse_handler = NULL;
1811 j = -1; /* j is the current active _group_ index */
1812 while (active_states_list != NULL) {
1813 state = (state_t *) pop(&active_states_list);
1814 handler = state->tension_handler;
1815 inverse_handler = state->inverse_tension_handler;
1816 group = state->tension_group;
1817 num_domains = state->num_domains;
1818 if (handler != old_handler || inverse_handler != old_inverse_handler || group != old_group) {
1819 j++; /* move to the next group */
1820 tdata = (tension_handler_data_t *) per_group_data[j];
1821 per_group_handlers[j] = handler;
1822 per_group_inverse_handlers[j] = inverse_handler;
1823 tdata->group_tension_params = NULL;
1825 tdata->persist = NULL; /* we might want to grab any appropriate data from the old persist here... */
1828 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);
1830 for (i=0; i<num_domains; i++) { /* add num_domains copies of the tension params */
1831 push(&tdata->group_tension_params, state->tension_params, 1);
1834 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);
1836 old_handler = handler;
1837 old_inverse_handler = inverse_handler;
1841 /* free old groups */
1842 if (*pPer_group_handlers != NULL)
1843 free(*pPer_group_handlers);
1844 if (*pPer_group_inverse_handlers != NULL)
1845 free(*pPer_group_inverse_handlers);
1846 if (*pPer_group_data != NULL) {
1847 for (i=0; i<*pNum_active_groups; i++)
1848 free((*pPer_group_data)[i]);
1849 free(*pPer_group_data);
1852 /* set new groups */
1854 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]);
1856 *pNum_active_groups = num_active_groups;
1857 *pPer_group_handlers = per_group_handlers;
1858 *pPer_group_inverse_handlers = per_group_inverse_handlers;
1859 *pPer_group_data = per_group_data;
1863 <<crosslink groups>>=
1864 void crosslink_groups(list_t *state_groups, list_t *transition_list)
1866 list_t *list, *out_trans=NULL, *associates=NULL;
1867 one_dim_fn_t *handler;
1870 state_groups = head(state_groups);
1871 while (state_groups != NULL) {
1872 /* find transitions leaving this state */
1874 fprintf(stderr, "finding transitions leaving %s\n", S(state_groups)->name);
1876 list = head(transition_list);
1877 while (list != NULL) {
1878 if (T(list)->initial_state == S(state_groups) && list_index(associates, T(list)->initial_state) != -1) {
1879 push(&out_trans, T(list), 1);
1883 n = list_length(out_trans);
1884 S(state_groups)->num_out_trans = n;
1885 S(state_groups)->out_trans = (transition_t **)malloc(sizeof(transition_t *)*n);
1886 assert(S(state_groups)->out_trans != NULL);
1887 for (i=0; i<n; i++) {
1888 S(state_groups)->out_trans[i] = (transition_t *)pop(&out_trans);
1891 /* find states grouped with this state */
1893 fprintf(stderr, "finding states grouped with %s\n", S(state_groups)->name);
1895 handler = S(state_groups)->tension_handler;
1896 group = S(state_groups)->tension_group;
1897 list = head(state_groups);
1898 while (list != NULL) {
1899 if (S(list)->tension_handler == handler && S(list)->tension_group == group && list != state_groups) {
1900 push(&associates, S(list), 1);
1904 n = list_length(associates);
1905 S(state_groups)->num_grouped_states = n;
1906 S(state_groups)->grouped_states = (state_t **)malloc(sizeof(state_t *)*n);
1907 assert(S(state_groups)->grouped_states != NULL);
1908 for (i=0; i<n; i++) {
1909 S(state_groups)->grouped_states[i] = (state_t *)pop(&associates);
1911 state_groups = state_groups->next;
1917 \section{String parsing}
1919 For handling command line arguments, we need parse delimited strings (e.g.\ [["param1,param2,param3..."]]).
1920 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1921 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1922 need to take care of parsing those parameters themselves.
1923 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1929 <<parse definitions>>
1930 <<parse declarations>>
1934 <<parse module makefile lines>>=
1935 NW_SAWSIM_MODS += parse
1936 CHECK_BINS += check_parse
1940 <<parse definitions>>=
1941 #define SEP ',' /* argument separator character */
1944 <<parse declarations>>=
1945 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1946 int *num, char ***string_array);
1947 extern void free_parsed_list(int num, char **string_array);
1950 [[parse_list_string]] allocates memory, don't forget to free it
1951 afterward with [[free_parsed_list]]. It does not alter the original.
1953 The string may start off with a [[deeper]] character
1954 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1955 [[x,y]], where the pointer is one character in on the copied string.
1956 However, when we go to free the memory, we need a pointer to the
1957 beginning of the string. In order to accommodate this for a string
1958 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1959 the first $N$ elements point to the separated arguments, and let the
1960 last element point to the start of the copied string regardless of
1962 <<parse delimited list string functions>>=
1963 /* TODO, split out into parse.hc */
1964 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1967 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1968 if (string[i] == deeper) {depth++;}
1969 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1975 void parse_list_string(char *string, char sep, char deeper, char shallower,
1976 int *num, char ***string_array)
1978 char *str=NULL, **ret=NULL;
1980 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
1982 *string_array = NULL;
1985 /* make a copy of the string, so we don't change the original */
1986 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1987 assert(str != NULL);
1988 strcpy(str, string); /* we know str is long enough */
1989 /* count the number of regions, so we can allocate pointers to them */
1992 n++; i++; /* move on to next argument */
1993 i += next_delim_index(str+i, sep, deeper, shallower);
1994 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1995 } while (str[i] != '\0');
1996 ret = (char **)malloc(sizeof(char *)*(n+1));
1997 assert(ret != NULL);
1998 /* replace the separators with '\0' & assign pointers */
1999 ret[n] = str; /* point to the front of the copied string */
2002 for(i=1; i<n; i++) {
2003 j += next_delim_index(str+j, sep, deeper, shallower);
2005 ret[i] = str+j; /* point to the separated arguments */
2007 /* strip off bounding braces, if any */
2008 for(i=0; i<n; i++) {
2009 if (ret[i][0] == deeper) {
2013 if (ret[i][strlen(ret[i])-1] == shallower)
2014 ret[i][strlen(ret[i])-1] = '\0';
2017 *string_array = ret;
2020 void free_parsed_list(int num, char **string_array)
2023 /* frees all items, since they were allocated as a single string*/
2024 free(string_array[num]);
2025 /* and free the array of pointers */
2031 \subsection{Parsing implementation}
2037 <<parse delimited list string functions>>
2041 #include <assert.h> /* assert() */
2042 #include <stdlib.h> /* NULL */
2043 #include <stdio.h> /* fprintf(), stdout *//*!!*/
2044 #include <string.h> /* strlen() */
2048 \subsection{Parsing unit tests}
2050 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2052 <<parse check includes>>
2053 <<string comparison definition and globals>>
2054 <<check relative error>>
2055 <<parse test suite>>
2056 <<main check program>>
2059 <<parse check includes>>=
2060 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2061 #include <stdio.h> /* printf() */
2062 #include <assert.h> /* assert() */
2063 #include <string.h> /* strlen() */
2068 <<parse test suite>>=
2069 <<parse list string tests>>
2070 <<parse suite definition>>
2073 <<parse suite definition>>=
2074 Suite *test_suite (void)
2076 Suite *s = suite_create ("k model");
2077 <<parse list string test case defs>>
2079 <<parse list string test case add>>
2084 <<parse list string tests>>=
2087 START_TEST(test_next_delim_index)
2089 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
2090 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
2091 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
2092 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
2093 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
2097 START_TEST(test_parse_list_null)
2101 parse_list_string(NULL, SEP, '{', '}',
2102 &num_param_args, ¶m_args);
2103 fail_unless(num_param_args == 0, NULL);
2104 fail_unless(param_args == NULL, NULL);
2107 START_TEST(test_parse_list_single_simple)
2111 parse_list_string("arg", SEP, '{', '}',
2112 &num_param_args, ¶m_args);
2113 fail_unless(num_param_args == 1, NULL);
2114 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
2117 START_TEST(test_parse_list_single_compound)
2121 parse_list_string("{x,y,z}", SEP, '{', '}',
2122 &num_param_args, ¶m_args);
2123 fail_unless(num_param_args == 1, NULL);
2124 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
2127 START_TEST(test_parse_list_double_simple)
2131 parse_list_string("abc,def", SEP, '{', '}',
2132 &num_param_args, ¶m_args);
2133 fail_unless(num_param_args == 2, NULL);
2134 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2135 fail_unless(STRMATCH(param_args[1],"def"), NULL);
2138 START_TEST(test_parse_list_compound)
2142 parse_list_string("abc,{def,ghi}", SEP, '{', '}',
2143 &num_param_args, ¶m_args);
2144 fail_unless(num_param_args == 2, NULL);
2145 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2146 fail_unless(STRMATCH(param_args[1],"def,ghi"), NULL);
2150 <<parse list string test case defs>>=
2151 TCase *tc_parse_list_string = tcase_create("parse list string");
2153 <<parse list string test case add>>=
2154 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
2155 tcase_add_test(tc_parse_list_string, test_parse_list_null);
2156 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
2157 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
2158 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
2159 tcase_add_test(tc_parse_list_string, test_parse_list_compound);
2160 suite_add_tcase(s, tc_parse_list_string);
2164 \section{Unit tests}
2166 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2173 <<check relative error>>
2176 <<main check program>>
2188 <<determine dt tests>>
2190 <<does domain transition tests>>
2191 <<randomly unfold domains tests>>
2192 <<suite definition>>
2195 <<suite definition>>=
2196 Suite *test_suite (void)
2198 Suite *s = suite_create ("sawsim");
2199 <<F test case defs>>
2200 <<determine dt test case defs>>
2201 <<happens test case defs>>
2202 <<does domain transition test case defs>>
2203 <<randomly unfold domains test case defs>>
2206 <<determine dt test case add>>
2207 <<happens test case add>>
2208 <<does domain transition test case add>>
2209 <<randomly unfold domains test case add>>
2212 tcase_add_checked_fixture(tc_strip_address,
2213 setup_strip_address,
2214 teardown_strip_address);
2220 <<main check program>>=
2224 Suite *s = test_suite();
2225 SRunner *sr = srunner_create(s);
2226 srunner_run_all(sr, CK_ENV);
2227 number_failed = srunner_ntests_failed(sr);
2229 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
2233 \subsection{F tests}
2235 <<worm-like chain tests>>
2237 <<F test case defs>>=
2238 <<worm-like chain test case def>>
2240 <<F test case add>>=
2241 <<worm-like chain test case add>>
2244 <<worm-like chain tests>>=
2245 <<worm-like chain function>>
2247 START_TEST(test_wlc_at_zero)
2249 double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
2250 fail_unless(abs(wlc(x, T, p, L)) < lim, \
2251 "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
2252 x, T, p, L, abs(wlc(x, T, p, L)), lim);
2256 START_TEST(test_wlc_at_half)
2258 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
2259 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
2260 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
2262 * linear term = x/L = 0.5
2263 * nonlinear + linear = 0.75 + 0.5 = 1.25
2264 * wlc = 10e21*1.25 = 12.5e21
2266 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
2267 "wlc(%g, %g, %g, %g) = %g != %g",
2268 x, T, p, L, wlc(x, T, p, L), 12.5e21);
2272 <<worm-like chain test case def>>=
2273 TCase *tc_wlc = tcase_create("WLC");
2276 <<worm-like chain test case add>>=
2277 tcase_add_test(tc_wlc, test_wlc_at_zero);
2278 tcase_add_test(tc_wlc, test_wlc_at_half);
2279 suite_add_tcase(s, tc_wlc);
2282 \subsection{Model tests}
2284 Check the searching with [[linear_k]].
2285 Check overwhelming force treatment with the heavyside-esque step [[?]].
2287 Hmm, I'm not really sure what I was doing here...
2288 <<determine dt tests>>=
2289 double linear_k(double F, environment_t *env, void *params)
2291 double Fnot = *(double *)params;
2295 START_TEST(test_determine_dt_linear_k)
2298 transition_t unfold={0};
2299 environment_t env = {0};
2300 double F[]={0,1,2,3,4,5,6};
2301 double dt_max=3.0, Fnot=3.0, x=1.0, max_p=0.1, max_dF=0.1, v=1.0;
2302 list_t *state_list=NULL, *transition_list=NULL;
2305 folded.tension_handler = &hooke_handler;
2306 folded.num_domains = 1;
2307 unfold.initial_state = &folded;
2308 unfold.k = linear_k;
2309 unfold.k_params = &Fnot;
2310 push(&state_list, &folded, 1);
2311 push(&transition_list, &unfold, 1);
2312 for( i=0; i < sizeof(F)/sizeof(double); i++) {
2313 //fail_unless(determine_dt(state_list, transition_list, &env, x, max_p, max_dF, dt_max, v)==1, NULL);
2318 <<determine dt test case defs>>=
2319 TCase *tc_determine_dt = tcase_create("Determine dt");
2321 <<determine dt test case add>>=
2322 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
2323 suite_add_tcase(s, tc_determine_dt);
2328 <<happens test case defs>>=
2330 <<happens test case add>>=
2333 <<does domain transition tests>>=
2335 <<does domain transition test case defs>>=
2337 <<does domain transition test case add>>=
2340 <<randomly unfold domains tests>>=
2342 <<randomly unfold domains test case defs>>=
2344 <<randomly unfold domains test case add>>=
2348 \section{Balancing group extensions}
2349 \label{sec.tension_balance}
2351 For a given total extension $x$ (of the piezo), the various domain
2352 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
2353 amounts, and we need to tweak the portion that each extends to balance
2354 the tension amongst the active groups. Since the tension balancing is
2355 separable from the bulk of the simulation, we break it out into a
2356 separate module. The interface is defined in [[tension_balance.h]],
2357 the implementation in [[tension_balance.c]], and the unit testing in
2358 [[check_tension_balance.c]]
2360 <<tension-balance.h>>=
2361 #ifndef TENSION_BALANCE_H
2362 #define TENSION_BALANCE_H
2364 <<tension balance function declaration>>
2365 #endif /* TENSION_BALANCE_H */
2368 <<tension balance functions>>=
2369 <<one dimensional bracket>>
2370 <<one dimensional solve>>
2372 <<group stiffness function>>
2373 <<chain stiffness function>>
2374 <<full chain stiffness function>>
2375 <<tension balance function>>
2378 <<tension balance module makefile lines>>=
2379 NW_SAWSIM_MODS += tension_balance
2380 CHECK_BINS += check_tension_balance
2381 check_tension_balance_MODS =
2384 The entire force balancing problem reduces to a solving two nested
2385 one-dimensional root-finding problems. First we define the one
2386 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
2387 populated group). $x(x_0)$ is also strictly monotonically increasing,
2388 so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
2389 stores the last successful value of $x$, since we expect to be taking
2390 small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
2391 indicates that the groups have changed.
2392 <<tension balance function declaration>>=
2393 double tension_balance(int num_tension_groups,
2394 one_dim_fn_t **tension_handlers,
2395 one_dim_fn_t **inverse_tension_handlers,
2396 void **params, /* array of pointers to tension_handler_data_t */
2401 <<tension balance function>>=
2402 double tension_balance(int num_tension_groups,
2403 one_dim_fn_t **tension_handlers,
2404 one_dim_fn_t **inverse_tension_handlers,
2405 void **params, /* array of pointers to tension_handler_data_t */
2410 static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL,NULL};
2411 double F, xo_guess, xo, lb, ub;
2412 double min_relx=1e-6, min_rely=1e-6;
2413 int max_steps=200, i;
2415 assert(num_tension_groups > 0);
2416 assert(tension_handlers != NULL);
2417 assert(params != NULL);
2418 assert(num_tension_groups > 0);
2420 if (last_x == -1 || x_xo_data.xi == NULL) { /* new group setup, reset x_xo_data */
2421 /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
2422 if (x_xo_data.xi != NULL)
2424 /* construct new x_xo_data */
2425 x_xo_data.n_groups = num_tension_groups;
2426 x_xo_data.tension_handler = tension_handlers;
2427 x_xo_data.inverse_tension_handler = inverse_tension_handlers;
2428 x_xo_data.handler_data = params;
2430 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);
2432 x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
2433 assert(x_xo_data.xi != NULL);
2434 for (i=0; i<num_tension_groups; i++)
2435 x_xo_data.xi[i] = last_x;
2438 if (num_tension_groups == 1) { /* only one group, no balancing required */
2440 x_xo_data.xi[0] = xo;
2442 fprintf(stderr, "%s:\tonly one active group, no balancing required\n", __FUNCTION__);
2446 fprintf(stderr, "%s:\tbalancing for x = %g with", __FUNCTION__, x);
2447 for(i=0; i<num_tension_groups; i++) fprintf(stderr, "\tx[%d]=%g", i, x_xo_data.xi[i]);
2448 fprintf(stderr, "\n");
2450 if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
2451 if (x == 0) xo_guess = length_scale;
2452 else xo_guess = x/num_tension_groups;
2454 fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
2456 oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
2457 } else { /* work off the last balanced point */
2459 fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
2462 lb = x_xo_data.xi[0];
2463 ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */
2464 } else if (x < last_x) {
2465 lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
2466 ub = x_xo_data.xi[0];
2467 } else { /* x == last_x */
2469 fprintf(stderr, "not moving\n");
2471 lb = x_xo_data.xi[0];
2472 ub = x_xo_data.xi[0];
2476 fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
2477 __FUNCTION__, x, lb, ub);
2479 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
2480 F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
2482 if (num_tension_groups > 1) {
2483 fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
2484 for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
2485 fprintf(stderr, "\n");
2490 F = (*tension_handlers[0])(xo, params[0]);
2492 if (stiffness != NULL) {
2493 *stiffness = chain_stiffness(&x_xo_data, min_relx);
2494 if (*stiffness == 0.0) { /* re-calculate based on full chain */
2495 *stiffness = full_chain_stiffness(num_tension_groups, tension_handlers,
2496 inverse_tension_handlers, params,
2497 last_x, x, min_relx, F);
2505 <<tension balance internal definitions>>=
2506 <<x of x0 definitions>>
2509 <<x of x0 definitions>>=
2510 typedef struct x_of_xo_data_struct {
2512 one_dim_fn_t **tension_handler; /* array of fn pointers */
2513 one_dim_fn_t **inverse_tension_handler; /* inverse fn pointers */
2514 void **handler_data; /* array of pointers to tension_handler_data_t */
2515 double *xi; /* array of group extensions */
2520 double x_of_xo(double xo, void *pdata)
2522 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
2523 double F, x=0, xi, xi_guess, lb, ub, slope;
2525 double min_relx=1e-6, min_rely=1e-6;
2527 assert(data->n_groups > 0);
2530 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);
2533 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
2535 fprintf(stderr, "%s:\tfinding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
2539 if (data->xi) data->xi[0] = xo;
2540 for (i=1; i < data->n_groups; i++) {
2541 if (data->inverse_tension_handler[i] != NULL) {
2542 xi = (*data->inverse_tension_handler[i])(F, data->handler_data[i]);
2543 } else { /* invert numerically */
2544 if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
2545 else xi_guess = data->xi[i];
2547 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]);
2549 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
2550 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
2555 fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
2559 fprintf(stderr, "%s:\tfound x(x0=%g) = %g\n", __FUNCTION__, xo, x);
2565 Determine the stiffness of a single tension group by taking a small
2566 step [[dx]] back from the position [[x]] for which we want the
2567 stiffness. The touchy part is determining a good [[dx]] and ensuring
2568 the whole interval is on [[x>=0]].
2569 <<group stiffness function>>=
2570 double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double x, double relx)
2572 double dx, xl, F, Fl, stiffness;
2574 assert(relx > 0 && relx < 1);
2575 if (x == 0 || relx == 0) {
2576 dx = 1e-12; /* HACK, how to get length scale? */
2586 F = (*tension_handler)(x, handler_data);
2587 Fl = (*tension_handler)(xl, handler_data);
2588 stiffness = (F-Fl)/dx;
2589 assert(stiffness >= 0);
2594 Attempt to calculate the chain stiffness by adding group stiffnesses
2595 as springs in series. In the case where a group has undefined
2596 stiffness (e.g. piston tension, Section \ref{sec.piston_tension}),
2597 this algorithm breaks down. In that case, [[tension_balance()]]
2598 (\ref{sec.tension_balance}) should use [[full_chain_stiffness()]]
2599 which uses the full chain's $F(x)$ rather than that of the individual
2600 domains'. We attempt the series approach first because, lacking the
2601 numerical inversion inside [[tension_balance()]], it is both faster
2602 and more accurate than the full-chain derivative.
2603 <<chain stiffness function>>=
2604 double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx)
2606 double stiffness, gstiff;
2608 /* add group stiffnesses in series */
2609 for (i=0; i < x_xo_data->n_groups; i++) {
2611 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);
2614 gstiff = group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx);
2617 stiffness += 1.0/gstiff;
2619 stiffness = 1.0 / stiffness;
2625 Determine the chain stiffness using only [[tension_balance()]]. This
2626 works around difficulties with tension models that have undefined
2627 stiffness (see discussion for [[chain_stiffness()]] above).
2628 <<full chain stiffness function>>=
2629 double full_chain_stiffness(int num_tension_groups,
2630 one_dim_fn_t **tension_handlers,
2631 one_dim_fn_t **inverse_tension_handlers,
2632 void **params, /* array of pointers to tension_handler_data_t */
2636 double F /* already evaluated F(x) */)
2638 double dx, xl, Fl, stiffness;
2641 assert(relx > 0 && relx < 1);
2643 /* Other option for dx that we ignore because of our poor tension_balance()
2644 * resolution (i.e. bad slopes if you zoom in too close):
2645 * if (last_x != -1.0 && last_x != x) // last_x limited
2646 * dx fabs(x-last_x);
2649 * if (1==1) { // constant max-value limited
2651 * dx = MIN(dx, dx_p);
2653 * if (x != 0 && relx != 0) { // relx limited
2655 * dx = MIN(dx, dx_p);
2657 * TODO, determine which of (min_relx, min_rely, max_steps) is actually
2658 * limiting tension_balance.
2660 dx = 1e-12; /* HACK, how to get length scale? */
2664 Fl = tension_balance(num_tension_groups, tension_handlers,
2665 inverse_tension_handlers, params,
2671 F = tension_balance(num_tension_groups, tension_handlers,
2672 inverse_tension_handlers, params,
2675 stiffness = (F-Fl)/dx;
2676 assert(stiffness >= 0);
2682 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
2683 number of steps for monotonically increasing $f(x)$. Simple
2684 bisection, so it's robust and converges fairly quickly. You can set
2685 the maximum number of steps to take, but if you have enough processor
2686 power, it's better to limit your precision with the relative limits
2687 [[min_relx]] and [[y]]. These ensure that the width of the bracket is
2688 small on the length scale of it's larger side. Note that \emph{both}
2689 relative limits must be satisfied for the search to stop.
2690 <<one dimensional function definition>>=
2691 /* equivalent to gsl_func_t */
2692 typedef double one_dim_fn_t(double x, void *params);
2694 <<one dimensional solve declaration>>=
2695 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2696 double min_relx, double min_rely, int max_steps,
2700 <<one dimensional solve>>=
2701 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2702 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2703 double min_relx, double min_rely, int max_steps,
2706 double yx, ylb, yub, x;
2709 ylb = (*f)(lb, params);
2710 yub = (*f)(ub, params);
2712 /* check some simple cases */
2714 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bracketed */
2715 else return lb; /* any x on the range [lb, ub] would work */
2717 if (ylb == y) { x=lb; yx=ylb; goto end; }
2718 if (yub == y) { x=ub; yx=yub; goto end; }
2721 fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2723 assert(ylb < y); assert(yub > y);
2724 for (i=0; i<max_steps; i++) {
2726 yx = (*f)(x, params);
2728 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);
2730 if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2731 else if (yx > y) { ub=x; yub=yx; }
2732 else /* yx < y */ { lb=x; ylb=yx; }
2737 fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2743 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2744 Generate bracketing $x$ values through bisection or exponential growth.
2745 <<one dimensional bracket declaration>>=
2746 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2749 <<one dimensional bracket>>=
2750 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2752 double yx, step, x=xguess;
2755 fprintf(stderr, "%s:\tbracketing %g (params %p)\n", __FUNCTION__, y, params);
2757 yx = (*f)(x, params);
2759 fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2766 fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2770 if (x == 0) x = length_scale/2.0; /* HACK */
2773 yx = (*f)(x, params);
2775 fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2780 fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2783 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2787 \subsection{Balancing implementation}
2789 <<tension-balance.c>>=
2792 <<tension balance includes>>
2793 #include "tension_balance.h"
2795 double length_scale = 1e-10; /* HACK */
2797 <<tension balance internal definitions>>
2798 <<tension balance functions>>
2801 <<tension balance includes>>=
2802 #include <assert.h> /* assert() */
2803 #include <stdlib.h> /* NULL */
2804 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2805 #include <stdio.h> /* fprintf(), stdout */
2809 \subsection{Balancing unit tests}
2811 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2812 <<check-tension-balance.c>>=
2813 <<tension balance check includes>>
2814 <<tension balance test suite>>
2815 <<main check program>>
2818 <<tension balance check includes>>=
2819 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2820 #include <assert.h> /* assert() */
2823 #include "tension_balance.h"
2826 <<tension balance test suite>>=
2827 <<tension balance function tests>>
2828 <<tension balance suite definition>>
2831 <<tension balance suite definition>>=
2832 Suite *test_suite (void)
2834 Suite *s = suite_create ("tension balance");
2835 <<tension balance function test case defs>>
2837 <<tension balance function test case adds>>
2842 <<tension balance function tests>>=
2843 <<check relative error>>
2845 double hooke(double x, void *pK)
2848 return *((double*)pK) * x;
2851 START_TEST(test_single_function)
2853 double k=5, x=3, last_x=2.0, F;
2854 one_dim_fn_t *handlers[] = {&hooke};
2855 one_dim_fn_t *inverse_handlers[] = {NULL};
2856 void *data[] = {&k};
2857 F = tension_balance(1, handlers, inverse_handlers, data, last_x, x, NULL);
2858 fail_unless(F = k*x, NULL);
2863 We can also test balancing two springs with different spring constants.
2864 The analytic solution for a general number of springs is given in Appendix \ref{sec.math_hooke}.
2865 <<tension balance function tests>>=
2866 START_TEST(test_double_hooke)
2868 double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
2869 one_dim_fn_t *handlers[] = {&hooke, &hooke};
2870 one_dim_fn_t *inverse_handlers[] = {NULL, NULL};
2871 void *data[] = {&k1, &k2};
2872 F = tension_balance(2, handlers, inverse_handlers, data, last_x, x, NULL);
2876 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2877 //CHECK_ERR(1e-6, x1e, xi[0]);
2878 //CHECK_ERR(1e-6, x2e, xi[1]);
2879 CHECK_ERR(1e-6, Fe, F);
2883 TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
2885 <<tension balance function test case defs>>=
2886 TCase *tc_tbfunc = tcase_create("tension balance function");
2889 <<tension balance function test case adds>>=
2890 tcase_add_test(tc_tbfunc, test_single_function);
2891 tcase_add_test(tc_tbfunc, test_double_hooke);
2892 suite_add_tcase(s, tc_tbfunc);
2898 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2899 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2900 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2906 <<list definitions>>
2907 <<list declarations>>
2911 <<list declarations>>=
2912 <<head and tail declarations>>
2913 <<list length declaration>>
2914 <<push declaration>>
2916 <<list destroy declaration>>
2917 <<list to array declaration>>
2918 <<move declaration>>
2919 <<sort declaration>>
2920 <<list index declaration>>
2921 <<list element declaration>>
2922 <<unique declaration>>
2923 <<list print declaration>>
2927 <<create and destroy node>>
2942 <<list module makefile lines>>=
2943 NW_SAWSIM_MODS += list
2944 CHECK_BINS += check_list
2948 <<list definitions>>=
2949 typedef struct list_struct {
2950 struct list_struct *next;
2951 struct list_struct *prev;
2955 <<comparison function typedef>>
2958 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2959 <<head and tail declarations>>=
2960 list_t *head(list_t *list);
2961 list_t *tail(list_t *list);
2964 list_t *head(list_t *list)
2968 while (list->prev != NULL) {
2974 list_t *tail(list_t *list)
2978 while (list->next != NULL) {
2985 <<list length declaration>>=
2986 int list_length(list_t *list);
2989 int list_length(list_t *list)
2996 while (list->next != NULL) {
3004 [[push]] inserts a node relative to the current position in the list
3005 without changing the current position.
3006 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
3007 so we set it to that of the pushed domain.
3008 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
3009 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
3010 <<push declaration>>=
3011 void push(list_t **pList, void *data, int below);
3014 void push(list_t **pList, void *data, int below)
3016 list_t *list, *new_node;
3017 assert(pList != NULL);
3019 new_node = create_node(data);
3022 else if (below==0) { /* insert above */
3023 if (list->prev != NULL)
3024 list->prev->next = new_node;
3025 new_node->prev = list->prev;
3026 list->prev = new_node;
3027 new_node->next = list;
3028 } else { /* insert below */
3029 if (list->next != NULL)
3030 list->next->prev = new_node;
3031 new_node->next = list->next;
3032 list->next = new_node;
3033 new_node->prev = list;
3038 [[pop]] removes the current domain node, moving the current position
3039 to the node after it, unless that node is [[NULL]], in which case move
3040 the current position to the node before it.
3041 <<pop declaration>>=
3042 void *pop(list_t **pList);
3045 void *pop(list_t **pList)
3047 list_t *list, *popped;
3049 assert(pList != NULL);
3051 assert(list != NULL); /* not an empty list */
3053 /* bypass the popped node */
3054 if (list->prev != NULL)
3055 list->prev->next = list->next;
3056 if (list->next != NULL) {
3057 list->next->prev = list->prev;
3058 *pList = list->next;
3060 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
3062 destroy_node(popped);
3067 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
3068 If the cleanup function is [[NULL]], no cleanup function is called.
3069 The nodes are not popped in any particular order.
3070 <<list destroy declaration>>=
3071 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
3074 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
3078 assert(pList != NULL);
3081 assert(list != NULL); /* not an empty list */
3082 while (list != NULL) {
3084 if (destroy != NULL)
3090 [[list_to_array]] converts a list to an array of pointers, preserving order.
3091 <<list to array declaration>>=
3092 void list_to_array(list_t *list, int *length, void ***array);
3095 void list_to_array(list_t *list, int *pLength, void ***pArray)
3099 assert(list != NULL);
3100 assert(pLength != NULL);
3101 assert(pArray != NULL);
3103 length = list_length(list);
3104 array = (void **)malloc(sizeof(void *)*length);
3105 assert(array != NULL);
3108 while (list != NULL) {
3109 array[i++] = list->d;
3117 [[move]] moves the current element along the list in either direction.
3118 <<move declaration>>=
3119 void move(list_t **pList, int downstream);
3122 void move(list_t **pList, int downstream)
3124 list_t *list, *flipper;
3126 assert(pList != NULL);
3127 list = *pList; /* a>B>c>d */
3128 assert(list != NULL); /* not an empty list */
3129 if (downstream == 0)
3130 flipper = list->prev; /* flipper is A */
3132 flipper = list->next; /* flipper is C */
3133 /* ensure there is somewhere to go in the direction we're heading */
3134 assert(flipper != NULL);
3135 /* remove the current element */
3136 data = pop(&list); /* data=B, a>C>d */
3137 /* and put it back in in the correct spot */
3139 if (downstream == 0) {
3140 push(&list, data, 0); /* b>A>c>d */
3141 list = list->prev; /* B>a>c>d */
3143 push(&list, data, 1); /* a>C>b>d */
3144 list = list->next; /* a>c>B>d */
3150 [[sort]] sorts a list from smallest to largest according to a user
3151 specified comparison.
3152 <<comparison function typedef>>=
3153 typedef int (comparison_fn_t)(void *A, void *B);
3156 <<int comparison function>>=
3157 int int_comparison(void *A, void *B)
3162 if (a > b) return 1;
3163 else if (a == b) return 0;
3168 <<sort declaration>>=
3169 void sort(list_t **pList, comparison_fn_t *comp);
3171 Here we use the bubble sort algorithm.
3173 void sort(list_t **pList, comparison_fn_t *comp)
3177 assert(pList != NULL);
3179 assert(list != NULL); /* not an empty list */
3180 while (stable == 0) {
3183 while (list->next != NULL) {
3184 if ((*comp)(list->d, list->next->d) > 0) {
3186 move(&list, 1 /* downstream */);
3196 [[list_index]] finds the location of [[data]] in [[list]]. Returns
3197 [[-1]] if [[data]] is not in [[list]].
3198 <<list index declaration>>=
3199 int list_index(list_t *list, void *data);
3203 int list_index(list_t *list, void *data)
3207 while (list != NULL) {
3208 if (list->d == data) return i;
3217 [[list_element]] returns the data bound to the $i^\text{th}$ node of the list.
3218 <<list element declaration>>=
3219 void *list_element(list_t *list, int i);
3222 void *list_element(list_t *list, int i)
3226 while (list != NULL) {
3227 if (i == j) return list->d;
3237 [[unique]] removes duplicates from a sorted list according to a user
3238 specified comparison. Don't do this unless you have other pointers
3239 any dynamically allocated data in your list, because [[unique]] just
3240 drops any non-unique data without freeing it.
3241 <<unique declaration>>=
3242 void unique(list_t **pList, comparison_fn_t *comp);
3245 void unique(list_t **pList, comparison_fn_t *comp)
3248 assert(pList != NULL);
3250 assert(list != NULL); /* not an empty list */
3252 while (list->next != NULL) {
3253 if ((*comp)(list->d, list->next->d) == 0) {
3262 [[list_print]] prints a list to a [[FILE *]] stream.
3263 <<list print declaration>>=
3264 void list_print(FILE *file, list_t *list, char *list_name);
3267 void list_print(FILE *file, list_t *list, char *list_name)
3269 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
3271 while (list != NULL) {
3272 fprintf(file, " %p", list->d);
3275 fprintf(file, "\n");
3279 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
3280 <<create and destroy node>>=
3281 list_t *create_node(void *data)
3283 list_t *ret = (list_t *)malloc(sizeof(list_t));
3284 assert(ret != NULL);
3291 void destroy_node(list_t *node)
3297 The user must free the data pointed to by the node on their own.
3299 \subsection{List implementation}
3309 #include <assert.h> /* assert() */
3310 #include <stdlib.h> /* malloc(), free(), rand() */
3311 #include <stdio.h> /* fprintf(), stdout, FILE */
3312 #include "global.h" /* destroy_data_func_t */
3315 \subsection{List unit tests}
3317 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3319 <<list check includes>>
3321 <<main check program>>
3324 <<list check includes>>=
3325 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3326 #include <stdio.h> /* FILE */
3332 <<list test suite>>=
3335 <<list suite definition>>
3338 <<list suite definition>>=
3339 Suite *test_suite (void)
3341 Suite *s = suite_create ("list");
3342 <<push test case defs>>
3343 <<pop test case defs>>
3345 <<push test case adds>>
3346 <<pop test case adds>>
3352 START_TEST(test_push)
3357 push(&list, (void *)i, 1);
3358 fail_unless(list == head(list), NULL);
3359 fail_unless( (int)list->d == 0 );
3360 for (i=0; i<n; i++) {
3361 /* we expect 0, n-1, n-2, ..., 1 */
3364 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
3370 <<push test case defs>>=
3371 TCase *tc_push = tcase_create("push");
3374 <<push test case adds>>=
3375 tcase_add_test(tc_push, test_push);
3376 suite_add_tcase(s, tc_push);
3381 <<pop test case defs>>=
3383 <<pop test case adds>>=
3386 \section{Function string evaluation}
3388 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).
3389 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
3390 The solution is to run a scripting language as a subprocess accessed via pipes.
3391 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
3393 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
3394 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
3395 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.
3396 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
3398 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.
3399 Then you can either statically or dynamically link to those hardcoded functions.
3400 While much less flexible, this approach would be a fairly simple-to-implement fallback.
3402 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
3403 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
3406 #ifndef STRING_EVAL_H
3407 #define STRING_EVAL_H
3409 <<string eval setup declaration>>
3410 <<string eval function declaration>>
3411 <<string eval teardown declaration>>
3412 #endif /* STRING_EVAL_H */
3415 <<string eval module makefile lines>>=
3416 NW_SAWSIM_MODS += string_eval
3417 CHECK_BINS += check_string_eval
3418 check_string_eval_MODS =
3421 For an introduction to POSIX process control, see\\
3422 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
3423 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
3424 and of course, the relavent [[man]] pages.
3426 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
3427 [[execvp]] replaces the calling process' program with a new program.
3428 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
3429 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
3430 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
3431 The new program needs command line arguments, just like it would if you were running it from a shell.
3432 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
3433 with the final array entry being a [[NULL]] pointer.
3435 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
3436 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
3437 <<string eval subprocess definitions>>=
3438 #define SUBPROCESS "python"
3439 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
3440 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
3441 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
3443 The [[i]] option lets Python know that it should run in interactive mode.
3444 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
3445 In interactive mode, python acts on every instruction as soon as it is recieved.
3446 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.
3447 %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.
3449 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
3450 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
3451 [[fork]] returns two copies of the same program, executing the original code.
3452 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
3453 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
3455 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.
3456 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
3457 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
3458 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
3459 <<string eval pipe definitions>>=
3460 #define PIPE_READ 0 /* the end you read from */
3461 #define PIPE_WRITE 1 /* the end you write to */
3463 #define STDIN 0 /* initial index of stdin pair */
3464 #define STDOUT 2 /* initial index of stdout pair */
3467 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
3469 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.
3471 <<string eval setup declaration>>=
3472 extern void string_eval_setup(FILE **pIn, FILE **pOut);
3474 <<string eval setup definition>>=
3475 void string_eval_setup(FILE **pIn, FILE **pOut)
3478 int pfd[4]; /* file descriptors for stdin and stdout */
3480 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
3481 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
3483 pid = fork(); /* split process into two copies */
3485 if (pid == -1) { /* fork error */
3486 perror("fork error");
3488 } else if (pid == 0) { /* child */
3489 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
3490 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
3491 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
3492 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
3493 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
3494 perror("exec error"); /* exec shouldn't return */
3496 } else { /* parent */
3497 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
3498 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
3499 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
3500 if ( *pIn == NULL ) {
3501 perror("fdopen (in)");
3504 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
3505 if ( *pOut == NULL ) {
3506 perror("fdopen (out)");
3513 To use the evaluating subprocess, we just pipe in our command, and read out the results.
3514 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
3515 <<string eval function declaration>>=
3516 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
3518 <<string eval function definition>>=
3519 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
3522 rc = fprintf(in, "%s", input);
3523 assert(rc == strlen(input));
3526 alarm(1); /* set a one second timeout on the read */
3527 assert( fgets(output, buflen, out) != NULL );
3528 alarm(0); /* cancel the timeout */
3529 //fprintf(stderr, "eval: %s ----> %s", input, output);
3532 The [[alarm]] calls set and clear a timeout on the returned output.
3533 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.
3534 This protects against invalid input for which a line of output is not printed to [[stdout]].
3535 Other invalid inputs (e.g.\ those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
3536 If you are getting strange results, check your python code seperately. TODO, better error handling.
3538 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
3539 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
3540 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
3541 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
3542 <<string eval teardown declaration>>=
3543 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
3546 <<string eval teardown definition>>=
3547 void string_eval_teardown(FILE **pIn, FILE **pOut)
3552 /* redirect Python's stderr */
3553 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
3557 assert( fclose(*pIn) == 0 );
3559 assert( fclose(*pOut) == 0 );
3562 /* wait for python to exit */
3564 pid = wait(&stat_loc);
3571 if (WIFEXITED(stat_loc)) {
3572 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
3573 } else if (WIFSIGNALED(stat_loc)) {
3574 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
3579 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
3581 \subsection{String evaluation implementation}
3585 <<string eval includes>>
3586 #include "string_eval.h"
3587 <<string eval internal definitions>>
3588 <<string eval functions>>
3591 <<string eval includes>>=
3592 #include <assert.h> /* assert() */
3593 #include <stdlib.h> /* NULL */
3594 #include <stdio.h> /* fprintf(), stdout, fdopen() */
3595 #include <string.h> /* strlen() */
3596 #include <sys/types.h> /* pid_t */
3597 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
3598 #include <wait.h> /* wait() */
3601 <<string eval internal definitions>>=
3602 <<string eval subprocess definitions>>
3603 <<string eval pipe definitions>>
3606 <<string eval functions>>=
3607 <<string eval setup definition>>
3608 <<string eval function definition>>
3609 <<string eval teardown definition>>
3612 \subsection{String evaluation unit tests}
3614 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3615 <<check-string-eval.c>>=
3616 <<string eval check includes>>
3617 <<string comparison definition and globals>>
3618 <<string eval test suite>>
3619 <<main check program>>
3622 <<string eval check includes>>=
3623 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3624 #include <stdio.h> /* printf() */
3625 #include <assert.h> /* assert() */
3626 #include <string.h> /* strlen() */
3627 #include <signal.h> /* SIGKILL */
3629 #include "string_eval.h"
3632 <<string eval test suite>>=
3633 <<string eval tests>>
3634 <<string eval suite definition>>
3637 <<string eval suite definition>>=
3638 Suite *test_suite (void)
3640 Suite *s = suite_create ("string eval");
3641 <<string eval test case defs>>
3643 <<string eval test case add>>
3648 <<string eval tests>>=
3649 START_TEST(test_setup_teardown)
3652 string_eval_setup(&in, &out);
3653 string_eval_teardown(&in, &out);
3656 START_TEST(test_invalid_command)
3659 char input[80], output[80]={};
3660 string_eval_setup(&in, &out);
3661 sprintf(input, "print ABCDefg\n");
3662 string_eval(in, out, input, 80, output);
3663 string_eval_teardown(&in, &out);
3666 START_TEST(test_simple_eval)
3669 char input[80], output[80]={};
3670 string_eval_setup(&in, &out);
3671 sprintf(input, "print 3+4*5\n");
3672 string_eval(in, out, input, 80, output);
3673 fail_unless(STRMATCH(output,"23\n"), NULL);
3674 string_eval_teardown(&in, &out);
3677 START_TEST(test_multiple_evals)
3680 char input[80], output[80]={};
3681 string_eval_setup(&in, &out);
3682 sprintf(input, "print 3+4*5\n");
3683 string_eval(in, out, input, 80, output);
3684 fail_unless(STRMATCH(output,"23\n"), NULL);
3685 sprintf(input, "print (3**2 + 4**2)**0.5\n");
3686 string_eval(in, out, input, 80, output);
3687 fail_unless(STRMATCH(output,"5.0\n"), NULL);
3688 string_eval_teardown(&in, &out);
3691 START_TEST(test_eval_with_state)
3694 char input[80], output[80]={};
3695 string_eval_setup(&in, &out);
3696 sprintf(input, "print 3+4*5\n");
3697 fprintf(in, "A = 3\n");
3698 sprintf(input, "print A*3\n");
3699 string_eval(in, out, input, 80, output);
3700 fail_unless(STRMATCH(output,"9\n"), NULL);
3701 string_eval_teardown(&in, &out);
3705 <<string eval test case defs>>=
3706 TCase *tc_string_eval = tcase_create("string_eval");
3708 <<string eval test case add>>=
3709 tcase_add_test(tc_string_eval, test_setup_teardown);
3710 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
3711 tcase_add_test(tc_string_eval, test_simple_eval);
3712 tcase_add_test(tc_string_eval, test_multiple_evals);
3713 tcase_add_test(tc_string_eval, test_eval_with_state);
3714 suite_add_tcase(s, tc_string_eval);
3718 \section{Accelerating function evaluation}
3720 My first version-0.3 code was running very slowly.
3721 With the options suggested in the help
3722 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
3723 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
3724 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
3725 That's only 15 calls per solution, so the search algorithm seems reasonable.
3726 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
3731 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
3733 #endif /* ACCEL_K_H */
3736 <<accel k module makefile lines>>=
3737 NW_SAWSIM_MODS += accel_k
3738 #CHECK_BINS += check_accel_k
3739 check_accel_k_MODS =
3743 #include <assert.h> /* assert() */
3744 #include <stdlib.h> /* realloc(), free(), NULL */
3745 #include "global.h" /* environment_t */
3746 #include "k_model.h" /* k_func_t */
3747 #include "interp.h" /* interp_* */
3748 #include "accel_k.h"
3750 typedef struct accel_k_struct {
3751 interp_table_t *itable;
3757 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3758 static int num_accels = 0;
3759 static accel_k_t *accels=NULL;
3761 /* Wrap k in the standard f(x) acceleration form */
3762 static double k_wrap(double F, void *params)
3764 accel_k_t *p = (accel_k_t *)params;
3765 return (*p->k)(F, p->env, p->params);
3768 static int k_tol(double FA, double kA, double FB, double kB)
3771 if (FB-FA > 1e-12) {
3772 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3773 return 1; /* unnacceptable */
3775 //printf("acceptable tol\n");
3776 return 0; /* acceptable */
3780 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3783 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3784 assert(accels != NULL);
3785 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3787 accels[i].env = env;
3788 accels[i].params = params;
3795 for (i=0; i<num_accels; i++)
3796 interp_table_free(accels[i].itable);
3798 if (accels) free(accels);
3802 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3805 for (i=0; i<num_accels; i++) {
3806 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3807 /* We've already setup this function.
3808 * WARNING: we're assuming param and env are the same. */
3809 return interp_table_eval(accels[i].itable, accels+i, F);
3812 /* set up a new acceleration environment */
3813 i = add_accel_k(k, env, params);
3814 return interp_table_eval(accels[i].itable, accels+i, F);
3818 \section{Tension models}
3819 \label{sec.tension_models}
3821 TODO: tension model intro.
3822 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.
3824 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3825 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]].
3827 <<tension-model.h>>=
3828 #ifndef TENSION_MODEL_H
3829 #define TENSION_MODEL_H
3831 <<tension handler types>>
3832 <<constant tension model declarations>>
3833 <<hooke tension model declarations>>
3834 <<worm-like chain tension model declarations>>
3835 <<freely-jointed chain tension model declarations>>
3836 <<piston tension model declarations>>
3837 <<find tension definitions>>
3838 <<tension model global declarations>>
3839 #endif /* TENSION_MODEL_H */
3842 <<tension model module makefile lines>>=
3843 NW_SAWSIM_MODS += tension_model
3844 #CHECK_BINS += check_tension_model
3845 check_tension_model_MODS =
3847 <<tension model utils makefile lines>>=
3848 TENSION_MODEL_MODS = tension_model parse list tension_balance
3849 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3850 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3851 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3852 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
3853 notangle -Rtension-model-utils.c $< > $@
3854 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
3855 gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3856 clean_tension_model_utils :
3857 rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3858 @ The pipe symbol [[|]] marks the prerequisits following it (in this
3859 case, the directories) as ‘order-only’ prerequisites. The timestamp
3860 on these prerequisits does not effect whether the rules are executed.
3861 This is appropriate for directories, where we don't need to recompile
3862 after an unrelated has been added to the directory, but only when the
3863 source prerequisites change. See the [[make]] documentation for more
3865 (\url{http://www.gnu.org/software/make/manual/html_node/Prerequisite-Types.html}).
3868 \label{sec.null_tension}
3870 For unstretchable domains.
3872 <<null tension model getopt>>=
3873 {NULL, NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3876 \subsection{Constant}
3877 \label{sec.const_tension}
3879 An infinitely stretchable domain providing a constant tension.
3880 Obviously this cannot be inverted, so you can't put this domain in
3881 series with any other domains. We include it mostly for testing
3884 <<constant tension functions>>=
3885 <<constant tension function>>
3886 <<constant tension parameter create/destroy functions>>
3889 <<constant tension model declarations>>=
3890 <<constant tension function declaration>>
3891 <<constant tension parameter create/destroy function declarations>>
3892 <<constant tension model global declarations>>
3896 <<constant tension function declaration>>=
3897 extern double const_tension_handler(double x, void *pdata);
3899 <<constant tension function>>=
3900 double const_tension_handler(double x, void *pdata)
3902 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3903 list_t *list = data->group_tension_params;
3908 assert(list != NULL); /* empty active group?! */
3909 F = ((const_tension_param_t *)list->d)->F;
3911 fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3913 while (list != NULL) {
3914 assert(((const_tension_param_t *)list->d)->F == F);
3922 <<constant tension structure>>=
3923 typedef struct const_tension_param_struct {
3924 double F; /* tension (force) in N */
3925 } const_tension_param_t;
3929 <<constant tension parameter create/destroy function declarations>>=
3930 extern void *string_create_const_tension_param_t(char **param_strings);
3931 extern void destroy_const_tension_param_t(void *p);
3934 <<constant tension parameter create/destroy functions>>=
3935 const_tension_param_t *create_const_tension_param_t(double F)
3937 const_tension_param_t *ret
3938 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3939 assert(ret != NULL);
3944 void *string_create_const_tension_param_t(char **param_strings)
3946 return create_const_tension_param_t(
3947 safe_strtod(param_strings[0],__FUNCTION__));
3950 void destroy_const_tension_param_t(void *p)
3958 <<constant tension model global declarations>>=
3959 extern int num_const_tension_params;
3960 extern char *const_tension_param_descriptions[];
3961 extern char const_tension_param_string[];
3964 <<constant tension model globals>>=
3965 int num_const_tension_params = 1;
3966 char *const_tension_param_descriptions[] = {"tension F, N"};
3967 char const_tension_param_string[] = "0";
3970 <<constant tension model getopt>>=
3971 {&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}
3977 <<hooke functions>>=
3978 <<internal hooke functions>>
3980 <<inverse hooke handler>>
3981 <<hooke parameter create/destroy functions>>
3984 <<hooke tension model declarations>>=
3985 <<hooke tension function declaration>>
3986 <<hooke tension parameter create/destroy function declarations>>
3987 <<hooke tension model global declarations>>
3991 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
3992 The behavior of a series of springs $k_i$ in series is given by
3994 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
3996 For a simple proof, see Appendix \ref{sec.math_hooke}.
3998 <<hooke tension function declaration>>=
3999 extern double hooke_handler(double x, void *pdata);
4000 extern double inverse_hooke_handler(double F, void *pdata);
4003 First we define a function that computes the effective spring constant
4004 (stored in a single [[hooke_param_t]]) for the entire group.
4005 <<internal hooke functions>>=
4006 static void hooke_param_grouper(tension_handler_data_t *data,
4007 hooke_param_t *param) {
4008 list_t *list = data->group_tension_params;
4012 assert(list != NULL); /* empty active group?! */
4013 while (list != NULL) {
4014 assert( ((hooke_param_t *)list->d)->k > 0 );
4015 k += 1.0/ ((hooke_param_t *)list->d)->k;
4017 fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
4018 __FUNCTION__, ((hooke_param_t *)list->d)->k);
4024 fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N (from data %p)\n",
4025 __FUNCTION__, k, x, k*x, data);
4032 Everyone knows Hooke's law: $F=kx$.
4034 double hooke_handler(double x, void *pdata)
4036 hooke_param_t param = {0};
4039 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
4045 Inverting Hooke's law gives $x=F/k$.
4046 <<inverse hooke handler>>=
4047 double inverse_hooke_handler(double F, void *pdata)
4049 hooke_param_t param = {0};
4052 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
4058 Not much to keep track of for a single effective spring, just the
4059 spring constant $k$.
4060 <<hooke structure>>=
4061 typedef struct hooke_param_struct {
4062 double k; /* spring constant in N/m */
4067 We wrap up our Hooke implementation with some book-keeping functions.
4068 <<hooke tension parameter create/destroy function declarations>>=
4069 extern void *string_create_hooke_param_t(char **param_strings);
4070 extern void destroy_hooke_param_t(void *p);
4074 <<hooke parameter create/destroy functions>>=
4075 hooke_param_t *create_hooke_param_t(double k)
4077 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
4078 assert(ret != NULL);
4083 void *string_create_hooke_param_t(char **param_strings)
4085 return create_hooke_param_t(safe_strtod(param_strings[0], __FUNCTION__));
4088 void destroy_hooke_param_t(void *p)
4095 <<hooke tension model global declarations>>=
4096 extern int num_hooke_params;
4097 extern char *hooke_param_descriptions[];
4098 extern char hooke_param_string[];
4101 <<hooke tension model globals>>=
4102 int num_hooke_params = 1;
4103 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
4104 char hooke_param_string[]="0.05";
4107 <<hooke tension model getopt>>=
4108 {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}
4111 \subsection{Worm-like chain}
4114 We can model several unfolded domains with a single worm-like chain.
4115 <<worm-like chain functions>>=
4116 <<internal worm-like chain functions>>
4117 <<worm-like chain handler>>
4118 <<inverse worm-like chain handler>>
4119 <<worm-like chain create/destroy functions>>
4122 <<worm-like chain tension model declarations>>=
4124 <<worm-like chain handler declaration>>
4125 <<inverse worm-like chain handler declaration>>
4126 <<worm-like chain create/destroy function declarations>>
4127 <<worm-like chain tension model global declarations>>
4131 <<internal worm-like chain functions>>=
4132 <<worm-like chain function>>
4133 <<inverse worm-like chain function>>
4134 <<worm-like chain parameter grouper>>
4137 The combination of all unfolded domains is modeled as a worm like chain,
4138 with the force $F_{\text{WLC}}$ approximately given by
4140 F_{\text{WLC}} \approx \frac{k_B T}{p}
4141 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
4143 where $x$ is the distance between the fixed ends,
4144 $k_B$ is Boltzmann's constant,
4145 $T$ is the absolute temperature,
4146 $p$ is the persistence length, and
4147 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
4148 <<worm-like chain function>>=
4149 static double wlc(double x, double T, double p, double L)
4151 double xL=0.0; /* uses global k_B */
4153 assert(T > 0); assert(p > 0); assert(L > 0);
4154 if (x >= L) return HUGE_VAL;
4155 xL = x/L; /* unitless */
4156 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
4157 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
4158 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
4163 Inverting the approximate WLC is fairly straightforward, if a bit tedious.
4164 Letting $x_L \equiv x/L$ and $F_T \equiv Fp/k_B T$ we have
4166 F_T &\approx \frac{1}{4} \p({\frac{1}{(1-x_L)^2} - 1}) + x_L \\
4167 F_T (1-x_L)^2 &\approx \frac{1}{4} - \frac{1}{4} (1-x_L)^2 + x_L(1-x_L)^2 \\
4168 0 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) (1-x_L)^2 + x_L(1-x_L)^2 \\
4169 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) \p({1-2x_L+x_L^2})
4170 + x_L - 2x_L^2 + x_L^3 \\
4171 &\approx \p[{\frac{1}{4} - \p({F_T + \frac{1}{4}})}]
4172 + x_L \p[{2F_T + \frac{1}{2} + 1}]
4173 + x_L^2 \p[{-\p({F_T + \frac{1}{4}}) - 2}]
4176 + x_L \p({2F_T + \frac{3}{2}})
4177 - x_L^2 \p({F_T + \frac{9}{4}})
4179 &\approx ax_L^3 + bx_L^2 + cx_L + d
4181 where $a \equiv 1$, $b \equiv -F_T - 9/4$, $c \equiv 2F_T + 3/2$, and
4183 % From \citet{wikipedia_cubic_function} the discriminant
4185 % \Delta &= -4b^3d + b^2c^2 - 4ac^3 + 18abcd - 27a^2d^2 \\
4186 % &= -4F_T\p({F_T + \frac{9}{4}})^3
4187 % + \p({F_T + \frac{9}{4}})^2 \cdot \p({2F_T + \frac{3}{2}})^2
4188 % - 4 \p({2F_T + \frac{3}{2}})^3
4189 % + 18F_T \p({F_T + \frac{9}{4}}) \cdot \p({2F_T + \frac{3}{2}})
4191 % &= -4F_T \p({F_T^3 + \frac{27}{4}F_T^2 +\frac{243}{16}F_T +\frac{729}{64}})
4192 %% a^3 + 3a^2b + 3ab^2 + b^3
4193 % + \p({F_T^2 + \frac{9}{2}F_T + \frac{81}{16}})
4194 % \cdot \p({4F_T^2 + 6F_T + \frac{9}{4}}) \\
4195 % &\qquad - 4 \p({8F_T^3 + 18F_T^2 + \frac{27}{2}F_T + \frac{27}{8}})
4196 %% a^3 + 3a^2b + 3ab^2 + b^3
4197 % + 18F_T \p({2F_T^2 + 6F_T + \frac{27}{8}})
4199 % &= -4F_T^4 - 27F_T^3 - \frac{243}{4}F_T^2 - \frac{729}{16}F_T
4200 % + 4F_T^4 + 24F_T^3 +\frac{99}{2}F_T^2 +\frac{81}{2}F_T +\frac{729}{64}\\
4201 % &\qquad - 32F_T^3 - 72F_T^2 - 54F_T - \frac{27}{2}
4202 % + 36F_T^3 + 108F_T^2 + \frac{243}{4}F_T
4204 % &= F_T^4 \p({-4 + 4}) + F_T^3 \p({-27 + 24 - 32 + 36})
4205 % + F_T^2 \p({-\frac{243}{4} + \frac{99}{2} - 72 + 108 -27}) \\
4206 % &\qquad + F_T \p({-\frac{729}{16} + \frac{81}{2} - 54 + \frac{243}{4}})
4207 % \p({\frac{729}{64} - \frac{27}{2}}) \\
4208 % &= F_T^3 - \frac{9}{4}F_T^2 + \frac{27}{16}F_T - \frac{135}{64}
4210 We can plug these coefficients into the GSL cubic solver to invert the
4211 WLC\citep{galassi05}. The applicable root is always the one which
4212 comes before the singularity, which will be the smallest real root.
4213 <<inverse worm-like chain function>>=
4214 static double inverse_wlc(double F, double T, double p, double L)
4216 double FT = F*p/(k_B*T); /* uses global k_B */
4217 double xL0, xL1, xL2;
4220 assert(T > 0); assert(p > 0); assert(L > 0);
4223 num_roots = gsl_poly_solve_cubic(-(FT+2.25),2*FT+1.5,-FT, &xL0, &xL1, &xL2);
4224 assert(num_roots > 0); assert(xL0 >= -DOUBLE_PRECISION); assert(xL0 < 1);
4225 if (xL0 < 0) xL0 = 0.0;
4231 First we define a function that computes the effective WLC parameters
4232 (stored in a single [[wlc_param_t]]) for the entire group.
4233 <<worm-like chain parameter grouper>>=
4234 static void wlc_param_grouper(tension_handler_data_t *data,
4235 wlc_param_t *param) {
4236 list_t *list = data->group_tension_params;
4237 double p=0.0, L=0.0;
4240 assert(list != NULL); /* empty active group?! */
4241 p = ((wlc_param_t *)list->d)->p;
4242 while (list != NULL) {
4243 /* enforce consistent p */
4244 assert( ((wlc_param_t *)list->d)->p == p);
4245 L += ((wlc_param_t *)list->d)->L;
4247 fprintf(stderr, "%s: WLC section %g m long in series\n",
4248 __FUNCTION__, ((wlc_param_t *)list->d)->L);
4253 fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
4254 __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
4262 <<worm-like chain handler declaration>>=
4263 extern double wlc_handler(double x, void *pdata);
4266 This model requires that each unfolded domain in the group have the
4267 same persistence length.
4268 <<worm-like chain handler>>=
4269 double wlc_handler(double x, void *pdata)
4271 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4272 wlc_param_t param = {0};
4275 wlc_param_grouper(data, ¶m);
4276 return wlc(x, data->env->T, param.p, param.L);
4281 <<inverse worm-like chain handler declaration>>=
4282 extern double inverse_wlc_handler(double F, void *pdata);
4285 <<inverse worm-like chain handler>>=
4286 double inverse_wlc_handler(double F, void *pdata)
4288 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4289 wlc_param_t param = {0};
4291 wlc_param_grouper(data, ¶m);
4292 return inverse_wlc(F, data->env->T, param.p, param.L);
4297 <<worm-like chain structure>>=
4298 typedef struct wlc_param_struct {
4299 double p; /* WLC persistence length (m) */
4300 double L; /* WLC contour length (m) */
4304 <<worm-like chain create/destroy function declarations>>=
4305 extern void *string_create_wlc_param_t(char **param_strings);
4306 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
4309 <<worm-like chain create/destroy functions>>=
4310 wlc_param_t *create_wlc_param_t(double p, double L)
4312 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
4313 assert(ret != NULL);
4319 void *string_create_wlc_param_t(char **param_strings)
4321 return create_wlc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4322 safe_strtod(param_strings[1], __FUNCTION__));
4325 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
4333 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.
4334 TODO (now we copy the string before parsing, but still do this for clarity).
4336 <<valid string write code>>=
4337 char string[] = "abc";
4340 <<valid string write code>>=
4341 char *string = "abc";
4345 <<worm-like chain tension model global declarations>>=
4346 extern int num_wlc_params;
4347 extern char *wlc_param_descriptions[];
4348 extern char wlc_param_string[];
4350 <<worm-like chain tension model globals>>=
4351 int num_wlc_params = 2;
4352 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
4353 char wlc_param_string[]="0.39e-9,28.4e-9";
4356 <<worm-like chain tension model getopt>>=
4357 {&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}
4359 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
4361 \subsection{Freely-jointed chain}
4364 <<freely-jointed chain functions>>=
4365 <<freely-jointed chain function>>
4366 <<freely-jointed chain handler>>
4367 <<freely-jointed chain create/destroy functions>>
4370 <<freely-jointed chain tension model declarations>>=
4371 <<freely-jointed chain handler declaration>>
4372 <<freely-jointed chain create/destroy function declarations>>
4373 <<freely-jointed chain tension model global declarations>>
4377 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
4378 The tension of a chain of $N$ such links is given by
4380 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
4382 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}.
4383 <<freely-jointed chain function>>=
4384 double langevin(double x, void *params)
4387 return 1.0/tanh(x) - 1/x;
4389 <<one dimensional bracket declaration>>
4390 <<one dimensional solve declaration>>
4391 double inv_langevin(double x)
4393 double lb=0.0, ub=1.0, ret;
4394 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
4395 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
4399 double fjc(double x, double T, double l, int N)
4401 double L = l*(double)N;
4402 if (x == 0) return 0;
4403 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
4404 return k_B*T/l * inv_langevin(x/L);
4408 <<freely-jointed chain handler declaration>>=
4409 extern double fjc_handler(double x, void *pdata);
4411 <<freely-jointed chain handler>>=
4412 double fjc_handler(double x, void *pdata)
4414 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4415 list_t *list = data->group_tension_params;
4420 assert(list != NULL); /* empty active group?! */
4421 l = ((fjc_param_t *)list->d)->l;
4422 while (list != NULL) {
4423 /* enforce consistent l */
4424 assert( ((fjc_param_t *)list->d)->l == l);
4425 N += ((fjc_param_t *)list->d)->N;
4427 fprintf(stderr, "%s: FJC section %d links long in series\n",
4428 __FUNCTION__, ((fjc_param_t *)list->d)->N);
4433 fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
4434 __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
4436 return fjc(x, data->env->T, l, N);
4440 <<freely-jointed chain structure>>=
4441 typedef struct fjc_param_struct {
4442 double l; /* FJC link length (m) */
4443 int N; /* FJC number of links */
4447 <<freely-jointed chain create/destroy function declarations>>=
4448 extern void *string_create_fjc_param_t(char **param_strings);
4449 extern void destroy_fjc_param_t(void *p);
4452 <<freely-jointed chain create/destroy functions>>=
4453 fjc_param_t *create_fjc_param_t(double l, double N)
4455 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
4456 assert(ret != NULL);
4462 void *string_create_fjc_param_t(char **param_strings)
4464 return create_fjc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4465 safe_strtod(param_strings[1], __FUNCTION__));
4468 void destroy_fjc_param_t(void *p)
4475 <<freely-jointed chain tension model global declarations>>=
4476 extern int num_fjc_params;
4477 extern char *fjc_param_descriptions[];
4478 extern char fjc_param_string[];
4481 <<freely-jointed chain tension model globals>>=
4482 int num_fjc_params = 2;
4483 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
4484 char fjc_param_string[]="0.5e-9,200";
4487 <<freely-jointed chain tension model getopt>>=
4488 {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}
4492 \label{sec.piston_tension}
4494 An infinitely stretchable domain with no tension for extensions less
4495 than a particular threshold $L$, and infinite tension for $x>L$. The
4496 tension at $x=L$ is undefined, but may be any positive value. The
4497 model is called the ``piston'' model because it resembles a
4498 frictionless piston sliding in a rigid cylinder of length $L$.
4501 0 & \text{if $x<L$}, \\
4502 \infty & \text{if $x>L$}.
4506 Note that because of it's infinte stiffness at $x=L$, fully extended
4507 domains with this tension model will be identical to completely rigid
4508 domains. However, there is a distinction between this model and rigid
4509 domains for $x<L$, because any reactions that occur at $F=0$
4510 (e.g. refolding) will have more time to occur.
4512 Because the tension at $x=L$ is undefined, the user must make sure
4513 domains with this tension model are never the initial active group,
4514 because it would break [[tension_balance()]] in [[find_tension()]]
4515 (see Section \ref{sec.tension_balance}).
4517 <<piston tension functions>>=
4518 <<piston tension parameter grouper>>
4519 <<piston tension handler>>
4520 <<inverse piston tension handler>>
4521 <<piston tension parameter create/destroy functions>>
4524 <<piston tension model declarations>>=
4525 <<piston tension handler declaration>>
4526 <<inverse piston tension handler declaration>>
4527 <<piston tension parameter create/destroy function declarations>>
4528 <<piston tension model global declarations>>
4532 First we define a function that computes the effective piston parameters
4533 (stored in a single [[piston_tension_param_t]]) for the entire group.
4534 <<piston tension parameter grouper>>=
4535 static void piston_tension_param_grouper(tension_handler_data_t *data,
4536 piston_tension_param_t *param) {
4537 list_t *list = data->group_tension_params;
4541 assert(list != NULL); /* empty active group?! */
4542 while (list != NULL) {
4543 L += ((piston_tension_param_t *)list->d)->L;
4549 <<piston tension handler declaration>>=
4550 extern double piston_tension_handler(double x, void *pdata);
4552 <<piston tension handler>>=
4553 double piston_tension_handler(double x, void *pdata)
4555 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4556 piston_tension_param_t param = {0};
4559 piston_tension_param_grouper(data, ¶m);
4560 if (x <= param.L) F = 0;
4563 fprintf(stderr, "%s: x %g, L %g, tension %g N\n", __FUNCTION__, x, param.L, F);
4570 We abuse [[x_of_xo()]] and [[oneD_solve()]] with this inverse
4571 definition (see Section \ref{sec.tension_balance}), because the
4572 $x(F=0)<=L$, but our function returns $x(F=0)=0$. This is fine when
4573 the total extension \emph{is} zero, but cheats a bit when there is
4574 some total extension. For any non-zero extension, the initial active
4575 group produces some tension (we hope. True for all our other tension
4576 models). This tension fully extends the piston group (of which there
4577 should be only one, since all piston states can get lumped into the
4578 same tension group.). If the total extension is $\le$ the target
4579 extension, the full extension is accurate, so the inverse was valid
4580 after all. If not, [[oneD_solve()]] will halve the extension of the
4581 first group, reducing the overall tension. For total extension less
4582 than the piston extension, this bisection continues for [[max_steps]],
4583 leaving $x_0$ reduced by a factor of $2^\text{max\_steps}$, a very
4584 small number. Because of this, the returned tension $F_0(x_0)$ is
4585 very small, even though the total extension found by [[x_of_xo()]]
4586 is still $>$ the piston length.
4588 While this works, it is clearly not the most elegant or robust
4589 solution. TODO: think of (and implement) a better procedure.
4591 <<inverse piston tension handler declaration>>=
4592 extern double inverse_piston_tension_handler(double x, void *pdata);
4594 <<inverse piston tension handler>>=
4595 double inverse_piston_tension_handler(double F, void *pdata)
4597 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4598 piston_tension_param_t param = {0};
4601 piston_tension_param_grouper(data, ¶m);
4602 if (F == 0.0) return 0;
4608 <<piston tension structure>>=
4609 typedef struct piston_tension_param_struct {
4610 double L; /* length of piston in m */
4611 } piston_tension_param_t;
4615 <<piston tension parameter create/destroy function declarations>>=
4616 extern void *string_create_piston_tension_param_t(char **param_strings);
4617 extern void destroy_piston_tension_param_t(void *p);
4621 <<piston tension parameter create/destroy functions>>=
4622 piston_tension_param_t *create_piston_tension_param_t(double L)
4624 piston_tension_param_t *ret
4625 = (piston_tension_param_t *) malloc(sizeof(piston_tension_param_t));
4626 assert(ret != NULL);
4631 void *string_create_piston_tension_param_t(char **param_strings)
4633 return create_piston_tension_param_t(
4634 safe_strtod(param_strings[0],__FUNCTION__));
4637 void destroy_piston_tension_param_t(void *p)
4645 <<piston tension model global declarations>>=
4646 extern int num_piston_tension_params;
4647 extern char *piston_tension_param_descriptions[];
4648 extern char piston_tension_param_string[];
4652 <<piston tension model globals>>=
4653 int num_piston_tension_params = 1;
4654 char *piston_tension_param_descriptions[] = {"length L, m"};
4655 char piston_tension_param_string[] = "0";
4659 <<piston tension model getopt>>=
4660 {&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}
4663 \subsection{Tension model implementation}
4665 <<tension-model.c>>=
4668 <<tension model includes>>
4669 #include "tension_model.h"
4670 <<tension model internal definitions>>
4671 <<tension model globals>>
4672 <<tension model functions>>
4675 <<tension model includes>>=
4676 #include <assert.h> /* assert() */
4677 #include <stdlib.h> /* NULL, strto*() */
4678 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4679 #include <stdio.h> /* fprintf(), stdout */
4680 #include <errno.h> /* errno, ERANGE */
4683 #include "tension_balance.h" /* oneD_*() */
4686 <<tension model internal definitions>>=
4687 <<constant tension structure>>
4689 <<worm-like chain structure>>
4690 <<freely-jointed chain structure>>
4691 <<piston tension structure>>
4694 <<tension model globals>>=
4695 <<tension handler array global>>
4696 <<constant tension model globals>>
4697 <<hooke tension model globals>>
4698 <<worm-like chain tension model globals>>
4699 <<freely-jointed chain tension model globals>>
4700 <<piston tension model globals>>
4703 <<tension model functions>>=
4705 <<constant tension functions>>
4707 <<worm-like chain functions>>
4708 <<freely-jointed chain functions>>
4709 <<piston tension functions>>
4712 \subsection{Utilities}
4714 The tension models can get complicated, and users may want to reassure
4715 themselves that this portion of the input physics are functioning properly. The
4716 stand-alone program [[tension_model_utils]] demonstrates the tension model
4717 interface without the overhead of having to understand a full simulation.
4718 [[tension_model_utils]] takes command line model arguments like the full
4719 simulation, and prints $F(x)$ data to the screen for the selected model over a
4722 <<tension-model-utils.c>>=
4724 <<tension model utility includes>>
4725 <<tension model utility definitions>>
4726 <<tension model utility getopt functions>>
4728 <<tension model utility main>>
4731 <<tension model utility main>>=
4732 int main(int argc, char **argv)
4734 <<tension model getopt array definition>>
4735 tension_model_getopt_t *model = NULL;
4739 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
4740 tension_handler_data_t tdata;
4741 int num_param_args; /* for INIT_MODEL() */
4742 char **param_args; /* for INIT_MODEL() */
4744 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model,
4745 &Fmax, &Xmax, &flags);
4747 if (flags & VFLAG) {
4748 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4750 INIT_MODEL("utility", model, model->params, params);
4751 tdata.group_tension_params = NULL;
4752 push(&tdata.group_tension_params, params, 1);
4754 tdata.persist = NULL;
4755 if (model->handler == NULL) {
4756 printf("No tension function for model %s\n", model->name);
4762 printf("#Distance (m)\tForce (N)\n");
4763 for (i=0; i<=N; i++) {
4764 x = Xmax*i/(double)N;
4765 F = (*model->handler)(x, &tdata);
4766 if (F < 0 || F > Fmax) break;
4767 printf("%g\t%g\n", x, F);
4769 if (flags & VFLAG && i <= N) {
4770 /* explain exit condition */
4772 printf("Impossible force %g\n", F);
4774 printf("Reached large force limit %g > %g\n", F, Fmax);
4777 params = pop(&tdata.group_tension_params);
4778 if (model->destructor)
4779 (*model->destructor)(params);
4784 <<tension model utility includes>>=
4787 #include <string.h> /* strlen() */
4788 #include <assert.h> /* assert() */
4789 #include <errno.h> /* errno, ERANGE */
4793 #include "tension_model.h"
4796 <<tension model utility definitions>>=
4797 <<version definition>>
4798 #define VFLAG 1 /* verbose */
4799 <<string comparison definition and globals>>
4800 <<tension model getopt definitions>>
4801 <<initialize model definition>>
4805 <<tension model utility getopt functions>>=
4807 <<index tension model>>
4808 <<tension model utility help>>
4809 <<tension model utility get options>>
4812 <<tension model utility help>>=
4813 void help(char *prog_name,
4815 int n_tension_models, tension_model_getopt_t *tension_models,
4816 int tension_model, double Fmax, double Xmax)
4819 printf("usage: %s [options]\n", prog_name);
4820 printf("Version %s\n\n", VERSION);
4821 printf("A utility for understanding the available tension models\n\n");
4822 printf("Environmental options:\n");
4823 printf("-T T\tTemperature T (currently %g K)\n", env->T);
4824 printf("-C T\tYou can also set the temperature T in Celsius\n");
4825 printf("Model options:\n");
4826 printf("-m model\tFolded tension model (currently %s)\n",
4827 tension_models[tension_model].name);
4828 printf("-a args\tFolded tension model argument string (currently %s)\n",
4829 tension_models[tension_model].params);
4830 printf("F(x) is calculated for a range of x and printed\n");
4831 printf("For example:\n");
4832 printf(" #Distance (m)\tForce (N)\n");
4833 printf(" 123.456\t7.89\n");
4835 printf("-F\tSet the maximum F value for the standard mode F(x) output (currently %g)\n", Fmax);
4836 printf("-X\tSet the maximum x value for the standard mode F(x) output (currently %g)\n", Xmax);
4837 printf("-V\tChange output to verbose mode\n");
4838 printf("-h\tPrint this help and exit\n");
4840 printf("Tension models:\n");
4841 for (i=0; i<n_tension_models; i++) {
4842 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
4843 for (j=0; j < tension_models[i].num_params ; j++)
4844 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
4845 printf("\t\tdefaults: %s\n", tension_models[i].params);
4850 <<tension model utility get options>>=
4851 void get_options(int argc, char **argv, environment_t *env,
4852 int n_tension_models, tension_model_getopt_t *tension_models,
4853 tension_model_getopt_t **model, double *Fmax, double *Xmax,
4854 unsigned int *flags)
4856 char *prog_name = NULL;
4857 char c, options[] = "T:C:m:a:F:X:Vh";
4858 int tension_model=0, num_strings;
4859 extern char *optarg;
4860 extern int optind, optopt, opterr;
4864 /* setup defaults */
4866 prog_name = argv[0];
4867 env->T = 300.0; /* K */
4868 *Fmax = 1e5; /* N */
4869 *Xmax = 1e-6; /* m */
4871 *model = tension_models;
4873 while ((c=getopt(argc, argv, options)) != -1) {
4875 case 'T': env->T = safe_strtod(optarg, "-T"); break;
4876 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
4878 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
4879 *model = tension_models+tension_model;
4882 tension_models[tension_model].params = optarg;
4884 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
4885 case 'X': *Xmax = safe_strtod(optarg, "-X"); break;
4886 case 'V': *flags |= VFLAG; break;
4888 fprintf(stderr, "unrecognized option '%c'\n", optopt);
4889 /* fall through to default case */
4891 help(prog_name, env, n_tension_models, tension_models, tension_model, *Fmax, *Xmax);
4900 \section{Unfolding rate models}
4901 \label{sec.k_models}
4903 We have two main choices for determining $k$: Bell's model, which assumes the
4904 folded domains exist in a single `folded' state and make instantaneous,
4905 stochastic transitions over a free energy barrier; and Kramers' model, which
4906 treats the unfolding as a thermalized diffusion process.
4907 We deal with these options like we dealt with the different tension models: by bundling all model-specific
4908 parameters into structures, with model specific functions for determining $k$.
4910 <<k func definition>>=
4911 typedef double k_func_t(double F, environment_t *env, void *params);
4914 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
4915 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]].
4917 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
4918 because \LaTeX\ doesn't like underscores outside math mode.
4919 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
4920 You could use spaces instead (`\verb+ +'),
4921 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
4922 which I don't like as much.
4928 <<k func definition>>
4929 <<null k declarations>>
4930 <<const k declarations>>
4931 <<bell k declarations>>
4932 <<kbell k declarations>>
4933 <<kramers k declarations>>
4934 <<kramers integ k declarations>>
4935 #endif /* K_MODEL_H */
4938 <<k model module makefile lines>>=
4939 NW_SAWSIM_MODS += k_model
4940 CHECK_BINS += check_k_model
4941 check_k_model_MODS = parse string_eval
4943 [[check_k_model]] also depends on the parse module.
4945 <<k model utils makefile lines>>=
4946 K_MODEL_MODS = k_model parse string_eval
4947 K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
4948 $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h)
4949 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
4950 notangle -Rk-model-utils.c $< > $@
4951 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
4952 gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
4953 clean_k_model_utils :
4954 rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
4960 For domains that are never folded.
4962 <<null k declarations>>=
4964 <<null k functions>>=
4969 <<null k model getopt>>=
4970 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
4973 \subsection{Constant rate model}
4976 This is a pretty boring model, but it is usefull for testing the engine.
4980 where $k_0$ is the constant unfolding rate.
4982 <<const k functions>>=
4983 <<const k function>>
4984 <<const k structure create/destroy functions>>
4987 <<const k declarations>>=
4988 <<const k function declaration>>
4989 <<const k structure create/destroy function declarations>>
4990 <<const k global declarations>>
4994 <<const k structure>>=
4995 typedef struct const_k_param_struct {
4996 double knot; /* transition rate k_0 (frac population per s) */
5000 <<const k function declaration>>=
5001 double const_k(double F, environment_t *env, void *const_k_params);
5003 <<const k function>>=
5004 double const_k(double F, environment_t *env, void *const_k_params)
5005 { /* Returns the rate constant k in frac pop per s. */
5006 const_k_param_t *p = (const_k_param_t *) const_k_params;
5008 assert(p->knot > 0);
5013 <<const k structure create/destroy function declarations>>=
5014 void *string_create_const_k_param_t(char **param_strings);
5015 void destroy_const_k_param_t(void *p);
5018 <<const k structure create/destroy functions>>=
5019 const_k_param_t *create_const_k_param_t(double knot)
5021 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
5022 assert(ret != NULL);
5027 void *string_create_const_k_param_t(char **param_strings)
5029 return create_const_k_param_t(safe_strtod(param_strings[0], __FUNCTION__));
5032 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
5039 <<const k global declarations>>=
5040 extern int num_const_k_params;
5041 extern char *const_k_param_descriptions[];
5042 extern char const_k_param_string[];
5045 <<const k globals>>=
5046 int num_const_k_params = 1;
5047 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
5048 char const_k_param_string[]="1";
5051 <<const k model getopt>>=
5052 {"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}
5055 \subsection{Bell's model}
5058 Quantitatively, Bell's model gives $k$ as
5060 k = k_0 \cdot e^{F dx / k_B T} \;,
5062 where $k_0$ is the unforced unfolding rate,
5063 $F$ is the applied unfolding force,
5064 $dx$ is the distance to the transition state, and
5065 $k_B T$ is the thermal energy\citep{schlierf06}.
5067 <<bell k functions>>=
5069 <<bell k structure create/destroy functions>>
5072 <<bell k declarations>>=
5073 <<bell k function declaration>>
5074 <<bell k structure create/destroy function declarations>>
5075 <<bell k global declarations>>
5079 <<bell k structure>>=
5080 typedef struct bell_param_struct {
5081 double knot; /* transition rate k_0 (frac population per s) */
5082 double dx; /* `distance to transition state' (nm) */
5086 <<bell k function declaration>>=
5087 double bell_k(double F, environment_t *env, void *bell_params);
5089 <<bell k function>>=
5090 double bell_k(double F, environment_t *env, void *bell_params)
5091 { /* Returns the rate constant k in frac pop per s.
5092 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5093 * uses global k_B in J/K */
5094 bell_param_t *p = (bell_param_t *) bell_params;
5095 assert(F >= 0); assert(env->T > 0);
5097 assert(p->knot > 0); //assert(p->dx >= 0); /* p < 0 for refolding reactions, etc. */
5099 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
5100 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
5101 p->knot * exp(F*p->dx / (k_B*env->T) ));
5102 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
5104 return p->knot * exp(F*p->dx / (k_B*env->T) );
5108 <<bell k structure create/destroy function declarations>>=
5109 void *string_create_bell_param_t(char **param_strings);
5110 void destroy_bell_param_t(void *p);
5113 <<bell k structure create/destroy functions>>=
5114 bell_param_t *create_bell_param_t(double knot, double dx)
5116 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
5117 assert(ret != NULL);
5123 void *string_create_bell_param_t(char **param_strings)
5125 return create_bell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5126 safe_strtod(param_strings[1], __FUNCTION__));
5129 void destroy_bell_param_t(void *p /* bell_param_t * */)
5136 <<bell k global declarations>>=
5137 extern int num_bell_params;
5138 extern char *bell_param_descriptions[];
5139 extern char bell_param_string[];
5143 int num_bell_params = 2;
5144 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
5145 char bell_param_string[]="3.3e-4,0.25e-9";
5148 <<bell k model getopt>>=
5149 {"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}
5151 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
5152 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5155 \subsection{Linker stiffness corrected Bell model}
5158 Walton et. al showed that the Bell model constant force approximation
5159 doesn't hold for biotin-streptavadin binding, and I extended their
5160 results to I27 for the 2009 Biophysical Society Annual
5161 Meeting\cite{walton08,king09}. More details to follow when I get done
5162 with the conference\ldots
5164 We adjust Bell's model to
5166 k = k_0 \cdot e^{\frac{F dx - \frac{1}{2} \kappa dx^2}{k_B T}} \;,
5168 where $k_0$ is the unforced unfolding rate,
5169 $F$ is the applied unfolding force,
5170 $\kappa$ is the effective spring constant,
5171 $dx$ is the distance to the transition state, and
5172 $k_B T$ is the thermal energy\citep{schlierf06}.
5174 <<kbell k functions>>=
5175 <<kbell k function>>
5176 <<kbell k structure create/destroy functions>>
5179 <<kbell k declarations>>=
5180 <<kbell k function declaration>>
5181 <<kbell k structure create/destroy function declarations>>
5182 <<kbell k global declarations>>
5186 <<kbell k structure>>=
5187 typedef struct kbell_param_struct {
5188 double knot; /* transition rate k_0 (frac population per s) */
5189 double dx; /* `distance to transition state' (nm) */
5193 <<kbell k function declaration>>=
5194 double kbell_k(double F, environment_t *env, void *kbell_params);
5196 <<kbell k function>>=
5197 double kbell_k(double F, environment_t *env, void *kbell_params)
5198 { /* Returns the rate constant k in frac pop per s.
5199 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5200 * uses global k_B in J/K */
5201 kbell_param_t *p = (kbell_param_t *) kbell_params;
5202 assert(F >= 0); assert(env->T > 0); assert(env->stiffness >= 0);
5204 assert(p->knot > 0); assert(p->dx >= 0);
5205 return p->knot * exp((F*p->dx - 0.5*env->stiffness*(p->dx)*(p->dx))/(k_B*env->T) );
5209 <<kbell k structure create/destroy function declarations>>=
5210 void *string_create_kbell_param_t(char **param_strings);
5211 void destroy_kbell_param_t(void *p);
5214 <<kbell k structure create/destroy functions>>=
5215 kbell_param_t *create_kbell_param_t(double knot, double dx)
5217 kbell_param_t *ret = (kbell_param_t *) malloc(sizeof(kbell_param_t));
5218 assert(ret != NULL);
5224 void *string_create_kbell_param_t(char **param_strings)
5226 return create_kbell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5227 safe_strtod(param_strings[1], __FUNCTION__));
5230 void destroy_kbell_param_t(void *p /* kbell_param_t * */)
5237 <<kbell k global declarations>>=
5238 extern int num_kbell_params;
5239 extern char *kbell_param_descriptions[];
5240 extern char kbell_param_string[];
5243 <<kbell k globals>>=
5244 int num_kbell_params = 2;
5245 char *kbell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
5246 char kbell_param_string[]="3.3e-4,0.25e-9";
5249 <<kbell k model getopt>>=
5250 {"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}
5252 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
5253 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5256 \subsection{Kramer's model}
5259 Kramer's model gives $k$ as
5261 % \frac{1}{k} = \frac{1}{D}
5262 % \int_{x_\text{min}}^{x_\text{max}}
5263 % e^{\frac{-U_F(x)}{k_B T}}
5264 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5267 %where $D$ is the diffusion constant,
5268 %$U_F(x)$ is the free energy along the unfolding coordinate, and
5269 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
5270 % before the minimum and after the maximum, respectively \citep{schlierf06}.
5272 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
5274 where $D$ is the diffusion constant,
5276 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
5278 are length scales where
5279 $x_c(F)$ is the location of the bound state and
5280 $x_{ts}(F)$ is the location of the transition state,
5281 $E(F,x)$ is the free energy, and
5282 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
5283 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
5284 \citet{evans97} Eqn.~3).
5286 <<kramers k functions>>=
5287 <<kramers k function>>
5288 <<kramers k structure create/destroy functions>>
5291 <<kramers k declarations>>=
5292 <<kramers k function declaration>>
5293 <<kramers k structure create/destroy function declarations>>
5294 <<kramers k global declarations>>
5298 <<kramers k structure>>=
5299 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
5300 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
5302 typedef struct kramers_param_struct {
5303 double D; /* diffusion rate (um/s) */
5310 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
5311 //kramers_x_func_t *xts; /* function returning transition x (nm) */
5312 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
5313 //kramers_E_func_t *E; /* function returning E (J) */
5314 //double *E_params; /* parametrize the energy landscape */
5315 //int n_E_params; /* length of E_params array */
5319 <<kramers k function declaration>>=
5320 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
5321 extern double kramers_k(double F, environment_t *env, void *kramers_params);
5323 <<kramers k function>>=
5324 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
5326 static char input[160]={0};
5327 static char output[80]={0};
5328 /* setup the environment */
5329 fprintf(in, "F = %g; T = %g\n", F, T);
5330 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
5331 string_eval(in, out, input, 80, output);
5332 return safe_strtod(output, __FUNCTION__);
5335 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
5337 static char input[160]={0};
5338 static char output[80]={0};
5339 /* setup the environment */
5340 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
5341 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
5342 string_eval(in, out, input, 80, output);
5343 return safe_strtod(output, __FUNCTION__);
5346 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
5348 kramers_param_t *p = (kramers_param_t *) kramers_params;
5349 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
5352 double kramers_k(double F, environment_t *env, void *kramers_params)
5353 { /* Returns the rate constant k in frac pop per s.
5354 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5355 * uses global k_B in J/K */
5356 kramers_param_t *p = (kramers_param_t *) kramers_params;
5357 double kbT, xc, xts, lc, lts, Eb;
5358 assert(F >= 0); assert(env->T > 0);
5361 assert(p->xc != NULL); assert(p->xts != NULL);
5362 assert(p->ddEdxx != NULL); assert(p->E != NULL);
5363 assert(p->in != NULL); assert(p->out != NULL);
5365 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
5366 if (xc == -1.0) return -1.0; /* error (Too much force) */
5367 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
5368 if (xc == -1.0) return -1.0; /* error (Too much force) */
5369 lc = sqrt(2.0*M_PI*kbT /
5370 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
5371 lts = sqrt(-2.0*M_PI*kbT/
5372 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
5373 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
5374 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
5375 //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));
5376 return p->D/(lc*lts) * exp(-Eb/kbT);
5380 <<kramers k structure create/destroy function declarations>>=
5381 void *string_create_kramers_param_t(char **param_strings);
5382 void destroy_kramers_param_t(void *p);
5385 <<kramers k structure create/destroy functions>>=
5386 kramers_param_t *create_kramers_param_t(double D,
5387 char *xc_fn, char *xts_fn,
5388 char *E_fn, char *ddEdxx_fn,
5389 char *global_define_string)
5390 // kramers_x_func_t *xc, kramers_x_func_t *xts,
5391 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
5392 // double *E_params, int n_E_params)
5394 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
5395 assert(ret != NULL);
5396 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5397 assert(ret->xc != NULL);
5398 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5399 assert(ret->xts != NULL);
5400 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
5401 assert(ret->E != NULL);
5402 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
5403 assert(ret->ddEdxx != NULL);
5405 strcpy(ret->xc, xc_fn);
5406 strcpy(ret->xts, xts_fn);
5407 strcpy(ret->E, E_fn);
5408 strcpy(ret->ddEdxx, ddEdxx_fn);
5409 //ret->E_params = E_params;
5410 //ret->n_E_params = n_E_params;
5411 string_eval_setup(&ret->in, &ret->out);
5412 fprintf(ret->in, "kB = %g\n", k_B);
5413 fprintf(ret->in, "%s\n", global_define_string);
5417 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
5418 void *string_create_kramers_param_t(char **param_strings)
5420 return create_kramers_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5428 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
5430 kramers_param_t *pk = (kramers_param_t *)p;
5432 string_eval_teardown(&pk->in, &pk->out);
5434 // free(pk->E_params);
5440 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
5441 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.
5442 We follow \citet{shillcock98} and use
5444 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
5445 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
5448 where TODO, variable meanings.%\citep{shillcock98}.
5449 The first and second derivatives of $E(F,E_b,x,x_s)$ are
5451 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\\
5452 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
5454 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
5455 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
5456 The roots are therefor at
5458 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
5459 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
5460 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
5463 <<kramers k global declarations>>=
5464 extern int num_kramers_params;
5465 extern char *kramers_param_descriptions[];
5466 extern char kramers_param_string[];
5469 <<kramers k globals>>=
5470 int num_kramers_params = 6;
5471 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"};
5472 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)}";
5474 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
5475 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
5476 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
5477 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.
5478 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
5479 It works so long as [[val_1]] is not [[false]].
5481 <<kramers k model getopt>>=
5482 {"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}
5485 \subsection{Kramer's model (natural cubic splines)}
5486 \label{sec.kramers_integ}
5488 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
5489 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
5490 \citet{schlierf06} Eqn.~2)
5492 \frac{1}{k} = \frac{1}{D}
5493 \int_{x_\text{min}}^{x_\text{max}}
5494 e^{\frac{U_F(x)}{k_B T}}
5495 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5498 where $D$ is the diffusion constant,
5499 $U_F(x)$ is the free energy along the unfolding coordinate, and
5500 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
5501 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
5503 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
5504 interpolating between them with cubic splines.
5505 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
5507 <<kramers integ k functions>>=
5508 <<spline functions>>
5509 <<kramers integ A k functions>>
5510 <<kramers integ B k functions>>
5511 <<kramers integ k function>>
5512 <<kramers integ k structure create/destroy functions>>
5515 <<kramers integ k declarations>>=
5516 <<kramers integ k function declaration>>
5517 <<kramers integ k structure create/destroy function declarations>>
5518 <<kramers integ k global declarations>>
5522 <<kramers integ k structure>>=
5523 typedef double func_t(double x, void *params);
5524 typedef struct kramers_integ_param_struct {
5525 double D; /* diffusion rate (um/s) */
5526 double x_min; /* integration bounds */
5528 func_t *E; /* function returning E (J) */
5529 void *E_params; /* parametrize the energy landscape */
5530 destroy_data_func_t *destroy_E_params;
5532 double F; /* for passing into GSL evaluations */
5534 } kramers_integ_param_t;
5537 <<spline param structure>>=
5538 typedef struct spline_param_struct {
5539 int n_knots; /* natural cubic spline knots for U(x) */
5542 gsl_spline *spline; /* GSL spline parameters */
5543 gsl_interp_accel *acc; /* GSL spline acceleration data */
5547 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
5549 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5551 <<kramers integ A k functions>>=
5552 double kramers_integ_k_integrandA(double x, void *params)
5554 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5555 double U = (*p->E)(x, p->E_params);
5556 double Fx = p->F * x;
5557 double kbT = k_B * p->env->T;
5558 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
5559 return exp(-(U-Fx)/kbT);
5561 double kramers_integ_k_integralA(double x, void *params)
5563 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5565 double result, abserr;
5566 size_t neval; /* for qng */
5567 static gsl_integration_workspace *w = NULL;
5568 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5569 f.function = &kramers_integ_k_integrandA;
5571 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5572 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5573 //fprintf(stderr, "integralA = %g\n", result);
5579 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
5580 \int_{x_\text{min}}^{x_\text{max}}
5581 e^{\frac{U_F(x)}{k_B T}}
5582 \text{Integral}_A(x)
5585 <<kramers integ B k functions>>=
5586 double kramers_integ_k_integrandB(double x, void *params)
5588 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5589 double U = (*p->E)(x, p->E_params);
5590 double Fx = p->F * x;
5591 double kbT = k_B * p->env->T;
5592 double intA = kramers_integ_k_integralA(x, params);
5593 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
5594 return exp((U-Fx)/kbT)*intA;
5596 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
5598 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5600 double result, abserr;
5601 size_t neval; /* for qng */
5602 static gsl_integration_workspace *w = NULL;
5603 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5604 f.function = &kramers_integ_k_integrandB;
5608 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5609 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5610 //fprintf(stderr, "integralB = %g\n", result);
5615 With the integrals taken care of, evaluating $k$ is simply
5617 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
5619 <<kramers integ k function declaration>>=
5620 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
5622 <<kramers integ k function>>=
5623 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
5624 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
5625 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
5626 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
5630 <<kramers integ k structure create/destroy function declarations>>=
5631 void *string_create_kramers_integ_param_t(char **param_strings);
5632 void destroy_kramers_integ_param_t(void *p);
5635 <<kramers integ k structure create/destroy functions>>=
5636 kramers_integ_param_t *create_kramers_integ_param_t(double D,
5637 double xmin, double xmax, func_t *E,
5639 destroy_data_func_t *destroy_E_params)
5641 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
5642 assert(ret != NULL);
5647 ret->E_params = E_params;
5648 ret->destroy_E_params = destroy_E_params;
5652 void *string_create_kramers_integ_param_t(char **param_strings)
5654 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
5655 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
5656 void *E_params = string_create_spline_param_t(param_strings+1);
5657 return create_kramers_integ_param_t(
5658 safe_strtod(param_strings[0], __FUNCTION__),
5659 safe_strtod(param_strings[2], __FUNCTION__),
5660 safe_strtod(param_strings[3], __FUNCTION__),
5661 &spline_eval, E_params, destroy_spline_param_t);
5664 void destroy_kramers_integ_param_t(void *params)
5666 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5669 (*p->destroy_E_params)(p->E_params);
5675 Finally we have the GSL spline wrappers:
5676 <<spline functions>>=
5678 <<create/destroy spline>>
5681 <<spline function>>=
5682 double spline_eval(double x, void *spline_params)
5684 spline_param_t *p = (spline_param_t *)(spline_params);
5685 return gsl_spline_eval(p->spline, x, p->acc);
5689 <<create/destroy spline>>=
5690 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
5692 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
5693 assert(ret != NULL);
5694 ret->n_knots = n_knots;
5697 ret->acc = gsl_interp_accel_alloc();
5698 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
5699 gsl_spline_init(ret->spline, x, y, n_knots);
5703 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
5704 void *string_create_spline_param_t(char **param_strings)
5708 double *x=NULL, *y=NULL;
5709 /* break into ordered pairs */
5710 parse_list_string(param_strings[0], SEP, '(', ')',
5711 &num_knots, &knot_coords);
5712 x = (double *)malloc(sizeof(double)*num_knots);
5714 y = (double *)malloc(sizeof(double)*num_knots);
5716 for (i=0; i<num_knots; i++) {
5717 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
5718 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
5721 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
5723 free_parsed_list(num_knots, knot_coords);
5724 return create_spline_param_t(num_knots, x, y);
5727 void destroy_spline_param_t(void *params)
5729 spline_param_t *p = (spline_param_t *)params;
5732 gsl_spline_free(p->spline);
5734 gsl_interp_accel_free(p->acc);
5744 <<kramers integ k global declarations>>=
5745 extern int num_kramers_integ_params;
5746 extern char *kramers_integ_param_descriptions[];
5747 extern char kramers_integ_param_string[];
5750 <<kramers integ k globals>>=
5751 int num_kramers_integ_params = 4;
5752 char *kramers_integ_param_descriptions[] = {
5753 "diffusion D, m^2 / s",
5754 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
5755 "starting integration bound x_min, m",
5756 "ending integration bound x_max, m"};
5757 //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";
5758 //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";
5759 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";
5760 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
5761 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
5762 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
5765 <<kramers integ k model getopt>>=
5766 {"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}
5768 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
5769 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5771 \subsection{Unfolding model implementation}
5775 <<k model includes>>
5776 #include "k_model.h"
5777 <<k model internal definitions>>
5779 <<k model functions>>
5782 <<k model includes>>=
5783 #include <assert.h> /* assert() */
5784 #include <stdlib.h> /* NULL, malloc(), strto*() */
5785 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
5786 #include <stdio.h> /* fprintf(), stdout */
5787 #include <string.h> /* strlen(), strcpy() */
5788 #include <errno.h> /* errno, ERANGE */
5789 #include <gsl/gsl_integration.h>
5790 #include <gsl/gsl_spline.h>
5795 <<k model internal definitions>>=
5796 <<const k structure>>
5797 <<bell k structure>>
5798 <<kbell k structure>>
5799 <<kramers k structure>>
5800 <<kramers integ k structure>>
5801 <<spline param structure>>
5804 <<k model globals>>=
5809 <<kramers k globals>>
5810 <<kramers integ k globals>>
5813 <<k model functions>>=
5815 <<null k functions>>
5816 <<const k functions>>
5817 <<bell k functions>>
5818 <<kbell k functions>>
5819 <<kramers k functions>>
5820 <<kramers integ k functions>>
5823 \subsection{Unfolding model unit tests}
5825 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
5826 <<check-k-model.c>>=
5827 <<k model check includes>>
5828 <<check relative error>>
5830 <<k model test suite>>
5831 <<main check program>>
5834 <<k model check includes>>=
5835 #include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
5836 #include <stdio.h> /* sprintf() */
5837 #include <assert.h> /* assert() */
5838 #include <math.h> /* exp() */
5839 #include <errno.h> /* errno, ERANGE */
5842 #include "k_model.h"
5845 <<k model test suite>>=
5849 <<k model suite definition>>
5852 <<k model suite definition>>=
5853 Suite *test_suite (void)
5855 Suite *s = suite_create ("k model");
5856 <<const k test case defs>>
5857 <<bell k test case defs>>
5859 <<const k test case adds>>
5860 <<bell k test case adds>>
5865 \subsubsection{Constant}
5867 <<const k test case defs>>=
5868 TCase *tc_const_k = tcase_create("Constant k");
5871 <<const k test case adds>>=
5872 tcase_add_test(tc_const_k, test_const_k_create_destroy);
5873 tcase_add_test(tc_const_k, test_const_k_over_range);
5874 suite_add_tcase(s, tc_const_k);
5879 START_TEST(test_const_k_create_destroy)
5881 char *knot[] = {"1","2","3","4","5","6"};
5882 char *params[] = {knot[0]};
5885 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5886 params[0] = knot[i];
5887 p = string_create_const_k_param_t(params);
5888 destroy_const_k_param_t(p);
5893 START_TEST(test_const_k_over_range)
5895 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
5896 char *knot[] = {"1","2","3","4","5","6"};
5897 char *params[] = {knot[0]};
5900 char param_string[80];
5903 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5904 params[0] = knot[i];
5905 p = string_create_const_k_param_t(params);
5906 for ( F=Fm; F<FM; F+=dF ) {
5907 fail_unless(const_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
5908 "const_k(%g, %g, &{%s}) = %g != %s",
5909 F, T, knot[i], const_k(F, &env, p), knot[i]);
5911 destroy_const_k_param_t(p);
5917 \subsubsection{Bell}
5919 <<bell k test case defs>>=
5920 TCase *tc_bell = tcase_create("Bell's k");
5923 <<bell k test case adds>>=
5924 tcase_add_test(tc_bell, test_bell_k_create_destroy);
5925 tcase_add_test(tc_bell, test_bell_k_at_zero);
5926 tcase_add_test(tc_bell, test_bell_k_at_one);
5927 suite_add_tcase(s, tc_bell);
5931 START_TEST(test_bell_k_create_destroy)
5933 char *knot[] = {"1","2","3","4","5","6"};
5935 char *params[] = {knot[0], dx};
5938 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5939 params[0] = knot[i];
5940 p = string_create_bell_param_t(params);
5941 destroy_bell_param_t(p);
5946 START_TEST(test_bell_k_at_zero)
5948 double F=0.0, T=300.0;
5949 char *knot[] = {"1","2","3","4","5","6"};
5951 char *params[] = {knot[0], dx};
5954 char param_string[80];
5957 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5958 params[0] = knot[i];
5959 p = string_create_bell_param_t(params);
5960 fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
5961 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
5962 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
5963 destroy_bell_param_t(p);
5968 START_TEST(test_bell_k_at_one)
5971 char *knot[] = {"1","2","3","4","5","6"};
5973 char *params[] = {knot[0], dx};
5974 double F= k_B*T/safe_strtod(dx, __FUNCTION__);
5977 char param_string[80];
5980 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5981 params[0] = knot[i];
5982 p = string_create_bell_param_t(params);
5983 CHECK_ERR(1e-6, safe_strtod(knot[i], __FUNCTION__)*exp(1.0), bell_k(F, &env, p));
5984 //fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__)*exp(1.0),
5985 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
5986 // F, T, knot[i], dx, bell_k(F, &env, p), safe_strtod(knot[i], __FUNCTION__)*exp(1.0));
5987 destroy_bell_param_t(p);
5993 <<kramers k tests>>=
5996 <<kramers k test case def>>=
5999 <<kramers k test case add>>=
6002 <<k model function tests>>=
6003 <<check relative error>>
6005 START_TEST(test_something)
6007 double k=5, x=3, last_x=2.0, F;
6008 one_dim_fn_t *handlers[] = {&hooke};
6009 void *data[] = {&k};
6010 double xi[] = {0.0};
6012 int new_active_groups = 1;
6013 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
6014 fail_unless(F = k*x, NULL);
6019 \subsection{Utilities}
6021 The unfolding models can get complicated, and are sometimes hard to explain to others.
6022 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
6023 the overhead of having to understand a full simulation.
6024 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
6025 or special mode, where the operation depends on the specific model selected.
6027 <<k-model-utils.c>>=
6029 <<k model utility includes>>
6030 <<k model utility definitions>>
6031 <<k model utility getopt functions>>
6032 <<k model utility multi model E>>
6033 <<k model utility main>>
6036 <<k model utility main>>=
6037 int main(int argc, char **argv)
6039 <<k model getopt array definition>>
6040 k_model_getopt_t *model = NULL;
6045 int num_param_args; /* for INIT_MODEL() */
6046 char **param_args; /* for INIT_MODEL() */
6048 double special_xmin;
6049 double special_xmax;
6050 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
6051 &Fmax, &special_xmin, &special_xmax, &flags);
6052 if (flags & VFLAG) {
6053 printf("#initializing model %s with parameters %s\n", model->name, model->params);
6055 INIT_MODEL("utility", model, model->params, params);
6058 if (model->k == NULL) {
6059 printf("No k function for model %s\n", model->name);
6063 printf("Fmax = %g <= 0\n", Fmax);
6069 printf("#F (N)\tk (%% pop. per s)\n");
6070 for (i=0; i<=N; i++) {
6071 F = Fmax*i/(double)N;
6072 k = (*model->k)(F, &env, params);
6074 printf("%g\t%g\n", F, k);
6079 if (model->k == NULL || model->k == &bell_k) {
6080 printf("No special mode for model %s\n", model->name);
6083 if (special_xmax <= special_xmin) {
6084 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
6088 double Xrng=(special_xmax-special_xmin), x, E;
6090 printf("#x (m)\tE (J)\n");
6091 for (i=0; i<=N; i++) {
6092 x = special_xmin + Xrng*i/(double)N;
6093 E = multi_model_E(model->k, params, &env, x);
6094 printf("%g\t%g\n", x, E);
6101 if (model->destructor)
6102 (*model->destructor)(params);
6107 <<k model utility multi model E>>=
6108 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
6110 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
6112 if (k == kramers_integ_k)
6113 E = (*p->E)(x, p->E_params);
6115 E = kramers_E(0, env, model_params, x);
6121 <<k model utility includes>>=
6124 #include <string.h> /* strlen() */
6125 #include <assert.h> /* assert() */
6126 #include <errno.h> /* errno, ERANGE */
6129 #include "string_eval.h"
6130 #include "k_model.h"
6133 <<k model utility definitions>>=
6134 <<version definition>>
6135 #define VFLAG 1 /* verbose */
6136 enum mode_t {M_K_OF_F, M_SPECIAL};
6137 <<string comparison definition and globals>>
6138 <<k model getopt definitions>>
6139 <<initialize model definition>>
6140 <<kramers k structure>>
6141 <<kramers integ k structure>>
6145 <<k model utility getopt functions>>=
6148 <<k model utility help>>
6149 <<k model utility get options>>
6152 <<k model utility help>>=
6153 void help(char *prog_name,
6155 int n_k_models, k_model_getopt_t *k_models,
6156 int k_model, double Fmax, double special_xmin, double special_xmax)
6159 printf("usage: %s [options]\n", prog_name);
6160 printf("Version %s\n\n", VERSION);
6161 printf("A utility for understanding the available unfolding models\n\n");
6162 printf("Environmental options:\n");
6163 printf("-T T\tTemperature T (currently %g K)\n", env->T);
6164 printf("-C T\tYou can also set the temperature T in Celsius\n");
6165 printf("Model options:\n");
6166 printf("-k model\tTransition rate model (currently %s)\n",
6167 k_models[k_model].name);
6168 printf("-K args\tTransition rate model argument string (currently %s)\n",
6169 k_models[k_model].params);
6170 printf("There are two output modes. In standard mode, k(F) is printed\n");
6171 printf("For example:\n");
6172 printf(" #Force (N)\tk (% pop. per s)\n");
6173 printf(" 123.456\t7.89\n");
6175 printf("In special mode, the output depends on the model.\n");
6176 printf("For models defining an energy function E(x), that function is printed\n");
6177 printf(" #Distance (m)\tE (J)\n");
6178 printf(" 0.0012\t0.0034\n");
6180 printf("-m\tChange output to standard mode\n");
6181 printf("-M\tChange output to special mode\n");
6182 printf("-F\tSet the maximum F value for the standard mode k(F) output (currently %g)\n", Fmax);
6183 printf("-x\tSet the minimum x value for the special mode E(x) output (currently %g)\n", special_xmin);
6184 printf("-X\tSet the maximum x value for the special mode E(x) output (currently %g)\n", special_xmax);
6185 printf("-V\tChange output to verbose mode\n");
6186 printf("-h\tPrint this help and exit\n");
6188 printf("Unfolding rate models:\n");
6189 for (i=0; i<n_k_models; i++) {
6190 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
6191 for (j=0; j < k_models[i].num_params ; j++)
6192 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
6193 printf("\t\tdefaults: %s\n", k_models[i].params);
6198 <<k model utility get options>>=
6199 void get_options(int argc, char **argv, environment_t *env,
6200 int n_k_models, k_model_getopt_t *k_models,
6201 enum mode_t *mode, k_model_getopt_t **model,
6202 double *Fmax, double *special_xmin, double *special_xmax,
6203 unsigned int *flags)
6205 char *prog_name = NULL;
6206 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
6208 extern char *optarg;
6209 extern int optind, optopt, opterr;
6213 /* setup defaults */
6215 prog_name = argv[0];
6216 env->T = 300.0; /* K */
6220 *Fmax = 1e-9; /* N */
6221 *special_xmax = 1e-8;
6222 *special_xmin = 0.0;
6224 while ((c=getopt(argc, argv, options)) != -1) {
6226 case 'T': env->T = safe_strtod(optarg, "-T"); break;
6227 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
6229 k_model = index_k_model(n_k_models, k_models, optarg);
6230 *model = k_models+k_model;
6233 k_models[k_model].params = optarg;
6235 case 'm': *mode = M_K_OF_F; break;
6236 case 'M': *mode = M_SPECIAL; break;
6237 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
6238 case 'x': *special_xmin = safe_strtod(optarg, "-x"); break;
6239 case 'X': *special_xmax = safe_strtod(optarg, "-X"); break;
6240 case 'V': *flags |= VFLAG; break;
6242 fprintf(stderr, "unrecognized option '%c'\n", optopt);
6243 /* fall through to default case */
6245 help(prog_name, env, n_k_models, k_models, k_model, *Fmax, *special_xmin, *special_xmax);
6254 \section{\LaTeX\ documentation}
6256 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
6257 The comment blocks are extracted (with nicely formatted code blocks), using
6258 <<latex makefile lines>>=
6259 $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6260 noweave -latex -index -delay $< > $@
6261 $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR)
6265 <<latex makefile lines>>=
6266 $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
6268 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
6269 $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
6270 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
6271 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
6272 mv $(BUILD_DIR)/sawsim.pdf $@
6274 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
6275 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
6276 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
6278 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
6279 <<latex makefile lines>>=
6281 rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
6282 $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
6283 $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.toc \
6284 $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib
6285 clean_latex : semi-clean_latex
6286 rm -f $(DOC_DIR)/sawsim.pdf
6291 [[make]] is a common utility on *nix systems for managing dependencies while building software.
6292 In this case, we have several \emph{targets} that we'd like to build:
6293 the main simulation program \prog;
6294 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
6295 the documentation [[sawsim.pdf]];
6296 and the [[Makefile]] itself.
6297 Besides the generated files, there is the utility target
6298 [[clean]] for removing all generated files except [[Makefile]].
6300 % [[dist]] for generating a distributable tar file.
6302 Extract the makefile with
6303 `[[notangle -Rmakefile src/sawsim.nw | sed 's/ /\t/' > Makefile]]'.
6304 The sed is needed because notangle replaces tabs with eight spaces,
6305 but [[make]] needs tabs.
6307 # Decide what directories to put things in
6312 # Note: Cannot use, for example, `./build', because make eats the `./'
6313 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
6315 # Modules (X.c, X.h) defined in the noweb file
6318 # Modules defined outside the noweb file
6319 FREE_SAWSIM_MODS = interp tavl
6321 <<list module makefile lines>>
6322 <<tension balance module makefile lines>>
6323 <<tension model module makefile lines>>
6324 <<k model module makefile lines>>
6325 <<parse module makefile lines>>
6326 <<string eval module makefile lines>>
6327 <<accel k module makefile lines>>
6329 SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS)
6331 # Everything needed for sawsim under one roof. sawsim.c must come first
6332 SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \
6333 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h)
6334 # Libraries needed to compile sawsim
6335 LIBS = gsl gslcblas m
6336 CHECK_LIBS = gsl gslcblas m check
6337 # The non-check binaries generated
6338 BINS = sawsim tension_model_utils k_model_utils sawsim_profile
6341 # Define the major targets
6342 all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
6344 view : $(DOC_DIR)/sawsim.pdf
6346 profile : $(BIN_DIR)/sawsim_profile | $(BIN_DIR)
6347 $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \
6348 -mnull -Mwlc -A0.39e-9,28e-9 -F8
6349 gprof $(BIN_DIR)/sawsim_profile gmon.out > $@
6350 check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim
6351 $(SHELL) -e -c 'for B in $^; do ./$$B; done'
6352 clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \
6353 clean_tension_model_utils \
6354 clean_k_model_utils clean_latex clean_check_sawsim
6355 rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \
6356 $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \
6357 $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
6358 $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
6359 $(BUILD_DIR)/global.h ./gmon.out
6360 $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
6362 # Various builds of sawsim
6363 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
6364 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6365 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
6366 gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6367 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
6368 gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6370 # Create the directories
6371 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
6374 # Copy over the source external to sawsim
6375 # Note: Cannot use, for example, `./build', because make eats the `./'
6376 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
6378 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6382 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6386 ## Basic source generated with noweb
6387 # The central files sawsim.c and global.h...
6388 $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6390 $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6391 notangle -Rglobal.h $< > $@
6392 # ... and the modules
6393 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6394 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6395 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6396 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6397 $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6398 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6399 # Note: I use `_' as a space in my file names, but noweb doesn't like
6400 # that so we use `-' to name the noweb roots and substitute here.
6402 ## Building the unit-test programs
6404 $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
6405 notangle -Rchecks $< > $@
6406 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
6407 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
6408 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR)
6409 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
6410 clean_check_sawsim :
6411 rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
6412 # ... and the modules
6414 $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
6415 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \
6416 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \
6417 $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6418 $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6419 $$(BUILD_DIR)/global.h | $(BIN_DIR)
6420 gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
6421 $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
6423 # todo: clean up the required modules to
6424 $(CHECK_BINS:%=clean_%) :
6425 rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
6427 # Cleaning up the modules
6429 $(SAWSIM_MODS:%=clean_%) :
6430 rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h)
6432 <<tension model utils makefile lines>>
6433 <<k model utils makefile lines>>
6434 <<latex makefile lines>>
6436 Makefile : $(SRC_DIR)/sawsim.nw
6437 notangle -Rmakefile $< | sed 's/ /\t/' > $@
6439 This is fairly self-explanatory. We compile a dynamically linked
6440 version ([[sawsim]]) and a statically linked version
6441 ([[sawsim_static]]). The static version is about 9 times as big, but
6442 it runs without needing dynamic GSL libraries to link to, so it's a
6443 better format for distributing to a cluster for parallel evaluation.
6447 \subsection{Hookean springs in series}
6448 \label{sec.math_hooke}
6451 The effective spring constant for $n$ springs $k_i$ in series is given by
6453 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6459 F &= k_i x_i &&\forall i \le n \\
6460 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
6461 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
6462 F &= k_1 x_1 = k_\text{eff} x \\
6463 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
6464 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6469 \addcontentsline{toc}{section}{References}
6470 \bibliography{sawsim}