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 \emph{major} 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$. I
159 dunno how that slipped by before. Now there are multi-domain tests in
160 testing to check for this sort of mistake.
162 <<version definition>>=
163 #define VERSION "0.8"
169 sawsim - program for simulating single molecule mechanical unfolding.
170 Copyright (C) 2008-2009, William Trevor King
172 This program is free software; you can redistribute it and/or
173 modify it under the terms of the GNU General Public License as
174 published by the Free Software Foundation; either version 3 of the
175 License, or (at your option) any later version.
177 This program is distributed in the hope that it will be useful, but
178 WITHOUT ANY WARRANTY; without even the implied warranty of
179 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
180 See the GNU General Public License for more details.
182 You should have received a copy of the GNU General Public License
183 along with this program; if not, write to the Free Software
184 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
187 The author may be contacted at <wking@drexel.edu> on the Internet, or
188 write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
189 Philadelphia PA 19104, USA.
202 The unfolding system is stored as a chain of \emph{domains} (Figure
203 \ref{fig.domain_chain}). Domains can be folded, globular protein
204 domains, unfolded protein linkers, AFM cantilevers, or other
205 stretchable system components. The system is modeled as graph of
206 possible domain states with transitions between them (Figure
207 \ref{fig.domain_states}).
211 \subfloat[][]{\label{fig.domain_chain}
212 \begin{tikzpicture}[->,node distance=2cm,font=\footnotesize]
213 \tikzstyle{every state}=[fill,draw=red!50,very thick,fill=red!20]
214 \node[state] (A) {domain 1};
215 \node[state] (B) [below of=A] {domain 2};
216 \node[state] (C) [below of=B] {.~.~.};
217 \node[state] (D) [below of=C] {domain $N$};
218 \node (S) [below of=D] {Surface};
219 \node (E) [above of=A] {};
221 \path[-] (A) edge (B)
222 (B) edge node [right] {Tension} (C)
225 \path[->,green] (A) edge node [right,black] {Extension} (E);
229 \subfloat[][]{\label{fig.domain_states}
230 \begin{tikzpicture}[->,node distance=2.5cm,shorten <=1pt,shorten >=1pt,font=\footnotesize]
231 \tikzstyle{every state}=[fill,draw=blue!50,very thick,fill=blue!20]
232 \node[state] (A) {cantilever};
233 \node[state] (C) [below of=A] {transition};
234 \node[state] (B) [left of=C] {folded};
235 \node[state] (D) [right of=C] {unfolded};
237 \path (B) edge [bend left] node [above] {$k_1$} (C)
238 (C) edge [bend left] node [below] {$k_1'$} (B)
239 edge [bend left] node [above] {$k_2$} (D)
240 (D) edge [bend left] node [below] {$k_2'$} (C);
243 \caption{\subref{fig.domain_chain} Extending a chain of domains. One end
244 of the chain is fixed, while the other is extended at a constant
245 speed. The domains are coupled with rigid linkers, so the domains
246 themselves must stretch to accomodate the extension.
247 \subref{fig.domain_states} Each domain exists in a discrete state. At
248 each timestep, it may transition into another state following a
249 user-defined state matrix such as this one, showing a metastable
250 transition state and an explicit ``cantilever'' domain.}
254 Each domain contributes to the total tension, which depends on the
255 chain extension. The particular model for the tension as a function
256 of extension is handled generally with function pointers. So far the
257 following models have been implemented:
259 \item Null (Appendix \ref{sec.null_tension}),
260 \item Constant (Appendix \ref{sec.const_tension}),
261 \item Hookean spring (Appendix \ref{sec.hooke}),
262 \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
263 \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
266 The tension model and parameters used for a particular domain depend
267 on the domain's current state. The transition rates between states
268 are also handled generally with function pointers, with
271 \item Null (Appendix \ref{sec.null_k}),
272 \item Bell's model (Appendix \ref{sec.bell}),
273 \item Spring constant corrected Bell model (Appendix \ref{sec.kbell}),
274 \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
275 \item Kramers' saddle point model (Appendix \ref{sec.kramers}).
278 The unfolding of the chain domains is modeled in two stages. First
279 the tension of the chain is calculated. Then the tension (assumed to
280 be constant over the length of the chain, see Section
281 \ref{sec.timescales}), is applied to each domain, and used to
282 calculate the probability of that domain changing state in the next
283 timestep $dt$. Then the time is stepped forward, and the process is
284 repeated. The simulation is complete when a pre-set time limit
285 [[t_max]] is reached or a pre-selected state [[stop_state]] is empty.
287 int main(int argc, char **argv)
289 <<tension model getopt array definition>>
290 <<k model getopt array definition>>
292 /* variables initialized in get_options() */
293 list_t *state_list=NULL, *transition_list=NULL;
294 environment_t env={0};
295 double P, t_max, dt_max, v, x_step;
296 state_t *stop_state=NULL;
298 /* variables used in the simulation loop */
299 double t=0, x=0, dt, F; /* time, distance, timestep, force */
300 int transition=1; /* boolean marking a transition for tension optimization */
302 get_options(argc, argv, &P, &t_max, &dt_max, &v, &x_step, &stop_state, &env,
303 NUM_TENSION_MODELS, tension_models, NUM_K_MODELS, k_models,
304 &state_list, &transition_list);
307 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Distance (m)\tForce (N)\tStiffness (N/m)\n");
308 if (flags & FLAG_OUTPUT_TRANSITION_FORCES) printf("#Force (N)\tInitial state\tFinal state\n");
309 while ((t_max > 0 && t <= t_max) || (stop_state != NULL && stop_state->num_domains > 0)) {
310 F = find_tension(state_list, &env, x, transition);
312 dt = determine_dt(state_list, transition_list, &env, x, P, dt_max, v);
314 dt = determine_dt(state_list, transition_list, &env, x, P, dt_max, 0);
315 transition=random_transitions(transition_list, F, dt, &env);
316 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness);
321 if (dt == dt_max) { /* step completed */
324 } else { /* still working on this step */
329 destroy_state_list(state_list);
330 destroy_transition_list(transition_list);
334 @ The meat of the simulation is bundled into the three functions
335 [[find_tension]], [[determine_dt]], and [[random_transitions]].
336 [[find_tension]] is discussed in Section \ref{sec.find_tension},
337 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
338 [[random_transitions]] in Section \ref{sec.transition_rate}.
340 The stretched distance is extended in one of two modes depending on
341 [[xstep]]. If [[xstep == 0]], then the pulling is continuous, at
342 least within the limits of the inherent discretization of a time
343 stepped approach. At any rate, the timestep is limited by the maximum
344 allowed timestep [[dt_max]] and the maximum allowed unfolding
345 probability in a given step [[P]]. In the second mode [[xstep > 0]],
346 and the end to end distance increases in discrete steps of that
347 length. The time between steps is chosen to maintain an average
348 unfolding velocity [[v]]. This approach is less work to simulate than
349 smooth pulling and also closer to the reality of many experiments, but
350 it is practically impossible to treat analytically. With both pulling
351 styles implemented, the effects of the discretization can be easily
352 calculated, bridging the gap between theory and experiment.
354 Environmental parameters important in determining reaction rates and
355 tensions (e.g.\ temperature, pH) are stored in a single structure to
356 facilitate extension to more complicated models in the future.
357 <<environment definition>>=
358 typedef struct environment_struct {
368 #define DOUBLE_PRECISION 1e-12
370 <<environment definition>>
371 <<one dimensional function definition>>
372 <<create/destroy definitions>>
374 #endif /* GLOBAL_H */
376 Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
377 included multiple times.
379 \section{Simulation functions}
381 <<simulation functions>>=
385 <<does domain transition>>
386 <<randomly transition domains>>
390 \label{sec.find_tension}
392 Because the stretched system may be made up of several parts (folded
393 domains, unfolded domains, spring-like cantilever, \ldots), we will
394 assign the domains to states with tension models and groups. The
395 models are used to determine the tension of all the domain in a given
396 state following some model (Hookean spring, WLC, \ldots). The domains
397 are assumed to be commutative, so ordering is ignored. The
398 interactions between the groups are assumed to be linear, but the
399 interactions between domains of the same group need not be. Each
400 model has a tension handler function, which gives the tension $F$
401 needed to stretch a given group of domains a distance $x$.
403 Groups represent collections of states which the model should treat as
404 a single entity. To understand the purpose of groups, consider a
405 system with two unfolded domain states (e.g.\ for two different
406 proteins) that were both modeled as WLCs. If both unfolded states had
407 the same persistence length, it would be wise to assign them to the
408 same group. This would reduce both the computational cost of
409 balancing the tension and the error associated with the inter-group
410 interaction linearity. Note that group numbers only matter
411 \emph{within} model classes, since grouping domains with seperate
412 models doesn't make sense. Within designated groups, the tension
413 parameters for each domain are still checked for compatibility, so if
414 you accidentally assign incompatible domains to the same group you
415 should get an appropriate error message.
417 <<find tension definitions>>=
418 #define NUM_TENSION_MODELS 5
422 <<tension handler array global declaration>>=
423 extern one_dim_fn_t *tension_handlers[];
426 <<tension handler array global>>=
427 one_dim_fn_t *tension_handlers[] = {
429 &const_tension_handler,
437 <<tension model global declarations>>=
438 <<tension handler array global declaration>>
441 <<tension handler types>>=
442 typedef struct tension_handler_data_struct {
443 list_t *group_tension_params; /* one item for each domain in the group */
446 } tension_handler_data_t;
450 After sorting the chain into separate groups $G_i$, with tension
451 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
452 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
453 \forall i,j$. Note that there may be several states within each
454 group. [[find_tension]] is basically a wrapper to feed proper input
455 into the [[tension_balance]] function. [[transition]] should be set
456 to 0 if there were no transitions since the previous call to
457 [[find_tension]] to support some optimizations that assume a small
458 increase in tension between steps (only true if no transition has
459 occured). While we're messing with the tension handlers, we also grab
460 a measure of the current linker stiffness for the environmental
461 variable, which is needed by some of the unfolding rate models
462 (e.g.\ the linker-stiffness corrected Bell model, Appendix
465 double find_tension(list_t *state_list, environment_t *env, double x, int transition)
467 static int num_active_groups=0;
468 static one_dim_fn_t **per_group_handlers = NULL;
469 static one_dim_fn_t **per_group_inverse_handlers = NULL;
470 static void **per_group_data = NULL; /* array of pointers to tension_handler_data_t */
471 static double last_x;
472 tension_handler_data_t *tdata;
477 fprintf(stderr, "%s:\tfinding tension at distance %g%s\n", __FUNCTION__, x, transition ? " after a transition" : "");
480 if (transition != 0 || num_active_groups == 0) { /* setup statics */
481 /* get new statics, freeing the old ones if they aren't NULL */
482 get_active_groups(state_list, &num_active_groups, &per_group_handlers, &per_group_inverse_handlers, &per_group_data);
484 /* fill in the group handlers and params */
486 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]);
488 for (i=0; i<num_active_groups; i++) {
489 tdata = (tension_handler_data_t *) per_group_data[i];
491 fprintf(stderr, "%s:\tactive group %d of %d (data %p, param list %p)", __FUNCTION__, i, num_active_groups, tdata, tdata->group_tension_params);
492 list_print(stderr, tdata->group_tension_params, " parameters");
495 /* tdata->persist continues unchanged */
498 } /* else continue with the current statics */
500 F = tension_balance(num_active_groups, per_group_handlers, per_group_inverse_handlers, (void **)per_group_data, last_x, x, &env->stiffness);
506 @ For the moment, we will restrict our group boundaries to a single
507 dimension, so $\sum_i x_i = x$, our total extension (it is this
508 restriction that causes the groups to interact linearly). We'll also
509 restrict our extensions to all be positive. With these restrictions,
510 the problem of balancing the tensions reduces to solving a system of
511 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
512 the number of active groups. In general this can be fairly
513 complicated, but without much loss of practicality we can restrict
514 ourselves to strictly monotonically increasing, non-negative tension
515 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
516 simpler. For example, it guarantees the existence of a unique, real
517 solution for finite forces. See Appendix \ref{sec.tension_balance}
518 for the tension balancing implementation.
520 Each group must define a parameter structure if the tension model is
522 a function to create the parameter structure,
523 a function to destroy the parameter structure, and
524 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as its data argument.
525 The tension-specific parameter structures are contained in the domain
526 group field. For optimization, a group may also define an inverse
527 tension function as an optimization. See the Section
528 \ref{sec.model_selection} for a discussion on the command line
529 framework and inverse function discussion. See the worm-like chain
530 implementation for an example (Section \ref{sec.wlc}).
532 The major limitation of our approach is that all of our tension models
533 are Markovian (memory-less), which is not really a problem for the
534 chains (see \citet{evans99} for WLC thermalization rate calculations).
536 \subsection{Transition rate}
537 \label{sec.transition_rate}
539 Each state transition is modeled as unimolecular, first order reaction
541 \text{State 1} \xrightarrow{k} \text{State 2}
543 With the rate constant $k$ defined by
545 \frac{d[\text{State 2}]}{dt} = -\frac{d[\text{State 1}]}{dt} = k [\text{State 1}].
547 So the probability of a given domain transitioning in a short time
553 For a group of $N$ identical domains, and therefore identicalunfolding
554 probabilities $P$, the probability of $n$ domains unfolding in a given
555 timestep is given by the binomial distribution
557 p(n) = \frac{N!}{n!(N-n)!}P^n(1-P)^(N-n).
559 The probability that \emph{none} of these domains unfold is then
563 so the probability that \emph{at least one} domain unfolds is
567 We take some time to discuss the meaning of $p(n>1)$
568 (i.e. multi-unfolding timesteps) in Section \ref{sec.timescales}.
570 <<does domain transition>>=
571 int domain_transitions(double F, double dt, environment_t *env, int num_domains, transition_t *transition)
572 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, number of 'reactant' domains, pointer to transition structure */
574 k = accel_k(transition->k, F, env, transition->k_params);
575 //(*transition->k)(F, env, domain->k_params);
576 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
577 return happens(1-pow((1.0-k*dt), num_domains)); /* N dice rolls for prob. k*dt event */
579 @ [[happens]] is a random decision making function defined in Appendix \ref{sec.utils}.
581 <<randomly transition domains>>=
582 int random_transitions(list_t *transition_list, double tension,
583 double dt, environment_t *env)
584 { /* returns 1 if there was a transition and 0 otherwise */
585 while (transition_list != NULL) {
586 if (T(transition_list)->initial_state->num_domains > 0
587 && domain_transitions(tension, dt, env, T(transition_list)->initial_state->num_domains, T(transition_list))) {
588 if (flags & FLAG_OUTPUT_TRANSITION_FORCES)
589 fprintf(stdout, "%g\t%s\t%s\n", tension, T(transition_list)->initial_state->name, T(transition_list)->final_state->name);
590 T(transition_list)->initial_state->num_domains--;
591 T(transition_list)->final_state->num_domains++;
592 return 1; /* our one domain has transitioned, stop looking for others */
594 transition_list = transition_list->next;
600 \subsection{Timescales}
601 \label{sec.timescales}
603 The simulation assumes that chain equilibration is instantaneous
604 relative to the simulation timestep $dt$, so that tension is uniform
605 along the chain. The quality of this assumption depends on your
606 particular chain. For example, a damped spring thermalizes on a
607 timescale of order $\tau = 1/\gamma$ where the damping ratio $\gamma
608 \equiv \omega_0(-\zeta \pm \sqrt{\zeta^2-1}$, the natural angular
609 frequency $omega_0 \equiv \sqrt{k/m}$, $\zeta \equiv b/2\sqrt{km}$,
610 and $b$ sets the drag $F_b = -bv$. For our cantilevers $\tau$ is
611 on the order of milliseconds, which is longer than a timestep.
612 However, the approximation is still reasonable, because there is
613 rarely unfolding during the cantilever return. You could, of course,
614 take cantilever, WLC, etc.\ response times into effect, but that
615 would significantly complicate a simulation that seems to work
616 well enough without bothering :p. For WLC and FJC relaxation times,
617 see the Appendix of \citet{evans99} and \citet{kroy07}.
619 Besides assuming our timestep is much greater than equilibration
620 times, we also force it to be short enough so that only one domain may
621 unfold in a given timestep. For an unfolding timescale of $dt_u =
622 1/k$, we adapt our timesteps to keep $dt \ll dt_u$, so the probability
623 of two domains unfolding in the same timestep is small (see
624 [[determine_dt]] in Section \ref{sec.adaptive_dt} and [[P]] in
625 [[main()]] in Section \ref{sec.main} set by [[-P]] in
626 [[get_options()]] in Section \ref{sec.get_options}). This approach
627 breaks down as the adaptive timestep scheme approaches the transition
628 timestep $dt_t$, but $dt_t \sim 1\U{$\mu$s}$ for Ig27-like proteins
629 \citep{klimov00}, so this shouldn't be a problem. To reassure
630 yourself, you can ask the simulator to print the smallest timestep
631 that was used with TODO.
633 Even if the likelyhood of two domains unfolding in the same timestep
634 is small, if you run long enough simulations it will eventually occur.
635 In this case, we assume that $dt_t \ll dt$, so even if two domains
636 would unfold if we stuck to our unfolding rate $k$ for an entire
637 timestep $dt$, in ``reality'' only one of those domains unfolds.
638 Because we do not know when in the timestep the unfolding took place,
639 we move on to the next timestep without checking for further
640 unfoldings. This ``unchecked timestep portion'' should not contribute
641 significant errors to the calculation, because if $dt$ was small
642 enough that unfolding was unlikely at the pre-unfolding force, it
643 should be even more unlikely at the post-unfolding force, especially
644 over only a fraction of the timestep. The only dangerous cases for
645 the current implementation are those where unfolding intermediates are
646 much less stable than their precursor states, in which case an
647 unfolding event during the remaining timestep may actually be likely.
648 A simple workaround for such cases is to pick the value for [[P]] to
649 be small enough that the $dt \ll dt_u$ assumption survives even under
650 a transition populating the unstable intermediate.
652 \subsection{Adaptive timesteps}
653 \label{sec.adaptive_dt}
655 We'd like to pick $dt$ so the probability of changing state
656 (e.g.\ unfolding) in the next timestep is small. If we don't adapt our
657 timestep, we also risk breaking our approximation that $F(x' \in [x,
658 x+v \cdot dt]) = F(x)$ or that only one domain changes state in a
659 given timestep. The simulation would have been accurate for
660 sufficiently small fixed timesteps, but adaptive timesteps will allow
661 us to move through low tension regions in fewer steps, leading to a
662 more efficient simulation.
664 Consider the two state folded $\rightarrow$ unfolded transition.
665 Because $F(x)$ is monotonically increasing between unfoldings,
666 excessively large timesteps will lead to erroneously small unfolding
667 rates (an thus, higher average unfolding force).
669 The actual adaptive timestep implementation is not particularly
670 interesting, since we are only required to reduce $dt$ to somewhere
671 below a set threshold, so I've removed it to Appendix
672 \ref{sec.adaptive_dt_implementation}.
678 The models provide the physics of the simulation, but the simulation
679 allows interchangeable models, and we are currently focusing on the
680 simulation itself, so we remove the actual model implementation to the
683 The tension models are defined in Appendix \ref{sec.tension_models}
684 and unfolding models are defined in Appendix \ref{sec.k_models}.
687 #define k_B 1.3806503e-23 /* J/K */
691 \section{Command line interface}
693 <<get option functions>>=
695 <<index tension model>>
700 \subsection{Model selection}
701 \label{sec.model_selection}
703 The main difficulty with the command line interface in \prog\ is
704 developing an intuitive method for accessing the possibly large number
705 of available models. We'll treat this generally by defining an array
706 of available models, containing their important parameters.
708 <<get options definitions>>=
709 <<model getopt definitions>>
712 <<create/destroy definitions>>=
713 typedef void *create_data_func_t(char **param_strings);
714 typedef void destroy_data_func_t(void *param_struct);
717 <<model getopt definitions>>=
718 <<tension model getopt definitions>>
719 <<k model getopt definitions>>
722 In order to choose models, we need a method of assembling a reaction
723 graph such as that in Figure \ref{fig.domain_states} and population
724 the graph with a set of domains. First we declare the domain states
727 [[-s <state-name>,<model>[,<group>] -S <model-params> [-N <initial-population]]] \\
731 [[-s unfolded,wlc,1 -S 1e-8,4e-10 -N 5]]
733 This creates the state named [[unfolded]], using the WLC model and one
734 as a group number (defaults to 0, see Section \ref{sec.find_tension}).
735 The tension parameters are then set to [[1e-8,4e-10]], which in the
736 case of the WLC are the contour and persistence lengths respectively.
737 Finally we populate the state by adding five domains to the state.
738 The name of the state is importent for identifying states later.
740 Once the domains have been created and populated, you can added
741 transition rates following
743 [[-k <initial-state>,<final-state>,<model> -K <model-params>]] \\
747 [[-k folded,unfolded,bell -K 3.3e-4,0.25e-9]]
749 This creates a reaction going from the [[folded]] state to the
750 [[unfolded]] state, with the rate constant given by the Bell model
751 (Appendix \ref{sec.bell}). The unfolding parameters are then set to
752 [[3.3e-4,0.25e-9]], which in the case of the Bell model are the
753 unforced rate and transition state distance respectively.
755 \subsubsection{Tension}
757 To access the various tension models from the command line, we define
758 a structure that contains all the neccessary information about the
760 <<tension model getopt definitions>>=
761 typedef struct tension_model_getopt_struct {
762 one_dim_fn_t *handler; /* fn implementing the model on a group */
763 one_dim_fn_t *inverse_handler; /* fn implementing inverse on group */
764 char *name; /* name string identifying the model */
765 char *description; /* a brief description of the model */
766 int num_params; /* number of extra params for the model */
767 char **param_descriptions; /* descriptions of extra params */
768 char *params; /* default values for the extra params */
769 create_data_func_t *creator; /* param-string -> model-data-struct fn */
770 destroy_data_func_t *destructor; /* fn to free the model data structure */
771 } tension_model_getopt_t;
772 @ Set [[inverse_handler]] to [[NULL]] if you wish to invert using
773 bisection. Obviously, this will be slower than computing an
774 analytical inverse and more prone to rounding errors, but it may be
775 easier if a closed-form inverse does not exist.
777 <<tension model getopt array definition>>=
778 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
779 <<null tension model getopt>>,
780 <<constant tension model getopt>>,
781 <<hooke tension model getopt>>,
782 <<worm-like chain tension model getopt>>,
783 <<freely-jointed chain tension model getopt>>
787 \subsubsection{Unfolding rate}
789 <<k model getopt definitions>>=
790 #define NUM_K_MODELS 6
792 typedef struct k_model_getopt_struct {
797 char **param_descriptions;
799 create_data_func_t *creator;
800 destroy_data_func_t *destructor;
804 <<k model getopt array definition>>=
805 k_model_getopt_t k_models[NUM_K_MODELS] = {
806 <<null k model getopt>>,
807 <<const k model getopt>>,
808 <<bell k model getopt>>,
809 <<kbell k model getopt>>,
810 <<kramers k model getopt>>,
811 <<kramers integ k model getopt>>
818 void help(char *prog_name, double P, double t_max, double dt_max, double v, double x_step,
819 state_t *stop_state, environment_t *env,
820 int n_tension_models, tension_model_getopt_t *tension_models,
821 int n_k_models, k_model_getopt_t *k_models,
822 int k_model, list_t *state_list)
827 if (state_list != NULL)
828 state = S(tail(state_list));
830 printf("usage: %s [options]\n", prog_name);
831 printf("Version %s\n\n", VERSION);
832 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
834 printf("Simulation options:\n");
836 printf("-q state\tStop simulation when this state empties (- for no limit, currently %s)\n", stop_state ? stop_state->name : "-" );
837 printf("-t tmax\tMaximum simulated time (-1 for no limit, currently %g s)\n", t_max);
838 printf("-d dt\tMaximum allowed timestep dt (currently %g s)\n", dt_max);
839 printf("-P P\tTarget probability for dt (currently %g)\n", P);
840 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
841 printf("-x dx\tPulling step size (currently %g m)\n", x_step);
842 printf("\tWhen dx = 0, the pulling is continuous.\n");
843 printf("\tWhen dx > 0, the pulling occurs in discrete steps.\n");
845 printf("Environmental options:\n");
846 printf("-T T\tTemperature T (currently %g K)\n", env->T);
847 printf("-C T\tYou can also set the temperature T in Celsius\n");
849 printf("State creation options:\n");
850 printf("-s name,model[[,group],params]\tCreate a new state.\n");
851 printf("Once you've set up the state, you may populate it with domains\n");
852 printf("-N n\tAdd n domains in the current state (%s)\n", state ? state->name : "-");
854 printf("Transition creation options:\n");
855 printf("-k initial_state,final_state,model[,params]\tCreate a new transition\n");
857 printf("To double check your reaction graph, you can produce graphviz dot output\n");
858 printf("describing the current states and transitions.\n");
859 printf("-D\tPrint dot description of states and transitions and exit.\n");
861 printf("Simulation output mode options:\n");
862 printf("There are two output modes. In standard mode, only the transition\n");
863 printf("events are printed. For example:\n");
864 printf(" #Force (N)\tInitial state\tFinal state\n");
865 printf(" 123.456e-12\tfolded\tunfolded\n");
867 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
868 printf(" #Distance (m)\tForce (N)\tStiffness (N/m)\n");
869 printf(" 0.001\t0.0023\t0.002\n");
871 printf("-V\tChange output to verbose mode.\n");
872 printf("-h\tPrint this help and exit.\n");
875 printf("Tension models:\n");
876 for (i=0; i<n_tension_models; i++) {
877 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
878 for (j=0; j < tension_models[i].num_params ; j++)
879 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
880 printf("\t\tdefaults: %s\n", tension_models[i].params);
882 printf("Unfolding rate models:\n");
883 for (i=0; i<n_k_models; i++) {
884 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
885 for (j=0; j < k_models[i].num_params ; j++)
886 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
887 printf("\t\tdefaults: %s\n", k_models[i].params);
890 printf("Examples:\n");
891 printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded. Stop once there are no more folded domains.\n");
892 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);
893 printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded. Stop after 0.25 seconds\n");
894 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);
898 \subsection{Get options}
899 \label{sec.get_options}
903 void get_options(int argc, char **argv,
904 double *pP, double *pT_max, double *pDt_max, double *pV,
905 double *pX_step, state_t **pStop_state, environment_t *env,
906 int n_tension_models, tension_model_getopt_t *tension_models,
907 int n_k_models, k_model_getopt_t *k_models,
908 list_t **pState_list, list_t **pTransition_list)
910 char *prog_name = NULL;
911 char c, options[] = "q:t:d:P:v:x:T:C:s:N:k:DVh";
912 int tension_model=0, k_model=0, initial_state=-1, final_state=-1;
913 char *state_name=NULL;
914 int i, n, tension_group=0;
915 list_t *temp_list=NULL;
918 void *model_params; /* for INIT_MODEL() */
919 int num_param_args; /* for INIT_MODEL() */
920 char **param_args; /* for INIT_MODEL() */
922 extern int optind, optopt, opterr;
927 for (i=0; i<n_tension_models; i++) {
928 fprintf(stderr, "tension model %i %s:\t model's handler %p\t inverse handler %p\n",
929 i, tension_models[i].name, tension_models[i].handler, tension_models[i].inverse_handler);
930 assert(tension_models[i].handler == tension_handlers[i]);
935 flags = FLAG_OUTPUT_TRANSITION_FORCES;
938 *pT_max = -1; /* s, -1 removes this limit */
939 *pDt_max = 0.001; /* s */
940 *pP = 1e-3; /* % pop per s */
941 *pV = 1e-6; /* m/s */
942 *pX_step = 0.0; /* m */
943 env->T = 300.0; /* K */
945 *pTransition_list = NULL;
947 while ((c=getopt(argc, argv, options)) != -1) {
950 if (STRMATCH(optarg, "-")) {
953 temp_list = head(*pState_list);
954 while (temp_list != NULL) {
955 if (STRMATCH(S(temp_list)->name, optarg)) {
956 *pStop_state = S(temp_list);
959 temp_list = temp_list->next;
963 case 't': *pT_max = safe_strtod(optarg, "-t"); break;
964 case 'd': *pDt_max = safe_strtod(optarg, "-d"); break;
965 case 'P': *pP = safe_strtod(optarg, "-P"); break;
966 case 'v': *pV = safe_strtod(optarg, "-v"); break;
967 case 'x': *pX_step = safe_strtod(optarg, "-x"); break;
968 case 'T': env->T = safe_strtod(optarg, "-T"); break;
969 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
971 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
972 assert(num_strings >= 2 && num_strings <= 4);
974 state_name = strings[0];
975 if (state_index(*pState_list, state_name) != -1) {
976 fprintf(stderr, "Invalid option '%s', repeated state name '%s'\n", optarg, state_name);
979 tension_model = index_tension_model(n_tension_models, tension_models, strings[1]);
980 if (num_strings == 4)
981 tension_group = safe_strtoi(strings[2], optarg);
984 if (num_strings >= 3)
985 tension_models[tension_model].params = strings[num_strings-1];
987 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);
989 INIT_MODEL(state_name, &tension_models[tension_model], tension_models[tension_model].params, model_params);
990 push(pState_list, create_state(state_name,
991 tension_models[tension_model].handler,
992 tension_models[tension_model].inverse_handler,
995 tension_models[tension_model].destructor,
997 *pState_list = tail(*pState_list);
999 free_parsed_list(num_strings, strings);
1002 n = safe_strtoi(optarg, "-N");
1004 fprintf(stderr, "%s:\tadding %d domains to %s\n", __FUNCTION__, n, S(*pState_list)->name);
1006 S(*pState_list)->num_domains += n;
1009 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
1010 assert(num_strings >= 3 && num_strings <= 4);
1011 initial_state = state_index(*pState_list, strings[0]);
1012 final_state = state_index(*pState_list, strings[1]);
1013 k_model = index_k_model(n_k_models, k_models, strings[2]);
1015 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);
1017 assert(initial_state != final_state);
1018 if (num_strings == 4)
1019 k_models[k_model].params = strings[3];
1020 INIT_MODEL(state_name, &k_models[k_model], k_models[k_model].params, model_params);
1021 push(pTransition_list, create_transition(list_element(*pState_list, initial_state),
1022 list_element(*pState_list, final_state),
1023 k_models[k_model].k,
1025 k_models[k_model].destructor), 1);
1026 free_parsed_list(num_strings, strings);
1029 print_state_graph(stdout, *pState_list, *pTransition_list);
1033 flags = FLAG_OUTPUT_FULL_CURVE;
1036 fprintf(stderr, "unrecognized option '%c'\n", optopt);
1037 /* fall through to default case */
1039 help(prog_name, *pP, *pT_max, *pDt_max, *pV, *pX_step, *pStop_state, env, n_tension_models, tension_models, n_k_models, k_models, k_model, *pState_list);
1043 /* check the arguments */
1044 assert(*pState_list != NULL); /* no domains! */
1045 assert(*pP > 0.0); assert(*pP < 1.0);
1046 assert(*pDt_max > 0.0);
1048 assert(env->T > 0.0);
1050 crosslink_groups(*pState_list, *pTransition_list);
1056 <<index tension model>>=
1057 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
1060 for (i=0; i<n_models; i++)
1061 if (STRMATCH(models[i].name, name))
1063 if (i == n_models) {
1064 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
1072 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
1075 for (i=0; i<n_models; i++)
1076 if (STRMATCH(models[i].name, name))
1078 if (i == n_models) {
1079 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
1086 <<initialize model definition>>=
1087 /* requires int num_param_args and char **param_args in the current scope
1089 * INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
1090 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
1092 #define INIT_MODEL(role, model, param_string, param_pointer) \
1094 parse_list_string((param_string), SEP, '{', '}', \
1095 &num_param_args, ¶m_args); \
1096 if (num_param_args != (model)->num_params) { \
1098 "%s model %s expected %d params,\n", \
1099 (role), (model)->name, (model)->num_params); \
1101 "not the %d params in '%s'\n", \
1102 num_param_args, (param_string)); \
1103 assert(num_param_args == (model)->num_params); \
1105 if ((model)->creator) \
1106 param_pointer = (*(model)->creator)(param_args); \
1108 param_pointer = NULL; \
1109 free_parsed_list(num_param_args, param_args); \
1113 First we define some safe string parsing functions. Currently these
1114 abort on any irregularity, but in the future they could print error
1115 messages, etc. [[static]] because the functions are currently
1116 defined in each [[*.c]] file that uses them, but perhaps they should
1117 be moved to [[utils.h/c]] or some such instead\ldots
1120 static int safe_strtoi(const char *s, const char *description) {
1124 ret = strtol(s, &endp, 10);
1125 if (endp[0] != '\0') { //strlen(endp) > 0) {
1126 fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1127 endp, s, description);
1128 assert(1==0); //strlen(endp) == 0);
1130 fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1133 } else if (errno == ERANGE) {
1134 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1135 assert(errno != ERANGE);
1140 static double safe_strtod(const char *s, const char *description) {
1144 ret = strtod(s, &endp);
1145 if (endp[0] != '\0') { //strlen(endp) > 0) {
1146 fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1147 endp, s, description);
1148 assert(1==0); //strlen(endp) == 0);
1150 fprintf(stderr, "Error: unparsable empty string '%s' for %s\n",
1153 } else if (errno == ERANGE) {
1154 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1155 assert(errno != ERANGE);
1164 \addcontentsline{toc}{section}{Appendicies}
1166 \section{sawsim.c details}
1170 The general layout of our simulation code is:
1181 We include [[math.h]], so don't forget to link to the libm with `-lm'.
1183 #include <assert.h> /* assert() */
1184 #include <stdlib.h> /* malloc(), free(), rand(), strto*() */
1185 #include <stdio.h> /* fprintf(), stdout */
1186 #include <string.h> /* strlen, strtok() */
1187 #include <math.h> /* exp(), M_PI, sqrt() */
1188 #include <sys/types.h> /* pid_t (returned by getpid()) */
1189 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
1190 #include <errno.h> /* errno, ERANGE (for safe_strto*()) */
1193 #include "tension_balance.h"
1194 #include "k_model.h"
1195 #include "tension_model.h"
1197 #include "accel_k.h"
1201 <<version definition>>
1202 <<flag definitions>>
1203 <<max/min definitions>>
1204 <<string comparison definition and globals>>
1205 <<initialize model definition>>
1206 <<get options definitions>>
1207 <<state definitions>>
1208 <<transition definitions>>
1217 <<create/destroy state>>
1218 <<destroy state list>>
1220 <<create/destroy transition>>
1221 <<destroy transition list>>
1222 <<print state graph>>
1224 <<simulation functions>>
1225 <<boilerplate functions>>
1228 <<boilerplate functions>>=
1230 <<get option functions>>
1236 srand(getpid()*time(NULL)); /* seed rand() */
1240 <<flag definitions>>=
1241 /* in octal b/c of prefixed '0' */
1242 #define FLAG_OUTPUT_FULL_CURVE 01
1243 #define FLAG_OUTPUT_TRANSITION_FORCES 02
1247 static unsigned long int flags = 0;
1250 \subsection{Utilities}
1253 <<max/min definitions>>=
1254 #define MAX(a,b) ((a)>(b) ? (a) : (b))
1255 #define MIN(a,b) ((a)<(b) ? (a) : (b))
1258 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
1259 <<string comparison definition and globals>>=
1260 // Check if two strings match, return 1 if they do
1261 static char *temp_string_A;
1262 static char *temp_string_B;
1263 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
1264 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
1265 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
1266 /* +1 to also compare the '\0' */
1269 We also define a macro for our [[check]] unit testing
1270 <<check relative error>>=
1271 #define CHECK_ERR(max_err, expected, received) \
1273 fail_unless( (received-expected)/expected < max_err, \
1274 "relative error %g >= %g in %s (Expected %g, received %g)", \
1275 (received-expected)/expected, max_err, #received, \
1276 expected, received); \
1277 fail_unless(-(received-expected)/expected < max_err, \
1278 "relative error %g >= %g in %s (Expected %g, received %g)", \
1279 -(received-expected)/expected, max_err, #received, \
1280 expected, received); \
1285 int happens(double probability)
1287 assert(probability >= 0.0); assert(probability <= 1.0);
1288 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*/
1293 \subsection{Adaptive timesteps}
1294 \label{sec.adaptive_dt_implementation}
1296 $F(x)$ increases with $x$, possibly exploding (e.g.\ in the worm-like
1297 chain model), so basing the timestep on the the unfolding probability
1298 at the current tension is dangerous, and we need to search for a $dt$
1299 for which $P(F(x+v*dt)) < P_\text{target}$. There are two cases to
1300 consider. In the most common, no domains have unfolded since the last
1301 step, and we expect the next step to be slightly shorter than the
1302 previous one. In the less common, domains did unfold in the last
1303 step, and we expect the next step to be considerably longer than the
1306 double search_dt(transition_t *transition,
1308 environment_t *env, double x,
1309 double target_prob, double max_dt, double v)
1310 { /* Returns a valid timestep dt in seconds for the given transition.
1311 * Takes a list of all states, a pointer env to the environmental data,
1312 * a starting separation x in m, a target_prob between 0 and 1,
1313 * max_dt in s, stretching velocity v in m/s.
1315 double F, k, dtCur, dtU, dtUCur, dtL, dt;
1317 /* get upper bound using the current position */
1318 F = find_tension(state_list, env, x, 0); /* BUG. repeated calculation */
1319 //printf("Start with x = %g (F = %g)\n", x, F);
1320 k = accel_k(transition->k, F, env, transition->k_params);
1321 //printf("x %g\tF %g\tk %g\n", x, F, k);
1322 dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */
1324 //printf("overshot max_dt\n");
1327 if (v == 0) /* discrete stepping, so force is temporarily constant */
1330 /* set a lower bound on dt too */
1333 /* The dt determined above may produce illegitimate forces or ks.
1334 * Reduce the upper bound until we have valid ks. */
1336 F = find_tension(state_list, env, x+v*dt, 0);
1337 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1340 F = find_tension(state_list, env, x+v*dt, 0);
1342 //printf("Try for dt = %g (F = %g)\n", dt, F);
1343 k = accel_k(transition->k, F, env, transition->k_params);
1344 /* returned k may be -1.0 */
1345 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1346 while (k == -1.0) { /* reduce step until we hit a valid k */
1348 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1349 F = find_tension(state_list, env, x+v*dt, 0);
1350 //printf("Try for dt = %g (F = %g)\n", dt, F);
1351 k = accel_k(transition->k, F, env, transition->k_params);
1352 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1354 assert(dtU > 1e-14); /* timestep to valid k too small */
1355 dtUCur = target_prob / k; /* safe timestep back from x+dtU */
1357 return dt; /* dtU is safe. */
1359 /* dtU wasn't safe, lets see what would be. */
1360 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1361 dt = (dtU + dtL) / 2.0;
1362 F = find_tension(state_list, env, x+v*dt, 0);
1363 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1364 k = accel_k(transition->k, F, env, transition->k_params);
1365 dtCur = target_prob / k;
1366 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1367 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1369 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1370 dtU = dt; /* ... stepping out only dtCur would be safe */
1373 break; /* dtCur = dt */
1375 return MAX(dtUCur, dtL);
1380 To determine $dt$ for an array of potentially different folded domains,
1381 we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
1384 double determine_dt(list_t *state_list, list_t *transition_list,
1385 environment_t *env, double x,
1386 double target_prob, double dt_max, double v)
1387 { /* Returns the timestep dt in seconds.
1388 * Takes the state and transition lists, a pointer to the environment,
1389 * an extension x in m, target_prob between 0 and 1, dt_max in s,
1391 double dt=dt_max, new_dt;
1392 assert(target_prob > 0.0); assert(target_prob < 1.0);
1393 assert(dt_max > 0.0);
1394 assert(state_list != NULL);
1395 assert(transition_list != NULL);
1396 transition_list = head(transition_list);
1397 while (transition_list != NULL) {
1398 if (T(transition_list)->initial_state->num_domains > 0) {
1399 new_dt = search_dt(T(transition_list), state_list, env, x, target_prob, dt, v);
1400 dt = MIN(dt, new_dt);
1402 transition_list = transition_list->next;
1409 \subsection{State data}
1410 \label{sec.state_data}
1412 The domains exist in one of an array of possible states. Each state
1413 has a tension model which determines it's force/extension curve
1414 (possibly [[null]] for rigid states, see Appendix
1415 \ref{sec.null_tension}).
1417 Groups are, for example, for two domain states with WLC tensions; one
1418 with $p=0.4\U{nm}$, $L=10\U{nm}$ and the other with $p=0.4\U{nm}$,
1419 $L=28\U{nm}$. Since the persistence lengths are the same, you would
1420 like to add the contour lengths of all the domains in \emph{both}
1421 states, and plug that total contour length into a single WLC
1422 evaluation (see Section \ref{sec.find_tension}).
1423 <<state definitions>>=
1424 typedef struct state_struct {
1425 char *name; /* name string */
1426 one_dim_fn_t *tension_handler; /* tension handler */
1427 one_dim_fn_t *inverse_tension_handler; /* inverse tension handler */
1428 int tension_group; /* for grouping similar tension states */
1429 void *tension_params; /* pointer to tension parameters */
1430 destroy_data_func_t *destroy_tension;
1431 int num_domains; /* number of domains currently in state */
1432 /* cached lookups generated from the above data */
1433 int num_out_trans; /* length of out_trans array */
1434 struct transition_struct **out_trans; /* transitions leaving this state */
1435 int num_grouped_states; /* length of grouped states array */
1436 struct state_struct **grouped_states; /* states grouped with this one (doesn't include self) */
1439 /* get the state data for the current list node */
1440 #define S(list) ((state_t *)(list)->d)
1442 @ [[tension_params]] is a pointer to the parameters used by the
1443 handler function when determining the tension. [[destroy_tension]]
1444 points to a function for freeing [[tension_params]]. We it in the
1445 state structure so that [[destroy_state]] and will have an easier job
1448 [[create_]] and [[destroy_state]] are simple wrappers around
1449 [[malloc]] and [[free]].
1450 <<create/destroy state>>=
1451 state_t *create_state(char *name,
1452 one_dim_fn_t *tension_handler,
1453 one_dim_fn_t *inverse_tension_handler,
1455 void *tension_params,
1456 destroy_data_func_t *destroy_tension,
1459 state_t *ret = (state_t *)malloc(sizeof(state_t));
1460 assert(ret != NULL);
1461 /* make a copy of the name, so we aren't dependent on the original */
1462 ret->name = (char *)malloc(sizeof(char)*(strlen(name)+1));
1463 assert(ret->name != NULL);
1464 strcpy(ret->name, name); /* we know ret->name is long enough */
1466 ret->tension_handler = tension_handler;
1467 ret->inverse_tension_handler = inverse_tension_handler;
1468 ret->tension_group = tension_group;
1469 ret->tension_params = tension_params;
1470 ret->destroy_tension = destroy_tension;
1471 ret->num_domains = num_domains;
1473 ret->num_out_trans = 0;
1474 ret->out_trans = NULL;
1475 ret->num_grouped_states = 0;
1476 ret->grouped_states = NULL;
1479 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);
1484 void destroy_state(state_t *state)
1488 fprintf(stderr, "freeing state %s (%p)\n", state->name, state);
1492 if (state->destroy_tension)
1493 (*state->destroy_tension)(state->tension_params);
1494 if (state->out_trans)
1495 free(state->out_trans);
1496 if (state->grouped_states)
1497 free(state->grouped_states);
1504 We'll be storing the states in a list (see Appendix \ref{sec.lists}),
1505 so we also define a convenience function for cleaning up.
1506 <<destroy state list>>=
1507 void destroy_state_list(list_t *state_list)
1509 state_list = head(state_list);
1510 while (state_list != NULL)
1511 destroy_state((state_t *) pop(&state_list));
1516 We can also define a convenient way to get a state index from it's name.
1518 int state_index(list_t *state_list, char *name)
1521 state_list = head(state_list);
1522 while (state_list != NULL) {
1523 if (STRMATCH(S(state_list)->name, name)) return i;
1524 state_list = state_list->next;
1532 \subsection{Transition data}
1533 \label{sec.transition_data}
1535 Once you've created states (Sections \ref{sec.main},
1536 \ref{sec.model_selection}, and \ref{sec.state_data}), you need to
1537 create transitions between them.
1538 <<transition definitions>>=
1539 typedef struct transition_struct {
1540 state_t *initial_state;
1541 state_t *final_state;
1542 k_func_t *k; /* transition rate model */
1543 void *k_params; /* pointer to k parameters */
1544 destroy_data_func_t *destroy_k;
1547 /* get the transition data for the current list node */
1548 #define T(list) ((transition_t *)(list)->d)
1549 @ [[k]] is a pointer to the function determining the transition rate
1550 for a given tension. [[k_params]] and [[destroy_k]] have similar
1551 roles with regards to [[k]] as [[tension_params]] and
1552 [[destroy_tension]] have with regards to [[tension_handler]] (see
1553 Section \ref{sec.state_data}).
1555 [[create_]] and [[destroy_transition]] are simple wrappers around
1556 [[malloc]] and [[free]].
1557 <<create/destroy transition>>=
1558 transition_t *create_transition(state_t *initial_state, state_t *final_state,
1561 destroy_data_func_t *destroy_k)
1563 transition_t *ret = (transition_t *)malloc(sizeof(transition_t));
1564 assert(ret != NULL);
1565 assert(initial_state != NULL);
1566 assert(final_state != NULL);
1567 ret->initial_state = initial_state;
1568 ret->final_state = final_state;
1570 ret->k_params = k_params;
1571 ret->destroy_k = destroy_k;
1573 fprintf(stderr, "generated transition %p from %s to %s, params %p\n", ret, ret->initial_state->name, ret->final_state->name, ret->k_params);
1578 void destroy_transition(transition_t *transition)
1582 fprintf(stderr, "freeing transition %p (%s to %s)\n", transition, transition->initial_state->name, transition->final_state->name);
1584 if (transition->destroy_k)
1585 (*transition->destroy_k)(transition->k_params);
1592 We'll be storing the transitions in a list (see Appendix
1593 \ref{sec.lists}), so we also define a convenience function for
1595 <<destroy transition list>>=
1596 void destroy_transition_list(list_t *transition_list)
1598 transition_list = head(transition_list);
1599 while (transition_list != NULL)
1600 destroy_transition((transition_t *) pop(&transition_list));
1605 \subsection{Graphviz output}
1606 \label{sec.graphviz_output}
1608 <<print state graph>>=
1609 void print_state_graph(FILE *file, list_t *state_list, list_t *transition_list)
1611 state_list = head(state_list);
1612 transition_list = head(transition_list);
1613 fprintf(file, "digraph states {\n node [shape=record,color=black,fillcolor=white,style=filled];\n\n");
1615 fprintf(file, " n%d [label=\"%s\"]\n", state_index(state_list, S(state_list)->name), S(state_list)->name);
1616 if (state_list->next == NULL) break;
1617 state_list = state_list->next;
1619 fprintf(file, "\n");
1620 while (transition_list != NULL) {
1621 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));
1622 transition_list = transition_list->next;
1624 fprintf(file, "}\n");
1628 \subsection{Domain model and group handling}
1630 <<group functions>>=
1631 <<crosslink groups>>
1632 <<get active groups>>
1635 Scan through a list of states and construct an array of all the
1637 <<get active groups>>=
1638 void get_active_groups(list_t *state_list,
1639 int *pNum_active_groups,
1640 one_dim_fn_t ***pPer_group_handlers,
1641 one_dim_fn_t ***pPer_group_inverse_handlers,
1642 void ***pPer_group_data /* pointer to array of pointers to tension_handler_data_t */)
1644 void **active_groups=NULL;
1645 one_dim_fn_t *handler, *old_handler, *inverse_handler, *old_inverse_handler;
1646 one_dim_fn_t **per_group_handlers=NULL, **per_group_inverse_handlers=NULL;
1647 void **per_group_data=NULL; /* array of pointers to tension_handler_data_t */
1648 int i, j, num_domains, group, old_group, num_active_groups=0;
1649 list_t *active_states_list=NULL;
1650 tension_handler_data_t *tdata=NULL;
1653 /* get all the active groups in a list */
1654 state_list = head(state_list);
1656 fprintf(stderr, "%s:\t", __FUNCTION__);
1657 list_print(stderr, state_list, "states");
1659 while (state_list != NULL) {
1660 state = S(state_list);
1661 handler = state->tension_handler;
1662 inverse_handler = state->inverse_tension_handler;
1663 group = state->tension_group;
1664 num_domains = state->num_domains;
1665 if (list_index(active_states_list, state) == -1) {
1666 /* we haven't added this state (or it's associated group) yet */
1667 if (num_domains > 0 && handler != NULL) { /* add the state */
1668 num_active_groups++;
1669 push(&active_states_list, state, 1);
1671 fprintf(stderr, "%s:\tactive %s state (%p) with %d domain(s)\n", __FUNCTION__, state->name, state, state->num_domains);
1673 for (i=0; i < state->num_grouped_states; i++) {
1674 if (state->grouped_states[i]->num_domains > 0) {
1675 /* add this related (and active) state now */
1676 assert(state->grouped_states[i]->tension_handler == handler);
1677 assert(state->grouped_states[i]->tension_group == group);
1678 push(&active_states_list, state->grouped_states[i], 1);
1680 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);
1684 } /* else empty state or NULL handler */
1685 } /* else state already added */
1686 state_list = state_list->next;
1689 fprintf(stderr, "%s:\t", __FUNCTION__);
1690 list_print(stderr, active_states_list, "active states");
1693 assert(num_active_groups <= list_length(active_states_list));
1694 assert(num_active_groups > 0);
1696 /* allocate the arrays we'll be returning */
1697 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1698 assert(per_group_handlers != NULL);
1699 per_group_inverse_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1700 assert(per_group_inverse_handlers != NULL);
1701 per_group_data = (void **)malloc(sizeof(void *)*num_active_groups);
1702 assert(per_group_data != NULL);
1704 fprintf(stderr, "%s:\tallocating per group data %p:", __FUNCTION__, per_group_data);
1706 for (i=0; i<num_active_groups; i++) { /* we might want to recycle any group data structs allocated in the previous round */
1707 per_group_data[i] = malloc(sizeof(tension_handler_data_t));
1708 assert(per_group_data[i] != NULL);
1710 fprintf(stderr, " %d,%p", i, per_group_data[i]);
1714 fprintf(stderr, "\n");
1717 /* populate the arrays we'll be returning */
1718 active_states_list = head(active_states_list);
1720 old_inverse_handler = NULL;
1721 j = -1; /* j is the current active _group_ index */
1722 while (active_states_list != NULL) {
1723 state = (state_t *) pop(&active_states_list);
1724 handler = state->tension_handler;
1725 inverse_handler = state->inverse_tension_handler;
1726 group = state->tension_group;
1727 num_domains = state->num_domains;
1728 if (handler != old_handler || inverse_handler != old_inverse_handler || group != old_group) {
1729 j++; /* move to the next group */
1730 tdata = (tension_handler_data_t *) per_group_data[j];
1731 per_group_handlers[j] = handler;
1732 per_group_inverse_handlers[j] = inverse_handler;
1733 tdata->group_tension_params = NULL;
1735 tdata->persist = NULL; /* we might want to grab any appropriate data from the old persist here... */
1738 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);
1740 for (i=0; i<num_domains; i++) { /* add num_domains copies of the tension params */
1741 push(&tdata->group_tension_params, state->tension_params, 1);
1744 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);
1746 old_handler = handler;
1747 old_inverse_handler = inverse_handler;
1751 /* free old groups */
1752 if (*pPer_group_handlers != NULL)
1753 free(*pPer_group_handlers);
1754 if (*pPer_group_inverse_handlers != NULL)
1755 free(*pPer_group_inverse_handlers);
1756 if (*pPer_group_data != NULL) {
1757 for (i=0; i<*pNum_active_groups; i++)
1758 free((*pPer_group_data)[i]);
1759 free(*pPer_group_data);
1762 /* set new groups */
1764 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]);
1766 *pNum_active_groups = num_active_groups;
1767 *pPer_group_handlers = per_group_handlers;
1768 *pPer_group_inverse_handlers = per_group_inverse_handlers;
1769 *pPer_group_data = per_group_data;
1773 <<crosslink groups>>=
1774 void crosslink_groups(list_t *state_groups, list_t *transition_list)
1776 list_t *list, *out_trans=NULL, *associates=NULL;
1777 one_dim_fn_t *handler;
1780 state_groups = head(state_groups);
1781 while (state_groups != NULL) {
1782 /* find transitions leaving this state */
1784 fprintf(stderr, "finding transitions leaving %s\n", S(state_groups)->name);
1786 list = head(transition_list);
1787 while (list != NULL) {
1788 if (T(list)->initial_state == S(state_groups) && list_index(associates, T(list)->initial_state) != -1) {
1789 push(&out_trans, T(list), 1);
1793 n = list_length(out_trans);
1794 S(state_groups)->num_out_trans = n;
1795 S(state_groups)->out_trans = (transition_t **)malloc(sizeof(transition_t *)*n);
1796 assert(S(state_groups)->out_trans != NULL);
1797 for (i=0; i<n; i++) {
1798 S(state_groups)->out_trans[i] = (transition_t *)pop(&out_trans);
1801 /* find states grouped with this state */
1803 fprintf(stderr, "finding states grouped with %s\n", S(state_groups)->name);
1805 handler = S(state_groups)->tension_handler;
1806 group = S(state_groups)->tension_group;
1807 list = head(state_groups);
1808 while (list != NULL) {
1809 if (S(list)->tension_handler == handler && S(list)->tension_group == group && list != state_groups) {
1810 push(&associates, S(list), 1);
1814 n = list_length(associates);
1815 S(state_groups)->num_grouped_states = n;
1816 S(state_groups)->grouped_states = (state_t **)malloc(sizeof(state_t *)*n);
1817 assert(S(state_groups)->grouped_states != NULL);
1818 for (i=0; i<n; i++) {
1819 S(state_groups)->grouped_states[i] = (state_t *)pop(&associates);
1821 state_groups = state_groups->next;
1827 \section{String parsing}
1829 For handling command line arguments, we need parse delimited strings (e.g.\ [["param1,param2,param3..."]]).
1830 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1831 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1832 need to take care of parsing those parameters themselves.
1833 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1839 <<parse definitions>>
1840 <<parse declarations>>
1844 <<parse module makefile lines>>=
1845 NW_SAWSIM_MODS += parse
1846 CHECK_BINS += check_parse
1850 <<parse definitions>>=
1851 #define SEP ',' /* argument separator character */
1854 <<parse declarations>>=
1855 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1856 int *num, char ***string_array);
1857 extern void free_parsed_list(int num, char **string_array);
1860 [[parse_list_string]] allocates memory, don't forget to free it
1861 afterward with [[free_parsed_list]]. It does not alter the original.
1863 The string may start off with a [[deeper]] character
1864 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1865 [[x,y]], where the pointer is one character in on the copied string.
1866 However, when we go to free the memory, we need a pointer to the
1867 beginning of the string. In order to accommodate this for a string
1868 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1869 the first $N$ elements point to the separated arguments, and let the
1870 last element point to the start of the copied string regardless of
1872 <<parse delimited list string functions>>=
1873 /* TODO, split out into parse.hc */
1874 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1877 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1878 if (string[i] == deeper) {depth++;}
1879 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1885 void parse_list_string(char *string, char sep, char deeper, char shallower,
1886 int *num, char ***string_array)
1888 char *str=NULL, **ret=NULL;
1890 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
1892 *string_array = NULL;
1895 /* make a copy of the string, so we don't change the original */
1896 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1897 assert(str != NULL);
1898 strcpy(str, string); /* we know str is long enough */
1899 /* count the number of regions, so we can allocate pointers to them */
1902 n++; i++; /* move on to next argument */
1903 i += next_delim_index(str+i, sep, deeper, shallower);
1904 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1905 } while (str[i] != '\0');
1906 ret = (char **)malloc(sizeof(char *)*(n+1));
1907 assert(ret != NULL);
1908 /* replace the separators with '\0' & assign pointers */
1909 ret[n] = str; /* point to the front of the copied string */
1912 for(i=1; i<n; i++) {
1913 j += next_delim_index(str+j, sep, deeper, shallower);
1915 ret[i] = str+j; /* point to the separated arguments */
1917 /* strip off bounding braces, if any */
1918 for(i=0; i<n; i++) {
1919 if (ret[i][0] == deeper) {
1923 if (ret[i][strlen(ret[i])-1] == shallower)
1924 ret[i][strlen(ret[i])-1] = '\0';
1927 *string_array = ret;
1930 void free_parsed_list(int num, char **string_array)
1933 /* frees all items, since they were allocated as a single string*/
1934 free(string_array[num]);
1935 /* and free the array of pointers */
1941 \subsection{Parsing implementation}
1947 <<parse delimited list string functions>>
1951 #include <assert.h> /* assert() */
1952 #include <stdlib.h> /* NULL */
1953 #include <stdio.h> /* fprintf(), stdout *//*!!*/
1954 #include <string.h> /* strlen() */
1958 \subsection{Parsing unit tests}
1960 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1962 <<parse check includes>>
1963 <<string comparison definition and globals>>
1964 <<check relative error>>
1965 <<parse test suite>>
1966 <<main check program>>
1969 <<parse check includes>>=
1970 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
1971 #include <stdio.h> /* printf() */
1972 #include <assert.h> /* assert() */
1973 #include <string.h> /* strlen() */
1978 <<parse test suite>>=
1979 <<parse list string tests>>
1980 <<parse suite definition>>
1983 <<parse suite definition>>=
1984 Suite *test_suite (void)
1986 Suite *s = suite_create ("k model");
1987 <<parse list string test case defs>>
1989 <<parse list string test case add>>
1994 <<parse list string tests>>=
1997 START_TEST(test_next_delim_index)
1999 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
2000 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
2001 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
2002 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
2003 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
2007 START_TEST(test_parse_list_null)
2011 parse_list_string(NULL, SEP, '{', '}',
2012 &num_param_args, ¶m_args);
2013 fail_unless(num_param_args == 0, NULL);
2014 fail_unless(param_args == NULL, NULL);
2017 START_TEST(test_parse_list_single_simple)
2021 parse_list_string("arg", SEP, '{', '}',
2022 &num_param_args, ¶m_args);
2023 fail_unless(num_param_args == 1, NULL);
2024 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
2027 START_TEST(test_parse_list_single_compound)
2031 parse_list_string("{x,y,z}", SEP, '{', '}',
2032 &num_param_args, ¶m_args);
2033 fail_unless(num_param_args == 1, NULL);
2034 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
2037 START_TEST(test_parse_list_double_simple)
2041 parse_list_string("abc,def", SEP, '{', '}',
2042 &num_param_args, ¶m_args);
2043 fail_unless(num_param_args == 2, NULL);
2044 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2045 fail_unless(STRMATCH(param_args[1],"def"), NULL);
2048 START_TEST(test_parse_list_compound)
2052 parse_list_string("abc,{def,ghi}", SEP, '{', '}',
2053 &num_param_args, ¶m_args);
2054 fail_unless(num_param_args == 2, NULL);
2055 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2056 fail_unless(STRMATCH(param_args[1],"def,ghi"), NULL);
2060 <<parse list string test case defs>>=
2061 TCase *tc_parse_list_string = tcase_create("parse list string");
2063 <<parse list string test case add>>=
2064 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
2065 tcase_add_test(tc_parse_list_string, test_parse_list_null);
2066 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
2067 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
2068 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
2069 tcase_add_test(tc_parse_list_string, test_parse_list_compound);
2070 suite_add_tcase(s, tc_parse_list_string);
2074 \section{Unit tests}
2076 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2083 <<check relative error>>
2086 <<main check program>>
2098 <<determine dt tests>>
2100 <<does domain transition tests>>
2101 <<randomly unfold domains tests>>
2102 <<suite definition>>
2105 <<suite definition>>=
2106 Suite *test_suite (void)
2108 Suite *s = suite_create ("sawsim");
2109 <<F test case defs>>
2110 <<determine dt test case defs>>
2111 <<happens test case defs>>
2112 <<does domain transition test case defs>>
2113 <<randomly unfold domains test case defs>>
2116 <<determine dt test case add>>
2117 <<happens test case add>>
2118 <<does domain transition test case add>>
2119 <<randomly unfold domains test case add>>
2122 tcase_add_checked_fixture(tc_strip_address,
2123 setup_strip_address,
2124 teardown_strip_address);
2130 <<main check program>>=
2134 Suite *s = test_suite();
2135 SRunner *sr = srunner_create(s);
2136 srunner_run_all(sr, CK_ENV);
2137 number_failed = srunner_ntests_failed(sr);
2139 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
2143 \subsection{F tests}
2145 <<worm-like chain tests>>
2147 <<F test case defs>>=
2148 <<worm-like chain test case def>>
2150 <<F test case add>>=
2151 <<worm-like chain test case add>>
2154 <<worm-like chain tests>>=
2155 <<worm-like chain function>>
2157 START_TEST(test_wlc_at_zero)
2159 double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
2160 fail_unless(abs(wlc(x, T, p, L)) < lim, \
2161 "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
2162 x, T, p, L, abs(wlc(x, T, p, L)), lim);
2166 START_TEST(test_wlc_at_half)
2168 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
2169 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
2170 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
2172 * linear term = x/L = 0.5
2173 * nonlinear + linear = 0.75 + 0.5 = 1.25
2174 * wlc = 10e21*1.25 = 12.5e21
2176 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
2177 "wlc(%g, %g, %g, %g) = %g != %g",
2178 x, T, p, L, wlc(x, T, p, L), 12.5e21);
2182 <<worm-like chain test case def>>=
2183 TCase *tc_wlc = tcase_create("WLC");
2186 <<worm-like chain test case add>>=
2187 tcase_add_test(tc_wlc, test_wlc_at_zero);
2188 tcase_add_test(tc_wlc, test_wlc_at_half);
2189 suite_add_tcase(s, tc_wlc);
2192 \subsection{Model tests}
2194 Check the searching with [[linear_k]].
2195 Check overwhelming force treatment with the heavyside-esque step [[?]].
2197 Hmm, I'm not really sure what I was doing here...
2198 <<determine dt tests>>=
2199 double linear_k(double F, environment_t *env, void *params)
2201 double Fnot = *(double *)params;
2205 START_TEST(test_determine_dt_linear_k)
2208 transition_t unfold={0};
2209 environment_t env = {0};
2210 double F[]={0,1,2,3,4,5,6};
2211 double dt_max=3.0, Fnot=3.0, x=1.0, targ_p=0.1, v=1.0;
2212 list_t *state_list=NULL, *transition_list=NULL;
2215 folded.tension_handler = &hooke_handler;
2216 folded.num_domains = 1;
2217 unfold.initial_state = &folded;
2218 unfold.k = linear_k;
2219 unfold.k_params = &Fnot;
2220 push(&state_list, &folded, 1);
2221 push(&transition_list, &unfold, 1);
2222 for( i=0; i < sizeof(F)/sizeof(double); i++) {
2223 //fail_unless(determine_dt(state_list, transition_list, &env, x, targ_p, dt_max, v)==1, NULL);
2228 <<determine dt test case defs>>=
2229 TCase *tc_determine_dt = tcase_create("Determine dt");
2231 <<determine dt test case add>>=
2232 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
2233 suite_add_tcase(s, tc_determine_dt);
2238 <<happens test case defs>>=
2240 <<happens test case add>>=
2243 <<does domain transition tests>>=
2245 <<does domain transition test case defs>>=
2247 <<does domain transition test case add>>=
2250 <<randomly unfold domains tests>>=
2252 <<randomly unfold domains test case defs>>=
2254 <<randomly unfold domains test case add>>=
2258 \section{Balancing group extensions}
2259 \label{sec.tension_balance}
2261 For a given total extension $x$ (of the piezo), the various domain
2262 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
2263 amounts, and we need to tweak the portion that each extends to balance
2264 the tension amongst the active groups. Since the tension balancing is
2265 separable from the bulk of the simulation, we break it out into a
2266 separate module. The interface is defined in [[tension_balance.h]],
2267 the implementation in [[tension_balance.c]], and the unit testing in
2268 [[check_tension_balance.c]]
2270 <<tension-balance.h>>=
2271 #ifndef TENSION_BALANCE_H
2272 #define TENSION_BALANCE_H
2274 <<tension balance function declaration>>
2275 #endif /* TENSION_BALANCE_H */
2278 <<tension balance functions>>=
2279 <<one dimensional bracket>>
2280 <<one dimensional solve>>
2282 <<group stiffness function>>
2283 <<chain stiffness function>>
2284 <<tension balance function>>
2287 <<tension balance module makefile lines>>=
2288 NW_SAWSIM_MODS += tension_balance
2289 CHECK_BINS += check_tension_balance
2290 check_tension_balance_MODS =
2293 The entire force balancing problem reduces to a solving two nested
2294 one-dimensional root-finding problems. First we define the one
2295 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
2296 populated group). $x(x_0)$ is also strictly monotonically increasing,
2297 so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
2298 stores the last successful value of $x$, since we expect to be taking
2299 small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
2300 indicates that the groups have changed.
2301 <<tension balance function declaration>>=
2302 double tension_balance(int num_tension_groups,
2303 one_dim_fn_t **tension_handlers,
2304 one_dim_fn_t **inverse_tension_handlers,
2305 void **params, /* array of pointers to tension_handler_data_t */
2310 <<tension balance function>>=
2311 double tension_balance(int num_tension_groups,
2312 one_dim_fn_t **tension_handlers,
2313 one_dim_fn_t **inverse_tension_handlers,
2314 void **params, /* array of pointers to tension_handler_data_t */
2319 static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL,NULL};
2320 double F, xo_guess, xo, lb, ub;
2321 double min_relx=1e-6, min_rely=1e-6;
2322 int max_steps=200, i;
2324 assert(num_tension_groups > 0);
2325 assert(tension_handlers != NULL);
2326 assert(params != NULL);
2327 assert(num_tension_groups > 0);
2329 if (last_x == -1 || x_xo_data.xi == NULL) { /* new group setup, reset x_xo_data */
2330 /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
2331 if (x_xo_data.xi != NULL)
2333 /* construct new x_xo_data */
2334 x_xo_data.n_groups = num_tension_groups;
2335 x_xo_data.tension_handler = tension_handlers;
2336 x_xo_data.inverse_tension_handler = inverse_tension_handlers;
2337 x_xo_data.handler_data = params;
2339 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);
2341 x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
2342 assert(x_xo_data.xi != NULL);
2343 for (i=0; i<num_tension_groups; i++)
2344 x_xo_data.xi[i] = last_x;
2347 if (num_tension_groups == 1) { /* only one group, no balancing required */
2349 x_xo_data.xi[0] = xo;
2351 fprintf(stderr, "%s:\tonly one active group, no balancing required\n", __FUNCTION__);
2355 fprintf(stderr, "%s:\tbalancing for x = %g with", __FUNCTION__, x);
2356 for(i=0; i<num_tension_groups; i++) fprintf(stderr, "\tx[%d]=%g", i, x_xo_data.xi[i]);
2357 fprintf(stderr, "\n");
2359 if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
2360 if (x == 0) xo_guess = length_scale;
2361 else xo_guess = x/num_tension_groups;
2363 fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
2365 oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
2366 } else { /* work off the last balanced point */
2368 fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
2371 lb = x_xo_data.xi[0];
2372 ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */
2373 } else if (x < last_x) {
2374 lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
2375 ub = x_xo_data.xi[0];
2376 } else { /* x == last_x */
2377 fprintf(stderr, "not moving\n");
2378 lb= x_xo_data.xi[0];
2379 ub= x_xo_data.xi[0];
2383 fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
2384 __FUNCTION__, x, lb, ub);
2386 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
2387 F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
2389 if (num_tension_groups > 1) {
2390 fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
2391 for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
2392 fprintf(stderr, "\n");
2397 F = (*tension_handlers[0])(xo, params[0]);
2399 if (stiffness != NULL)
2400 *stiffness = chain_stiffness(&x_xo_data, min_relx);
2407 <<tension balance internal definitions>>=
2408 <<x of x0 definitions>>
2411 <<x of x0 definitions>>=
2412 typedef struct x_of_xo_data_struct {
2414 one_dim_fn_t **tension_handler; /* array of fn pointers */
2415 one_dim_fn_t **inverse_tension_handler; /* inverse fn pointers */
2416 void **handler_data; /* array of pointers to tension_handler_data_t */
2417 double *xi; /* array of group extensions */
2422 double x_of_xo(double xo, void *pdata)
2424 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
2425 double F, x=0, xi, xi_guess, lb, ub, slope;
2427 double min_relx=1e-6, min_rely=1e-6;
2429 assert(data->n_groups > 0);
2432 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);
2435 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
2437 fprintf(stderr, "%s:\tfinding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
2441 if (data->xi) data->xi[0] = xo;
2442 for (i=1; i < data->n_groups; i++) {
2443 if (data->inverse_tension_handler[i] != NULL) {
2444 xi = (*data->inverse_tension_handler[i])(F, data->handler_data[i]);
2445 } else { /* invert numerically */
2446 if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
2447 else xi_guess = data->xi[i];
2449 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]);
2451 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
2452 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
2457 fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
2461 fprintf(stderr, "%s:\tfound x(x0=%g) = %g\n", __FUNCTION__, xo, x);
2467 <<group stiffness function>>=
2468 double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double x, double relx)
2470 double dx, xl, F, Fl, stiffness;
2472 assert(relx > 0 && relx < 1);
2473 if (x == 0 || relx == 0) {
2474 dx = 1e-12; /* HACK, how to get length scale? */
2484 F = (*tension_handler)(x, handler_data);
2485 Fl = (*tension_handler)(xl, handler_data);
2486 stiffness = (F-Fl)/dx;
2487 assert(stiffness > 0);
2492 <<chain stiffness function>>=
2493 double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx)
2497 /* add group stiffnesses in series */
2498 for (i=0; i < x_xo_data->n_groups; i++) {
2500 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);
2503 stiffness += 1.0/group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx);
2505 stiffness = 1.0 / stiffness;
2511 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
2512 number of steps for monotonically increasing $f(x)$. Simple
2513 bisection, so it's robust and converges fairly quickly. You can set
2514 the maximum number of steps to take, but if you have enough processor
2515 power, it's better to limit your precision with the relative limits
2516 [[min_relx]] and [[y]]. These ensure that the width of the bracket is
2517 small on the length scale of it's larger side. Note that \emph{both}
2518 relative limits must be satisfied for the search to stop.
2519 <<one dimensional function definition>>=
2520 /* equivalent to gsl_func_t */
2521 typedef double one_dim_fn_t(double x, void *params);
2523 <<one dimensional solve declaration>>=
2524 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2525 double min_relx, double min_rely, int max_steps,
2529 <<one dimensional solve>>=
2530 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2531 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2532 double min_relx, double min_rely, int max_steps,
2535 double yx, ylb, yub, x;
2538 ylb = (*f)(lb, params);
2539 yub = (*f)(ub, params);
2541 /* check some simple cases */
2543 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bracketed */
2544 else return lb; /* any x on the range [lb, ub] would work */
2546 if (ylb == y) { x=lb; yx=ylb; goto end; }
2547 if (yub == y) { x=ub; yx=yub; goto end; }
2550 fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2552 assert(ylb < y); assert(yub > y);
2553 for (i=0; i<max_steps; i++) {
2555 yx = (*f)(x, params);
2557 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);
2559 if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2560 else if (yx > y) { ub=x; yub=yx; }
2561 else /* yx < y */ { lb=x; ylb=yx; }
2566 fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2572 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2573 Generate bracketing $x$ values through bisection or exponential growth.
2574 <<one dimensional bracket declaration>>=
2575 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2578 <<one dimensional bracket>>=
2579 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2581 double yx, step, x=xguess;
2584 fprintf(stderr, "%s:\tbracketing %g (params %p)\n", __FUNCTION__, y, params);
2586 yx = (*f)(x, params);
2588 fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2595 fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2599 if (x == 0) x = length_scale/2.0; /* HACK */
2602 yx = (*f)(x, params);
2604 fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2609 fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2612 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2616 \subsection{Balancing implementation}
2618 <<tension-balance.c>>=
2621 <<tension balance includes>>
2622 #include "tension_balance.h"
2624 double length_scale = 1e-10; /* HACK */
2626 <<tension balance internal definitions>>
2627 <<tension balance functions>>
2630 <<tension balance includes>>=
2631 #include <assert.h> /* assert() */
2632 #include <stdlib.h> /* NULL */
2633 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2634 #include <stdio.h> /* fprintf(), stdout */
2638 \subsection{Balancing unit tests}
2640 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2641 <<check-tension-balance.c>>=
2642 <<tension balance check includes>>
2643 <<tension balance test suite>>
2644 <<main check program>>
2647 <<tension balance check includes>>=
2648 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2649 #include <assert.h> /* assert() */
2652 #include "tension_balance.h"
2655 <<tension balance test suite>>=
2656 <<tension balance function tests>>
2657 <<tension balance suite definition>>
2660 <<tension balance suite definition>>=
2661 Suite *test_suite (void)
2663 Suite *s = suite_create ("tension balance");
2664 <<tension balance function test case defs>>
2666 <<tension balance function test case adds>>
2671 <<tension balance function tests>>=
2672 <<check relative error>>
2674 double hooke(double x, void *pK)
2677 return *((double*)pK) * x;
2680 START_TEST(test_single_function)
2682 double k=5, x=3, last_x=2.0, F;
2683 one_dim_fn_t *handlers[] = {&hooke};
2684 one_dim_fn_t *inverse_handlers[] = {NULL};
2685 void *data[] = {&k};
2686 F = tension_balance(1, handlers, inverse_handlers, data, last_x, x, NULL);
2687 fail_unless(F = k*x, NULL);
2692 We can also test balancing two springs with different spring constants.
2693 The analytic solution for a general number of springs is given in Appendix \ref{sec.math_hooke}.
2694 <<tension balance function tests>>=
2695 START_TEST(test_double_hooke)
2697 double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
2698 one_dim_fn_t *handlers[] = {&hooke, &hooke};
2699 one_dim_fn_t *inverse_handlers[] = {NULL, NULL};
2700 void *data[] = {&k1, &k2};
2701 F = tension_balance(2, handlers, inverse_handlers, data, last_x, x, NULL);
2705 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2706 //CHECK_ERR(1e-6, x1e, xi[0]);
2707 //CHECK_ERR(1e-6, x2e, xi[1]);
2708 CHECK_ERR(1e-6, Fe, F);
2712 TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
2714 <<tension balance function test case defs>>=
2715 TCase *tc_tbfunc = tcase_create("tension balance function");
2718 <<tension balance function test case adds>>=
2719 tcase_add_test(tc_tbfunc, test_single_function);
2720 tcase_add_test(tc_tbfunc, test_double_hooke);
2721 suite_add_tcase(s, tc_tbfunc);
2727 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2728 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2729 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2735 <<list definitions>>
2736 <<list declarations>>
2740 <<list declarations>>=
2741 <<head and tail declarations>>
2742 <<list length declaration>>
2743 <<push declaration>>
2745 <<list destroy declaration>>
2746 <<list to array declaration>>
2747 <<move declaration>>
2748 <<sort declaration>>
2749 <<list index declaration>>
2750 <<list element declaration>>
2751 <<unique declaration>>
2752 <<list print declaration>>
2756 <<create and destroy node>>
2771 <<list module makefile lines>>=
2772 NW_SAWSIM_MODS += list
2773 CHECK_BINS += check_list
2777 <<list definitions>>=
2778 typedef struct list_struct {
2779 struct list_struct *next;
2780 struct list_struct *prev;
2784 <<comparison function typedef>>
2787 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2788 <<head and tail declarations>>=
2789 list_t *head(list_t *list);
2790 list_t *tail(list_t *list);
2793 list_t *head(list_t *list)
2797 while (list->prev != NULL) {
2803 list_t *tail(list_t *list)
2807 while (list->next != NULL) {
2814 <<list length declaration>>=
2815 int list_length(list_t *list);
2818 int list_length(list_t *list)
2825 while (list->next != NULL) {
2833 [[push]] inserts a node relative to the current position in the list
2834 without changing the current position.
2835 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2836 so we set it to that of the pushed domain.
2837 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2838 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2839 <<push declaration>>=
2840 void push(list_t **pList, void *data, int below);
2843 void push(list_t **pList, void *data, int below)
2845 list_t *list, *new_node;
2846 assert(pList != NULL);
2848 new_node = create_node(data);
2851 else if (below==0) { /* insert above */
2852 if (list->prev != NULL)
2853 list->prev->next = new_node;
2854 new_node->prev = list->prev;
2855 list->prev = new_node;
2856 new_node->next = list;
2857 } else { /* insert below */
2858 if (list->next != NULL)
2859 list->next->prev = new_node;
2860 new_node->next = list->next;
2861 list->next = new_node;
2862 new_node->prev = list;
2867 [[pop]] removes the current domain node, moving the current position
2868 to the node after it, unless that node is [[NULL]], in which case move
2869 the current position to the node before it.
2870 <<pop declaration>>=
2871 void *pop(list_t **pList);
2874 void *pop(list_t **pList)
2876 list_t *list, *popped;
2878 assert(pList != NULL);
2880 assert(list != NULL); /* not an empty list */
2882 /* bypass the popped node */
2883 if (list->prev != NULL)
2884 list->prev->next = list->next;
2885 if (list->next != NULL) {
2886 list->next->prev = list->prev;
2887 *pList = list->next;
2889 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
2891 destroy_node(popped);
2896 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
2897 If the cleanup function is [[NULL]], no cleanup function is called.
2898 The nodes are not popped in any particular order.
2899 <<list destroy declaration>>=
2900 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
2903 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
2907 assert(pList != NULL);
2910 assert(list != NULL); /* not an empty list */
2911 while (list != NULL) {
2913 if (destroy != NULL)
2919 [[list_to_array]] converts a list to an array of pointers, preserving order.
2920 <<list to array declaration>>=
2921 void list_to_array(list_t *list, int *length, void ***array);
2924 void list_to_array(list_t *list, int *pLength, void ***pArray)
2928 assert(list != NULL);
2929 assert(pLength != NULL);
2930 assert(pArray != NULL);
2932 length = list_length(list);
2933 array = (void **)malloc(sizeof(void *)*length);
2934 assert(array != NULL);
2937 while (list != NULL) {
2938 array[i++] = list->d;
2946 [[move]] moves the current element along the list in either direction.
2947 <<move declaration>>=
2948 void move(list_t **pList, int downstream);
2951 void move(list_t **pList, int downstream)
2953 list_t *list, *flipper;
2955 assert(pList != NULL);
2956 list = *pList; /* a>B>c>d */
2957 assert(list != NULL); /* not an empty list */
2958 if (downstream == 0)
2959 flipper = list->prev; /* flipper is A */
2961 flipper = list->next; /* flipper is C */
2962 /* ensure there is somewhere to go in the direction we're heading */
2963 assert(flipper != NULL);
2964 /* remove the current element */
2965 data = pop(&list); /* data=B, a>C>d */
2966 /* and put it back in in the correct spot */
2968 if (downstream == 0) {
2969 push(&list, data, 0); /* b>A>c>d */
2970 list = list->prev; /* B>a>c>d */
2972 push(&list, data, 1); /* a>C>b>d */
2973 list = list->next; /* a>c>B>d */
2979 [[sort]] sorts a list from smallest to largest according to a user
2980 specified comparison.
2981 <<comparison function typedef>>=
2982 typedef int (comparison_fn_t)(void *A, void *B);
2985 <<int comparison function>>=
2986 int int_comparison(void *A, void *B)
2991 if (a > b) return 1;
2992 else if (a == b) return 0;
2997 <<sort declaration>>=
2998 void sort(list_t **pList, comparison_fn_t *comp);
3000 Here we use the bubble sort algorithm.
3002 void sort(list_t **pList, comparison_fn_t *comp)
3006 assert(pList != NULL);
3008 assert(list != NULL); /* not an empty list */
3009 while (stable == 0) {
3012 while (list->next != NULL) {
3013 if ((*comp)(list->d, list->next->d) > 0) {
3015 move(&list, 1 /* downstream */);
3025 [[list_index]] finds the location of [[data]] in [[list]]. Returns
3026 [[-1]] if [[data]] is not in [[list]].
3027 <<list index declaration>>=
3028 int list_index(list_t *list, void *data);
3032 int list_index(list_t *list, void *data)
3036 while (list != NULL) {
3037 if (list->d == data) return i;
3046 [[list_element]] returns the data bound to the $i^\text{th}$ node of the list.
3047 <<list element declaration>>=
3048 void *list_element(list_t *list, int i);
3051 void *list_element(list_t *list, int i)
3055 while (list != NULL) {
3056 if (i == j) return list->d;
3066 [[unique]] removes duplicates from a sorted list according to a user
3067 specified comparison. Don't do this unless you have other pointers
3068 any dynamically allocated data in your list, because [[unique]] just
3069 drops any non-unique data without freeing it.
3070 <<unique declaration>>=
3071 void unique(list_t **pList, comparison_fn_t *comp);
3074 void unique(list_t **pList, comparison_fn_t *comp)
3077 assert(pList != NULL);
3079 assert(list != NULL); /* not an empty list */
3081 while (list->next != NULL) {
3082 if ((*comp)(list->d, list->next->d) == 0) {
3091 [[list_print]] prints a list to a [[FILE *]] stream.
3092 <<list print declaration>>=
3093 void list_print(FILE *file, list_t *list, char *list_name);
3096 void list_print(FILE *file, list_t *list, char *list_name)
3098 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
3100 while (list != NULL) {
3101 fprintf(file, " %p", list->d);
3104 fprintf(file, "\n");
3108 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
3109 <<create and destroy node>>=
3110 list_t *create_node(void *data)
3112 list_t *ret = (list_t *)malloc(sizeof(list_t));
3113 assert(ret != NULL);
3120 void destroy_node(list_t *node)
3126 The user must free the data pointed to by the node on their own.
3128 \subsection{List implementation}
3138 #include <assert.h> /* assert() */
3139 #include <stdlib.h> /* malloc(), free(), rand() */
3140 #include <stdio.h> /* fprintf(), stdout, FILE */
3141 #include "global.h" /* destroy_data_func_t */
3144 \subsection{List unit tests}
3146 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3148 <<list check includes>>
3150 <<main check program>>
3153 <<list check includes>>=
3154 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3155 #include <stdio.h> /* FILE */
3161 <<list test suite>>=
3164 <<list suite definition>>
3167 <<list suite definition>>=
3168 Suite *test_suite (void)
3170 Suite *s = suite_create ("list");
3171 <<push test case defs>>
3172 <<pop test case defs>>
3174 <<push test case adds>>
3175 <<pop test case adds>>
3181 START_TEST(test_push)
3186 push(&list, (void *)i, 1);
3187 fail_unless(list == head(list), NULL);
3188 fail_unless( (int)list->d == 0 );
3189 for (i=0; i<n; i++) {
3190 /* we expect 0, n-1, n-2, ..., 1 */
3193 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
3199 <<push test case defs>>=
3200 TCase *tc_push = tcase_create("push");
3203 <<push test case adds>>=
3204 tcase_add_test(tc_push, test_push);
3205 suite_add_tcase(s, tc_push);
3210 <<pop test case defs>>=
3212 <<pop test case adds>>=
3215 \section{Function string evaluation}
3217 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).
3218 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
3219 The solution is to run a scripting language as a subprocess accessed via pipes.
3220 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
3222 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
3223 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
3224 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.
3225 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
3227 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.
3228 Then you can either statically or dynamically link to those hardcoded functions.
3229 While much less flexible, this approach would be a fairly simple-to-implement fallback.
3231 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
3232 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
3235 #ifndef STRING_EVAL_H
3236 #define STRING_EVAL_H
3238 <<string eval setup declaration>>
3239 <<string eval function declaration>>
3240 <<string eval teardown declaration>>
3241 #endif /* STRING_EVAL_H */
3244 <<string eval module makefile lines>>=
3245 NW_SAWSIM_MODS += string_eval
3246 CHECK_BINS += check_string_eval
3247 check_string_eval_MODS =
3250 For an introduction to POSIX process control, see\\
3251 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
3252 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
3253 and of course, the relavent [[man]] pages.
3255 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
3256 [[execvp]] replaces the calling process' program with a new program.
3257 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
3258 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
3259 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
3260 The new program needs command line arguments, just like it would if you were running it from a shell.
3261 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
3262 with the final array entry being a [[NULL]] pointer.
3264 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
3265 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
3266 <<string eval subprocess definitions>>=
3267 #define SUBPROCESS "python"
3268 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
3269 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
3270 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
3272 The [[i]] option lets Python know that it should run in interactive mode.
3273 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
3274 In interactive mode, python acts on every instruction as soon as it is recieved.
3275 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.
3276 %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.
3278 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
3279 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
3280 [[fork]] returns two copies of the same program, executing the original code.
3281 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
3282 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
3284 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.
3285 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
3286 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
3287 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
3288 <<string eval pipe definitions>>=
3289 #define PIPE_READ 0 /* the end you read from */
3290 #define PIPE_WRITE 1 /* the end you write to */
3292 #define STDIN 0 /* initial index of stdin pair */
3293 #define STDOUT 2 /* initial index of stdout pair */
3296 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
3298 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.
3300 <<string eval setup declaration>>=
3301 extern void string_eval_setup(FILE **pIn, FILE **pOut);
3303 <<string eval setup definition>>=
3304 void string_eval_setup(FILE **pIn, FILE **pOut)
3307 int pfd[4]; /* file descriptors for stdin and stdout */
3309 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
3310 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
3312 pid = fork(); /* split process into two copies */
3314 if (pid == -1) { /* fork error */
3315 perror("fork error");
3317 } else if (pid == 0) { /* child */
3318 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
3319 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
3320 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
3321 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
3322 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
3323 perror("exec error"); /* exec shouldn't return */
3325 } else { /* parent */
3326 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
3327 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
3328 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
3329 if ( *pIn == NULL ) {
3330 perror("fdopen (in)");
3333 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
3334 if ( *pOut == NULL ) {
3335 perror("fdopen (out)");
3342 To use the evaluating subprocess, we just pipe in our command, and read out the results.
3343 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
3344 <<string eval function declaration>>=
3345 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
3347 <<string eval function definition>>=
3348 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
3351 rc = fprintf(in, "%s", input);
3352 assert(rc == strlen(input));
3355 alarm(1); /* set a one second timeout on the read */
3356 assert( fgets(output, buflen, out) != NULL );
3357 alarm(0); /* cancel the timeout */
3358 //fprintf(stderr, "eval: %s ----> %s", input, output);
3361 The [[alarm]] calls set and clear a timeout on the returned output.
3362 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.
3363 This protects against invalid input for which a line of output is not printed to [[stdout]].
3364 Other invalid inputs (e.g.\ those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
3365 If you are getting strange results, check your python code seperately. TODO, better error handling.
3367 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
3368 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
3369 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
3370 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
3371 <<string eval teardown declaration>>=
3372 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
3375 <<string eval teardown definition>>=
3376 void string_eval_teardown(FILE **pIn, FILE **pOut)
3381 /* redirect Python's stderr */
3382 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
3386 assert( fclose(*pIn) == 0 );
3388 assert( fclose(*pOut) == 0 );
3391 /* wait for python to exit */
3393 pid = wait(&stat_loc);
3400 if (WIFEXITED(stat_loc)) {
3401 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
3402 } else if (WIFSIGNALED(stat_loc)) {
3403 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
3408 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
3410 \subsection{String evaluation implementation}
3414 <<string eval includes>>
3415 #include "string_eval.h"
3416 <<string eval internal definitions>>
3417 <<string eval functions>>
3420 <<string eval includes>>=
3421 #include <assert.h> /* assert() */
3422 #include <stdlib.h> /* NULL */
3423 #include <stdio.h> /* fprintf(), stdout, fdopen() */
3424 #include <string.h> /* strlen() */
3425 #include <sys/types.h> /* pid_t */
3426 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
3427 #include <wait.h> /* wait() */
3430 <<string eval internal definitions>>=
3431 <<string eval subprocess definitions>>
3432 <<string eval pipe definitions>>
3435 <<string eval functions>>=
3436 <<string eval setup definition>>
3437 <<string eval function definition>>
3438 <<string eval teardown definition>>
3441 \subsection{String evaluation unit tests}
3443 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3444 <<check-string-eval.c>>=
3445 <<string eval check includes>>
3446 <<string comparison definition and globals>>
3447 <<string eval test suite>>
3448 <<main check program>>
3451 <<string eval check includes>>=
3452 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3453 #include <stdio.h> /* printf() */
3454 #include <assert.h> /* assert() */
3455 #include <string.h> /* strlen() */
3456 #include <signal.h> /* SIGKILL */
3458 #include "string_eval.h"
3461 <<string eval test suite>>=
3462 <<string eval tests>>
3463 <<string eval suite definition>>
3466 <<string eval suite definition>>=
3467 Suite *test_suite (void)
3469 Suite *s = suite_create ("string eval");
3470 <<string eval test case defs>>
3472 <<string eval test case add>>
3477 <<string eval tests>>=
3478 START_TEST(test_setup_teardown)
3481 string_eval_setup(&in, &out);
3482 string_eval_teardown(&in, &out);
3485 START_TEST(test_invalid_command)
3488 char input[80], output[80]={};
3489 string_eval_setup(&in, &out);
3490 sprintf(input, "print ABCDefg\n");
3491 string_eval(in, out, input, 80, output);
3492 string_eval_teardown(&in, &out);
3495 START_TEST(test_simple_eval)
3498 char input[80], output[80]={};
3499 string_eval_setup(&in, &out);
3500 sprintf(input, "print 3+4*5\n");
3501 string_eval(in, out, input, 80, output);
3502 fail_unless(STRMATCH(output,"23\n"), NULL);
3503 string_eval_teardown(&in, &out);
3506 START_TEST(test_multiple_evals)
3509 char input[80], output[80]={};
3510 string_eval_setup(&in, &out);
3511 sprintf(input, "print 3+4*5\n");
3512 string_eval(in, out, input, 80, output);
3513 fail_unless(STRMATCH(output,"23\n"), NULL);
3514 sprintf(input, "print (3**2 + 4**2)**0.5\n");
3515 string_eval(in, out, input, 80, output);
3516 fail_unless(STRMATCH(output,"5.0\n"), NULL);
3517 string_eval_teardown(&in, &out);
3520 START_TEST(test_eval_with_state)
3523 char input[80], output[80]={};
3524 string_eval_setup(&in, &out);
3525 sprintf(input, "print 3+4*5\n");
3526 fprintf(in, "A = 3\n");
3527 sprintf(input, "print A*3\n");
3528 string_eval(in, out, input, 80, output);
3529 fail_unless(STRMATCH(output,"9\n"), NULL);
3530 string_eval_teardown(&in, &out);
3534 <<string eval test case defs>>=
3535 TCase *tc_string_eval = tcase_create("string_eval");
3537 <<string eval test case add>>=
3538 tcase_add_test(tc_string_eval, test_setup_teardown);
3539 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
3540 tcase_add_test(tc_string_eval, test_simple_eval);
3541 tcase_add_test(tc_string_eval, test_multiple_evals);
3542 tcase_add_test(tc_string_eval, test_eval_with_state);
3543 suite_add_tcase(s, tc_string_eval);
3547 \section{Accelerating function evaluation}
3549 My first version-0.3 code was running very slowly.
3550 With the options suggested in the help
3551 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
3552 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
3553 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
3554 That's only 15 calls per solution, so the search algorithm seems reasonable.
3555 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
3560 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
3562 #endif /* ACCEL_K_H */
3565 <<accel k module makefile lines>>=
3566 NW_SAWSIM_MODS += accel_k
3567 #CHECK_BINS += check_accel_k
3568 check_accel_k_MODS =
3572 #include <assert.h> /* assert() */
3573 #include <stdlib.h> /* realloc(), free(), NULL */
3574 #include "global.h" /* environment_t */
3575 #include "k_model.h" /* k_func_t */
3576 #include "interp.h" /* interp_* */
3577 #include "accel_k.h"
3579 typedef struct accel_k_struct {
3580 interp_table_t *itable;
3586 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3587 static int num_accels = 0;
3588 static accel_k_t *accels=NULL;
3590 /* Wrap k in the standard f(x) acceleration form */
3591 static double k_wrap(double F, void *params)
3593 accel_k_t *p = (accel_k_t *)params;
3594 return (*p->k)(F, p->env, p->params);
3597 static int k_tol(double FA, double kA, double FB, double kB)
3600 if (FB-FA > 1e-12) {
3601 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3602 return 1; /* unnacceptable */
3604 //printf("acceptable tol\n");
3605 return 0; /* acceptable */
3609 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3612 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3613 assert(accels != NULL);
3614 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3616 accels[i].env = env;
3617 accels[i].params = params;
3624 for (i=0; i<num_accels; i++)
3625 interp_table_free(accels[i].itable);
3627 if (accels) free(accels);
3631 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3634 for (i=0; i<num_accels; i++) {
3635 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3636 /* We've already setup this function.
3637 * WARNING: we're assuming param and env are the same. */
3638 return interp_table_eval(accels[i].itable, accels+i, F);
3641 /* set up a new acceleration environment */
3642 i = add_accel_k(k, env, params);
3643 return interp_table_eval(accels[i].itable, accels+i, F);
3647 \section{Tension models}
3648 \label{sec.tension_models}
3650 TODO: tension model intro.
3651 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.
3653 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3654 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]].
3656 <<tension-model.h>>=
3657 #ifndef TENSION_MODEL_H
3658 #define TENSION_MODEL_H
3660 <<tension handler types>>
3661 <<constant tension model declarations>>
3662 <<hooke tension model declarations>>
3663 <<worm-like chain tension model declarations>>
3664 <<freely-jointed chain tension model declarations>>
3665 <<find tension definitions>>
3666 <<tension model global declarations>>
3667 #endif /* TENSION_MODEL_H */
3670 <<tension model module makefile lines>>=
3671 NW_SAWSIM_MODS += tension_model
3672 #CHECK_BINS += check_tension_model
3673 check_tension_model_MODS =
3675 <<tension model utils makefile lines>>=
3676 TENSION_MODEL_MODS = tension_model parse list tension_balance
3677 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3678 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3679 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3680 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
3681 notangle -Rtension-model-utils.c $< > $@
3682 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
3683 gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3684 clean_tension_model_utils :
3685 rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3686 @ The pipe symbol [[|]] marks the prerequisits following it (in this
3687 case, the directories) as ‘order-only’ prerequisites. The timestamp
3688 on these prerequisits does not effect whether the rules are executed.
3689 This is appropriate for directories, where we don't need to recompile
3690 after an unrelated has been added to the directory, but only when the
3691 source prerequisites change. See the [[make]] documentation for more
3693 (\url{http://www.gnu.org/software/make/manual/html_node/Prerequisite-Types.html}).
3696 \label{sec.null_tension}
3698 For unstretchable domains.
3700 <<null tension model getopt>>=
3701 {NULL, NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3704 \subsection{Constant}
3705 \label{sec.const_tension}
3707 An infinitely stretchable domain providing a constant tension.
3708 Obviously this cannot be inverted, so you can't put this domain in
3709 series with any other domains. We include it mostly for testing
3712 <<constant tension functions>>=
3713 <<constant tension function>>
3714 <<constant tension parameter create/destroy functions>>
3717 <<constant tension model declarations>>=
3718 <<constant tension function declaration>>
3719 <<constant tension parameter create/destroy function declarations>>
3720 <<constant tension model global declarations>>
3724 <<constant tension function declaration>>=
3725 extern double const_tension_handler(double x, void *pdata);
3727 <<constant tension function>>=
3728 double const_tension_handler(double x, void *pdata)
3730 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3731 list_t *list = data->group_tension_params;
3736 assert(list != NULL); /* empty active group?! */
3737 F = ((const_tension_param_t *)list->d)->F;
3739 fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3741 while (list != NULL) {
3742 assert(((const_tension_param_t *)list->d)->F == F);
3750 <<constant tension structure>>=
3751 typedef struct const_tension_param_struct {
3752 double F; /* tension (force) in N */
3753 } const_tension_param_t;
3757 <<constant tension parameter create/destroy function declarations>>=
3758 extern void *string_create_const_tension_param_t(char **param_strings);
3759 extern void destroy_const_tension_param_t(void *p);
3762 <<constant tension parameter create/destroy functions>>=
3763 const_tension_param_t *create_const_tension_param_t(double F)
3765 const_tension_param_t *ret
3766 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3767 assert(ret != NULL);
3772 void *string_create_const_tension_param_t(char **param_strings)
3774 return create_const_tension_param_t(
3775 safe_strtod(param_strings[0],__FUNCTION__));
3778 void destroy_const_tension_param_t(void *p)
3786 <<constant tension model global declarations>>=
3787 extern int num_const_tension_params;
3788 extern char *const_tension_param_descriptions[];
3789 extern char const_tension_param_string[];
3792 <<constant tension model globals>>=
3793 int num_const_tension_params = 1;
3794 char *const_tension_param_descriptions[] = {"tension F, N"};
3795 char const_tension_param_string[] = "0";
3798 <<constant tension model getopt>>=
3799 {&const_tension_handler, NULL, "const", "an infinitely stretchable domain with constant tension", 1, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t}
3805 <<hooke functions>>=
3806 <<internal hooke functions>>
3808 <<inverse hooke handler>>
3809 <<hooke parameter create/destroy functions>>
3812 <<hooke tension model declarations>>=
3813 <<hooke tension function declaration>>
3814 <<hooke tension parameter create/destroy function declarations>>
3815 <<hooke tension model global declarations>>
3819 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
3820 The behavior of a series of springs $k_i$ in series is given by
3822 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
3824 For a simple proof, see Appendix \ref{sec.math_hooke}.
3826 <<hooke tension function declaration>>=
3827 extern double hooke_handler(double x, void *pdata);
3828 extern double inverse_hooke_handler(double F, void *pdata);
3831 First we define a function that computes the effective spring constant
3832 (stored in a single [[hooke_param_t]]) for the entire group.
3833 <<internal hooke functions>>=
3834 static void hooke_param_grouper(tension_handler_data_t *data,
3835 hooke_param_t *param) {
3836 list_t *list = data->group_tension_params;
3840 assert(list != NULL); /* empty active group?! */
3841 while (list != NULL) {
3842 assert( ((hooke_param_t *)list->d)->k > 0 );
3843 k += 1.0/ ((hooke_param_t *)list->d)->k;
3845 fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
3846 __FUNCTION__, ((hooke_param_t *)list->d)->k);
3852 fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N (from data %p)\n",
3853 __FUNCTION__, k, x, k*x, data);
3860 Everyone knows Hooke's law: $F=kx$.
3862 double hooke_handler(double x, void *pdata)
3864 hooke_param_t param = {0};
3867 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
3873 Inverting Hooke's law gives $x=F/k$.
3874 <<inverse hooke handler>>=
3875 double inverse_hooke_handler(double F, void *pdata)
3877 hooke_param_t param = {0};
3880 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
3886 Not much to keep track of for a single effective spring, just the
3887 spring constant $k$.
3888 <<hooke structure>>=
3889 typedef struct hooke_param_struct {
3890 double k; /* spring constant in N/m */
3895 We wrap up our Hooke implementation with some book-keeping functions.
3896 <<hooke tension parameter create/destroy function declarations>>=
3897 extern void *string_create_hooke_param_t(char **param_strings);
3898 extern void destroy_hooke_param_t(void *p);
3902 <<hooke parameter create/destroy functions>>=
3903 hooke_param_t *create_hooke_param_t(double k)
3905 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
3906 assert(ret != NULL);
3911 void *string_create_hooke_param_t(char **param_strings)
3913 return create_hooke_param_t(safe_strtod(param_strings[0], __FUNCTION__));
3916 void destroy_hooke_param_t(void *p)
3923 <<hooke tension model global declarations>>=
3924 extern int num_hooke_params;
3925 extern char *hooke_param_descriptions[];
3926 extern char hooke_param_string[];
3929 <<hooke tension model globals>>=
3930 int num_hooke_params = 1;
3931 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
3932 char hooke_param_string[]="0.05";
3935 <<hooke tension model getopt>>=
3936 {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}
3939 \subsection{Worm-like chain}
3942 We can model several unfolded domains with a single worm-like chain.
3943 <<worm-like chain functions>>=
3944 <<internal worm-like chain functions>>
3945 <<worm-like chain handler>>
3946 <<inverse worm-like chain handler>>
3947 <<worm-like chain create/destroy functions>>
3950 <<worm-like chain tension model declarations>>=
3952 <<worm-like chain handler declaration>>
3953 <<inverse worm-like chain handler declaration>>
3954 <<worm-like chain create/destroy function declarations>>
3955 <<worm-like chain tension model global declarations>>
3959 <<internal worm-like chain functions>>=
3960 <<worm-like chain function>>
3961 <<inverse worm-like chain function>>
3962 <<worm-like chain parameter grouper>>
3965 The combination of all unfolded domains is modeled as a worm like chain,
3966 with the force $F_{\text{WLC}}$ approximately given by
3968 F_{\text{WLC}} \approx \frac{k_B T}{p}
3969 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
3971 where $x$ is the distance between the fixed ends,
3972 $k_B$ is Boltzmann's constant,
3973 $T$ is the absolute temperature,
3974 $p$ is the persistence length, and
3975 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
3976 <<worm-like chain function>>=
3977 static double wlc(double x, double T, double p, double L)
3979 double xL=0.0; /* uses global k_B */
3981 assert(T > 0); assert(p > 0); assert(L > 0);
3982 if (x >= L) return HUGE_VAL;
3983 xL = x/L; /* unitless */
3984 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
3985 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
3986 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
3991 Inverting the approximate WLC is fairly straightforward, if a bit tedious.
3992 Letting $x_L \equiv x/L$ and $F_T \equiv Fp/k_B T$ we have
3994 F_T &\approx \frac{1}{4} \p({\frac{1}{(1-x_L)^2} - 1}) + x_L \\
3995 F_T (1-x_L)^2 &\approx \frac{1}{4} - \frac{1}{4} (1-x_L)^2 + x_L(1-x_L)^2 \\
3996 0 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) (1-x_L)^2 + x_L(1-x_L)^2 \\
3997 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) \p({1-2x_L+x_L^2})
3998 + x_L - 2x_L^2 + x_L^3 \\
3999 &\approx \p[{\frac{1}{4} - \p({F_T + \frac{1}{4}})}]
4000 + x_L \p[{2F_T + \frac{1}{2} + 1}]
4001 + x_L^2 \p[{-\p({F_T + \frac{1}{4}}) - 2}]
4004 + x_L \p({2F_T + \frac{3}{2}})
4005 - x_L^2 \p({F_T + \frac{9}{4}})
4007 &\approx ax_L^3 + bx_L^2 + cx_L + d
4009 where $a \equiv 1$, $b \equiv -F_T - 9/4$, $c \equiv 2F_T + 3/2$, and
4011 % From \citet{wikipedia_cubic_function} the discriminant
4013 % \Delta &= -4b^3d + b^2c^2 - 4ac^3 + 18abcd - 27a^2d^2 \\
4014 % &= -4F_T\p({F_T + \frac{9}{4}})^3
4015 % + \p({F_T + \frac{9}{4}})^2 \cdot \p({2F_T + \frac{3}{2}})^2
4016 % - 4 \p({2F_T + \frac{3}{2}})^3
4017 % + 18F_T \p({F_T + \frac{9}{4}}) \cdot \p({2F_T + \frac{3}{2}})
4019 % &= -4F_T \p({F_T^3 + \frac{27}{4}F_T^2 +\frac{243}{16}F_T +\frac{729}{64}})
4020 %% a^3 + 3a^2b + 3ab^2 + b^3
4021 % + \p({F_T^2 + \frac{9}{2}F_T + \frac{81}{16}})
4022 % \cdot \p({4F_T^2 + 6F_T + \frac{9}{4}}) \\
4023 % &\qquad - 4 \p({8F_T^3 + 18F_T^2 + \frac{27}{2}F_T + \frac{27}{8}})
4024 %% a^3 + 3a^2b + 3ab^2 + b^3
4025 % + 18F_T \p({2F_T^2 + 6F_T + \frac{27}{8}})
4027 % &= -4F_T^4 - 27F_T^3 - \frac{243}{4}F_T^2 - \frac{729}{16}F_T
4028 % + 4F_T^4 + 24F_T^3 +\frac{99}{2}F_T^2 +\frac{81}{2}F_T +\frac{729}{64}\\
4029 % &\qquad - 32F_T^3 - 72F_T^2 - 54F_T - \frac{27}{2}
4030 % + 36F_T^3 + 108F_T^2 + \frac{243}{4}F_T
4032 % &= F_T^4 \p({-4 + 4}) + F_T^3 \p({-27 + 24 - 32 + 36})
4033 % + F_T^2 \p({-\frac{243}{4} + \frac{99}{2} - 72 + 108 -27}) \\
4034 % &\qquad + F_T \p({-\frac{729}{16} + \frac{81}{2} - 54 + \frac{243}{4}})
4035 % \p({\frac{729}{64} - \frac{27}{2}}) \\
4036 % &= F_T^3 - \frac{9}{4}F_T^2 + \frac{27}{16}F_T - \frac{135}{64}
4038 We can plug these coefficients into the GSL cubic solver to invert the
4039 WLC\citep{galassi05}. The applicable root is always the one which
4040 comes before the singularity, which will be the smallest real root.
4041 <<inverse worm-like chain function>>=
4042 static double inverse_wlc(double F, double T, double p, double L)
4044 double FT = F*p/(k_B*T); /* uses global k_B */
4045 double xL0, xL1, xL2;
4048 assert(T > 0); assert(p > 0); assert(L > 0);
4051 num_roots = gsl_poly_solve_cubic(-(FT+2.25),2*FT+1.5,-FT, &xL0, &xL1, &xL2);
4052 assert(num_roots > 0); assert(xL0 >= -DOUBLE_PRECISION); assert(xL0 < 1);
4053 if (xL0 < 0) xL0 = 0.0;
4059 First we define a function that computes the effective WLC parameters
4060 (stored in a single [[wlc_param_t]]) for the entire group.
4061 <<worm-like chain parameter grouper>>=
4062 static void wlc_param_grouper(tension_handler_data_t *data,
4063 wlc_param_t *param) {
4064 list_t *list = data->group_tension_params;
4065 double p=0.0, L=0.0;
4068 assert(list != NULL); /* empty active group?! */
4069 p = ((wlc_param_t *)list->d)->p;
4070 while (list != NULL) {
4071 /* enforce consistent p */
4072 assert( ((wlc_param_t *)list->d)->p == p);
4073 L += ((wlc_param_t *)list->d)->L;
4075 fprintf(stderr, "%s: WLC section %g m long in series\n",
4076 __FUNCTION__, ((wlc_param_t *)list->d)->L);
4081 fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
4082 __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
4090 <<worm-like chain handler declaration>>=
4091 extern double wlc_handler(double x, void *pdata);
4094 This model requires that each unfolded domain in the group have the
4095 same persistence length.
4096 <<worm-like chain handler>>=
4097 double wlc_handler(double x, void *pdata)
4099 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4100 wlc_param_t param = {0};
4102 wlc_param_grouper(data, ¶m);
4103 return wlc(x, data->env->T, param.p, param.L);
4108 <<inverse worm-like chain handler declaration>>=
4109 extern double inverse_wlc_handler(double F, void *pdata);
4112 <<inverse worm-like chain handler>>=
4113 double inverse_wlc_handler(double F, void *pdata)
4115 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4116 wlc_param_t param = {0};
4118 wlc_param_grouper(data, ¶m);
4119 return inverse_wlc(F, data->env->T, param.p, param.L);
4124 <<worm-like chain structure>>=
4125 typedef struct wlc_param_struct {
4126 double p; /* WLC persistence length (m) */
4127 double L; /* WLC contour length (m) */
4131 <<worm-like chain create/destroy function declarations>>=
4132 extern void *string_create_wlc_param_t(char **param_strings);
4133 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
4136 <<worm-like chain create/destroy functions>>=
4137 wlc_param_t *create_wlc_param_t(double p, double L)
4139 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
4140 assert(ret != NULL);
4146 void *string_create_wlc_param_t(char **param_strings)
4148 return create_wlc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4149 safe_strtod(param_strings[1], __FUNCTION__));
4152 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
4160 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.
4161 TODO (now we copy the string before parsing, but still do this for clarity).
4163 <<valid string write code>>=
4164 char string[] = "abc";
4167 <<valid string write code>>=
4168 char *string = "abc";
4172 <<worm-like chain tension model global declarations>>=
4173 extern int num_wlc_params;
4174 extern char *wlc_param_descriptions[];
4175 extern char wlc_param_string[];
4177 <<worm-like chain tension model globals>>=
4178 int num_wlc_params = 2;
4179 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
4180 char wlc_param_string[]="0.39e-9,28.4e-9";
4183 <<worm-like chain tension model getopt>>=
4184 {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}
4186 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
4188 \subsection{Freely-jointed chain}
4191 <<freely-jointed chain functions>>=
4192 <<freely-jointed chain function>>
4193 <<freely-jointed chain handler>>
4194 <<freely-jointed chain create/destroy functions>>
4197 <<freely-jointed chain tension model declarations>>=
4198 <<freely-jointed chain handler declaration>>
4199 <<freely-jointed chain create/destroy function declarations>>
4200 <<freely-jointed chain tension model global declarations>>
4204 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
4205 The tension of a chain of $N$ such links is given by
4207 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
4209 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}.
4210 <<freely-jointed chain function>>=
4211 double langevin(double x, void *params)
4214 return 1.0/tanh(x) - 1/x;
4216 <<one dimensional bracket declaration>>
4217 <<one dimensional solve declaration>>
4218 double inv_langevin(double x)
4220 double lb=0.0, ub=1.0, ret;
4221 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
4222 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
4226 double fjc(double x, double T, double l, int N)
4228 double L = l*(double)N;
4229 if (x == 0) return 0;
4230 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
4231 return k_B*T/l * inv_langevin(x/L);
4235 <<freely-jointed chain handler declaration>>=
4236 extern double fjc_handler(double x, void *pdata);
4238 <<freely-jointed chain handler>>=
4239 double fjc_handler(double x, void *pdata)
4241 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4242 list_t *list = data->group_tension_params;
4247 assert(list != NULL); /* empty active group?! */
4248 l = ((fjc_param_t *)list->d)->l;
4249 while (list != NULL) {
4250 /* enforce consistent l */
4251 assert( ((fjc_param_t *)list->d)->l == l);
4252 N += ((fjc_param_t *)list->d)->N;
4254 fprintf(stderr, "%s: FJC section %d links long in series\n",
4255 __FUNCTION__, ((fjc_param_t *)list->d)->N);
4260 fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
4261 __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
4263 return fjc(x, data->env->T, l, N);
4267 <<freely-jointed chain structure>>=
4268 typedef struct fjc_param_struct {
4269 double l; /* FJC link length (m) */
4270 int N; /* FJC number of links */
4274 <<freely-jointed chain create/destroy function declarations>>=
4275 extern void *string_create_fjc_param_t(char **param_strings);
4276 extern void destroy_fjc_param_t(void *p);
4279 <<freely-jointed chain create/destroy functions>>=
4280 fjc_param_t *create_fjc_param_t(double l, double N)
4282 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
4283 assert(ret != NULL);
4289 void *string_create_fjc_param_t(char **param_strings)
4291 return create_fjc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4292 safe_strtod(param_strings[1], __FUNCTION__));
4295 void destroy_fjc_param_t(void *p)
4302 <<freely-jointed chain tension model global declarations>>=
4303 extern int num_fjc_params;
4304 extern char *fjc_param_descriptions[];
4305 extern char fjc_param_string[];
4308 <<freely-jointed chain tension model globals>>=
4309 int num_fjc_params = 2;
4310 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
4311 char fjc_param_string[]="0.5e-9,200";
4314 <<freely-jointed chain tension model getopt>>=
4315 {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}
4318 \subsection{Tension model implementation}
4320 <<tension-model.c>>=
4323 <<tension model includes>>
4324 #include "tension_model.h"
4325 <<tension model internal definitions>>
4326 <<tension model globals>>
4327 <<tension model functions>>
4330 <<tension model includes>>=
4331 #include <assert.h> /* assert() */
4332 #include <stdlib.h> /* NULL, strto*() */
4333 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4334 #include <stdio.h> /* fprintf(), stdout */
4335 #include <errno.h> /* errno, ERANGE */
4338 #include "tension_balance.h" /* oneD_*() */
4341 <<tension model internal definitions>>=
4342 <<constant tension structure>>
4344 <<worm-like chain structure>>
4345 <<freely-jointed chain structure>>
4348 <<tension model globals>>=
4349 <<tension handler array global>>
4350 <<constant tension model globals>>
4351 <<hooke tension model globals>>
4352 <<worm-like chain tension model globals>>
4353 <<freely-jointed chain tension model globals>>
4356 <<tension model functions>>=
4358 <<constant tension functions>>
4360 <<worm-like chain functions>>
4361 <<freely-jointed chain functions>>
4365 \subsection{Utilities}
4367 The tension models can get complicated, and users may want to reassure
4368 themselves that this portion of the input physics are functioning properly. The
4369 stand-alone program [[tension_model_utils]] demonstrates the tension model
4370 interface without the overhead of having to understand a full simulation.
4371 [[tension_model_utils]] takes command line model arguments like the full
4372 simulation, and prints $F(x)$ data to the screen for the selected model over a
4375 <<tension-model-utils.c>>=
4377 <<tension model utility includes>>
4378 <<tension model utility definitions>>
4379 <<tension model utility getopt functions>>
4381 <<tension model utility main>>
4384 <<tension model utility main>>=
4385 int main(int argc, char **argv)
4387 <<tension model getopt array definition>>
4388 tension_model_getopt_t *model = NULL;
4392 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
4393 tension_handler_data_t tdata;
4394 int num_param_args; /* for INIT_MODEL() */
4395 char **param_args; /* for INIT_MODEL() */
4397 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model,
4398 &Fmax, &Xmax, &flags);
4400 if (flags & VFLAG) {
4401 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4403 INIT_MODEL("utility", model, model->params, params);
4404 tdata.group_tension_params = NULL;
4405 push(&tdata.group_tension_params, params, 1);
4407 tdata.persist = NULL;
4408 if (model->handler == NULL) {
4409 printf("No tension function for model %s\n", model->name);
4415 printf("#Distance (m)\tForce (N)\n");
4416 for (i=0; i<=N; i++) {
4417 x = Xmax*i/(double)N;
4418 F = (*model->handler)(x, &tdata);
4419 if (F < 0 || F > Fmax) break;
4420 printf("%g\t%g\n", x, F);
4422 if (flags & VFLAG && i <= N) {
4423 /* explain exit condition */
4425 printf("Impossible force %g\n", F);
4427 printf("Reached large force limit %g > %g\n", F, Fmax);
4430 params = pop(&tdata.group_tension_params);
4431 if (model->destructor)
4432 (*model->destructor)(params);
4437 <<tension model utility includes>>=
4440 #include <string.h> /* strlen() */
4441 #include <assert.h> /* assert() */
4442 #include <errno.h> /* errno, ERANGE */
4446 #include "tension_model.h"
4449 <<tension model utility definitions>>=
4450 <<version definition>>
4451 #define VFLAG 1 /* verbose */
4452 <<string comparison definition and globals>>
4453 <<tension model getopt definitions>>
4454 <<initialize model definition>>
4458 <<tension model utility getopt functions>>=
4460 <<index tension model>>
4461 <<tension model utility help>>
4462 <<tension model utility get options>>
4465 <<tension model utility help>>=
4466 void help(char *prog_name,
4468 int n_tension_models, tension_model_getopt_t *tension_models,
4469 int tension_model, double Fmax, double Xmax)
4472 printf("usage: %s [options]\n", prog_name);
4473 printf("Version %s\n\n", VERSION);
4474 printf("A utility for understanding the available tension models\n\n");
4475 printf("Environmental options:\n");
4476 printf("-T T\tTemperature T (currently %g K)\n", env->T);
4477 printf("-C T\tYou can also set the temperature T in Celsius\n");
4478 printf("Model options:\n");
4479 printf("-m model\tFolded tension model (currently %s)\n",
4480 tension_models[tension_model].name);
4481 printf("-a args\tFolded tension model argument string (currently %s)\n",
4482 tension_models[tension_model].params);
4483 printf("F(x) is calculated for a range of x and printed\n");
4484 printf("For example:\n");
4485 printf(" #Distance (m)\tForce (N)\n");
4486 printf(" 123.456\t7.89\n");
4488 printf("-F\tSet the maximum F value for the standard mode F(x) output (currently %g)\n", Fmax);
4489 printf("-X\tSet the maximum x value for the standard mode F(x) output (currently %g)\n", Xmax);
4490 printf("-V\tChange output to verbose mode\n");
4491 printf("-h\tPrint this help and exit\n");
4493 printf("Tension models:\n");
4494 for (i=0; i<n_tension_models; i++) {
4495 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
4496 for (j=0; j < tension_models[i].num_params ; j++)
4497 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
4498 printf("\t\tdefaults: %s\n", tension_models[i].params);
4503 <<tension model utility get options>>=
4504 void get_options(int argc, char **argv, environment_t *env,
4505 int n_tension_models, tension_model_getopt_t *tension_models,
4506 tension_model_getopt_t **model, double *Fmax, double *Xmax,
4507 unsigned int *flags)
4509 char *prog_name = NULL;
4510 char c, options[] = "T:C:m:a:F:X:Vh";
4511 int tension_model=0, num_strings;
4512 extern char *optarg;
4513 extern int optind, optopt, opterr;
4517 /* setup defaults */
4519 prog_name = argv[0];
4520 env->T = 300.0; /* K */
4521 *Fmax = 1e5; /* N */
4522 *Xmax = 1e-6; /* m */
4524 *model = tension_models;
4526 while ((c=getopt(argc, argv, options)) != -1) {
4528 case 'T': env->T = safe_strtod(optarg, "-T"); break;
4529 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
4531 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
4532 *model = tension_models+tension_model;
4535 tension_models[tension_model].params = optarg;
4537 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
4538 case 'X': *Xmax = safe_strtod(optarg, "-X"); break;
4539 case 'V': *flags |= VFLAG; break;
4541 fprintf(stderr, "unrecognized option '%c'\n", optopt);
4542 /* fall through to default case */
4544 help(prog_name, env, n_tension_models, tension_models, tension_model, *Fmax, *Xmax);
4553 \section{Unfolding rate models}
4554 \label{sec.k_models}
4556 We have two main choices for determining $k$: Bell's model, which assumes the
4557 folded domains exist in a single `folded' state and make instantaneous,
4558 stochastic transitions over a free energy barrier; and Kramers' model, which
4559 treats the unfolding as a thermalized diffusion process.
4560 We deal with these options like we dealt with the different tension models: by bundling all model-specific
4561 parameters into structures, with model specific functions for determining $k$.
4563 <<k func definition>>=
4564 typedef double k_func_t(double F, environment_t *env, void *params);
4567 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
4568 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]].
4570 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
4571 because \LaTeX\ doesn't like underscores outside math mode.
4572 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
4573 You could use spaces instead (`\verb+ +'),
4574 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
4575 which I don't like as much.
4581 <<k func definition>>
4582 <<null k declarations>>
4583 <<const k declarations>>
4584 <<bell k declarations>>
4585 <<kbell k declarations>>
4586 <<kramers k declarations>>
4587 <<kramers integ k declarations>>
4588 #endif /* K_MODEL_H */
4591 <<k model module makefile lines>>=
4592 NW_SAWSIM_MODS += k_model
4593 CHECK_BINS += check_k_model
4594 check_k_model_MODS = parse string_eval
4596 [[check_k_model]] also depends on the parse module.
4598 <<k model utils makefile lines>>=
4599 K_MODEL_MODS = k_model parse string_eval
4600 K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
4601 $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h)
4602 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
4603 notangle -Rk-model-utils.c $< > $@
4604 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
4605 gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
4606 clean_k_model_utils :
4607 rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
4613 For domains that are never folded.
4615 <<null k declarations>>=
4617 <<null k functions>>=
4622 <<null k model getopt>>=
4623 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
4626 \subsection{Constant rate model}
4629 This is a pretty boring model, but it is usefull for testing the engine.
4633 where $k_0$ is the constant unfolding rate.
4635 <<const k functions>>=
4636 <<const k function>>
4637 <<const k structure create/destroy functions>>
4640 <<const k declarations>>=
4641 <<const k function declaration>>
4642 <<const k structure create/destroy function declarations>>
4643 <<const k global declarations>>
4647 <<const k structure>>=
4648 typedef struct const_k_param_struct {
4649 double knot; /* transition rate k_0 (frac population per s) */
4653 <<const k function declaration>>=
4654 double const_k(double F, environment_t *env, void *const_k_params);
4656 <<const k function>>=
4657 double const_k(double F, environment_t *env, void *const_k_params)
4658 { /* Returns the rate constant k in frac pop per s. */
4659 const_k_param_t *p = (const_k_param_t *) const_k_params;
4661 assert(p->knot > 0);
4666 <<const k structure create/destroy function declarations>>=
4667 void *string_create_const_k_param_t(char **param_strings);
4668 void destroy_const_k_param_t(void *p);
4671 <<const k structure create/destroy functions>>=
4672 const_k_param_t *create_const_k_param_t(double knot)
4674 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
4675 assert(ret != NULL);
4680 void *string_create_const_k_param_t(char **param_strings)
4682 return create_const_k_param_t(safe_strtod(param_strings[0], __FUNCTION__));
4685 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
4692 <<const k global declarations>>=
4693 extern int num_const_k_params;
4694 extern char *const_k_param_descriptions[];
4695 extern char const_k_param_string[];
4698 <<const k globals>>=
4699 int num_const_k_params = 1;
4700 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
4701 char const_k_param_string[]="1";
4704 <<const k model getopt>>=
4705 {"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}
4708 \subsection{Bell's model}
4711 Quantitatively, Bell's model gives $k$ as
4713 k = k_0 \cdot e^{F dx / k_B T} \;,
4715 where $k_0$ is the unforced unfolding rate,
4716 $F$ is the applied unfolding force,
4717 $dx$ is the distance to the transition state, and
4718 $k_B T$ is the thermal energy\citep{schlierf06}.
4720 <<bell k functions>>=
4722 <<bell k structure create/destroy functions>>
4725 <<bell k declarations>>=
4726 <<bell k function declaration>>
4727 <<bell k structure create/destroy function declarations>>
4728 <<bell k global declarations>>
4732 <<bell k structure>>=
4733 typedef struct bell_param_struct {
4734 double knot; /* transition rate k_0 (frac population per s) */
4735 double dx; /* `distance to transition state' (nm) */
4739 <<bell k function declaration>>=
4740 double bell_k(double F, environment_t *env, void *bell_params);
4742 <<bell k function>>=
4743 double bell_k(double F, environment_t *env, void *bell_params)
4744 { /* Returns the rate constant k in frac pop per s.
4745 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4746 * uses global k_B in J/K */
4747 bell_param_t *p = (bell_param_t *) bell_params;
4748 assert(F >= 0); assert(env->T > 0);
4750 assert(p->knot > 0); //assert(p->dx >= 0); /* p < 0 for refolding reactions, etc. */
4752 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
4753 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
4754 p->knot * exp(F*p->dx / (k_B*env->T) ));
4755 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
4757 return p->knot * exp(F*p->dx / (k_B*env->T) );
4761 <<bell k structure create/destroy function declarations>>=
4762 void *string_create_bell_param_t(char **param_strings);
4763 void destroy_bell_param_t(void *p);
4766 <<bell k structure create/destroy functions>>=
4767 bell_param_t *create_bell_param_t(double knot, double dx)
4769 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
4770 assert(ret != NULL);
4776 void *string_create_bell_param_t(char **param_strings)
4778 return create_bell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4779 safe_strtod(param_strings[1], __FUNCTION__));
4782 void destroy_bell_param_t(void *p /* bell_param_t * */)
4789 <<bell k global declarations>>=
4790 extern int num_bell_params;
4791 extern char *bell_param_descriptions[];
4792 extern char bell_param_string[];
4796 int num_bell_params = 2;
4797 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4798 char bell_param_string[]="3.3e-4,0.25e-9";
4801 <<bell k model getopt>>=
4802 {"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}
4804 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4805 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4808 \subsection{Linker stiffness corrected Bell model}
4811 Walton et. al showed that the Bell model constant force approximation
4812 doesn't hold for biotin-streptavadin binding, and I extended their
4813 results to I27 for the 2009 Biophysical Society Annual
4814 Meeting\cite{walton08,king09}. More details to follow when I get done
4815 with the conference\ldots
4817 We adjust Bell's model to
4819 k = k_0 \cdot e^{\frac{F dx - \frac{1}{2} \kappa dx^2}{k_B T}} \;,
4821 where $k_0$ is the unforced unfolding rate,
4822 $F$ is the applied unfolding force,
4823 $\kappa$ is the effective spring constant,
4824 $dx$ is the distance to the transition state, and
4825 $k_B T$ is the thermal energy\citep{schlierf06}.
4827 <<kbell k functions>>=
4828 <<kbell k function>>
4829 <<kbell k structure create/destroy functions>>
4832 <<kbell k declarations>>=
4833 <<kbell k function declaration>>
4834 <<kbell k structure create/destroy function declarations>>
4835 <<kbell k global declarations>>
4839 <<kbell k structure>>=
4840 typedef struct kbell_param_struct {
4841 double knot; /* transition rate k_0 (frac population per s) */
4842 double dx; /* `distance to transition state' (nm) */
4846 <<kbell k function declaration>>=
4847 double kbell_k(double F, environment_t *env, void *kbell_params);
4849 <<kbell k function>>=
4850 double kbell_k(double F, environment_t *env, void *kbell_params)
4851 { /* Returns the rate constant k in frac pop per s.
4852 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4853 * uses global k_B in J/K */
4854 kbell_param_t *p = (kbell_param_t *) kbell_params;
4855 assert(F >= 0); assert(env->T > 0); assert(env->stiffness >= 0);
4857 assert(p->knot > 0); assert(p->dx >= 0);
4858 return p->knot * exp((F*p->dx - 0.5*env->stiffness*(p->dx)*(p->dx))/(k_B*env->T) );
4862 <<kbell k structure create/destroy function declarations>>=
4863 void *string_create_kbell_param_t(char **param_strings);
4864 void destroy_kbell_param_t(void *p);
4867 <<kbell k structure create/destroy functions>>=
4868 kbell_param_t *create_kbell_param_t(double knot, double dx)
4870 kbell_param_t *ret = (kbell_param_t *) malloc(sizeof(kbell_param_t));
4871 assert(ret != NULL);
4877 void *string_create_kbell_param_t(char **param_strings)
4879 return create_kbell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4880 safe_strtod(param_strings[1], __FUNCTION__));
4883 void destroy_kbell_param_t(void *p /* kbell_param_t * */)
4890 <<kbell k global declarations>>=
4891 extern int num_kbell_params;
4892 extern char *kbell_param_descriptions[];
4893 extern char kbell_param_string[];
4896 <<kbell k globals>>=
4897 int num_kbell_params = 2;
4898 char *kbell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4899 char kbell_param_string[]="3.3e-4,0.25e-9";
4902 <<kbell k model getopt>>=
4903 {"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}
4905 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4906 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4909 \subsection{Kramer's model}
4912 Kramer's model gives $k$ as
4914 % \frac{1}{k} = \frac{1}{D}
4915 % \int_{x_\text{min}}^{x_\text{max}}
4916 % e^{\frac{-U_F(x)}{k_B T}}
4917 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4920 %where $D$ is the diffusion constant,
4921 %$U_F(x)$ is the free energy along the unfolding coordinate, and
4922 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4923 % before the minimum and after the maximum, respectively \citep{schlierf06}.
4925 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
4927 where $D$ is the diffusion constant,
4929 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
4931 are length scales where
4932 $x_c(F)$ is the location of the bound state and
4933 $x_{ts}(F)$ is the location of the transition state,
4934 $E(F,x)$ is the free energy, and
4935 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
4936 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
4937 \citet{evans97} Eqn.~3).
4939 <<kramers k functions>>=
4940 <<kramers k function>>
4941 <<kramers k structure create/destroy functions>>
4944 <<kramers k declarations>>=
4945 <<kramers k function declaration>>
4946 <<kramers k structure create/destroy function declarations>>
4947 <<kramers k global declarations>>
4951 <<kramers k structure>>=
4952 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
4953 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
4955 typedef struct kramers_param_struct {
4956 double D; /* diffusion rate (um/s) */
4963 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
4964 //kramers_x_func_t *xts; /* function returning transition x (nm) */
4965 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
4966 //kramers_E_func_t *E; /* function returning E (J) */
4967 //double *E_params; /* parametrize the energy landscape */
4968 //int n_E_params; /* length of E_params array */
4972 <<kramers k function declaration>>=
4973 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
4974 extern double kramers_k(double F, environment_t *env, void *kramers_params);
4976 <<kramers k function>>=
4977 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
4979 static char input[160]={0};
4980 static char output[80]={0};
4981 /* setup the environment */
4982 fprintf(in, "F = %g; T = %g\n", F, T);
4983 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4984 string_eval(in, out, input, 80, output);
4985 return safe_strtod(output, __FUNCTION__);
4988 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
4990 static char input[160]={0};
4991 static char output[80]={0};
4992 /* setup the environment */
4993 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
4994 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4995 string_eval(in, out, input, 80, output);
4996 return safe_strtod(output, __FUNCTION__);
4999 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
5001 kramers_param_t *p = (kramers_param_t *) kramers_params;
5002 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
5005 double kramers_k(double F, environment_t *env, void *kramers_params)
5006 { /* Returns the rate constant k in frac pop per s.
5007 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5008 * uses global k_B in J/K */
5009 kramers_param_t *p = (kramers_param_t *) kramers_params;
5010 double kbT, xc, xts, lc, lts, Eb;
5011 assert(F >= 0); assert(env->T > 0);
5014 assert(p->xc != NULL); assert(p->xts != NULL);
5015 assert(p->ddEdxx != NULL); assert(p->E != NULL);
5016 assert(p->in != NULL); assert(p->out != NULL);
5018 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
5019 if (xc == -1.0) return -1.0; /* error (Too much force) */
5020 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
5021 if (xc == -1.0) return -1.0; /* error (Too much force) */
5022 lc = sqrt(2.0*M_PI*kbT /
5023 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
5024 lts = sqrt(-2.0*M_PI*kbT/
5025 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
5026 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
5027 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
5028 //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));
5029 return p->D/(lc*lts) * exp(-Eb/kbT);
5033 <<kramers k structure create/destroy function declarations>>=
5034 void *string_create_kramers_param_t(char **param_strings);
5035 void destroy_kramers_param_t(void *p);
5038 <<kramers k structure create/destroy functions>>=
5039 kramers_param_t *create_kramers_param_t(double D,
5040 char *xc_fn, char *xts_fn,
5041 char *E_fn, char *ddEdxx_fn,
5042 char *global_define_string)
5043 // kramers_x_func_t *xc, kramers_x_func_t *xts,
5044 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
5045 // double *E_params, int n_E_params)
5047 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
5048 assert(ret != NULL);
5049 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5050 assert(ret->xc != NULL);
5051 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5052 assert(ret->xts != NULL);
5053 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
5054 assert(ret->E != NULL);
5055 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
5056 assert(ret->ddEdxx != NULL);
5058 strcpy(ret->xc, xc_fn);
5059 strcpy(ret->xts, xts_fn);
5060 strcpy(ret->E, E_fn);
5061 strcpy(ret->ddEdxx, ddEdxx_fn);
5062 //ret->E_params = E_params;
5063 //ret->n_E_params = n_E_params;
5064 string_eval_setup(&ret->in, &ret->out);
5065 fprintf(ret->in, "kB = %g\n", k_B);
5066 fprintf(ret->in, "%s\n", global_define_string);
5070 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
5071 void *string_create_kramers_param_t(char **param_strings)
5073 return create_kramers_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5081 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
5083 kramers_param_t *pk = (kramers_param_t *)p;
5085 string_eval_teardown(&pk->in, &pk->out);
5087 // free(pk->E_params);
5093 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
5094 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.
5095 We follow \citet{shillcock98} and use
5097 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
5098 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
5101 where TODO, variable meanings.%\citep{shillcock98}.
5102 The first and second derivatives of $E(F,E_b,x,x_s)$ are
5104 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\\
5105 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
5107 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
5108 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
5109 The roots are therefor at
5111 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
5112 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
5113 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
5116 <<kramers k global declarations>>=
5117 extern int num_kramers_params;
5118 extern char *kramers_param_descriptions[];
5119 extern char kramers_param_string[];
5122 <<kramers k globals>>=
5123 int num_kramers_params = 6;
5124 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"};
5125 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)}";
5127 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
5128 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
5129 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
5130 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.
5131 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
5132 It works so long as [[val_1]] is not [[false]].
5134 <<kramers k model getopt>>=
5135 {"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}
5138 \subsection{Kramer's model (natural cubic splines)}
5139 \label{sec.kramers_integ}
5141 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
5142 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
5143 \citet{schlierf06} Eqn.~2)
5145 \frac{1}{k} = \frac{1}{D}
5146 \int_{x_\text{min}}^{x_\text{max}}
5147 e^{\frac{U_F(x)}{k_B T}}
5148 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5151 where $D$ is the diffusion constant,
5152 $U_F(x)$ is the free energy along the unfolding coordinate, and
5153 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
5154 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
5156 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
5157 interpolating between them with cubic splines.
5158 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
5160 <<kramers integ k functions>>=
5161 <<spline functions>>
5162 <<kramers integ A k functions>>
5163 <<kramers integ B k functions>>
5164 <<kramers integ k function>>
5165 <<kramers integ k structure create/destroy functions>>
5168 <<kramers integ k declarations>>=
5169 <<kramers integ k function declaration>>
5170 <<kramers integ k structure create/destroy function declarations>>
5171 <<kramers integ k global declarations>>
5175 <<kramers integ k structure>>=
5176 typedef double func_t(double x, void *params);
5177 typedef struct kramers_integ_param_struct {
5178 double D; /* diffusion rate (um/s) */
5179 double x_min; /* integration bounds */
5181 func_t *E; /* function returning E (J) */
5182 void *E_params; /* parametrize the energy landscape */
5183 destroy_data_func_t *destroy_E_params;
5185 double F; /* for passing into GSL evaluations */
5187 } kramers_integ_param_t;
5190 <<spline param structure>>=
5191 typedef struct spline_param_struct {
5192 int n_knots; /* natural cubic spline knots for U(x) */
5195 gsl_spline *spline; /* GSL spline parameters */
5196 gsl_interp_accel *acc; /* GSL spline acceleration data */
5200 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
5202 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5204 <<kramers integ A k functions>>=
5205 double kramers_integ_k_integrandA(double x, void *params)
5207 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5208 double U = (*p->E)(x, p->E_params);
5209 double Fx = p->F * x;
5210 double kbT = k_B * p->env->T;
5211 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
5212 return exp(-(U-Fx)/kbT);
5214 double kramers_integ_k_integralA(double x, void *params)
5216 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5218 double result, abserr;
5219 size_t neval; /* for qng */
5220 static gsl_integration_workspace *w = NULL;
5221 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5222 f.function = &kramers_integ_k_integrandA;
5224 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5225 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5226 //fprintf(stderr, "integralA = %g\n", result);
5232 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
5233 \int_{x_\text{min}}^{x_\text{max}}
5234 e^{\frac{U_F(x)}{k_B T}}
5235 \text{Integral}_A(x)
5238 <<kramers integ B k functions>>=
5239 double kramers_integ_k_integrandB(double x, void *params)
5241 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5242 double U = (*p->E)(x, p->E_params);
5243 double Fx = p->F * x;
5244 double kbT = k_B * p->env->T;
5245 double intA = kramers_integ_k_integralA(x, params);
5246 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
5247 return exp((U-Fx)/kbT)*intA;
5249 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
5251 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5253 double result, abserr;
5254 size_t neval; /* for qng */
5255 static gsl_integration_workspace *w = NULL;
5256 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5257 f.function = &kramers_integ_k_integrandB;
5261 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5262 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5263 //fprintf(stderr, "integralB = %g\n", result);
5268 With the integrals taken care of, evaluating $k$ is simply
5270 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
5272 <<kramers integ k function declaration>>=
5273 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
5275 <<kramers integ k function>>=
5276 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
5277 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
5278 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
5279 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
5283 <<kramers integ k structure create/destroy function declarations>>=
5284 void *string_create_kramers_integ_param_t(char **param_strings);
5285 void destroy_kramers_integ_param_t(void *p);
5288 <<kramers integ k structure create/destroy functions>>=
5289 kramers_integ_param_t *create_kramers_integ_param_t(double D,
5290 double xmin, double xmax, func_t *E,
5292 destroy_data_func_t *destroy_E_params)
5294 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
5295 assert(ret != NULL);
5300 ret->E_params = E_params;
5301 ret->destroy_E_params = destroy_E_params;
5305 void *string_create_kramers_integ_param_t(char **param_strings)
5307 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
5308 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
5309 void *E_params = string_create_spline_param_t(param_strings+1);
5310 return create_kramers_integ_param_t(
5311 safe_strtod(param_strings[0], __FUNCTION__),
5312 safe_strtod(param_strings[2], __FUNCTION__),
5313 safe_strtod(param_strings[3], __FUNCTION__),
5314 &spline_eval, E_params, destroy_spline_param_t);
5317 void destroy_kramers_integ_param_t(void *params)
5319 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5322 (*p->destroy_E_params)(p->E_params);
5328 Finally we have the GSL spline wrappers:
5329 <<spline functions>>=
5331 <<create/destroy spline>>
5334 <<spline function>>=
5335 double spline_eval(double x, void *spline_params)
5337 spline_param_t *p = (spline_param_t *)(spline_params);
5338 return gsl_spline_eval(p->spline, x, p->acc);
5342 <<create/destroy spline>>=
5343 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
5345 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
5346 assert(ret != NULL);
5347 ret->n_knots = n_knots;
5350 ret->acc = gsl_interp_accel_alloc();
5351 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
5352 gsl_spline_init(ret->spline, x, y, n_knots);
5356 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
5357 void *string_create_spline_param_t(char **param_strings)
5361 double *x=NULL, *y=NULL;
5362 /* break into ordered pairs */
5363 parse_list_string(param_strings[0], SEP, '(', ')',
5364 &num_knots, &knot_coords);
5365 x = (double *)malloc(sizeof(double)*num_knots);
5367 y = (double *)malloc(sizeof(double)*num_knots);
5369 for (i=0; i<num_knots; i++) {
5370 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
5371 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
5374 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
5376 free_parsed_list(num_knots, knot_coords);
5377 return create_spline_param_t(num_knots, x, y);
5380 void destroy_spline_param_t(void *params)
5382 spline_param_t *p = (spline_param_t *)params;
5385 gsl_spline_free(p->spline);
5387 gsl_interp_accel_free(p->acc);
5397 <<kramers integ k global declarations>>=
5398 extern int num_kramers_integ_params;
5399 extern char *kramers_integ_param_descriptions[];
5400 extern char kramers_integ_param_string[];
5403 <<kramers integ k globals>>=
5404 int num_kramers_integ_params = 4;
5405 char *kramers_integ_param_descriptions[] = {
5406 "diffusion D, m^2 / s",
5407 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
5408 "starting integration bound x_min, m",
5409 "ending integration bound x_max, m"};
5410 //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";
5411 //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";
5412 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";
5413 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
5414 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
5415 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
5418 <<kramers integ k model getopt>>=
5419 {"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}
5421 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
5422 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5424 \subsection{Unfolding model implementation}
5428 <<k model includes>>
5429 #include "k_model.h"
5430 <<k model internal definitions>>
5432 <<k model functions>>
5435 <<k model includes>>=
5436 #include <assert.h> /* assert() */
5437 #include <stdlib.h> /* NULL, malloc(), strto*() */
5438 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
5439 #include <stdio.h> /* fprintf(), stdout */
5440 #include <string.h> /* strlen(), strcpy() */
5441 #include <errno.h> /* errno, ERANGE */
5442 #include <gsl/gsl_integration.h>
5443 #include <gsl/gsl_spline.h>
5448 <<k model internal definitions>>=
5449 <<const k structure>>
5450 <<bell k structure>>
5451 <<kbell k structure>>
5452 <<kramers k structure>>
5453 <<kramers integ k structure>>
5454 <<spline param structure>>
5457 <<k model globals>>=
5462 <<kramers k globals>>
5463 <<kramers integ k globals>>
5466 <<k model functions>>=
5468 <<null k functions>>
5469 <<const k functions>>
5470 <<bell k functions>>
5471 <<kbell k functions>>
5472 <<kramers k functions>>
5473 <<kramers integ k functions>>
5476 \subsection{Unfolding model unit tests}
5478 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
5479 <<check-k-model.c>>=
5480 <<k model check includes>>
5481 <<check relative error>>
5483 <<k model test suite>>
5484 <<main check program>>
5487 <<k model check includes>>=
5488 #include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
5489 #include <stdio.h> /* sprintf() */
5490 #include <assert.h> /* assert() */
5491 #include <math.h> /* exp() */
5492 #include <errno.h> /* errno, ERANGE */
5495 #include "k_model.h"
5498 <<k model test suite>>=
5502 <<k model suite definition>>
5505 <<k model suite definition>>=
5506 Suite *test_suite (void)
5508 Suite *s = suite_create ("k model");
5509 <<const k test case defs>>
5510 <<bell k test case defs>>
5512 <<const k test case adds>>
5513 <<bell k test case adds>>
5518 \subsubsection{Constant}
5520 <<const k test case defs>>=
5521 TCase *tc_const_k = tcase_create("Constant k");
5524 <<const k test case adds>>=
5525 tcase_add_test(tc_const_k, test_const_k_create_destroy);
5526 tcase_add_test(tc_const_k, test_const_k_over_range);
5527 suite_add_tcase(s, tc_const_k);
5532 START_TEST(test_const_k_create_destroy)
5534 char *knot[] = {"1","2","3","4","5","6"};
5535 char *params[] = {knot[0]};
5538 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5539 params[0] = knot[i];
5540 p = string_create_const_k_param_t(params);
5541 destroy_const_k_param_t(p);
5546 START_TEST(test_const_k_over_range)
5548 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
5549 char *knot[] = {"1","2","3","4","5","6"};
5550 char *params[] = {knot[0]};
5553 char param_string[80];
5556 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5557 params[0] = knot[i];
5558 p = string_create_const_k_param_t(params);
5559 for ( F=Fm; F<FM; F+=dF ) {
5560 fail_unless(const_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
5561 "const_k(%g, %g, &{%s}) = %g != %s",
5562 F, T, knot[i], const_k(F, &env, p), knot[i]);
5564 destroy_const_k_param_t(p);
5570 \subsubsection{Bell}
5572 <<bell k test case defs>>=
5573 TCase *tc_bell = tcase_create("Bell's k");
5576 <<bell k test case adds>>=
5577 tcase_add_test(tc_bell, test_bell_k_create_destroy);
5578 tcase_add_test(tc_bell, test_bell_k_at_zero);
5579 tcase_add_test(tc_bell, test_bell_k_at_one);
5580 suite_add_tcase(s, tc_bell);
5584 START_TEST(test_bell_k_create_destroy)
5586 char *knot[] = {"1","2","3","4","5","6"};
5588 char *params[] = {knot[0], dx};
5591 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5592 params[0] = knot[i];
5593 p = string_create_bell_param_t(params);
5594 destroy_bell_param_t(p);
5599 START_TEST(test_bell_k_at_zero)
5601 double F=0.0, T=300.0;
5602 char *knot[] = {"1","2","3","4","5","6"};
5604 char *params[] = {knot[0], dx};
5607 char param_string[80];
5610 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5611 params[0] = knot[i];
5612 p = string_create_bell_param_t(params);
5613 fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
5614 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
5615 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
5616 destroy_bell_param_t(p);
5621 START_TEST(test_bell_k_at_one)
5624 char *knot[] = {"1","2","3","4","5","6"};
5626 char *params[] = {knot[0], dx};
5627 double F= k_B*T/safe_strtod(dx, __FUNCTION__);
5630 char param_string[80];
5633 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5634 params[0] = knot[i];
5635 p = string_create_bell_param_t(params);
5636 CHECK_ERR(1e-6, safe_strtod(knot[i], __FUNCTION__)*exp(1.0), bell_k(F, &env, p));
5637 //fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__)*exp(1.0),
5638 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
5639 // F, T, knot[i], dx, bell_k(F, &env, p), safe_strtod(knot[i], __FUNCTION__)*exp(1.0));
5640 destroy_bell_param_t(p);
5646 <<kramers k tests>>=
5649 <<kramers k test case def>>=
5652 <<kramers k test case add>>=
5655 <<k model function tests>>=
5656 <<check relative error>>
5658 START_TEST(test_something)
5660 double k=5, x=3, last_x=2.0, F;
5661 one_dim_fn_t *handlers[] = {&hooke};
5662 void *data[] = {&k};
5663 double xi[] = {0.0};
5665 int new_active_groups = 1;
5666 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
5667 fail_unless(F = k*x, NULL);
5672 \subsection{Utilities}
5674 The unfolding models can get complicated, and are sometimes hard to explain to others.
5675 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
5676 the overhead of having to understand a full simulation.
5677 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
5678 or special mode, where the operation depends on the specific model selected.
5680 <<k-model-utils.c>>=
5682 <<k model utility includes>>
5683 <<k model utility definitions>>
5684 <<k model utility getopt functions>>
5685 <<k model utility multi model E>>
5686 <<k model utility main>>
5689 <<k model utility main>>=
5690 int main(int argc, char **argv)
5692 <<k model getopt array definition>>
5693 k_model_getopt_t *model = NULL;
5698 int num_param_args; /* for INIT_MODEL() */
5699 char **param_args; /* for INIT_MODEL() */
5701 double special_xmin;
5702 double special_xmax;
5703 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
5704 &Fmax, &special_xmin, &special_xmax, &flags);
5705 if (flags & VFLAG) {
5706 printf("#initializing model %s with parameters %s\n", model->name, model->params);
5708 INIT_MODEL("utility", model, model->params, params);
5711 if (model->k == NULL) {
5712 printf("No k function for model %s\n", model->name);
5716 printf("Fmax = %g <= 0\n", Fmax);
5722 printf("#F (N)\tk (%% pop. per s)\n");
5723 for (i=0; i<=N; i++) {
5724 F = Fmax*i/(double)N;
5725 k = (*model->k)(F, &env, params);
5727 printf("%g\t%g\n", F, k);
5732 if (model->k == NULL || model->k == &bell_k) {
5733 printf("No special mode for model %s\n", model->name);
5736 if (special_xmax <= special_xmin) {
5737 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
5741 double Xrng=(special_xmax-special_xmin), x, E;
5743 printf("#x (m)\tE (J)\n");
5744 for (i=0; i<=N; i++) {
5745 x = special_xmin + Xrng*i/(double)N;
5746 E = multi_model_E(model->k, params, &env, x);
5747 printf("%g\t%g\n", x, E);
5754 if (model->destructor)
5755 (*model->destructor)(params);
5760 <<k model utility multi model E>>=
5761 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
5763 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
5765 if (k == kramers_integ_k)
5766 E = (*p->E)(x, p->E_params);
5768 E = kramers_E(0, env, model_params, x);
5774 <<k model utility includes>>=
5777 #include <string.h> /* strlen() */
5778 #include <assert.h> /* assert() */
5779 #include <errno.h> /* errno, ERANGE */
5782 #include "string_eval.h"
5783 #include "k_model.h"
5786 <<k model utility definitions>>=
5787 <<version definition>>
5788 #define VFLAG 1 /* verbose */
5789 enum mode_t {M_K_OF_F, M_SPECIAL};
5790 <<string comparison definition and globals>>
5791 <<k model getopt definitions>>
5792 <<initialize model definition>>
5793 <<kramers k structure>>
5794 <<kramers integ k structure>>
5798 <<k model utility getopt functions>>=
5801 <<k model utility help>>
5802 <<k model utility get options>>
5805 <<k model utility help>>=
5806 void help(char *prog_name,
5808 int n_k_models, k_model_getopt_t *k_models,
5809 int k_model, double Fmax, double special_xmin, double special_xmax)
5812 printf("usage: %s [options]\n", prog_name);
5813 printf("Version %s\n\n", VERSION);
5814 printf("A utility for understanding the available unfolding models\n\n");
5815 printf("Environmental options:\n");
5816 printf("-T T\tTemperature T (currently %g K)\n", env->T);
5817 printf("-C T\tYou can also set the temperature T in Celsius\n");
5818 printf("Model options:\n");
5819 printf("-k model\tTransition rate model (currently %s)\n",
5820 k_models[k_model].name);
5821 printf("-K args\tTransition rate model argument string (currently %s)\n",
5822 k_models[k_model].params);
5823 printf("There are two output modes. In standard mode, k(F) is printed\n");
5824 printf("For example:\n");
5825 printf(" #Force (N)\tk (% pop. per s)\n");
5826 printf(" 123.456\t7.89\n");
5828 printf("In special mode, the output depends on the model.\n");
5829 printf("For models defining an energy function E(x), that function is printed\n");
5830 printf(" #Distance (m)\tE (J)\n");
5831 printf(" 0.0012\t0.0034\n");
5833 printf("-m\tChange output to standard mode\n");
5834 printf("-M\tChange output to special mode\n");
5835 printf("-F\tSet the maximum F value for the standard mode k(F) output (currently %g)\n", Fmax);
5836 printf("-x\tSet the minimum x value for the special mode E(x) output (currently %g)\n", special_xmin);
5837 printf("-X\tSet the maximum x value for the special mode E(x) output (currently %g)\n", special_xmax);
5838 printf("-V\tChange output to verbose mode\n");
5839 printf("-h\tPrint this help and exit\n");
5841 printf("Unfolding rate models:\n");
5842 for (i=0; i<n_k_models; i++) {
5843 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
5844 for (j=0; j < k_models[i].num_params ; j++)
5845 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
5846 printf("\t\tdefaults: %s\n", k_models[i].params);
5851 <<k model utility get options>>=
5852 void get_options(int argc, char **argv, environment_t *env,
5853 int n_k_models, k_model_getopt_t *k_models,
5854 enum mode_t *mode, k_model_getopt_t **model,
5855 double *Fmax, double *special_xmin, double *special_xmax,
5856 unsigned int *flags)
5858 char *prog_name = NULL;
5859 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
5861 extern char *optarg;
5862 extern int optind, optopt, opterr;
5866 /* setup defaults */
5868 prog_name = argv[0];
5869 env->T = 300.0; /* K */
5873 *Fmax = 1e-9; /* N */
5874 *special_xmax = 1e-8;
5875 *special_xmin = 0.0;
5877 while ((c=getopt(argc, argv, options)) != -1) {
5879 case 'T': env->T = safe_strtod(optarg, "-T"); break;
5880 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
5882 k_model = index_k_model(n_k_models, k_models, optarg);
5883 *model = k_models+k_model;
5886 k_models[k_model].params = optarg;
5888 case 'm': *mode = M_K_OF_F; break;
5889 case 'M': *mode = M_SPECIAL; break;
5890 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
5891 case 'x': *special_xmin = safe_strtod(optarg, "-x"); break;
5892 case 'X': *special_xmax = safe_strtod(optarg, "-X"); break;
5893 case 'V': *flags |= VFLAG; break;
5895 fprintf(stderr, "unrecognized option '%c'\n", optopt);
5896 /* fall through to default case */
5898 help(prog_name, env, n_k_models, k_models, k_model, *Fmax, *special_xmin, *special_xmax);
5907 \section{\LaTeX\ documentation}
5909 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
5910 The comment blocks are extracted (with nicely formatted code blocks), using
5911 <<latex makefile lines>>=
5912 $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5913 noweave -latex -index -delay $< > $@
5914 $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR)
5918 <<latex makefile lines>>=
5919 $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
5921 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5922 $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
5923 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5924 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5925 mv $(BUILD_DIR)/sawsim.pdf $@
5927 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
5928 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
5929 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
5931 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
5932 <<latex makefile lines>>=
5934 rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
5935 $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
5936 $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.toc \
5937 $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib
5938 clean_latex : semi-clean_latex
5939 rm -f $(DOC_DIR)/sawsim.pdf
5944 [[make]] is a common utility on *nix systems for managing dependencies while building software.
5945 In this case, we have several \emph{targets} that we'd like to build:
5946 the main simulation program \prog;
5947 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
5948 the documentation [[sawsim.pdf]];
5949 and the [[Makefile]] itself.
5950 Besides the generated files, there is the utility target
5951 [[clean]] for removing all generated files except [[Makefile]].
5953 % [[dist]] for generating a distributable tar file.
5955 Extract the makefile with
5956 `[[notangle -Rmakefile src/sawsim.nw | sed 's/ /\t/' > Makefile]]'.
5957 The sed is needed because notangle replaces tabs with eight spaces,
5958 but [[make]] needs tabs.
5960 # Decide what directories to put things in
5965 # Note: Cannot use, for example, `./build', because make eats the `./'
5966 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5968 # Modules (X.c, X.h) defined in the noweb file
5971 # Modules defined outside the noweb file
5972 FREE_SAWSIM_MODS = interp tavl
5974 <<list module makefile lines>>
5975 <<tension balance module makefile lines>>
5976 <<tension model module makefile lines>>
5977 <<k model module makefile lines>>
5978 <<parse module makefile lines>>
5979 <<string eval module makefile lines>>
5980 <<accel k module makefile lines>>
5982 SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS)
5984 # Everything needed for sawsim under one roof. sawsim.c must come first
5985 SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \
5986 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h)
5987 # Libraries needed to compile sawsim
5988 LIBS = gsl gslcblas m
5989 CHECK_LIBS = gsl gslcblas m check
5990 # The non-check binaries generated
5991 BINS = sawsim tension_model_utils k_model_utils sawsim_profile
5994 # Define the major targets
5995 all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
5997 view : $(DOC_DIR)/sawsim.pdf
5999 profile : $(BIN_DIR)/sawsim_profile | $(BIN_DIR)
6000 $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \
6001 -mnull -Mwlc -A0.39e-9,28e-9 -F8
6002 gprof $(BIN_DIR)/sawsim_profile gmon.out > $@
6003 check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim
6004 $(SHELL) -e -c 'for B in $^; do ./$$B; done'
6005 clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \
6006 clean_tension_model_utils \
6007 clean_k_model_utils clean_latex clean_check_sawsim
6008 rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \
6009 $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \
6010 $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
6011 $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
6012 $(BUILD_DIR)/global.h ./gmon.out
6013 $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
6015 # Various builds of sawsim
6016 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
6017 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6018 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
6019 gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6020 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
6021 gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6023 # Create the directories
6024 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
6027 # Copy over the source external to sawsim
6028 # Note: Cannot use, for example, `./build', because make eats the `./'
6029 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
6031 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6035 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6039 ## Basic source generated with noweb
6040 # The central files sawsim.c and global.h...
6041 $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6043 $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6044 notangle -Rglobal.h $< > $@
6045 # ... and the modules
6046 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6047 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6048 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6049 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6050 $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6051 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6052 # Note: I use `_' as a space in my file names, but noweb doesn't like
6053 # that so we use `-' to name the noweb roots and substitute here.
6055 ## Building the unit-test programs
6057 $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
6058 notangle -Rchecks $< > $@
6059 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
6060 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
6061 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR)
6062 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
6063 clean_check_sawsim :
6064 rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
6065 # ... and the modules
6067 $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
6068 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \
6069 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \
6070 $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6071 $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6072 $$(BUILD_DIR)/global.h | $(BIN_DIR)
6073 gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
6074 $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
6076 # todo: clean up the required modules to
6077 $(CHECK_BINS:%=clean_%) :
6078 rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
6080 # Cleaning up the modules
6082 $(SAWSIM_MODS:%=clean_%) :
6083 rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h)
6085 <<tension model utils makefile lines>>
6086 <<k model utils makefile lines>>
6087 <<latex makefile lines>>
6089 Makefile : $(SRC_DIR)/sawsim.nw
6090 notangle -Rmakefile $< | sed 's/ /\t/' > $@
6092 This is fairly self-explanatory. We compile a dynamically linked
6093 version ([[sawsim]]) and a statically linked version
6094 ([[sawsim_static]]). The static version is about 9 times as big, but
6095 it runs without needing dynamic GSL libraries to link to, so it's a
6096 better format for distributing to a cluster for parallel evaluation.
6100 \subsection{Hookean springs in series}
6101 \label{sec.math_hooke}
6104 The effective spring constant for $n$ springs $k_i$ in series is given by
6106 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6112 F &= k_i x_i &&\forall i \le n \\
6113 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
6114 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
6115 F &= k_1 x_1 = k_\text{eff} x \\
6116 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
6117 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6122 \addcontentsline{toc}{section}{References}
6123 \bibliography{sawsim}