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);
1129 } else if (errno == ERANGE) {
1130 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1131 assert(errno != ERANGE);
1136 static double safe_strtod(const char *s, const char *description) {
1140 ret = strtod(s, &endp);
1141 if (endp[0] != '\0') { //strlen(endp) > 0) {
1142 fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n",
1143 endp, s, description);
1144 assert(1==0); //strlen(endp) == 0);
1145 } else if (errno == ERANGE) {
1146 fprintf(stderr, "Error: '%s' out of range for %s\n", s, description);
1147 assert(errno != ERANGE);
1156 \addcontentsline{toc}{section}{Appendicies}
1158 \section{sawsim.c details}
1162 The general layout of our simulation code is:
1173 We include [[math.h]], so don't forget to link to the libm with `-lm'.
1175 #include <assert.h> /* assert() */
1176 #include <stdlib.h> /* malloc(), free(), rand(), strto*() */
1177 #include <stdio.h> /* fprintf(), stdout */
1178 #include <string.h> /* strlen, strtok() */
1179 #include <math.h> /* exp(), M_PI, sqrt() */
1180 #include <sys/types.h> /* pid_t (returned by getpid()) */
1181 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
1182 #include <errno.h> /* errno, ERANGE (for safe_strto*()) */
1185 #include "tension_balance.h"
1186 #include "k_model.h"
1187 #include "tension_model.h"
1189 #include "accel_k.h"
1193 <<version definition>>
1194 <<flag definitions>>
1195 <<max/min definitions>>
1196 <<string comparison definition and globals>>
1197 <<initialize model definition>>
1198 <<get options definitions>>
1199 <<state definitions>>
1200 <<transition definitions>>
1209 <<create/destroy state>>
1210 <<destroy state list>>
1212 <<create/destroy transition>>
1213 <<destroy transition list>>
1214 <<print state graph>>
1216 <<simulation functions>>
1217 <<boilerplate functions>>
1220 <<boilerplate functions>>=
1222 <<get option functions>>
1228 srand(getpid()*time(NULL)); /* seed rand() */
1232 <<flag definitions>>=
1233 /* in octal b/c of prefixed '0' */
1234 #define FLAG_OUTPUT_FULL_CURVE 01
1235 #define FLAG_OUTPUT_TRANSITION_FORCES 02
1239 static unsigned long int flags = 0;
1242 \subsection{Utilities}
1245 <<max/min definitions>>=
1246 #define MAX(a,b) ((a)>(b) ? (a) : (b))
1247 #define MIN(a,b) ((a)<(b) ? (a) : (b))
1250 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
1251 <<string comparison definition and globals>>=
1252 // Check if two strings match, return 1 if they do
1253 static char *temp_string_A;
1254 static char *temp_string_B;
1255 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
1256 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
1257 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
1258 /* +1 to also compare the '\0' */
1261 We also define a macro for our [[check]] unit testing
1262 <<check relative error>>=
1263 #define CHECK_ERR(max_err, expected, received) \
1265 fail_unless( (received-expected)/expected < max_err, \
1266 "relative error %g >= %g in %s (Expected %g, received %g)", \
1267 (received-expected)/expected, max_err, #received, \
1268 expected, received); \
1269 fail_unless(-(received-expected)/expected < max_err, \
1270 "relative error %g >= %g in %s (Expected %g, received %g)", \
1271 -(received-expected)/expected, max_err, #received, \
1272 expected, received); \
1277 int happens(double probability)
1279 assert(probability >= 0.0); assert(probability <= 1.0);
1280 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*/
1285 \subsection{Adaptive timesteps}
1286 \label{sec.adaptive_dt_implementation}
1288 $F(x)$ increases with $x$, possibly exploding (e.g.\ in the worm-like
1289 chain model), so basing the timestep on the the unfolding probability
1290 at the current tension is dangerous, and we need to search for a $dt$
1291 for which $P(F(x+v*dt)) < P_\text{target}$. There are two cases to
1292 consider. In the most common, no domains have unfolded since the last
1293 step, and we expect the next step to be slightly shorter than the
1294 previous one. In the less common, domains did unfold in the last
1295 step, and we expect the next step to be considerably longer than the
1298 double search_dt(transition_t *transition,
1300 environment_t *env, double x,
1301 double target_prob, double max_dt, double v)
1302 { /* Returns a valid timestep dt in seconds for the given transition.
1303 * Takes a list of all states, a pointer env to the environmental data,
1304 * a starting separation x in m, a target_prob between 0 and 1,
1305 * max_dt in s, stretching velocity v in m/s.
1307 double F, k, dtCur, dtU, dtUCur, dtL, dt;
1309 /* get upper bound using the current position */
1310 F = find_tension(state_list, env, x, 0); /* BUG. repeated calculation */
1311 //printf("Start with x = %g (F = %g)\n", x, F);
1312 k = accel_k(transition->k, F, env, transition->k_params);
1313 //printf("x %g\tF %g\tk %g\n", x, F, k);
1314 dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */
1316 //printf("overshot max_dt\n");
1319 if (v == 0) /* discrete stepping, so force is temporarily constant */
1322 /* set a lower bound on dt too */
1325 /* The dt determined above may produce illegitimate forces or ks.
1326 * Reduce the upper bound until we have valid ks. */
1328 F = find_tension(state_list, env, x+v*dt, 0);
1329 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1332 F = find_tension(state_list, env, x+v*dt, 0);
1334 //printf("Try for dt = %g (F = %g)\n", dt, F);
1335 k = accel_k(transition->k, F, env, transition->k_params);
1336 /* returned k may be -1.0 */
1337 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1338 while (k == -1.0) { /* reduce step until we hit a valid k */
1340 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1341 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 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1346 assert(dtU > 1e-14); /* timestep to valid k too small */
1347 dtUCur = target_prob / k; /* safe timestep back from x+dtU */
1349 return dt; /* dtU is safe. */
1351 /* dtU wasn't safe, lets see what would be. */
1352 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1353 dt = (dtU + dtL) / 2.0;
1354 F = find_tension(state_list, env, x+v*dt, 0);
1355 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1356 k = accel_k(transition->k, F, env, transition->k_params);
1357 dtCur = target_prob / k;
1358 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1359 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1361 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1362 dtU = dt; /* ... stepping out only dtCur would be safe */
1365 break; /* dtCur = dt */
1367 return MAX(dtUCur, dtL);
1372 To determine $dt$ for an array of potentially different folded domains,
1373 we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
1376 double determine_dt(list_t *state_list, list_t *transition_list,
1377 environment_t *env, double x,
1378 double target_prob, double dt_max, double v)
1379 { /* Returns the timestep dt in seconds.
1380 * Takes the state and transition lists, a pointer to the environment,
1381 * an extension x in m, target_prob between 0 and 1, dt_max in s,
1383 double dt=dt_max, new_dt;
1384 assert(target_prob > 0.0); assert(target_prob < 1.0);
1385 assert(dt_max > 0.0);
1386 assert(state_list != NULL);
1387 assert(transition_list != NULL);
1388 transition_list = head(transition_list);
1389 while (transition_list != NULL) {
1390 if (T(transition_list)->initial_state->num_domains > 0) {
1391 new_dt = search_dt(T(transition_list), state_list, env, x, target_prob, dt, v);
1392 dt = MIN(dt, new_dt);
1394 transition_list = transition_list->next;
1401 \subsection{State data}
1402 \label{sec.state_data}
1404 The domains exist in one of an array of possible states. Each state
1405 has a tension model which determines it's force/extension curve
1406 (possibly [[null]] for rigid states, see Appendix
1407 \ref{sec.null_tension}).
1409 Groups are, for example, for two domain states with WLC tensions; one
1410 with $p=0.4\U{nm}$, $L=10\U{nm}$ and the other with $p=0.4\U{nm}$,
1411 $L=28\U{nm}$. Since the persistence lengths are the same, you would
1412 like to add the contour lengths of all the domains in \emph{both}
1413 states, and plug that total contour length into a single WLC
1414 evaluation (see Section \ref{sec.find_tension}).
1415 <<state definitions>>=
1416 typedef struct state_struct {
1417 char *name; /* name string */
1418 one_dim_fn_t *tension_handler; /* tension handler */
1419 one_dim_fn_t *inverse_tension_handler; /* inverse tension handler */
1420 int tension_group; /* for grouping similar tension states */
1421 void *tension_params; /* pointer to tension parameters */
1422 destroy_data_func_t *destroy_tension;
1423 int num_domains; /* number of domains currently in state */
1424 /* cached lookups generated from the above data */
1425 int num_out_trans; /* length of out_trans array */
1426 struct transition_struct **out_trans; /* transitions leaving this state */
1427 int num_grouped_states; /* length of grouped states array */
1428 struct state_struct **grouped_states; /* states grouped with this one (doesn't include self) */
1431 /* get the state data for the current list node */
1432 #define S(list) ((state_t *)(list)->d)
1434 @ [[tension_params]] is a pointer to the parameters used by the
1435 handler function when determining the tension. [[destroy_tension]]
1436 points to a function for freeing [[tension_params]]. We it in the
1437 state structure so that [[destroy_state]] and will have an easier job
1440 [[create_]] and [[destroy_state]] are simple wrappers around
1441 [[malloc]] and [[free]].
1442 <<create/destroy state>>=
1443 state_t *create_state(char *name,
1444 one_dim_fn_t *tension_handler,
1445 one_dim_fn_t *inverse_tension_handler,
1447 void *tension_params,
1448 destroy_data_func_t *destroy_tension,
1451 state_t *ret = (state_t *)malloc(sizeof(state_t));
1452 assert(ret != NULL);
1453 /* make a copy of the name, so we aren't dependent on the original */
1454 ret->name = (char *)malloc(sizeof(char)*(strlen(name)+1));
1455 assert(ret->name != NULL);
1456 strcpy(ret->name, name); /* we know ret->name is long enough */
1458 ret->tension_handler = tension_handler;
1459 ret->inverse_tension_handler = inverse_tension_handler;
1460 ret->tension_group = tension_group;
1461 ret->tension_params = tension_params;
1462 ret->destroy_tension = destroy_tension;
1463 ret->num_domains = num_domains;
1465 ret->num_out_trans = 0;
1466 ret->out_trans = NULL;
1467 ret->num_grouped_states = 0;
1468 ret->grouped_states = NULL;
1471 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);
1476 void destroy_state(state_t *state)
1480 fprintf(stderr, "freeing state %s (%p)\n", state->name, state);
1484 if (state->destroy_tension)
1485 (*state->destroy_tension)(state->tension_params);
1486 if (state->out_trans)
1487 free(state->out_trans);
1488 if (state->grouped_states)
1489 free(state->grouped_states);
1496 We'll be storing the states in a list (see Appendix \ref{sec.lists}),
1497 so we also define a convenience function for cleaning up.
1498 <<destroy state list>>=
1499 void destroy_state_list(list_t *state_list)
1501 state_list = head(state_list);
1502 while (state_list != NULL)
1503 destroy_state((state_t *) pop(&state_list));
1508 We can also define a convenient way to get a state index from it's name.
1510 int state_index(list_t *state_list, char *name)
1513 state_list = head(state_list);
1514 while (state_list != NULL) {
1515 if (STRMATCH(S(state_list)->name, name)) return i;
1516 state_list = state_list->next;
1524 \subsection{Transition data}
1525 \label{sec.transition_data}
1527 Once you've created states (Sections \ref{sec.main},
1528 \ref{sec.model_selection}, and \ref{sec.state_data}), you need to
1529 create transitions between them.
1530 <<transition definitions>>=
1531 typedef struct transition_struct {
1532 state_t *initial_state;
1533 state_t *final_state;
1534 k_func_t *k; /* transition rate model */
1535 void *k_params; /* pointer to k parameters */
1536 destroy_data_func_t *destroy_k;
1539 /* get the transition data for the current list node */
1540 #define T(list) ((transition_t *)(list)->d)
1541 @ [[k]] is a pointer to the function determining the transition rate
1542 for a given tension. [[k_params]] and [[destroy_k]] have similar
1543 roles with regards to [[k]] as [[tension_params]] and
1544 [[destroy_tension]] have with regards to [[tension_handler]] (see
1545 Section \ref{sec.state_data}).
1547 [[create_]] and [[destroy_transition]] are simple wrappers around
1548 [[malloc]] and [[free]].
1549 <<create/destroy transition>>=
1550 transition_t *create_transition(state_t *initial_state, state_t *final_state,
1553 destroy_data_func_t *destroy_k)
1555 transition_t *ret = (transition_t *)malloc(sizeof(transition_t));
1556 assert(ret != NULL);
1557 assert(initial_state != NULL);
1558 assert(final_state != NULL);
1559 ret->initial_state = initial_state;
1560 ret->final_state = final_state;
1562 ret->k_params = k_params;
1563 ret->destroy_k = destroy_k;
1565 fprintf(stderr, "generated transition %p from %s to %s, params %p\n", ret, ret->initial_state->name, ret->final_state->name, ret->k_params);
1570 void destroy_transition(transition_t *transition)
1574 fprintf(stderr, "freeing transition %p (%s to %s)\n", transition, transition->initial_state->name, transition->final_state->name);
1576 if (transition->destroy_k)
1577 (*transition->destroy_k)(transition->k_params);
1584 We'll be storing the transitions in a list (see Appendix
1585 \ref{sec.lists}), so we also define a convenience function for
1587 <<destroy transition list>>=
1588 void destroy_transition_list(list_t *transition_list)
1590 transition_list = head(transition_list);
1591 while (transition_list != NULL)
1592 destroy_transition((transition_t *) pop(&transition_list));
1597 \subsection{Graphviz output}
1598 \label{sec.graphviz_output}
1600 <<print state graph>>=
1601 void print_state_graph(FILE *file, list_t *state_list, list_t *transition_list)
1603 state_list = head(state_list);
1604 transition_list = head(transition_list);
1605 fprintf(file, "digraph states {\n node [shape=record,color=black,fillcolor=white,style=filled];\n\n");
1607 fprintf(file, " n%d [label=\"%s\"]\n", state_index(state_list, S(state_list)->name), S(state_list)->name);
1608 if (state_list->next == NULL) break;
1609 state_list = state_list->next;
1611 fprintf(file, "\n");
1612 while (transition_list != NULL) {
1613 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));
1614 transition_list = transition_list->next;
1616 fprintf(file, "}\n");
1620 \subsection{Domain model and group handling}
1622 <<group functions>>=
1623 <<crosslink groups>>
1624 <<get active groups>>
1627 Scan through a list of states and construct an array of all the
1629 <<get active groups>>=
1630 void get_active_groups(list_t *state_list,
1631 int *pNum_active_groups,
1632 one_dim_fn_t ***pPer_group_handlers,
1633 one_dim_fn_t ***pPer_group_inverse_handlers,
1634 void ***pPer_group_data /* pointer to array of pointers to tension_handler_data_t */)
1636 void **active_groups=NULL;
1637 one_dim_fn_t *handler, *old_handler, *inverse_handler, *old_inverse_handler;
1638 one_dim_fn_t **per_group_handlers=NULL, **per_group_inverse_handlers=NULL;
1639 void **per_group_data=NULL; /* array of pointers to tension_handler_data_t */
1640 int i, j, num_domains, group, old_group, num_active_groups=0;
1641 list_t *active_states_list=NULL;
1642 tension_handler_data_t *tdata=NULL;
1645 /* get all the active groups in a list */
1646 state_list = head(state_list);
1648 fprintf(stderr, "%s:\t", __FUNCTION__);
1649 list_print(stderr, state_list, "states");
1651 while (state_list != NULL) {
1652 state = S(state_list);
1653 handler = state->tension_handler;
1654 inverse_handler = state->inverse_tension_handler;
1655 group = state->tension_group;
1656 num_domains = state->num_domains;
1657 if (list_index(active_states_list, state) == -1) {
1658 /* we haven't added this state (or it's associated group) yet */
1659 if (num_domains > 0 && handler != NULL) { /* add the state */
1660 num_active_groups++;
1661 push(&active_states_list, state, 1);
1663 fprintf(stderr, "%s:\tactive %s state (%p) with %d domain(s)\n", __FUNCTION__, state->name, state, state->num_domains);
1665 for (i=0; i < state->num_grouped_states; i++) {
1666 if (state->grouped_states[i]->num_domains > 0) {
1667 /* add this related (and active) state now */
1668 assert(state->grouped_states[i]->tension_handler == handler);
1669 assert(state->grouped_states[i]->tension_group == group);
1670 push(&active_states_list, state->grouped_states[i], 1);
1672 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);
1676 } /* else empty state or NULL handler */
1677 } /* else state already added */
1678 state_list = state_list->next;
1681 fprintf(stderr, "%s:\t", __FUNCTION__);
1682 list_print(stderr, active_states_list, "active states");
1685 assert(num_active_groups <= list_length(active_states_list));
1686 assert(num_active_groups > 0);
1688 /* allocate the arrays we'll be returning */
1689 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1690 assert(per_group_handlers != NULL);
1691 per_group_inverse_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1692 assert(per_group_inverse_handlers != NULL);
1693 per_group_data = (void **)malloc(sizeof(void *)*num_active_groups);
1694 assert(per_group_data != NULL);
1696 fprintf(stderr, "%s:\tallocating per group data %p:", __FUNCTION__, per_group_data);
1698 for (i=0; i<num_active_groups; i++) { /* we might want to recycle any group data structs allocated in the previous round */
1699 per_group_data[i] = malloc(sizeof(tension_handler_data_t));
1700 assert(per_group_data[i] != NULL);
1702 fprintf(stderr, " %d,%p", i, per_group_data[i]);
1706 fprintf(stderr, "\n");
1709 /* populate the arrays we'll be returning */
1710 active_states_list = head(active_states_list);
1712 old_inverse_handler = NULL;
1713 j = -1; /* j is the current active _group_ index */
1714 while (active_states_list != NULL) {
1715 state = (state_t *) pop(&active_states_list);
1716 handler = state->tension_handler;
1717 inverse_handler = state->inverse_tension_handler;
1718 group = state->tension_group;
1719 num_domains = state->num_domains;
1720 if (handler != old_handler || inverse_handler != old_inverse_handler || group != old_group) {
1721 j++; /* move to the next group */
1722 tdata = (tension_handler_data_t *) per_group_data[j];
1723 per_group_handlers[j] = handler;
1724 per_group_inverse_handlers[j] = inverse_handler;
1725 tdata->group_tension_params = NULL;
1727 tdata->persist = NULL; /* we might want to grab any appropriate data from the old persist here... */
1730 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);
1732 for (i=0; i<num_domains; i++) { /* add num_domains copies of the tension params */
1733 push(&tdata->group_tension_params, state->tension_params, 1);
1736 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);
1738 old_handler = handler;
1739 old_inverse_handler = inverse_handler;
1743 /* free old groups */
1744 if (*pPer_group_handlers != NULL)
1745 free(*pPer_group_handlers);
1746 if (*pPer_group_inverse_handlers != NULL)
1747 free(*pPer_group_inverse_handlers);
1748 if (*pPer_group_data != NULL) {
1749 for (i=0; i<*pNum_active_groups; i++)
1750 free((*pPer_group_data)[i]);
1751 free(*pPer_group_data);
1754 /* set new groups */
1756 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]);
1758 *pNum_active_groups = num_active_groups;
1759 *pPer_group_handlers = per_group_handlers;
1760 *pPer_group_inverse_handlers = per_group_inverse_handlers;
1761 *pPer_group_data = per_group_data;
1765 <<crosslink groups>>=
1766 void crosslink_groups(list_t *state_groups, list_t *transition_list)
1768 list_t *list, *out_trans=NULL, *associates=NULL;
1769 one_dim_fn_t *handler;
1772 state_groups = head(state_groups);
1773 while (state_groups != NULL) {
1774 /* find transitions leaving this state */
1776 fprintf(stderr, "finding transitions leaving %s\n", S(state_groups)->name);
1778 list = head(transition_list);
1779 while (list != NULL) {
1780 if (T(list)->initial_state == S(state_groups) && list_index(associates, T(list)->initial_state) != -1) {
1781 push(&out_trans, T(list), 1);
1785 n = list_length(out_trans);
1786 S(state_groups)->num_out_trans = n;
1787 S(state_groups)->out_trans = (transition_t **)malloc(sizeof(transition_t *)*n);
1788 assert(S(state_groups)->out_trans != NULL);
1789 for (i=0; i<n; i++) {
1790 S(state_groups)->out_trans[i] = (transition_t *)pop(&out_trans);
1793 /* find states grouped with this state */
1795 fprintf(stderr, "finding states grouped with %s\n", S(state_groups)->name);
1797 handler = S(state_groups)->tension_handler;
1798 group = S(state_groups)->tension_group;
1799 list = head(state_groups);
1800 while (list != NULL) {
1801 if (S(list)->tension_handler == handler && S(list)->tension_group == group && list != state_groups) {
1802 push(&associates, S(list), 1);
1806 n = list_length(associates);
1807 S(state_groups)->num_grouped_states = n;
1808 S(state_groups)->grouped_states = (state_t **)malloc(sizeof(state_t *)*n);
1809 assert(S(state_groups)->grouped_states != NULL);
1810 for (i=0; i<n; i++) {
1811 S(state_groups)->grouped_states[i] = (state_t *)pop(&associates);
1813 state_groups = state_groups->next;
1819 \section{String parsing}
1821 For handling command line arguments, we need parse delimited strings (e.g.\ [["param1,param2,param3..."]]).
1822 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1823 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1824 need to take care of parsing those parameters themselves.
1825 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1831 <<parse definitions>>
1832 <<parse declarations>>
1836 <<parse module makefile lines>>=
1837 NW_SAWSIM_MODS += parse
1838 CHECK_BINS += check_parse
1842 <<parse definitions>>=
1843 #define SEP ',' /* argument separator character */
1846 <<parse declarations>>=
1847 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1848 int *num, char ***string_array);
1849 extern void free_parsed_list(int num, char **string_array);
1852 [[parse_list_string]] allocates memory, don't forget to free it
1853 afterward with [[free_parsed_list]]. It does not alter the original.
1855 The string may start off with a [[deeper]] character
1856 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1857 [[x,y]], where the pointer is one character in on the copied string.
1858 However, when we go to free the memory, we need a pointer to the
1859 beginning of the string. In order to accommodate this for a string
1860 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1861 the first $N$ elements point to the separated arguments, and let the
1862 last element point to the start of the copied string regardless of
1864 <<parse delimited list string functions>>=
1865 /* TODO, split out into parse.hc */
1866 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1869 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1870 if (string[i] == deeper) {depth++;}
1871 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1877 void parse_list_string(char *string, char sep, char deeper, char shallower,
1878 int *num, char ***string_array)
1880 char *str=NULL, **ret=NULL;
1882 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
1884 *string_array = NULL;
1887 /* make a copy of the string, so we don't change the original */
1888 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1889 assert(str != NULL);
1890 strcpy(str, string); /* we know str is long enough */
1891 /* count the number of regions, so we can allocate pointers to them */
1894 n++; i++; /* move on to next argument */
1895 i += next_delim_index(str+i, sep, deeper, shallower);
1896 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1897 } while (str[i] != '\0');
1898 ret = (char **)malloc(sizeof(char *)*(n+1));
1899 assert(ret != NULL);
1900 /* replace the separators with '\0' & assign pointers */
1901 ret[n] = str; /* point to the front of the copied string */
1904 for(i=1; i<n; i++) {
1905 j += next_delim_index(str+j, sep, deeper, shallower);
1907 ret[i] = str+j; /* point to the separated arguments */
1909 /* strip off bounding braces, if any */
1910 for(i=0; i<n; i++) {
1911 if (ret[i][0] == deeper) {
1915 if (ret[i][strlen(ret[i])-1] == shallower)
1916 ret[i][strlen(ret[i])-1] = '\0';
1919 *string_array = ret;
1922 void free_parsed_list(int num, char **string_array)
1925 /* frees all items, since they were allocated as a single string*/
1926 free(string_array[num]);
1927 /* and free the array of pointers */
1933 \subsection{Parsing implementation}
1939 <<parse delimited list string functions>>
1943 #include <assert.h> /* assert() */
1944 #include <stdlib.h> /* NULL */
1945 #include <stdio.h> /* fprintf(), stdout *//*!!*/
1946 #include <string.h> /* strlen() */
1950 \subsection{Parsing unit tests}
1952 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1954 <<parse check includes>>
1955 <<string comparison definition and globals>>
1956 <<check relative error>>
1957 <<parse test suite>>
1958 <<main check program>>
1961 <<parse check includes>>=
1962 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
1963 #include <stdio.h> /* printf() */
1964 #include <assert.h> /* assert() */
1965 #include <string.h> /* strlen() */
1970 <<parse test suite>>=
1971 <<parse list string tests>>
1972 <<parse suite definition>>
1975 <<parse suite definition>>=
1976 Suite *test_suite (void)
1978 Suite *s = suite_create ("k model");
1979 <<parse list string test case defs>>
1981 <<parse list string test case add>>
1986 <<parse list string tests>>=
1989 START_TEST(test_next_delim_index)
1991 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
1992 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
1993 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
1994 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
1995 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
1999 START_TEST(test_parse_list_null)
2003 parse_list_string(NULL, SEP, '{', '}',
2004 &num_param_args, ¶m_args);
2005 fail_unless(num_param_args == 0, NULL);
2006 fail_unless(param_args == NULL, NULL);
2009 START_TEST(test_parse_list_single_simple)
2013 parse_list_string("arg", SEP, '{', '}',
2014 &num_param_args, ¶m_args);
2015 fail_unless(num_param_args == 1, NULL);
2016 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
2019 START_TEST(test_parse_list_single_compound)
2023 parse_list_string("{x,y,z}", SEP, '{', '}',
2024 &num_param_args, ¶m_args);
2025 fail_unless(num_param_args == 1, NULL);
2026 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
2029 START_TEST(test_parse_list_double_simple)
2033 parse_list_string("abc,def", SEP, '{', '}',
2034 &num_param_args, ¶m_args);
2035 fail_unless(num_param_args == 2, NULL);
2036 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2037 fail_unless(STRMATCH(param_args[1],"def"), NULL);
2040 START_TEST(test_parse_list_compound)
2044 parse_list_string("abc,{def,ghi}", SEP, '{', '}',
2045 &num_param_args, ¶m_args);
2046 fail_unless(num_param_args == 2, NULL);
2047 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
2048 fail_unless(STRMATCH(param_args[1],"def,ghi"), NULL);
2052 <<parse list string test case defs>>=
2053 TCase *tc_parse_list_string = tcase_create("parse list string");
2055 <<parse list string test case add>>=
2056 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
2057 tcase_add_test(tc_parse_list_string, test_parse_list_null);
2058 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
2059 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
2060 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
2061 tcase_add_test(tc_parse_list_string, test_parse_list_compound);
2062 suite_add_tcase(s, tc_parse_list_string);
2066 \section{Unit tests}
2068 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2075 <<check relative error>>
2078 <<main check program>>
2090 <<determine dt tests>>
2092 <<does domain transition tests>>
2093 <<randomly unfold domains tests>>
2094 <<suite definition>>
2097 <<suite definition>>=
2098 Suite *test_suite (void)
2100 Suite *s = suite_create ("sawsim");
2101 <<F test case defs>>
2102 <<determine dt test case defs>>
2103 <<happens test case defs>>
2104 <<does domain transition test case defs>>
2105 <<randomly unfold domains test case defs>>
2108 <<determine dt test case add>>
2109 <<happens test case add>>
2110 <<does domain transition test case add>>
2111 <<randomly unfold domains test case add>>
2114 tcase_add_checked_fixture(tc_strip_address,
2115 setup_strip_address,
2116 teardown_strip_address);
2122 <<main check program>>=
2126 Suite *s = test_suite();
2127 SRunner *sr = srunner_create(s);
2128 srunner_run_all(sr, CK_ENV);
2129 number_failed = srunner_ntests_failed(sr);
2131 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
2135 \subsection{F tests}
2137 <<worm-like chain tests>>
2139 <<F test case defs>>=
2140 <<worm-like chain test case def>>
2142 <<F test case add>>=
2143 <<worm-like chain test case add>>
2146 <<worm-like chain tests>>=
2147 <<worm-like chain function>>
2149 START_TEST(test_wlc_at_zero)
2151 double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
2152 fail_unless(abs(wlc(x, T, p, L)) < lim, \
2153 "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
2154 x, T, p, L, abs(wlc(x, T, p, L)), lim);
2158 START_TEST(test_wlc_at_half)
2160 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
2161 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
2162 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
2164 * linear term = x/L = 0.5
2165 * nonlinear + linear = 0.75 + 0.5 = 1.25
2166 * wlc = 10e21*1.25 = 12.5e21
2168 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
2169 "wlc(%g, %g, %g, %g) = %g != %g",
2170 x, T, p, L, wlc(x, T, p, L), 12.5e21);
2174 <<worm-like chain test case def>>=
2175 TCase *tc_wlc = tcase_create("WLC");
2178 <<worm-like chain test case add>>=
2179 tcase_add_test(tc_wlc, test_wlc_at_zero);
2180 tcase_add_test(tc_wlc, test_wlc_at_half);
2181 suite_add_tcase(s, tc_wlc);
2184 \subsection{Model tests}
2186 Check the searching with [[linear_k]].
2187 Check overwhelming force treatment with the heavyside-esque step [[?]].
2189 Hmm, I'm not really sure what I was doing here...
2190 <<determine dt tests>>=
2191 double linear_k(double F, environment_t *env, void *params)
2193 double Fnot = *(double *)params;
2197 START_TEST(test_determine_dt_linear_k)
2200 transition_t unfold={0};
2201 environment_t env = {0};
2202 double F[]={0,1,2,3,4,5,6};
2203 double dt_max=3.0, Fnot=3.0, x=1.0, targ_p=0.1, v=1.0;
2204 list_t *state_list=NULL, *transition_list=NULL;
2207 folded.tension_handler = &hooke_handler;
2208 folded.num_domains = 1;
2209 unfold.initial_state = &folded;
2210 unfold.k = linear_k;
2211 unfold.k_params = &Fnot;
2212 push(&state_list, &folded, 1);
2213 push(&transition_list, &unfold, 1);
2214 for( i=0; i < sizeof(F)/sizeof(double); i++) {
2215 //fail_unless(determine_dt(state_list, transition_list, &env, x, targ_p, dt_max, v)==1, NULL);
2220 <<determine dt test case defs>>=
2221 TCase *tc_determine_dt = tcase_create("Determine dt");
2223 <<determine dt test case add>>=
2224 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
2225 suite_add_tcase(s, tc_determine_dt);
2230 <<happens test case defs>>=
2232 <<happens test case add>>=
2235 <<does domain transition tests>>=
2237 <<does domain transition test case defs>>=
2239 <<does domain transition test case add>>=
2242 <<randomly unfold domains tests>>=
2244 <<randomly unfold domains test case defs>>=
2246 <<randomly unfold domains test case add>>=
2250 \section{Balancing group extensions}
2251 \label{sec.tension_balance}
2253 For a given total extension $x$ (of the piezo), the various domain
2254 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
2255 amounts, and we need to tweak the portion that each extends to balance
2256 the tension amongst the active groups. Since the tension balancing is
2257 separable from the bulk of the simulation, we break it out into a
2258 separate module. The interface is defined in [[tension_balance.h]],
2259 the implementation in [[tension_balance.c]], and the unit testing in
2260 [[check_tension_balance.c]]
2262 <<tension-balance.h>>=
2263 #ifndef TENSION_BALANCE_H
2264 #define TENSION_BALANCE_H
2266 <<tension balance function declaration>>
2267 #endif /* TENSION_BALANCE_H */
2270 <<tension balance functions>>=
2271 <<one dimensional bracket>>
2272 <<one dimensional solve>>
2274 <<group stiffness function>>
2275 <<chain stiffness function>>
2276 <<tension balance function>>
2279 <<tension balance module makefile lines>>=
2280 NW_SAWSIM_MODS += tension_balance
2281 CHECK_BINS += check_tension_balance
2282 check_tension_balance_MODS =
2285 The entire force balancing problem reduces to a solving two nested
2286 one-dimensional root-finding problems. First we define the one
2287 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
2288 populated group). $x(x_0)$ is also strictly monotonically increasing,
2289 so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
2290 stores the last successful value of $x$, since we expect to be taking
2291 small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
2292 indicates that the groups have changed.
2293 <<tension balance function declaration>>=
2294 double tension_balance(int num_tension_groups,
2295 one_dim_fn_t **tension_handlers,
2296 one_dim_fn_t **inverse_tension_handlers,
2297 void **params, /* array of pointers to tension_handler_data_t */
2302 <<tension balance function>>=
2303 double tension_balance(int num_tension_groups,
2304 one_dim_fn_t **tension_handlers,
2305 one_dim_fn_t **inverse_tension_handlers,
2306 void **params, /* array of pointers to tension_handler_data_t */
2311 static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL,NULL};
2312 double F, xo_guess, xo, lb, ub;
2313 double min_relx=1e-6, min_rely=1e-6;
2314 int max_steps=200, i;
2316 assert(num_tension_groups > 0);
2317 assert(tension_handlers != NULL);
2318 assert(params != NULL);
2319 assert(num_tension_groups > 0);
2321 if (last_x == -1 || x_xo_data.xi == NULL) { /* new group setup, reset x_xo_data */
2322 /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
2323 if (x_xo_data.xi != NULL)
2325 /* construct new x_xo_data */
2326 x_xo_data.n_groups = num_tension_groups;
2327 x_xo_data.tension_handler = tension_handlers;
2328 x_xo_data.inverse_tension_handler = inverse_tension_handlers;
2329 x_xo_data.handler_data = params;
2331 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);
2333 x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
2334 assert(x_xo_data.xi != NULL);
2335 for (i=0; i<num_tension_groups; i++)
2336 x_xo_data.xi[i] = last_x;
2339 if (num_tension_groups == 1) { /* only one group, no balancing required */
2341 x_xo_data.xi[0] = xo;
2343 fprintf(stderr, "%s:\tonly one active group, no balancing required\n", __FUNCTION__);
2347 fprintf(stderr, "%s:\tbalancing for x = %g with", __FUNCTION__, x);
2348 for(i=0; i<num_tension_groups; i++) fprintf(stderr, "\tx[%d]=%g", i, x_xo_data.xi[i]);
2349 fprintf(stderr, "\n");
2351 if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
2352 if (x == 0) xo_guess = length_scale;
2353 else xo_guess = x/num_tension_groups;
2355 fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
2357 oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
2358 } else { /* work off the last balanced point */
2360 fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
2363 lb = x_xo_data.xi[0];
2364 ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */
2365 } else if (x < last_x) {
2366 lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
2367 ub = x_xo_data.xi[0];
2368 } else { /* x == last_x */
2369 fprintf(stderr, "not moving\n");
2370 lb= x_xo_data.xi[0];
2371 ub= x_xo_data.xi[0];
2375 fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
2376 __FUNCTION__, x, lb, ub);
2378 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
2379 F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
2381 if (num_tension_groups > 1) {
2382 fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
2383 for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
2384 fprintf(stderr, "\n");
2389 F = (*tension_handlers[0])(xo, params[0]);
2391 if (stiffness != NULL)
2392 *stiffness = chain_stiffness(&x_xo_data, min_relx);
2399 <<tension balance internal definitions>>=
2400 <<x of x0 definitions>>
2403 <<x of x0 definitions>>=
2404 typedef struct x_of_xo_data_struct {
2406 one_dim_fn_t **tension_handler; /* array of fn pointers */
2407 one_dim_fn_t **inverse_tension_handler; /* inverse fn pointers */
2408 void **handler_data; /* array of pointers to tension_handler_data_t */
2409 double *xi; /* array of group extensions */
2414 double x_of_xo(double xo, void *pdata)
2416 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
2417 double F, x=0, xi, xi_guess, lb, ub, slope;
2419 double min_relx=1e-6, min_rely=1e-6;
2421 assert(data->n_groups > 0);
2424 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);
2427 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
2429 fprintf(stderr, "%s:\tfinding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
2433 if (data->xi) data->xi[0] = xo;
2434 for (i=1; i < data->n_groups; i++) {
2435 if (data->inverse_tension_handler[i] != NULL) {
2436 xi = (*data->inverse_tension_handler[i])(F, data->handler_data[i]);
2437 } else { /* invert numerically */
2438 if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
2439 else xi_guess = data->xi[i];
2441 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]);
2443 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
2444 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
2449 fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
2453 fprintf(stderr, "%s:\tfound x(x0=%g) = %g\n", __FUNCTION__, xo, x);
2459 <<group stiffness function>>=
2460 double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double x, double relx)
2462 double dx, xl, F, Fl, stiffness;
2464 assert(relx > 0 && relx < 1);
2465 if (x == 0 || relx == 0) {
2466 dx = 1e-12; /* HACK, how to get length scale? */
2476 F = (*tension_handler)(x, handler_data);
2477 Fl = (*tension_handler)(xl, handler_data);
2478 stiffness = (F-Fl)/dx;
2479 assert(stiffness > 0);
2484 <<chain stiffness function>>=
2485 double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx)
2489 /* add group stiffnesses in series */
2490 for (i=0; i < x_xo_data->n_groups; i++) {
2492 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);
2495 stiffness += 1.0/group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx);
2497 stiffness = 1.0 / stiffness;
2503 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
2504 number of steps for monotonically increasing $f(x)$. Simple
2505 bisection, so it's robust and converges fairly quickly. You can set
2506 the maximum number of steps to take, but if you have enough processor
2507 power, it's better to limit your precision with the relative limits
2508 [[min_relx]] and [[y]]. These ensure that the width of the bracket is
2509 small on the length scale of it's larger side. Note that \emph{both}
2510 relative limits must be satisfied for the search to stop.
2511 <<one dimensional function definition>>=
2512 /* equivalent to gsl_func_t */
2513 typedef double one_dim_fn_t(double x, void *params);
2515 <<one dimensional solve declaration>>=
2516 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2517 double min_relx, double min_rely, int max_steps,
2521 <<one dimensional solve>>=
2522 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2523 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2524 double min_relx, double min_rely, int max_steps,
2527 double yx, ylb, yub, x;
2530 ylb = (*f)(lb, params);
2531 yub = (*f)(ub, params);
2533 /* check some simple cases */
2535 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bracketed */
2536 else return lb; /* any x on the range [lb, ub] would work */
2538 if (ylb == y) { x=lb; yx=ylb; goto end; }
2539 if (yub == y) { x=ub; yx=yub; goto end; }
2542 fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2544 assert(ylb < y); assert(yub > y);
2545 for (i=0; i<max_steps; i++) {
2547 yx = (*f)(x, params);
2549 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);
2551 if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2552 else if (yx > y) { ub=x; yub=yx; }
2553 else /* yx < y */ { lb=x; ylb=yx; }
2558 fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2564 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2565 Generate bracketing $x$ values through bisection or exponential growth.
2566 <<one dimensional bracket declaration>>=
2567 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2570 <<one dimensional bracket>>=
2571 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2573 double yx, step, x=xguess;
2576 fprintf(stderr, "%s:\tbracketing %g (params %p)\n", __FUNCTION__, y, params);
2578 yx = (*f)(x, params);
2580 fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2587 fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2591 if (x == 0) x = length_scale/2.0; /* HACK */
2594 yx = (*f)(x, params);
2596 fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2601 fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2604 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2608 \subsection{Balancing implementation}
2610 <<tension-balance.c>>=
2613 <<tension balance includes>>
2614 #include "tension_balance.h"
2616 double length_scale = 1e-10; /* HACK */
2618 <<tension balance internal definitions>>
2619 <<tension balance functions>>
2622 <<tension balance includes>>=
2623 #include <assert.h> /* assert() */
2624 #include <stdlib.h> /* NULL */
2625 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2626 #include <stdio.h> /* fprintf(), stdout */
2630 \subsection{Balancing unit tests}
2632 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2633 <<check-tension-balance.c>>=
2634 <<tension balance check includes>>
2635 <<tension balance test suite>>
2636 <<main check program>>
2639 <<tension balance check includes>>=
2640 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2641 #include <assert.h> /* assert() */
2644 #include "tension_balance.h"
2647 <<tension balance test suite>>=
2648 <<tension balance function tests>>
2649 <<tension balance suite definition>>
2652 <<tension balance suite definition>>=
2653 Suite *test_suite (void)
2655 Suite *s = suite_create ("tension balance");
2656 <<tension balance function test case defs>>
2658 <<tension balance function test case adds>>
2663 <<tension balance function tests>>=
2664 <<check relative error>>
2666 double hooke(double x, void *pK)
2669 return *((double*)pK) * x;
2672 START_TEST(test_single_function)
2674 double k=5, x=3, last_x=2.0, F;
2675 one_dim_fn_t *handlers[] = {&hooke};
2676 one_dim_fn_t *inverse_handlers[] = {NULL};
2677 void *data[] = {&k};
2678 F = tension_balance(1, handlers, inverse_handlers, data, last_x, x, NULL);
2679 fail_unless(F = k*x, NULL);
2684 We can also test balancing two springs with different spring constants.
2685 The analytic solution for a general number of springs is given in Appendix \ref{sec.math_hooke}.
2686 <<tension balance function tests>>=
2687 START_TEST(test_double_hooke)
2689 double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
2690 one_dim_fn_t *handlers[] = {&hooke, &hooke};
2691 one_dim_fn_t *inverse_handlers[] = {NULL, NULL};
2692 void *data[] = {&k1, &k2};
2693 F = tension_balance(2, handlers, inverse_handlers, data, last_x, x, NULL);
2697 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2698 //CHECK_ERR(1e-6, x1e, xi[0]);
2699 //CHECK_ERR(1e-6, x2e, xi[1]);
2700 CHECK_ERR(1e-6, Fe, F);
2704 TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
2706 <<tension balance function test case defs>>=
2707 TCase *tc_tbfunc = tcase_create("tension balance function");
2710 <<tension balance function test case adds>>=
2711 tcase_add_test(tc_tbfunc, test_single_function);
2712 tcase_add_test(tc_tbfunc, test_double_hooke);
2713 suite_add_tcase(s, tc_tbfunc);
2719 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2720 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2721 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2727 <<list definitions>>
2728 <<list declarations>>
2732 <<list declarations>>=
2733 <<head and tail declarations>>
2734 <<list length declaration>>
2735 <<push declaration>>
2737 <<list destroy declaration>>
2738 <<list to array declaration>>
2739 <<move declaration>>
2740 <<sort declaration>>
2741 <<list index declaration>>
2742 <<list element declaration>>
2743 <<unique declaration>>
2744 <<list print declaration>>
2748 <<create and destroy node>>
2763 <<list module makefile lines>>=
2764 NW_SAWSIM_MODS += list
2765 CHECK_BINS += check_list
2769 <<list definitions>>=
2770 typedef struct list_struct {
2771 struct list_struct *next;
2772 struct list_struct *prev;
2776 <<comparison function typedef>>
2779 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2780 <<head and tail declarations>>=
2781 list_t *head(list_t *list);
2782 list_t *tail(list_t *list);
2785 list_t *head(list_t *list)
2789 while (list->prev != NULL) {
2795 list_t *tail(list_t *list)
2799 while (list->next != NULL) {
2806 <<list length declaration>>=
2807 int list_length(list_t *list);
2810 int list_length(list_t *list)
2817 while (list->next != NULL) {
2825 [[push]] inserts a node relative to the current position in the list
2826 without changing the current position.
2827 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2828 so we set it to that of the pushed domain.
2829 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2830 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2831 <<push declaration>>=
2832 void push(list_t **pList, void *data, int below);
2835 void push(list_t **pList, void *data, int below)
2837 list_t *list, *new_node;
2838 assert(pList != NULL);
2840 new_node = create_node(data);
2843 else if (below==0) { /* insert above */
2844 if (list->prev != NULL)
2845 list->prev->next = new_node;
2846 new_node->prev = list->prev;
2847 list->prev = new_node;
2848 new_node->next = list;
2849 } else { /* insert below */
2850 if (list->next != NULL)
2851 list->next->prev = new_node;
2852 new_node->next = list->next;
2853 list->next = new_node;
2854 new_node->prev = list;
2859 [[pop]] removes the current domain node, moving the current position
2860 to the node after it, unless that node is [[NULL]], in which case move
2861 the current position to the node before it.
2862 <<pop declaration>>=
2863 void *pop(list_t **pList);
2866 void *pop(list_t **pList)
2868 list_t *list, *popped;
2870 assert(pList != NULL);
2872 assert(list != NULL); /* not an empty list */
2874 /* bypass the popped node */
2875 if (list->prev != NULL)
2876 list->prev->next = list->next;
2877 if (list->next != NULL) {
2878 list->next->prev = list->prev;
2879 *pList = list->next;
2881 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
2883 destroy_node(popped);
2888 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
2889 If the cleanup function is [[NULL]], no cleanup function is called.
2890 The nodes are not popped in any particular order.
2891 <<list destroy declaration>>=
2892 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
2895 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
2899 assert(pList != NULL);
2902 assert(list != NULL); /* not an empty list */
2903 while (list != NULL) {
2905 if (destroy != NULL)
2911 [[list_to_array]] converts a list to an array of pointers, preserving order.
2912 <<list to array declaration>>=
2913 void list_to_array(list_t *list, int *length, void ***array);
2916 void list_to_array(list_t *list, int *pLength, void ***pArray)
2920 assert(list != NULL);
2921 assert(pLength != NULL);
2922 assert(pArray != NULL);
2924 length = list_length(list);
2925 array = (void **)malloc(sizeof(void *)*length);
2926 assert(array != NULL);
2929 while (list != NULL) {
2930 array[i++] = list->d;
2938 [[move]] moves the current element along the list in either direction.
2939 <<move declaration>>=
2940 void move(list_t **pList, int downstream);
2943 void move(list_t **pList, int downstream)
2945 list_t *list, *flipper;
2947 assert(pList != NULL);
2948 list = *pList; /* a>B>c>d */
2949 assert(list != NULL); /* not an empty list */
2950 if (downstream == 0)
2951 flipper = list->prev; /* flipper is A */
2953 flipper = list->next; /* flipper is C */
2954 /* ensure there is somewhere to go in the direction we're heading */
2955 assert(flipper != NULL);
2956 /* remove the current element */
2957 data = pop(&list); /* data=B, a>C>d */
2958 /* and put it back in in the correct spot */
2960 if (downstream == 0) {
2961 push(&list, data, 0); /* b>A>c>d */
2962 list = list->prev; /* B>a>c>d */
2964 push(&list, data, 1); /* a>C>b>d */
2965 list = list->next; /* a>c>B>d */
2971 [[sort]] sorts a list from smallest to largest according to a user
2972 specified comparison.
2973 <<comparison function typedef>>=
2974 typedef int (comparison_fn_t)(void *A, void *B);
2977 <<int comparison function>>=
2978 int int_comparison(void *A, void *B)
2983 if (a > b) return 1;
2984 else if (a == b) return 0;
2989 <<sort declaration>>=
2990 void sort(list_t **pList, comparison_fn_t *comp);
2992 Here we use the bubble sort algorithm.
2994 void sort(list_t **pList, comparison_fn_t *comp)
2998 assert(pList != NULL);
3000 assert(list != NULL); /* not an empty list */
3001 while (stable == 0) {
3004 while (list->next != NULL) {
3005 if ((*comp)(list->d, list->next->d) > 0) {
3007 move(&list, 1 /* downstream */);
3017 [[list_index]] finds the location of [[data]] in [[list]]. Returns
3018 [[-1]] if [[data]] is not in [[list]].
3019 <<list index declaration>>=
3020 int list_index(list_t *list, void *data);
3024 int list_index(list_t *list, void *data)
3028 while (list != NULL) {
3029 if (list->d == data) return i;
3038 [[list_element]] returns the data bound to the $i^\text{th}$ node of the list.
3039 <<list element declaration>>=
3040 void *list_element(list_t *list, int i);
3043 void *list_element(list_t *list, int i)
3047 while (list != NULL) {
3048 if (i == j) return list->d;
3058 [[unique]] removes duplicates from a sorted list according to a user
3059 specified comparison. Don't do this unless you have other pointers
3060 any dynamically allocated data in your list, because [[unique]] just
3061 drops any non-unique data without freeing it.
3062 <<unique declaration>>=
3063 void unique(list_t **pList, comparison_fn_t *comp);
3066 void unique(list_t **pList, comparison_fn_t *comp)
3069 assert(pList != NULL);
3071 assert(list != NULL); /* not an empty list */
3073 while (list->next != NULL) {
3074 if ((*comp)(list->d, list->next->d) == 0) {
3083 [[list_print]] prints a list to a [[FILE *]] stream.
3084 <<list print declaration>>=
3085 void list_print(FILE *file, list_t *list, char *list_name);
3088 void list_print(FILE *file, list_t *list, char *list_name)
3090 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
3092 while (list != NULL) {
3093 fprintf(file, " %p", list->d);
3096 fprintf(file, "\n");
3100 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
3101 <<create and destroy node>>=
3102 list_t *create_node(void *data)
3104 list_t *ret = (list_t *)malloc(sizeof(list_t));
3105 assert(ret != NULL);
3112 void destroy_node(list_t *node)
3118 The user must free the data pointed to by the node on their own.
3120 \subsection{List implementation}
3130 #include <assert.h> /* assert() */
3131 #include <stdlib.h> /* malloc(), free(), rand() */
3132 #include <stdio.h> /* fprintf(), stdout, FILE */
3133 #include "global.h" /* destroy_data_func_t */
3136 \subsection{List unit tests}
3138 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3140 <<list check includes>>
3142 <<main check program>>
3145 <<list check includes>>=
3146 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3147 #include <stdio.h> /* FILE */
3153 <<list test suite>>=
3156 <<list suite definition>>
3159 <<list suite definition>>=
3160 Suite *test_suite (void)
3162 Suite *s = suite_create ("list");
3163 <<push test case defs>>
3164 <<pop test case defs>>
3166 <<push test case adds>>
3167 <<pop test case adds>>
3173 START_TEST(test_push)
3178 push(&list, (void *)i, 1);
3179 fail_unless(list == head(list), NULL);
3180 fail_unless( (int)list->d == 0 );
3181 for (i=0; i<n; i++) {
3182 /* we expect 0, n-1, n-2, ..., 1 */
3185 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
3191 <<push test case defs>>=
3192 TCase *tc_push = tcase_create("push");
3195 <<push test case adds>>=
3196 tcase_add_test(tc_push, test_push);
3197 suite_add_tcase(s, tc_push);
3202 <<pop test case defs>>=
3204 <<pop test case adds>>=
3207 \section{Function string evaluation}
3209 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).
3210 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
3211 The solution is to run a scripting language as a subprocess accessed via pipes.
3212 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
3214 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
3215 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
3216 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.
3217 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
3219 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.
3220 Then you can either statically or dynamically link to those hardcoded functions.
3221 While much less flexible, this approach would be a fairly simple-to-implement fallback.
3223 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
3224 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
3227 #ifndef STRING_EVAL_H
3228 #define STRING_EVAL_H
3230 <<string eval setup declaration>>
3231 <<string eval function declaration>>
3232 <<string eval teardown declaration>>
3233 #endif /* STRING_EVAL_H */
3236 <<string eval module makefile lines>>=
3237 NW_SAWSIM_MODS += string_eval
3238 CHECK_BINS += check_string_eval
3239 check_string_eval_MODS =
3242 For an introduction to POSIX process control, see\\
3243 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
3244 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
3245 and of course, the relavent [[man]] pages.
3247 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
3248 [[execvp]] replaces the calling process' program with a new program.
3249 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
3250 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
3251 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
3252 The new program needs command line arguments, just like it would if you were running it from a shell.
3253 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
3254 with the final array entry being a [[NULL]] pointer.
3256 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
3257 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
3258 <<string eval subprocess definitions>>=
3259 #define SUBPROCESS "python"
3260 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
3261 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
3262 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
3264 The [[i]] option lets Python know that it should run in interactive mode.
3265 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
3266 In interactive mode, python acts on every instruction as soon as it is recieved.
3267 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.
3268 %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.
3270 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
3271 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
3272 [[fork]] returns two copies of the same program, executing the original code.
3273 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
3274 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
3276 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.
3277 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
3278 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
3279 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
3280 <<string eval pipe definitions>>=
3281 #define PIPE_READ 0 /* the end you read from */
3282 #define PIPE_WRITE 1 /* the end you write to */
3284 #define STDIN 0 /* initial index of stdin pair */
3285 #define STDOUT 2 /* initial index of stdout pair */
3288 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
3290 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.
3292 <<string eval setup declaration>>=
3293 extern void string_eval_setup(FILE **pIn, FILE **pOut);
3295 <<string eval setup definition>>=
3296 void string_eval_setup(FILE **pIn, FILE **pOut)
3299 int pfd[4]; /* file descriptors for stdin and stdout */
3301 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
3302 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
3304 pid = fork(); /* split process into two copies */
3306 if (pid == -1) { /* fork error */
3307 perror("fork error");
3309 } else if (pid == 0) { /* child */
3310 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
3311 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
3312 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
3313 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
3314 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
3315 perror("exec error"); /* exec shouldn't return */
3317 } else { /* parent */
3318 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
3319 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
3320 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
3321 if ( *pIn == NULL ) {
3322 perror("fdopen (in)");
3325 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
3326 if ( *pOut == NULL ) {
3327 perror("fdopen (out)");
3334 To use the evaluating subprocess, we just pipe in our command, and read out the results.
3335 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
3336 <<string eval function declaration>>=
3337 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
3339 <<string eval function definition>>=
3340 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
3343 rc = fprintf(in, "%s", input);
3344 assert(rc == strlen(input));
3347 alarm(1); /* set a one second timeout on the read */
3348 assert( fgets(output, buflen, out) != NULL );
3349 alarm(0); /* cancel the timeout */
3350 //fprintf(stderr, "eval: %s ----> %s", input, output);
3353 The [[alarm]] calls set and clear a timeout on the returned output.
3354 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.
3355 This protects against invalid input for which a line of output is not printed to [[stdout]].
3356 Other invalid inputs (e.g.\ those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
3357 If you are getting strange results, check your python code seperately. TODO, better error handling.
3359 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
3360 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
3361 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
3362 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
3363 <<string eval teardown declaration>>=
3364 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
3367 <<string eval teardown definition>>=
3368 void string_eval_teardown(FILE **pIn, FILE **pOut)
3373 /* redirect Python's stderr */
3374 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
3378 assert( fclose(*pIn) == 0 );
3380 assert( fclose(*pOut) == 0 );
3383 /* wait for python to exit */
3385 pid = wait(&stat_loc);
3392 if (WIFEXITED(stat_loc)) {
3393 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
3394 } else if (WIFSIGNALED(stat_loc)) {
3395 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
3400 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
3402 \subsection{String evaluation implementation}
3406 <<string eval includes>>
3407 #include "string_eval.h"
3408 <<string eval internal definitions>>
3409 <<string eval functions>>
3412 <<string eval includes>>=
3413 #include <assert.h> /* assert() */
3414 #include <stdlib.h> /* NULL */
3415 #include <stdio.h> /* fprintf(), stdout, fdopen() */
3416 #include <string.h> /* strlen() */
3417 #include <sys/types.h> /* pid_t */
3418 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
3419 #include <wait.h> /* wait() */
3422 <<string eval internal definitions>>=
3423 <<string eval subprocess definitions>>
3424 <<string eval pipe definitions>>
3427 <<string eval functions>>=
3428 <<string eval setup definition>>
3429 <<string eval function definition>>
3430 <<string eval teardown definition>>
3433 \subsection{String evaluation unit tests}
3435 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3436 <<check-string-eval.c>>=
3437 <<string eval check includes>>
3438 <<string comparison definition and globals>>
3439 <<string eval test suite>>
3440 <<main check program>>
3443 <<string eval check includes>>=
3444 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3445 #include <stdio.h> /* printf() */
3446 #include <assert.h> /* assert() */
3447 #include <string.h> /* strlen() */
3448 #include <signal.h> /* SIGKILL */
3450 #include "string_eval.h"
3453 <<string eval test suite>>=
3454 <<string eval tests>>
3455 <<string eval suite definition>>
3458 <<string eval suite definition>>=
3459 Suite *test_suite (void)
3461 Suite *s = suite_create ("string eval");
3462 <<string eval test case defs>>
3464 <<string eval test case add>>
3469 <<string eval tests>>=
3470 START_TEST(test_setup_teardown)
3473 string_eval_setup(&in, &out);
3474 string_eval_teardown(&in, &out);
3477 START_TEST(test_invalid_command)
3480 char input[80], output[80]={};
3481 string_eval_setup(&in, &out);
3482 sprintf(input, "print ABCDefg\n");
3483 string_eval(in, out, input, 80, output);
3484 string_eval_teardown(&in, &out);
3487 START_TEST(test_simple_eval)
3490 char input[80], output[80]={};
3491 string_eval_setup(&in, &out);
3492 sprintf(input, "print 3+4*5\n");
3493 string_eval(in, out, input, 80, output);
3494 fail_unless(STRMATCH(output,"23\n"), NULL);
3495 string_eval_teardown(&in, &out);
3498 START_TEST(test_multiple_evals)
3501 char input[80], output[80]={};
3502 string_eval_setup(&in, &out);
3503 sprintf(input, "print 3+4*5\n");
3504 string_eval(in, out, input, 80, output);
3505 fail_unless(STRMATCH(output,"23\n"), NULL);
3506 sprintf(input, "print (3**2 + 4**2)**0.5\n");
3507 string_eval(in, out, input, 80, output);
3508 fail_unless(STRMATCH(output,"5.0\n"), NULL);
3509 string_eval_teardown(&in, &out);
3512 START_TEST(test_eval_with_state)
3515 char input[80], output[80]={};
3516 string_eval_setup(&in, &out);
3517 sprintf(input, "print 3+4*5\n");
3518 fprintf(in, "A = 3\n");
3519 sprintf(input, "print A*3\n");
3520 string_eval(in, out, input, 80, output);
3521 fail_unless(STRMATCH(output,"9\n"), NULL);
3522 string_eval_teardown(&in, &out);
3526 <<string eval test case defs>>=
3527 TCase *tc_string_eval = tcase_create("string_eval");
3529 <<string eval test case add>>=
3530 tcase_add_test(tc_string_eval, test_setup_teardown);
3531 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
3532 tcase_add_test(tc_string_eval, test_simple_eval);
3533 tcase_add_test(tc_string_eval, test_multiple_evals);
3534 tcase_add_test(tc_string_eval, test_eval_with_state);
3535 suite_add_tcase(s, tc_string_eval);
3539 \section{Accelerating function evaluation}
3541 My first version-0.3 code was running very slowly.
3542 With the options suggested in the help
3543 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
3544 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
3545 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
3546 That's only 15 calls per solution, so the search algorithm seems reasonable.
3547 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
3552 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
3554 #endif /* ACCEL_K_H */
3557 <<accel k module makefile lines>>=
3558 NW_SAWSIM_MODS += accel_k
3559 #CHECK_BINS += check_accel_k
3560 check_accel_k_MODS =
3564 #include <assert.h> /* assert() */
3565 #include <stdlib.h> /* realloc(), free(), NULL */
3566 #include "global.h" /* environment_t */
3567 #include "k_model.h" /* k_func_t */
3568 #include "interp.h" /* interp_* */
3569 #include "accel_k.h"
3571 typedef struct accel_k_struct {
3572 interp_table_t *itable;
3578 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3579 static int num_accels = 0;
3580 static accel_k_t *accels=NULL;
3582 /* Wrap k in the standard f(x) acceleration form */
3583 static double k_wrap(double F, void *params)
3585 accel_k_t *p = (accel_k_t *)params;
3586 return (*p->k)(F, p->env, p->params);
3589 static int k_tol(double FA, double kA, double FB, double kB)
3592 if (FB-FA > 1e-12) {
3593 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3594 return 1; /* unnacceptable */
3596 //printf("acceptable tol\n");
3597 return 0; /* acceptable */
3601 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3604 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3605 assert(accels != NULL);
3606 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3608 accels[i].env = env;
3609 accels[i].params = params;
3616 for (i=0; i<num_accels; i++)
3617 interp_table_free(accels[i].itable);
3619 if (accels) free(accels);
3623 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3626 for (i=0; i<num_accels; i++) {
3627 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3628 /* We've already setup this function.
3629 * WARNING: we're assuming param and env are the same. */
3630 return interp_table_eval(accels[i].itable, accels+i, F);
3633 /* set up a new acceleration environment */
3634 i = add_accel_k(k, env, params);
3635 return interp_table_eval(accels[i].itable, accels+i, F);
3639 \section{Tension models}
3640 \label{sec.tension_models}
3642 TODO: tension model intro.
3643 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.
3645 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3646 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]].
3648 <<tension-model.h>>=
3649 #ifndef TENSION_MODEL_H
3650 #define TENSION_MODEL_H
3652 <<tension handler types>>
3653 <<constant tension model declarations>>
3654 <<hooke tension model declarations>>
3655 <<worm-like chain tension model declarations>>
3656 <<freely-jointed chain tension model declarations>>
3657 <<find tension definitions>>
3658 <<tension model global declarations>>
3659 #endif /* TENSION_MODEL_H */
3662 <<tension model module makefile lines>>=
3663 NW_SAWSIM_MODS += tension_model
3664 #CHECK_BINS += check_tension_model
3665 check_tension_model_MODS =
3667 <<tension model utils makefile lines>>=
3668 TENSION_MODEL_MODS = tension_model parse list tension_balance
3669 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3670 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3671 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3672 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
3673 notangle -Rtension-model-utils.c $< > $@
3674 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
3675 gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3676 clean_tension_model_utils :
3677 rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3678 @ The pipe symbol [[|]] marks the prerequisits following it (in this
3679 case, the directories) as ‘order-only’ prerequisites. The timestamp
3680 on these prerequisits does not effect whether the rules are executed.
3681 This is appropriate for directories, where we don't need to recompile
3682 after an unrelated has been added to the directory, but only when the
3683 source prerequisites change. See the [[make]] documentation for more
3685 (\url{http://www.gnu.org/software/make/manual/html_node/Prerequisite-Types.html}).
3688 \label{sec.null_tension}
3690 For unstretchable domains.
3692 <<null tension model getopt>>=
3693 {NULL, NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3696 \subsection{Constant}
3697 \label{sec.const_tension}
3699 An infinitely stretchable domain providing a constant tension.
3700 Obviously this cannot be inverted, so you can't put this domain in
3701 series with any other domains. We include it mostly for testing
3704 <<constant tension functions>>=
3705 <<constant tension function>>
3706 <<constant tension parameter create/destroy functions>>
3709 <<constant tension model declarations>>=
3710 <<constant tension function declaration>>
3711 <<constant tension parameter create/destroy function declarations>>
3712 <<constant tension model global declarations>>
3716 <<constant tension function declaration>>=
3717 extern double const_tension_handler(double x, void *pdata);
3719 <<constant tension function>>=
3720 double const_tension_handler(double x, void *pdata)
3722 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3723 list_t *list = data->group_tension_params;
3728 assert(list != NULL); /* empty active group?! */
3729 F = ((const_tension_param_t *)list->d)->F;
3731 fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3733 while (list != NULL) {
3734 assert(((const_tension_param_t *)list->d)->F == F);
3742 <<constant tension structure>>=
3743 typedef struct const_tension_param_struct {
3744 double F; /* tension (force) in N */
3745 } const_tension_param_t;
3749 <<constant tension parameter create/destroy function declarations>>=
3750 extern void *string_create_const_tension_param_t(char **param_strings);
3751 extern void destroy_const_tension_param_t(void *p);
3754 <<constant tension parameter create/destroy functions>>=
3755 const_tension_param_t *create_const_tension_param_t(double F)
3757 const_tension_param_t *ret
3758 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3759 assert(ret != NULL);
3764 void *string_create_const_tension_param_t(char **param_strings)
3766 return create_const_tension_param_t(
3767 safe_strtod(param_strings[0],__FUNCTION__));
3770 void destroy_const_tension_param_t(void *p)
3778 <<constant tension model global declarations>>=
3779 extern int num_const_tension_params;
3780 extern char *const_tension_param_descriptions[];
3781 extern char const_tension_param_string[];
3784 <<constant tension model globals>>=
3785 int num_const_tension_params = 1;
3786 char *const_tension_param_descriptions[] = {"tension F, N"};
3787 char const_tension_param_string[] = "0";
3790 <<constant tension model getopt>>=
3791 {&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}
3797 <<hooke functions>>=
3798 <<internal hooke functions>>
3800 <<inverse hooke handler>>
3801 <<hooke parameter create/destroy functions>>
3804 <<hooke tension model declarations>>=
3805 <<hooke tension function declaration>>
3806 <<hooke tension parameter create/destroy function declarations>>
3807 <<hooke tension model global declarations>>
3811 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
3812 The behavior of a series of springs $k_i$ in series is given by
3814 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
3816 For a simple proof, see Appendix \ref{sec.math_hooke}.
3818 <<hooke tension function declaration>>=
3819 extern double hooke_handler(double x, void *pdata);
3820 extern double inverse_hooke_handler(double F, void *pdata);
3823 First we define a function that computes the effective spring constant
3824 (stored in a single [[hooke_param_t]]) for the entire group.
3825 <<internal hooke functions>>=
3826 static void hooke_param_grouper(tension_handler_data_t *data,
3827 hooke_param_t *param) {
3828 list_t *list = data->group_tension_params;
3832 assert(list != NULL); /* empty active group?! */
3833 while (list != NULL) {
3834 assert( ((hooke_param_t *)list->d)->k > 0 );
3835 k += 1.0/ ((hooke_param_t *)list->d)->k;
3837 fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
3838 __FUNCTION__, ((hooke_param_t *)list->d)->k);
3844 fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N (from data %p)\n",
3845 __FUNCTION__, k, x, k*x, data);
3852 Everyone knows Hooke's law: $F=kx$.
3854 double hooke_handler(double x, void *pdata)
3856 hooke_param_t param = {0};
3859 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
3865 Inverting Hooke's law gives $x=F/k$.
3866 <<inverse hooke handler>>=
3867 double inverse_hooke_handler(double F, void *pdata)
3869 hooke_param_t param = {0};
3872 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
3878 Not much to keep track of for a single effective spring, just the
3879 spring constant $k$.
3880 <<hooke structure>>=
3881 typedef struct hooke_param_struct {
3882 double k; /* spring constant in N/m */
3887 We wrap up our Hooke implementation with some book-keeping functions.
3888 <<hooke tension parameter create/destroy function declarations>>=
3889 extern void *string_create_hooke_param_t(char **param_strings);
3890 extern void destroy_hooke_param_t(void *p);
3894 <<hooke parameter create/destroy functions>>=
3895 hooke_param_t *create_hooke_param_t(double k)
3897 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
3898 assert(ret != NULL);
3903 void *string_create_hooke_param_t(char **param_strings)
3905 return create_hooke_param_t(safe_strtod(param_strings[0], __FUNCTION__));
3908 void destroy_hooke_param_t(void *p)
3915 <<hooke tension model global declarations>>=
3916 extern int num_hooke_params;
3917 extern char *hooke_param_descriptions[];
3918 extern char hooke_param_string[];
3921 <<hooke tension model globals>>=
3922 int num_hooke_params = 1;
3923 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
3924 char hooke_param_string[]="0.05";
3927 <<hooke tension model getopt>>=
3928 {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}
3931 \subsection{Worm-like chain}
3934 We can model several unfolded domains with a single worm-like chain.
3935 <<worm-like chain functions>>=
3936 <<internal worm-like chain functions>>
3937 <<worm-like chain handler>>
3938 <<inverse worm-like chain handler>>
3939 <<worm-like chain create/destroy functions>>
3942 <<worm-like chain tension model declarations>>=
3944 <<worm-like chain handler declaration>>
3945 <<inverse worm-like chain handler declaration>>
3946 <<worm-like chain create/destroy function declarations>>
3947 <<worm-like chain tension model global declarations>>
3951 <<internal worm-like chain functions>>=
3952 <<worm-like chain function>>
3953 <<inverse worm-like chain function>>
3954 <<worm-like chain parameter grouper>>
3957 The combination of all unfolded domains is modeled as a worm like chain,
3958 with the force $F_{\text{WLC}}$ approximately given by
3960 F_{\text{WLC}} \approx \frac{k_B T}{p}
3961 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
3963 where $x$ is the distance between the fixed ends,
3964 $k_B$ is Boltzmann's constant,
3965 $T$ is the absolute temperature,
3966 $p$ is the persistence length, and
3967 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
3968 <<worm-like chain function>>=
3969 static double wlc(double x, double T, double p, double L)
3971 double xL=0.0; /* uses global k_B */
3973 assert(T > 0); assert(p > 0); assert(L > 0);
3974 if (x >= L) return HUGE_VAL;
3975 xL = x/L; /* unitless */
3976 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
3977 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
3978 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
3983 Inverting the approximate WLC is fairly straightforward, if a bit tedious.
3984 Letting $x_L \equiv x/L$ and $F_T \equiv Fp/k_B T$ we have
3986 F_T &\approx \frac{1}{4} \p({\frac{1}{(1-x_L)^2} - 1}) + x_L \\
3987 F_T (1-x_L)^2 &\approx \frac{1}{4} - \frac{1}{4} (1-x_L)^2 + x_L(1-x_L)^2 \\
3988 0 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) (1-x_L)^2 + x_L(1-x_L)^2 \\
3989 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) \p({1-2x_L+x_L^2})
3990 + x_L - 2x_L^2 + x_L^3 \\
3991 &\approx \p[{\frac{1}{4} - \p({F_T + \frac{1}{4}})}]
3992 + x_L \p[{2F_T + \frac{1}{2} + 1}]
3993 + x_L^2 \p[{-\p({F_T + \frac{1}{4}}) - 2}]
3996 + x_L \p({2F_T + \frac{3}{2}})
3997 - x_L^2 \p({F_T + \frac{9}{4}})
3999 &\approx ax_L^3 + bx_L^2 + cx_L + d
4001 where $a \equiv 1$, $b \equiv -F_T - 9/4$, $c \equiv 2F_T + 3/2$, and
4003 % From \citet{wikipedia_cubic_function} the discriminant
4005 % \Delta &= -4b^3d + b^2c^2 - 4ac^3 + 18abcd - 27a^2d^2 \\
4006 % &= -4F_T\p({F_T + \frac{9}{4}})^3
4007 % + \p({F_T + \frac{9}{4}})^2 \cdot \p({2F_T + \frac{3}{2}})^2
4008 % - 4 \p({2F_T + \frac{3}{2}})^3
4009 % + 18F_T \p({F_T + \frac{9}{4}}) \cdot \p({2F_T + \frac{3}{2}})
4011 % &= -4F_T \p({F_T^3 + \frac{27}{4}F_T^2 +\frac{243}{16}F_T +\frac{729}{64}})
4012 %% a^3 + 3a^2b + 3ab^2 + b^3
4013 % + \p({F_T^2 + \frac{9}{2}F_T + \frac{81}{16}})
4014 % \cdot \p({4F_T^2 + 6F_T + \frac{9}{4}}) \\
4015 % &\qquad - 4 \p({8F_T^3 + 18F_T^2 + \frac{27}{2}F_T + \frac{27}{8}})
4016 %% a^3 + 3a^2b + 3ab^2 + b^3
4017 % + 18F_T \p({2F_T^2 + 6F_T + \frac{27}{8}})
4019 % &= -4F_T^4 - 27F_T^3 - \frac{243}{4}F_T^2 - \frac{729}{16}F_T
4020 % + 4F_T^4 + 24F_T^3 +\frac{99}{2}F_T^2 +\frac{81}{2}F_T +\frac{729}{64}\\
4021 % &\qquad - 32F_T^3 - 72F_T^2 - 54F_T - \frac{27}{2}
4022 % + 36F_T^3 + 108F_T^2 + \frac{243}{4}F_T
4024 % &= F_T^4 \p({-4 + 4}) + F_T^3 \p({-27 + 24 - 32 + 36})
4025 % + F_T^2 \p({-\frac{243}{4} + \frac{99}{2} - 72 + 108 -27}) \\
4026 % &\qquad + F_T \p({-\frac{729}{16} + \frac{81}{2} - 54 + \frac{243}{4}})
4027 % \p({\frac{729}{64} - \frac{27}{2}}) \\
4028 % &= F_T^3 - \frac{9}{4}F_T^2 + \frac{27}{16}F_T - \frac{135}{64}
4030 We can plug these coefficients into the GSL cubic solver to invert the
4031 WLC\citep{galassi05}. The applicable root is always the one which
4032 comes before the singularity, which will be the smallest real root.
4033 <<inverse worm-like chain function>>=
4034 static double inverse_wlc(double F, double T, double p, double L)
4036 double FT = F*p/(k_B*T); /* uses global k_B */
4037 double xL0, xL1, xL2;
4040 assert(T > 0); assert(p > 0); assert(L > 0);
4043 num_roots = gsl_poly_solve_cubic(-(FT+2.25),2*FT+1.5,-FT, &xL0, &xL1, &xL2);
4044 assert(num_roots > 0); assert(xL0 >= -DOUBLE_PRECISION); assert(xL0 < 1);
4045 if (xL0 < 0) xL0 = 0.0;
4051 First we define a function that computes the effective WLC parameters
4052 (stored in a single [[wlc_param_t]]) for the entire group.
4053 <<worm-like chain parameter grouper>>=
4054 static void wlc_param_grouper(tension_handler_data_t *data,
4055 wlc_param_t *param) {
4056 list_t *list = data->group_tension_params;
4057 double p=0.0, L=0.0;
4060 assert(list != NULL); /* empty active group?! */
4061 p = ((wlc_param_t *)list->d)->p;
4062 while (list != NULL) {
4063 /* enforce consistent p */
4064 assert( ((wlc_param_t *)list->d)->p == p);
4065 L += ((wlc_param_t *)list->d)->L;
4067 fprintf(stderr, "%s: WLC section %g m long in series\n",
4068 __FUNCTION__, ((wlc_param_t *)list->d)->L);
4073 fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
4074 __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
4082 <<worm-like chain handler declaration>>=
4083 extern double wlc_handler(double x, void *pdata);
4086 This model requires that each unfolded domain in the group have the
4087 same persistence length.
4088 <<worm-like chain handler>>=
4089 double wlc_handler(double x, void *pdata)
4091 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4092 wlc_param_t param = {0};
4094 wlc_param_grouper(data, ¶m);
4095 return wlc(x, data->env->T, param.p, param.L);
4100 <<inverse worm-like chain handler declaration>>=
4101 extern double inverse_wlc_handler(double F, void *pdata);
4104 <<inverse worm-like chain handler>>=
4105 double inverse_wlc_handler(double F, void *pdata)
4107 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4108 wlc_param_t param = {0};
4110 wlc_param_grouper(data, ¶m);
4111 return inverse_wlc(F, data->env->T, param.p, param.L);
4116 <<worm-like chain structure>>=
4117 typedef struct wlc_param_struct {
4118 double p; /* WLC persistence length (m) */
4119 double L; /* WLC contour length (m) */
4123 <<worm-like chain create/destroy function declarations>>=
4124 extern void *string_create_wlc_param_t(char **param_strings);
4125 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
4128 <<worm-like chain create/destroy functions>>=
4129 wlc_param_t *create_wlc_param_t(double p, double L)
4131 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
4132 assert(ret != NULL);
4138 void *string_create_wlc_param_t(char **param_strings)
4140 return create_wlc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4141 safe_strtod(param_strings[1], __FUNCTION__));
4144 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
4152 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.
4153 TODO (now we copy the string before parsing, but still do this for clarity).
4155 <<valid string write code>>=
4156 char string[] = "abc";
4159 <<valid string write code>>=
4160 char *string = "abc";
4164 <<worm-like chain tension model global declarations>>=
4165 extern int num_wlc_params;
4166 extern char *wlc_param_descriptions[];
4167 extern char wlc_param_string[];
4169 <<worm-like chain tension model globals>>=
4170 int num_wlc_params = 2;
4171 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
4172 char wlc_param_string[]="0.39e-9,28.4e-9";
4175 <<worm-like chain tension model getopt>>=
4176 {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}
4178 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
4180 \subsection{Freely-jointed chain}
4183 <<freely-jointed chain functions>>=
4184 <<freely-jointed chain function>>
4185 <<freely-jointed chain handler>>
4186 <<freely-jointed chain create/destroy functions>>
4189 <<freely-jointed chain tension model declarations>>=
4190 <<freely-jointed chain handler declaration>>
4191 <<freely-jointed chain create/destroy function declarations>>
4192 <<freely-jointed chain tension model global declarations>>
4196 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
4197 The tension of a chain of $N$ such links is given by
4199 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
4201 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}.
4202 <<freely-jointed chain function>>=
4203 double langevin(double x, void *params)
4206 return 1.0/tanh(x) - 1/x;
4208 <<one dimensional bracket declaration>>
4209 <<one dimensional solve declaration>>
4210 double inv_langevin(double x)
4212 double lb=0.0, ub=1.0, ret;
4213 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
4214 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
4218 double fjc(double x, double T, double l, int N)
4220 double L = l*(double)N;
4221 if (x == 0) return 0;
4222 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
4223 return k_B*T/l * inv_langevin(x/L);
4227 <<freely-jointed chain handler declaration>>=
4228 extern double fjc_handler(double x, void *pdata);
4230 <<freely-jointed chain handler>>=
4231 double fjc_handler(double x, void *pdata)
4233 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4234 list_t *list = data->group_tension_params;
4239 assert(list != NULL); /* empty active group?! */
4240 l = ((fjc_param_t *)list->d)->l;
4241 while (list != NULL) {
4242 /* enforce consistent l */
4243 assert( ((fjc_param_t *)list->d)->l == l);
4244 N += ((fjc_param_t *)list->d)->N;
4246 fprintf(stderr, "%s: FJC section %d links long in series\n",
4247 __FUNCTION__, ((fjc_param_t *)list->d)->N);
4252 fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
4253 __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
4255 return fjc(x, data->env->T, l, N);
4259 <<freely-jointed chain structure>>=
4260 typedef struct fjc_param_struct {
4261 double l; /* FJC link length (m) */
4262 int N; /* FJC number of links */
4266 <<freely-jointed chain create/destroy function declarations>>=
4267 extern void *string_create_fjc_param_t(char **param_strings);
4268 extern void destroy_fjc_param_t(void *p);
4271 <<freely-jointed chain create/destroy functions>>=
4272 fjc_param_t *create_fjc_param_t(double l, double N)
4274 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
4275 assert(ret != NULL);
4281 void *string_create_fjc_param_t(char **param_strings)
4283 return create_fjc_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4284 safe_strtod(param_strings[1], __FUNCTION__));
4287 void destroy_fjc_param_t(void *p)
4294 <<freely-jointed chain tension model global declarations>>=
4295 extern int num_fjc_params;
4296 extern char *fjc_param_descriptions[];
4297 extern char fjc_param_string[];
4300 <<freely-jointed chain tension model globals>>=
4301 int num_fjc_params = 2;
4302 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
4303 char fjc_param_string[]="0.5e-9,200";
4306 <<freely-jointed chain tension model getopt>>=
4307 {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}
4310 \subsection{Tension model implementation}
4312 <<tension-model.c>>=
4315 <<tension model includes>>
4316 #include "tension_model.h"
4317 <<tension model internal definitions>>
4318 <<tension model globals>>
4319 <<tension model functions>>
4322 <<tension model includes>>=
4323 #include <assert.h> /* assert() */
4324 #include <stdlib.h> /* NULL, strto*() */
4325 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4326 #include <stdio.h> /* fprintf(), stdout */
4327 #include <errno.h> /* errno, ERANGE */
4330 #include "tension_balance.h" /* oneD_*() */
4333 <<tension model internal definitions>>=
4334 <<constant tension structure>>
4336 <<worm-like chain structure>>
4337 <<freely-jointed chain structure>>
4340 <<tension model globals>>=
4341 <<tension handler array global>>
4342 <<constant tension model globals>>
4343 <<hooke tension model globals>>
4344 <<worm-like chain tension model globals>>
4345 <<freely-jointed chain tension model globals>>
4348 <<tension model functions>>=
4350 <<constant tension functions>>
4352 <<worm-like chain functions>>
4353 <<freely-jointed chain functions>>
4357 \subsection{Utilities}
4359 The tension models can get complicated, and users may want to reassure
4360 themselves that this portion of the input physics are functioning properly. The
4361 stand-alone program [[tension_model_utils]] demonstrates the tension model
4362 interface without the overhead of having to understand a full simulation.
4363 [[tension_model_utils]] takes command line model arguments like the full
4364 simulation, and prints $F(x)$ data to the screen for the selected model over a
4367 <<tension-model-utils.c>>=
4369 <<tension model utility includes>>
4370 <<tension model utility definitions>>
4371 <<tension model utility getopt functions>>
4373 <<tension model utility main>>
4376 <<tension model utility main>>=
4377 int main(int argc, char **argv)
4379 <<tension model getopt array definition>>
4380 tension_model_getopt_t *model = NULL;
4384 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
4385 tension_handler_data_t tdata;
4386 int num_param_args; /* for INIT_MODEL() */
4387 char **param_args; /* for INIT_MODEL() */
4389 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model,
4390 &Fmax, &Xmax, &flags);
4392 if (flags & VFLAG) {
4393 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4395 INIT_MODEL("utility", model, model->params, params);
4396 tdata.group_tension_params = NULL;
4397 push(&tdata.group_tension_params, params, 1);
4399 tdata.persist = NULL;
4400 if (model->handler == NULL) {
4401 printf("No tension function for model %s\n", model->name);
4407 printf("#Distance (m)\tForce (N)\n");
4408 for (i=0; i<=N; i++) {
4409 x = Xmax*i/(double)N;
4410 F = (*model->handler)(x, &tdata);
4411 if (F < 0 || F > Fmax) break;
4412 printf("%g\t%g\n", x, F);
4414 if (flags & VFLAG && i <= N) {
4415 /* explain exit condition */
4417 printf("Impossible force %g\n", F);
4419 printf("Reached large force limit %g > %g\n", F, Fmax);
4422 params = pop(&tdata.group_tension_params);
4423 if (model->destructor)
4424 (*model->destructor)(params);
4429 <<tension model utility includes>>=
4432 #include <string.h> /* strlen() */
4433 #include <assert.h> /* assert() */
4434 #include <errno.h> /* errno, ERANGE */
4438 #include "tension_model.h"
4441 <<tension model utility definitions>>=
4442 <<version definition>>
4443 #define VFLAG 1 /* verbose */
4444 <<string comparison definition and globals>>
4445 <<tension model getopt definitions>>
4446 <<initialize model definition>>
4450 <<tension model utility getopt functions>>=
4452 <<index tension model>>
4453 <<tension model utility help>>
4454 <<tension model utility get options>>
4457 <<tension model utility help>>=
4458 void help(char *prog_name,
4460 int n_tension_models, tension_model_getopt_t *tension_models,
4461 int tension_model, double Fmax, double Xmax)
4464 printf("usage: %s [options]\n", prog_name);
4465 printf("Version %s\n\n", VERSION);
4466 printf("A utility for understanding the available tension models\n\n");
4467 printf("Environmental options:\n");
4468 printf("-T T\tTemperature T (currently %g K)\n", env->T);
4469 printf("-C T\tYou can also set the temperature T in Celsius\n");
4470 printf("Model options:\n");
4471 printf("-m model\tFolded tension model (currently %s)\n",
4472 tension_models[tension_model].name);
4473 printf("-a args\tFolded tension model argument string (currently %s)\n",
4474 tension_models[tension_model].params);
4475 printf("F(x) is calculated for a range of x and printed\n");
4476 printf("For example:\n");
4477 printf(" #Distance (m)\tForce (N)\n");
4478 printf(" 123.456\t7.89\n");
4480 printf("-F\tSet the maximum F value for the standard mode F(x) output (currently %g)\n", Fmax);
4481 printf("-X\tSet the maximum x value for the standard mode F(x) output (currently %g)\n", Xmax);
4482 printf("-V\tChange output to verbose mode\n");
4483 printf("-h\tPrint this help and exit\n");
4485 printf("Tension models:\n");
4486 for (i=0; i<n_tension_models; i++) {
4487 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
4488 for (j=0; j < tension_models[i].num_params ; j++)
4489 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
4490 printf("\t\tdefaults: %s\n", tension_models[i].params);
4495 <<tension model utility get options>>=
4496 void get_options(int argc, char **argv, environment_t *env,
4497 int n_tension_models, tension_model_getopt_t *tension_models,
4498 tension_model_getopt_t **model, double *Fmax, double *Xmax,
4499 unsigned int *flags)
4501 char *prog_name = NULL;
4502 char c, options[] = "T:C:m:a:F:X:Vh";
4503 int tension_model=0, num_strings;
4504 extern char *optarg;
4505 extern int optind, optopt, opterr;
4509 /* setup defaults */
4511 prog_name = argv[0];
4512 env->T = 300.0; /* K */
4513 *Fmax = 1e5; /* N */
4514 *Xmax = 1e-6; /* m */
4516 *model = tension_models;
4518 while ((c=getopt(argc, argv, options)) != -1) {
4520 case 'T': env->T = safe_strtod(optarg, "-T"); break;
4521 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
4523 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
4524 *model = tension_models+tension_model;
4527 tension_models[tension_model].params = optarg;
4529 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
4530 case 'X': *Xmax = safe_strtod(optarg, "-X"); break;
4531 case 'V': *flags |= VFLAG; break;
4533 fprintf(stderr, "unrecognized option '%c'\n", optopt);
4534 /* fall through to default case */
4536 help(prog_name, env, n_tension_models, tension_models, tension_model, *Fmax, *Xmax);
4545 \section{Unfolding rate models}
4546 \label{sec.k_models}
4548 We have two main choices for determining $k$: Bell's model, which assumes the
4549 folded domains exist in a single `folded' state and make instantaneous,
4550 stochastic transitions over a free energy barrier; and Kramers' model, which
4551 treats the unfolding as a thermalized diffusion process.
4552 We deal with these options like we dealt with the different tension models: by bundling all model-specific
4553 parameters into structures, with model specific functions for determining $k$.
4555 <<k func definition>>=
4556 typedef double k_func_t(double F, environment_t *env, void *params);
4559 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
4560 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]].
4562 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
4563 because \LaTeX\ doesn't like underscores outside math mode.
4564 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
4565 You could use spaces instead (`\verb+ +'),
4566 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
4567 which I don't like as much.
4573 <<k func definition>>
4574 <<null k declarations>>
4575 <<const k declarations>>
4576 <<bell k declarations>>
4577 <<kbell k declarations>>
4578 <<kramers k declarations>>
4579 <<kramers integ k declarations>>
4580 #endif /* K_MODEL_H */
4583 <<k model module makefile lines>>=
4584 NW_SAWSIM_MODS += k_model
4585 CHECK_BINS += check_k_model
4586 check_k_model_MODS = parse string_eval
4588 [[check_k_model]] also depends on the parse module.
4590 <<k model utils makefile lines>>=
4591 K_MODEL_MODS = k_model parse string_eval
4592 K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
4593 $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h)
4594 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
4595 notangle -Rk-model-utils.c $< > $@
4596 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
4597 gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
4598 clean_k_model_utils :
4599 rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
4605 For domains that are never folded.
4607 <<null k declarations>>=
4609 <<null k functions>>=
4614 <<null k model getopt>>=
4615 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
4618 \subsection{Constant rate model}
4621 This is a pretty boring model, but it is usefull for testing the engine.
4625 where $k_0$ is the constant unfolding rate.
4627 <<const k functions>>=
4628 <<const k function>>
4629 <<const k structure create/destroy functions>>
4632 <<const k declarations>>=
4633 <<const k function declaration>>
4634 <<const k structure create/destroy function declarations>>
4635 <<const k global declarations>>
4639 <<const k structure>>=
4640 typedef struct const_k_param_struct {
4641 double knot; /* transition rate k_0 (frac population per s) */
4645 <<const k function declaration>>=
4646 double const_k(double F, environment_t *env, void *const_k_params);
4648 <<const k function>>=
4649 double const_k(double F, environment_t *env, void *const_k_params)
4650 { /* Returns the rate constant k in frac pop per s. */
4651 const_k_param_t *p = (const_k_param_t *) const_k_params;
4653 assert(p->knot > 0);
4658 <<const k structure create/destroy function declarations>>=
4659 void *string_create_const_k_param_t(char **param_strings);
4660 void destroy_const_k_param_t(void *p);
4663 <<const k structure create/destroy functions>>=
4664 const_k_param_t *create_const_k_param_t(double knot)
4666 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
4667 assert(ret != NULL);
4672 void *string_create_const_k_param_t(char **param_strings)
4674 return create_const_k_param_t(safe_strtod(param_strings[0], __FUNCTION__));
4677 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
4684 <<const k global declarations>>=
4685 extern int num_const_k_params;
4686 extern char *const_k_param_descriptions[];
4687 extern char const_k_param_string[];
4690 <<const k globals>>=
4691 int num_const_k_params = 1;
4692 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
4693 char const_k_param_string[]="1";
4696 <<const k model getopt>>=
4697 {"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}
4700 \subsection{Bell's model}
4703 Quantitatively, Bell's model gives $k$ as
4705 k = k_0 \cdot e^{F dx / k_B T} \;,
4707 where $k_0$ is the unforced unfolding rate,
4708 $F$ is the applied unfolding force,
4709 $dx$ is the distance to the transition state, and
4710 $k_B T$ is the thermal energy\citep{schlierf06}.
4712 <<bell k functions>>=
4714 <<bell k structure create/destroy functions>>
4717 <<bell k declarations>>=
4718 <<bell k function declaration>>
4719 <<bell k structure create/destroy function declarations>>
4720 <<bell k global declarations>>
4724 <<bell k structure>>=
4725 typedef struct bell_param_struct {
4726 double knot; /* transition rate k_0 (frac population per s) */
4727 double dx; /* `distance to transition state' (nm) */
4731 <<bell k function declaration>>=
4732 double bell_k(double F, environment_t *env, void *bell_params);
4734 <<bell k function>>=
4735 double bell_k(double F, environment_t *env, void *bell_params)
4736 { /* Returns the rate constant k in frac pop per s.
4737 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4738 * uses global k_B in J/K */
4739 bell_param_t *p = (bell_param_t *) bell_params;
4740 assert(F >= 0); assert(env->T > 0);
4742 assert(p->knot > 0); //assert(p->dx >= 0); /* p < 0 for refolding reactions, etc. */
4744 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
4745 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
4746 p->knot * exp(F*p->dx / (k_B*env->T) ));
4747 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
4749 return p->knot * exp(F*p->dx / (k_B*env->T) );
4753 <<bell k structure create/destroy function declarations>>=
4754 void *string_create_bell_param_t(char **param_strings);
4755 void destroy_bell_param_t(void *p);
4758 <<bell k structure create/destroy functions>>=
4759 bell_param_t *create_bell_param_t(double knot, double dx)
4761 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
4762 assert(ret != NULL);
4768 void *string_create_bell_param_t(char **param_strings)
4770 return create_bell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4771 safe_strtod(param_strings[1], __FUNCTION__));
4774 void destroy_bell_param_t(void *p /* bell_param_t * */)
4781 <<bell k global declarations>>=
4782 extern int num_bell_params;
4783 extern char *bell_param_descriptions[];
4784 extern char bell_param_string[];
4788 int num_bell_params = 2;
4789 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4790 char bell_param_string[]="3.3e-4,0.25e-9";
4793 <<bell k model getopt>>=
4794 {"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}
4796 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4797 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4800 \subsection{Linker stiffness corrected Bell model}
4803 Walton et. al showed that the Bell model constant force approximation
4804 doesn't hold for biotin-streptavadin binding, and I extended their
4805 results to I27 for the 2009 Biophysical Society Annual
4806 Meeting\cite{walton08,king09}. More details to follow when I get done
4807 with the conference\ldots
4809 We adjust Bell's model to
4811 k = k_0 \cdot e^{\frac{F dx - \frac{1}{2} \kappa dx^2}{k_B T}} \;,
4813 where $k_0$ is the unforced unfolding rate,
4814 $F$ is the applied unfolding force,
4815 $\kappa$ is the effective spring constant,
4816 $dx$ is the distance to the transition state, and
4817 $k_B T$ is the thermal energy\citep{schlierf06}.
4819 <<kbell k functions>>=
4820 <<kbell k function>>
4821 <<kbell k structure create/destroy functions>>
4824 <<kbell k declarations>>=
4825 <<kbell k function declaration>>
4826 <<kbell k structure create/destroy function declarations>>
4827 <<kbell k global declarations>>
4831 <<kbell k structure>>=
4832 typedef struct kbell_param_struct {
4833 double knot; /* transition rate k_0 (frac population per s) */
4834 double dx; /* `distance to transition state' (nm) */
4838 <<kbell k function declaration>>=
4839 double kbell_k(double F, environment_t *env, void *kbell_params);
4841 <<kbell k function>>=
4842 double kbell_k(double F, environment_t *env, void *kbell_params)
4843 { /* Returns the rate constant k in frac pop per s.
4844 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4845 * uses global k_B in J/K */
4846 kbell_param_t *p = (kbell_param_t *) kbell_params;
4847 assert(F >= 0); assert(env->T > 0); assert(env->stiffness >= 0);
4849 assert(p->knot > 0); assert(p->dx >= 0);
4850 return p->knot * exp((F*p->dx - 0.5*env->stiffness*(p->dx)*(p->dx))/(k_B*env->T) );
4854 <<kbell k structure create/destroy function declarations>>=
4855 void *string_create_kbell_param_t(char **param_strings);
4856 void destroy_kbell_param_t(void *p);
4859 <<kbell k structure create/destroy functions>>=
4860 kbell_param_t *create_kbell_param_t(double knot, double dx)
4862 kbell_param_t *ret = (kbell_param_t *) malloc(sizeof(kbell_param_t));
4863 assert(ret != NULL);
4869 void *string_create_kbell_param_t(char **param_strings)
4871 return create_kbell_param_t(safe_strtod(param_strings[0], __FUNCTION__),
4872 safe_strtod(param_strings[1], __FUNCTION__));
4875 void destroy_kbell_param_t(void *p /* kbell_param_t * */)
4882 <<kbell k global declarations>>=
4883 extern int num_kbell_params;
4884 extern char *kbell_param_descriptions[];
4885 extern char kbell_param_string[];
4888 <<kbell k globals>>=
4889 int num_kbell_params = 2;
4890 char *kbell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4891 char kbell_param_string[]="3.3e-4,0.25e-9";
4894 <<kbell k model getopt>>=
4895 {"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}
4897 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4898 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4901 \subsection{Kramer's model}
4904 Kramer's model gives $k$ as
4906 % \frac{1}{k} = \frac{1}{D}
4907 % \int_{x_\text{min}}^{x_\text{max}}
4908 % e^{\frac{-U_F(x)}{k_B T}}
4909 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4912 %where $D$ is the diffusion constant,
4913 %$U_F(x)$ is the free energy along the unfolding coordinate, and
4914 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4915 % before the minimum and after the maximum, respectively \citep{schlierf06}.
4917 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
4919 where $D$ is the diffusion constant,
4921 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
4923 are length scales where
4924 $x_c(F)$ is the location of the bound state and
4925 $x_{ts}(F)$ is the location of the transition state,
4926 $E(F,x)$ is the free energy, and
4927 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
4928 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
4929 \citet{evans97} Eqn.~3).
4931 <<kramers k functions>>=
4932 <<kramers k function>>
4933 <<kramers k structure create/destroy functions>>
4936 <<kramers k declarations>>=
4937 <<kramers k function declaration>>
4938 <<kramers k structure create/destroy function declarations>>
4939 <<kramers k global declarations>>
4943 <<kramers k structure>>=
4944 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
4945 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
4947 typedef struct kramers_param_struct {
4948 double D; /* diffusion rate (um/s) */
4955 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
4956 //kramers_x_func_t *xts; /* function returning transition x (nm) */
4957 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
4958 //kramers_E_func_t *E; /* function returning E (J) */
4959 //double *E_params; /* parametrize the energy landscape */
4960 //int n_E_params; /* length of E_params array */
4964 <<kramers k function declaration>>=
4965 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
4966 extern double kramers_k(double F, environment_t *env, void *kramers_params);
4968 <<kramers k function>>=
4969 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
4971 static char input[160]={0};
4972 static char output[80]={0};
4973 /* setup the environment */
4974 fprintf(in, "F = %g; T = %g\n", F, T);
4975 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4976 string_eval(in, out, input, 80, output);
4977 return safe_strtod(output, __FUNCTION__);
4980 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
4982 static char input[160]={0};
4983 static char output[80]={0};
4984 /* setup the environment */
4985 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
4986 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4987 string_eval(in, out, input, 80, output);
4988 return safe_strtod(output, __FUNCTION__);
4991 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
4993 kramers_param_t *p = (kramers_param_t *) kramers_params;
4994 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
4997 double kramers_k(double F, environment_t *env, void *kramers_params)
4998 { /* Returns the rate constant k in frac pop per s.
4999 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
5000 * uses global k_B in J/K */
5001 kramers_param_t *p = (kramers_param_t *) kramers_params;
5002 double kbT, xc, xts, lc, lts, Eb;
5003 assert(F >= 0); assert(env->T > 0);
5006 assert(p->xc != NULL); assert(p->xts != NULL);
5007 assert(p->ddEdxx != NULL); assert(p->E != NULL);
5008 assert(p->in != NULL); assert(p->out != NULL);
5010 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
5011 if (xc == -1.0) return -1.0; /* error (Too much force) */
5012 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
5013 if (xc == -1.0) return -1.0; /* error (Too much force) */
5014 lc = sqrt(2.0*M_PI*kbT /
5015 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
5016 lts = sqrt(-2.0*M_PI*kbT/
5017 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
5018 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
5019 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
5020 //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));
5021 return p->D/(lc*lts) * exp(-Eb/kbT);
5025 <<kramers k structure create/destroy function declarations>>=
5026 void *string_create_kramers_param_t(char **param_strings);
5027 void destroy_kramers_param_t(void *p);
5030 <<kramers k structure create/destroy functions>>=
5031 kramers_param_t *create_kramers_param_t(double D,
5032 char *xc_fn, char *xts_fn,
5033 char *E_fn, char *ddEdxx_fn,
5034 char *global_define_string)
5035 // kramers_x_func_t *xc, kramers_x_func_t *xts,
5036 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
5037 // double *E_params, int n_E_params)
5039 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
5040 assert(ret != NULL);
5041 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5042 assert(ret->xc != NULL);
5043 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
5044 assert(ret->xts != NULL);
5045 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
5046 assert(ret->E != NULL);
5047 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
5048 assert(ret->ddEdxx != NULL);
5050 strcpy(ret->xc, xc_fn);
5051 strcpy(ret->xts, xts_fn);
5052 strcpy(ret->E, E_fn);
5053 strcpy(ret->ddEdxx, ddEdxx_fn);
5054 //ret->E_params = E_params;
5055 //ret->n_E_params = n_E_params;
5056 string_eval_setup(&ret->in, &ret->out);
5057 fprintf(ret->in, "kB = %g\n", k_B);
5058 fprintf(ret->in, "%s\n", global_define_string);
5062 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
5063 void *string_create_kramers_param_t(char **param_strings)
5065 return create_kramers_param_t(safe_strtod(param_strings[0], __FUNCTION__),
5073 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
5075 kramers_param_t *pk = (kramers_param_t *)p;
5077 string_eval_teardown(&pk->in, &pk->out);
5079 // free(pk->E_params);
5085 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
5086 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.
5087 We follow \citet{shillcock98} and use
5089 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
5090 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
5093 where TODO, variable meanings.%\citep{shillcock98}.
5094 The first and second derivatives of $E(F,E_b,x,x_s)$ are
5096 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\\
5097 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
5099 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
5100 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
5101 The roots are therefor at
5103 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
5104 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
5105 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
5108 <<kramers k global declarations>>=
5109 extern int num_kramers_params;
5110 extern char *kramers_param_descriptions[];
5111 extern char kramers_param_string[];
5114 <<kramers k globals>>=
5115 int num_kramers_params = 6;
5116 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"};
5117 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)}";
5119 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
5120 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
5121 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
5122 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.
5123 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
5124 It works so long as [[val_1]] is not [[false]].
5126 <<kramers k model getopt>>=
5127 {"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}
5130 \subsection{Kramer's model (natural cubic splines)}
5131 \label{sec.kramers_integ}
5133 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
5134 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
5135 \citet{schlierf06} Eqn.~2)
5137 \frac{1}{k} = \frac{1}{D}
5138 \int_{x_\text{min}}^{x_\text{max}}
5139 e^{\frac{U_F(x)}{k_B T}}
5140 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5143 where $D$ is the diffusion constant,
5144 $U_F(x)$ is the free energy along the unfolding coordinate, and
5145 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
5146 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
5148 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
5149 interpolating between them with cubic splines.
5150 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
5152 <<kramers integ k functions>>=
5153 <<spline functions>>
5154 <<kramers integ A k functions>>
5155 <<kramers integ B k functions>>
5156 <<kramers integ k function>>
5157 <<kramers integ k structure create/destroy functions>>
5160 <<kramers integ k declarations>>=
5161 <<kramers integ k function declaration>>
5162 <<kramers integ k structure create/destroy function declarations>>
5163 <<kramers integ k global declarations>>
5167 <<kramers integ k structure>>=
5168 typedef double func_t(double x, void *params);
5169 typedef struct kramers_integ_param_struct {
5170 double D; /* diffusion rate (um/s) */
5171 double x_min; /* integration bounds */
5173 func_t *E; /* function returning E (J) */
5174 void *E_params; /* parametrize the energy landscape */
5175 destroy_data_func_t *destroy_E_params;
5177 double F; /* for passing into GSL evaluations */
5179 } kramers_integ_param_t;
5182 <<spline param structure>>=
5183 typedef struct spline_param_struct {
5184 int n_knots; /* natural cubic spline knots for U(x) */
5187 gsl_spline *spline; /* GSL spline parameters */
5188 gsl_interp_accel *acc; /* GSL spline acceleration data */
5192 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
5194 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5196 <<kramers integ A k functions>>=
5197 double kramers_integ_k_integrandA(double x, void *params)
5199 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5200 double U = (*p->E)(x, p->E_params);
5201 double Fx = p->F * x;
5202 double kbT = k_B * p->env->T;
5203 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
5204 return exp(-(U-Fx)/kbT);
5206 double kramers_integ_k_integralA(double x, void *params)
5208 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5210 double result, abserr;
5211 size_t neval; /* for qng */
5212 static gsl_integration_workspace *w = NULL;
5213 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5214 f.function = &kramers_integ_k_integrandA;
5216 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5217 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5218 //fprintf(stderr, "integralA = %g\n", result);
5224 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
5225 \int_{x_\text{min}}^{x_\text{max}}
5226 e^{\frac{U_F(x)}{k_B T}}
5227 \text{Integral}_A(x)
5230 <<kramers integ B k functions>>=
5231 double kramers_integ_k_integrandB(double x, void *params)
5233 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5234 double U = (*p->E)(x, p->E_params);
5235 double Fx = p->F * x;
5236 double kbT = k_B * p->env->T;
5237 double intA = kramers_integ_k_integralA(x, params);
5238 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
5239 return exp((U-Fx)/kbT)*intA;
5241 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
5243 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5245 double result, abserr;
5246 size_t neval; /* for qng */
5247 static gsl_integration_workspace *w = NULL;
5248 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5249 f.function = &kramers_integ_k_integrandB;
5253 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5254 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5255 //fprintf(stderr, "integralB = %g\n", result);
5260 With the integrals taken care of, evaluating $k$ is simply
5262 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
5264 <<kramers integ k function declaration>>=
5265 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
5267 <<kramers integ k function>>=
5268 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
5269 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
5270 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
5271 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
5275 <<kramers integ k structure create/destroy function declarations>>=
5276 void *string_create_kramers_integ_param_t(char **param_strings);
5277 void destroy_kramers_integ_param_t(void *p);
5280 <<kramers integ k structure create/destroy functions>>=
5281 kramers_integ_param_t *create_kramers_integ_param_t(double D,
5282 double xmin, double xmax, func_t *E,
5284 destroy_data_func_t *destroy_E_params)
5286 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
5287 assert(ret != NULL);
5292 ret->E_params = E_params;
5293 ret->destroy_E_params = destroy_E_params;
5297 void *string_create_kramers_integ_param_t(char **param_strings)
5299 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
5300 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
5301 void *E_params = string_create_spline_param_t(param_strings+1);
5302 return create_kramers_integ_param_t(
5303 safe_strtod(param_strings[0], __FUNCTION__),
5304 safe_strtod(param_strings[2], __FUNCTION__),
5305 safe_strtod(param_strings[3], __FUNCTION__),
5306 &spline_eval, E_params, destroy_spline_param_t);
5309 void destroy_kramers_integ_param_t(void *params)
5311 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5314 (*p->destroy_E_params)(p->E_params);
5320 Finally we have the GSL spline wrappers:
5321 <<spline functions>>=
5323 <<create/destroy spline>>
5326 <<spline function>>=
5327 double spline_eval(double x, void *spline_params)
5329 spline_param_t *p = (spline_param_t *)(spline_params);
5330 return gsl_spline_eval(p->spline, x, p->acc);
5334 <<create/destroy spline>>=
5335 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
5337 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
5338 assert(ret != NULL);
5339 ret->n_knots = n_knots;
5342 ret->acc = gsl_interp_accel_alloc();
5343 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
5344 gsl_spline_init(ret->spline, x, y, n_knots);
5348 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
5349 void *string_create_spline_param_t(char **param_strings)
5353 double *x=NULL, *y=NULL;
5354 /* break into ordered pairs */
5355 parse_list_string(param_strings[0], SEP, '(', ')',
5356 &num_knots, &knot_coords);
5357 x = (double *)malloc(sizeof(double)*num_knots);
5359 y = (double *)malloc(sizeof(double)*num_knots);
5361 for (i=0; i<num_knots; i++) {
5362 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
5363 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
5366 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
5368 free_parsed_list(num_knots, knot_coords);
5369 return create_spline_param_t(num_knots, x, y);
5372 void destroy_spline_param_t(void *params)
5374 spline_param_t *p = (spline_param_t *)params;
5377 gsl_spline_free(p->spline);
5379 gsl_interp_accel_free(p->acc);
5389 <<kramers integ k global declarations>>=
5390 extern int num_kramers_integ_params;
5391 extern char *kramers_integ_param_descriptions[];
5392 extern char kramers_integ_param_string[];
5395 <<kramers integ k globals>>=
5396 int num_kramers_integ_params = 4;
5397 char *kramers_integ_param_descriptions[] = {
5398 "diffusion D, m^2 / s",
5399 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
5400 "starting integration bound x_min, m",
5401 "ending integration bound x_max, m"};
5402 //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";
5403 //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";
5404 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";
5405 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
5406 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
5407 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
5410 <<kramers integ k model getopt>>=
5411 {"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}
5413 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
5414 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5416 \subsection{Unfolding model implementation}
5420 <<k model includes>>
5421 #include "k_model.h"
5422 <<k model internal definitions>>
5424 <<k model functions>>
5427 <<k model includes>>=
5428 #include <assert.h> /* assert() */
5429 #include <stdlib.h> /* NULL, malloc(), strto*() */
5430 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
5431 #include <stdio.h> /* fprintf(), stdout */
5432 #include <string.h> /* strlen(), strcpy() */
5433 #include <errno.h> /* errno, ERANGE */
5434 #include <gsl/gsl_integration.h>
5435 #include <gsl/gsl_spline.h>
5440 <<k model internal definitions>>=
5441 <<const k structure>>
5442 <<bell k structure>>
5443 <<kbell k structure>>
5444 <<kramers k structure>>
5445 <<kramers integ k structure>>
5446 <<spline param structure>>
5449 <<k model globals>>=
5454 <<kramers k globals>>
5455 <<kramers integ k globals>>
5458 <<k model functions>>=
5460 <<null k functions>>
5461 <<const k functions>>
5462 <<bell k functions>>
5463 <<kbell k functions>>
5464 <<kramers k functions>>
5465 <<kramers integ k functions>>
5468 \subsection{Unfolding model unit tests}
5470 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
5471 <<check-k-model.c>>=
5472 <<k model check includes>>
5473 <<check relative error>>
5475 <<k model test suite>>
5476 <<main check program>>
5479 <<k model check includes>>=
5480 #include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
5481 #include <stdio.h> /* sprintf() */
5482 #include <assert.h> /* assert() */
5483 #include <math.h> /* exp() */
5484 #include <errno.h> /* errno, ERANGE */
5487 #include "k_model.h"
5490 <<k model test suite>>=
5494 <<k model suite definition>>
5497 <<k model suite definition>>=
5498 Suite *test_suite (void)
5500 Suite *s = suite_create ("k model");
5501 <<const k test case defs>>
5502 <<bell k test case defs>>
5504 <<const k test case adds>>
5505 <<bell k test case adds>>
5510 \subsubsection{Constant}
5512 <<const k test case defs>>=
5513 TCase *tc_const_k = tcase_create("Constant k");
5516 <<const k test case adds>>=
5517 tcase_add_test(tc_const_k, test_const_k_create_destroy);
5518 tcase_add_test(tc_const_k, test_const_k_over_range);
5519 suite_add_tcase(s, tc_const_k);
5524 START_TEST(test_const_k_create_destroy)
5526 char *knot[] = {"1","2","3","4","5","6"};
5527 char *params[] = {knot[0]};
5530 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5531 params[0] = knot[i];
5532 p = string_create_const_k_param_t(params);
5533 destroy_const_k_param_t(p);
5538 START_TEST(test_const_k_over_range)
5540 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
5541 char *knot[] = {"1","2","3","4","5","6"};
5542 char *params[] = {knot[0]};
5545 char param_string[80];
5548 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5549 params[0] = knot[i];
5550 p = string_create_const_k_param_t(params);
5551 for ( F=Fm; F<FM; F+=dF ) {
5552 fail_unless(const_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
5553 "const_k(%g, %g, &{%s}) = %g != %s",
5554 F, T, knot[i], const_k(F, &env, p), knot[i]);
5556 destroy_const_k_param_t(p);
5562 \subsubsection{Bell}
5564 <<bell k test case defs>>=
5565 TCase *tc_bell = tcase_create("Bell's k");
5568 <<bell k test case adds>>=
5569 tcase_add_test(tc_bell, test_bell_k_create_destroy);
5570 tcase_add_test(tc_bell, test_bell_k_at_zero);
5571 tcase_add_test(tc_bell, test_bell_k_at_one);
5572 suite_add_tcase(s, tc_bell);
5576 START_TEST(test_bell_k_create_destroy)
5578 char *knot[] = {"1","2","3","4","5","6"};
5580 char *params[] = {knot[0], dx};
5583 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5584 params[0] = knot[i];
5585 p = string_create_bell_param_t(params);
5586 destroy_bell_param_t(p);
5591 START_TEST(test_bell_k_at_zero)
5593 double F=0.0, T=300.0;
5594 char *knot[] = {"1","2","3","4","5","6"};
5596 char *params[] = {knot[0], dx};
5599 char param_string[80];
5602 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5603 params[0] = knot[i];
5604 p = string_create_bell_param_t(params);
5605 fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
5606 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
5607 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
5608 destroy_bell_param_t(p);
5613 START_TEST(test_bell_k_at_one)
5616 char *knot[] = {"1","2","3","4","5","6"};
5618 char *params[] = {knot[0], dx};
5619 double F= k_B*T/safe_strtod(dx, __FUNCTION__);
5622 char param_string[80];
5625 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5626 params[0] = knot[i];
5627 p = string_create_bell_param_t(params);
5628 CHECK_ERR(1e-6, safe_strtod(knot[i], __FUNCTION__)*exp(1.0), bell_k(F, &env, p));
5629 //fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__)*exp(1.0),
5630 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
5631 // F, T, knot[i], dx, bell_k(F, &env, p), safe_strtod(knot[i], __FUNCTION__)*exp(1.0));
5632 destroy_bell_param_t(p);
5638 <<kramers k tests>>=
5641 <<kramers k test case def>>=
5644 <<kramers k test case add>>=
5647 <<k model function tests>>=
5648 <<check relative error>>
5650 START_TEST(test_something)
5652 double k=5, x=3, last_x=2.0, F;
5653 one_dim_fn_t *handlers[] = {&hooke};
5654 void *data[] = {&k};
5655 double xi[] = {0.0};
5657 int new_active_groups = 1;
5658 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
5659 fail_unless(F = k*x, NULL);
5664 \subsection{Utilities}
5666 The unfolding models can get complicated, and are sometimes hard to explain to others.
5667 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
5668 the overhead of having to understand a full simulation.
5669 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
5670 or special mode, where the operation depends on the specific model selected.
5672 <<k-model-utils.c>>=
5674 <<k model utility includes>>
5675 <<k model utility definitions>>
5676 <<k model utility getopt functions>>
5677 <<k model utility multi model E>>
5678 <<k model utility main>>
5681 <<k model utility main>>=
5682 int main(int argc, char **argv)
5684 <<k model getopt array definition>>
5685 k_model_getopt_t *model = NULL;
5690 int num_param_args; /* for INIT_MODEL() */
5691 char **param_args; /* for INIT_MODEL() */
5693 double special_xmin;
5694 double special_xmax;
5695 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
5696 &Fmax, &special_xmin, &special_xmax, &flags);
5697 if (flags & VFLAG) {
5698 printf("#initializing model %s with parameters %s\n", model->name, model->params);
5700 INIT_MODEL("utility", model, model->params, params);
5703 if (model->k == NULL) {
5704 printf("No k function for model %s\n", model->name);
5708 printf("Fmax = %g <= 0\n", Fmax);
5714 printf("#F (N)\tk (%% pop. per s)\n");
5715 for (i=0; i<=N; i++) {
5716 F = Fmax*i/(double)N;
5717 k = (*model->k)(F, &env, params);
5719 printf("%g\t%g\n", F, k);
5724 if (model->k == NULL || model->k == &bell_k) {
5725 printf("No special mode for model %s\n", model->name);
5728 if (special_xmax <= special_xmin) {
5729 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
5733 double Xrng=(special_xmax-special_xmin), x, E;
5735 printf("#x (m)\tE (J)\n");
5736 for (i=0; i<=N; i++) {
5737 x = special_xmin + Xrng*i/(double)N;
5738 E = multi_model_E(model->k, params, &env, x);
5739 printf("%g\t%g\n", x, E);
5746 if (model->destructor)
5747 (*model->destructor)(params);
5752 <<k model utility multi model E>>=
5753 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
5755 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
5757 if (k == kramers_integ_k)
5758 E = (*p->E)(x, p->E_params);
5760 E = kramers_E(0, env, model_params, x);
5766 <<k model utility includes>>=
5769 #include <string.h> /* strlen() */
5770 #include <assert.h> /* assert() */
5771 #include <errno.h> /* errno, ERANGE */
5774 #include "string_eval.h"
5775 #include "k_model.h"
5778 <<k model utility definitions>>=
5779 <<version definition>>
5780 #define VFLAG 1 /* verbose */
5781 enum mode_t {M_K_OF_F, M_SPECIAL};
5782 <<string comparison definition and globals>>
5783 <<k model getopt definitions>>
5784 <<initialize model definition>>
5785 <<kramers k structure>>
5786 <<kramers integ k structure>>
5790 <<k model utility getopt functions>>=
5793 <<k model utility help>>
5794 <<k model utility get options>>
5797 <<k model utility help>>=
5798 void help(char *prog_name,
5800 int n_k_models, k_model_getopt_t *k_models,
5801 int k_model, double Fmax, double special_xmin, double special_xmax)
5804 printf("usage: %s [options]\n", prog_name);
5805 printf("Version %s\n\n", VERSION);
5806 printf("A utility for understanding the available unfolding models\n\n");
5807 printf("Environmental options:\n");
5808 printf("-T T\tTemperature T (currently %g K)\n", env->T);
5809 printf("-C T\tYou can also set the temperature T in Celsius\n");
5810 printf("Model options:\n");
5811 printf("-k model\tTransition rate model (currently %s)\n",
5812 k_models[k_model].name);
5813 printf("-K args\tTransition rate model argument string (currently %s)\n",
5814 k_models[k_model].params);
5815 printf("There are two output modes. In standard mode, k(F) is printed\n");
5816 printf("For example:\n");
5817 printf(" #Force (N)\tk (% pop. per s)\n");
5818 printf(" 123.456\t7.89\n");
5820 printf("In special mode, the output depends on the model.\n");
5821 printf("For models defining an energy function E(x), that function is printed\n");
5822 printf(" #Distance (m)\tE (J)\n");
5823 printf(" 0.0012\t0.0034\n");
5825 printf("-m\tChange output to standard mode\n");
5826 printf("-M\tChange output to special mode\n");
5827 printf("-F\tSet the maximum F value for the standard mode k(F) output (currently %g)\n", Fmax);
5828 printf("-x\tSet the minimum x value for the special mode E(x) output (currently %g)\n", special_xmin);
5829 printf("-X\tSet the maximum x value for the special mode E(x) output (currently %g)\n", special_xmax);
5830 printf("-V\tChange output to verbose mode\n");
5831 printf("-h\tPrint this help and exit\n");
5833 printf("Unfolding rate models:\n");
5834 for (i=0; i<n_k_models; i++) {
5835 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
5836 for (j=0; j < k_models[i].num_params ; j++)
5837 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
5838 printf("\t\tdefaults: %s\n", k_models[i].params);
5843 <<k model utility get options>>=
5844 void get_options(int argc, char **argv, environment_t *env,
5845 int n_k_models, k_model_getopt_t *k_models,
5846 enum mode_t *mode, k_model_getopt_t **model,
5847 double *Fmax, double *special_xmin, double *special_xmax,
5848 unsigned int *flags)
5850 char *prog_name = NULL;
5851 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
5853 extern char *optarg;
5854 extern int optind, optopt, opterr;
5858 /* setup defaults */
5860 prog_name = argv[0];
5861 env->T = 300.0; /* K */
5865 *Fmax = 1e-9; /* N */
5866 *special_xmax = 1e-8;
5867 *special_xmin = 0.0;
5869 while ((c=getopt(argc, argv, options)) != -1) {
5871 case 'T': env->T = safe_strtod(optarg, "-T"); break;
5872 case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break;
5874 k_model = index_k_model(n_k_models, k_models, optarg);
5875 *model = k_models+k_model;
5878 k_models[k_model].params = optarg;
5880 case 'm': *mode = M_K_OF_F; break;
5881 case 'M': *mode = M_SPECIAL; break;
5882 case 'F': *Fmax = safe_strtod(optarg, "-F"); break;
5883 case 'x': *special_xmin = safe_strtod(optarg, "-x"); break;
5884 case 'X': *special_xmax = safe_strtod(optarg, "-X"); break;
5885 case 'V': *flags |= VFLAG; break;
5887 fprintf(stderr, "unrecognized option '%c'\n", optopt);
5888 /* fall through to default case */
5890 help(prog_name, env, n_k_models, k_models, k_model, *Fmax, *special_xmin, *special_xmax);
5899 \section{\LaTeX\ documentation}
5901 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
5902 The comment blocks are extracted (with nicely formatted code blocks), using
5903 <<latex makefile lines>>=
5904 $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5905 noweave -latex -index -delay $< > $@
5906 $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR)
5910 <<latex makefile lines>>=
5911 $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
5913 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5914 $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
5915 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5916 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5917 mv $(BUILD_DIR)/sawsim.pdf $@
5919 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
5920 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
5921 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
5923 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
5924 <<latex makefile lines>>=
5926 rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
5927 $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
5928 $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.toc \
5929 $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib
5930 clean_latex : semi-clean_latex
5931 rm -f $(DOC_DIR)/sawsim.pdf
5936 [[make]] is a common utility on *nix systems for managing dependencies while building software.
5937 In this case, we have several \emph{targets} that we'd like to build:
5938 the main simulation program \prog;
5939 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
5940 the documentation [[sawsim.pdf]];
5941 and the [[Makefile]] itself.
5942 Besides the generated files, there is the utility target
5943 [[clean]] for removing all generated files except [[Makefile]].
5945 % [[dist]] for generating a distributable tar file.
5947 Extract the makefile with
5948 `[[notangle -Rmakefile src/sawsim.nw | sed 's/ /\t/' > Makefile]]'.
5949 The sed is needed because notangle replaces tabs with eight spaces,
5950 but [[make]] needs tabs.
5952 # Decide what directories to put things in
5957 # Note: Cannot use, for example, `./build', because make eats the `./'
5958 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5960 # Modules (X.c, X.h) defined in the noweb file
5963 # Modules defined outside the noweb file
5964 FREE_SAWSIM_MODS = interp tavl
5966 <<list module makefile lines>>
5967 <<tension balance module makefile lines>>
5968 <<tension model module makefile lines>>
5969 <<k model module makefile lines>>
5970 <<parse module makefile lines>>
5971 <<string eval module makefile lines>>
5972 <<accel k module makefile lines>>
5974 SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS)
5976 # Everything needed for sawsim under one roof. sawsim.c must come first
5977 SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \
5978 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h)
5979 # Libraries needed to compile sawsim
5980 LIBS = gsl gslcblas m
5981 CHECK_LIBS = gsl gslcblas m check
5982 # The non-check binaries generated
5983 BINS = sawsim tension_model_utils k_model_utils sawsim_profile
5986 # Define the major targets
5987 all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
5989 view : $(DOC_DIR)/sawsim.pdf
5991 profile : $(BIN_DIR)/sawsim_profile | $(BIN_DIR)
5992 $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \
5993 -mnull -Mwlc -A0.39e-9,28e-9 -F8
5994 gprof $(BIN_DIR)/sawsim_profile gmon.out > $@
5995 check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim
5996 $(SHELL) -e -c 'for B in $^; do ./$$B; done'
5997 clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \
5998 clean_tension_model_utils \
5999 clean_k_model_utils clean_latex clean_check_sawsim
6000 rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \
6001 $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \
6002 $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
6003 $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
6004 $(BUILD_DIR)/global.h ./gmon.out
6005 $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
6007 # Various builds of sawsim
6008 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
6009 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6010 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
6011 gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6012 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
6013 gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
6015 # Create the directories
6016 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
6019 # Copy over the source external to sawsim
6020 # Note: Cannot use, for example, `./build', because make eats the `./'
6021 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
6023 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6027 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
6031 ## Basic source generated with noweb
6032 # The central files sawsim.c and global.h...
6033 $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6035 $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6036 notangle -Rglobal.h $< > $@
6037 # ... and the modules
6038 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6039 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6040 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6041 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6042 $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
6043 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
6044 # Note: I use `_' as a space in my file names, but noweb doesn't like
6045 # that so we use `-' to name the noweb roots and substitute here.
6047 ## Building the unit-test programs
6049 $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
6050 notangle -Rchecks $< > $@
6051 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
6052 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
6053 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR)
6054 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
6055 clean_check_sawsim :
6056 rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
6057 # ... and the modules
6059 $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
6060 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \
6061 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \
6062 $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6063 $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
6064 $$(BUILD_DIR)/global.h | $(BIN_DIR)
6065 gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
6066 $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
6068 # todo: clean up the required modules to
6069 $(CHECK_BINS:%=clean_%) :
6070 rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
6072 # Cleaning up the modules
6074 $(SAWSIM_MODS:%=clean_%) :
6075 rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h)
6077 <<tension model utils makefile lines>>
6078 <<k model utils makefile lines>>
6079 <<latex makefile lines>>
6081 Makefile : $(SRC_DIR)/sawsim.nw
6082 notangle -Rmakefile $< | sed 's/ /\t/' > $@
6084 This is fairly self-explanatory. We compile a dynamically linked
6085 version ([[sawsim]]) and a statically linked version
6086 ([[sawsim_static]]). The static version is about 9 times as big, but
6087 it runs without needing dynamic GSL libraries to link to, so it's a
6088 better format for distributing to a cluster for parallel evaluation.
6092 \subsection{Hookean springs in series}
6093 \label{sec.math_hooke}
6096 The effective spring constant for $n$ springs $k_i$ in series is given by
6098 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6104 F &= k_i x_i &&\forall i \le n \\
6105 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
6106 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
6107 F &= k_1 x_1 = k_\text{eff} x \\
6108 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
6109 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6114 \addcontentsline{toc}{section}{References}
6115 \bibliography{sawsim}