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.
151 <<version definition>>=
152 #define VERSION "0.6"
158 sawsim - program for simulating single molecule mechanical unfolding.
159 Copyright (C) 2008-2009, William Trevor King
161 This program is free software; you can redistribute it and/or
162 modify it under the terms of the GNU General Public License as
163 published by the Free Software Foundation; either version 3 of the
164 License, or (at your option) any later version.
166 This program is distributed in the hope that it will be useful, but
167 WITHOUT ANY WARRANTY; without even the implied warranty of
168 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
169 See the GNU General Public License for more details.
171 You should have received a copy of the GNU General Public License
172 along with this program; if not, write to the Free Software
173 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
176 The author may be contacted at <wking@drexel.edu> on the Internet, or
177 write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
178 Philadelphia PA 19104, USA.
190 The unfolding system is stored as a chain of \emph{domains} (Figure
191 \ref{fig.domain_chain}). Domains can be folded, globular protein
192 domains, unfolded protein linkers, AFM cantilevers, or other
193 stretchable system components. The system is modeled as graph of
194 possible domain states with transitions between them (Figure
195 \ref{fig.domain_states}).
199 \subfloat[][]{\label{fig.domain_chain}
200 \begin{tikzpicture}[->,node distance=2cm,font=\footnotesize]
201 \tikzstyle{every state}=[fill,draw=red!50,very thick,fill=red!20]
202 \node[state] (A) {domain 1};
203 \node[state] (B) [below of=A] {domain 2};
204 \node[state] (C) [below of=B] {.~.~.};
205 \node[state] (D) [below of=C] {domain $N$};
206 \node (S) [below of=D] {Surface};
207 \node (E) [above of=A] {};
209 \path[-] (A) edge (B)
210 (B) edge node [right] {Tension} (C)
213 \path[->,green] (A) edge node [right,black] {Extension} (E);
217 \subfloat[][]{\label{fig.domain_states}
218 \begin{tikzpicture}[->,node distance=2.5cm,shorten <=1pt,shorten >=1pt,font=\footnotesize]
219 \tikzstyle{every state}=[fill,draw=blue!50,very thick,fill=blue!20]
220 \node[state] (A) {cantilever};
221 \node[state] (C) [below of=A] {transition};
222 \node[state] (B) [left of=C] {folded};
223 \node[state] (D) [right of=C] {unfolded};
225 \path (B) edge [bend left] node [above] {$k_1$} (C)
226 (C) edge [bend left] node [below] {$k_1'$} (B)
227 edge [bend left] node [above] {$k_2$} (D)
228 (D) edge [bend left] node [below] {$k_2'$} (C);
231 \caption{\subref{fig.domain_chain} Extending a chain of domains. One end
232 of the chain is fixed, while the other is extended at a constant
233 speed. The domains are coupled with rigid linkers, so the domains
234 themselves must stretch to accomodate the extension.
235 \subref{fig.domain_states} Each domain exists in a discrete state. At
236 each timestep, it may transition into another state following a
237 user-defined state matrix such as this one, showing a metastable
238 transition state and an explicit ``cantilever'' domain.}
242 Each domain contributes to the total tension, which depends on the
243 chain extension. The particular model for the tension as a function
244 of extension is handled generally with function pointers. So far the
245 following models have been implemented:
247 \item Null (Appendix \ref{sec.null_tension}),
248 \item Constant (Appendix \ref{sec.const_tension}),
249 \item Hookean spring (Appendix \ref{sec.hooke}),
250 \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
251 \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
254 The tension model and parameters used for a particular domain depend
255 on whether the domain's current state. The transition rates between
256 states are also handled generally with function pointers, with
259 \item Null (Appendix \ref{sec.null_k}),
260 \item Bell's model (Appendix \ref{sec.bell}),
261 \item Spring constant corrected Bell model (Appendix \ref{sec.kbell}),
262 \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
263 \item Kramers' saddle point model (Appendix \ref{sec.kramers}).
266 The unfolding of the chain domains is modeled in two stages. First
267 the tension of the chain is calculated. Then the tension (assumed to
268 be constant over the length of the chain, see Section
269 \ref{sec.timescales}), is applied to each folded domain, and
270 used to calculate the probability of that subdomain unfolding in the
271 next timestep $dt$. Then the time is stepped forward, and the process
275 int main(int argc, char **argv)
277 <<tension model getopt array definition>>
278 <<k model getopt array definition>>
279 list_t *domain_list=NULL; /* variables initialized in get_options() */
280 environment_t env={0};
281 double P, dt_max, v, xstep;
282 int num_folded=0, unfolding;
283 double x=0, dt, F; /* variables used in the simulation loop */
284 get_options(argc, argv, &P, &dt_max, &v, &xstep, &env, NUM_TENSION_MODELS,
285 tension_models, NUM_K_MODELS, k_models, &domain_list, &num_folded);
287 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Distance (m)\tForce (N)\tStiffness (N/m)\n");
288 if (flags & FLAG_OUTPUT_UNFOLDING_FORCES) printf("#Unfolding Force (N)\n");
289 while (num_folded > 0) {
290 F = find_tension(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, unfolding);
292 dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, v);
294 dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, 0);
295 unfolding=random_unfoldings(domain_list, F, dt, &env, &num_folded);
296 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness);
300 if (dt == dt_max) { /* step completed */
303 } else { /* still working on this step */
308 destroy_domain_list(domain_list);
312 @ The meat of the simulation is bundled into the three functions
313 [[find_tension]], [[determine_dt]], and [[random_unfoldings]].
314 [[find_tension]] is discussed in Section \ref{sec.find_tension},
315 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
316 [[random_unfoldings]] in Section \ref{sec.unfolding_rate}.
318 The stretched distance is extended in one of two modes depending on
319 [[xstep]]. If [[xstep == 0]], then the pulling is continuous, at
320 least within the limits of the inherent discretization of a time
321 stepped approach. At any rate, the timestep is limited by the maximum
322 allowed timestep [[dt_max]] and the maximum allowed unfolding
323 probability in a given step [[P]]. In the second mode [[xstep > 0]],
324 and the end to end distance increases in discrete steps of that
325 length. The time between steps is chosen to maintain an average
326 unfolding velocity [[v]]. This approach is less work to simulate than
327 smooth pulling and also closer to the reality of many experiments, but
328 it is practically impossible to treat analytically. With both pulling
329 styles implemented, the effects of the discretization can be easily
330 calculated, bridging the gap between theory and experiment.
332 Environmental parameters important in determining reaction rates and
333 tensions (e.g. temperature, pH) are stored in a single structure to
334 facilitate extension to more complicated models in the future.
335 <<environment definition>>=
336 typedef struct environment_struct {
345 <<environment definition>>
346 <<one dimensional function definition>>
347 <<create/destroy definitions>>
349 #endif /* GLOBAL_H */
351 Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
352 included multiple times.
354 \section{Simulation functions}
356 <<simulation functions>>=
360 <<does domain unfold>>
361 <<randomly unfold domains>>
365 \label{sec.find_tension}
367 Because the stretched system may be made up of several parts (folded
368 domains, unfolded domains, spring-like cantilever, \ldots), we will
369 assign the domains to models and groups. The models are used to
370 determine a domain's tension (Hookean spring, WLC, \ldots). Groups
371 will represent collections of domains which the model should treat as
372 a single entity. A domain may belong to different groups or models
373 depending on its state. For example, a domain might be modeled as a
374 freely-jointed chain when it is folded and as a worm-like chain when
375 it is unfolded. The domains are assumed to be commutative, so
376 ordering is ignored. The interactions between the groups are assumed
377 to be linear, but the interactions between domains of the same group
378 need not be. This allows for non-linear group models such as th
379 worm-like or freely-jointed chains. Each model has a tension handler
380 function, which gives the tension $F$ needed to stretch a given group
381 of domains a distance $x$.
383 To understand the purpose of groups, consider a chain of proteins
384 where two folded proteins $A$ and $B$ are modeled as WLCs with
385 persistence length $p_f$ and two unfolded proteins $C$ and $D$ are
386 modeled as WLCs with persistence length $p_u$. The proteins are
387 attached to a cantilever $E$ qof spring constant $κ$. The string
388 would get split into three lists:
390 \begin{tabular}{llll}
391 Model & Group & List & Tension \\
392 WLC & 0 & $[A,B]$ & $F_\text{WLC}(x, p_f, L_A+L_B)$ \\
393 WLC & 1 & $[C,D]$ & $F_\text{WLC}(x, p_u, L_C+L_D)$ \\
394 Hooke & 0 & $[E]$ & $F_\text{Hooke}(x, κ)$ \\
397 Note that group numbers only matter \emph{within} model classes, since
398 grouping domains with seperate models doesn't make sense.
400 <<find tension definitions>>=
401 #define NUM_TENSION_MODELS 5
405 <<tension handler array global declaration>>=
406 extern one_dim_fn_t *tension_handlers[];
409 <<tension handler array global>>=
410 one_dim_fn_t *tension_handlers[] = {
412 &const_tension_handler,
420 <<tension model global declarations>>=
421 <<tension handler array global declaration>>
424 <<tension handler types>>=
425 typedef struct tension_handler_data_struct {
426 /* int num_domains; */
427 /* domain_t *domains; */
431 } tension_handler_data_t;
435 After sorting the chain into separate groups $G_i$, with tension
436 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
437 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
438 \forall i,j$. Note that there may be several groups within each model
439 class. [[find_tension]] is basically a wrapper to feed proper input
440 into the [[tension_balance]] function. [[unfolding]] should be set to
441 0 if there were no unfoldings since the previous call to
442 [[find_tension]]. While were messing with the tension handlers, we
443 also grab a measure of the current linker stiffness for the
444 environmental variable, which is needed by some of the unfolding rate
445 models (e.g. the linker-stiffness corrected Bell model, Appendix
448 double find_tension(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list, environment_t *env, double x, int unfolding)
450 static int num_active_groups=0;
451 static one_dim_fn_t **per_group_handlers = NULL;
452 static void **per_group_params = NULL;
453 static double last_x;
454 tension_handler_data_t *tdata;
459 fprintf(stderr, "finding tension at distance %g%s\n", x, unfolding ? " after an unfolding" : "");
460 list_print(stderr, domain_list, "domain list");
463 if (unfolding != 0 || num_active_groups == 0) { /* setup statics */
464 /* free old statics */
465 if (per_group_handlers != NULL)
466 free(per_group_handlers);
467 if (per_group_params != NULL) {
468 for (i=0; i < num_active_groups; i++) {
469 assert(per_group_params[i] != NULL);
470 free(per_group_params[i]);
472 free(per_group_params);
475 /* get new statics */
476 get_active_groups(domain_list, num_tension_handlers, tension_handlers, &num_active_groups, &per_group_handlers, &per_group_params);
478 /* fill in the group handlers and params */
479 for (i=0; i<num_active_groups; i++) {
480 tdata = (tension_handler_data_t *) per_group_params[i];
482 fprintf(stderr, "active group %d of %d", i, num_active_groups);
483 list_print(stderr, tdata->group, " parameters");
486 /* tdata->persist continues unchanged */
489 } /* else roll with the current statics */
491 F = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_params, last_x, x);
498 fprintf(stderr, "finding linker stiffness at distance %g, tension %g\n", x, F);
504 if (dx == 0) { /* e.g. if x == 0 */
508 Flow = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_params, last_x, xlow);
509 env->stiffness = (F-Flow)/dx;
512 fprintf(stderr, " stiffness (%g-%g)/%g = %g\t(dx = %g = %g-%g)\n", F, Flow, dx, env->stiffness, dx, x, xlow);
520 @ For the moment, we will restrict our group boundaries to a single
521 dimension, so $\sum_i x_i = x$, our total extension (it is this
522 restriction that causes the groups to interact linearly). We'll also
523 restrict our extensions to all be positive. With these restrictions,
524 the problem of balancing the tensions reduces to solving a system of
525 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
526 the number of active groups. In general this can be fairly
527 complicated, but without much loss of practicality we can restrict
528 ourselves to strictly monotonically increasing, non-negative tension
529 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
530 simpler. For example, it guarantees the existence of a unique, real
531 solution for finite forces. See Appendix \ref{app.tension_balance}
532 for the tension balancing implementation.
534 Each group must define a parameter structure if the tension model is
536 a function to create the parameter structure,
537 a function to destroy the parameter structure, and
538 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as it's data argument.
539 The tension-specific parameter structures are contained in the domain
540 group field. See the Section \ref{app.model_selection} for a
541 discussion on the command line framework. See the worm-like chain
542 implementation for an example (Section \ref{sec.wlc}).
544 The major limitation of our approach is that all of our tension models
545 are Markovian (memory-less), which is not really a problem for the
546 chains (see \citet{evans99} for WLC thermalization rate calculations).
548 \subsection{Unfolding rate}
549 \label{sec.unfolding_rate}
551 Each folded domain is modeled as unimolecular, first order reaction
553 \text{Folded} \xrightarrow{k} % in tex: X atop Y
556 With the rate constant $k$ defined by
558 \frac{d[\text{Folded}]}{dt} = -k [\text{Folded}].
560 So the probability of a given protein unfolding in a short time $dt$
566 <<does domain unfold>>=
567 int domain_unfolds(double F, double dt, environment_t *env, domain_t *domain)
568 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to a folded domain */
570 k = accel_k(domain->k, F, env, domain->k_params);
571 //(*domain->k)(F, env, domain->k_params);
572 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
573 return happens(k*dt); /* dice roll for prob. k*dt event */
575 @ [[happens]] is a random decision making function defined in Appendix \ref{app.utils}.
577 <<randomly unfold domains>>=
578 int random_unfoldings(list_t *domain_list, double tension,
579 double dt, environment_t *env,
580 int *num_folded_domains)
581 { /* returns 1 if there was an unfolding and 0 otherwise */
582 while (domain_list != NULL) {
583 if (D(domain_list)->state == DS_FOLDED
584 && domain_unfolds(tension, dt, env, D(domain_list))) {
585 if (flags & FLAG_OUTPUT_UNFOLDING_FORCES)
586 fprintf(stdout, "%g\n", tension);
587 D(domain_list)->state = DS_UNFOLDED;
588 (*num_folded_domains)--;
589 return 1; /* our one domain has unfolded, stop looking for others */
591 domain_list = domain_list->next;
597 \subsection{Timescales}
598 \label{sec.timescales}
600 The simulation assumes that chain equilibration is instantaneous
601 relative to the simulation timestep $dt$, so that tension is uniform
602 along the chain. The quality of this assumption depends on your
603 particular chain. For example, a damped spring thermalizes on a
604 timescale of order $\tau = 1/\gamma$ where the damping ratio $\gamma
605 \equiv \omega_0(-\zeta \pm \sqrt{\zeta^2-1}$, the natural angular
606 frequency $omega_0 \equiv \sqrt{k/m}$, $\zeta \equiv b/2\sqrt{km}$,
607 and $b$ sets the drag $F_b = -bv$. For our cantilevers $\tau$ is
608 on the order of milliseconds, which is longer than a timestep.
609 However, the approximation is still reasonable, because there is
610 rarely unfolding during the cantilever return. You could, of course,
611 take cantilever, WLC, etc.\ response times into effect, but that
612 would significantly complicate a simulation that seems to work
613 well enough without bothering :p. For WLC and FJC relaxation times,
614 see the Appendix of \citet{evans99} and \citet{kroy07}.
616 Besides assuming our timestep is much greater than equilibration
617 times, we also force it to be short enough so that only one domain may
618 unfold in a given timestep. For an unfolding timescale of $dt_u =
619 1/k$ we adapt our timesteps to keep $dt \ll dt_u$, so the probability
620 of two domains unfolding in the same timestep is negligible. This
621 approach breaks down as the adaptive timestep scheme approaches $dt
622 \sim dt_u$, but $dt_u \sim 1\U{$\mu$s}$ for Ig27-like proteins
623 \citep{klimov00}, so this shouldn't be a problem. To reassure
624 yourself, you can ask the simulator to print the smallest timestep
625 that was used with TODO.
627 \subsection{Adaptive timesteps}
628 \label{sec.adaptive_dt}
630 We'd like to pick $dt$ so the probability of changing state
631 (e.g. unfolding) in the next timestep is small. If we don't adapt our
632 timestep, we also risk breaking our approximation that $F(x' \in [x,
633 x+v \cdot dt]) = F(x)$ or that only one domain changes state in a
634 given timestep. The simulation would have been accurate for
635 sufficiently small fixed timesteps, but adaptive timesteps will allow
636 us to move through low tension regions in fewer steps, leading to a
637 more efficient simulation.
639 Consider the two state folded $\rightarrow$ unfolded transition.
640 Because $F(x)$ is monotonically increasing between unfoldings,
641 excessively large timesteps will lead to erroneously small unfolding
642 rates (an thus, higher average unfolding force).
644 The actual adaptive timestep implementation is not particularly
645 interesting, since we are only required to reduce $dt$ to somewhere
646 below a set threshold, so I've removed it to Appendix
647 \ref{app.adaptive_dt}.
653 The models provide the physics of the simulation, but the simulation
654 allows interchangeable models, and we are currently focusing on the
655 simulation itself, so we remove the actual model implementation to the
658 The tension models are defined in Appendix \ref{sec.tension_models}
659 and unfolding models are defined in Appendix \ref{sec.k_models}.
662 #define k_B 1.3806503e-23 /* J/K */
666 \section{Command line interface}
668 <<get option functions>>=
670 <<index tension model>>
676 \subsection{Model selection}
677 \label{app.model_selection}
679 The main difficulty with the command line interface in \prog\ is
680 developing an intuitive method for accessing the possibly large number
681 of available models. We'll treat this generally by defining an array
682 of available models, containing their important parameters.
684 <<get options definitions>>=
685 <<model getopt definitions>>
688 <<create/destroy definitions>>=
689 typedef void *create_data_func_t(char **param_strings);
690 typedef void destroy_data_func_t(void *param_struct);
693 <<model getopt definitions>>=
694 <<tension model getopt definitions>>
695 <<k model getopt definitions>>
699 \subsubsection{Tension}
701 To access the various tension models from the command line, we define
702 a structure that contains all the neccessary information about the
704 <<tension model getopt definitions>>=
705 typedef struct tension_model_getopt_struct {
706 one_dim_fn_t *handler; /* fn implementing the model on a group */
707 char *name; /* name string identifying the model */
708 char *description; /* a brief description of the model */
709 int num_params; /* number of extra params for the model */
710 char **param_descriptions; /* descriptions of extra params */
711 char *params; /* default values for the extra params */
712 create_data_func_t *creator; /* param-string -> model-data-struct fn */
713 destroy_data_func_t *destructor; /* fn to free the model data structure */
714 } tension_model_getopt_t;
717 <<tension model getopt array definition>>=
718 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
719 <<null tension model getopt>>,
720 <<constant tension model getopt>>,
721 <<hooke tension model getopt>>,
722 <<worm-like chain tension model getopt>>,
723 <<freely-jointed chain tension model getopt>>
727 \subsubsection{Unfolding rate}
729 <<k model getopt definitions>>=
730 #define NUM_K_MODELS 6
732 typedef struct k_model_getopt_struct {
737 char **param_descriptions;
739 create_data_func_t *creator;
740 destroy_data_func_t *destructor;
744 <<k model getopt array definition>>=
745 k_model_getopt_t k_models[NUM_K_MODELS] = {
746 <<null k model getopt>>,
747 <<const k model getopt>>,
748 <<bell k model getopt>>,
749 <<kbell k model getopt>>,
750 <<kramers k model getopt>>,
751 <<kramers integ k model getopt>>
758 void help(char *prog_name, double P, double dt_max, double v, double xstep,
760 int n_tension_models, tension_model_getopt_t *tension_models,
761 int folded_tension_model, int unfolded_tension_model,
762 int n_k_models, k_model_getopt_t *k_models,
766 printf("usage: %s [options]\n", prog_name);
767 printf("Version %s\n\n", VERSION);
768 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
769 printf("Simulation options:\n");
770 printf("-P P\tTarget probability for dt (currently %g)\n", P);
771 printf("-t dt\tMaximum allowed timestep dt (currently %g)\n", dt_max);
772 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
773 printf("-x dx\tPulling step size (currently %g m)\n", xstep);
774 printf("\tWhen dx = 0, the pulling is continuous\n");
775 printf("\tWhen dx > 0, the pulling occurs in discrete steps\n");
776 printf("Environmental options:\n");
777 printf("-T T\tTemperature T (currently %g K)\n", env->T);
778 printf("-C T\tYou can also set the temperature T in Celsius\n");
779 printf("Model options:\n");
780 printf("The domains exist in either the folded or unfolded state\n");
781 printf("The following options change the default behavior in each state.\n");
782 printf("-m model[,group]\tFolded tension model (currently %s)\n",
783 tension_models[folded_tension_model].name);
784 printf("-a args\tFolded tension model argument string (currently %s)\n",
785 tension_models[folded_tension_model].params);
786 printf("-M model[,group]\tUnfolded tension model (currently %s)\n",
787 tension_models[unfolded_tension_model].name);
788 printf("-A args\tUnfolded tension model argument string (currently %s)\n",
789 tension_models[unfolded_tension_model].params);
790 printf("The following options change the unfolding rate.\n");
791 printf("-k model\tTransition rate model (currently %s)\n",
792 k_models[k_model].name);
793 printf("-K args\tTransition rate model argument string (currently %s)\n",
794 k_models[k_model].params);
795 printf("Domain creation options:\n");
796 printf("Once you've set up the appropriate default models, you need to add the domains\n");
797 printf("-F n\tAdd n folded domains with the current models\n");
798 printf("-U n\tAdd n unfolded domains with the current models\n");
799 printf("Output mode options:\n");
800 printf("There are two output modes. In standard mode, only the unfolding\n");
801 printf("events are printed. For example:\n");
802 printf(" #Unfolding Force (N)\n");
803 printf(" 123.456e-12\n");
805 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
806 printf(" #Distance (m)\tForce (N)\tStiffness (N/m)\n");
807 printf(" 0.001\t0.0023\n");
809 printf("-V\tChange output to verbose mode\n");
810 printf("-h\tPrint this help and exit\n");
812 printf("Tension models:\n");
813 for (i=0; i<n_tension_models; i++) {
814 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
815 for (j=0; j < tension_models[i].num_params ; j++)
816 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
817 printf("\t\tdefaults: %s\n", tension_models[i].params);
819 printf("Unfolding rate models:\n");
820 for (i=0; i<n_k_models; i++) {
821 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
822 for (j=0; j < k_models[i].num_params ; j++)
823 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
824 printf("\t\tdefaults: %s\n", k_models[i].params);
826 printf("Examples:\n");
827 printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded\n");
828 printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8\n", prog_name);
829 printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded\n");
830 printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K5e-5,.225e-9 -mwlc -a1e-8,4e-10 -Mwlc,1 -A0.4e-9,28.1e-9 -F8\n", prog_name);
834 \subsection{Get options}
837 void get_options(int argc, char **argv,
838 double *pP, double *pDt_max, double *pV, double *pXstep,
840 int n_tension_models, tension_model_getopt_t *tension_models,
841 int n_k_models, k_model_getopt_t *k_models,
842 list_t **pDomain_list, int *pNum_folded)
844 char *prog_name = NULL;
845 char c, options[] = "P:t:v:x:T:C:m:a:M:A:k:K:F:U:Vh";
846 int ftension_model=0, utension_model=0, k_model=0;
847 char *ftension_params=NULL, *utension_params=NULL;
848 int i, n, ftension_group=0, utension_group=0;
852 extern int optind, optopt, opterr;
857 for (i=0; i<n_tension_models; i++) {
858 fprintf(stderr, "tension model %i %s:\t model's handler %p\n",
859 i, tension_models[i].name, tension_models[i].handler);
860 assert(tension_models[i].handler == tension_handlers[i]);
865 flags = FLAG_OUTPUT_UNFOLDING_FORCES;
867 *pP = 1e-3; /* % pop per s */
868 *pDt_max = 0.001; /* s */
869 *pV = 1e-6; /* m/s */
870 *pXstep = 0.0; /* m */
871 env->T = 300.0; /* K */
873 while ((c=getopt(argc, argv, options)) != -1) {
875 case 'P': *pP = atof(optarg); break;
876 case 't': *pDt_max = atof(optarg); break;
877 case 'v': *pV = atof(optarg); break;
878 case 'x': *pXstep = atof(optarg); break;
879 case 'T': env->T = atof(optarg); break;
880 case 'C': env->T = atof(optarg)+273.15; break;
882 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
883 assert(num_strings == 1 || num_strings == 2);
884 if (num_strings == 1)
887 ftension_group = atoi(strings[1]);
888 ftension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
889 ftension_params = tension_models[ftension_model].params;
890 free_parsed_list(num_strings, strings);
893 ftension_params = optarg;
896 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
897 assert(num_strings == 1 || num_strings == 2);
898 if (num_strings == 1)
901 utension_group = atoi(strings[1]);
902 utension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
903 utension_params = tension_models[utension_model].params;
904 free_parsed_list(num_strings, strings);
907 utension_params = optarg;
910 k_model = index_k_model(n_k_models, k_models, optarg);
913 k_models[k_model].params = optarg;
918 for (i=0; i<n; i++) {
919 push(pDomain_list, generate_domain(DS_FOLDED,
920 tension_models+ftension_model,
923 tension_models+utension_model,
926 k_models+k_model), 1);
933 for (i=0; i<n; i++) {
934 push(pDomain_list, generate_domain(DS_UNFOLDED,
935 tension_models+ftension_model,
938 tension_models+utension_model,
941 k_models+k_model), 1);
945 flags = FLAG_OUTPUT_FULL_CURVE;
948 fprintf(stderr, "unrecognized option '%c'\n", optopt);
949 /* fall through to default case */
951 help(prog_name, *pP, *pDt_max, *pV, *pXstep, env, n_tension_models, tension_models, ftension_model, utension_model, n_k_models, k_models, k_model);
955 /* check the arguments */
956 assert(*pDomain_list != NULL); /* no domains! */
957 assert(*pP > 0.0); assert(*pP < 1.0);
958 assert(*pDt_max > 0.0);
960 assert(env->T > 0.0);
965 <<index tension model>>=
966 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
969 for (i=0; i<n_models; i++)
970 if (STRMATCH(models[i].name, name))
973 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
981 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
984 for (i=0; i<n_models; i++)
985 if (STRMATCH(models[i].name, name))
988 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
995 <<initialize model definition>>=
996 /* requires int num_param_args and char **param_args in the current scope
998 * INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
999 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
1001 #define INIT_MODEL(role, model, param_string, param_pointer) \
1003 parse_list_string(param_string, SEP, '{', '}', \
1004 &num_param_args, ¶m_args); \
1005 if (num_param_args != model->num_params) { \
1007 "%s model %s expected %d params,\n", \
1008 role, model->name, model->num_params); \
1010 "not the %d params in '%s'\n", \
1011 num_param_args, param_string); \
1012 assert(num_param_args == model->num_params); \
1014 if (model->creator) \
1015 param_pointer = (*model->creator)(param_args); \
1017 param_pointer = NULL; \
1018 free_parsed_list(num_param_args, param_args); \
1022 <<generate domain>>=
1023 void *generate_domain(enum domain_state_t state,
1024 tension_model_getopt_t *folded_model,
1026 char *folded_param_string,
1027 tension_model_getopt_t *unfolded_model,
1029 char *unfolded_param_string,
1030 k_model_getopt_t *k_model)
1032 void *folded_params, *unfolded_params, *k_params;
1033 int num_param_args; /* for INIT_MODEL() */
1034 char **param_args; /* for INIT_MODEL() */
1037 fprintf(stderr, "generating %s ", state==DS_FOLDED ? "folded" : "unfolded");
1038 fprintf(stderr, "domain (k: \"%s\", \"%s\", f: \"%s\" (%p) group %d \"%s\", u: \"%s\" (%p) group %d \"%s\")\n",
1039 k_model->name, k_model->params,
1040 folded_model->name, folded_model->handler, folded_group, folded_param_string,
1041 unfolded_model->name, unfolded_model->handler, unfolded_group, unfolded_param_string);
1043 INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
1044 INIT_MODEL("unfolded", unfolded_model, unfolded_param_string, unfolded_params);
1045 INIT_MODEL("k", k_model, k_model->params, k_params);
1046 return create_domain(state,
1047 k_model->k, k_params, k_model->destructor,
1048 folded_model->handler, folded_group, folded_params, folded_model->destructor,
1049 unfolded_model->handler, unfolded_group, unfolded_params, unfolded_model->destructor
1056 \addcontentsline{toc}{section}{Appendicies}
1058 \section{sawsim.c details}
1062 The general layout of our simulation code is:
1072 We include [[math.h]], so don't forget to link to the libm with `-lm'.
1074 #include <assert.h> /* assert() */
1075 #include <stdlib.h> /* malloc(), free(), rand() */
1076 #include <stdio.h> /* fprintf(), stdout */
1077 #include <string.h> /* strlen, strtok() */
1078 #include <math.h> /* exp(), M_PI, sqrt() */
1079 #include <sys/types.h> /* pid_t (returned by getpid()) */
1080 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
1083 #include "tension_balance.h"
1084 #include "k_model.h"
1085 #include "tension_model.h"
1087 #include "accel_k.h"
1091 <<version definition>>
1092 <<flag definitions>>
1093 <<max/min definitions>>
1094 <<string comparison definition and globals>>
1095 <<initialize model definition>>
1096 <<get options definitions>>
1097 <<domain definitions>>
1106 <<create/destroy domain>>
1107 <<destroy domain list>>
1109 <<simulation functions>>
1110 <<boilerplate functions>>
1113 <<boilerplate functions>>=
1115 <<get option functions>>
1121 srand(getpid()*time(NULL)); /* seed rand() */
1125 <<flag definitions>>=
1126 /* in octal b/c of prefixed '0' */
1127 #define FLAG_OUTPUT_FULL_CURVE 01
1128 #define FLAG_OUTPUT_UNFOLDING_FORCES 02
1132 static unsigned long int flags = 0;
1135 \subsection{Utilities}
1138 <<max/min definitions>>=
1139 #define MAX(a,b) ((a)>(b) ? (a) : (b))
1140 #define MIN(a,b) ((a)<(b) ? (a) : (b))
1143 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
1144 <<string comparison definition and globals>>=
1145 // Check if two strings match, return 1 if they do
1146 static char *temp_string_A;
1147 static char *temp_string_B;
1148 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
1149 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
1150 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
1151 /* +1 to also compare the '\0' */
1154 We also define a macro for our [[check]] unit testing
1155 <<check relative error>>=
1156 #define CHECK_ERR(max_err, expected, received) \
1158 fail_unless( (received-expected)/expected < max_err, \
1159 "relative error %g >= %g in %s (Expected %g, received %g)", \
1160 (received-expected)/expected, max_err, #received, \
1161 expected, received); \
1162 fail_unless(-(received-expected)/expected < max_err, \
1163 "relative error %g >= %g in %s (Expected %g, received %g)", \
1164 -(received-expected)/expected, max_err, #received, \
1165 expected, received); \
1170 int happens(double probability)
1172 assert(probability >= 0.0); assert(probability <= 1.0);
1173 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*/
1178 \subsection{Adaptive timesteps}
1179 \label{app.adaptive_dt}
1181 $F(x)$ increases with $x$, possibly exploding, as in the worm-like chain model,
1182 so basing the timestep on the the unfolding probability at the current tension
1183 is dangerous, and we need to search for a $dt$ for which
1184 $P(F(x+v*dt)) < P_\text{target}$.
1185 There are two cases to consider.
1186 In the most common, no domains have unfolded since the last step,
1187 and we expect the next step to be slightly shorter than the previous one.
1188 In the less common, domains did unfold in the last step,
1189 and we expect the next step to be considerably longer than the previous one.
1191 double search_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers,
1192 list_t *domain_list,
1193 environment_t *env, double x,
1194 double target_prob, double max_dt, double v)
1195 { /* Returns the timestep dt in seconds for the current folded domain.
1196 * Takes a list of tension handlers, the list of domains,
1197 * a pointer env to the environmental data, a starting separation x in m,
1198 * a target_prob between 0 and 1,
1199 * max_dt in s, stretching velocity v in m/s.
1201 double F, k, dtCur, dtU, dtUCur, dtL, dt;
1203 /* get upper bound using the current position */
1204 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x, 0); /* BUG. repeated calculation */
1205 //printf("Start with x = %g (F = %g)\n", x, F);
1206 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1207 //printf("x %g\tF %g\tk %g\n", x, F, k);
1208 dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */
1210 //printf("overshot max_dt\n");
1213 if (v == 0) /* discrete stepping, so force is temporatily constant */
1216 /* set a lower bound on dt too */
1219 /* The dt determined above may produce illegitimate forces or ks.
1220 * Reduce the upper bound until we have valid ks. */
1222 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1223 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1226 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1228 //printf("Try for dt = %g (F = %g)\n", dt, F);
1229 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1230 /* returned k may be -1.0 */
1231 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1232 while (k == -1.0) { /* reduce step until we hit a valid k */
1234 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1235 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1236 //printf("Try for dt = %g (F = %g)\n", dt, F);
1237 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1238 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1240 assert(dtU > 1e-14); /* timestep to valid k too small */
1241 dtUCur = target_prob / k; /* safe timestep back from x+dtU */
1243 return dt; /* dtU is safe. */
1245 /* dtU wasn't safe, lets see what would be. */
1246 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1247 dt = (dtU + dtL) / 2.0;
1248 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1249 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1250 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1251 dtCur = target_prob / k;
1252 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1253 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1255 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1256 dtU = dt; /* ... stepping out only dtCur would be safe */
1259 break; /* dtCur = dt */
1261 return MAX(dtUCur, dtL);
1265 To determine $dt$ for an array of potentially different folded domains,
1266 we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
1269 double determine_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list,
1270 environment_t *env, double x,
1271 double target_prob, double dt_max, double v)
1272 { /* Returns the timestep dt in seconds.
1273 * Takes the list of folded domains, target_prob between 0 and 1,
1274 * F in N, and T in K. */
1275 double dt=dt_max, new_dt;
1276 assert(target_prob > 0.0); assert(target_prob < 1.0);
1277 assert(dt_max > 0.0);
1279 while (domain_list != NULL) {
1280 if (D(domain_list)->state == DS_FOLDED) {
1281 new_dt = search_dt(num_tension_handlers, tension_handlers, domain_list, env, x, target_prob, dt, v);
1282 dt = MIN(dt, new_dt);
1284 domain_list = domain_list->next;
1290 \subsection{Domain data}
1292 Currently domains exist in two states, folded and unfolded, and the
1293 only allowed transitions are folded $\rightarrow$ unfolded. Of
1294 course, it wouldn't be too complicated to extent this to a multi-state
1295 system, with an array containing the domains group for each possible
1296 state, and a matrix of transition-rate-calculating functions.
1297 However, at this point such generality seems unnecessary at this
1299 <<domain definitions>>=
1300 enum domain_state_t {DS_FOLDED,
1304 typedef struct domain_struct {
1305 enum domain_state_t state;
1306 one_dim_fn_t *folded_handler;
1308 one_dim_fn_t *unfolded_handler;
1310 k_func_t *k; /* function returning unfolding rate */
1311 void *folded_params; /* pointer to folded parameters */
1312 void *unfolded_params; /* pointer to unfolded parameters */
1313 void *k_params; /* pointer to k parameters */
1314 destroy_data_func_t *destroy_folded;
1315 destroy_data_func_t *destroy_unfolded;
1316 destroy_data_func_t *destroy_k;
1319 /* get the domain data for the current list node */
1320 #define D(list) ((domain_t *)(list)->d)
1321 /* get the tension params for the current list node */
1322 #define D_TP(list) (((domain_t *)(list)->d)->state == DS_FOLDED \
1323 ? ((domain_t *)(list)->d)->folded_params \
1324 : ((domain_t *)(list)->d)->unfolded_params)
1326 [[k]] is a pointer to the function determining the unfolding rate for a given tension.
1327 [[folded_params]] is a pointer to the parameters used by the function pointed to by [[k]].
1328 [[unfolded_params]] is a pointer to the parameters used by the group-appropriate handler function when determining the tension.
1329 The [[destroy_*]] pointers point to functions for freeing the memory [[*_params]].
1330 We store them with the domain data so that [[destroy_domain]] doesn't have to know which type of domain it's cleaning up after.
1332 [[create_]] and [[destroy_domain]] are simple wrappers around [[malloc]] and [[free]].
1333 <<create/destroy domain>>=
1334 domain_t *create_domain(enum domain_state_t state,
1337 destroy_data_func_t *destroy_k,
1338 one_dim_fn_t *folded_handler,
1340 void *folded_params,
1341 destroy_data_func_t *destroy_folded,
1342 one_dim_fn_t *unfolded_handler,
1344 void *unfolded_params,
1345 destroy_data_func_t *destroy_unfolded)
1347 domain_t *ret = (domain_t *)malloc(sizeof(domain_t));
1348 assert(ret != NULL);
1349 if (state == DS_FOLDED) {
1350 assert(k != NULL); /* the pointer points somewhere valid */
1351 assert(*k != NULL); /* and there is something useful there */
1353 assert(state == DS_UNFOLDED);
1355 ret->folded_handler = folded_handler;
1356 ret->folded_group = folded_group;
1357 ret->unfolded_handler = unfolded_handler;
1358 ret->unfolded_group = unfolded_group;
1360 ret->folded_params = folded_params;
1361 ret->unfolded_params = unfolded_params;
1362 ret->k_params = k_params;
1363 ret->destroy_folded = destroy_folded;
1364 ret->destroy_unfolded = destroy_unfolded;
1365 ret->destroy_k = destroy_k;
1367 fprintf(stderr, "generated domain %p. state %d\n", ret, ret->state);
1372 void destroy_domain(domain_t *domain)
1375 //printf("domain %p & %p\n", *domain, domain);
1376 if (domain->destroy_folded)
1377 (*domain->destroy_folded)(domain->folded_params);
1378 if (domain->destroy_unfolded)
1379 (*domain->destroy_unfolded)(domain->unfolded_params);
1380 if (domain->destroy_k)
1381 (*domain->destroy_k)(domain->k_params);
1387 <<destroy domain list>>=
1388 void destroy_domain_list(list_t *domain_list)
1390 domain_list = head(domain_list);
1391 while (domain_list != NULL)
1392 destroy_domain((domain_t *) pop(&domain_list));
1396 \subsection{Domain model and group handling}
1398 <<group functions>>=
1399 <<get tension handler>>
1401 <<int comparison function>>
1402 <<find possible groups>>
1405 <<get active groups>>
1408 <<get tension handler>>=
1409 one_dim_fn_t *get_tension_handler(domain_t *domain)
1411 if (domain->state == DS_FOLDED)
1412 return domain->folded_handler;
1414 if (domain->state != DS_UNFOLDED)
1415 fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
1416 assert(domain->state == DS_UNFOLDED);
1417 return domain->unfolded_handler;
1423 int get_group(domain_t *domain)
1425 if (domain->state == DS_FOLDED)
1426 return domain->folded_group;
1428 if (domain->state != DS_UNFOLDED)
1429 fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
1430 assert(domain->state == DS_UNFOLDED);
1431 return domain->unfolded_group;
1436 We already know all possible tension classes, since they are all
1437 hardcoded into \prog. However, there may be any number of possible
1438 groups. We define a function [[find_possible_groups]] to search for
1439 possible (not neccessarily active) groups. It can search for
1440 subgroups of a single [[handler]], or by setting [[handler]] to
1441 [[NULL]], subgroups of \emph{any} handler. It returns a list of group
1443 <<find possible groups>>=
1444 list_t *find_possible_groups(list_t *list, one_dim_fn_t *handler) {
1447 while (list != NULL) {
1448 if (handler == NULL || get_tension_handler(D(list)) == handler) {
1449 push(&ret, &D(list)->folded_group, 1);
1450 push(&ret, &D(list)->unfolded_group, 1);
1456 return ret; /* no groups with this handler, no processing needed */
1458 /* sort the ret list, and remove duplicates */
1459 sort(&ret, &int_comparison);
1460 unique(&ret, &int_comparison);
1462 fprintf(stderr, "handler %p possible groups:", handler);
1464 while (list != NULL) {
1465 fprintf(stderr, " %d", *((int *)list->d));
1468 fprintf(stderr, "\n");
1474 Search a [[list]] of domains, and determine whether a particular model
1475 class and group number combination is active.
1476 <<is group active>>=
1477 int is_group_active(list_t *list, one_dim_fn_t *handler, int group)
1480 while (list != NULL) {
1481 if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group)
1489 Search a [[list]] of domains, and return all domains matching a
1490 particular model class and group number.
1492 list_t *get_group_list(list_t *list, one_dim_fn_t *handler, int group)
1496 while (list != NULL) {
1497 if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group) {
1498 push(&ret, D_TP(list), 1); /* add a pointer to the appropriate tension parameters to our new list. */
1500 fprintf(stderr, "push domain %p %s tension parameters %p onto active group %p %d list %p\n", D(list), D(list)->state == DS_FOLDED ? "folded" : "unfolded", D_TP(list), handler, group, ret);
1508 Because all the node data in lists returned by [[get_group_list]] is
1509 also in the main domain list, you shouldn't destroy the node data
1510 popped off when destroying the group lists. It will all get cleaned
1511 up when the main domain list is destroyed.
1513 Put all this together to scan through a list of domains and construct
1514 an array of all the active groups.
1515 <<get active groups>>=
1516 void get_active_groups(list_t *list,
1517 int num_tension_handlers, one_dim_fn_t **tension_handlers,
1518 int *pNum_active_groups, one_dim_fn_t ***pPer_group_handlers, void ***pActive_groups)
1520 void **active_groups=NULL;
1521 one_dim_fn_t *handler, **per_group_handlers=NULL;
1522 int i, num_active_groups, *group;
1523 list_t *possible_group_numbers, *active_group_list, *active_groups_list=NULL, *per_group_handlers_list=NULL;
1525 /* get all the active groups in a list */
1526 for (i=0; i<num_tension_handlers; i++) {
1527 handler = tension_handlers[i];
1528 if (handler == NULL) continue; /* tensionless handler */
1529 possible_group_numbers = head(find_possible_groups(list, handler));
1530 while (possible_group_numbers != NULL) {
1531 group = (int *)pop(&possible_group_numbers);
1532 if (is_group_active(list, handler, *group) == 1) {
1533 active_group_list = get_group_list(list, handler, *group);
1534 push(&active_groups_list, active_group_list, 1);
1535 push(&per_group_handlers_list, handler, 1);
1537 fprintf(stderr, "active group %p for %p %d headed by tension parameters %p\n", active_group_list, handler, *group, head(active_group_list)->d);
1538 list_print(stderr, active_group_list, "active group");
1544 /* allocate the array we'll be returning */
1545 num_active_groups = list_length(active_groups_list);
1546 active_groups = (void **) malloc(sizeof(void *)*num_active_groups);
1547 assert(active_groups != NULL);
1548 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1549 assert(per_group_handlers != NULL);
1551 /* populate the array we'll be returning */
1552 active_groups_list = head(active_groups_list);
1553 for (i=0; i<num_active_groups; i++) {
1554 handler = pop(&per_group_handlers_list);
1555 per_group_handlers[i] = handler;
1557 active_groups[i] = (void *) malloc(sizeof(tension_handler_data_t));
1558 assert(active_groups[i] != NULL);
1559 active_group_list = pop(&active_groups_list);
1560 ((tension_handler_data_t *)active_groups[i])->group = active_group_list;
1561 ((tension_handler_data_t *)active_groups[i])->env = NULL;
1562 ((tension_handler_data_t *)active_groups[i])->persist = NULL;
1564 assert(active_groups_list == NULL);
1565 assert(per_group_handlers_list == NULL);
1567 *pNum_active_groups = num_active_groups;
1568 *pPer_group_handlers = per_group_handlers;
1569 *pActive_groups = active_groups;
1574 \section{String parsing}
1576 For handling command line arguments, we need parse delimited strings (e.g. [["param1,param2,param3..."]]).
1577 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1578 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1579 need to take care of parsing those parameters themselves.
1580 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1586 <<parse definitions>>
1587 <<parse declarations>>
1591 <<parse module makefile lines>>=
1592 NW_SAWSIM_MODS += parse
1593 CHECK_BINS += check_parse
1597 <<parse definitions>>=
1598 #define SEP ',' /* argument separator character */
1601 <<parse declarations>>=
1602 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1603 int *num, char ***string_array);
1604 extern void free_parsed_list(int num, char **string_array);
1607 [[parse_list_string]] allocates memory, don't forget to free it
1608 afterward with [[free_parsed_list]]. It does not alter the original.
1610 The string may start off with a [[deeper]] character
1611 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1612 [[x,y]], where the pointer is one character in on the copied string.
1613 However, when we go to free the memory, we need a pointer to the
1614 beginning of the string. In order to accommodate this for a string
1615 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1616 the first $N$ elements point to the separated arguments, and let the
1617 last element point to the start of the copied string regardless of
1619 <<parse delimited list string functions>>=
1620 /* TODO, split out into parse.hc */
1621 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1624 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1625 if (string[i] == deeper) {depth++;}
1626 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1632 void parse_list_string(char *string, char sep, char deeper, char shallower,
1633 int *num, char ***string_array)
1635 char *str=NULL, **ret=NULL;
1637 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
1639 *string_array = NULL;
1642 /* make a copy of the string, so we don't change the original */
1643 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1644 assert(str != NULL);
1645 strcpy(str, string); /* we know str is long enough */
1646 /* count the number of regions, so we can allocate pointers to them */
1649 n++; i++; /* move on to next argument */
1650 i += next_delim_index(str+i, sep, deeper, shallower);
1651 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1652 } while (str[i] != '\0');
1653 ret = (char **)malloc(sizeof(char *)*(n+1));
1654 assert(ret != NULL);
1655 /* replace the separators with '\0' & assign pointers */
1656 ret[n] = str; /* point to the front of the copied string */
1659 for(i=1; i<n; i++) {
1660 j += next_delim_index(str+j, sep, deeper, shallower);
1662 ret[i] = str+j; /* point to the separated arguments */
1664 /* strip off bounding braces, if any */
1665 for(i=0; i<n; i++) {
1666 if (ret[i][0] == deeper) {
1670 if (ret[i][strlen(ret[i])-1] == shallower)
1671 ret[i][strlen(ret[i])-1] = '\0';
1674 *string_array = ret;
1677 void free_parsed_list(int num, char **string_array)
1680 /* frees all items, since they were allocated as a single string*/
1681 free(string_array[num]);
1682 /* and free the array of pointers */
1688 \subsection{Parsing implementation}
1694 <<parse delimited list string functions>>
1698 #include <assert.h> /* assert() */
1699 #include <stdlib.h> /* NULL */
1700 #include <stdio.h> /* fprintf(), stdout *//*!!*/
1701 #include <string.h> /* strlen() */
1705 \subsection{Parsing unit tests}
1707 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1709 <<parse check includes>>
1710 <<string comparison definition and globals>>
1711 <<check relative error>>
1712 <<parse test suite>>
1713 <<main check program>>
1716 <<parse check includes>>=
1717 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
1718 #include <stdio.h> /* printf() */
1719 #include <assert.h> /* assert() */
1720 #include <string.h> /* strlen() */
1725 <<parse test suite>>=
1726 <<parse list string tests>>
1727 <<parse suite definition>>
1730 <<parse suite definition>>=
1731 Suite *test_suite (void)
1733 Suite *s = suite_create ("k model");
1734 <<parse list string test case defs>>
1736 <<parse list string test case add>>
1741 <<parse list string tests>>=
1744 START_TEST(test_next_delim_index)
1746 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
1747 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
1748 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
1749 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
1750 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
1754 START_TEST(test_parse_list_null)
1758 parse_list_string(NULL, SEP, '{', '}',
1759 &num_param_args, ¶m_args);
1760 fail_unless(num_param_args == 0, NULL);
1761 fail_unless(param_args == NULL, NULL);
1764 START_TEST(test_parse_list_single_simple)
1768 parse_list_string("arg", SEP, '{', '}',
1769 &num_param_args, ¶m_args);
1770 fail_unless(num_param_args == 1, NULL);
1771 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
1774 START_TEST(test_parse_list_single_compound)
1778 parse_list_string("{x,y,z}", SEP, '{', '}',
1779 &num_param_args, ¶m_args);
1780 fail_unless(num_param_args == 1, NULL);
1781 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
1784 START_TEST(test_parse_list_double_simple)
1788 parse_list_string("abc,def", SEP, '{', '}',
1789 &num_param_args, ¶m_args);
1790 fail_unless(num_param_args == 2, NULL);
1791 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
1792 fail_unless(STRMATCH(param_args[1],"def"), NULL);
1796 <<parse list string test case defs>>=
1797 TCase *tc_parse_list_string = tcase_create("parse list string");
1799 <<parse list string test case add>>=
1800 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
1801 tcase_add_test(tc_parse_list_string, test_parse_list_null);
1802 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
1803 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
1804 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
1805 suite_add_tcase(s, tc_parse_list_string);
1809 \section{Unit tests}
1811 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1818 <<check relative error>>
1821 <<main check program>>
1833 <<determine dt tests>>
1835 <<does domain unfold tests>>
1836 <<randomly unfold domains tests>>
1837 <<suite definition>>
1840 <<suite definition>>=
1841 Suite *test_suite (void)
1843 Suite *s = suite_create ("sawsim");
1844 <<F test case defs>>
1845 <<determine dt test case defs>>
1846 <<happens test case defs>>
1847 <<does domain unfold test case defs>>
1848 <<randomly unfold domains test case defs>>
1851 <<determine dt test case add>>
1852 <<happens test case add>>
1853 <<does domain unfold test case add>>
1854 <<randomly unfold domains test case add>>
1857 tcase_add_checked_fixture(tc_strip_address,
1858 setup_strip_address,
1859 teardown_strip_address);
1865 <<main check program>>=
1869 Suite *s = test_suite();
1870 SRunner *sr = srunner_create(s);
1871 srunner_run_all(sr, CK_ENV);
1872 number_failed = srunner_ntests_failed(sr);
1874 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
1878 \subsection{F tests}
1880 <<worm-like chain tests>>
1882 <<F test case defs>>=
1883 <<worm-like chain test case def>>
1885 <<F test case add>>=
1886 <<worm-like chain test case add>>
1889 <<worm-like chain tests>>=
1890 extern double wlc(double x, double T, double p, double L);
1891 START_TEST(test_wlc_at_zero)
1893 double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
1894 fail_unless(abs(wlc(x, T, p, L)) < lim, \
1895 "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
1896 x, T, p, L, abs(wlc(x, T, p, L)), lim);
1900 START_TEST(test_wlc_at_half)
1902 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
1903 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
1904 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
1906 * linear term = x/L = 0.5
1907 * nonlinear + linear = 0.75 + 0.5 = 1.25
1908 * wlc = 10e21*1.25 = 12.5e21
1910 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
1911 "wlc(%g, %g, %g, %g) = %g != %g",
1912 x, T, p, L, wlc(x, T, p, L), 12.5e21);
1916 <<worm-like chain test case def>>=
1917 TCase *tc_wlc = tcase_create("WLC");
1920 <<worm-like chain test case add>>=
1921 tcase_add_test(tc_wlc, test_wlc_at_zero);
1922 tcase_add_test(tc_wlc, test_wlc_at_half);
1923 suite_add_tcase(s, tc_wlc);
1926 \subsection{Model tests}
1928 Check the searching with [[linear_k]].
1929 Check overwhelming force treatment with the heavyside-esque step [[?]].
1930 <<determine dt tests>>=
1931 double linear_k(double F, environment_t *env, void *params)
1933 double Fnot = *(double *)params;
1937 START_TEST(test_determine_dt_linear_k)
1940 double dt_max=3.0, Fnot=3.0;
1941 double F[]={0,1,2,3,4,5,6};
1942 domain_t dom; /* use both parts at once for folded/unfolded */
1946 dom->next = dom->prev = NULL;
1947 dom->k_func_t = linear_k;
1948 dom->folded_params = &Fnot;
1949 dom->unfolded_params = !!!!!!!!!
1950 dom->destroy_folded = dom->destroy_unfolded = NULL;
1951 for( i=0; i < sizeof(F)/sizeof(double); i++) {
1952 fail_unless(determine_dt(folded, unfolded, targ_p, dt_max, v, x, T)==1, NULL);
1958 <<determine dt test case defs>>=
1959 TCase *tc_determine_dt = tcase_create("Determine dt");
1961 <<determine dt test case add>>=
1962 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
1963 suite_add_tcase(s, tc_determine_dt);
1968 <<happens test case defs>>=
1970 <<happens test case add>>=
1973 <<does domain unfold tests>>=
1975 <<does domain unfold test case defs>>=
1977 <<does domain unfold test case add>>=
1980 <<randomly unfold domains tests>>=
1982 <<randomly unfold domains test case defs>>=
1984 <<randomly unfold domains test case add>>=
1988 \section{Balancing group extensions}
1989 \label{app.tension_balance}
1991 For a given total extension $x$ (of the piezo), the various domain
1992 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
1993 amounts, and we need to tweak the portion that each extends to balance
1994 the tension amongst the active groups. Since the tension balancing is
1995 separable from the bulk of the simulation, we break it out into a
1996 separate module. The interface is defined in [[tension_balance.h]],
1997 the implementation in [[tension_balance.c]], and the unit testing in
1998 [[check_tension_balance.c]]
2000 <<tension-balance.h>>=
2001 #ifndef TENSION_BALANCE_H
2002 #define TENSION_BALANCE_H
2004 <<tension balance function declaration>>
2005 #endif /* TENSION_BALANCE_H */
2008 <<tension balance functions>>=
2009 <<one dimensional bracket>>
2010 <<one dimensional solve>>
2012 <<tension balance function>>
2015 <<tension balance module makefile lines>>=
2016 NW_SAWSIM_MODS += tension_balance
2017 CHECK_BINS += check_tension_balance
2018 check_tension_balance_MODS =
2021 The entire force balancing problem reduces to a solving two nested
2022 one-dimensional root-finding problems. First we define the one
2023 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
2024 non-empty group). $x(x_0)$ is also strictly monotonically increasing,
2025 so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
2026 stores the last successful value of $x$, since we expect to be taking
2027 small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
2028 indicates that the groups have changed.
2029 <<tension balance function declaration>>=
2030 double tension_balance(int num_tension_groups,
2031 one_dim_fn_t **tension_handlers,
2036 <<tension balance function>>=
2037 double tension_balance(int num_tension_groups,
2038 one_dim_fn_t **tension_handlers,
2043 static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL};
2044 double F, xo_guess, xo, lb, ub;
2045 double min_relx=1e-6, min_rely=1e-6;
2046 int max_steps=200, i;
2048 assert(num_tension_groups > 0);
2049 assert(tension_handlers != NULL);
2050 assert(params != NULL);
2051 assert(num_tension_groups > 0);
2053 if (num_tension_groups == 1) { /* only one group, no balancing required */
2056 //fprintf(stderr, "balancing for x = %g with ", x);
2057 //for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
2058 //fprintf(stderr, "\n");
2059 if (last_x == -1) { /* new group setup, reset x_xo_data */
2060 /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
2061 if (x_xo_data.xi != NULL)
2063 /* construct new x_xo_data */
2064 x_xo_data.n_groups = num_tension_groups;
2065 x_xo_data.tension_handler = tension_handlers;
2066 x_xo_data.handler_data = params;
2067 x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
2068 for (i=0; i<num_tension_groups; i++)
2069 x_xo_data.xi[i] = -1.0;
2071 if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
2072 if (x == 0) xo_guess = length_scale;
2073 else xo_guess = x/num_tension_groups;
2075 fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
2077 oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
2078 } else { /* work off the last balanced point */
2080 fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
2083 lb = x_xo_data.xi[0];
2084 ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */
2085 } else if (x < last_x) {
2086 lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
2087 ub = x_xo_data.xi[0];
2088 } else { /* x == last_x */
2089 fprintf(stderr, "not moving\n");
2090 lb= x_xo_data.xi[0];
2091 ub= x_xo_data.xi[0];
2095 fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
2096 __FUNCTION__, x, lb, ub);
2098 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
2099 F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
2101 if (num_tension_groups > 1) {
2102 fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
2103 for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
2104 fprintf(stderr, "\n");
2108 F = (*tension_handlers[0])(xo, params[0]);
2113 <<tension balance internal definitions>>=
2114 <<x of x0 definitions>>
2117 <<x of x0 definitions>>=
2118 typedef struct x_of_xo_data_struct {
2120 one_dim_fn_t **tension_handler; /* array of fn pointers */
2121 void **handler_data; /* array of void* pointers */
2127 double x_of_xo(double xo, void *pdata)
2129 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
2130 double F, x=0, xi, xi_guess, lb, ub;
2132 double min_relx=1e-6, min_rely=1e-6;
2134 assert(data->n_groups > 0);
2136 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
2138 fprintf(stderr, "%s: finding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
2141 if (data->xi) data->xi[0] = xo;
2142 for (i=1; i < data->n_groups; i++) {
2143 if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
2144 else xi_guess = data->xi[i];
2146 fprintf(stderr, "%s: searching for proper x[%d] for F[%d](x[%d]) = %g N\n", __FUNCTION__, i, i, i, F);
2148 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
2149 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
2151 fprintf(stderr, "%s: found F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
2155 if (data->xi) data->xi[i] = xi;
2158 fprintf(stderr, "%s: found x(x0=%g) = %g\n", __FUNCTION__, xo, x);
2164 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
2165 number of steps for monotonically increasing $f(x)$. Simple
2166 bisection, so it's robust and converges fairly quickly. You can set
2167 the maximum number of steps to take, but if you have enough processor
2168 power, it's better to limit your precision with the relative limits
2169 [[min_relx]] and [[y]]. These ensure that the width of the bracket is
2170 small on the length scale of it's larger side. Note that \emph{both}
2171 relative limits must be satisfied for the search to stop.
2172 <<one dimensional function definition>>=
2173 /* equivalent to gsl_func_t */
2174 typedef double one_dim_fn_t(double x, void *params);
2176 <<one dimensional solve declaration>>=
2177 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2178 double min_relx, double min_rely, int max_steps,
2182 <<one dimensional solve>>=
2183 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2184 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2185 double min_relx, double min_rely, int max_steps,
2188 double yx, ylb, yub, x;
2191 ylb = (*f)(lb, params);
2192 yub = (*f)(ub, params);
2194 /* check some simple cases */
2196 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bounded */
2197 else return lb; /* any x on the range [lb, ub] would work */
2199 if (ylb == y) { x=lb; yx=ylb; goto end; }
2200 if (yub == y) { x=ub; yx=yub; goto end; }
2203 fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2205 assert(ylb < y); assert(yub > y);
2206 for (i=0; i<max_steps; i++) {
2208 yx = (*f)(x, params);
2210 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);
2212 if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2213 else if (yx > y) { ub=x; yub=yx; }
2214 else /* yx < y */ { lb=x; ylb=yx; }
2219 fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2225 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2226 Generate bracketing $x$ values through bisection or exponential growth.
2227 <<one dimensional bracket declaration>>=
2228 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2231 <<one dimensional bracket>>=
2232 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2234 double yx, step, x=xguess;
2236 yx = (*f)(x, params);
2238 fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2245 fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2249 if (x == 0) x = length_scale/2.0; /* HACK */
2252 yx = (*f)(x, params);
2254 fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2259 fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2262 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2266 \subsection{Balancing implementation}
2268 <<tension-balance.c>>=
2271 <<tension balance includes>>
2272 #include "tension_balance.h"
2274 double length_scale = 1e-10; /* HACK */
2276 <<tension balance internal definitions>>
2277 <<tension balance functions>>
2280 <<tension balance includes>>=
2281 #include <assert.h> /* assert() */
2282 #include <stdlib.h> /* NULL */
2283 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2284 #include <stdio.h> /* fprintf(), stdout */
2288 \subsection{Balancing unit tests}
2290 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2291 <<check-tension-balance.c>>=
2292 <<tension balance check includes>>
2293 <<tension balance test suite>>
2294 <<main check program>>
2297 <<tension balance check includes>>=
2298 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2299 #include <assert.h> /* assert() */
2302 #include "tension_balance.h"
2305 <<tension balance test suite>>=
2306 <<tension balance function tests>>
2307 <<tension balance suite definition>>
2310 <<tension balance suite definition>>=
2311 Suite *test_suite (void)
2313 Suite *s = suite_create ("tension balance");
2314 <<tension balance function test case defs>>
2316 <<tension balance function test case adds>>
2321 <<tension balance function tests>>=
2322 <<check relative error>>
2324 double hooke(double x, void *pK)
2327 return *((double*)pK) * x;
2330 START_TEST(test_single_function)
2332 double k=5, x=3, last_x=2.0, F;
2333 one_dim_fn_t *handlers[] = {&hooke};
2334 void *data[] = {&k};
2335 F = tension_balance(1, handlers, data, last_x, x);
2336 fail_unless(F = k*x, NULL);
2341 We can also test balancing two springs with different spring constants.
2342 The analytic solution for a general number of springs is given in Appendix \ref{app.math_hooke}.
2343 <<tension balance function tests>>=
2344 START_TEST(test_double_hooke)
2346 double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
2347 one_dim_fn_t *handlers[] = {&hooke, &hooke};
2348 void *data[] = {&k1, &k2};
2349 F = tension_balance(2, handlers, data, last_x, x);
2353 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2354 //CHECK_ERR(1e-6, x1e, xi[0]);
2355 //CHECK_ERR(1e-6, x2e, xi[1]);
2356 CHECK_ERR(1e-6, Fe, F);
2360 TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
2362 <<tension balance function test case defs>>=
2363 TCase *tc_tbfunc = tcase_create("tension balance function");
2366 <<tension balance function test case adds>>=
2367 tcase_add_test(tc_tbfunc, test_single_function);
2368 tcase_add_test(tc_tbfunc, test_double_hooke);
2369 suite_add_tcase(s, tc_tbfunc);
2374 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2375 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2376 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2382 <<list definitions>>
2383 <<list declarations>>
2387 <<list declarations>>=
2388 <<head and tail declarations>>
2389 <<list length declaration>>
2390 <<push declaration>>
2392 <<list destroy declaration>>
2393 <<list to array declaration>>
2394 <<move declaration>>
2395 <<sort declaration>>
2396 <<unique declaration>>
2397 <<list print declaration>>
2401 <<create and destroy node>>
2414 <<list module makefile lines>>=
2415 NW_SAWSIM_MODS += list
2416 CHECK_BINS += check_list
2420 <<list definitions>>=
2421 typedef struct list_struct {
2422 struct list_struct *next;
2423 struct list_struct *prev;
2427 <<comparison function typedef>>
2430 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2431 <<head and tail declarations>>=
2432 list_t *head(list_t *list);
2433 list_t *tail(list_t *list);
2436 list_t *head(list_t *list)
2440 while (list->prev != NULL) {
2446 list_t *tail(list_t *list)
2450 while (list->next != NULL) {
2457 <<list length declaration>>=
2458 int list_length(list_t *list);
2461 int list_length(list_t *list)
2468 while (list->next != NULL) {
2476 [[push]] inserts a node relative to the current position in the list
2477 without changing the current position.
2478 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2479 so we set it to that of the pushed domain.
2480 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2481 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2482 <<push declaration>>=
2483 void push(list_t **pList, void *data, int below);
2486 void push(list_t **pList, void *data, int below)
2488 list_t *list, *new_node;
2489 assert(pList != NULL);
2491 new_node = create_node(data);
2494 else if (below==0) { /* insert above */
2495 if (list->prev != NULL)
2496 list->prev->next = new_node;
2497 new_node->prev = list->prev;
2498 list->prev = new_node;
2499 new_node->next = list;
2500 } else { /* insert below */
2501 if (list->next != NULL)
2502 list->next->prev = new_node;
2503 new_node->next = list->next;
2504 list->next = new_node;
2505 new_node->prev = list;
2510 [[pop]] removes the current domain node, moving the current position
2511 to the node after it, unless that node is [[NULL]], in which case move
2512 the current position to the node before it.
2513 <<pop declaration>>=
2514 void *pop(list_t **pList);
2517 void *pop(list_t **pList)
2519 list_t *list, *popped;
2521 assert(pList != NULL);
2523 assert(list != NULL); /* not an empty list */
2525 /* bypass the popped node */
2526 if (list->prev != NULL)
2527 list->prev->next = list->next;
2528 if (list->next != NULL) {
2529 list->next->prev = list->prev;
2530 *pList = list->next;
2532 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
2534 destroy_node(popped);
2539 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
2540 If the cleanup function is [[NULL]], no cleanup function is called.
2541 The nodes are not popped in any particular order.
2542 <<list destroy declaration>>=
2543 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
2546 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
2550 assert(pList != NULL);
2553 assert(list != NULL); /* not an empty list */
2554 while (list != NULL) {
2556 if (destroy != NULL)
2562 [[list_to_array]] converts a list to an array of pointers, preserving order.
2563 <<list to array declaration>>=
2564 void list_to_array(list_t *list, int *length, void ***array);
2567 void list_to_array(list_t *list, int *pLength, void ***pArray)
2571 assert(list != NULL);
2572 assert(pLength != NULL);
2573 assert(pArray != NULL);
2575 length = list_length(list);
2576 array = (void **)malloc(sizeof(void *)*length);
2577 assert(array != NULL);
2580 while (list != NULL) {
2581 array[i++] = list->d;
2589 [[move]] moves the current element along the list in either direction.
2590 <<move declaration>>=
2591 void move(list_t **pList, int downstream);
2594 void move(list_t **pList, int downstream)
2596 list_t *list, *flipper;
2598 assert(pList != NULL);
2599 list = *pList; /* a>B>c>d */
2600 assert(list != NULL); /* not an empty list */
2601 if (downstream == 0)
2602 flipper = list->prev; /* flipper is A */
2604 flipper = list->next; /* flipper is C */
2605 /* ensure there is somewhere to go in the direction we're heading */
2606 assert(flipper != NULL);
2607 /* remove the current element */
2608 data = pop(&list); /* data=B, a>C>d */
2609 /* and put it back in in the correct spot */
2611 if (downstream == 0) {
2612 push(&list, data, 0); /* b>A>c>d */
2613 list = list->prev; /* B>a>c>d */
2615 push(&list, data, 1); /* a>C>b>d */
2616 list = list->next; /* a>c>B>d */
2622 [[sort]] sorts a list from smallest to largest according to a user
2623 specified comparison.
2624 <<comparison function typedef>>=
2625 typedef int (comparison_fn_t)(void *A, void *B);
2628 <<int comparison function>>=
2629 int int_comparison(void *A, void *B)
2634 if (a > b) return 1;
2635 else if (a == b) return 0;
2639 <<sort declaration>>=
2640 void sort(list_t **pList, comparison_fn_t *comp);
2642 Here we use the bubble sort algorithm.
2644 void sort(list_t **pList, comparison_fn_t *comp)
2648 assert(pList != NULL);
2650 assert(list != NULL); /* not an empty list */
2651 while (stable == 0) {
2654 while (list->next != NULL) {
2655 if ((*comp)(list->d, list->next->d) > 0) {
2657 move(&list, 1 /* downstream */);
2666 [[unique]] removes duplicates from a sorted list according to a user
2667 specified comparison. Don't do this unless you have other pointers
2668 any dynamically allocated data in your list, because [[unique]] just
2669 drops any non-unique data without freeing it.
2670 <<unique declaration>>=
2671 void unique(list_t **pList, comparison_fn_t *comp);
2674 void unique(list_t **pList, comparison_fn_t *comp)
2677 assert(pList != NULL);
2679 assert(list != NULL); /* not an empty list */
2681 while (list->next != NULL) {
2682 if ((*comp)(list->d, list->next->d) == 0) {
2691 [[list_print]] prints a list to a [[FILE *]] stream.
2692 <<list print declaration>>=
2693 void list_print(FILE *file, list_t *list, char *list_name);
2696 void list_print(FILE *file, list_t *list, char *list_name)
2698 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
2700 while (list != NULL) {
2701 fprintf(file, " %p", list->d);
2704 fprintf(file, "\n");
2708 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
2709 <<create and destroy node>>=
2710 list_t *create_node(void *data)
2712 list_t *ret = (list_t *)malloc(sizeof(list_t));
2713 assert(ret != NULL);
2720 void destroy_node(list_t *node)
2726 The user must free the data pointed to by the node on their own.
2728 \subsection{List implementation}
2738 #include <assert.h> /* assert() */
2739 #include <stdlib.h> /* malloc(), free(), rand() */
2740 #include <stdio.h> /* fprintf(), stdout, FILE */
2741 #include "global.h" /* destroy_data_func_t */
2744 \subsection{List unit tests}
2746 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2748 <<list check includes>>
2750 <<main check program>>
2753 <<list check includes>>=
2754 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2755 #include <stdio.h> /* FILE */
2761 <<list test suite>>=
2764 <<list suite definition>>
2767 <<list suite definition>>=
2768 Suite *test_suite (void)
2770 Suite *s = suite_create ("list");
2771 <<push test case defs>>
2772 <<pop test case defs>>
2774 <<push test case adds>>
2775 <<pop test case adds>>
2781 START_TEST(test_push)
2786 push(&list, (void *)i, 1);
2787 fail_unless(list == head(list), NULL);
2788 fail_unless( (int)list->d == 0 );
2789 for (i=0; i<n; i++) {
2790 /* we expect 0, n-1, n-2, ..., 1 */
2793 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
2799 <<push test case defs>>=
2800 TCase *tc_push = tcase_create("push");
2803 <<push test case adds>>=
2804 tcase_add_test(tc_push, test_push);
2805 suite_add_tcase(s, tc_push);
2810 <<pop test case defs>>=
2812 <<pop test case adds>>=
2815 \section{Function string evaluation}
2817 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).
2818 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
2819 The solution is to run a scripting language as a subprocess accessed via pipes.
2820 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
2822 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
2823 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
2824 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.
2825 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
2827 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.
2828 Then you can either statically or dynamically link to those hardcoded functions.
2829 While much less flexible, this approach would be a fairly simple-to-implement fallback.
2831 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
2832 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
2835 #ifndef STRING_EVAL_H
2836 #define STRING_EVAL_H
2838 <<string eval setup declaration>>
2839 <<string eval function declaration>>
2840 <<string eval teardown declaration>>
2841 #endif /* STRING_EVAL_H */
2844 <<string eval module makefile lines>>=
2845 NW_SAWSIM_MODS += string_eval
2846 CHECK_BINS += check_string_eval
2847 check_string_eval_MODS =
2850 For an introduction to POSIX process control, see\\
2851 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
2852 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
2853 and of course, the relavent [[man]] pages.
2855 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
2856 [[execvp]] replaces the calling process' program with a new program.
2857 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
2858 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
2859 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
2860 The new program needs command line arguments, just like it would if you were running it from a shell.
2861 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
2862 with the final array entry being a [[NULL]] pointer.
2864 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
2865 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
2866 <<string eval subprocess definitions>>=
2867 #define SUBPROCESS "python"
2868 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
2869 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
2870 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
2872 The [[i]] option lets Python know that it should run in interactive mode.
2873 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
2874 In interactive mode, python acts on every instruction as soon as it is recieved.
2875 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.
2876 %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.
2878 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
2879 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
2880 [[fork]] returns two copies of the same program, executing the original code.
2881 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
2882 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
2884 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.
2885 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
2886 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
2887 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
2888 <<string eval pipe definitions>>=
2889 #define PIPE_READ 0 /* the end you read from */
2890 #define PIPE_WRITE 1 /* the end you write to */
2892 #define STDIN 0 /* initial index of stdin pair */
2893 #define STDOUT 2 /* initial index of stdout pair */
2896 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
2898 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.
2900 <<string eval setup declaration>>=
2901 extern void string_eval_setup(FILE **pIn, FILE **pOut);
2903 <<string eval setup definition>>=
2904 void string_eval_setup(FILE **pIn, FILE **pOut)
2907 int pfd[4]; /* file descriptors for stdin and stdout */
2909 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
2910 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
2912 pid = fork(); /* split process into two copies */
2914 if (pid == -1) { /* fork error */
2915 perror("fork error");
2917 } else if (pid == 0) { /* child */
2918 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
2919 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
2920 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
2921 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
2922 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
2923 perror("exec error"); /* exec shouldn't return */
2925 } else { /* parent */
2926 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
2927 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
2928 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
2929 if ( *pIn == NULL ) {
2930 perror("fdopen (in)");
2933 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
2934 if ( *pOut == NULL ) {
2935 perror("fdopen (out)");
2942 To use the evaluating subprocess, we just pipe in our command, and read out the results.
2943 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
2944 <<string eval function declaration>>=
2945 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
2947 <<string eval function definition>>=
2948 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
2951 rc = fprintf(in, "%s", input);
2952 assert(rc == strlen(input));
2955 alarm(1); /* set a one second timeout on the read */
2956 assert( fgets(output, buflen, out) != NULL );
2957 alarm(0); /* cancel the timeout */
2958 //fprintf(stderr, "eval: %s ----> %s", input, output);
2961 The [[alarm]] calls set and clear a timeout on the returned output.
2962 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.
2963 This protects against invalid input for which a line of output is not printed to [[stdout]].
2964 Other invalid inputs (e.g. those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
2965 If you are getting strange results, check your python code seperately. TODO, better error handling.
2967 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
2968 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
2969 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
2970 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
2971 <<string eval teardown declaration>>=
2972 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
2975 <<string eval teardown definition>>=
2976 void string_eval_teardown(FILE **pIn, FILE **pOut)
2981 /* redirect Python's stderr */
2982 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
2986 assert( fclose(*pIn) == 0 );
2988 assert( fclose(*pOut) == 0 );
2991 /* wait for python to exit */
2993 pid = wait(&stat_loc);
3000 if (WIFEXITED(stat_loc)) {
3001 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
3002 } else if (WIFSIGNALED(stat_loc)) {
3003 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
3008 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
3010 \subsection{String evaluation implementation}
3014 <<string eval includes>>
3015 #include "string_eval.h"
3016 <<string eval internal definitions>>
3017 <<string eval functions>>
3020 <<string eval includes>>=
3021 #include <assert.h> /* assert() */
3022 #include <stdlib.h> /* NULL */
3023 #include <stdio.h> /* fprintf(), stdout, fdopen() */
3024 #include <string.h> /* strlen() */
3025 #include <sys/types.h> /* pid_t */
3026 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
3027 #include <wait.h> /* wait() */
3030 <<string eval internal definitions>>=
3031 <<string eval subprocess definitions>>
3032 <<string eval pipe definitions>>
3035 <<string eval functions>>=
3036 <<string eval setup definition>>
3037 <<string eval function definition>>
3038 <<string eval teardown definition>>
3041 \subsection{String evaluation unit tests}
3043 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3044 <<check-string-eval.c>>=
3045 <<string eval check includes>>
3046 <<string comparison definition and globals>>
3047 <<string eval test suite>>
3048 <<main check program>>
3051 <<string eval check includes>>=
3052 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
3053 #include <stdio.h> /* printf() */
3054 #include <assert.h> /* assert() */
3055 #include <string.h> /* strlen() */
3056 #include <signal.h> /* SIGKILL */
3058 #include "string_eval.h"
3061 <<string eval test suite>>=
3062 <<string eval tests>>
3063 <<string eval suite definition>>
3066 <<string eval suite definition>>=
3067 Suite *test_suite (void)
3069 Suite *s = suite_create ("string eval");
3070 <<string eval test case defs>>
3072 <<string eval test case add>>
3077 <<string eval tests>>=
3078 START_TEST(test_setup_teardown)
3081 string_eval_setup(&in, &out);
3082 string_eval_teardown(&in, &out);
3085 START_TEST(test_invalid_command)
3088 char input[80], output[80]={};
3089 string_eval_setup(&in, &out);
3090 sprintf(input, "print ABCDefg\n");
3091 string_eval(in, out, input, 80, output);
3092 string_eval_teardown(&in, &out);
3095 START_TEST(test_simple_eval)
3098 char input[80], output[80]={};
3099 string_eval_setup(&in, &out);
3100 sprintf(input, "print 3+4*5\n");
3101 string_eval(in, out, input, 80, output);
3102 fail_unless(STRMATCH(output,"23\n"), NULL);
3103 string_eval_teardown(&in, &out);
3106 START_TEST(test_multiple_evals)
3109 char input[80], output[80]={};
3110 string_eval_setup(&in, &out);
3111 sprintf(input, "print 3+4*5\n");
3112 string_eval(in, out, input, 80, output);
3113 fail_unless(STRMATCH(output,"23\n"), NULL);
3114 sprintf(input, "print (3**2 + 4**2)**0.5\n");
3115 string_eval(in, out, input, 80, output);
3116 fail_unless(STRMATCH(output,"5.0\n"), NULL);
3117 string_eval_teardown(&in, &out);
3120 START_TEST(test_eval_with_state)
3123 char input[80], output[80]={};
3124 string_eval_setup(&in, &out);
3125 sprintf(input, "print 3+4*5\n");
3126 fprintf(in, "A = 3\n");
3127 sprintf(input, "print A*3\n");
3128 string_eval(in, out, input, 80, output);
3129 fail_unless(STRMATCH(output,"9\n"), NULL);
3130 string_eval_teardown(&in, &out);
3134 <<string eval test case defs>>=
3135 TCase *tc_string_eval = tcase_create("string_eval");
3137 <<string eval test case add>>=
3138 tcase_add_test(tc_string_eval, test_setup_teardown);
3139 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
3140 tcase_add_test(tc_string_eval, test_simple_eval);
3141 tcase_add_test(tc_string_eval, test_multiple_evals);
3142 tcase_add_test(tc_string_eval, test_eval_with_state);
3143 suite_add_tcase(s, tc_string_eval);
3147 \section{Accelerating function evaluation}
3149 My first version-0.3 code was running very slowly.
3150 With the options suggested in the help
3151 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
3152 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
3153 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
3154 That's only 15 calls per solution, so the search algorithm seems reasonable.
3155 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
3160 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
3162 #endif /* ACCEL_K_H */
3165 <<accel k module makefile lines>>=
3166 NW_SAWSIM_MODS += accel_k
3167 #CHECK_BINS += check_accel_k
3168 check_accel_k_MODS =
3172 #include <assert.h> /* assert() */
3173 #include <stdlib.h> /* realloc(), free(), NULL */
3174 #include "global.h" /* environment_t */
3175 #include "k_model.h" /* k_func_t */
3176 #include "interp.h" /* interp_* */
3177 #include "accel_k.h"
3179 typedef struct accel_k_struct {
3180 interp_table_t *itable;
3186 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3187 static int num_accels = 0;
3188 static accel_k_t *accels=NULL;
3190 /* Wrap k in the standard f(x) acceleration form */
3191 static double k_wrap(double F, void *params)
3193 accel_k_t *p = (accel_k_t *)params;
3194 return (*p->k)(F, p->env, p->params);
3197 static int k_tol(double FA, double kA, double FB, double kB)
3200 if (FB-FA > 1e-12) {
3201 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3202 return 1; /* unnacceptable */
3204 //printf("acceptable tol\n");
3205 return 0; /* acceptable */
3209 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3212 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3213 assert(accels != NULL);
3214 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3216 accels[i].env = env;
3217 accels[i].params = params;
3224 for (i=0; i<num_accels; i++)
3225 interp_table_free(accels[i].itable);
3227 if (accels) free(accels);
3231 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3234 for (i=0; i<num_accels; i++) {
3235 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3236 /* We've already setup this function.
3237 * WARNING: we're assuming param and env are the same. */
3238 return interp_table_eval(accels[i].itable, accels+i, F);
3241 /* set up a new acceleration environment */
3242 i = add_accel_k(k, env, params);
3243 return interp_table_eval(accels[i].itable, accels+i, F);
3247 \section{Tension models}
3248 \label{sec.tension_models}
3250 TODO: tension model intro.
3251 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.
3253 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3254 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]].
3256 <<tension-model.h>>=
3257 #ifndef TENSION_MODEL_H
3258 #define TENSION_MODEL_H
3260 <<tension handler types>>
3261 <<constant tension model declarations>>
3262 <<hooke tension model declarations>>
3263 <<worm-like chain tension model declarations>>
3264 <<freely-jointed chain tension model declarations>>
3265 <<find tension definitions>>
3266 <<tension model global declarations>>
3267 #endif /* TENSION_MODEL_H */
3270 <<tension model module makefile lines>>=
3271 NW_SAWSIM_MODS += tension_model
3272 #CHECK_BINS += check_tension_model
3273 check_tension_model_MODS =
3275 <<tension model utils makefile lines>>=
3276 TENSION_MODEL_MODS = tension_model parse list tension_balance
3277 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3278 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3279 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3280 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
3281 notangle -Rtension-model-utils.c $< > $@
3282 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
3283 gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3284 clean_tension_model_utils :
3285 rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3286 @ The pipe symbol [[|]] marks the prerequisits following it (in this
3287 case, the directories) as ‘order-only’ prerequisites. The timestamp
3288 on these prerequisits does not effect whether the rules are executed.
3289 This is appropriate for directories, where we don't need to recompile
3290 after an unrelated has been added to the directory, but only when the
3291 source prerequisites change. See the [[make]] documentation for more
3293 (\url{http://www.gnu.org/software/make/manual/html_node/Prerequisite-Types.html}).
3296 \label{sec.null_tension}
3298 For unstretchable domains.
3300 <<null tension model getopt>>=
3301 {NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3304 \subsection{Constant}
3305 \label{sec.const_tension}
3307 <<constant tension functions>>=
3308 <<constant tension function>>
3309 <<constant tension parameter create/destroy functions>>
3312 <<constant tension model declarations>>=
3313 <<constant tension function declaration>>
3314 <<constant tension parameter create/destroy function declarations>>
3315 <<constant tension model global declarations>>
3319 An infinitely stretchable domain providing a constant tension.
3321 <<constant tension function declaration>>=
3322 extern double const_tension_handler(double x, void *pdata);
3324 <<constant tension function>>=
3325 double const_tension_handler(double x, void *pdata)
3327 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3328 list_t *list = data->group;
3333 assert(list != NULL); /* empty active group?! */
3334 F = ((const_tension_param_t *)list->d)->F;
3336 fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3338 while (list != NULL) {
3339 assert(((const_tension_param_t *)list->d)->F == F);
3347 <<constant tension structure>>=
3348 typedef struct const_tension_param_struct {
3349 double F; /* tension (force) in N */
3350 } const_tension_param_t;
3354 <<constant tension parameter create/destroy function declarations>>=
3355 extern void *string_create_const_tension_param_t(char **param_strings);
3356 extern void destroy_const_tension_param_t(void *p);
3359 <<constant tension parameter create/destroy functions>>=
3360 const_tension_param_t *create_const_tension_param_t(double F)
3362 const_tension_param_t *ret
3363 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3364 assert(ret != NULL);
3369 void *string_create_const_tension_param_t(char **param_strings)
3371 return create_const_tension_param_t(atof(param_strings[0]));
3374 void destroy_const_tension_param_t(void *p)
3382 <<constant tension model global declarations>>=
3383 extern int num_const_tension_params;
3384 extern char *const_tension_param_descriptions[];
3385 extern char const_tension_param_string[];
3388 <<constant tension model globals>>=
3389 int num_const_tension_params = 1;
3390 char *const_tension_param_descriptions[] = {"tension F, N"};
3391 char const_tension_param_string[] = "0";
3394 <<constant tension model getopt>>=
3395 {&const_tension_handler, "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}
3401 <<hooke functions>>=
3403 <<hooke parameter create/destroy functions>>
3406 <<hooke tension model declarations>>=
3407 <<hooke tension function declaration>>
3408 <<hooke tension parameter create/destroy function declarations>>
3409 <<hooke tension model global declarations>>
3413 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
3414 The behavior of a series of springs $k_i$ in series is given by
3416 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
3418 For a simple proof, see Appendix \ref{app.math_hooke}.
3420 <<hooke tension function declaration>>=
3421 extern double hooke_handler(double x, void *pdata);
3425 double hooke_handler(double x, void *pdata)
3427 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3428 list_t *list = data->group;
3433 assert(list != NULL); /* empty active group?! */
3434 while (list != NULL) {
3435 assert( ((hooke_param_t *)list->d)->k > 0 );
3436 k += 1.0/ ((hooke_param_t *)list->d)->k;
3438 fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
3439 __FUNCTION__, ((hooke_param_t *)list->d)->k);
3445 fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N\n",
3446 __FUNCTION__, k, x, k*x);
3452 <<hooke structure>>=
3453 typedef struct hooke_param_struct {
3454 double k; /* spring constant in N/m */
3458 <<hooke tension parameter create/destroy function declarations>>=
3459 extern void *string_create_hooke_param_t(char **param_strings);
3460 extern void destroy_hooke_param_t(void *p);
3463 <<hooke parameter create/destroy functions>>=
3464 hooke_param_t *create_hooke_param_t(double k)
3466 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
3467 assert(ret != NULL);
3472 void *string_create_hooke_param_t(char **param_strings)
3474 return create_hooke_param_t(atof(param_strings[0]));
3477 void destroy_hooke_param_t(void *p)
3484 <<hooke tension model global declarations>>=
3485 extern int num_hooke_params;
3486 extern char *hooke_param_descriptions[];
3487 extern char hooke_param_string[];
3490 <<hooke tension model globals>>=
3491 int num_hooke_params = 1;
3492 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
3493 char hooke_param_string[]="0.05";
3496 <<hooke tension model getopt>>=
3497 {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}
3500 \subsection{Worm-like chain}
3503 We can model several unfolded domains with a single worm-like chain.
3504 <<worm-like chain functions>>=
3505 <<worm-like chain function>>
3506 <<worm-like chain handler>>
3507 <<worm-like chain create/destroy functions>>
3510 <<worm-like chain tension model declarations>>=
3511 <<worm-like chain handler declaration>>
3512 <<worm-like chain create/destroy function declarations>>
3513 <<worm-like chain tension model global declarations>>
3517 The combination of all unfolded domains is modeled as a worm like chain,
3518 with the force $F_{\text{WLC}}$ approximately given by
3520 F_{\text{WLC}} \approx \frac{k_B T}{p}
3521 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
3523 where $x$ is the distance between the fixed ends,
3524 $k_B$ is Boltzmann's constant,
3525 $T$ is the absolute temperature,
3526 $p$ is the persistence length, and
3527 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
3528 <<worm-like chain function>>=
3529 double wlc(double x, double T, double p, double L)
3531 double xL=0.0; /* uses global k_B */
3533 assert(T > 0); assert(p > 0); assert(L > 0);
3534 if (x >= L) return HUGE_VAL;
3535 xL = x/L; /* unitless */
3536 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
3537 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
3538 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
3541 This model assumes that each unfolded domain has the same persistence length.
3543 <<worm-like chain handler declaration>>=
3544 extern double wlc_handler(double x, void *pdata);
3547 <<worm-like chain handler>>=
3548 double wlc_handler(double x, void *pdata)
3550 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3551 list_t *list = data->group;
3555 assert(list != NULL); /* empty active group?! */
3556 p = ((wlc_param_t *)list->d)->p;
3557 while (list != NULL) {
3558 /* enforce consistent p */
3559 assert( ((wlc_param_t *)list->d)->p == p);
3560 L += ((wlc_param_t *)list->d)->L;
3562 fprintf(stderr, "%s: WLC section %g m long in series\n",
3563 __FUNCTION__, ((wlc_param_t *)list->d)->L);
3568 fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
3569 __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
3571 return wlc(x, data->env->T, p, L);
3575 <<worm-like chain structure>>=
3576 typedef struct wlc_param_struct {
3577 double p; /* WLC persistence length (m) */
3578 double L; /* WLC contour length (m) */
3582 <<worm-like chain create/destroy function declarations>>=
3583 extern void *string_create_wlc_param_t(char **param_strings);
3584 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
3587 <<worm-like chain create/destroy functions>>=
3588 wlc_param_t *create_wlc_param_t(double p, double L)
3590 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
3591 assert(ret != NULL);
3597 void *string_create_wlc_param_t(char **param_strings)
3599 return create_wlc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3602 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
3609 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.
3610 TODO (now we copy the string before parsing, but still do this for clarity).
3612 <<valid string write code>>=
3613 char string[] = "abc";
3616 <<valid string write code>>=
3617 char *string = "abc";
3621 <<worm-like chain tension model global declarations>>=
3622 extern int num_wlc_params;
3623 extern char *wlc_param_descriptions[];
3624 extern char wlc_param_string[];
3626 <<worm-like chain tension model globals>>=
3627 int num_wlc_params = 2;
3628 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
3629 char wlc_param_string[]="0.39e-9,28.4e-9";
3632 <<worm-like chain tension model getopt>>=
3633 {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}
3635 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
3637 \subsection{Freely-jointed chain}
3640 <<freely-jointed chain functions>>=
3641 <<freely-jointed chain function>>
3642 <<freely-jointed chain handler>>
3643 <<freely-jointed chain create/destroy functions>>
3646 <<freely-jointed chain tension model declarations>>=
3647 <<freely-jointed chain handler declaration>>
3648 <<freely-jointed chain create/destroy function declarations>>
3649 <<freely-jointed chain tension model global declarations>>
3653 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
3654 The tension of a chain of $N$ such links is given by
3656 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
3658 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}.
3659 <<freely-jointed chain function>>=
3660 double langevin(double x, void *params)
3663 return 1.0/tanh(x) - 1/x;
3665 <<one dimensional bracket declaration>>
3666 <<one dimensional solve declaration>>
3667 double inv_langevin(double x)
3669 double lb=0.0, ub=1.0, ret;
3670 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
3671 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
3675 double fjc(double x, double T, double l, int N)
3677 double L = l*(double)N;
3678 if (x == 0) return 0;
3679 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
3680 return k_B*T/l * inv_langevin(x/L);
3684 <<freely-jointed chain handler declaration>>=
3685 extern double fjc_handler(double x, void *pdata);
3687 <<freely-jointed chain handler>>=
3688 double fjc_handler(double x, void *pdata)
3690 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3691 list_t *list = data->group;
3696 assert(list != NULL); /* empty active group?! */
3697 l = ((fjc_param_t *)list->d)->l;
3698 while (list != NULL) {
3699 /* enforce consistent l */
3700 assert( ((fjc_param_t *)list->d)->l == l);
3701 N += ((fjc_param_t *)list->d)->N;
3703 fprintf(stderr, "%s: FJC section %d links long in series\n",
3704 __FUNCTION__, ((fjc_param_t *)list->d)->N);
3709 fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
3710 __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
3712 return fjc(x, data->env->T, l, N);
3716 <<freely-jointed chain structure>>=
3717 typedef struct fjc_param_struct {
3718 double l; /* FJC link length (m) */
3719 int N; /* FJC number of links */
3723 <<freely-jointed chain create/destroy function declarations>>=
3724 extern void *string_create_fjc_param_t(char **param_strings);
3725 extern void destroy_fjc_param_t(void *p);
3728 <<freely-jointed chain create/destroy functions>>=
3729 fjc_param_t *create_fjc_param_t(double l, double N)
3731 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
3732 assert(ret != NULL);
3738 void *string_create_fjc_param_t(char **param_strings)
3740 return create_fjc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3743 void destroy_fjc_param_t(void *p)
3750 <<freely-jointed chain tension model global declarations>>=
3751 extern int num_fjc_params;
3752 extern char *fjc_param_descriptions[];
3753 extern char fjc_param_string[];
3756 <<freely-jointed chain tension model globals>>=
3757 int num_fjc_params = 2;
3758 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
3759 char fjc_param_string[]="0.5e-9,200";
3762 <<freely-jointed chain tension model getopt>>=
3763 {fjc_handler, "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}
3766 \subsection{Tension model implementation}
3768 <<tension-model.c>>=
3771 <<tension model includes>>
3772 #include "tension_model.h"
3773 <<tension model internal definitions>>
3774 <<tension model globals>>
3775 <<tension model functions>>
3778 <<tension model includes>>=
3779 #include <assert.h> /* assert() */
3780 #include <stdlib.h> /* NULL */
3781 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
3782 #include <stdio.h> /* fprintf(), stdout */
3785 #include "tension_balance.h" /* oneD_*() */
3788 <<tension model internal definitions>>=
3789 <<constant tension structure>>
3791 <<worm-like chain structure>>
3792 <<freely-jointed chain structure>>
3795 <<tension model globals>>=
3796 <<tension handler array global>>
3797 <<constant tension model globals>>
3798 <<hooke tension model globals>>
3799 <<worm-like chain tension model globals>>
3800 <<freely-jointed chain tension model globals>>
3803 <<tension model functions>>=
3804 <<constant tension functions>>
3806 <<worm-like chain functions>>
3807 <<freely-jointed chain functions>>
3811 \subsection{Utilities}
3813 The tension models can get complicated, and users may want to reassure
3814 themselves that this portion of the input physics are functioning properly. The
3815 stand-alone program [[tension_model_utils]] demonstrates the tension model
3816 interface without the overhead of having to understand a full simulation.
3817 [[tension_model_utils]] takes command line model arguments like the full
3818 simulation, and prints $F(x)$ data to the screen for the selected model over a
3821 <<tension-model-utils.c>>=
3823 <<tension model utility includes>>
3824 <<tension model utility definitions>>
3825 <<tension model utility getopt functions>>
3827 <<tension model utility main>>
3830 <<tension model utility main>>=
3831 int main(int argc, char **argv)
3833 <<tension model getopt array definition>>
3834 tension_model_getopt_t *model = NULL;
3838 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
3839 tension_handler_data_t tdata;
3840 int num_param_args; /* for INIT_MODEL() */
3841 char **param_args; /* for INIT_MODEL() */
3843 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model,
3844 &Fmax, &Xmax, &flags);
3846 if (flags & VFLAG) {
3847 printf("#initializing model %s with parameters %s\n", model->name, model->params);
3849 INIT_MODEL("utility", model, model->params, params);
3851 push(&tdata.group, params, 1);
3853 tdata.persist = NULL;
3854 if (model->handler == NULL) {
3855 printf("No tension function for model %s\n", model->name);
3861 printf("#Distance (m)\tForce (N)\n");
3862 for (i=0; i<=N; i++) {
3863 x = Xmax*i/(double)N;
3864 F = (*model->handler)(x, &tdata);
3865 if (F < 0 || F > Fmax) break;
3866 printf("%g\t%g\n", x, F);
3868 if (flags & VFLAG && i <= N) {
3869 /* explain exit condition */
3871 printf("Impossible force %g\n", F);
3873 printf("Reached large force limit %g > %g\n", F, Fmax);
3876 params = pop(&tdata.group);
3877 if (model->destructor)
3878 (*model->destructor)(params);
3883 <<tension model utility includes>>=
3886 #include <string.h> /* strlen() */
3887 #include <assert.h> /* assert() */
3891 #include "tension_model.h"
3894 <<tension model utility definitions>>=
3895 <<version definition>>
3896 #define VFLAG 1 /* verbose */
3897 <<string comparison definition and globals>>
3898 <<tension model getopt definitions>>
3899 <<initialize model definition>>
3903 <<tension model utility getopt functions>>=
3904 <<index tension model>>
3905 <<tension model utility help>>
3906 <<tension model utility get options>>
3909 <<tension model utility help>>=
3910 void help(char *prog_name,
3912 int n_tension_models, tension_model_getopt_t *tension_models,
3913 int tension_model, double Fmax, double Xmax)
3916 printf("usage: %s [options]\n", prog_name);
3917 printf("Version %s\n\n", VERSION);
3918 printf("A utility for understanding the available tension models\n\n");
3919 printf("Environmental options:\n");
3920 printf("-T T\tTemperature T (currently %g K)\n", env->T);
3921 printf("-C T\tYou can also set the temperature T in Celsius\n");
3922 printf("Model options:\n");
3923 printf("-m model\tFolded tension model (currently %s)\n",
3924 tension_models[tension_model].name);
3925 printf("-a args\tFolded tension model argument string (currently %s)\n",
3926 tension_models[tension_model].params);
3927 printf("F(x) is calculated for a range of x and printed\n");
3928 printf("For example:\n");
3929 printf(" #Distance (m)\tForce (N)\n");
3930 printf(" 123.456\t7.89\n");
3932 printf("-F\tSet the maximum F value for the standard mode F(x) output (currently %g)\n", Fmax);
3933 printf("-X\tSet the maximum x value for the standard mode F(x) output (currently %g)\n", Xmax);
3934 printf("-V\tChange output to verbose mode\n");
3935 printf("-h\tPrint this help and exit\n");
3937 printf("Tension models:\n");
3938 for (i=0; i<n_tension_models; i++) {
3939 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
3940 for (j=0; j < tension_models[i].num_params ; j++)
3941 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
3942 printf("\t\tdefaults: %s\n", tension_models[i].params);
3947 <<tension model utility get options>>=
3948 void get_options(int argc, char **argv, environment_t *env,
3949 int n_tension_models, tension_model_getopt_t *tension_models,
3950 tension_model_getopt_t **model, double *Fmax, double *Xmax,
3951 unsigned int *flags)
3953 char *prog_name = NULL;
3954 char c, options[] = "T:C:m:a:F:X:Vh";
3955 int tension_model=0, num_strings;
3956 extern char *optarg;
3957 extern int optind, optopt, opterr;
3961 /* setup defaults */
3963 prog_name = argv[0];
3964 env->T = 300.0; /* K */
3968 *model = tension_models;
3970 while ((c=getopt(argc, argv, options)) != -1) {
3972 case 'T': env->T = atof(optarg); break;
3973 case 'C': env->T = atof(optarg)+273.15; break;
3975 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
3976 *model = tension_models+tension_model;
3979 tension_models[tension_model].params = optarg;
3981 case 'F': *Fmax = atof(optarg); break;
3982 case 'X': *Xmax = atof(optarg); break;
3983 case 'V': *flags |= VFLAG; break;
3985 fprintf(stderr, "unrecognized option '%c'\n", optopt);
3986 /* fall through to default case */
3988 help(prog_name, env, n_tension_models, tension_models, tension_model, *Fmax, *Xmax);
3997 \section{Unfolding rate models}
3998 \label{sec.k_models}
4000 We have two main choices for determining $k$: Bell's model, which assumes the
4001 folded domains exist in a single `folded' state and make instantaneous,
4002 stochastic transitions over a free energy barrier; and Kramers' model, which
4003 treats the unfolding as a thermalized diffusion process.
4004 We deal with these options like we dealt with the different tension models: by bundling all model-specific
4005 parameters into structures, with model specific functions for determining $k$.
4007 <<k func definition>>=
4008 typedef double k_func_t(double F, environment_t *env, void *params);
4011 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
4012 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]].
4014 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
4015 because \LaTeX\ doesn't like underscores outside math mode.
4016 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
4017 You could use spaces instead (`\verb+ +'),
4018 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
4019 which I don't like as much.
4025 <<k func definition>>
4026 <<null k declarations>>
4027 <<const k declarations>>
4028 <<bell k declarations>>
4029 <<kbell k declarations>>
4030 <<kramers k declarations>>
4031 <<kramers integ k declarations>>
4032 #endif /* K_MODEL_H */
4035 <<k model module makefile lines>>=
4036 NW_SAWSIM_MODS += k_model
4037 CHECK_BINS += check_k_model
4038 check_k_model_MODS = parse string_eval
4040 [[check_k_model]] also depends on the parse module.
4042 <<k model utils makefile lines>>=
4043 K_MODEL_MODS = k_model parse string_eval
4044 K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
4045 $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h)
4046 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
4047 notangle -Rk-model-utils.c $< > $@
4048 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
4049 gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
4050 clean_k_model_utils :
4051 rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
4057 For domains that are never folded.
4059 <<null k declarations>>=
4061 <<null k functions>>=
4066 <<null k model getopt>>=
4067 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
4070 \subsection{Constant rate model}
4073 This is a pretty boring model, but it is usefull for testing the engine.
4077 where $k_0$ is the constant unfolding rate.
4079 <<const k functions>>=
4080 <<const k function>>
4081 <<const k structure create/destroy functions>>
4084 <<const k declarations>>=
4085 <<const k function declaration>>
4086 <<const k structure create/destroy function declarations>>
4087 <<const k global declarations>>
4091 <<const k structure>>=
4092 typedef struct const_k_param_struct {
4093 double knot; /* transition rate k_0 (frac population per s) */
4097 <<const k function declaration>>=
4098 double const_k(double F, environment_t *env, void *const_k_params);
4100 <<const k function>>=
4101 double const_k(double F, environment_t *env, void *const_k_params)
4102 { /* Returns the rate constant k in frac pop per s. */
4103 const_k_param_t *p = (const_k_param_t *) const_k_params;
4105 assert(p->knot > 0);
4110 <<const k structure create/destroy function declarations>>=
4111 void *string_create_const_k_param_t(char **param_strings);
4112 void destroy_const_k_param_t(void *p);
4115 <<const k structure create/destroy functions>>=
4116 const_k_param_t *create_const_k_param_t(double knot)
4118 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
4119 assert(ret != NULL);
4124 void *string_create_const_k_param_t(char **param_strings)
4126 return create_const_k_param_t(atof(param_strings[0]));
4129 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
4136 <<const k global declarations>>=
4137 extern int num_const_k_params;
4138 extern char *const_k_param_descriptions[];
4139 extern char const_k_param_string[];
4142 <<const k globals>>=
4143 int num_const_k_params = 1;
4144 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
4145 char const_k_param_string[]="1";
4148 <<const k model getopt>>=
4149 {"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}
4152 \subsection{Bell's model}
4155 Quantitatively, Bell's model gives $k$ as
4157 k = k_0 \cdot e^{F dx / k_B T} \;,
4159 where $k_0$ is the unforced unfolding rate,
4160 $F$ is the applied unfolding force,
4161 $dx$ is the distance to the transition state, and
4162 $k_B T$ is the thermal energy\citep{schlierf06}.
4164 <<bell k functions>>=
4166 <<bell k structure create/destroy functions>>
4169 <<bell k declarations>>=
4170 <<bell k function declaration>>
4171 <<bell k structure create/destroy function declarations>>
4172 <<bell k global declarations>>
4176 <<bell k structure>>=
4177 typedef struct bell_param_struct {
4178 double knot; /* transition rate k_0 (frac population per s) */
4179 double dx; /* `distance to transition state' (nm) */
4183 <<bell k function declaration>>=
4184 double bell_k(double F, environment_t *env, void *bell_params);
4186 <<bell k function>>=
4187 double bell_k(double F, environment_t *env, void *bell_params)
4188 { /* Returns the rate constant k in frac pop per s.
4189 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4190 * uses global k_B in J/K */
4191 bell_param_t *p = (bell_param_t *) bell_params;
4192 assert(F >= 0); assert(env->T > 0);
4194 assert(p->knot > 0); assert(p->dx >= 0);
4196 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
4197 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
4198 p->knot * exp(F*p->dx / (k_B*env->T) ));
4199 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
4201 return p->knot * exp(F*p->dx / (k_B*env->T) );
4205 <<bell k structure create/destroy function declarations>>=
4206 void *string_create_bell_param_t(char **param_strings);
4207 void destroy_bell_param_t(void *p);
4210 <<bell k structure create/destroy functions>>=
4211 bell_param_t *create_bell_param_t(double knot, double dx)
4213 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
4214 assert(ret != NULL);
4220 void *string_create_bell_param_t(char **param_strings)
4222 return create_bell_param_t(atof(param_strings[0]), atof(param_strings[1]));
4225 void destroy_bell_param_t(void *p /* bell_param_t * */)
4232 <<bell k global declarations>>=
4233 extern int num_bell_params;
4234 extern char *bell_param_descriptions[];
4235 extern char bell_param_string[];
4239 int num_bell_params = 2;
4240 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4241 char bell_param_string[]="3.3e-4,0.25e-9";
4244 <<bell k model getopt>>=
4245 {"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}
4247 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4248 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4251 \subsection{Linker stiffness corrected Bell model}
4254 Walton et. al showed that the Bell model constant force approximation
4255 doesn't hold for biotin-streptavadin binding, and I extended their
4256 results to I27 for the 2009 Biophysical Society Annual
4257 Meeting\cite{walton08,king09}. More details to follow when I get done
4258 with the conference...
4260 We adjust Bell's model to
4262 k = k_0 \cdot e^{\frac{F dx - \frac{1}{2} \kappa dx^2}{k_B T}} \;,
4264 where $k_0$ is the unforced unfolding rate,
4265 $F$ is the applied unfolding force,
4266 $\kappa$ is the effective spring constant,
4267 $dx$ is the distance to the transition state, and
4268 $k_B T$ is the thermal energy\citep{schlierf06}.
4270 <<kbell k functions>>=
4271 <<kbell k function>>
4272 <<kbell k structure create/destroy functions>>
4275 <<kbell k declarations>>=
4276 <<kbell k function declaration>>
4277 <<kbell k structure create/destroy function declarations>>
4278 <<kbell k global declarations>>
4282 <<kbell k structure>>=
4283 typedef struct kbell_param_struct {
4284 double knot; /* transition rate k_0 (frac population per s) */
4285 double dx; /* `distance to transition state' (nm) */
4289 <<kbell k function declaration>>=
4290 double kbell_k(double F, environment_t *env, void *kbell_params);
4292 <<kbell k function>>=
4293 double kbell_k(double F, environment_t *env, void *kbell_params)
4294 { /* Returns the rate constant k in frac pop per s.
4295 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4296 * uses global k_B in J/K */
4297 kbell_param_t *p = (kbell_param_t *) kbell_params;
4298 assert(F >= 0); assert(env->T > 0); assert(env->stiffness >= 0);
4300 assert(p->knot > 0); assert(p->dx >= 0);
4301 return p->knot * exp((F*p->dx - 0.5*env->stiffness*(p->dx)*(p->dx))/(k_B*env->T) );
4305 <<kbell k structure create/destroy function declarations>>=
4306 void *string_create_kbell_param_t(char **param_strings);
4307 void destroy_kbell_param_t(void *p);
4310 <<kbell k structure create/destroy functions>>=
4311 kbell_param_t *create_kbell_param_t(double knot, double dx)
4313 kbell_param_t *ret = (kbell_param_t *) malloc(sizeof(kbell_param_t));
4314 assert(ret != NULL);
4320 void *string_create_kbell_param_t(char **param_strings)
4322 return create_kbell_param_t(atof(param_strings[0]), atof(param_strings[1]));
4325 void destroy_kbell_param_t(void *p /* kbell_param_t * */)
4332 <<kbell k global declarations>>=
4333 extern int num_kbell_params;
4334 extern char *kbell_param_descriptions[];
4335 extern char kbell_param_string[];
4338 <<kbell k globals>>=
4339 int num_kbell_params = 2;
4340 char *kbell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4341 char kbell_param_string[]="3.3e-4,0.25e-9";
4344 <<kbell k model getopt>>=
4345 {"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}
4347 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4348 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4351 \subsection{Kramer's model}
4354 Kramer's model gives $k$ as
4356 % \frac{1}{k} = \frac{1}{D}
4357 % \int_{x_\text{min}}^{x_\text{max}}
4358 % e^{\frac{-U_F(x)}{k_B T}}
4359 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4362 %where $D$ is the diffusion constant,
4363 %$U_F(x)$ is the free energy along the unfolding coordinate, and
4364 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4365 % before the minimum and after the maximum, respectively \citep{schlierf06}.
4367 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
4369 where $D$ is the diffusion constant,
4371 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
4373 are length scales where
4374 $x_c(F)$ is the location of the bound state and
4375 $x_{ts}(F)$ is the location of the transition state,
4376 $E(F,x)$ is the free energy, and
4377 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
4378 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
4379 \citet{evans97} Eqn.~3).
4381 <<kramers k functions>>=
4382 <<kramers k function>>
4383 <<kramers k structure create/destroy functions>>
4386 <<kramers k declarations>>=
4387 <<kramers k function declaration>>
4388 <<kramers k structure create/destroy function declarations>>
4389 <<kramers k global declarations>>
4393 <<kramers k structure>>=
4394 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
4395 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
4397 typedef struct kramers_param_struct {
4398 double D; /* diffusion rate (um/s) */
4405 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
4406 //kramers_x_func_t *xts; /* function returning transition x (nm) */
4407 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
4408 //kramers_E_func_t *E; /* function returning E (J) */
4409 //double *E_params; /* parametrize the energy landscape */
4410 //int n_E_params; /* length of E_params array */
4414 <<kramers k function declaration>>=
4415 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
4416 extern double kramers_k(double F, environment_t *env, void *kramers_params);
4418 <<kramers k function>>=
4419 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
4421 static char input[160]={0};
4422 static char output[80]={0};
4423 /* setup the environment */
4424 fprintf(in, "F = %g; T = %g\n", F, T);
4425 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4426 string_eval(in, out, input, 80, output);
4427 return atof(output);
4430 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
4432 static char input[160]={0};
4433 static char output[80]={0};
4434 /* setup the environment */
4435 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
4436 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4437 string_eval(in, out, input, 80, output);
4438 return atof(output);
4441 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
4443 kramers_param_t *p = (kramers_param_t *) kramers_params;
4444 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
4447 double kramers_k(double F, environment_t *env, void *kramers_params)
4448 { /* Returns the rate constant k in frac pop per s.
4449 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4450 * uses global k_B in J/K */
4451 kramers_param_t *p = (kramers_param_t *) kramers_params;
4452 double kbT, xc, xts, lc, lts, Eb;
4453 assert(F >= 0); assert(env->T > 0);
4456 assert(p->xc != NULL); assert(p->xts != NULL);
4457 assert(p->ddEdxx != NULL); assert(p->E != NULL);
4458 assert(p->in != NULL); assert(p->out != NULL);
4460 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
4461 if (xc == -1.0) return -1.0; /* error (Too much force) */
4462 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
4463 if (xc == -1.0) return -1.0; /* error (Too much force) */
4464 lc = sqrt(2.0*M_PI*kbT /
4465 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
4466 lts = sqrt(-2.0*M_PI*kbT/
4467 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
4468 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
4469 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
4470 //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));
4471 return p->D/(lc*lts) * exp(-Eb/kbT);
4475 <<kramers k structure create/destroy function declarations>>=
4476 void *string_create_kramers_param_t(char **param_strings);
4477 void destroy_kramers_param_t(void *p);
4480 <<kramers k structure create/destroy functions>>=
4481 kramers_param_t *create_kramers_param_t(double D,
4482 char *xc_fn, char *xts_fn,
4483 char *E_fn, char *ddEdxx_fn,
4484 char *global_define_string)
4485 // kramers_x_func_t *xc, kramers_x_func_t *xts,
4486 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
4487 // double *E_params, int n_E_params)
4489 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
4490 assert(ret != NULL);
4491 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4492 assert(ret->xc != NULL);
4493 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4494 assert(ret->xts != NULL);
4495 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
4496 assert(ret->E != NULL);
4497 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
4498 assert(ret->ddEdxx != NULL);
4500 strcpy(ret->xc, xc_fn);
4501 strcpy(ret->xts, xts_fn);
4502 strcpy(ret->E, E_fn);
4503 strcpy(ret->ddEdxx, ddEdxx_fn);
4504 //ret->E_params = E_params;
4505 //ret->n_E_params = n_E_params;
4506 string_eval_setup(&ret->in, &ret->out);
4507 fprintf(ret->in, "kB = %g\n", k_B);
4508 fprintf(ret->in, "%s\n", global_define_string);
4512 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
4513 void *string_create_kramers_param_t(char **param_strings)
4515 return create_kramers_param_t(atof(param_strings[0]),
4523 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
4525 kramers_param_t *pk = (kramers_param_t *)p;
4527 string_eval_teardown(&pk->in, &pk->out);
4529 // free(pk->E_params);
4535 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
4536 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.
4537 We follow \citet{shillcock98} and use
4539 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
4540 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
4543 where TODO, variable meanings.%\citep{shillcock98}.
4544 The first and second derivatives of $E(F,E_b,x,x_s)$ are
4546 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\\
4547 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
4549 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
4550 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
4551 The roots are therefor at
4553 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
4554 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
4555 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
4558 <<kramers k global declarations>>=
4559 extern int num_kramers_params;
4560 extern char *kramers_param_descriptions[];
4561 extern char kramers_param_string[];
4564 <<kramers k globals>>=
4565 int num_kramers_params = 6;
4566 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"};
4567 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)}";
4569 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
4570 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
4571 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
4572 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{app.adaptive_dt} reduces it's timestep appropriately.
4573 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
4574 It works so long as [[val_1]] is not [[false]].
4576 <<kramers k model getopt>>=
4577 {"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}
4580 \subsection{Kramer's model (natural cubic splines)}
4581 \label{sec.kramers_integ}
4583 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
4584 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
4585 \citet{schlierf06} Eqn.~2)
4587 \frac{1}{k} = \frac{1}{D}
4588 \int_{x_\text{min}}^{x_\text{max}}
4589 e^{\frac{U_F(x)}{k_B T}}
4590 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4593 where $D$ is the diffusion constant,
4594 $U_F(x)$ is the free energy along the unfolding coordinate, and
4595 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4596 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
4598 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
4599 interpolating between them with cubic splines.
4600 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
4602 <<kramers integ k functions>>=
4603 <<spline functions>>
4604 <<kramers integ A k functions>>
4605 <<kramers integ B k functions>>
4606 <<kramers integ k function>>
4607 <<kramers integ k structure create/destroy functions>>
4610 <<kramers integ k declarations>>=
4611 <<kramers integ k function declaration>>
4612 <<kramers integ k structure create/destroy function declarations>>
4613 <<kramers integ k global declarations>>
4617 <<kramers integ k structure>>=
4618 typedef double func_t(double x, void *params);
4619 typedef struct kramers_integ_param_struct {
4620 double D; /* diffusion rate (um/s) */
4621 double x_min; /* integration bounds */
4623 func_t *E; /* function returning E (J) */
4624 void *E_params; /* parametrize the energy landscape */
4625 destroy_data_func_t *destroy_E_params;
4627 double F; /* for passing into GSL evaluations */
4629 } kramers_integ_param_t;
4632 <<spline param structure>>=
4633 typedef struct spline_param_struct {
4634 int n_knots; /* natural cubic spline knots for U(x) */
4637 gsl_spline *spline; /* GSL spline parameters */
4638 gsl_interp_accel *acc; /* GSL spline acceleration data */
4642 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
4644 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4646 <<kramers integ A k functions>>=
4647 double kramers_integ_k_integrandA(double x, void *params)
4649 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4650 double U = (*p->E)(x, p->E_params);
4651 double Fx = p->F * x;
4652 double kbT = k_B * p->env->T;
4653 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
4654 return exp(-(U-Fx)/kbT);
4656 double kramers_integ_k_integralA(double x, void *params)
4658 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4660 double result, abserr;
4661 size_t neval; /* for qng */
4662 static gsl_integration_workspace *w = NULL;
4663 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4664 f.function = &kramers_integ_k_integrandA;
4666 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4667 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4668 //fprintf(stderr, "integralA = %g\n", result);
4674 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
4675 \int_{x_\text{min}}^{x_\text{max}}
4676 e^{\frac{U_F(x)}{k_B T}}
4677 \text{Integral}_A(x)
4680 <<kramers integ B k functions>>=
4681 double kramers_integ_k_integrandB(double x, void *params)
4683 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4684 double U = (*p->E)(x, p->E_params);
4685 double Fx = p->F * x;
4686 double kbT = k_B * p->env->T;
4687 double intA = kramers_integ_k_integralA(x, params);
4688 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
4689 return exp((U-Fx)/kbT)*intA;
4691 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
4693 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4695 double result, abserr;
4696 size_t neval; /* for qng */
4697 static gsl_integration_workspace *w = NULL;
4698 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4699 f.function = &kramers_integ_k_integrandB;
4703 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4704 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4705 //fprintf(stderr, "integralB = %g\n", result);
4710 With the integrals taken care of, evaluating $k$ is simply
4712 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
4714 <<kramers integ k function declaration>>=
4715 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
4717 <<kramers integ k function>>=
4718 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
4719 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
4720 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
4721 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
4725 <<kramers integ k structure create/destroy function declarations>>=
4726 void *string_create_kramers_integ_param_t(char **param_strings);
4727 void destroy_kramers_integ_param_t(void *p);
4730 <<kramers integ k structure create/destroy functions>>=
4731 kramers_integ_param_t *create_kramers_integ_param_t(double D,
4732 double xmin, double xmax, func_t *E,
4734 destroy_data_func_t *destroy_E_params)
4736 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
4737 assert(ret != NULL);
4742 ret->E_params = E_params;
4743 ret->destroy_E_params = destroy_E_params;
4747 void *string_create_kramers_integ_param_t(char **param_strings)
4749 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
4750 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
4751 void *E_params = string_create_spline_param_t(param_strings+1);
4752 return create_kramers_integ_param_t(atof(param_strings[0]),
4753 atof(param_strings[2]),
4754 atof(param_strings[3]),
4755 &spline_eval, E_params,
4756 destroy_spline_param_t);
4759 void destroy_kramers_integ_param_t(void *params)
4761 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4764 (*p->destroy_E_params)(p->E_params);
4770 Finally we have the GSL spline wrappers:
4771 <<spline functions>>=
4773 <<create/destroy spline>>
4776 <<spline function>>=
4777 double spline_eval(double x, void *spline_params)
4779 spline_param_t *p = (spline_param_t *)(spline_params);
4780 return gsl_spline_eval(p->spline, x, p->acc);
4784 <<create/destroy spline>>=
4785 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
4787 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
4788 assert(ret != NULL);
4789 ret->n_knots = n_knots;
4792 ret->acc = gsl_interp_accel_alloc();
4793 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
4794 gsl_spline_init(ret->spline, x, y, n_knots);
4798 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
4799 void *string_create_spline_param_t(char **param_strings)
4803 double *x=NULL, *y=NULL;
4804 /* break into ordered pairs */
4805 parse_list_string(param_strings[0], SEP, '(', ')',
4806 &num_knots, &knot_coords);
4807 x = (double *)malloc(sizeof(double)*num_knots);
4809 y = (double *)malloc(sizeof(double)*num_knots);
4811 for (i=0; i<num_knots; i++) {
4812 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
4813 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
4816 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
4818 free_parsed_list(num_knots, knot_coords);
4819 return create_spline_param_t(num_knots, x, y);
4822 void destroy_spline_param_t(void *params)
4824 spline_param_t *p = (spline_param_t *)params;
4827 gsl_spline_free(p->spline);
4829 gsl_interp_accel_free(p->acc);
4839 <<kramers integ k global declarations>>=
4840 extern int num_kramers_integ_params;
4841 extern char *kramers_integ_param_descriptions[];
4842 extern char kramers_integ_param_string[];
4845 <<kramers integ k globals>>=
4846 int num_kramers_integ_params = 4;
4847 char *kramers_integ_param_descriptions[] = {
4848 "diffusion D, m^2 / s",
4849 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
4850 "starting integration bound x_min, m",
4851 "ending integration bound x_max, m"};
4852 //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";
4853 //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";
4854 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";
4855 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
4856 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
4857 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
4860 <<kramers integ k model getopt>>=
4861 {"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}
4863 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4864 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4866 \subsection{Unfolding model implementation}
4870 <<k model includes>>
4871 #include "k_model.h"
4872 <<k model internal definitions>>
4874 <<k model functions>>
4877 <<k model includes>>=
4878 #include <assert.h> /* assert() */
4879 #include <stdlib.h> /* NULL, malloc() */
4880 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4881 #include <stdio.h> /* fprintf(), stdout */
4882 #include <string.h> /* strlen(), strcpy() */
4883 #include <gsl/gsl_integration.h>
4884 #include <gsl/gsl_spline.h>
4889 <<k model internal definitions>>=
4890 <<const k structure>>
4891 <<bell k structure>>
4892 <<kbell k structure>>
4893 <<kramers k structure>>
4894 <<kramers integ k structure>>
4895 <<spline param structure>>
4898 <<k model globals>>=
4903 <<kramers k globals>>
4904 <<kramers integ k globals>>
4907 <<k model functions>>=
4908 <<null k functions>>
4909 <<const k functions>>
4910 <<bell k functions>>
4911 <<kbell k functions>>
4912 <<kramers k functions>>
4913 <<kramers integ k functions>>
4916 \subsection{Unfolding model unit tests}
4918 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
4919 <<check-k-model.c>>=
4920 <<k model check includes>>
4921 <<check relative error>>
4923 <<k model test suite>>
4924 <<main check program>>
4927 <<k model check includes>>=
4928 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
4929 #include <stdio.h> /* sprintf() */
4930 #include <assert.h> /* assert() */
4931 #include <math.h> /* exp() */
4934 #include "k_model.h"
4937 <<k model test suite>>=
4940 <<k model suite definition>>
4943 <<k model suite definition>>=
4944 Suite *test_suite (void)
4946 Suite *s = suite_create ("k model");
4947 <<const k test case defs>>
4948 <<bell k test case defs>>
4950 <<const k test case adds>>
4951 <<bell k test case adds>>
4956 \subsubsection{Constant}
4958 <<const k test case defs>>=
4959 TCase *tc_const_k = tcase_create("Constant k");
4962 <<const k test case adds>>=
4963 tcase_add_test(tc_const_k, test_const_k_create_destroy);
4964 tcase_add_test(tc_const_k, test_const_k_over_range);
4965 suite_add_tcase(s, tc_const_k);
4970 START_TEST(test_const_k_create_destroy)
4972 char *knot[] = {"1","2","3","4","5","6"};
4973 char *params[] = {knot[0]};
4976 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4977 params[0] = knot[i];
4978 p = string_create_const_k_param_t(params);
4979 destroy_const_k_param_t(p);
4984 START_TEST(test_const_k_over_range)
4986 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
4987 char *knot[] = {"1","2","3","4","5","6"};
4988 char *params[] = {knot[0]};
4991 char param_string[80];
4994 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4995 params[0] = knot[i];
4996 p = string_create_const_k_param_t(params);
4997 for ( F=Fm; F<FM; F+=dF ) {
4998 fail_unless(const_k(F, &env, p)==atof(knot[i]),
4999 "const_k(%g, %g, &{%s}) = %g != %s",
5000 F, T, knot[i], const_k(F, &env, p), knot[i]);
5002 destroy_const_k_param_t(p);
5008 \subsubsection{Bell}
5010 <<bell k test case defs>>=
5011 TCase *tc_bell = tcase_create("Bell's k");
5014 <<bell k test case adds>>=
5015 tcase_add_test(tc_bell, test_bell_k_create_destroy);
5016 tcase_add_test(tc_bell, test_bell_k_at_zero);
5017 tcase_add_test(tc_bell, test_bell_k_at_one);
5018 suite_add_tcase(s, tc_bell);
5022 START_TEST(test_bell_k_create_destroy)
5024 char *knot[] = {"1","2","3","4","5","6"};
5026 char *params[] = {knot[0], dx};
5029 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5030 params[0] = knot[i];
5031 p = string_create_bell_param_t(params);
5032 destroy_bell_param_t(p);
5037 START_TEST(test_bell_k_at_zero)
5039 double F=0.0, T=300.0;
5040 char *knot[] = {"1","2","3","4","5","6"};
5042 char *params[] = {knot[0], dx};
5045 char param_string[80];
5048 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5049 params[0] = knot[i];
5050 p = string_create_bell_param_t(params);
5051 fail_unless(bell_k(F, &env, p)==atof(knot[i]),
5052 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
5053 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
5054 destroy_bell_param_t(p);
5059 START_TEST(test_bell_k_at_one)
5062 char *knot[] = {"1","2","3","4","5","6"};
5064 char *params[] = {knot[0], dx};
5065 double F= k_B*T/atof(dx);
5068 char param_string[80];
5071 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5072 params[0] = knot[i];
5073 p = string_create_bell_param_t(params);
5074 CHECK_ERR(1e-6, atof(knot[i])*exp(1.0), bell_k(F, &env, p));
5075 //fail_unless(bell_k(F, &env, p)==atof(knot[i])*exp(1.0),
5076 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
5077 // F, T, knot[i], dx, bell_k(F, &env, p), atof(knot[i])*exp(1.0));
5078 destroy_bell_param_t(p);
5084 <<kramers k tests>>=
5087 <<kramers k test case def>>=
5090 <<kramers k test case add>>=
5093 <<k model function tests>>=
5094 <<check relative error>>
5096 START_TEST(test_something)
5098 double k=5, x=3, last_x=2.0, F;
5099 one_dim_fn_t *handlers[] = {&hooke};
5100 void *data[] = {&k};
5101 double xi[] = {0.0};
5103 int new_active_groups = 1;
5104 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
5105 fail_unless(F = k*x, NULL);
5110 \subsection{Utilities}
5112 The unfolding models can get complicated, and are sometimes hard to explain to others.
5113 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
5114 the overhead of having to understand a full simulation.
5115 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
5116 or special mode, where the operation depends on the specific model selected.
5118 <<k-model-utils.c>>=
5120 <<k model utility includes>>
5121 <<k model utility definitions>>
5122 <<k model utility getopt functions>>
5123 <<k model utility multi model E>>
5124 <<k model utility main>>
5127 <<k model utility main>>=
5128 int main(int argc, char **argv)
5130 <<k model getopt array definition>>
5131 k_model_getopt_t *model = NULL;
5136 int num_param_args; /* for INIT_MODEL() */
5137 char **param_args; /* for INIT_MODEL() */
5139 double special_xmin;
5140 double special_xmax;
5141 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
5142 &Fmax, &special_xmin, &special_xmax, &flags);
5143 if (flags & VFLAG) {
5144 printf("#initializing model %s with parameters %s\n", model->name, model->params);
5146 INIT_MODEL("utility", model, model->params, params);
5149 if (model->k == NULL) {
5150 printf("No k function for model %s\n", model->name);
5154 printf("Fmax = %g <= 0\n", Fmax);
5160 printf("#F (N)\tk (%% pop. per s)\n");
5161 for (i=0; i<=N; i++) {
5162 F = Fmax*i/(double)N;
5163 k = (*model->k)(F, &env, params);
5165 printf("%g\t%g\n", F, k);
5170 if (model->k == NULL || model->k == &bell_k) {
5171 printf("No special mode for model %s\n", model->name);
5174 if (special_xmax <= special_xmin) {
5175 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
5179 double Xrng=(special_xmax-special_xmin), x, E;
5181 printf("#x (m)\tE (J)\n");
5182 for (i=0; i<=N; i++) {
5183 x = special_xmin + Xrng*i/(double)N;
5184 E = multi_model_E(model->k, params, &env, x);
5185 printf("%g\t%g\n", x, E);
5192 if (model->destructor)
5193 (*model->destructor)(params);
5198 <<k model utility multi model E>>=
5199 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
5201 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
5203 if (k == kramers_integ_k)
5204 E = (*p->E)(x, p->E_params);
5206 E = kramers_E(0, env, model_params, x);
5212 <<k model utility includes>>=
5215 #include <string.h> /* strlen() */
5216 #include <assert.h> /* assert() */
5219 #include "string_eval.h"
5220 #include "k_model.h"
5223 <<k model utility definitions>>=
5224 <<version definition>>
5225 #define VFLAG 1 /* verbose */
5226 enum mode_t {M_K_OF_F, M_SPECIAL};
5227 <<string comparison definition and globals>>
5228 <<k model getopt definitions>>
5229 <<initialize model definition>>
5230 <<kramers k structure>>
5231 <<kramers integ k structure>>
5235 <<k model utility getopt functions>>=
5237 <<k model utility help>>
5238 <<k model utility get options>>
5241 <<k model utility help>>=
5242 void help(char *prog_name,
5244 int n_k_models, k_model_getopt_t *k_models,
5245 int k_model, double Fmax, double special_xmin, double special_xmax)
5248 printf("usage: %s [options]\n", prog_name);
5249 printf("Version %s\n\n", VERSION);
5250 printf("A utility for understanding the available unfolding models\n\n");
5251 printf("Environmental options:\n");
5252 printf("-T T\tTemperature T (currently %g K)\n", env->T);
5253 printf("-C T\tYou can also set the temperature T in Celsius\n");
5254 printf("Model options:\n");
5255 printf("-k model\tTransition rate model (currently %s)\n",
5256 k_models[k_model].name);
5257 printf("-K args\tTransition rate model argument string (currently %s)\n",
5258 k_models[k_model].params);
5259 printf("There are two output modes. In standard mode, k(F) is printed\n");
5260 printf("For example:\n");
5261 printf(" #Force (N)\tk (% pop. per s)\n");
5262 printf(" 123.456\t7.89\n");
5264 printf("In special mode, the output depends on the model.\n");
5265 printf("For models defining an energy function E(x), that function is printed\n");
5266 printf(" #Distance (m)\tE (J)\n");
5267 printf(" 0.0012\t0.0034\n");
5269 printf("-m\tChange output to standard mode\n");
5270 printf("-M\tChange output to special mode\n");
5271 printf("-F\tSet the maximum F value for the standard mode k(F) output (currently %g)\n", Fmax);
5272 printf("-x\tSet the minimum x value for the special mode E(x) output (currently %g)\n", special_xmin);
5273 printf("-X\tSet the maximum x value for the special mode E(x) output (currently %g)\n", special_xmax);
5274 printf("-V\tChange output to verbose mode\n");
5275 printf("-h\tPrint this help and exit\n");
5277 printf("Unfolding rate models:\n");
5278 for (i=0; i<n_k_models; i++) {
5279 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
5280 for (j=0; j < k_models[i].num_params ; j++)
5281 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
5282 printf("\t\tdefaults: %s\n", k_models[i].params);
5287 <<k model utility get options>>=
5288 void get_options(int argc, char **argv, environment_t *env,
5289 int n_k_models, k_model_getopt_t *k_models,
5290 enum mode_t *mode, k_model_getopt_t **model,
5291 double *Fmax, double *special_xmin, double *special_xmax,
5292 unsigned int *flags)
5294 char *prog_name = NULL;
5295 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
5297 extern char *optarg;
5298 extern int optind, optopt, opterr;
5302 /* setup defaults */
5304 prog_name = argv[0];
5305 env->T = 300.0; /* K */
5310 *special_xmax = 1e-8;
5311 *special_xmin = 0.0;
5313 while ((c=getopt(argc, argv, options)) != -1) {
5315 case 'T': env->T = atof(optarg); break;
5316 case 'C': env->T = atof(optarg)+273.15; break;
5318 k_model = index_k_model(n_k_models, k_models, optarg);
5319 *model = k_models+k_model;
5322 k_models[k_model].params = optarg;
5324 case 'm': *mode = M_K_OF_F; break;
5325 case 'M': *mode = M_SPECIAL; break;
5326 case 'F': *Fmax = atof(optarg); break;
5327 case 'x': *special_xmin = atof(optarg); break;
5328 case 'X': *special_xmax = atof(optarg); break;
5329 case 'V': *flags |= VFLAG; break;
5331 fprintf(stderr, "unrecognized option '%c'\n", optopt);
5332 /* fall through to default case */
5334 help(prog_name, env, n_k_models, k_models, k_model, *Fmax, *special_xmin, *special_xmax);
5343 \section{\LaTeX\ documentation}
5345 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
5346 The comment blocks are extracted (with nicely formatted code blocks), using
5347 <<latex makefile lines>>=
5348 $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5349 noweave -latex -index -delay $< > $@
5350 $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR)
5354 <<latex makefile lines>>=
5355 $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
5357 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5358 $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
5359 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5360 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5361 mv $(BUILD_DIR)/sawsim.pdf $@
5363 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
5364 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
5365 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
5367 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
5368 <<latex makefile lines>>=
5370 rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
5371 $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
5372 $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.tex \
5373 $(BUILD_DIR)/sawsim.bib
5374 clean_latex : semi-clean_latex
5375 rm -f $(DOC_DIR)/sawsim.pdf
5380 [[make]] is a common utility on *nix systems for managing dependencies while building software.
5381 In this case, we have several \emph{targets} that we'd like to build:
5382 the main simulation program \prog;
5383 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
5384 the documentation [[sawsim.pdf]];
5385 and the [[Makefile]] itself.
5386 Besides the generated files, there is the utility target
5387 [[clean]] for removing all generated files except [[Makefile]].
5389 % [[dist]] for generating a distributable tar file.
5391 Extract the makefile with
5392 `[[notangle -Rmakefile src/sawsim.nw | sed 's/ /\t/' > Makefile]]'.
5393 The sed is needed because notangle replaces tabs with eight spaces,
5394 but [[make]] needs tabs.
5396 # Decide what directories to put things in
5401 # Note: Cannot use, for example, `./build', because make eats the `./'
5402 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5404 # Modules (X.c, X.h) defined in the noweb file
5407 # Modules defined outside the noweb file
5408 FREE_SAWSIM_MODS = interp tavl
5410 <<list module makefile lines>>
5411 <<tension balance module makefile lines>>
5412 <<tension model module makefile lines>>
5413 <<k model module makefile lines>>
5414 <<parse module makefile lines>>
5415 <<string eval module makefile lines>>
5416 <<accel k module makefile lines>>
5418 SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS)
5420 # Everything needed for sawsim under one roof. sawsim.c must come first
5421 SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \
5422 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h)
5423 # Libraries needed to compile sawsim
5424 LIBS = gsl gslcblas m
5425 CHECK_LIBS = gsl gslcblas m check
5426 # The non-check binaries generated
5427 BINS = sawsim tension_model_utils k_model_utils sawsim_profile
5430 # Define the major targets
5431 all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
5433 view : $(DOC_DIR)/sawsim.pdf
5435 profile : $(BIN_DIR)/sawsim_profile | $(BIN_DIR)
5436 $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \
5437 -mnull -Mwlc -A0.39e-9,28e-9 -F8
5438 gprof $(BIN_DIR)/sawsim_profile gmon.out > $@
5439 check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim
5440 $(SHELL) -e -c 'for B in $^; do ./$$B; done'
5441 clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \
5442 clean_tension_model_utils \
5443 clean_k_model_utils clean_latex clean_check_sawsim
5444 rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \
5445 $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \
5446 $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
5447 $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
5448 $(BUILD_DIR)/global.h ./gmon.out
5449 $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
5451 # Various builds of sawsim
5452 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
5453 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5454 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
5455 gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5456 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
5457 gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5459 # Create the directories
5460 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
5463 # Copy over the source external to sawsim
5464 # Note: Cannot use, for example, `./build', because make eats the `./'
5465 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5467 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
5471 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
5475 ## Basic source generated with noweb
5476 # The central files sawsim.c and global.h...
5477 $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5479 $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5480 notangle -Rglobal.h $< > $@
5481 # ... and the modules
5482 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5483 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5484 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5485 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5486 $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5487 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5488 # Note: I use `_' as a space in my file names, but noweb doesn't like
5489 # that so we use `-' to name the noweb roots and substitute here.
5491 ## Building the unit-test programs
5493 $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
5494 notangle -Rchecks $< > $@
5495 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
5496 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
5497 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR)
5498 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
5499 clean_check_sawsim :
5500 rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
5501 # ... and the modules
5503 $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
5504 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \
5505 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \
5506 $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
5507 $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
5508 $$(BUILD_DIR)/global.h | $(BIN_DIR)
5509 gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
5510 $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
5512 # todo: clean up the required modules to
5513 $(CHECK_BINS:%=clean_%) :
5514 rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
5516 # Cleaning up the modules
5518 $(SAWSIM_MODS:%=clean_%) :
5519 rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h)
5521 <<tension model utils makefile lines>>
5522 <<k model utils makefile lines>>
5523 <<latex makefile lines>>
5525 Makefile : $(SRC_DIR)/sawsim.nw
5526 notangle -Rmakefile $< | sed 's/ /\t/' > $@
5528 This is fairly self-explanatory. We compile a dynamically linked
5529 version ([[sawsim]]) and a statically linked version
5530 ([[sawsim_static]]). The static version is about 9 times as big, but
5531 it runs without needing dynamic GSL libraries to link to, so it's a
5532 better format for distributing to a cluster for parallel evaluation.
5536 \subsection{Hookean springs in series}
5537 \label{app.math_hooke}
5540 The effective spring constant for $n$ springs $k_i$ in series is given by
5542 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5548 F &= k_i x_i &&\forall i \le n \\
5549 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
5550 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
5551 F &= k_1 x_1 = k_\text{eff} x \\
5552 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
5553 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5558 \addcontentsline{toc}{section}{References}
5559 \bibliography{sawsim}