1 %%%%%%%%%%%%%%%%%%%%%%%%% Start LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%
2 % we have our own preamble,
3 % use `noweave -delay` to not print noweb's automatic one
4 % -*- mode: Noweb; noweb-code-mode: c-mode -*-
5 \documentclass[letterpaper, 11pt]{article}
8 \noweboptions{smallcode}
10 \usepackage{url} % needed for natbib for underscores in urls?
11 \usepackage[breaklinks,colorlinks]{hyperref} % for internal links, \phantomsection
12 % breaklinks breaks long links
13 % colorlinks colors link text instead of boxing it
14 \usepackage{amsmath} % for \text{}, \xrightarrow[sub]{super}
15 \usepackage{amsthm} % for theorem and proof environments
16 \newtheorem{theorem}{Theorem}
18 \usepackage{subfig} % Support side-by-side figures
19 \usepackage{pgf} % Fancy graphics
20 \usepackage{tikz} % A nice, inline PGF frontend
21 \usetikzlibrary{automata} % Graph-theory library
23 \usepackage{doc} % BibTeX
24 \usepackage[super,sort&compress]{natbib} % fancy citation extensions
25 % super selects citations in superscript mode
26 % sort&compress automatically sorts and compresses compound citations (\cite{a,b,...})
28 \usepackage[mathletters]{ucs} % use math fonts to allow Greek UTF-8 in the text
29 \usepackage[utf8x]{inputenc} % allow UTF-8 characters
31 %\bibliographystyle{ieeetr} % pick the bibliography style, short and sweet
32 %\bibliographystyle{plain} % pick the bibliography style, includes dates
33 \bibliographystyle{plainnat}
34 \defcitealias{sw:noweb}{{\tt noweb}}
35 \defcitealias{galassi05}{GSL}
36 \defcitealias{sw:check}{{\tt check}}
37 \defcitealias{sw:python}{Python}
42 \textheight 9.5in % leave a bit of extra room for page numbers (at top of page)
45 \setlength{\parindent}{0pt}
46 \setlength{\parskip}{5pt}
48 % For some reason, the \maketitle wants to go on it's own page
49 % so we'll just hardcode the following in our first page.
50 %\title{Sawsim: a sawtooth protein unfolding simulator}
51 %\author{W.~Trevor King}
54 \newcommand{\prog}{[[sawsim]]}
55 \newcommand{\lang}{[[C]]}
56 \newcommand{\p}[3]{\left#1 #2 \right#3} % e.g. \p({complicated expr})
57 \newcommand{\Th}{\ensuremath{\sup{\text{th}}}} % ordinal th
58 \newcommand{\U}[1]{\ensuremath{\text{ #1}}} % units: $1\U{m/s}$
60 % single spaced lists, from Alvin Alexander
61 % http://www.devdaily.com/blog/post/latex/control-line-spacing-in-itemize-enumerate-tags/
62 \newenvironment{packed_item}{
64 \setlength{\itemsep}{1pt}
65 \setlength{\parskip}{0pt}
66 \setlength{\parsep}{0pt}
70 \nwfilename{sawsim.nw}
75 Sawsim: a sawtooth protein unfolding simulator \\
83 %%%%%%%%%%%%%%%%%%%%%%%%%%%% End LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%%
85 \section{Introduction}
87 The unfolding of multi-globular protein strings is a stochastic
88 process. In the \prog\ program, we use Monte Carlo methods to
89 simulate this unfolding through an explicit time-stepping approach.
90 We develop a program in \lang\ to simulate probability of unfolding
91 under a constant extension velocity of a chain of globular domains.
93 In order to extract physical meaning from many mechanical unfolding
94 simulations, it is necessary to compare the experimental results
95 against the results of simulations using different models and
96 parameters. The model and parameters that best match the experimental
97 data provide the ‘best’ model for that experiment.
99 Previous mechanical unfolding studies have rolled their own
100 simulations for modelling their experiment, and there has been little
101 exchange and standardization of these simulations between groups.
102 The \prog\ program attempts to fill this gap, providing a flexible,
103 extensible, \emph{community} program, to make it easier for groups
104 to share the burden and increase the transparency of data analysis.
106 ‘Sawsim’ is blend of ‘sawtooth simulation’.
108 \subsection{Related work}
112 \subsection{About this document}
114 This document was written in \citetalias{sw:noweb} with code blocks in
115 \lang\ and comment blocks in \LaTeX.
117 \subsection{System Requirements and Dependencies}
120 \subsection{Change Log}
122 Version 0.0 used only the unfolded domain WLC to determine the tension
123 and had a fixed timestep $dt$.
125 Version 0.1 added adaptive timesteps, adjusting so that the unfolding
126 probability for a given domain was constant.
128 Version 0.2 added Kramers' model unfolding rate calculations, and a
129 switch to select Bell's or Kramers' model. This induced a major
130 rewrite, introducing the tension group abstraction, and a split to
131 multiple \lang\ source files. Also added a unit-test suites for
132 sanity checking, and switched to SI units throughout.
134 Version 0.3 added integrated Kramers' model unfolding rate
135 calculations. Also replaced some of my hand-coded routines with GNU
136 Scientific Library \citepalias{galassi05} calls.
138 Version 0.4 added Python evaluation for the saddle-point Kramers'
139 model unfolding rate calculations, so that the model functions could
140 be controlled from the command line. Also added interpolating lookup
141 tables to accelerate slow unfolding rate model evaluation, and command
142 line control of tension groups.
144 Version 0.5 added the stiffness environmental parameter and it's
145 associated unfolding models.
147 Version 0.6 generalized from two state unfolding to arbitrary
148 $N$-state rate reactions. This allows simulations of
149 unfolding-refolding, metastable intermediates, etc., but required
150 reorganizing the command line interface.
152 Version 0.7 added tension model inverses, which seems to reduce
153 computation time by about a factor of 2. I was expecting more, but
154 I'll take what I can get.
156 Version 0.8 fixed a minor bug in unfolding probability for
157 multi-domain groups. The probability of at least one domain unfolding
158 had been calculated as $P_N=NP_1$, rather than $P_N=1-(1-P_1)^N$.
159 However, for small $P$ the two are equivalent.
161 Version 0.9 added the [[-P]] option to sawsim, setting the target
162 probability for the per-state transition $P_N$, not the per-domain
165 Version 0.10 added the [[-F]] option to sawsim, setting a limit on the
166 force change $dF$ in a single timestep for continuous pulling
167 simulations. It also added the piston tension model (Section
168 \ref{sec.piston_tension}), and adjusted the stiffness calculations to
169 handle domains with undefined stiffness.
171 <<version definition>>=
172 #define VERSION "0.10"
178 sawsim - program for simulating single molecule mechanical unfolding.
179 Copyright (C) 2008-2009, William Trevor King
181 This program is free software; you can redistribute it and/or
182 modify it under the terms of the GNU General Public License as
183 published by the Free Software Foundation; either version 3 of the
184 License, or (at your option) any later version.
186 This program is distributed in the hope that it will be useful, but
187 WITHOUT ANY WARRANTY; without even the implied warranty of
188 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
189 See the GNU General Public License for more details.
191 You should have received a copy of the GNU General Public License
192 along with this program; if not, write to the Free Software
193 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
196 The author may be contacted at <wking@drexel.edu> on the Internet, or
197 write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
198 Philadelphia PA 19104, USA.
211 The unfolding system is stored as a chain of \emph{domains} (Figure
212 \ref{fig.domain_chain}). Domains can be folded, globular protein
213 domains, unfolded protein linkers, AFM cantilevers, or other
214 stretchable system components. The system is modeled as graph of
215 possible domain states with transitions between them (Figure
216 \ref{fig.domain_states}).
220 \subfloat[][]{\label{fig.domain_chain}
221 \begin{tikzpicture}[->,node distance=2cm,font=\footnotesize]
222 \tikzstyle{every state}=[fill,draw=red!50,very thick,fill=red!20]
223 \node[state] (A) {domain 1};
224 \node[state] (B) [below of=A] {domain 2};
225 \node[state] (C) [below of=B] {.~.~.};
226 \node[state] (D) [below of=C] {domain $N$};
227 \node (S) [below of=D] {Surface};
228 \node (E) [above of=A] {};
230 \path[-] (A) edge (B)
231 (B) edge node [right] {Tension} (C)
234 \path[->,green] (A) edge node [right,black] {Extension} (E);
238 \subfloat[][]{\label{fig.domain_states}
239 \begin{tikzpicture}[->,node distance=2.5cm,shorten <=1pt,shorten >=1pt,font=\footnotesize]
240 \tikzstyle{every state}=[fill,draw=blue!50,very thick,fill=blue!20]
241 \node[state] (A) {cantilever};
242 \node[state] (C) [below of=A] {transition};
243 \node[state] (B) [left of=C] {folded};
244 \node[state] (D) [right of=C] {unfolded};
246 \path (B) edge [bend left] node [above] {$k_1$} (C)
247 (C) edge [bend left] node [below] {$k_1'$} (B)
248 edge [bend left] node [above] {$k_2$} (D)
249 (D) edge [bend left] node [below] {$k_2'$} (C);
252 \caption{\subref{fig.domain_chain} Extending a chain of domains. One end
253 of the chain is fixed, while the other is extended at a constant
254 speed. The domains are coupled with rigid linkers, so the domains
255 themselves must stretch to accomodate the extension.
256 \subref{fig.domain_states} Each domain exists in a discrete state. At
257 each timestep, it may transition into another state following a
258 user-defined state matrix such as this one, showing a metastable
259 transition state and an explicit ``cantilever'' domain.}
263 Each domain contributes to the total tension, which depends on the
264 chain extension. The particular model for the tension as a function
265 of extension is handled generally with function pointers. So far the
266 following models have been implemented:
268 \item Null (Appendix \ref{sec.null_tension}),
269 \item Constant (Appendix \ref{sec.const_tension}),
270 \item Hookean spring (Appendix \ref{sec.hooke}),
271 \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
272 \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
273 \item Piston (Appendix \ref{sec.piston_tension}),
276 The tension model and parameters used for a particular domain depend
277 on the domain's current state. The transition rates between states
278 are also handled generally with function pointers, with
281 \item Null (Appendix \ref{sec.null_k}),
282 \item Bell's model (Appendix \ref{sec.bell}),
283 \item Spring constant corrected Bell model (Appendix \ref{sec.kbell}),
284 \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
285 \item Kramers' saddle point model (Appendix \ref{sec.kramers}).
288 The unfolding of the chain domains is modeled in two stages. First
289 the tension of the chain is calculated. Then the tension (assumed to
290 be constant over the length of the chain, see Section
291 \ref{sec.timescales}), is applied to each domain, and used to
292 calculate the probability of that domain changing state in the next
293 timestep $dt$. Then the time is stepped forward, and the process is
294 repeated. The simulation is complete when a pre-set time limit
295 [[t_max]] is reached or a pre-selected state [[stop_state]] is empty.
297 int main(int argc, char **argv)
299 <<tension model getopt array definition>>
300 <<k model getopt array definition>>
302 /* variables initialized in get_options() */
303 list_t *state_list=NULL, *transition_list=NULL;
304 environment_t env={0};
305 double P, dF, t_max, dt_max, v, x_step;
306 state_t *stop_state=NULL;
308 /* variables used in the simulation loop */
309 double t=0, x=0, dt, F; /* time, distance, timestep, force */
310 int transition=1; /* boolean marking a transition for tension and timestep optimization */
312 get_options(argc, argv, &P, &dF, &t_max, &dt_max, &v, &x_step,
313 &stop_state, &env, NUM_TENSION_MODELS, tension_models,
314 NUM_K_MODELS, k_models, &state_list, &transition_list);
317 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Distance (m)\tForce (N)\tStiffness (N/m)\n");
318 if (flags & FLAG_OUTPUT_TRANSITION_FORCES) printf("#Force (N)\tInitial state\tFinal state\n");
319 while ((t_max > 0 && t <= t_max) || (stop_state != NULL && stop_state->num_domains > 0)) {
320 F = find_tension(state_list, &env, x, transition, 0);
322 dt = determine_dt(state_list, transition_list, &env, F, x, dF, P, dt_max, v, transition);
324 dt = determine_dt(state_list, transition_list, &env, F, x, dF, P, dt_max, 0, transition);
325 transition=random_transitions(transition_list, F, dt, &env);
326 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness);
331 if (dt == dt_max) { /* step completed */
334 } else { /* still working on this step */
339 destroy_state_list(state_list);
340 destroy_transition_list(transition_list);
344 @ The meat of the simulation is bundled into the three functions
345 [[find_tension]], [[determine_dt]], and [[random_transitions]].
346 [[find_tension]] is discussed in Section \ref{sec.find_tension},
347 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
348 [[random_transitions]] in Section \ref{sec.transition_rate}.
350 The stretched distance is extended in one of two modes depending on
351 [[xstep]]. If [[xstep == 0]], then the pulling is continuous, at
352 least within the limits of the inherent discretization of a time
353 stepped approach. At any rate, the timestep is limited by the maximum
354 allowed timestep [[dt_max]] and the maximum allowed unfolding
355 probability in a given step [[P]]. In the second mode [[xstep > 0]],
356 and the end to end distance increases in discrete steps of that
357 length. The time between steps is chosen to maintain an average
358 unfolding velocity [[v]]. This approach is less work to simulate than
359 smooth pulling and also closer to the reality of many experiments, but
360 it is practically impossible to treat analytically. With both pulling
361 styles implemented, the effects of the discretization can be easily
362 calculated, bridging the gap between theory and experiment.
364 Environmental parameters important in determining reaction rates and
365 tensions (e.g.\ temperature, pH) are stored in a single structure to
366 facilitate extension to more complicated models in the future.
367 <<environment definition>>=
368 typedef struct environment_struct {
378 #define DOUBLE_PRECISION 1e-12
380 <<environment definition>>
381 <<one dimensional function definition>>
382 <<create/destroy definitions>>
384 #endif /* GLOBAL_H */
386 Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
387 included multiple times.
389 \section{Simulation functions}
391 <<simulation functions>>=
395 <<does domain transition>>
396 <<randomly transition domains>>
400 \label{sec.find_tension}
402 Because the stretched system may be made up of several parts (folded
403 domains, unfolded domains, spring-like cantilever, \ldots), we will
404 assign the domains to states with tension models and groups. The
405 models are used to determine the tension of all the domain in a given
406 state following some model (Hookean spring, WLC, \ldots). The domains
407 are assumed to be commutative, so ordering is ignored. The
408 interactions between the groups are assumed to be linear, but the
409 interactions between domains of the same group need not be. Each
410 model has a tension handler function, which gives the tension $F$
411 needed to stretch a given group of domains a distance $x$.
413 Groups represent collections of states which the model should treat as
414 a single entity. To understand the purpose of groups, consider a
415 system with two unfolded domain states (e.g.\ for two different
416 proteins) that were both modeled as WLCs. If both unfolded states had
417 the same persistence length, it would be wise to assign them to the
418 same group. This would reduce both the computational cost of
419 balancing the tension and the error associated with the inter-group
420 interaction linearity. Note that group numbers only matter
421 \emph{within} model classes, since grouping domains with seperate
422 models doesn't make sense. Within designated groups, the tension
423 parameters for each domain are still checked for compatibility, so if
424 you accidentally assign incompatible domains to the same group you
425 should get an appropriate error message.
427 <<find tension definitions>>=
428 #define NUM_TENSION_MODELS 6
432 <<tension handler array global declaration>>=
433 extern one_dim_fn_t *tension_handlers[];
436 <<tension handler array global>>=
437 one_dim_fn_t *tension_handlers[] = {
439 &const_tension_handler,
447 <<tension model global declarations>>=
448 <<tension handler array global declaration>>
451 <<tension handler types>>=
452 typedef struct tension_handler_data_struct {
453 list_t *group_tension_params; /* one item for each domain in the group */
456 } tension_handler_data_t;
460 After sorting the chain into separate groups $G_i$, with tension
461 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
462 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
463 \forall i,j$. Note that there may be several states within each
464 group. [[find_tension]] is basically a wrapper to feed proper input
465 into the [[tension_balance]] function. [[transition]] should be set
466 to 0 if there were no transitions since the previous call to
467 [[find_tension]] to support some optimizations that assume a small
468 increase in tension between steps (only true if no transition has
469 occured). While we're messing with the tension handlers and if
470 [[const_env==0]], we also grab a measure of the current linker
471 stiffness for the environmental variable, which is needed by some of
472 the unfolding rate models (e.g.\ the linker-stiffness corrected Bell
473 model, Appendix \ref{sec.kbell}).
475 double find_tension(list_t *state_list, environment_t *env, double x,
476 int transition, int const_env)
477 { // TODO: !cont_env -> step_call, only alter env or statics if step_call==1
478 static int num_active_groups=0;
479 static one_dim_fn_t **per_group_handlers = NULL;
480 static one_dim_fn_t **per_group_inverse_handlers = NULL;
481 static void **per_group_data = NULL; /* array of pointers to tension_handler_data_t */
482 static double last_x;
483 tension_handler_data_t *tdata;
484 double F, *pStiffness;
488 fprintf(stderr, "%s:\tfinding tension at distance %g%s\n", __FUNCTION__, x, transition ? " after a transition" : "");
491 if (transition != 0 || num_active_groups == 0) { /* setup statics */
492 /* get new statics, freeing the old ones if they aren't NULL */
493 get_active_groups(state_list, &num_active_groups, &per_group_handlers, &per_group_inverse_handlers, &per_group_data);
495 /* fill in the group handlers and params */
497 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]);
499 for (i=0; i<num_active_groups; i++) {
500 tdata = (tension_handler_data_t *) per_group_data[i];
502 fprintf(stderr, "%s:\tactive group %d of %d (data %p, param list %p)", __FUNCTION__, i, num_active_groups, tdata, tdata->group_tension_params);
503 list_print(stderr, tdata->group_tension_params, " parameters");
506 /* tdata->persist continues unchanged */
509 } /* else continue with the current statics */
514 pStiffness = &env->stiffness;
516 F = tension_balance(num_active_groups, per_group_handlers, per_group_inverse_handlers, (void **)per_group_data, last_x, x, pStiffness);
522 @ For the moment, we will restrict our group boundaries to a single
523 dimension, so $\sum_i x_i = x$, our total extension (it is this
524 restriction that causes the groups to interact linearly). We'll also
525 restrict our extensions to all be positive. With these restrictions,
526 the problem of balancing the tensions reduces to solving a system of
527 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
528 the number of active groups. In general this can be fairly
529 complicated, but without much loss of practicality we can restrict
530 ourselves to strictly monotonically increasing, non-negative tension
531 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
532 simpler. For example, it guarantees the existence of a unique, real
533 solution for finite forces. See Appendix \ref{sec.tension_balance}
534 for the tension balancing implementation.
536 Each group must define a parameter structure if the tension model is
538 a function to create the parameter structure,
539 a function to destroy the parameter structure, and
540 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as its data argument.
541 The tension-specific parameter structures are contained in the domain
542 group field. For optimization, a group may also define an inverse
543 tension function as an optimization. See the Section
544 \ref{sec.model_selection} for a discussion on the command line
545 framework and inverse function discussion. See the worm-like chain
546 implementation for an example (Section \ref{sec.wlc}).
548 The major limitation of our approach is that all of our tension models
549 are Markovian (memory-less), which is not really a problem for the
550 chains (see \citet{evans99} for WLC thermalization rate calculations).
552 \subsection{Transition rate}
553 \label{sec.transition_rate}
555 Each state transition is modeled as unimolecular, first order reaction
557 \text{State 1} \xrightarrow{k} \text{State 2}
559 With the rate constant $k$ defined by
561 \frac{d[\text{State 2}]}{dt} = -\frac{d[\text{State 1}]}{dt} = k [\text{State 1}].
563 So the probability of a given domain transitioning in a short time
569 For a group of $N$ identical domains, and therefore identical
570 unfolding probabilities $P_1$, the probability of $n$ domains
571 unfolding in a given timestep is given by the binomial distribution
573 P(n) = \frac{N!}{n!(N-n)!}P_1^n(1-P_1)^{(N-n)}.
575 The probability that \emph{none} of these domains unfold is then
579 so the probability that \emph{at least one} domain unfolds is
583 Note that for small $P$,
585 p(n>0) = 1-(1-NP_1+\mathcal{O}(P_1^2)) = NP_1 - \mathcal{O}(P^2)
588 This conversion from $P_1$ to $P_N$ is handled by the [[pN]] macro
590 /* find multi-domain transition rate from the single domain rate */
591 #define PN(P,N) (1.0-pow(1.0-(P), (N)))
595 We can also define a macro to find the approprate $dt$ to achieve a
596 target $P_N$ with a given $k$ and $N$.
598 P_N &= 1-(1-P_1)^N \\
599 1-P_1 &= (1-P_N)^{1/N} \\
600 P_1 &= 1-(1-P_N)^{1/N}
603 #define P1(PN,N) (1.0-pow(1.0-(PN), 1.0/(N)))
606 We take some time to discuss the meaning of $p(n>1)$
607 (i.e. multi-unfolding timesteps) in Section \ref{sec.timescales}.
609 <<does domain transition>>=
610 int domain_transitions(double F, double dt, environment_t *env, transition_t *transition)
611 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to transition structure */
613 k = accel_k(transition->k, F, env, transition->k_params);
614 //(*transition->k)(F, env, domain->k_params);
615 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
616 return happens(PN(k*dt, N_INIT(transition)));
618 @ [[happens]] is a random decision making function defined in Appendix \ref{sec.utils}.
620 <<randomly transition domains>>=
621 int random_transitions(list_t *transition_list, double tension,
622 double dt, environment_t *env)
623 { /* returns 1 if there was a transition and 0 otherwise */
624 while (transition_list != NULL) {
625 if (T(transition_list)->initial_state->num_domains > 0
626 && domain_transitions(tension, dt, env, T(transition_list))) {
627 if (flags & FLAG_OUTPUT_TRANSITION_FORCES)
628 fprintf(stdout, "%g\t%s\t%s\n", tension, T(transition_list)->initial_state->name, T(transition_list)->final_state->name);
629 T(transition_list)->initial_state->num_domains--;
630 T(transition_list)->final_state->num_domains++;
631 return 1; /* our one domain has transitioned, stop looking for others */
633 transition_list = transition_list->next;
639 \subsection{Timescales}
640 \label{sec.timescales}
642 The simulation assumes that chain equilibration is instantaneous
643 relative to the simulation timestep $dt$, so that tension is uniform
644 along the chain. The quality of this assumption depends on your
645 particular chain. For example, a damped spring thermalizes on a
646 timescale of order $\tau = 1/\gamma$ where the damping ratio $\gamma
647 \equiv \omega_0(-\zeta \pm \sqrt{\zeta^2-1}$, the natural angular
648 frequency $omega_0 \equiv \sqrt{k/m}$, $\zeta \equiv b/2\sqrt{km}$,
649 and $b$ sets the drag $F_b = -bv$. For our cantilevers $\tau$ is
650 on the order of milliseconds, which is longer than a timestep.
651 However, the approximation is still reasonable, because there is
652 rarely unfolding during the cantilever return. You could, of course,
653 take cantilever, WLC, etc.\ response times into effect, but that
654 would significantly complicate a simulation that seems to work
655 well enough without bothering :p. For WLC and FJC relaxation times,
656 see the Appendix of \citet{evans99} and \citet{kroy07}.
658 Besides assuming our timestep is much greater than equilibration
659 times, we also force it to be short enough so that only one domain may
660 unfold in a given timestep. For an unfolding timescale of $dt_u =
661 1/k$, we adapt our timesteps to keep $dt \ll dt_u$, so the probability
662 of two domains unfolding in the same timestep is small (see
663 [[determine_dt]] in Section \ref{sec.adaptive_dt} and [[P]] in
664 [[main()]] in Section \ref{sec.main} set by [[-P]] in
665 [[get_options()]] in Section \ref{sec.get_options}). This approach
666 breaks down as the adaptive timestep scheme approaches the transition
667 timestep $dt_t$, but $dt_t \sim 1\U{$\mu$s}$ for Ig27-like proteins
668 \citep{klimov00}, so this shouldn't be a problem. To reassure
669 yourself, you can ask the simulator to print the smallest timestep
670 that was used with TODO.
672 Even if the likelyhood of two domains unfolding in the same timestep
673 is small, if you run long enough simulations it will eventually occur.
674 In this case, we assume that $dt_t \ll dt$, so even if two domains
675 would unfold if we stuck to our unfolding rate $k$ for an entire
676 timestep $dt$, in ``reality'' only one of those domains unfolds.
677 Because we do not know when in the timestep the unfolding took place,
678 we move on to the next timestep without checking for further
679 unfoldings. This ``unchecked timestep portion'' should not contribute
680 significant errors to the calculation, because if $dt$ was small
681 enough that unfolding was unlikely at the pre-unfolding force, it
682 should be even more unlikely at the post-unfolding force, especially
683 over only a fraction of the timestep. The only dangerous cases for
684 the current implementation are those where unfolding intermediates are
685 much less stable than their precursor states, in which case an
686 unfolding event during the remaining timestep may actually be likely.
687 A simple workaround for such cases is to pick the value for [[P]] to
688 be small enough that the $dt \ll dt_u$ assumption survives even under
689 a transition populating the unstable intermediate.
691 \subsection{Adaptive timesteps}
692 \label{sec.adaptive_dt}
694 We'd like to pick $dt$ so the probability of changing state
695 (e.g.\ unfolding) in the next timestep is small. If we don't adapt our
696 timestep, we also risk breaking our approximation that $F(x' \in [x,
697 x+v \cdot dt]) = F(x)$ or that only one domain changes state in a
698 given timestep. The simulation would have been accurate for
699 sufficiently small fixed timesteps, but adaptive timesteps will allow
700 us to move through low tension regions in fewer steps, leading to a
701 more efficient simulation.
703 Consider the two state folded $\rightarrow$ unfolded transition.
704 Because $F(x)$ is monotonically increasing between unfoldings,
705 excessively large timesteps will lead to erroneously small unfolding
706 rates (an thus, higher average unfolding force).
708 The actual adaptive timestep implementation is not particularly
709 interesting, since we are only required to reduce $dt$ to somewhere
710 below a set threshold, so I've removed it to Appendix
711 \ref{sec.adaptive_dt_implementation}.
717 The models provide the physics of the simulation, but the simulation
718 allows interchangeable models, and we are currently focusing on the
719 simulation itself, so we remove the actual model implementation to the
722 The tension models are defined in Appendix \ref{sec.tension_models}
723 and unfolding models are defined in Appendix \ref{sec.k_models}.
726 #define k_B 1.3806503e-23 /* J/K */
730 \section{Command line interface}
732 <<get option functions>>=
734 <<index tension model>>
739 \subsection{Model selection}
740 \label{sec.model_selection}
742 The main difficulty with the command line interface in \prog\ is
743 developing an intuitive method for accessing the possibly large number
744 of available models. We'll treat this generally by defining an array
745 of available models, containing their important parameters.
747 <<get options definitions>>=
748 <<model getopt definitions>>
751 <<create/destroy definitions>>=
752 typedef void *create_data_func_t(char **param_strings);
753 typedef void destroy_data_func_t(void *param_struct);
756 <<model getopt definitions>>=
757 <<tension model getopt definitions>>
758 <<k model getopt definitions>>
761 In order to choose models, we need a method of assembling a reaction
762 graph such as that in Figure \ref{fig.domain_states} and population
763 the graph with a set of domains. First we declare the domain states
766 [[-s <state-name>,<model>[,<group>] -S <model-params> [-N <initial-population]]] \\
770 [[-s unfolded,wlc,1 -S 1e-8,4e-10 -N 5]]
772 This creates the state named [[unfolded]], using the WLC model and one
773 as a group number (defaults to 0, see Section \ref{sec.find_tension}).
774 The tension parameters are then set to [[1e-8,4e-10]], which in the
775 case of the WLC are the contour and persistence lengths respectively.
776 Finally we populate the state by adding five domains to the state.
777 The name of the state is importent for identifying states later.
779 Once the domains have been created and populated, you can added
780 transition rates following
782 [[-k <initial-state>,<final-state>,<model> -K <model-params>]] \\
786 [[-k folded,unfolded,bell -K 3.3e-4,0.25e-9]]
788 This creates a reaction going from the [[folded]] state to the
789 [[unfolded]] state, with the rate constant given by the Bell model
790 (Appendix \ref{sec.bell}). The unfolding parameters are then set to
791 [[3.3e-4,0.25e-9]], which in the case of the Bell model are the
792 unforced rate and transition state distance respectively.
794 \subsubsection{Tension}
796 To access the various tension models from the command line, we define
797 a structure that contains all the neccessary information about the
799 <<tension model getopt definitions>>=
800 typedef struct tension_model_getopt_struct {
801 one_dim_fn_t *handler; /* fn implementing the model on a group */
802 one_dim_fn_t *inverse_handler; /* fn implementing inverse on group */
803 char *name; /* name string identifying the model */
804 char *description; /* a brief description of the model */
805 int num_params; /* number of extra params for the model */
806 char **param_descriptions; /* descriptions of extra params */
807 char *params; /* default values for the extra params */
808 create_data_func_t *creator; /* param-string -> model-data-struct fn */
809 destroy_data_func_t *destructor; /* fn to free the model data structure */
810 } tension_model_getopt_t;
811 @ Set [[inverse_handler]] to [[NULL]] if you wish to invert using
812 bisection. Obviously, this will be slower than computing an
813 analytical inverse and more prone to rounding errors, but it may be
814 easier if a closed-form inverse does not exist.
816 <<tension model getopt array definition>>=
817 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
818 <<null tension model getopt>>,
819 <<constant tension model getopt>>,
820 <<hooke tension model getopt>>,
821 <<worm-like chain tension model getopt>>,
822 <<freely-jointed chain tension model getopt>>,
823 <<piston tension model getopt>>
827 \subsubsection{Unfolding rate}
829 <<k model getopt definitions>>=
830 #define NUM_K_MODELS 6
832 typedef struct k_model_getopt_struct {
837 char **param_descriptions;
839 create_data_func_t *creator;
840 destroy_data_func_t *destructor;
844 <<k model getopt array definition>>=
845 k_model_getopt_t k_models[NUM_K_MODELS] = {
846 <<null k model getopt>>,
847 <<const k model getopt>>,
848 <<bell k model getopt>>,
849 <<kbell k model getopt>>,
850 <<kramers k model getopt>>,
851 <<kramers integ k model getopt>>
858 void help(char *prog_name, double P, double dF,
859 double t_max, double dt_max, double v, double x_step,
860 state_t *stop_state, environment_t *env,
861 int n_tension_models, tension_model_getopt_t *tension_models,
862 int n_k_models, k_model_getopt_t *k_models,
863 int k_model, list_t *state_list)
868 if (state_list != NULL)
869 state = S(tail(state_list));
871 printf("usage: %s [options]\n", prog_name);
872 printf("Version %s\n\n", VERSION);
873 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
875 printf("Simulation options:\n");
877 printf("-q state\tStop simulation when this state empties (- for no limit, currently %s)\n", stop_state ? stop_state->name : "-" );
878 printf("-t tmax\tMaximum simulated time (-1 for no limit, currently %g s)\n", t_max);
879 printf("-d dt\tMaximum allowed timestep dt (currently %g s)\n", dt_max);
880 printf("-P P\tMaximum transition probability for dt (currently %g)\n", P);
881 printf("-F dF\tMaximum change in force over dt. Only relevant if dx > 0. (currently %g N)\n", dF);
882 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
883 printf("-x dx\tPulling step size (currently %g m)\n", x_step);
884 printf("\tWhen dx = 0, the pulling is continuous.\n");
885 printf("\tWhen dx > 0, the pulling occurs in discrete steps.\n");
887 printf("Environmental options:\n");
888 printf("-T T\tTemperature T (currently %g K)\n", env->T);
889 printf("-C T\tYou can also set the temperature T in Celsius\n");
891 printf("State creation options:\n");
892 printf("-s name,model[[,group],params]\tCreate a new state.\n");
893 printf("Once you've set up the state, you may populate it with domains\n");
894 printf("-N n\tAdd n domains in the current state (%s)\n", state ? state->name : "-");
896 printf("Transition creation options:\n");
897 printf("-k initial_state,final_state,model[,params]\tCreate a new transition\n");
899 printf("To double check your reaction graph, you can produce graphviz dot output\n");
900 printf("describing the current states and transitions.\n");
901 printf("-D\tPrint dot description of states and transitions and exit.\n");
903 printf("Simulation output mode options:\n");
904 printf("There are two output modes. In standard mode, only the transition\n");
905 printf("events are printed. For example:\n");
906 printf(" #Force (N)\tInitial state\tFinal state\n");
907 printf(" 123.456e-12\tfolded\tunfolded\n");
909 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
910 printf(" #Distance (m)\tForce (N)\tStiffness (N/m)\n");
911 printf(" 0.001\t0.0023\t0.002\n");
913 printf("-V\tChange output to verbose mode.\n");
914 printf("-h\tPrint this help and exit.\n");
917 printf("Tension models:\n");
918 for (i=0; i<n_tension_models; i++) {
919 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
920 for (j=0; j < tension_models[i].num_params ; j++)
921 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
922 printf("\t\tdefaults: %s\n", tension_models[i].params);
924 printf("Unfolding rate models:\n");
925 for (i=0; i<n_k_models; i++) {
926 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
927 for (j=0; j < k_models[i].num_params ; j++)
928 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
929 printf("\t\tdefaults: %s\n", k_models[i].params);
932 printf("Examples:\n");
933 printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded. Stop once there are no more folded domains.\n");
934 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);
935 printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded. Stop after 0.25 seconds\n");
936 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);
940 \subsection{Get options}
941 \label{sec.get_options}
945 void get_options(int argc, char **argv, double *pP, double *pDF,
946 double *pT_max, double *pDt_max, double *pV,
947 double *pX_step, state_t **pStop_state, environment_t *env,
948 int n_tension_models, tension_model_getopt_t *tension_models,
949 int n_k_models, k_model_getopt_t *k_models,
950 list_t **pState_list, list_t **pTransition_list)
952 char *prog_name = NULL;
953 char c, options[] = "q:t:d:P:F:v:x:T:C:s:N:k:DVh";
954 int tension_model=0, k_model=0, initial_state=-1, final_state=-1;
955 char *state_name=NULL;
956 int i, n, tension_group=0;
957 list_t *temp_list=NULL;
960 void *model_params; /* for INIT_MODEL() */
961 int num_param_args; /* for INIT_MODEL() */
962 char **param_args; /* for INIT_MODEL() */
964 extern int optind, optopt, opterr;
969 for (i=0; i<n_tension_models; i++) {
970 fprintf(stderr, "tension model %i %s:\t model's handler %p\t inverse handler %p\n",
971 i, tension_models[i].name, tension_models[i].handler, tension_models[i].inverse_handler);
972 assert(tension_models[i].handler == tension_handlers[i]);
977 flags = FLAG_OUTPUT_TRANSITION_FORCES;
980 *pT_max = -1; /* s, -1 removes this limit */
981 *pDt_max = 0.001; /* s */
982 *pP = 1e-3; /* % pop per s */
983 *pDF = 1e-12; /* N */
984 *pV = 1e-6; /* m/s */
985 *pX_step = 0.0; /* m */
986 env->T = 300.0; /* K */
988 *pTransition_list = NULL;
990 while ((c=getopt(argc, argv, options)) != -1) {
993 if (STRMATCH(optarg, "-")) {
996 temp_list = head(*pState_list);
997 while (temp_list != NULL) {
998 if (STRMATCH(S(temp_list)->name, optarg)) {
999 *pStop_state = S(temp_list);
1002 temp_list = temp_list->next;
1006 case 't': *pT_max = safe_strtod(optarg, "-t"); break;
1007 case 'd': *pDt_max = safe_strtod(optarg, "-d"); break;
1008 case 'P': *pP = safe_strtod(optarg, "-P"); break;
1009 case 'F': *pDF = safe_strtod(optarg, "-F"); break;
1010 case 'v': *pV = safe_strtod(optarg, "-v"); break;
1011 case 'x': *pX_step = safe_strtod(optarg, "-x"); break;
1012 case 'T': env->T = safe_strtod(optarg, "-T"); break;
1013 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
1015 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
1016 assert(num_strings >= 2 && num_strings <= 4);
1018 state_name = strings[0];
1019 if (state_index(*pState_list, state_name) != -1) {
1020 fprintf(stderr, "Invalid option '%s', repeated state name '%s'\n", optarg, state_name);
1023 tension_model = index_tension_model(n_tension_models, tension_models, strings[1]);
1024 if (num_strings == 4)
1025 tension_group = safe_strtoi(strings[2], optarg);
1028 if (num_strings >= 3)
1029 tension_models[tension_model].params = strings[num_strings-1];
1031 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);
1033 INIT_MODEL(state_name, &tension_models[tension_model], tension_models[tension_model].params, model_params);
1034 push(pState_list, create_state(state_name,
1035 tension_models[tension_model].handler,
1036 tension_models[tension_model].inverse_handler,
1039 tension_models[tension_model].destructor,
1041 *pState_list = tail(*pState_list);
1043 free_parsed_list(num_strings, strings);
1046 n = safe_strtoi(optarg, "-N");
1048 fprintf(stderr, "%s:\tadding %d domains to %s\n", __FUNCTION__, n, S(*pState_list)->name);
1050 S(*pState_list)->num_domains += n;
1053 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
1054 assert(num_strings >= 3 && num_strings <= 4);
1055 initial_state = state_index(*pState_list, strings[0]);
1056 final_state = state_index(*pState_list, strings[1]);
1057 k_model = index_k_model(n_k_models, k_models, strings[2]);
1059 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);
1061 assert(initial_state != final_state);
1062 if (num_strings == 4)
1063 k_models[k_model].params = strings[3];
1064 INIT_MODEL(state_name, &k_models[k_model], k_models[k_model].params, model_params);
1065 push(pTransition_list, create_transition(list_element(*pState_list, initial_state),
1066 list_element(*pState_list, final_state),
1067 k_models[k_model].k,
1069 k_models[k_model].destructor), 1);
1070 free_parsed_list(num_strings, strings);
1073 print_state_graph(stdout, *pState_list, *pTransition_list);
1077 flags = FLAG_OUTPUT_FULL_CURVE;
1080 fprintf(stderr, "unrecognized option '%c'\n", optopt);
1081 /* fall through to default case */
1083 help(prog_name, *pP, *pDF, *pT_max, *pDt_max, *pV, *pX_step,
1084 *pStop_state, env, n_tension_models, tension_models,
1085 n_k_models, k_models, k_model, *pState_list);
1089 /* check the arguments */
1090 assert(*pState_list != NULL); /* no domains! */
1091 assert(*pP > 0.0); assert(*pP < 1.0);
1092 assert(*pDt_max > 0.0);
1094 assert(env->T > 0.0);
1096 crosslink_groups(*pState_list, *pTransition_list);
1102 <<index tension model>>=
1103 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
1106 for (i=0; i<n_models; i++)
1107 if (STRMATCH(models[i].name, name))
1109 if (i == n_models) {
1110 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
1118 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
1121 for (i=0; i<n_models; i++)
1122 if (STRMATCH(models[i].name, name))
1124 if (i == n_models) {
1125 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
1132 <<initialize model definition>>=
1133 /* requires int num_param_args and char **param_args in the current scope
1135 * INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
1136 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
1138 #define INIT_MODEL(role, model, param_string, param_pointer) \
1140 parse_list_string((param_string), SEP, '{', '}', \
1141 &num_param_args, ¶m_args); \
1142 if (num_param_args != (model)->num_params) { \
1144 "%s model %s expected %d params,\n", \
1145 (role), (model)->name, (model)->num_params); \
1147 "not the %d params in '%s'\n", \
1148 num_param_args, (param_string)); \
1149 assert(num_param_args == (model)->num_params); \
1151 if ((model)->creator) \
1152 param_pointer = (*(model)->creator)(param_args); \
1154 param_pointer = NULL; \
1155 free_parsed_list(num_param_args, param_args); \
1159 First we define some safe string parsing functions. Currently these
1160 abort on any irregularity, but in the future they could print error
1161 messages, etc. [[static]] because the functions are currently
1162 defined in each [[*.c]] file that uses them, but perhaps they should
1163 be moved to [[utils.h/c]] or some such instead\ldots
1166 static int safe_strtoi(const char *s, const char *description) {
1170 ret = strtol(s, &endp, 10);
1171 if (endp[0] != '\0') { //strlen(endp) > 0) {
1172 fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1173 endp, s, description);
1174 assert(1==0); //strlen(endp) == 0);
1176 fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1179 } else if (errno == ERANGE) {
1180 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1181 assert(errno != ERANGE);
1186 static double safe_strtod(const char *s, const char *description) {
1190 ret = strtod(s, &endp);
1191 if (endp[0] != '\0') { //strlen(endp) > 0) {
1192 fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1193 endp, s, description);
1194 assert(1==0); //strlen(endp) == 0);
1196 fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1199 } else if (errno == ERANGE) {
1200 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1201 assert(errno != ERANGE);
1210 \addcontentsline{toc}{section}{Appendicies}
1212 \section{sawsim.c details}
1216 The general layout of our simulation code is:
1227 We include [[math.h]], so don't forget to link to the libm with `-lm'.
1229 #include <assert.h> /* assert() */
1230 #include <stdlib.h> /* malloc(), free(), rand(), strto*() */
1231 #include <stdio.h> /* fprintf(), stdout */
1232 #include <string.h> /* strlen, strtok() */
1233 #include <math.h> /* exp(), M_PI, sqrt() */
1234 #include <sys/types.h> /* pid_t (returned by getpid()) */
1235 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
1236 #include <errno.h> /* errno, ERANGE (for safe_strto*()) */
1239 #include "tension_balance.h"
1240 #include "k_model.h"
1241 #include "tension_model.h"
1243 #include "accel_k.h"
1247 <<version definition>>
1248 <<flag definitions>>
1249 <<max/min definitions>>
1250 <<string comparison definition and globals>>
1251 <<initialize model definition>>
1252 <<get options definitions>>
1253 <<state definitions>>
1254 <<transition definitions>>
1263 <<create/destroy state>>
1264 <<destroy state list>>
1266 <<create/destroy transition>>
1267 <<destroy transition list>>
1268 <<print state graph>>
1270 <<simulation functions>>
1271 <<boilerplate functions>>
1274 <<boilerplate functions>>=
1276 <<get option functions>>
1282 srand(getpid()*time(NULL)); /* seed rand() */
1286 <<flag definitions>>=
1287 /* in octal b/c of prefixed '0' */
1288 #define FLAG_OUTPUT_FULL_CURVE 01
1289 #define FLAG_OUTPUT_TRANSITION_FORCES 02
1293 static unsigned long int flags = 0;
1296 \subsection{Utilities}
1299 <<max/min definitions>>=
1300 #define MAX(a,b) ((a)>(b) ? (a) : (b))
1301 #define MIN(a,b) ((a)<(b) ? (a) : (b))
1304 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
1305 <<string comparison definition and globals>>=
1306 // Check if two strings match, return 1 if they do
1307 static char *temp_string_A;
1308 static char *temp_string_B;
1309 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
1310 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
1311 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
1312 /* +1 to also compare the '\0' */
1315 We also define a macro for our [[check]] unit testing
1316 <<check relative error>>=
1317 #define CHECK_ERR(max_err, expected, received) \
1319 fail_unless( (received-expected)/expected < max_err, \
1320 "relative error %g >= %g in %s (Expected %g, received %g)", \
1321 (received-expected)/expected, max_err, #received, \
1322 expected, received); \
1323 fail_unless(-(received-expected)/expected < max_err, \
1324 "relative error %g >= %g in %s (Expected %g, received %g)", \
1325 -(received-expected)/expected, max_err, #received, \
1326 expected, received); \
1331 int happens(double probability)
1333 assert(probability >= 0.0); assert(probability <= 1.0);
1334 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*/
1339 \subsection{Adaptive timesteps}
1340 \label{sec.adaptive_dt_implementation}
1342 $F(x)$ increases with $x$, possibly exploding (e.g.\ in the worm-like
1343 chain model), so basing the timestep on the unfolding probability at
1344 the current tension is dangerous, and we need to search for a $dt$ for
1345 which $P_N(F(x+v*dt))<P_\text{max}$ (See Section
1346 \ref{sec.transition_rate} for a discussion of $P_N$). There are two
1347 cases to consider. In the most common, no domains have unfolded since
1348 the last step, and we expect the next step to be slightly shorter than
1349 the previous one. In the less common, domains did unfold in the last
1350 step, and we expect the next step to be considerably longer than the
1353 double search_dt(transition_t *transition,
1355 environment_t *env, double x,
1356 double max_prob, double max_dt, double v)
1357 { /* Returns a valid timestep dt in seconds for the given transition.
1358 * Takes a list of all states, a pointer env to the environmental data,
1359 * a starting separation x in m, a max_prob between 0 and 1,
1360 * max_dt in s, stretching velocity v in m/s.
1362 double F, k, P, dtCur, dtU, dtUCur, dtL, dt;
1364 /* get upper bound using the current position */
1365 F = find_tension(state_list, env, x, 0, 1); /* BUG. repeated calculation */
1366 //printf("Start with x = %g (F = %g)\n", x, F);
1367 k = accel_k(transition->k, F, env, transition->k_params);
1368 //printf("x %g\tF %g\tk %g\n", x, F, k);
1369 dtU = P1(max_prob, N_INIT(transition)) / k; /* upper bound on dt */
1371 //printf("overshot max_dt\n");
1374 if (v == 0) /* discrete stepping, so force is temporarily constant */
1377 /* set a lower bound on dt too */
1380 /* The dt determined above may produce illegitimate forces or ks.
1381 * Reduce the upper bound until we have valid ks. */
1383 F = find_tension(state_list, env, x+v*dt, 0, 1);
1384 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1387 F = find_tension(state_list, env, x+v*dt, 0, 1);
1389 //printf("Try for dt = %g (F = %g)\n", dt, F);
1390 k = accel_k(transition->k, F, env, transition->k_params);
1391 /* returned k may be -1.0 */
1392 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1393 while (k == -1.0) { /* reduce step until we hit a valid k */
1395 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1396 F = find_tension(state_list, env, x+v*dt, 0, 1);
1397 //printf("Try for dt = %g (F = %g)\n", dt, F);
1398 k = accel_k(transition->k, F, env, transition->k_params);
1399 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1401 assert(dtU > 1e-14); /* timestep to valid k too small */
1402 /* safe timestep back from x+dtU */
1403 dtUCur = P1(max_prob, N_INIT(transition)) / k;
1405 return dt; /* dtU is safe. */
1407 /* dtU wasn't safe, lets see what would be. */
1408 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1409 dt = (dtU + dtL) / 2.0;
1410 F = find_tension(state_list, env, x+v*dt, 0, 1);
1411 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1412 k = accel_k(transition->k, F, env, transition->k_params);
1413 dtCur = P1(max_prob, N_INIT(transition)) / k;
1414 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1415 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1417 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1418 dtU = dt; /* ... stepping out only dtCur would be safe */
1421 break; /* dtCur = dt */
1423 return MAX(dtUCur, dtL);
1428 Determine $dt$ for an array of potentially different folded domains.
1429 We apply [[search_dt()]] to each populated domain to find the maximum
1430 $dt$ the domain can handle. The final $dt$ is less than the
1431 individual domain $dt$s (to satisfy [[search_dt()]], see above). If
1432 $v>0$ (continuous pulling), $dt$ is also reduced to satisfy
1433 $|F(x+v*dt)-F(x)|<dF_\text{max}$, which limits errors due to our
1434 transition rate assumption that the force does not change appreciably
1435 over a single timestep.
1438 double determine_dt(list_t *state_list, list_t *transition_list,
1439 environment_t *env, double F, double x, double max_dF,
1440 double max_prob, double dt_max, double v, int transition)
1441 { /* Returns the timestep dt in seconds.
1442 * Takes the state and transition lists, a pointer to the environment,
1443 * the force F(x) in N, an extension x in m, max_dF in N,
1444 * max_prob between 0 and 1, dt_max in s, v in m/s, and transition in {0,1}*/
1445 double dt=dt_max, new_dt, F_new;
1446 assert(max_dF > 0.0);
1447 assert(max_prob > 0.0); assert(max_prob < 1.0);
1448 assert(dt_max > 0.0);
1449 assert(state_list != NULL);
1450 assert(transition_list != NULL);
1451 transition_list = head(transition_list);
1454 /* Ensure |dF| < max_dF */
1455 F_new = find_tension(state_list, env, x+v*dt, transition, 1);
1456 while (fabs(F_new - F) > max_dF) {
1458 F_new = find_tension(state_list, env, x+v*dt, transition, 1);
1462 /* Ensure all individual domains pass search_dt() */
1463 while (transition_list != NULL) {
1464 if (T(transition_list)->initial_state->num_domains > 0) {
1465 new_dt = search_dt(T(transition_list),state_list,env,x,max_prob,dt,v);
1466 dt = MIN(dt, new_dt);
1468 transition_list = transition_list->next;
1475 \subsection{State data}
1476 \label{sec.state_data}
1478 The domains exist in one of an array of possible states. Each state
1479 has a tension model which determines it's force/extension curve
1480 (possibly [[null]] for rigid states, see Appendix
1481 \ref{sec.null_tension}).
1483 Groups are, for example, for two domain states with WLC tensions; one
1484 with $p=0.4\U{nm}$, $L=10\U{nm}$ and the other with $p=0.4\U{nm}$,
1485 $L=28\U{nm}$. Since the persistence lengths are the same, you would
1486 like to add the contour lengths of all the domains in \emph{both}
1487 states, and plug that total contour length into a single WLC
1488 evaluation (see Section \ref{sec.find_tension}).
1489 <<state definitions>>=
1490 typedef struct state_struct {
1491 char *name; /* name string */
1492 one_dim_fn_t *tension_handler; /* tension handler */
1493 one_dim_fn_t *inverse_tension_handler; /* inverse tension handler */
1494 int tension_group; /* for grouping similar tension states */
1495 void *tension_params; /* pointer to tension parameters */
1496 destroy_data_func_t *destroy_tension;
1497 int num_domains; /* number of domains currently in state */
1498 /* cached lookups generated from the above data */
1499 int num_out_trans; /* length of out_trans array */
1500 struct transition_struct **out_trans; /* transitions leaving this state */
1501 int num_grouped_states; /* length of grouped states array */
1502 struct state_struct **grouped_states; /* states grouped with this one (doesn't include self) */
1505 /* get the state data for the current list node */
1506 #define S(list) ((state_t *)(list)->d)
1508 @ [[tension_params]] is a pointer to the parameters used by the
1509 handler function when determining the tension. [[destroy_tension]]
1510 points to a function for freeing [[tension_params]]. We it in the
1511 state structure so that [[destroy_state]] and will have an easier job
1514 [[create_]] and [[destroy_state]] are simple wrappers around
1515 [[malloc]] and [[free]].
1516 <<create/destroy state>>=
1517 state_t *create_state(char *name,
1518 one_dim_fn_t *tension_handler,
1519 one_dim_fn_t *inverse_tension_handler,
1521 void *tension_params,
1522 destroy_data_func_t *destroy_tension,
1525 state_t *ret = (state_t *)malloc(sizeof(state_t));
1526 assert(ret != NULL);
1527 /* make a copy of the name, so we aren't dependent on the original */
1528 ret->name = (char *)malloc(sizeof(char)*(strlen(name)+1));
1529 assert(ret->name != NULL);
1530 strcpy(ret->name, name); /* we know ret->name is long enough */
1532 ret->tension_handler = tension_handler;
1533 ret->inverse_tension_handler = inverse_tension_handler;
1534 ret->tension_group = tension_group;
1535 ret->tension_params = tension_params;
1536 ret->destroy_tension = destroy_tension;
1537 ret->num_domains = num_domains;
1539 ret->num_out_trans = 0;
1540 ret->out_trans = NULL;
1541 ret->num_grouped_states = 0;
1542 ret->grouped_states = NULL;
1545 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);
1550 void destroy_state(state_t *state)
1554 fprintf(stderr, "freeing state %s (%p)\n", state->name, state);
1558 if (state->destroy_tension)
1559 (*state->destroy_tension)(state->tension_params);
1560 if (state->out_trans)
1561 free(state->out_trans);
1562 if (state->grouped_states)
1563 free(state->grouped_states);
1570 We'll be storing the states in a list (see Appendix \ref{sec.lists}),
1571 so we also define a convenience function for cleaning up.
1572 <<destroy state list>>=
1573 void destroy_state_list(list_t *state_list)
1575 state_list = head(state_list);
1576 while (state_list != NULL)
1577 destroy_state((state_t *) pop(&state_list));
1582 We can also define a convenient way to get a state index from it's name.
1584 int state_index(list_t *state_list, char *name)
1587 state_list = head(state_list);
1588 while (state_list != NULL) {
1589 if (STRMATCH(S(state_list)->name, name)) return i;
1590 state_list = state_list->next;
1598 \subsection{Transition data}
1599 \label{sec.transition_data}
1601 Once you've created states (Sections \ref{sec.main},
1602 \ref{sec.model_selection}, and \ref{sec.state_data}), you need to
1603 create transitions between them.
1604 <<transition definitions>>=
1605 typedef struct transition_struct {
1606 state_t *initial_state;
1607 state_t *final_state;
1608 k_func_t *k; /* transition rate model */
1609 void *k_params; /* pointer to k parameters */
1610 destroy_data_func_t *destroy_k;
1613 /* get the transition data for the current list node */
1614 #define T(list) ((transition_t *)(list)->d)
1616 /* get the number of initial state population for a transition state */
1617 #define N_INIT(transition) ((transition)->initial_state->num_domains)
1621 @ [[k]] is a pointer to the function determining the transition rate
1622 for a given tension. [[k_params]] and [[destroy_k]] have similar
1623 roles with regards to [[k]] as [[tension_params]] and
1624 [[destroy_tension]] have with regards to [[tension_handler]] (see
1625 Section \ref{sec.state_data}).
1627 [[create_]] and [[destroy_transition]] are simple wrappers around
1628 [[malloc]] and [[free]].
1629 <<create/destroy transition>>=
1630 transition_t *create_transition(state_t *initial_state, state_t *final_state,
1633 destroy_data_func_t *destroy_k)
1635 transition_t *ret = (transition_t *)malloc(sizeof(transition_t));
1636 assert(ret != NULL);
1637 assert(initial_state != NULL);
1638 assert(final_state != NULL);
1639 ret->initial_state = initial_state;
1640 ret->final_state = final_state;
1642 ret->k_params = k_params;
1643 ret->destroy_k = destroy_k;
1645 fprintf(stderr, "generated transition %p from %s to %s, params %p\n", ret, ret->initial_state->name, ret->final_state->name, ret->k_params);
1650 void destroy_transition(transition_t *transition)
1654 fprintf(stderr, "freeing transition %p (%s to %s)\n", transition, transition->initial_state->name, transition->final_state->name);
1656 if (transition->destroy_k)
1657 (*transition->destroy_k)(transition->k_params);
1664 We'll be storing the transitions in a list (see Appendix
1665 \ref{sec.lists}), so we also define a convenience function for
1667 <<destroy transition list>>=
1668 void destroy_transition_list(list_t *transition_list)
1670 transition_list = head(transition_list);
1671 while (transition_list != NULL)
1672 destroy_transition((transition_t *) pop(&transition_list));
1677 \subsection{Graphviz output}
1678 \label{sec.graphviz_output}
1680 <<print state graph>>=
1681 void print_state_graph(FILE *file, list_t *state_list, list_t *transition_list)
1683 state_list = head(state_list);
1684 transition_list = head(transition_list);
1685 fprintf(file, "digraph states {\n node [shape=record,color=black,fillcolor=white,style=filled];\n\n");
1687 fprintf(file, " n%d [label=\"%s\"]\n", state_index(state_list, S(state_list)->name), S(state_list)->name);
1688 if (state_list->next == NULL) break;
1689 state_list = state_list->next;
1691 fprintf(file, "\n");
1692 while (transition_list != NULL) {
1693 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));
1694 transition_list = transition_list->next;
1696 fprintf(file, "}\n");
1700 \subsection{Domain model and group handling}
1702 <<group functions>>=
1703 <<crosslink groups>>
1704 <<get active groups>>
1707 Scan through a list of states and construct an array of all the
1709 <<get active groups>>=
1710 void get_active_groups(list_t *state_list,
1711 int *pNum_active_groups,
1712 one_dim_fn_t ***pPer_group_handlers,
1713 one_dim_fn_t ***pPer_group_inverse_handlers,
1714 void ***pPer_group_data /* pointer to array of pointers to tension_handler_data_t */)
1716 void **active_groups=NULL;
1717 one_dim_fn_t *handler, *old_handler, *inverse_handler, *old_inverse_handler;
1718 one_dim_fn_t **per_group_handlers=NULL, **per_group_inverse_handlers=NULL;
1719 void **per_group_data=NULL; /* array of pointers to tension_handler_data_t */
1720 int i, j, num_domains, group, old_group, num_active_groups=0;
1721 list_t *active_states_list=NULL;
1722 tension_handler_data_t *tdata=NULL;
1725 /* get all the active groups in a list */
1726 state_list = head(state_list);
1728 fprintf(stderr, "%s:\t", __FUNCTION__);
1729 list_print(stderr, state_list, "states");
1731 while (state_list != NULL) {
1732 state = S(state_list);
1733 handler = state->tension_handler;
1734 inverse_handler = state->inverse_tension_handler;
1735 group = state->tension_group;
1736 num_domains = state->num_domains;
1737 if (list_index(active_states_list, state) == -1) {
1738 /* we haven't added this state (or it's associated group) yet */
1739 if (num_domains > 0 && handler != NULL) { /* add the state */
1740 num_active_groups++;
1741 push(&active_states_list, state, 1);
1743 fprintf(stderr, "%s:\tactive %s state (%p) with %d domain(s)\n", __FUNCTION__, state->name, state, state->num_domains);
1745 for (i=0; i < state->num_grouped_states; i++) {
1746 if (state->grouped_states[i]->num_domains > 0) {
1747 /* add this related (and active) state now */
1748 assert(state->grouped_states[i]->tension_handler == handler);
1749 assert(state->grouped_states[i]->tension_group == group);
1750 push(&active_states_list, state->grouped_states[i], 1);
1752 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);
1756 } /* else empty state or NULL handler */
1757 } /* else state already added */
1758 state_list = state_list->next;
1761 fprintf(stderr, "%s:\t", __FUNCTION__);
1762 list_print(stderr, active_states_list, "active states");
1765 assert(num_active_groups <= list_length(active_states_list));
1766 assert(num_active_groups > 0);
1768 /* allocate the arrays we'll be returning */
1769 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1770 assert(per_group_handlers != NULL);
1771 per_group_inverse_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1772 assert(per_group_inverse_handlers != NULL);
1773 per_group_data = (void **)malloc(sizeof(void *)*num_active_groups);
1774 assert(per_group_data != NULL);
1776 fprintf(stderr, "%s:\tallocating per group data %p:", __FUNCTION__, per_group_data);
1778 for (i=0; i<num_active_groups; i++) { /* we might want to recycle any group data structs allocated in the previous round */
1779 per_group_data[i] = malloc(sizeof(tension_handler_data_t));
1780 assert(per_group_data[i] != NULL);
1782 fprintf(stderr, " %d,%p", i, per_group_data[i]);
1786 fprintf(stderr, "\n");
1789 /* populate the arrays we'll be returning */
1790 active_states_list = head(active_states_list);
1792 old_inverse_handler = NULL;
1793 j = -1; /* j is the current active _group_ index */
1794 while (active_states_list != NULL) {
1795 state = (state_t *) pop(&active_states_list);
1796 handler = state->tension_handler;
1797 inverse_handler = state->inverse_tension_handler;
1798 group = state->tension_group;
1799 num_domains = state->num_domains;
1800 if (handler != old_handler || inverse_handler != old_inverse_handler || group != old_group) {
1801 j++; /* move to the next group */
1802 tdata = (tension_handler_data_t *) per_group_data[j];
1803 per_group_handlers[j] = handler;
1804 per_group_inverse_handlers[j] = inverse_handler;
1805 tdata->group_tension_params = NULL;
1807 tdata->persist = NULL; /* we might want to grab any appropriate data from the old persist here... */
1810 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);
1812 for (i=0; i<num_domains; i++) { /* add num_domains copies of the tension params */
1813 push(&tdata->group_tension_params, state->tension_params, 1);
1816 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);
1818 old_handler = handler;
1819 old_inverse_handler = inverse_handler;
1823 /* free old groups */
1824 if (*pPer_group_handlers != NULL)
1825 free(*pPer_group_handlers);
1826 if (*pPer_group_inverse_handlers != NULL)
1827 free(*pPer_group_inverse_handlers);
1828 if (*pPer_group_data != NULL) {
1829 for (i=0; i<*pNum_active_groups; i++)
1830 free((*pPer_group_data)[i]);
1831 free(*pPer_group_data);
1834 /* set new groups */
1836 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]);
1838 *pNum_active_groups = num_active_groups;
1839 *pPer_group_handlers = per_group_handlers;
1840 *pPer_group_inverse_handlers = per_group_inverse_handlers;
1841 *pPer_group_data = per_group_data;
1845 <<crosslink groups>>=
1846 void crosslink_groups(list_t *state_groups, list_t *transition_list)
1848 list_t *list, *out_trans=NULL, *associates=NULL;
1849 one_dim_fn_t *handler;
1852 state_groups = head(state_groups);
1853 while (state_groups != NULL) {
1854 /* find transitions leaving this state */
1856 fprintf(stderr, "finding transitions leaving %s\n", S(state_groups)->name);
1858 list = head(transition_list);
1859 while (list != NULL) {
1860 if (T(list)->initial_state == S(state_groups) && list_index(associates, T(list)->initial_state) != -1) {
1861 push(&out_trans, T(list), 1);
1865 n = list_length(out_trans);
1866 S(state_groups)->num_out_trans = n;
1867 S(state_groups)->out_trans = (transition_t **)malloc(sizeof(transition_t *)*n);
1868 assert(S(state_groups)->out_trans != NULL);
1869 for (i=0; i<n; i++) {
1870 S(state_groups)->out_trans[i] = (transition_t *)pop(&out_trans);
1873 /* find states grouped with this state */
1875 fprintf(stderr, "finding states grouped with %s\n", S(state_groups)->name);
1877 handler = S(state_groups)->tension_handler;
1878 group = S(state_groups)->tension_group;
1879 list = head(state_groups);
1880 while (list != NULL) {
1881 if (S(list)->tension_handler == handler && S(list)->tension_group == group && list != state_groups) {
1882 push(&associates, S(list), 1);
1886 n = list_length(associates);
1887 S(state_groups)->num_grouped_states = n;
1888 S(state_groups)->grouped_states = (state_t **)malloc(sizeof(state_t *)*n);
1889 assert(S(state_groups)->grouped_states != NULL);
1890 for (i=0; i<n; i++) {
1891 S(state_groups)->grouped_states[i] = (state_t *)pop(&associates);
1893 state_groups = state_groups->next;
1899 \section{String parsing}
1901 For handling command line arguments, we need parse delimited strings (e.g.\ [["param1,param2,param3..."]]).
1902 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1903 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1904 need to take care of parsing those parameters themselves.
1905 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1911 <<parse definitions>>
1912 <<parse declarations>>
1916 <<parse module makefile lines>>=
1917 NW_SAWSIM_MODS += parse
1918 CHECK_BINS += check_parse
1922 <<parse definitions>>=
1923 #define SEP ',' /* argument separator character */
1926 <<parse declarations>>=
1927 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1928 int *num, char ***string_array);
1929 extern void free_parsed_list(int num, char **string_array);
1932 [[parse_list_string]] allocates memory, don't forget to free it
1933 afterward with [[free_parsed_list]]. It does not alter the original.
1935 The string may start off with a [[deeper]] character
1936 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1937 [[x,y]], where the pointer is one character in on the copied string.
1938 However, when we go to free the memory, we need a pointer to the
1939 beginning of the string. In order to accommodate this for a string
1940 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1941 the first $N$ elements point to the separated arguments, and let the
1942 last element point to the start of the copied string regardless of
1944 <<parse delimited list string functions>>=
1945 /* TODO, split out into parse.hc */
1946 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1949 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1950 if (string[i] == deeper) {depth++;}
1951 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1957 void parse_list_string(char *string, char sep, char deeper, char shallower,
1958 int *num, char ***string_array)
1960 char *str=NULL, **ret=NULL;
1962 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
1964 *string_array = NULL;
1967 /* make a copy of the string, so we don't change the original */
1968 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1969 assert(str != NULL);
1970 strcpy(str, string); /* we know str is long enough */
1971 /* count the number of regions, so we can allocate pointers to them */
1974 n++; i++; /* move on to next argument */
1975 i += next_delim_index(str+i, sep, deeper, shallower);
1976 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1977 } while (str[i] != '\0');
1978 ret = (char **)malloc(sizeof(char *)*(n+1));
1979 assert(ret != NULL);
1980 /* replace the separators with '\0' & assign pointers */
1981 ret[n] = str; /* point to the front of the copied string */
1984 for(i=1; i<n; i++) {
1985 j += next_delim_index(str+j, sep, deeper, shallower);
1987 ret[i] = str+j; /* point to the separated arguments */
1989 /* strip off bounding braces, if any */
1990 for(i=0; i<n; i++) {
1991 if (ret[i][0] == deeper) {
1995 if (ret[i][strlen(ret[i])-1] == shallower)
1996 ret[i][strlen(ret[i])-1] = '\0';
1999 *string_array = ret;
2002 void free_parsed_list(int num, char **string_array)
2005 /* frees all items, since they were allocated as a single string*/
2006 free(string_array[num]);
2007 /* and free the array of pointers */
2013 \subsection{Parsing implementation}
2019 <<parse delimited list string functions>>
2023 #include <assert.h> /* assert() */
2024 #include <stdlib.h> /* NULL */
2025 #include <stdio.h> /* fprintf(), stdout *//*!!*/
2026 #include <string.h> /* strlen() */
2030 \subsection{Parsing unit tests}
2032 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2034 <<parse check includes>>
2035 <<string comparison definition and globals>>
2036 <<check relative error>>
2037 <<parse test suite>>
2038 <<main check program>>
2041 <<parse check includes>>=
2042 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2043 #include <stdio.h> /* printf() */
2044 #include <assert.h> /* assert() */
2045 #include <string.h> /* strlen() */
2050 <<parse test suite>>=
2051 <<parse list string tests>>
2052 <<parse suite definition>>
2055 <<parse suite definition>>=
2056 Suite *test_suite (void)
2058 Suite *s = suite_create ("k model");
2059 <<parse list string test case defs>>
2061 <<parse list string test case add>>
2066 <<parse list string tests>>=
2069 START_TEST(test_next_delim_index)
2071 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
2072 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
2073 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
2074 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
2075 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
2079 START_TEST(test_parse_list_null)
2083 parse_list_string(NULL, SEP, '{', '}',
2084 &num_param_args, ¶m_args);
2085 fail_unless(num_param_args == 0, NULL);
2086 fail_unless(param_args == NULL, NULL);
2089 START_TEST(test_parse_list_single_simple)
2093 parse_list_string("arg", SEP, '{', '}',
2094 &num_param_args, ¶m_args);
2095 fail_unless(num_param_args == 1, NULL);
2096 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
2099 START_TEST(test_parse_list_single_compound)
2103 parse_list_string("{x,y,z}", SEP, '{', '}',
2104 &num_param_args, ¶m_args);
2105 fail_unless(num_param_args == 1, NULL);
2106 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
2109 START_TEST(test_parse_list_double_simple)
2113 parse_list_string("abc,def", SEP, '{', '}',
2114 &num_param_args, ¶m_args);
2115 fail_unless(num_param_args == 2, NULL);
2116 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2117 fail_unless(STRMATCH(param_args[1],"def"), NULL);
2120 START_TEST(test_parse_list_compound)
2124 parse_list_string("abc,{def,ghi}", SEP, '{', '}',
2125 &num_param_args, ¶m_args);
2126 fail_unless(num_param_args == 2, NULL);
2127 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2128 fail_unless(STRMATCH(param_args[1],"def,ghi"), NULL);
2132 <<parse list string test case defs>>=
2133 TCase *tc_parse_list_string = tcase_create("parse list string");
2135 <<parse list string test case add>>=
2136 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
2137 tcase_add_test(tc_parse_list_string, test_parse_list_null);
2138 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
2139 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
2140 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
2141 tcase_add_test(tc_parse_list_string, test_parse_list_compound);
2142 suite_add_tcase(s, tc_parse_list_string);
2146 \section{Unit tests}
2148 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2155 <<check relative error>>
2158 <<main check program>>
2170 <<determine dt tests>>
2172 <<does domain transition tests>>
2173 <<randomly unfold domains tests>>
2174 <<suite definition>>
2177 <<suite definition>>=
2178 Suite *test_suite (void)
2180 Suite *s = suite_create ("sawsim");
2181 <<F test case defs>>
2182 <<determine dt test case defs>>
2183 <<happens test case defs>>
2184 <<does domain transition test case defs>>
2185 <<randomly unfold domains test case defs>>
2188 <<determine dt test case add>>
2189 <<happens test case add>>
2190 <<does domain transition test case add>>
2191 <<randomly unfold domains test case add>>
2194 tcase_add_checked_fixture(tc_strip_address,
2195 setup_strip_address,
2196 teardown_strip_address);
2202 <<main check program>>=
2206 Suite *s = test_suite();
2207 SRunner *sr = srunner_create(s);
2208 srunner_run_all(sr, CK_ENV);
2209 number_failed = srunner_ntests_failed(sr);
2211 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
2215 \subsection{F tests}
2217 <<worm-like chain tests>>
2219 <<F test case defs>>=
2220 <<worm-like chain test case def>>
2222 <<F test case add>>=
2223 <<worm-like chain test case add>>
2226 <<worm-like chain tests>>=
2227 <<worm-like chain function>>
2229 START_TEST(test_wlc_at_zero)
2231 double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
2232 fail_unless(abs(wlc(x, T, p, L)) < lim, \
2233 "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
2234 x, T, p, L, abs(wlc(x, T, p, L)), lim);
2238 START_TEST(test_wlc_at_half)
2240 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
2241 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
2242 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
2244 * linear term = x/L = 0.5
2245 * nonlinear + linear = 0.75 + 0.5 = 1.25
2246 * wlc = 10e21*1.25 = 12.5e21
2248 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
2249 "wlc(%g, %g, %g, %g) = %g != %g",
2250 x, T, p, L, wlc(x, T, p, L), 12.5e21);
2254 <<worm-like chain test case def>>=
2255 TCase *tc_wlc = tcase_create("WLC");
2258 <<worm-like chain test case add>>=
2259 tcase_add_test(tc_wlc, test_wlc_at_zero);
2260 tcase_add_test(tc_wlc, test_wlc_at_half);
2261 suite_add_tcase(s, tc_wlc);
2264 \subsection{Model tests}
2266 Check the searching with [[linear_k]].
2267 Check overwhelming force treatment with the heavyside-esque step [[?]].
2269 Hmm, I'm not really sure what I was doing here...
2270 <<determine dt tests>>=
2271 double linear_k(double F, environment_t *env, void *params)
2273 double Fnot = *(double *)params;
2277 START_TEST(test_determine_dt_linear_k)
2280 transition_t unfold={0};
2281 environment_t env = {0};
2282 double F[]={0,1,2,3,4,5,6};
2283 double dt_max=3.0, Fnot=3.0, x=1.0, max_p=0.1, max_dF=0.1, v=1.0;
2284 list_t *state_list=NULL, *transition_list=NULL;
2287 folded.tension_handler = &hooke_handler;
2288 folded.num_domains = 1;
2289 unfold.initial_state = &folded;
2290 unfold.k = linear_k;
2291 unfold.k_params = &Fnot;
2292 push(&state_list, &folded, 1);
2293 push(&transition_list, &unfold, 1);
2294 for( i=0; i < sizeof(F)/sizeof(double); i++) {
2295 //fail_unless(determine_dt(state_list, transition_list, &env, x, max_p, max_dF, dt_max, v)==1, NULL);
2300 <<determine dt test case defs>>=
2301 TCase *tc_determine_dt = tcase_create("Determine dt");
2303 <<determine dt test case add>>=
2304 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
2305 suite_add_tcase(s, tc_determine_dt);
2310 <<happens test case defs>>=
2312 <<happens test case add>>=
2315 <<does domain transition tests>>=
2317 <<does domain transition test case defs>>=
2319 <<does domain transition test case add>>=
2322 <<randomly unfold domains tests>>=
2324 <<randomly unfold domains test case defs>>=
2326 <<randomly unfold domains test case add>>=
2330 \section{Balancing group extensions}
2331 \label{sec.tension_balance}
2333 For a given total extension $x$ (of the piezo), the various domain
2334 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
2335 amounts, and we need to tweak the portion that each extends to balance
2336 the tension amongst the active groups. Since the tension balancing is
2337 separable from the bulk of the simulation, we break it out into a
2338 separate module. The interface is defined in [[tension_balance.h]],
2339 the implementation in [[tension_balance.c]], and the unit testing in
2340 [[check_tension_balance.c]]
2342 <<tension-balance.h>>=
2343 #ifndef TENSION_BALANCE_H
2344 #define TENSION_BALANCE_H
2346 <<tension balance function declaration>>
2347 #endif /* TENSION_BALANCE_H */
2350 <<tension balance functions>>=
2351 <<one dimensional bracket>>
2352 <<one dimensional solve>>
2354 <<group stiffness function>>
2355 <<chain stiffness function>>
2356 <<full chain stiffness function>>
2357 <<tension balance function>>
2360 <<tension balance module makefile lines>>=
2361 NW_SAWSIM_MODS += tension_balance
2362 CHECK_BINS += check_tension_balance
2363 check_tension_balance_MODS =
2366 The entire force balancing problem reduces to a solving two nested
2367 one-dimensional root-finding problems. First we define the one
2368 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
2369 populated group). $x(x_0)$ is also strictly monotonically increasing,
2370 so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
2371 stores the last successful value of $x$, since we expect to be taking
2372 small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
2373 indicates that the groups have changed.
2374 <<tension balance function declaration>>=
2375 double tension_balance(int num_tension_groups,
2376 one_dim_fn_t **tension_handlers,
2377 one_dim_fn_t **inverse_tension_handlers,
2378 void **params, /* array of pointers to tension_handler_data_t */
2383 <<tension balance function>>=
2384 double tension_balance(int num_tension_groups,
2385 one_dim_fn_t **tension_handlers,
2386 one_dim_fn_t **inverse_tension_handlers,
2387 void **params, /* array of pointers to tension_handler_data_t */
2392 static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL,NULL};
2393 double F, xo_guess, xo, lb, ub;
2394 double min_relx=1e-6, min_rely=1e-6;
2395 int max_steps=200, i;
2397 assert(num_tension_groups > 0);
2398 assert(tension_handlers != NULL);
2399 assert(params != NULL);
2400 assert(num_tension_groups > 0);
2402 if (last_x == -1 || x_xo_data.xi == NULL) { /* new group setup, reset x_xo_data */
2403 /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
2404 if (x_xo_data.xi != NULL)
2406 /* construct new x_xo_data */
2407 x_xo_data.n_groups = num_tension_groups;
2408 x_xo_data.tension_handler = tension_handlers;
2409 x_xo_data.inverse_tension_handler = inverse_tension_handlers;
2410 x_xo_data.handler_data = params;
2412 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);
2414 x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
2415 assert(x_xo_data.xi != NULL);
2416 for (i=0; i<num_tension_groups; i++)
2417 x_xo_data.xi[i] = last_x;
2420 if (num_tension_groups == 1) { /* only one group, no balancing required */
2422 x_xo_data.xi[0] = xo;
2424 fprintf(stderr, "%s:\tonly one active group, no balancing required\n", __FUNCTION__);
2428 fprintf(stderr, "%s:\tbalancing for x = %g with", __FUNCTION__, x);
2429 for(i=0; i<num_tension_groups; i++) fprintf(stderr, "\tx[%d]=%g", i, x_xo_data.xi[i]);
2430 fprintf(stderr, "\n");
2432 if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
2433 if (x == 0) xo_guess = length_scale;
2434 else xo_guess = x/num_tension_groups;
2436 fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
2438 oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
2439 } else { /* work off the last balanced point */
2441 fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
2444 lb = x_xo_data.xi[0];
2445 ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */
2446 } else if (x < last_x) {
2447 lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
2448 ub = x_xo_data.xi[0];
2449 } else { /* x == last_x */
2451 fprintf(stderr, "not moving\n");
2453 lb = x_xo_data.xi[0];
2454 ub = x_xo_data.xi[0];
2458 fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
2459 __FUNCTION__, x, lb, ub);
2461 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
2462 F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
2464 if (num_tension_groups > 1) {
2465 fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
2466 for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
2467 fprintf(stderr, "\n");
2472 F = (*tension_handlers[0])(xo, params[0]);
2474 if (stiffness != NULL) {
2475 *stiffness = chain_stiffness(&x_xo_data, min_relx);
2476 if (*stiffness == 0.0) { /* re-calculate based on full chain */
2477 *stiffness = full_chain_stiffness(num_tension_groups, tension_handlers,
2478 inverse_tension_handlers, params,
2479 last_x, x, min_relx, F);
2487 <<tension balance internal definitions>>=
2488 <<x of x0 definitions>>
2491 <<x of x0 definitions>>=
2492 typedef struct x_of_xo_data_struct {
2494 one_dim_fn_t **tension_handler; /* array of fn pointers */
2495 one_dim_fn_t **inverse_tension_handler; /* inverse fn pointers */
2496 void **handler_data; /* array of pointers to tension_handler_data_t */
2497 double *xi; /* array of group extensions */
2502 double x_of_xo(double xo, void *pdata)
2504 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
2505 double F, x=0, xi, xi_guess, lb, ub, slope;
2507 double min_relx=1e-6, min_rely=1e-6;
2509 assert(data->n_groups > 0);
2512 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);
2515 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
2517 fprintf(stderr, "%s:\tfinding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
2521 if (data->xi) data->xi[0] = xo;
2522 for (i=1; i < data->n_groups; i++) {
2523 if (data->inverse_tension_handler[i] != NULL) {
2524 xi = (*data->inverse_tension_handler[i])(F, data->handler_data[i]);
2525 } else { /* invert numerically */
2526 if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
2527 else xi_guess = data->xi[i];
2529 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]);
2531 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
2532 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
2537 fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
2541 fprintf(stderr, "%s:\tfound x(x0=%g) = %g\n", __FUNCTION__, xo, x);
2547 Determine the stiffness of a single tension group by taking a small
2548 step [[dx]] back from the position [[x]] for which we want the
2549 stiffness. The touchy part is determining a good [[dx]] and ensuring
2550 the whole interval is on [[x>=0]].
2551 <<group stiffness function>>=
2552 double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double x, double relx)
2554 double dx, xl, F, Fl, stiffness;
2556 assert(relx > 0 && relx < 1);
2557 if (x == 0 || relx == 0) {
2558 dx = 1e-12; /* HACK, how to get length scale? */
2568 F = (*tension_handler)(x, handler_data);
2569 Fl = (*tension_handler)(xl, handler_data);
2570 stiffness = (F-Fl)/dx;
2571 assert(stiffness >= 0);
2576 Attempt to calculate the chain stiffness by adding group stiffnesses
2577 as springs in series. In the case where a group has undefined
2578 stiffness (e.g. piston tension, Section \ref{sec.piston_tension}),
2579 this algorithm breaks down. In that case, [[tension_balance()]]
2580 (\ref{sec.tension_balance}) should use [[full_chain_stiffness()]]
2581 which uses the full chain's $F(x)$ rather than that of the individual
2582 domains'. We attempt the series approach first because, lacking the
2583 numerical inversion inside [[tension_balance()]], it is both faster
2584 and more accurate than the full-chain derivative.
2585 <<chain stiffness function>>=
2586 double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx)
2588 double stiffness, gstiff;
2590 /* add group stiffnesses in series */
2591 for (i=0; i < x_xo_data->n_groups; i++) {
2593 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);
2596 gstiff = group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx);
2599 stiffness += 1.0/gstiff;
2601 stiffness = 1.0 / stiffness;
2607 Determine the chain stiffness using only [[tension_balance()]]. This
2608 works around difficulties with tension models that have undefined
2609 stiffness (see discussion for [[chain_stiffness()]] above).
2610 <<full chain stiffness function>>=
2611 double full_chain_stiffness(int num_tension_groups,
2612 one_dim_fn_t **tension_handlers,
2613 one_dim_fn_t **inverse_tension_handlers,
2614 void **params, /* array of pointers to tension_handler_data_t */
2618 double F /* already evaluated F(x) */)
2620 double dx, xl, Fl, stiffness;
2623 assert(relx > 0 && relx < 1);
2625 /* Other option for dx that we ignore because of our poor tension_balance()
2626 * resolution (i.e. bad slopes if you zoom in too close):
2627 * if (last_x != -1.0 && last_x != x) // last_x limited
2628 * dx fabs(x-last_x);
2631 * if (1==1) { // constant max-value limited
2633 * dx = MIN(dx, dx_p);
2635 * if (x != 0 && relx != 0) { // relx limited
2637 * dx = MIN(dx, dx_p);
2639 * TODO, determine which of (min_relx, min_rely, max_steps) is actually
2640 * limiting tension_balance.
2642 dx = 1e-12; /* HACK, how to get length scale? */
2646 Fl = tension_balance(num_tension_groups, tension_handlers,
2647 inverse_tension_handlers, params,
2653 F = tension_balance(num_tension_groups, tension_handlers,
2654 inverse_tension_handlers, params,
2657 stiffness = (F-Fl)/dx;
2658 assert(stiffness >= 0);
2664 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
2665 number of steps for monotonically increasing $f(x)$. Simple
2666 bisection, so it's robust and converges fairly quickly. You can set
2667 the maximum number of steps to take, but if you have enough processor
2668 power, it's better to limit your precision with the relative limits
2669 [[min_relx]] and [[y]]. These ensure that the width of the bracket is
2670 small on the length scale of it's larger side. Note that \emph{both}
2671 relative limits must be satisfied for the search to stop.
2672 <<one dimensional function definition>>=
2673 /* equivalent to gsl_func_t */
2674 typedef double one_dim_fn_t(double x, void *params);
2676 <<one dimensional solve declaration>>=
2677 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2678 double min_relx, double min_rely, int max_steps,
2682 <<one dimensional solve>>=
2683 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2684 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2685 double min_relx, double min_rely, int max_steps,
2688 double yx, ylb, yub, x;
2691 ylb = (*f)(lb, params);
2692 yub = (*f)(ub, params);
2694 /* check some simple cases */
2696 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bracketed */
2697 else return lb; /* any x on the range [lb, ub] would work */
2699 if (ylb == y) { x=lb; yx=ylb; goto end; }
2700 if (yub == y) { x=ub; yx=yub; goto end; }
2703 fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2705 assert(ylb < y); assert(yub > y);
2706 for (i=0; i<max_steps; i++) {
2708 yx = (*f)(x, params);
2710 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);
2712 if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2713 else if (yx > y) { ub=x; yub=yx; }
2714 else /* yx < y */ { lb=x; ylb=yx; }
2719 fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2725 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2726 Generate bracketing $x$ values through bisection or exponential growth.
2727 <<one dimensional bracket declaration>>=
2728 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2731 <<one dimensional bracket>>=
2732 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2734 double yx, step, x=xguess;
2737 fprintf(stderr, "%s:\tbracketing %g (params %p)\n", __FUNCTION__, y, params);
2739 yx = (*f)(x, params);
2741 fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2748 fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2752 if (x == 0) x = length_scale/2.0; /* HACK */
2755 yx = (*f)(x, params);
2757 fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2762 fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2765 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2769 \subsection{Balancing implementation}
2771 <<tension-balance.c>>=
2774 <<tension balance includes>>
2775 #include "tension_balance.h"
2777 double length_scale = 1e-10; /* HACK */
2779 <<tension balance internal definitions>>
2780 <<tension balance functions>>
2783 <<tension balance includes>>=
2784 #include <assert.h> /* assert() */
2785 #include <stdlib.h> /* NULL */
2786 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2787 #include <stdio.h> /* fprintf(), stdout */
2791 \subsection{Balancing unit tests}
2793 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2794 <<check-tension-balance.c>>=
2795 <<tension balance check includes>>
2796 <<tension balance test suite>>
2797 <<main check program>>
2800 <<tension balance check includes>>=
2801 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2802 #include <assert.h> /* assert() */
2805 #include "tension_balance.h"
2808 <<tension balance test suite>>=
2809 <<tension balance function tests>>
2810 <<tension balance suite definition>>
2813 <<tension balance suite definition>>=
2814 Suite *test_suite (void)
2816 Suite *s = suite_create ("tension balance");
2817 <<tension balance function test case defs>>
2819 <<tension balance function test case adds>>
2824 <<tension balance function tests>>=
2825 <<check relative error>>
2827 double hooke(double x, void *pK)
2830 return *((double*)pK) * x;
2833 START_TEST(test_single_function)
2835 double k=5, x=3, last_x=2.0, F;
2836 one_dim_fn_t *handlers[] = {&hooke};
2837 one_dim_fn_t *inverse_handlers[] = {NULL};
2838 void *data[] = {&k};
2839 F = tension_balance(1, handlers, inverse_handlers, data, last_x, x, NULL);
2840 fail_unless(F = k*x, NULL);
2845 We can also test balancing two springs with different spring constants.
2846 The analytic solution for a general number of springs is given in Appendix \ref{sec.math_hooke}.
2847 <<tension balance function tests>>=
2848 START_TEST(test_double_hooke)
2850 double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
2851 one_dim_fn_t *handlers[] = {&hooke, &hooke};
2852 one_dim_fn_t *inverse_handlers[] = {NULL, NULL};
2853 void *data[] = {&k1, &k2};
2854 F = tension_balance(2, handlers, inverse_handlers, data, last_x, x, NULL);
2858 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2859 //CHECK_ERR(1e-6, x1e, xi[0]);
2860 //CHECK_ERR(1e-6, x2e, xi[1]);
2861 CHECK_ERR(1e-6, Fe, F);
2865 TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
2867 <<tension balance function test case defs>>=
2868 TCase *tc_tbfunc = tcase_create("tension balance function");
2871 <<tension balance function test case adds>>=
2872 tcase_add_test(tc_tbfunc, test_single_function);
2873 tcase_add_test(tc_tbfunc, test_double_hooke);
2874 suite_add_tcase(s, tc_tbfunc);
2880 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2881 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2882 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2888 <<list definitions>>
2889 <<list declarations>>
2893 <<list declarations>>=
2894 <<head and tail declarations>>
2895 <<list length declaration>>
2896 <<push declaration>>
2898 <<list destroy declaration>>
2899 <<list to array declaration>>
2900 <<move declaration>>
2901 <<sort declaration>>
2902 <<list index declaration>>
2903 <<list element declaration>>
2904 <<unique declaration>>
2905 <<list print declaration>>
2909 <<create and destroy node>>
2924 <<list module makefile lines>>=
2925 NW_SAWSIM_MODS += list
2926 CHECK_BINS += check_list
2930 <<list definitions>>=
2931 typedef struct list_struct {
2932 struct list_struct *next;
2933 struct list_struct *prev;
2937 <<comparison function typedef>>
2940 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2941 <<head and tail declarations>>=
2942 list_t *head(list_t *list);
2943 list_t *tail(list_t *list);
2946 list_t *head(list_t *list)
2950 while (list->prev != NULL) {
2956 list_t *tail(list_t *list)
2960 while (list->next != NULL) {
2967 <<list length declaration>>=
2968 int list_length(list_t *list);
2971 int list_length(list_t *list)
2978 while (list->next != NULL) {
2986 [[push]] inserts a node relative to the current position in the list
2987 without changing the current position.
2988 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2989 so we set it to that of the pushed domain.
2990 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2991 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2992 <<push declaration>>=
2993 void push(list_t **pList, void *data, int below);
2996 void push(list_t **pList, void *data, int below)
2998 list_t *list, *new_node;
2999 assert(pList != NULL);
3001 new_node = create_node(data);
3004 else if (below==0) { /* insert above */
3005 if (list->prev != NULL)
3006 list->prev->next = new_node;
3007 new_node->prev = list->prev;
3008 list->prev = new_node;
3009 new_node->next = list;
3010 } else { /* insert below */
3011 if (list->next != NULL)
3012 list->next->prev = new_node;
3013 new_node->next = list->next;
3014 list->next = new_node;
3015 new_node->prev = list;
3020 [[pop]] removes the current domain node, moving the current position
3021 to the node after it, unless that node is [[NULL]], in which case move
3022 the current position to the node before it.
3023 <<pop declaration>>=
3024 void *pop(list_t **pList);
3027 void *pop(list_t **pList)
3029 list_t *list, *popped;
3031 assert(pList != NULL);
3033 assert(list != NULL); /* not an empty list */
3035 /* bypass the popped node */
3036 if (list->prev != NULL)
3037 list->prev->next = list->next;
3038 if (list->next != NULL) {
3039 list->next->prev = list->prev;
3040 *pList = list->next;
3042 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
3044 destroy_node(popped);
3049 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
3050 If the cleanup function is [[NULL]], no cleanup function is called.
3051 The nodes are not popped in any particular order.
3052 <<list destroy declaration>>=
3053 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
3056 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
3060 assert(pList != NULL);
3063 assert(list != NULL); /* not an empty list */
3064 while (list != NULL) {
3066 if (destroy != NULL)
3072 [[list_to_array]] converts a list to an array of pointers, preserving order.
3073 <<list to array declaration>>=
3074 void list_to_array(list_t *list, int *length, void ***array);
3077 void list_to_array(list_t *list, int *pLength, void ***pArray)
3081 assert(list != NULL);
3082 assert(pLength != NULL);
3083 assert(pArray != NULL);
3085 length = list_length(list);
3086 array = (void **)malloc(sizeof(void *)*length);
3087 assert(array != NULL);
3090 while (list != NULL) {
3091 array[i++] = list->d;
3099 [[move]] moves the current element along the list in either direction.
3100 <<move declaration>>=
3101 void move(list_t **pList, int downstream);
3104 void move(list_t **pList, int downstream)
3106 list_t *list, *flipper;
3108 assert(pList != NULL);
3109 list = *pList; /* a>B>c>d */
3110 assert(list != NULL); /* not an empty list */
3111 if (downstream == 0)
3112 flipper = list->prev; /* flipper is A */
3114 flipper = list->next; /* flipper is C */
3115 /* ensure there is somewhere to go in the direction we're heading */
3116 assert(flipper != NULL);
3117 /* remove the current element */
3118 data = pop(&list); /* data=B, a>C>d */
3119 /* and put it back in in the correct spot */
3121 if (downstream == 0) {
3122 push(&list, data, 0); /* b>A>c>d */
3123 list = list->prev; /* B>a>c>d */
3125 push(&list, data, 1); /* a>C>b>d */
3126 list = list->next; /* a>c>B>d */
3132 [[sort]] sorts a list from smallest to largest according to a user
3133 specified comparison.
3134 <<comparison function typedef>>=
3135 typedef int (comparison_fn_t)(void *A, void *B);
3138 <<int comparison function>>=
3139 int int_comparison(void *A, void *B)
3144 if (a > b) return 1;
3145 else if (a == b) return 0;
3150 <<sort declaration>>=
3151 void sort(list_t **pList, comparison_fn_t *comp);
3153 Here we use the bubble sort algorithm.
3155 void sort(list_t **pList, comparison_fn_t *comp)
3159 assert(pList != NULL);
3161 assert(list != NULL); /* not an empty list */
3162 while (stable == 0) {
3165 while (list->next != NULL) {
3166 if ((*comp)(list->d, list->next->d) > 0) {
3168 move(&list, 1 /* downstream */);
3178 [[list_index]] finds the location of [[data]] in [[list]]. Returns
3179 [[-1]] if [[data]] is not in [[list]].
3180 <<list index declaration>>=
3181 int list_index(list_t *list, void *data);
3185 int list_index(list_t *list, void *data)
3189 while (list != NULL) {
3190 if (list->d == data) return i;
3199 [[list_element]] returns the data bound to the $i^\text{th}$ node of the list.
3200 <<list element declaration>>=
3201 void *list_element(list_t *list, int i);
3204 void *list_element(list_t *list, int i)
3208 while (list != NULL) {
3209 if (i == j) return list->d;
3219 [[unique]] removes duplicates from a sorted list according to a user
3220 specified comparison. Don't do this unless you have other pointers
3221 any dynamically allocated data in your list, because [[unique]] just
3222 drops any non-unique data without freeing it.
3223 <<unique declaration>>=
3224 void unique(list_t **pList, comparison_fn_t *comp);
3227 void unique(list_t **pList, comparison_fn_t *comp)
3230 assert(pList != NULL);
3232 assert(list != NULL); /* not an empty list */
3234 while (list->next != NULL) {
3235 if ((*comp)(list->d, list->next->d) == 0) {
3244 [[list_print]] prints a list to a [[FILE *]] stream.
3245 <<list print declaration>>=
3246 void list_print(FILE *file, list_t *list, char *list_name);
3249 void list_print(FILE *file, list_t *list, char *list_name)
3251 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
3253 while (list != NULL) {
3254 fprintf(file, " %p", list->d);
3257 fprintf(file, "\n");
3261 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
3262 <<create and destroy node>>=
3263 list_t *create_node(void *data)
3265 list_t *ret = (list_t *)malloc(sizeof(list_t));
3266 assert(ret != NULL);
3273 void destroy_node(list_t *node)
3279 The user must free the data pointed to by the node on their own.
3281 \subsection{List implementation}
3291 #include <assert.h> /* assert() */
3292 #include <stdlib.h> /* malloc(), free(), rand() */
3293 #include <stdio.h> /* fprintf(), stdout, FILE */
3294 #include "global.h" /* destroy_data_func_t */
3297 \subsection{List unit tests}
3299 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3301 <<list check includes>>
3303 <<main check program>>
3306 <<list check includes>>=
3307 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3308 #include <stdio.h> /* FILE */
3314 <<list test suite>>=
3317 <<list suite definition>>
3320 <<list suite definition>>=
3321 Suite *test_suite (void)
3323 Suite *s = suite_create ("list");
3324 <<push test case defs>>
3325 <<pop test case defs>>
3327 <<push test case adds>>
3328 <<pop test case adds>>
3334 START_TEST(test_push)
3339 push(&list, (void *)i, 1);
3340 fail_unless(list == head(list), NULL);
3341 fail_unless( (int)list->d == 0 );
3342 for (i=0; i<n; i++) {
3343 /* we expect 0, n-1, n-2, ..., 1 */
3346 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
3352 <<push test case defs>>=
3353 TCase *tc_push = tcase_create("push");
3356 <<push test case adds>>=
3357 tcase_add_test(tc_push, test_push);
3358 suite_add_tcase(s, tc_push);
3363 <<pop test case defs>>=
3365 <<pop test case adds>>=
3368 \section{Function string evaluation}
3370 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).
3371 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
3372 The solution is to run a scripting language as a subprocess accessed via pipes.
3373 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
3375 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
3376 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
3377 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.
3378 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
3380 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.
3381 Then you can either statically or dynamically link to those hardcoded functions.
3382 While much less flexible, this approach would be a fairly simple-to-implement fallback.
3384 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
3385 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
3388 #ifndef STRING_EVAL_H
3389 #define STRING_EVAL_H
3391 <<string eval setup declaration>>
3392 <<string eval function declaration>>
3393 <<string eval teardown declaration>>
3394 #endif /* STRING_EVAL_H */
3397 <<string eval module makefile lines>>=
3398 NW_SAWSIM_MODS += string_eval
3399 CHECK_BINS += check_string_eval
3400 check_string_eval_MODS =
3403 For an introduction to POSIX process control, see\\
3404 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
3405 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
3406 and of course, the relavent [[man]] pages.
3408 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
3409 [[execvp]] replaces the calling process' program with a new program.
3410 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
3411 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
3412 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
3413 The new program needs command line arguments, just like it would if you were running it from a shell.
3414 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
3415 with the final array entry being a [[NULL]] pointer.
3417 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
3418 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
3419 <<string eval subprocess definitions>>=
3420 #define SUBPROCESS "python"
3421 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
3422 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
3423 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
3425 The [[i]] option lets Python know that it should run in interactive mode.
3426 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
3427 In interactive mode, python acts on every instruction as soon as it is recieved.
3428 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.
3429 %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.
3431 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
3432 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
3433 [[fork]] returns two copies of the same program, executing the original code.
3434 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
3435 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
3437 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.
3438 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
3439 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
3440 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
3441 <<string eval pipe definitions>>=
3442 #define PIPE_READ 0 /* the end you read from */
3443 #define PIPE_WRITE 1 /* the end you write to */
3445 #define STDIN 0 /* initial index of stdin pair */
3446 #define STDOUT 2 /* initial index of stdout pair */
3449 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
3451 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.
3453 <<string eval setup declaration>>=
3454 extern void string_eval_setup(FILE **pIn, FILE **pOut);
3456 <<string eval setup definition>>=
3457 void string_eval_setup(FILE **pIn, FILE **pOut)
3460 int pfd[4]; /* file descriptors for stdin and stdout */
3462 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
3463 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
3465 pid = fork(); /* split process into two copies */
3467 if (pid == -1) { /* fork error */
3468 perror("fork error");
3470 } else if (pid == 0) { /* child */
3471 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
3472 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
3473 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
3474 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
3475 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
3476 perror("exec error"); /* exec shouldn't return */
3478 } else { /* parent */
3479 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
3480 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
3481 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
3482 if ( *pIn == NULL ) {
3483 perror("fdopen (in)");
3486 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
3487 if ( *pOut == NULL ) {
3488 perror("fdopen (out)");
3495 To use the evaluating subprocess, we just pipe in our command, and read out the results.
3496 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
3497 <<string eval function declaration>>=
3498 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
3500 <<string eval function definition>>=
3501 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
3504 rc = fprintf(in, "%s", input);
3505 assert(rc == strlen(input));
3508 alarm(1); /* set a one second timeout on the read */
3509 assert( fgets(output, buflen, out) != NULL );
3510 alarm(0); /* cancel the timeout */
3511 //fprintf(stderr, "eval: %s ----> %s", input, output);
3514 The [[alarm]] calls set and clear a timeout on the returned output.
3515 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.
3516 This protects against invalid input for which a line of output is not printed to [[stdout]].
3517 Other invalid inputs (e.g.\ those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
3518 If you are getting strange results, check your python code seperately. TODO, better error handling.
3520 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
3521 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
3522 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
3523 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
3524 <<string eval teardown declaration>>=
3525 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
3528 <<string eval teardown definition>>=
3529 void string_eval_teardown(FILE **pIn, FILE **pOut)
3534 /* redirect Python's stderr */
3535 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
3539 assert( fclose(*pIn) == 0 );
3541 assert( fclose(*pOut) == 0 );
3544 /* wait for python to exit */
3546 pid = wait(&stat_loc);
3553 if (WIFEXITED(stat_loc)) {
3554 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
3555 } else if (WIFSIGNALED(stat_loc)) {
3556 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
3561 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
3563 \subsection{String evaluation implementation}
3567 <<string eval includes>>
3568 #include "string_eval.h"
3569 <<string eval internal definitions>>
3570 <<string eval functions>>
3573 <<string eval includes>>=
3574 #include <assert.h> /* assert() */
3575 #include <stdlib.h> /* NULL */
3576 #include <stdio.h> /* fprintf(), stdout, fdopen() */
3577 #include <string.h> /* strlen() */
3578 #include <sys/types.h> /* pid_t */
3579 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
3580 #include <wait.h> /* wait() */
3583 <<string eval internal definitions>>=
3584 <<string eval subprocess definitions>>
3585 <<string eval pipe definitions>>
3588 <<string eval functions>>=
3589 <<string eval setup definition>>
3590 <<string eval function definition>>
3591 <<string eval teardown definition>>
3594 \subsection{String evaluation unit tests}
3596 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3597 <<check-string-eval.c>>=
3598 <<string eval check includes>>
3599 <<string comparison definition and globals>>
3600 <<string eval test suite>>
3601 <<main check program>>
3604 <<string eval check includes>>=
3605 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3606 #include <stdio.h> /* printf() */
3607 #include <assert.h> /* assert() */
3608 #include <string.h> /* strlen() */
3609 #include <signal.h> /* SIGKILL */
3611 #include "string_eval.h"
3614 <<string eval test suite>>=
3615 <<string eval tests>>
3616 <<string eval suite definition>>
3619 <<string eval suite definition>>=
3620 Suite *test_suite (void)
3622 Suite *s = suite_create ("string eval");
3623 <<string eval test case defs>>
3625 <<string eval test case add>>
3630 <<string eval tests>>=
3631 START_TEST(test_setup_teardown)
3634 string_eval_setup(&in, &out);
3635 string_eval_teardown(&in, &out);
3638 START_TEST(test_invalid_command)
3641 char input[80], output[80]={};
3642 string_eval_setup(&in, &out);
3643 sprintf(input, "print ABCDefg\n");
3644 string_eval(in, out, input, 80, output);
3645 string_eval_teardown(&in, &out);
3648 START_TEST(test_simple_eval)
3651 char input[80], output[80]={};
3652 string_eval_setup(&in, &out);
3653 sprintf(input, "print 3+4*5\n");
3654 string_eval(in, out, input, 80, output);
3655 fail_unless(STRMATCH(output,"23\n"), NULL);
3656 string_eval_teardown(&in, &out);
3659 START_TEST(test_multiple_evals)
3662 char input[80], output[80]={};
3663 string_eval_setup(&in, &out);
3664 sprintf(input, "print 3+4*5\n");
3665 string_eval(in, out, input, 80, output);
3666 fail_unless(STRMATCH(output,"23\n"), NULL);
3667 sprintf(input, "print (3**2 + 4**2)**0.5\n");
3668 string_eval(in, out, input, 80, output);
3669 fail_unless(STRMATCH(output,"5.0\n"), NULL);
3670 string_eval_teardown(&in, &out);
3673 START_TEST(test_eval_with_state)
3676 char input[80], output[80]={};
3677 string_eval_setup(&in, &out);
3678 sprintf(input, "print 3+4*5\n");
3679 fprintf(in, "A = 3\n");
3680 sprintf(input, "print A*3\n");
3681 string_eval(in, out, input, 80, output);
3682 fail_unless(STRMATCH(output,"9\n"), NULL);
3683 string_eval_teardown(&in, &out);
3687 <<string eval test case defs>>=
3688 TCase *tc_string_eval = tcase_create("string_eval");
3690 <<string eval test case add>>=
3691 tcase_add_test(tc_string_eval, test_setup_teardown);
3692 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
3693 tcase_add_test(tc_string_eval, test_simple_eval);
3694 tcase_add_test(tc_string_eval, test_multiple_evals);
3695 tcase_add_test(tc_string_eval, test_eval_with_state);
3696 suite_add_tcase(s, tc_string_eval);
3700 \section{Accelerating function evaluation}
3702 My first version-0.3 code was running very slowly.
3703 With the options suggested in the help
3704 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
3705 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
3706 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
3707 That's only 15 calls per solution, so the search algorithm seems reasonable.
3708 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
3713 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
3715 #endif /* ACCEL_K_H */
3718 <<accel k module makefile lines>>=
3719 NW_SAWSIM_MODS += accel_k
3720 #CHECK_BINS += check_accel_k
3721 check_accel_k_MODS =
3725 #include <assert.h> /* assert() */
3726 #include <stdlib.h> /* realloc(), free(), NULL */
3727 #include "global.h" /* environment_t */
3728 #include "k_model.h" /* k_func_t */
3729 #include "interp.h" /* interp_* */
3730 #include "accel_k.h"
3732 typedef struct accel_k_struct {
3733 interp_table_t *itable;
3739 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3740 static int num_accels = 0;
3741 static accel_k_t *accels=NULL;
3743 /* Wrap k in the standard f(x) acceleration form */
3744 static double k_wrap(double F, void *params)
3746 accel_k_t *p = (accel_k_t *)params;
3747 return (*p->k)(F, p->env, p->params);
3750 static int k_tol(double FA, double kA, double FB, double kB)
3753 if (FB-FA > 1e-12) {
3754 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3755 return 1; /* unnacceptable */
3757 //printf("acceptable tol\n");
3758 return 0; /* acceptable */
3762 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3765 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3766 assert(accels != NULL);
3767 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3769 accels[i].env = env;
3770 accels[i].params = params;
3777 for (i=0; i<num_accels; i++)
3778 interp_table_free(accels[i].itable);
3780 if (accels) free(accels);
3784 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3787 for (i=0; i<num_accels; i++) {
3788 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3789 /* We've already setup this function.
3790 * WARNING: we're assuming param and env are the same. */
3791 return interp_table_eval(accels[i].itable, accels+i, F);
3794 /* set up a new acceleration environment */
3795 i = add_accel_k(k, env, params);
3796 return interp_table_eval(accels[i].itable, accels+i, F);
3800 \section{Tension models}
3801 \label{sec.tension_models}
3803 TODO: tension model intro.
3804 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.
3806 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3807 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]].
3809 <<tension-model.h>>=
3810 #ifndef TENSION_MODEL_H
3811 #define TENSION_MODEL_H
3813 <<tension handler types>>
3814 <<constant tension model declarations>>
3815 <<hooke tension model declarations>>
3816 <<worm-like chain tension model declarations>>
3817 <<freely-jointed chain tension model declarations>>
3818 <<piston tension model declarations>>
3819 <<find tension definitions>>
3820 <<tension model global declarations>>
3821 #endif /* TENSION_MODEL_H */
3824 <<tension model module makefile lines>>=
3825 NW_SAWSIM_MODS += tension_model
3826 #CHECK_BINS += check_tension_model
3827 check_tension_model_MODS =
3829 <<tension model utils makefile lines>>=
3830 TENSION_MODEL_MODS = tension_model parse list tension_balance
3831 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3832 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3833 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3834 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
3835 notangle -Rtension-model-utils.c $< > $@
3836 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
3837 gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3838 clean_tension_model_utils :
3839 rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3840 @ The pipe symbol [[|]] marks the prerequisits following it (in this
3841 case, the directories) as ‘order-only’ prerequisites. The timestamp
3842 on these prerequisits does not effect whether the rules are executed.
3843 This is appropriate for directories, where we don't need to recompile
3844 after an unrelated has been added to the directory, but only when the
3845 source prerequisites change. See the [[make]] documentation for more
3847 (\url{http://www.gnu.org/software/make/manual/html_node/Prerequisite-Types.html}).
3850 \label{sec.null_tension}
3852 For unstretchable domains.
3854 <<null tension model getopt>>=
3855 {NULL, NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3858 \subsection{Constant}
3859 \label{sec.const_tension}
3861 An infinitely stretchable domain providing a constant tension.
3862 Obviously this cannot be inverted, so you can't put this domain in
3863 series with any other domains. We include it mostly for testing
3866 <<constant tension functions>>=
3867 <<constant tension function>>
3868 <<constant tension parameter create/destroy functions>>
3871 <<constant tension model declarations>>=
3872 <<constant tension function declaration>>
3873 <<constant tension parameter create/destroy function declarations>>
3874 <<constant tension model global declarations>>
3878 <<constant tension function declaration>>=
3879 extern double const_tension_handler(double x, void *pdata);
3881 <<constant tension function>>=
3882 double const_tension_handler(double x, void *pdata)
3884 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3885 list_t *list = data->group_tension_params;
3890 assert(list != NULL); /* empty active group?! */
3891 F = ((const_tension_param_t *)list->d)->F;
3893 fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3895 while (list != NULL) {
3896 assert(((const_tension_param_t *)list->d)->F == F);
3904 <<constant tension structure>>=
3905 typedef struct const_tension_param_struct {
3906 double F; /* tension (force) in N */
3907 } const_tension_param_t;
3911 <<constant tension parameter create/destroy function declarations>>=
3912 extern void *string_create_const_tension_param_t(char **param_strings);
3913 extern void destroy_const_tension_param_t(void *p);
3916 <<constant tension parameter create/destroy functions>>=
3917 const_tension_param_t *create_const_tension_param_t(double F)
3919 const_tension_param_t *ret
3920 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3921 assert(ret != NULL);
3926 void *string_create_const_tension_param_t(char **param_strings)
3928 return create_const_tension_param_t(
3929 safe_strtod(param_strings[0],__FUNCTION__));
3932 void destroy_const_tension_param_t(void *p)
3940 <<constant tension model global declarations>>=
3941 extern int num_const_tension_params;
3942 extern char *const_tension_param_descriptions[];
3943 extern char const_tension_param_string[];
3946 <<constant tension model globals>>=
3947 int num_const_tension_params = 1;
3948 char *const_tension_param_descriptions[] = {"tension F, N"};
3949 char const_tension_param_string[] = "0";
3952 <<constant tension model getopt>>=
3953 {&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}
3959 <<hooke functions>>=
3960 <<internal hooke functions>>
3962 <<inverse hooke handler>>
3963 <<hooke parameter create/destroy functions>>
3966 <<hooke tension model declarations>>=
3967 <<hooke tension function declaration>>
3968 <<hooke tension parameter create/destroy function declarations>>
3969 <<hooke tension model global declarations>>
3973 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
3974 The behavior of a series of springs $k_i$ in series is given by
3976 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
3978 For a simple proof, see Appendix \ref{sec.math_hooke}.
3980 <<hooke tension function declaration>>=
3981 extern double hooke_handler(double x, void *pdata);
3982 extern double inverse_hooke_handler(double F, void *pdata);
3985 First we define a function that computes the effective spring constant
3986 (stored in a single [[hooke_param_t]]) for the entire group.
3987 <<internal hooke functions>>=
3988 static void hooke_param_grouper(tension_handler_data_t *data,
3989 hooke_param_t *param) {
3990 list_t *list = data->group_tension_params;
3994 assert(list != NULL); /* empty active group?! */
3995 while (list != NULL) {
3996 assert( ((hooke_param_t *)list->d)->k > 0 );
3997 k += 1.0/ ((hooke_param_t *)list->d)->k;
3999 fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
4000 __FUNCTION__, ((hooke_param_t *)list->d)->k);
4006 fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N (from data %p)\n",
4007 __FUNCTION__, k, x, k*x, data);
4014 Everyone knows Hooke's law: $F=kx$.
4016 double hooke_handler(double x, void *pdata)
4018 hooke_param_t param = {0};
4021 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
4027 Inverting Hooke's law gives $x=F/k$.
4028 <<inverse hooke handler>>=
4029 double inverse_hooke_handler(double F, void *pdata)
4031 hooke_param_t param = {0};
4034 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
4040 Not much to keep track of for a single effective spring, just the
4041 spring constant $k$.
4042 <<hooke structure>>=
4043 typedef struct hooke_param_struct {
4044 double k; /* spring constant in N/m */
4049 We wrap up our Hooke implementation with some book-keeping functions.
4050 <<hooke tension parameter create/destroy function declarations>>=
4051 extern void *string_create_hooke_param_t(char **param_strings);
4052 extern void destroy_hooke_param_t(void *p);
4056 <<hooke parameter create/destroy functions>>=
4057 hooke_param_t *create_hooke_param_t(double k)
4059 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
4060 assert(ret != NULL);
4065 void *string_create_hooke_param_t(char **param_strings)
4067 return create_hooke_param_t(safe_strtod(param_strings[0], __FUNCTION__));
4070 void destroy_hooke_param_t(void *p)
4077 <<hooke tension model global declarations>>=
4078 extern int num_hooke_params;
4079 extern char *hooke_param_descriptions[];
4080 extern char hooke_param_string[];
4083 <<hooke tension model globals>>=
4084 int num_hooke_params = 1;
4085 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
4086 char hooke_param_string[]="0.05";
4089 <<hooke tension model getopt>>=
4090 {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}
4093 \subsection{Worm-like chain}
4096 We can model several unfolded domains with a single worm-like chain.
4097 <<worm-like chain functions>>=
4098 <<internal worm-like chain functions>>
4099 <<worm-like chain handler>>
4100 <<inverse worm-like chain handler>>
4101 <<worm-like chain create/destroy functions>>
4104 <<worm-like chain tension model declarations>>=
4106 <<worm-like chain handler declaration>>
4107 <<inverse worm-like chain handler declaration>>
4108 <<worm-like chain create/destroy function declarations>>
4109 <<worm-like chain tension model global declarations>>
4113 <<internal worm-like chain functions>>=
4114 <<worm-like chain function>>
4115 <<inverse worm-like chain function>>
4116 <<worm-like chain parameter grouper>>
4119 The combination of all unfolded domains is modeled as a worm like chain,
4120 with the force $F_{\text{WLC}}$ approximately given by
4122 F_{\text{WLC}} \approx \frac{k_B T}{p}
4123 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
4125 where $x$ is the distance between the fixed ends,
4126 $k_B$ is Boltzmann's constant,
4127 $T$ is the absolute temperature,
4128 $p$ is the persistence length, and
4129 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
4130 <<worm-like chain function>>=
4131 static double wlc(double x, double T, double p, double L)
4133 double xL=0.0; /* uses global k_B */
4135 assert(T > 0); assert(p > 0); assert(L > 0);
4136 if (x >= L) return HUGE_VAL;
4137 xL = x/L; /* unitless */
4138 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
4139 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
4140 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
4145 Inverting the approximate WLC is fairly straightforward, if a bit tedious.
4146 Letting $x_L \equiv x/L$ and $F_T \equiv Fp/k_B T$ we have
4148 F_T &\approx \frac{1}{4} \p({\frac{1}{(1-x_L)^2} - 1}) + x_L \\
4149 F_T (1-x_L)^2 &\approx \frac{1}{4} - \frac{1}{4} (1-x_L)^2 + x_L(1-x_L)^2 \\
4150 0 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) (1-x_L)^2 + x_L(1-x_L)^2 \\
4151 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) \p({1-2x_L+x_L^2})
4152 + x_L - 2x_L^2 + x_L^3 \\
4153 &\approx \p[{\frac{1}{4} - \p({F_T + \frac{1}{4}})}]
4154 + x_L \p[{2F_T + \frac{1}{2} + 1}]
4155 + x_L^2 \p[{-\p({F_T + \frac{1}{4}}) - 2}]
4158 + x_L \p({2F_T + \frac{3}{2}})
4159 - x_L^2 \p({F_T + \frac{9}{4}})
4161 &\approx ax_L^3 + bx_L^2 + cx_L + d
4163 where $a \equiv 1$, $b \equiv -F_T - 9/4$, $c \equiv 2F_T + 3/2$, and
4165 % From \citet{wikipedia_cubic_function} the discriminant
4167 % \Delta &= -4b^3d + b^2c^2 - 4ac^3 + 18abcd - 27a^2d^2 \\
4168 % &= -4F_T\p({F_T + \frac{9}{4}})^3
4169 % + \p({F_T + \frac{9}{4}})^2 \cdot \p({2F_T + \frac{3}{2}})^2
4170 % - 4 \p({2F_T + \frac{3}{2}})^3
4171 % + 18F_T \p({F_T + \frac{9}{4}}) \cdot \p({2F_T + \frac{3}{2}})
4173 % &= -4F_T \p({F_T^3 + \frac{27}{4}F_T^2 +\frac{243}{16}F_T +\frac{729}{64}})
4174 %% a^3 + 3a^2b + 3ab^2 + b^3
4175 % + \p({F_T^2 + \frac{9}{2}F_T + \frac{81}{16}})
4176 % \cdot \p({4F_T^2 + 6F_T + \frac{9}{4}}) \\
4177 % &\qquad - 4 \p({8F_T^3 + 18F_T^2 + \frac{27}{2}F_T + \frac{27}{8}})
4178 %% a^3 + 3a^2b + 3ab^2 + b^3
4179 % + 18F_T \p({2F_T^2 + 6F_T + \frac{27}{8}})
4181 % &= -4F_T^4 - 27F_T^3 - \frac{243}{4}F_T^2 - \frac{729}{16}F_T
4182 % + 4F_T^4 + 24F_T^3 +\frac{99}{2}F_T^2 +\frac{81}{2}F_T +\frac{729}{64}\\
4183 % &\qquad - 32F_T^3 - 72F_T^2 - 54F_T - \frac{27}{2}
4184 % + 36F_T^3 + 108F_T^2 + \frac{243}{4}F_T
4186 % &= F_T^4 \p({-4 + 4}) + F_T^3 \p({-27 + 24 - 32 + 36})
4187 % + F_T^2 \p({-\frac{243}{4} + \frac{99}{2} - 72 + 108 -27}) \\
4188 % &\qquad + F_T \p({-\frac{729}{16} + \frac{81}{2} - 54 + \frac{243}{4}})
4189 % \p({\frac{729}{64} - \frac{27}{2}}) \\
4190 % &= F_T^3 - \frac{9}{4}F_T^2 + \frac{27}{16}F_T - \frac{135}{64}
4192 We can plug these coefficients into the GSL cubic solver to invert the
4193 WLC\citep{galassi05}. The applicable root is always the one which
4194 comes before the singularity, which will be the smallest real root.
4195 <<inverse worm-like chain function>>=
4196 static double inverse_wlc(double F, double T, double p, double L)
4198 double FT = F*p/(k_B*T); /* uses global k_B */
4199 double xL0, xL1, xL2;
4202 assert(T > 0); assert(p > 0); assert(L > 0);
4205 num_roots = gsl_poly_solve_cubic(-(FT+2.25),2*FT+1.5,-FT, &xL0, &xL1, &xL2);
4206 assert(num_roots > 0); assert(xL0 >= -DOUBLE_PRECISION); assert(xL0 < 1);
4207 if (xL0 < 0) xL0 = 0.0;
4213 First we define a function that computes the effective WLC parameters
4214 (stored in a single [[wlc_param_t]]) for the entire group.
4215 <<worm-like chain parameter grouper>>=
4216 static void wlc_param_grouper(tension_handler_data_t *data,
4217 wlc_param_t *param) {
4218 list_t *list = data->group_tension_params;
4219 double p=0.0, L=0.0;
4222 assert(list != NULL); /* empty active group?! */
4223 p = ((wlc_param_t *)list->d)->p;
4224 while (list != NULL) {
4225 /* enforce consistent p */
4226 assert( ((wlc_param_t *)list->d)->p == p);
4227 L += ((wlc_param_t *)list->d)->L;
4229 fprintf(stderr, "%s: WLC section %g m long in series\n",
4230 __FUNCTION__, ((wlc_param_t *)list->d)->L);
4235 fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
4236 __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
4244 <<worm-like chain handler declaration>>=
4245 extern double wlc_handler(double x, void *pdata);
4248 This model requires that each unfolded domain in the group have the
4249 same persistence length.
4250 <<worm-like chain handler>>=
4251 double wlc_handler(double x, void *pdata)
4253 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4254 wlc_param_t param = {0};
4257 wlc_param_grouper(data, ¶m);
4258 return wlc(x, data->env->T, param.p, param.L);
4263 <<inverse worm-like chain handler declaration>>=
4264 extern double inverse_wlc_handler(double F, void *pdata);
4267 <<inverse worm-like chain handler>>=
4268 double inverse_wlc_handler(double F, void *pdata)
4270 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4271 wlc_param_t param = {0};
4273 wlc_param_grouper(data, ¶m);
4274 return inverse_wlc(F, data->env->T, param.p, param.L);
4279 <<worm-like chain structure>>=
4280 typedef struct wlc_param_struct {
4281 double p; /* WLC persistence length (m) */
4282 double L; /* WLC contour length (m) */
4286 <<worm-like chain create/destroy function declarations>>=
4287 extern void *string_create_wlc_param_t(char **param_strings);
4288 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
4291 <<worm-like chain create/destroy functions>>=
4292 wlc_param_t *create_wlc_param_t(double p, double L)
4294 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
4295 assert(ret != NULL);
4301 void *string_create_wlc_param_t(char **param_strings)
4303 return create_wlc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4304 safe_strtod(param_strings[1], __FUNCTION__));
4307 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
4315 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.
4316 TODO (now we copy the string before parsing, but still do this for clarity).
4318 <<valid string write code>>=
4319 char string[] = "abc";
4322 <<valid string write code>>=
4323 char *string = "abc";
4327 <<worm-like chain tension model global declarations>>=
4328 extern int num_wlc_params;
4329 extern char *wlc_param_descriptions[];
4330 extern char wlc_param_string[];
4332 <<worm-like chain tension model globals>>=
4333 int num_wlc_params = 2;
4334 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
4335 char wlc_param_string[]="0.39e-9,28.4e-9";
4338 <<worm-like chain tension model getopt>>=
4339 {&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}
4341 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
4343 \subsection{Freely-jointed chain}
4346 <<freely-jointed chain functions>>=
4347 <<freely-jointed chain function>>
4348 <<freely-jointed chain handler>>
4349 <<freely-jointed chain create/destroy functions>>
4352 <<freely-jointed chain tension model declarations>>=
4353 <<freely-jointed chain handler declaration>>
4354 <<freely-jointed chain create/destroy function declarations>>
4355 <<freely-jointed chain tension model global declarations>>
4359 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
4360 The tension of a chain of $N$ such links is given by
4362 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
4364 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}.
4365 <<freely-jointed chain function>>=
4366 double langevin(double x, void *params)
4369 return 1.0/tanh(x) - 1/x;
4371 <<one dimensional bracket declaration>>
4372 <<one dimensional solve declaration>>
4373 double inv_langevin(double x)
4375 double lb=0.0, ub=1.0, ret;
4376 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
4377 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
4381 double fjc(double x, double T, double l, int N)
4383 double L = l*(double)N;
4384 if (x == 0) return 0;
4385 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
4386 return k_B*T/l * inv_langevin(x/L);
4390 <<freely-jointed chain handler declaration>>=
4391 extern double fjc_handler(double x, void *pdata);
4393 <<freely-jointed chain handler>>=
4394 double fjc_handler(double x, void *pdata)
4396 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4397 list_t *list = data->group_tension_params;
4402 assert(list != NULL); /* empty active group?! */
4403 l = ((fjc_param_t *)list->d)->l;
4404 while (list != NULL) {
4405 /* enforce consistent l */
4406 assert( ((fjc_param_t *)list->d)->l == l);
4407 N += ((fjc_param_t *)list->d)->N;
4409 fprintf(stderr, "%s: FJC section %d links long in series\n",
4410 __FUNCTION__, ((fjc_param_t *)list->d)->N);
4415 fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
4416 __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
4418 return fjc(x, data->env->T, l, N);
4422 <<freely-jointed chain structure>>=
4423 typedef struct fjc_param_struct {
4424 double l; /* FJC link length (m) */
4425 int N; /* FJC number of links */
4429 <<freely-jointed chain create/destroy function declarations>>=
4430 extern void *string_create_fjc_param_t(char **param_strings);
4431 extern void destroy_fjc_param_t(void *p);
4434 <<freely-jointed chain create/destroy functions>>=
4435 fjc_param_t *create_fjc_param_t(double l, double N)
4437 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
4438 assert(ret != NULL);
4444 void *string_create_fjc_param_t(char **param_strings)
4446 return create_fjc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4447 safe_strtod(param_strings[1], __FUNCTION__));
4450 void destroy_fjc_param_t(void *p)
4457 <<freely-jointed chain tension model global declarations>>=
4458 extern int num_fjc_params;
4459 extern char *fjc_param_descriptions[];
4460 extern char fjc_param_string[];
4463 <<freely-jointed chain tension model globals>>=
4464 int num_fjc_params = 2;
4465 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
4466 char fjc_param_string[]="0.5e-9,200";
4469 <<freely-jointed chain tension model getopt>>=
4470 {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}
4474 \label{sec.piston_tension}
4476 An infinitely stretchable domain with no tension for extensions less
4477 than a particular threshold $L$, and infinite tension for $x>L$. The
4478 tension at $x=L$ is undefined, but may be any positive value. The
4479 model is called the ``piston'' model because it resembles a
4480 frictionless piston sliding in a rigid cylinder of length $L$.
4483 0 & \text{if $x<L$}, \\
4484 \infty & \text{if $x>L$}.
4488 Note that because of it's infinte stiffness at $x=L$, fully extended
4489 domains with this tension model will be identical to completely rigid
4490 domains. However, there is a distinction between this model and rigid
4491 domains for $x<L$, because any reactions that occur at $F=0$
4492 (e.g. refolding) will have more time to occur.
4494 Because the tension at $x=L$ is undefined, the user must make sure
4495 domains with this tension model are never the initial active group,
4496 because it would break [[tension_balance()]] in [[find_tension()]]
4497 (see Section \ref{sec.tension_balance}).
4499 <<piston tension functions>>=
4500 <<piston tension parameter grouper>>
4501 <<piston tension handler>>
4502 <<inverse piston tension handler>>
4503 <<piston tension parameter create/destroy functions>>
4506 <<piston tension model declarations>>=
4507 <<piston tension handler declaration>>
4508 <<inverse piston tension handler declaration>>
4509 <<piston tension parameter create/destroy function declarations>>
4510 <<piston tension model global declarations>>
4514 First we define a function that computes the effective piston parameters
4515 (stored in a single [[piston_tension_param_t]]) for the entire group.
4516 <<piston tension parameter grouper>>=
4517 static void piston_tension_param_grouper(tension_handler_data_t *data,
4518 piston_tension_param_t *param) {
4519 list_t *list = data->group_tension_params;
4523 assert(list != NULL); /* empty active group?! */
4524 while (list != NULL) {
4525 L += ((piston_tension_param_t *)list->d)->L;
4531 <<piston tension handler declaration>>=
4532 extern double piston_tension_handler(double x, void *pdata);
4534 <<piston tension handler>>=
4535 double piston_tension_handler(double x, void *pdata)
4537 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4538 piston_tension_param_t param = {0};
4541 piston_tension_param_grouper(data, ¶m);
4542 if (x <= param.L) F = 0;
4545 fprintf(stderr, "%s: x %g, L %g, tension %g N\n", __FUNCTION__, x, param.L, F);
4552 We abuse [[x_of_xo()]] and [[oneD_solve()]] with this inverse
4553 definition (see Section \ref{sec.tension_balance}), because the
4554 $x(F=0)<=L$, but our function returns $x(F=0)=0$. This is fine when
4555 the total extension \emph{is} zero, but cheats a bit when there is
4556 some total extension. For any non-zero extension, the initial active
4557 group produces some tension (we hope. True for all our other tension
4558 models). This tension fully extends the piston group (of which there
4559 should be only one, since all piston states can get lumped into the
4560 same tension group.). If the total extension is $\le$ the target
4561 extension, the full extension is accurate, so the inverse was valid
4562 after all. If not, [[oneD_solve()]] will halve the extension of the
4563 first group, reducing the overall tension. For total extension less
4564 than the piston extension, this bisection continues for [[max_steps]],
4565 leaving $x_0$ reduced by a factor of $2^\text{max\_steps}$, a very
4566 small number. Because of this, the returned tension $F_0(x_0)$ is
4567 very small, even though the total extension found by [[x_of_xo()]]
4568 is still $>$ the piston length.
4570 While this works, it is clearly not the most elegant or robust
4571 solution. TODO: think of (and implement) a better procedure.
4573 <<inverse piston tension handler declaration>>=
4574 extern double inverse_piston_tension_handler(double x, void *pdata);
4576 <<inverse piston tension handler>>=
4577 double inverse_piston_tension_handler(double F, void *pdata)
4579 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4580 piston_tension_param_t param = {0};
4583 piston_tension_param_grouper(data, ¶m);
4584 if (F == 0.0) return 0;
4590 <<piston tension structure>>=
4591 typedef struct piston_tension_param_struct {
4592 double L; /* length of piston in m */
4593 } piston_tension_param_t;
4597 <<piston tension parameter create/destroy function declarations>>=
4598 extern void *string_create_piston_tension_param_t(char **param_strings);
4599 extern void destroy_piston_tension_param_t(void *p);
4603 <<piston tension parameter create/destroy functions>>=
4604 piston_tension_param_t *create_piston_tension_param_t(double L)
4606 piston_tension_param_t *ret
4607 = (piston_tension_param_t *) malloc(sizeof(piston_tension_param_t));
4608 assert(ret != NULL);
4613 void *string_create_piston_tension_param_t(char **param_strings)
4615 return create_piston_tension_param_t(
4616 safe_strtod(param_strings[0],__FUNCTION__));
4619 void destroy_piston_tension_param_t(void *p)
4627 <<piston tension model global declarations>>=
4628 extern int num_piston_tension_params;
4629 extern char *piston_tension_param_descriptions[];
4630 extern char piston_tension_param_string[];
4634 <<piston tension model globals>>=
4635 int num_piston_tension_params = 1;
4636 char *piston_tension_param_descriptions[] = {"length L, m"};
4637 char piston_tension_param_string[] = "0";
4641 <<piston tension model getopt>>=
4642 {&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}
4645 \subsection{Tension model implementation}
4647 <<tension-model.c>>=
4650 <<tension model includes>>
4651 #include "tension_model.h"
4652 <<tension model internal definitions>>
4653 <<tension model globals>>
4654 <<tension model functions>>
4657 <<tension model includes>>=
4658 #include <assert.h> /* assert() */
4659 #include <stdlib.h> /* NULL, strto*() */
4660 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4661 #include <stdio.h> /* fprintf(), stdout */
4662 #include <errno.h> /* errno, ERANGE */
4665 #include "tension_balance.h" /* oneD_*() */
4668 <<tension model internal definitions>>=
4669 <<constant tension structure>>
4671 <<worm-like chain structure>>
4672 <<freely-jointed chain structure>>
4673 <<piston tension structure>>
4676 <<tension model globals>>=
4677 <<tension handler array global>>
4678 <<constant tension model globals>>
4679 <<hooke tension model globals>>
4680 <<worm-like chain tension model globals>>
4681 <<freely-jointed chain tension model globals>>
4682 <<piston tension model globals>>
4685 <<tension model functions>>=
4687 <<constant tension functions>>
4689 <<worm-like chain functions>>
4690 <<freely-jointed chain functions>>
4691 <<piston tension functions>>
4694 \subsection{Utilities}
4696 The tension models can get complicated, and users may want to reassure
4697 themselves that this portion of the input physics are functioning properly. The
4698 stand-alone program [[tension_model_utils]] demonstrates the tension model
4699 interface without the overhead of having to understand a full simulation.
4700 [[tension_model_utils]] takes command line model arguments like the full
4701 simulation, and prints $F(x)$ data to the screen for the selected model over a
4704 <<tension-model-utils.c>>=
4706 <<tension model utility includes>>
4707 <<tension model utility definitions>>
4708 <<tension model utility getopt functions>>
4710 <<tension model utility main>>
4713 <<tension model utility main>>=
4714 int main(int argc, char **argv)
4716 <<tension model getopt array definition>>
4717 tension_model_getopt_t *model = NULL;
4721 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
4722 tension_handler_data_t tdata;
4723 int num_param_args; /* for INIT_MODEL() */
4724 char **param_args; /* for INIT_MODEL() */
4726 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model,
4727 &Fmax, &Xmax, &flags);
4729 if (flags & VFLAG) {
4730 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4732 INIT_MODEL("utility", model, model->params, params);
4733 tdata.group_tension_params = NULL;
4734 push(&tdata.group_tension_params, params, 1);
4736 tdata.persist = NULL;
4737 if (model->handler == NULL) {
4738 printf("No tension function for model %s\n", model->name);
4744 printf("#Distance (m)\tForce (N)\n");
4745 for (i=0; i<=N; i++) {
4746 x = Xmax*i/(double)N;
4747 F = (*model->handler)(x, &tdata);
4748 if (F < 0 || F > Fmax) break;
4749 printf("%g\t%g\n", x, F);
4751 if (flags & VFLAG && i <= N) {
4752 /* explain exit condition */
4754 printf("Impossible force %g\n", F);
4756 printf("Reached large force limit %g > %g\n", F, Fmax);
4759 params = pop(&tdata.group_tension_params);
4760 if (model->destructor)
4761 (*model->destructor)(params);
4766 <<tension model utility includes>>=
4769 #include <string.h> /* strlen() */
4770 #include <assert.h> /* assert() */
4771 #include <errno.h> /* errno, ERANGE */
4775 #include "tension_model.h"
4778 <<tension model utility definitions>>=
4779 <<version definition>>
4780 #define VFLAG 1 /* verbose */
4781 <<string comparison definition and globals>>
4782 <<tension model getopt definitions>>
4783 <<initialize model definition>>
4787 <<tension model utility getopt functions>>=
4789 <<index tension model>>
4790 <<tension model utility help>>
4791 <<tension model utility get options>>
4794 <<tension model utility help>>=
4795 void help(char *prog_name,
4797 int n_tension_models, tension_model_getopt_t *tension_models,
4798 int tension_model, double Fmax, double Xmax)
4801 printf("usage: %s [options]\n", prog_name);
4802 printf("Version %s\n\n", VERSION);
4803 printf("A utility for understanding the available tension models\n\n");
4804 printf("Environmental options:\n");
4805 printf("-T T\tTemperature T (currently %g K)\n", env->T);
4806 printf("-C T\tYou can also set the temperature T in Celsius\n");
4807 printf("Model options:\n");
4808 printf("-m model\tFolded tension model (currently %s)\n",
4809 tension_models[tension_model].name);
4810 printf("-a args\tFolded tension model argument string (currently %s)\n",
4811 tension_models[tension_model].params);
4812 printf("F(x) is calculated for a range of x and printed\n");
4813 printf("For example:\n");
4814 printf(" #Distance (m)\tForce (N)\n");
4815 printf(" 123.456\t7.89\n");
4817 printf("-F\tSet the maximum F value for the standard mode F(x) output (currently %g)\n", Fmax);
4818 printf("-X\tSet the maximum x value for the standard mode F(x) output (currently %g)\n", Xmax);
4819 printf("-V\tChange output to verbose mode\n");
4820 printf("-h\tPrint this help and exit\n");
4822 printf("Tension models:\n");
4823 for (i=0; i<n_tension_models; i++) {
4824 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
4825 for (j=0; j < tension_models[i].num_params ; j++)
4826 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
4827 printf("\t\tdefaults: %s\n", tension_models[i].params);
4832 <<tension model utility get options>>=
4833 void get_options(int argc, char **argv, environment_t *env,
4834 int n_tension_models, tension_model_getopt_t *tension_models,
4835 tension_model_getopt_t **model, double *Fmax, double *Xmax,
4836 unsigned int *flags)
4838 char *prog_name = NULL;
4839 char c, options[] = "T:C:m:a:F:X:Vh";
4840 int tension_model=0, num_strings;
4841 extern char *optarg;
4842 extern int optind, optopt, opterr;
4846 /* setup defaults */
4848 prog_name = argv[0];
4849 env->T = 300.0; /* K */
4850 *Fmax = 1e5; /* N */
4851 *Xmax = 1e-6; /* m */
4853 *model = tension_models;
4855 while ((c=getopt(argc, argv, options)) != -1) {
4857 case 'T': env->T = safe_strtod(optarg, "-T"); break;
4858 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
4860 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
4861 *model = tension_models+tension_model;
4864 tension_models[tension_model].params = optarg;
4866 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
4867 case 'X': *Xmax = safe_strtod(optarg, "-X"); break;
4868 case 'V': *flags |= VFLAG; break;
4870 fprintf(stderr, "unrecognized option '%c'\n", optopt);
4871 /* fall through to default case */
4873 help(prog_name, env, n_tension_models, tension_models, tension_model, *Fmax, *Xmax);
4882 \section{Unfolding rate models}
4883 \label{sec.k_models}
4885 We have two main choices for determining $k$: Bell's model, which assumes the
4886 folded domains exist in a single `folded' state and make instantaneous,
4887 stochastic transitions over a free energy barrier; and Kramers' model, which
4888 treats the unfolding as a thermalized diffusion process.
4889 We deal with these options like we dealt with the different tension models: by bundling all model-specific
4890 parameters into structures, with model specific functions for determining $k$.
4892 <<k func definition>>=
4893 typedef double k_func_t(double F, environment_t *env, void *params);
4896 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
4897 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]].
4899 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
4900 because \LaTeX\ doesn't like underscores outside math mode.
4901 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
4902 You could use spaces instead (`\verb+ +'),
4903 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
4904 which I don't like as much.
4910 <<k func definition>>
4911 <<null k declarations>>
4912 <<const k declarations>>
4913 <<bell k declarations>>
4914 <<kbell k declarations>>
4915 <<kramers k declarations>>
4916 <<kramers integ k declarations>>
4917 #endif /* K_MODEL_H */
4920 <<k model module makefile lines>>=
4921 NW_SAWSIM_MODS += k_model
4922 CHECK_BINS += check_k_model
4923 check_k_model_MODS = parse string_eval
4925 [[check_k_model]] also depends on the parse module.
4927 <<k model utils makefile lines>>=
4928 K_MODEL_MODS = k_model parse string_eval
4929 K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
4930 $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h)
4931 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
4932 notangle -Rk-model-utils.c $< > $@
4933 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
4934 gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
4935 clean_k_model_utils :
4936 rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
4942 For domains that are never folded.
4944 <<null k declarations>>=
4946 <<null k functions>>=
4951 <<null k model getopt>>=
4952 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
4955 \subsection{Constant rate model}
4958 This is a pretty boring model, but it is usefull for testing the engine.
4962 where $k_0$ is the constant unfolding rate.
4964 <<const k functions>>=
4965 <<const k function>>
4966 <<const k structure create/destroy functions>>
4969 <<const k declarations>>=
4970 <<const k function declaration>>
4971 <<const k structure create/destroy function declarations>>
4972 <<const k global declarations>>
4976 <<const k structure>>=
4977 typedef struct const_k_param_struct {
4978 double knot; /* transition rate k_0 (frac population per s) */
4982 <<const k function declaration>>=
4983 double const_k(double F, environment_t *env, void *const_k_params);
4985 <<const k function>>=
4986 double const_k(double F, environment_t *env, void *const_k_params)
4987 { /* Returns the rate constant k in frac pop per s. */
4988 const_k_param_t *p = (const_k_param_t *) const_k_params;
4990 assert(p->knot > 0);
4995 <<const k structure create/destroy function declarations>>=
4996 void *string_create_const_k_param_t(char **param_strings);
4997 void destroy_const_k_param_t(void *p);
5000 <<const k structure create/destroy functions>>=
5001 const_k_param_t *create_const_k_param_t(double knot)
5003 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
5004 assert(ret != NULL);
5009 void *string_create_const_k_param_t(char **param_strings)
5011 return create_const_k_param_t(safe_strtod(param_strings[0], __FUNCTION__));
5014 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
5021 <<const k global declarations>>=
5022 extern int num_const_k_params;
5023 extern char *const_k_param_descriptions[];
5024 extern char const_k_param_string[];
5027 <<const k globals>>=
5028 int num_const_k_params = 1;
5029 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
5030 char const_k_param_string[]="1";
5033 <<const k model getopt>>=
5034 {"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}
5037 \subsection{Bell's model}
5040 Quantitatively, Bell's model gives $k$ as
5042 k = k_0 \cdot e^{F dx / k_B T} \;,
5044 where $k_0$ is the unforced unfolding rate,
5045 $F$ is the applied unfolding force,
5046 $dx$ is the distance to the transition state, and
5047 $k_B T$ is the thermal energy\citep{schlierf06}.
5049 <<bell k functions>>=
5051 <<bell k structure create/destroy functions>>
5054 <<bell k declarations>>=
5055 <<bell k function declaration>>
5056 <<bell k structure create/destroy function declarations>>
5057 <<bell k global declarations>>
5061 <<bell k structure>>=
5062 typedef struct bell_param_struct {
5063 double knot; /* transition rate k_0 (frac population per s) */
5064 double dx; /* `distance to transition state' (nm) */
5068 <<bell k function declaration>>=
5069 double bell_k(double F, environment_t *env, void *bell_params);
5071 <<bell k function>>=
5072 double bell_k(double F, environment_t *env, void *bell_params)
5073 { /* Returns the rate constant k in frac pop per s.
5074 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5075 * uses global k_B in J/K */
5076 bell_param_t *p = (bell_param_t *) bell_params;
5077 assert(F >= 0); assert(env->T > 0);
5079 assert(p->knot > 0); //assert(p->dx >= 0); /* p < 0 for refolding reactions, etc. */
5081 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
5082 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
5083 p->knot * exp(F*p->dx / (k_B*env->T) ));
5084 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
5086 return p->knot * exp(F*p->dx / (k_B*env->T) );
5090 <<bell k structure create/destroy function declarations>>=
5091 void *string_create_bell_param_t(char **param_strings);
5092 void destroy_bell_param_t(void *p);
5095 <<bell k structure create/destroy functions>>=
5096 bell_param_t *create_bell_param_t(double knot, double dx)
5098 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
5099 assert(ret != NULL);
5105 void *string_create_bell_param_t(char **param_strings)
5107 return create_bell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5108 safe_strtod(param_strings[1], __FUNCTION__));
5111 void destroy_bell_param_t(void *p /* bell_param_t * */)
5118 <<bell k global declarations>>=
5119 extern int num_bell_params;
5120 extern char *bell_param_descriptions[];
5121 extern char bell_param_string[];
5125 int num_bell_params = 2;
5126 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
5127 char bell_param_string[]="3.3e-4,0.25e-9";
5130 <<bell k model getopt>>=
5131 {"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}
5133 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
5134 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5137 \subsection{Linker stiffness corrected Bell model}
5140 Walton et. al showed that the Bell model constant force approximation
5141 doesn't hold for biotin-streptavadin binding, and I extended their
5142 results to I27 for the 2009 Biophysical Society Annual
5143 Meeting\cite{walton08,king09}. More details to follow when I get done
5144 with the conference\ldots
5146 We adjust Bell's model to
5148 k = k_0 \cdot e^{\frac{F dx - \frac{1}{2} \kappa dx^2}{k_B T}} \;,
5150 where $k_0$ is the unforced unfolding rate,
5151 $F$ is the applied unfolding force,
5152 $\kappa$ is the effective spring constant,
5153 $dx$ is the distance to the transition state, and
5154 $k_B T$ is the thermal energy\citep{schlierf06}.
5156 <<kbell k functions>>=
5157 <<kbell k function>>
5158 <<kbell k structure create/destroy functions>>
5161 <<kbell k declarations>>=
5162 <<kbell k function declaration>>
5163 <<kbell k structure create/destroy function declarations>>
5164 <<kbell k global declarations>>
5168 <<kbell k structure>>=
5169 typedef struct kbell_param_struct {
5170 double knot; /* transition rate k_0 (frac population per s) */
5171 double dx; /* `distance to transition state' (nm) */
5175 <<kbell k function declaration>>=
5176 double kbell_k(double F, environment_t *env, void *kbell_params);
5178 <<kbell k function>>=
5179 double kbell_k(double F, environment_t *env, void *kbell_params)
5180 { /* Returns the rate constant k in frac pop per s.
5181 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5182 * uses global k_B in J/K */
5183 kbell_param_t *p = (kbell_param_t *) kbell_params;
5184 assert(F >= 0); assert(env->T > 0); assert(env->stiffness >= 0);
5186 assert(p->knot > 0); assert(p->dx >= 0);
5187 return p->knot * exp((F*p->dx - 0.5*env->stiffness*(p->dx)*(p->dx))/(k_B*env->T) );
5191 <<kbell k structure create/destroy function declarations>>=
5192 void *string_create_kbell_param_t(char **param_strings);
5193 void destroy_kbell_param_t(void *p);
5196 <<kbell k structure create/destroy functions>>=
5197 kbell_param_t *create_kbell_param_t(double knot, double dx)
5199 kbell_param_t *ret = (kbell_param_t *) malloc(sizeof(kbell_param_t));
5200 assert(ret != NULL);
5206 void *string_create_kbell_param_t(char **param_strings)
5208 return create_kbell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5209 safe_strtod(param_strings[1], __FUNCTION__));
5212 void destroy_kbell_param_t(void *p /* kbell_param_t * */)
5219 <<kbell k global declarations>>=
5220 extern int num_kbell_params;
5221 extern char *kbell_param_descriptions[];
5222 extern char kbell_param_string[];
5225 <<kbell k globals>>=
5226 int num_kbell_params = 2;
5227 char *kbell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
5228 char kbell_param_string[]="3.3e-4,0.25e-9";
5231 <<kbell k model getopt>>=
5232 {"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}
5234 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
5235 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5238 \subsection{Kramer's model}
5241 Kramer's model gives $k$ as
5243 % \frac{1}{k} = \frac{1}{D}
5244 % \int_{x_\text{min}}^{x_\text{max}}
5245 % e^{\frac{-U_F(x)}{k_B T}}
5246 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5249 %where $D$ is the diffusion constant,
5250 %$U_F(x)$ is the free energy along the unfolding coordinate, and
5251 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
5252 % before the minimum and after the maximum, respectively \citep{schlierf06}.
5254 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
5256 where $D$ is the diffusion constant,
5258 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
5260 are length scales where
5261 $x_c(F)$ is the location of the bound state and
5262 $x_{ts}(F)$ is the location of the transition state,
5263 $E(F,x)$ is the free energy, and
5264 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
5265 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
5266 \citet{evans97} Eqn.~3).
5268 <<kramers k functions>>=
5269 <<kramers k function>>
5270 <<kramers k structure create/destroy functions>>
5273 <<kramers k declarations>>=
5274 <<kramers k function declaration>>
5275 <<kramers k structure create/destroy function declarations>>
5276 <<kramers k global declarations>>
5280 <<kramers k structure>>=
5281 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
5282 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
5284 typedef struct kramers_param_struct {
5285 double D; /* diffusion rate (um/s) */
5292 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
5293 //kramers_x_func_t *xts; /* function returning transition x (nm) */
5294 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
5295 //kramers_E_func_t *E; /* function returning E (J) */
5296 //double *E_params; /* parametrize the energy landscape */
5297 //int n_E_params; /* length of E_params array */
5301 <<kramers k function declaration>>=
5302 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
5303 extern double kramers_k(double F, environment_t *env, void *kramers_params);
5305 <<kramers k function>>=
5306 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
5308 static char input[160]={0};
5309 static char output[80]={0};
5310 /* setup the environment */
5311 fprintf(in, "F = %g; T = %g\n", F, T);
5312 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
5313 string_eval(in, out, input, 80, output);
5314 return safe_strtod(output, __FUNCTION__);
5317 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
5319 static char input[160]={0};
5320 static char output[80]={0};
5321 /* setup the environment */
5322 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
5323 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
5324 string_eval(in, out, input, 80, output);
5325 return safe_strtod(output, __FUNCTION__);
5328 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
5330 kramers_param_t *p = (kramers_param_t *) kramers_params;
5331 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
5334 double kramers_k(double F, environment_t *env, void *kramers_params)
5335 { /* Returns the rate constant k in frac pop per s.
5336 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5337 * uses global k_B in J/K */
5338 kramers_param_t *p = (kramers_param_t *) kramers_params;
5339 double kbT, xc, xts, lc, lts, Eb;
5340 assert(F >= 0); assert(env->T > 0);
5343 assert(p->xc != NULL); assert(p->xts != NULL);
5344 assert(p->ddEdxx != NULL); assert(p->E != NULL);
5345 assert(p->in != NULL); assert(p->out != NULL);
5347 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
5348 if (xc == -1.0) return -1.0; /* error (Too much force) */
5349 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
5350 if (xc == -1.0) return -1.0; /* error (Too much force) */
5351 lc = sqrt(2.0*M_PI*kbT /
5352 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
5353 lts = sqrt(-2.0*M_PI*kbT/
5354 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
5355 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
5356 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
5357 //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));
5358 return p->D/(lc*lts) * exp(-Eb/kbT);
5362 <<kramers k structure create/destroy function declarations>>=
5363 void *string_create_kramers_param_t(char **param_strings);
5364 void destroy_kramers_param_t(void *p);
5367 <<kramers k structure create/destroy functions>>=
5368 kramers_param_t *create_kramers_param_t(double D,
5369 char *xc_fn, char *xts_fn,
5370 char *E_fn, char *ddEdxx_fn,
5371 char *global_define_string)
5372 // kramers_x_func_t *xc, kramers_x_func_t *xts,
5373 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
5374 // double *E_params, int n_E_params)
5376 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
5377 assert(ret != NULL);
5378 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5379 assert(ret->xc != NULL);
5380 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5381 assert(ret->xts != NULL);
5382 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
5383 assert(ret->E != NULL);
5384 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
5385 assert(ret->ddEdxx != NULL);
5387 strcpy(ret->xc, xc_fn);
5388 strcpy(ret->xts, xts_fn);
5389 strcpy(ret->E, E_fn);
5390 strcpy(ret->ddEdxx, ddEdxx_fn);
5391 //ret->E_params = E_params;
5392 //ret->n_E_params = n_E_params;
5393 string_eval_setup(&ret->in, &ret->out);
5394 fprintf(ret->in, "kB = %g\n", k_B);
5395 fprintf(ret->in, "%s\n", global_define_string);
5399 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
5400 void *string_create_kramers_param_t(char **param_strings)
5402 return create_kramers_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5410 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
5412 kramers_param_t *pk = (kramers_param_t *)p;
5414 string_eval_teardown(&pk->in, &pk->out);
5416 // free(pk->E_params);
5422 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
5423 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.
5424 We follow \citet{shillcock98} and use
5426 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
5427 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
5430 where TODO, variable meanings.%\citep{shillcock98}.
5431 The first and second derivatives of $E(F,E_b,x,x_s)$ are
5433 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\\
5434 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
5436 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
5437 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
5438 The roots are therefor at
5440 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
5441 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
5442 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
5445 <<kramers k global declarations>>=
5446 extern int num_kramers_params;
5447 extern char *kramers_param_descriptions[];
5448 extern char kramers_param_string[];
5451 <<kramers k globals>>=
5452 int num_kramers_params = 6;
5453 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"};
5454 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)}";
5456 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
5457 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
5458 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
5459 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.
5460 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
5461 It works so long as [[val_1]] is not [[false]].
5463 <<kramers k model getopt>>=
5464 {"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}
5467 \subsection{Kramer's model (natural cubic splines)}
5468 \label{sec.kramers_integ}
5470 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
5471 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
5472 \citet{schlierf06} Eqn.~2)
5474 \frac{1}{k} = \frac{1}{D}
5475 \int_{x_\text{min}}^{x_\text{max}}
5476 e^{\frac{U_F(x)}{k_B T}}
5477 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5480 where $D$ is the diffusion constant,
5481 $U_F(x)$ is the free energy along the unfolding coordinate, and
5482 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
5483 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
5485 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
5486 interpolating between them with cubic splines.
5487 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
5489 <<kramers integ k functions>>=
5490 <<spline functions>>
5491 <<kramers integ A k functions>>
5492 <<kramers integ B k functions>>
5493 <<kramers integ k function>>
5494 <<kramers integ k structure create/destroy functions>>
5497 <<kramers integ k declarations>>=
5498 <<kramers integ k function declaration>>
5499 <<kramers integ k structure create/destroy function declarations>>
5500 <<kramers integ k global declarations>>
5504 <<kramers integ k structure>>=
5505 typedef double func_t(double x, void *params);
5506 typedef struct kramers_integ_param_struct {
5507 double D; /* diffusion rate (um/s) */
5508 double x_min; /* integration bounds */
5510 func_t *E; /* function returning E (J) */
5511 void *E_params; /* parametrize the energy landscape */
5512 destroy_data_func_t *destroy_E_params;
5514 double F; /* for passing into GSL evaluations */
5516 } kramers_integ_param_t;
5519 <<spline param structure>>=
5520 typedef struct spline_param_struct {
5521 int n_knots; /* natural cubic spline knots for U(x) */
5524 gsl_spline *spline; /* GSL spline parameters */
5525 gsl_interp_accel *acc; /* GSL spline acceleration data */
5529 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
5531 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5533 <<kramers integ A k functions>>=
5534 double kramers_integ_k_integrandA(double x, void *params)
5536 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5537 double U = (*p->E)(x, p->E_params);
5538 double Fx = p->F * x;
5539 double kbT = k_B * p->env->T;
5540 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
5541 return exp(-(U-Fx)/kbT);
5543 double kramers_integ_k_integralA(double x, void *params)
5545 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5547 double result, abserr;
5548 size_t neval; /* for qng */
5549 static gsl_integration_workspace *w = NULL;
5550 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5551 f.function = &kramers_integ_k_integrandA;
5553 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5554 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5555 //fprintf(stderr, "integralA = %g\n", result);
5561 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
5562 \int_{x_\text{min}}^{x_\text{max}}
5563 e^{\frac{U_F(x)}{k_B T}}
5564 \text{Integral}_A(x)
5567 <<kramers integ B k functions>>=
5568 double kramers_integ_k_integrandB(double x, void *params)
5570 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5571 double U = (*p->E)(x, p->E_params);
5572 double Fx = p->F * x;
5573 double kbT = k_B * p->env->T;
5574 double intA = kramers_integ_k_integralA(x, params);
5575 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
5576 return exp((U-Fx)/kbT)*intA;
5578 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
5580 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5582 double result, abserr;
5583 size_t neval; /* for qng */
5584 static gsl_integration_workspace *w = NULL;
5585 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5586 f.function = &kramers_integ_k_integrandB;
5590 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5591 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5592 //fprintf(stderr, "integralB = %g\n", result);
5597 With the integrals taken care of, evaluating $k$ is simply
5599 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
5601 <<kramers integ k function declaration>>=
5602 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
5604 <<kramers integ k function>>=
5605 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
5606 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
5607 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
5608 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
5612 <<kramers integ k structure create/destroy function declarations>>=
5613 void *string_create_kramers_integ_param_t(char **param_strings);
5614 void destroy_kramers_integ_param_t(void *p);
5617 <<kramers integ k structure create/destroy functions>>=
5618 kramers_integ_param_t *create_kramers_integ_param_t(double D,
5619 double xmin, double xmax, func_t *E,
5621 destroy_data_func_t *destroy_E_params)
5623 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
5624 assert(ret != NULL);
5629 ret->E_params = E_params;
5630 ret->destroy_E_params = destroy_E_params;
5634 void *string_create_kramers_integ_param_t(char **param_strings)
5636 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
5637 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
5638 void *E_params = string_create_spline_param_t(param_strings+1);
5639 return create_kramers_integ_param_t(
5640 safe_strtod(param_strings[0], __FUNCTION__),
5641 safe_strtod(param_strings[2], __FUNCTION__),
5642 safe_strtod(param_strings[3], __FUNCTION__),
5643 &spline_eval, E_params, destroy_spline_param_t);
5646 void destroy_kramers_integ_param_t(void *params)
5648 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5651 (*p->destroy_E_params)(p->E_params);
5657 Finally we have the GSL spline wrappers:
5658 <<spline functions>>=
5660 <<create/destroy spline>>
5663 <<spline function>>=
5664 double spline_eval(double x, void *spline_params)
5666 spline_param_t *p = (spline_param_t *)(spline_params);
5667 return gsl_spline_eval(p->spline, x, p->acc);
5671 <<create/destroy spline>>=
5672 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
5674 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
5675 assert(ret != NULL);
5676 ret->n_knots = n_knots;
5679 ret->acc = gsl_interp_accel_alloc();
5680 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
5681 gsl_spline_init(ret->spline, x, y, n_knots);
5685 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
5686 void *string_create_spline_param_t(char **param_strings)
5690 double *x=NULL, *y=NULL;
5691 /* break into ordered pairs */
5692 parse_list_string(param_strings[0], SEP, '(', ')',
5693 &num_knots, &knot_coords);
5694 x = (double *)malloc(sizeof(double)*num_knots);
5696 y = (double *)malloc(sizeof(double)*num_knots);
5698 for (i=0; i<num_knots; i++) {
5699 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
5700 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
5703 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
5705 free_parsed_list(num_knots, knot_coords);
5706 return create_spline_param_t(num_knots, x, y);
5709 void destroy_spline_param_t(void *params)
5711 spline_param_t *p = (spline_param_t *)params;
5714 gsl_spline_free(p->spline);
5716 gsl_interp_accel_free(p->acc);
5726 <<kramers integ k global declarations>>=
5727 extern int num_kramers_integ_params;
5728 extern char *kramers_integ_param_descriptions[];
5729 extern char kramers_integ_param_string[];
5732 <<kramers integ k globals>>=
5733 int num_kramers_integ_params = 4;
5734 char *kramers_integ_param_descriptions[] = {
5735 "diffusion D, m^2 / s",
5736 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
5737 "starting integration bound x_min, m",
5738 "ending integration bound x_max, m"};
5739 //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";
5740 //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";
5741 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";
5742 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
5743 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
5744 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
5747 <<kramers integ k model getopt>>=
5748 {"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}
5750 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
5751 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5753 \subsection{Unfolding model implementation}
5757 <<k model includes>>
5758 #include "k_model.h"
5759 <<k model internal definitions>>
5761 <<k model functions>>
5764 <<k model includes>>=
5765 #include <assert.h> /* assert() */
5766 #include <stdlib.h> /* NULL, malloc(), strto*() */
5767 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
5768 #include <stdio.h> /* fprintf(), stdout */
5769 #include <string.h> /* strlen(), strcpy() */
5770 #include <errno.h> /* errno, ERANGE */
5771 #include <gsl/gsl_integration.h>
5772 #include <gsl/gsl_spline.h>
5777 <<k model internal definitions>>=
5778 <<const k structure>>
5779 <<bell k structure>>
5780 <<kbell k structure>>
5781 <<kramers k structure>>
5782 <<kramers integ k structure>>
5783 <<spline param structure>>
5786 <<k model globals>>=
5791 <<kramers k globals>>
5792 <<kramers integ k globals>>
5795 <<k model functions>>=
5797 <<null k functions>>
5798 <<const k functions>>
5799 <<bell k functions>>
5800 <<kbell k functions>>
5801 <<kramers k functions>>
5802 <<kramers integ k functions>>
5805 \subsection{Unfolding model unit tests}
5807 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
5808 <<check-k-model.c>>=
5809 <<k model check includes>>
5810 <<check relative error>>
5812 <<k model test suite>>
5813 <<main check program>>
5816 <<k model check includes>>=
5817 #include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
5818 #include <stdio.h> /* sprintf() */
5819 #include <assert.h> /* assert() */
5820 #include <math.h> /* exp() */
5821 #include <errno.h> /* errno, ERANGE */
5824 #include "k_model.h"
5827 <<k model test suite>>=
5831 <<k model suite definition>>
5834 <<k model suite definition>>=
5835 Suite *test_suite (void)
5837 Suite *s = suite_create ("k model");
5838 <<const k test case defs>>
5839 <<bell k test case defs>>
5841 <<const k test case adds>>
5842 <<bell k test case adds>>
5847 \subsubsection{Constant}
5849 <<const k test case defs>>=
5850 TCase *tc_const_k = tcase_create("Constant k");
5853 <<const k test case adds>>=
5854 tcase_add_test(tc_const_k, test_const_k_create_destroy);
5855 tcase_add_test(tc_const_k, test_const_k_over_range);
5856 suite_add_tcase(s, tc_const_k);
5861 START_TEST(test_const_k_create_destroy)
5863 char *knot[] = {"1","2","3","4","5","6"};
5864 char *params[] = {knot[0]};
5867 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5868 params[0] = knot[i];
5869 p = string_create_const_k_param_t(params);
5870 destroy_const_k_param_t(p);
5875 START_TEST(test_const_k_over_range)
5877 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
5878 char *knot[] = {"1","2","3","4","5","6"};
5879 char *params[] = {knot[0]};
5882 char param_string[80];
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 for ( F=Fm; F<FM; F+=dF ) {
5889 fail_unless(const_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
5890 "const_k(%g, %g, &{%s}) = %g != %s",
5891 F, T, knot[i], const_k(F, &env, p), knot[i]);
5893 destroy_const_k_param_t(p);
5899 \subsubsection{Bell}
5901 <<bell k test case defs>>=
5902 TCase *tc_bell = tcase_create("Bell's k");
5905 <<bell k test case adds>>=
5906 tcase_add_test(tc_bell, test_bell_k_create_destroy);
5907 tcase_add_test(tc_bell, test_bell_k_at_zero);
5908 tcase_add_test(tc_bell, test_bell_k_at_one);
5909 suite_add_tcase(s, tc_bell);
5913 START_TEST(test_bell_k_create_destroy)
5915 char *knot[] = {"1","2","3","4","5","6"};
5917 char *params[] = {knot[0], dx};
5920 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5921 params[0] = knot[i];
5922 p = string_create_bell_param_t(params);
5923 destroy_bell_param_t(p);
5928 START_TEST(test_bell_k_at_zero)
5930 double F=0.0, T=300.0;
5931 char *knot[] = {"1","2","3","4","5","6"};
5933 char *params[] = {knot[0], dx};
5936 char param_string[80];
5939 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5940 params[0] = knot[i];
5941 p = string_create_bell_param_t(params);
5942 fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
5943 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
5944 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
5945 destroy_bell_param_t(p);
5950 START_TEST(test_bell_k_at_one)
5953 char *knot[] = {"1","2","3","4","5","6"};
5955 char *params[] = {knot[0], dx};
5956 double F= k_B*T/safe_strtod(dx, __FUNCTION__);
5959 char param_string[80];
5962 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5963 params[0] = knot[i];
5964 p = string_create_bell_param_t(params);
5965 CHECK_ERR(1e-6, safe_strtod(knot[i], __FUNCTION__)*exp(1.0), bell_k(F, &env, p));
5966 //fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__)*exp(1.0),
5967 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
5968 // F, T, knot[i], dx, bell_k(F, &env, p), safe_strtod(knot[i], __FUNCTION__)*exp(1.0));
5969 destroy_bell_param_t(p);
5975 <<kramers k tests>>=
5978 <<kramers k test case def>>=
5981 <<kramers k test case add>>=
5984 <<k model function tests>>=
5985 <<check relative error>>
5987 START_TEST(test_something)
5989 double k=5, x=3, last_x=2.0, F;
5990 one_dim_fn_t *handlers[] = {&hooke};
5991 void *data[] = {&k};
5992 double xi[] = {0.0};
5994 int new_active_groups = 1;
5995 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
5996 fail_unless(F = k*x, NULL);
6001 \subsection{Utilities}
6003 The unfolding models can get complicated, and are sometimes hard to explain to others.
6004 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
6005 the overhead of having to understand a full simulation.
6006 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
6007 or special mode, where the operation depends on the specific model selected.
6009 <<k-model-utils.c>>=
6011 <<k model utility includes>>
6012 <<k model utility definitions>>
6013 <<k model utility getopt functions>>
6014 <<k model utility multi model E>>
6015 <<k model utility main>>
6018 <<k model utility main>>=
6019 int main(int argc, char **argv)
6021 <<k model getopt array definition>>
6022 k_model_getopt_t *model = NULL;
6027 int num_param_args; /* for INIT_MODEL() */
6028 char **param_args; /* for INIT_MODEL() */
6030 double special_xmin;
6031 double special_xmax;
6032 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
6033 &Fmax, &special_xmin, &special_xmax, &flags);
6034 if (flags & VFLAG) {
6035 printf("#initializing model %s with parameters %s\n", model->name, model->params);
6037 INIT_MODEL("utility", model, model->params, params);
6040 if (model->k == NULL) {
6041 printf("No k function for model %s\n", model->name);
6045 printf("Fmax = %g <= 0\n", Fmax);
6051 printf("#F (N)\tk (%% pop. per s)\n");
6052 for (i=0; i<=N; i++) {
6053 F = Fmax*i/(double)N;
6054 k = (*model->k)(F, &env, params);
6056 printf("%g\t%g\n", F, k);
6061 if (model->k == NULL || model->k == &bell_k) {
6062 printf("No special mode for model %s\n", model->name);
6065 if (special_xmax <= special_xmin) {
6066 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
6070 double Xrng=(special_xmax-special_xmin), x, E;
6072 printf("#x (m)\tE (J)\n");
6073 for (i=0; i<=N; i++) {
6074 x = special_xmin + Xrng*i/(double)N;
6075 E = multi_model_E(model->k, params, &env, x);
6076 printf("%g\t%g\n", x, E);
6083 if (model->destructor)
6084 (*model->destructor)(params);
6089 <<k model utility multi model E>>=
6090 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
6092 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
6094 if (k == kramers_integ_k)
6095 E = (*p->E)(x, p->E_params);
6097 E = kramers_E(0, env, model_params, x);
6103 <<k model utility includes>>=
6106 #include <string.h> /* strlen() */
6107 #include <assert.h> /* assert() */
6108 #include <errno.h> /* errno, ERANGE */
6111 #include "string_eval.h"
6112 #include "k_model.h"
6115 <<k model utility definitions>>=
6116 <<version definition>>
6117 #define VFLAG 1 /* verbose */
6118 enum mode_t {M_K_OF_F, M_SPECIAL};
6119 <<string comparison definition and globals>>
6120 <<k model getopt definitions>>
6121 <<initialize model definition>>
6122 <<kramers k structure>>
6123 <<kramers integ k structure>>
6127 <<k model utility getopt functions>>=
6130 <<k model utility help>>
6131 <<k model utility get options>>
6134 <<k model utility help>>=
6135 void help(char *prog_name,
6137 int n_k_models, k_model_getopt_t *k_models,
6138 int k_model, double Fmax, double special_xmin, double special_xmax)
6141 printf("usage: %s [options]\n", prog_name);
6142 printf("Version %s\n\n", VERSION);
6143 printf("A utility for understanding the available unfolding models\n\n");
6144 printf("Environmental options:\n");
6145 printf("-T T\tTemperature T (currently %g K)\n", env->T);
6146 printf("-C T\tYou can also set the temperature T in Celsius\n");
6147 printf("Model options:\n");
6148 printf("-k model\tTransition rate model (currently %s)\n",
6149 k_models[k_model].name);
6150 printf("-K args\tTransition rate model argument string (currently %s)\n",
6151 k_models[k_model].params);
6152 printf("There are two output modes. In standard mode, k(F) is printed\n");
6153 printf("For example:\n");
6154 printf(" #Force (N)\tk (% pop. per s)\n");
6155 printf(" 123.456\t7.89\n");
6157 printf("In special mode, the output depends on the model.\n");
6158 printf("For models defining an energy function E(x), that function is printed\n");
6159 printf(" #Distance (m)\tE (J)\n");
6160 printf(" 0.0012\t0.0034\n");
6162 printf("-m\tChange output to standard mode\n");
6163 printf("-M\tChange output to special mode\n");
6164 printf("-F\tSet the maximum F value for the standard mode k(F) output (currently %g)\n", Fmax);
6165 printf("-x\tSet the minimum x value for the special mode E(x) output (currently %g)\n", special_xmin);
6166 printf("-X\tSet the maximum x value for the special mode E(x) output (currently %g)\n", special_xmax);
6167 printf("-V\tChange output to verbose mode\n");
6168 printf("-h\tPrint this help and exit\n");
6170 printf("Unfolding rate models:\n");
6171 for (i=0; i<n_k_models; i++) {
6172 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
6173 for (j=0; j < k_models[i].num_params ; j++)
6174 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
6175 printf("\t\tdefaults: %s\n", k_models[i].params);
6180 <<k model utility get options>>=
6181 void get_options(int argc, char **argv, environment_t *env,
6182 int n_k_models, k_model_getopt_t *k_models,
6183 enum mode_t *mode, k_model_getopt_t **model,
6184 double *Fmax, double *special_xmin, double *special_xmax,
6185 unsigned int *flags)
6187 char *prog_name = NULL;
6188 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
6190 extern char *optarg;
6191 extern int optind, optopt, opterr;
6195 /* setup defaults */
6197 prog_name = argv[0];
6198 env->T = 300.0; /* K */
6202 *Fmax = 1e-9; /* N */
6203 *special_xmax = 1e-8;
6204 *special_xmin = 0.0;
6206 while ((c=getopt(argc, argv, options)) != -1) {
6208 case 'T': env->T = safe_strtod(optarg, "-T"); break;
6209 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
6211 k_model = index_k_model(n_k_models, k_models, optarg);
6212 *model = k_models+k_model;
6215 k_models[k_model].params = optarg;
6217 case 'm': *mode = M_K_OF_F; break;
6218 case 'M': *mode = M_SPECIAL; break;
6219 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
6220 case 'x': *special_xmin = safe_strtod(optarg, "-x"); break;
6221 case 'X': *special_xmax = safe_strtod(optarg, "-X"); break;
6222 case 'V': *flags |= VFLAG; break;
6224 fprintf(stderr, "unrecognized option '%c'\n", optopt);
6225 /* fall through to default case */
6227 help(prog_name, env, n_k_models, k_models, k_model, *Fmax, *special_xmin, *special_xmax);
6236 \section{\LaTeX\ documentation}
6238 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
6239 The comment blocks are extracted (with nicely formatted code blocks), using
6240 <<latex makefile lines>>=
6241 $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6242 noweave -latex -index -delay $< > $@
6243 $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR)
6247 <<latex makefile lines>>=
6248 $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
6250 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
6251 $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
6252 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
6253 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
6254 mv $(BUILD_DIR)/sawsim.pdf $@
6256 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
6257 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
6258 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
6260 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
6261 <<latex makefile lines>>=
6263 rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
6264 $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
6265 $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.toc \
6266 $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib
6267 clean_latex : semi-clean_latex
6268 rm -f $(DOC_DIR)/sawsim.pdf
6273 [[make]] is a common utility on *nix systems for managing dependencies while building software.
6274 In this case, we have several \emph{targets} that we'd like to build:
6275 the main simulation program \prog;
6276 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
6277 the documentation [[sawsim.pdf]];
6278 and the [[Makefile]] itself.
6279 Besides the generated files, there is the utility target
6280 [[clean]] for removing all generated files except [[Makefile]].
6282 % [[dist]] for generating a distributable tar file.
6284 Extract the makefile with
6285 `[[notangle -Rmakefile src/sawsim.nw | sed 's/ /\t/' > Makefile]]'.
6286 The sed is needed because notangle replaces tabs with eight spaces,
6287 but [[make]] needs tabs.
6289 # Decide what directories to put things in
6294 # Note: Cannot use, for example, `./build', because make eats the `./'
6295 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
6297 # Modules (X.c, X.h) defined in the noweb file
6300 # Modules defined outside the noweb file
6301 FREE_SAWSIM_MODS = interp tavl
6303 <<list module makefile lines>>
6304 <<tension balance module makefile lines>>
6305 <<tension model module makefile lines>>
6306 <<k model module makefile lines>>
6307 <<parse module makefile lines>>
6308 <<string eval module makefile lines>>
6309 <<accel k module makefile lines>>
6311 SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS)
6313 # Everything needed for sawsim under one roof. sawsim.c must come first
6314 SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \
6315 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h)
6316 # Libraries needed to compile sawsim
6317 LIBS = gsl gslcblas m
6318 CHECK_LIBS = gsl gslcblas m check
6319 # The non-check binaries generated
6320 BINS = sawsim tension_model_utils k_model_utils sawsim_profile
6323 # Define the major targets
6324 all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
6326 view : $(DOC_DIR)/sawsim.pdf
6328 profile : $(BIN_DIR)/sawsim_profile | $(BIN_DIR)
6329 $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \
6330 -mnull -Mwlc -A0.39e-9,28e-9 -F8
6331 gprof $(BIN_DIR)/sawsim_profile gmon.out > $@
6332 check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim
6333 $(SHELL) -e -c 'for B in $^; do ./$$B; done'
6334 clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \
6335 clean_tension_model_utils \
6336 clean_k_model_utils clean_latex clean_check_sawsim
6337 rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \
6338 $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \
6339 $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
6340 $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
6341 $(BUILD_DIR)/global.h ./gmon.out
6342 $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
6344 # Various builds of sawsim
6345 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
6346 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6347 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
6348 gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6349 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
6350 gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6352 # Create the directories
6353 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
6356 # Copy over the source external to sawsim
6357 # Note: Cannot use, for example, `./build', because make eats the `./'
6358 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
6360 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6364 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6368 ## Basic source generated with noweb
6369 # The central files sawsim.c and global.h...
6370 $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6372 $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6373 notangle -Rglobal.h $< > $@
6374 # ... and the modules
6375 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6376 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6377 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6378 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6379 $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6380 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6381 # Note: I use `_' as a space in my file names, but noweb doesn't like
6382 # that so we use `-' to name the noweb roots and substitute here.
6384 ## Building the unit-test programs
6386 $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
6387 notangle -Rchecks $< > $@
6388 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
6389 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
6390 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR)
6391 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
6392 clean_check_sawsim :
6393 rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
6394 # ... and the modules
6396 $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
6397 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \
6398 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \
6399 $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6400 $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6401 $$(BUILD_DIR)/global.h | $(BIN_DIR)
6402 gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
6403 $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
6405 # todo: clean up the required modules to
6406 $(CHECK_BINS:%=clean_%) :
6407 rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
6409 # Cleaning up the modules
6411 $(SAWSIM_MODS:%=clean_%) :
6412 rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h)
6414 <<tension model utils makefile lines>>
6415 <<k model utils makefile lines>>
6416 <<latex makefile lines>>
6418 Makefile : $(SRC_DIR)/sawsim.nw
6419 notangle -Rmakefile $< | sed 's/ /\t/' > $@
6421 This is fairly self-explanatory. We compile a dynamically linked
6422 version ([[sawsim]]) and a statically linked version
6423 ([[sawsim_static]]). The static version is about 9 times as big, but
6424 it runs without needing dynamic GSL libraries to link to, so it's a
6425 better format for distributing to a cluster for parallel evaluation.
6429 \subsection{Hookean springs in series}
6430 \label{sec.math_hooke}
6433 The effective spring constant for $n$ springs $k_i$ in series is given by
6435 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6441 F &= k_i x_i &&\forall i \le n \\
6442 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
6443 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
6444 F &= k_1 x_1 = k_\text{eff} x \\
6445 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
6446 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6451 \addcontentsline{toc}{section}{References}
6452 \bibliography{sawsim}