1 %%%%%%%%%%%%%%%%%%%%%%%%% Start LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%
2 % we have our own preamble,
3 % use `noweave -delay` to not print noweb's automatic one
4 % -*- mode: Noweb; noweb-code-mode: c-mode -*-
5 \documentclass[letterpaper, 11pt]{article}
8 \noweboptions{smallcode}
10 \usepackage{url} % needed for natbib for underscores in urls?
11 \usepackage[breaklinks,colorlinks]{hyperref} % for internal links, \phantomsection
12 % breaklinks breaks long links
13 % colorlinks colors link text instead of boxing it
14 \usepackage{amsmath} % for \text{}, \xrightarrow[sub]{super}
15 \usepackage{amsthm} % for theorem and proof environments
16 \newtheorem{theorem}{Theorem}
18 \usepackage{subfig} % Support side-by-side figures
19 \usepackage{pgf} % Fancy graphics
20 \usepackage{tikz} % A nice, inline PGF frontend
21 \usetikzlibrary{automata} % Graph-theory library
23 \usepackage{doc} % BibTeX
24 \usepackage[super,sort&compress]{natbib} % fancy citation extensions
25 % super selects citations in superscript mode
26 % sort&compress automatically sorts and compresses compound citations (\cite{a,b,...})
28 \usepackage[mathletters]{ucs} % use math fonts to allow Greek UTF-8 in the text
29 \usepackage[utf8x]{inputenc} % allow UTF-8 characters
31 %\bibliographystyle{ieeetr} % pick the bibliography style, short and sweet
32 %\bibliographystyle{plain} % pick the bibliography style, includes dates
33 \bibliographystyle{plainnat}
34 \defcitealias{sw:noweb}{{\tt noweb}}
35 \defcitealias{galassi05}{GSL}
36 \defcitealias{sw:check}{{\tt check}}
37 \defcitealias{sw:python}{Python}
42 \textheight 9.5in % leave a bit of extra room for page numbers (at top of page)
45 \setlength{\parindent}{0pt}
46 \setlength{\parskip}{5pt}
48 % For some reason, the \maketitle wants to go on it's own page
49 % so we'll just hardcode the following in our first page.
50 %\title{Sawsim: a sawtooth protein unfolding simulator}
51 %\author{W.~Trevor King}
54 \newcommand{\prog}{[[sawsim]]}
55 \newcommand{\lang}{[[C]]}
56 \newcommand{\p}[3]{\left#1 #2 \right#3} % e.g. \p({complicated expr})
57 \newcommand{\Th}{\ensuremath{\sup{\text{th}}}} % ordinal th
58 \newcommand{\U}[1]{\ensuremath{\text{ #1}}} % units: $1\U{m/s}$
60 % single spaced lists, from Alvin Alexander
61 % http://www.devdaily.com/blog/post/latex/control-line-spacing-in-itemize-enumerate-tags/
62 \newenvironment{packed_item}{
64 \setlength{\itemsep}{1pt}
65 \setlength{\parskip}{0pt}
66 \setlength{\parsep}{0pt}
70 \nwfilename{sawsim.nw}
75 Sawsim: a sawtooth protein unfolding simulator \\
83 %%%%%%%%%%%%%%%%%%%%%%%%%%%% End LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%%
85 \section{Introduction}
87 The unfolding of multi-globular protein strings is a stochastic
88 process. In the \prog\ program, we use Monte Carlo methods to
89 simulate this unfolding through an explicit time-stepping approach.
90 We develop a program in \lang\ to simulate probability of unfolding
91 under a constant extension velocity of a chain of globular domains.
93 In order to extract physical meaning from many mechanical unfolding
94 simulations, it is necessary to compare the experimental results
95 against the results of simulations using different models and
96 parameters. The model and parameters that best match the experimental
97 data provide the ‘best’ model for that experiment.
99 Previous mechanical unfolding studies have rolled their own
100 simulations for modelling their experiment, and there has been little
101 exchange and standardization of these simulations between groups.
102 The \prog\ program attempts to fill this gap, providing a flexible,
103 extensible, \emph{community} program, to make it easier for groups
104 to share the burden and increase the transparency of data analysis.
106 ‘Sawsim’ is blend of ‘sawtooth simulation’.
108 \subsection{Related work}
112 \subsection{About this document}
114 This document was written in \citetalias{sw:noweb} with code blocks in
115 \lang\ and comment blocks in \LaTeX.
117 \subsection{System Requirements and Dependencies}
120 \subsection{Change Log}
122 Version 0.0 used only the unfolded domain WLC to determine the tension
123 and had a fixed timestep $dt$.
125 Version 0.1 added adaptive timesteps, adjusting so that the unfolding
126 probability for a given domain was constant.
128 Version 0.2 added Kramers' model unfolding rate calculations, and a
129 switch to select Bell's or Kramers' model. This induced a major
130 rewrite, introducing the tension group abstraction, and a split to
131 multiple \lang\ source files. Also added a unit-test suites for
132 sanity checking, and switched to SI units throughout.
134 Version 0.3 added integrated Kramers' model unfolding rate
135 calculations. Also replaced some of my hand-coded routines with GNU
136 Scientific Library \citepalias{galassi05} calls.
138 Version 0.4 added Python evaluation for the saddle-point Kramers'
139 model unfolding rate calculations, so that the model functions could
140 be controlled from the command line. Also added interpolating lookup
141 tables to accelerate slow unfolding rate model evaluation, and command
142 line control of tension groups.
144 Version 0.5 added the stiffness environmental parameter and it's
145 associated unfolding models.
147 Version 0.6 generalized from two state unfolding to arbitrary
148 $N$-state rate reactions. This allows simulations of
149 unfolding-refolding, metastable intermediates, etc., but required
150 reorganizing the command line interface.
152 Version 0.7 added tension model inverses, which seems to reduce
153 computation time by about a factor of 2. I was expecting more, but
154 I'll take what I can get.
156 <<version definition>>=
157 #define VERSION "0.7"
163 sawsim - program for simulating single molecule mechanical unfolding.
164 Copyright (C) 2008-2009, William Trevor King
166 This program is free software; you can redistribute it and/or
167 modify it under the terms of the GNU General Public License as
168 published by the Free Software Foundation; either version 3 of the
169 License, or (at your option) any later version.
171 This program is distributed in the hope that it will be useful, but
172 WITHOUT ANY WARRANTY; without even the implied warranty of
173 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
174 See the GNU General Public License for more details.
176 You should have received a copy of the GNU General Public License
177 along with this program; if not, write to the Free Software
178 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
181 The author may be contacted at <wking@drexel.edu> on the Internet, or
182 write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
183 Philadelphia PA 19104, USA.
196 The unfolding system is stored as a chain of \emph{domains} (Figure
197 \ref{fig.domain_chain}). Domains can be folded, globular protein
198 domains, unfolded protein linkers, AFM cantilevers, or other
199 stretchable system components. The system is modeled as graph of
200 possible domain states with transitions between them (Figure
201 \ref{fig.domain_states}).
205 \subfloat[][]{\label{fig.domain_chain}
206 \begin{tikzpicture}[->,node distance=2cm,font=\footnotesize]
207 \tikzstyle{every state}=[fill,draw=red!50,very thick,fill=red!20]
208 \node[state] (A) {domain 1};
209 \node[state] (B) [below of=A] {domain 2};
210 \node[state] (C) [below of=B] {.~.~.};
211 \node[state] (D) [below of=C] {domain $N$};
212 \node (S) [below of=D] {Surface};
213 \node (E) [above of=A] {};
215 \path[-] (A) edge (B)
216 (B) edge node [right] {Tension} (C)
219 \path[->,green] (A) edge node [right,black] {Extension} (E);
223 \subfloat[][]{\label{fig.domain_states}
224 \begin{tikzpicture}[->,node distance=2.5cm,shorten <=1pt,shorten >=1pt,font=\footnotesize]
225 \tikzstyle{every state}=[fill,draw=blue!50,very thick,fill=blue!20]
226 \node[state] (A) {cantilever};
227 \node[state] (C) [below of=A] {transition};
228 \node[state] (B) [left of=C] {folded};
229 \node[state] (D) [right of=C] {unfolded};
231 \path (B) edge [bend left] node [above] {$k_1$} (C)
232 (C) edge [bend left] node [below] {$k_1'$} (B)
233 edge [bend left] node [above] {$k_2$} (D)
234 (D) edge [bend left] node [below] {$k_2'$} (C);
237 \caption{\subref{fig.domain_chain} Extending a chain of domains. One end
238 of the chain is fixed, while the other is extended at a constant
239 speed. The domains are coupled with rigid linkers, so the domains
240 themselves must stretch to accomodate the extension.
241 \subref{fig.domain_states} Each domain exists in a discrete state. At
242 each timestep, it may transition into another state following a
243 user-defined state matrix such as this one, showing a metastable
244 transition state and an explicit ``cantilever'' domain.}
248 Each domain contributes to the total tension, which depends on the
249 chain extension. The particular model for the tension as a function
250 of extension is handled generally with function pointers. So far the
251 following models have been implemented:
253 \item Null (Appendix \ref{sec.null_tension}),
254 \item Constant (Appendix \ref{sec.const_tension}),
255 \item Hookean spring (Appendix \ref{sec.hooke}),
256 \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
257 \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
260 The tension model and parameters used for a particular domain depend
261 on whether the domain's current state. The transition rates between
262 states are also handled generally with function pointers, with
265 \item Null (Appendix \ref{sec.null_k}),
266 \item Bell's model (Appendix \ref{sec.bell}),
267 \item Spring constant corrected Bell model (Appendix \ref{sec.kbell}),
268 \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
269 \item Kramers' saddle point model (Appendix \ref{sec.kramers}).
272 The unfolding of the chain domains is modeled in two stages. First
273 the tension of the chain is calculated. Then the tension (assumed to
274 be constant over the length of the chain, see Section
275 \ref{sec.timescales}), is applied to each domain, and used to
276 calculate the probability of that domain changing state in the next
277 timestep $dt$. Then the time is stepped forward, and the process is
278 repeated. The simulation is complete when a pre-set time limit
279 [[t_max]] is reached or a pre-selected state [[stop_state]] is empty.
281 int main(int argc, char **argv)
283 <<tension model getopt array definition>>
284 <<k model getopt array definition>>
286 /* variables initialized in get_options() */
287 list_t *state_list=NULL, *transition_list=NULL;
288 environment_t env={0};
289 double P, t_max, dt_max, v, x_step;
290 state_t *stop_state=NULL;
292 /* variables used in the simulation loop */
293 double t=0, x=0, dt, F; /* time, distance, timestep, force */
294 int transition=1; /* boolean marking a transition for tension optimization */
296 get_options(argc, argv, &P, &t_max, &dt_max, &v, &x_step, &stop_state, &env,
297 NUM_TENSION_MODELS, tension_models, NUM_K_MODELS, k_models,
298 &state_list, &transition_list);
301 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Distance (m)\tForce (N)\tStiffness (N/m)\n");
302 if (flags & FLAG_OUTPUT_TRANSITION_FORCES) printf("#Force (N)\tInitial state\tFinal state\n");
303 while ((t_max > 0 && t <= t_max) || (stop_state != NULL && stop_state->num_domains > 0)) {
304 F = find_tension(state_list, &env, x, transition);
306 dt = determine_dt(state_list, transition_list, &env, x, P, dt_max, v);
308 dt = determine_dt(state_list, transition_list, &env, x, P, dt_max, 0);
309 transition=random_transitions(transition_list, F, dt, &env);
310 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness);
315 if (dt == dt_max) { /* step completed */
318 } else { /* still working on this step */
323 destroy_state_list(state_list);
324 destroy_transition_list(transition_list);
328 @ The meat of the simulation is bundled into the three functions
329 [[find_tension]], [[determine_dt]], and [[random_transitions]].
330 [[find_tension]] is discussed in Section \ref{sec.find_tension},
331 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
332 [[random_transitions]] in Section \ref{sec.transition_rate}.
334 The stretched distance is extended in one of two modes depending on
335 [[xstep]]. If [[xstep == 0]], then the pulling is continuous, at
336 least within the limits of the inherent discretization of a time
337 stepped approach. At any rate, the timestep is limited by the maximum
338 allowed timestep [[dt_max]] and the maximum allowed unfolding
339 probability in a given step [[P]]. In the second mode [[xstep > 0]],
340 and the end to end distance increases in discrete steps of that
341 length. The time between steps is chosen to maintain an average
342 unfolding velocity [[v]]. This approach is less work to simulate than
343 smooth pulling and also closer to the reality of many experiments, but
344 it is practically impossible to treat analytically. With both pulling
345 styles implemented, the effects of the discretization can be easily
346 calculated, bridging the gap between theory and experiment.
348 Environmental parameters important in determining reaction rates and
349 tensions (e.g.\ temperature, pH) are stored in a single structure to
350 facilitate extension to more complicated models in the future.
351 <<environment definition>>=
352 typedef struct environment_struct {
362 #define DOUBLE_PRECISION 1e-12
364 <<environment definition>>
365 <<one dimensional function definition>>
366 <<create/destroy definitions>>
368 #endif /* GLOBAL_H */
370 Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
371 included multiple times.
373 \section{Simulation functions}
375 <<simulation functions>>=
379 <<does domain transition>>
380 <<randomly transition domains>>
384 \label{sec.find_tension}
386 Because the stretched system may be made up of several parts (folded
387 domains, unfolded domains, spring-like cantilever, \ldots), we will
388 assign the domains to states with tension models and groups. The
389 models are used to determine the tension of all the domain in a given
390 state following some model (Hookean spring, WLC, \ldots). The domains
391 are assumed to be commutative, so ordering is ignored. The
392 interactions between the groups are assumed to be linear, but the
393 interactions between domains of the same group need not be. Each
394 model has a tension handler function, which gives the tension $F$
395 needed to stretch a given group of domains a distance $x$.
397 GROUPS ARE CURRENTLY DISABLED.
399 Groups represent collections of states which the model should treat as
400 a single entity. To understand the purpose of groups, consider a
401 system with two unfolded domain states (e.g.\ for two different
402 proteins) that were both modeled as WLCs. If both unfolded states had
403 the same persistence length, it would be wise to assign them to the
404 same group. This would reduce both the computational cost of
405 balancing the tension and the error associated with the inter-group
406 interaction linearity. Note that group numbers only matter
407 \emph{within} model classes, since grouping domains with seperate
408 models doesn't make sense.
410 <<find tension definitions>>=
411 #define NUM_TENSION_MODELS 5
415 <<tension handler array global declaration>>=
416 extern one_dim_fn_t *tension_handlers[];
419 <<tension handler array global>>=
420 one_dim_fn_t *tension_handlers[] = {
422 &const_tension_handler,
430 <<tension model global declarations>>=
431 <<tension handler array global declaration>>
434 <<tension handler types>>=
435 typedef struct tension_handler_data_struct {
436 list_t *group_tension_params; /* one item for each domain in the group */
439 } tension_handler_data_t;
443 After sorting the chain into separate groups $G_i$, with tension
444 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
445 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
446 \forall i,j$. Note that there may be several states within each
447 group. [[find_tension]] is basically a wrapper to feed proper input
448 into the [[tension_balance]] function. [[transition]] should be set
449 to 0 if there were no transitions since the previous call to
450 [[find_tension]] to support some optimizations that assume a small
451 increase in tension between steps (only true if no transition has
452 occured). While were messing with the tension handlers, we also grab
453 a measure of the current linker stiffness for the environmental
454 variable, which is needed by some of the unfolding rate models
455 (e.g.\ the linker-stiffness corrected Bell model, Appendix
458 double find_tension(list_t *state_list, environment_t *env, double x, int transition)
460 static int num_active_groups=0;
461 static one_dim_fn_t **per_group_handlers = NULL;
462 static one_dim_fn_t **per_group_inverse_handlers = NULL;
463 static void **per_group_data = NULL; /* array of pointers to tension_handler_data_t */
464 static double last_x;
465 tension_handler_data_t *tdata;
470 fprintf(stderr, "%s:\tfinding tension at distance %g%s\n", __FUNCTION__, x, transition ? " after a transition" : "");
473 if (transition != 0 || num_active_groups == 0) { /* setup statics */
474 /* get new statics, freeing the old ones if they aren't NULL */
475 get_active_groups(state_list, &num_active_groups, &per_group_handlers, &per_group_inverse_handlers, &per_group_data);
477 /* fill in the group handlers and params */
479 fprintf(stderr, "%s:\t%d groups, tension handler array %p, inverse tension handler array %p, handler data array %p, (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_inverse_handlers, per_group_data, per_group_data[0]);
481 for (i=0; i<num_active_groups; i++) {
482 tdata = (tension_handler_data_t *) per_group_data[i];
484 fprintf(stderr, "%s:\tactive group %d of %d (data %p, param list %p)", __FUNCTION__, i, num_active_groups, tdata, tdata->group_tension_params);
485 list_print(stderr, tdata->group_tension_params, " parameters");
488 /* tdata->persist continues unchanged */
491 } /* else continue with the current statics */
493 F = tension_balance(num_active_groups, per_group_handlers, per_group_inverse_handlers, (void **)per_group_data, last_x, x, &env->stiffness);
499 @ For the moment, we will restrict our group boundaries to a single
500 dimension, so $\sum_i x_i = x$, our total extension (it is this
501 restriction that causes the groups to interact linearly). We'll also
502 restrict our extensions to all be positive. With these restrictions,
503 the problem of balancing the tensions reduces to solving a system of
504 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
505 the number of active groups. In general this can be fairly
506 complicated, but without much loss of practicality we can restrict
507 ourselves to strictly monotonically increasing, non-negative tension
508 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
509 simpler. For example, it guarantees the existence of a unique, real
510 solution for finite forces. See Appendix \ref{sec.tension_balance}
511 for the tension balancing implementation.
513 Each group must define a parameter structure if the tension model is
515 a function to create the parameter structure,
516 a function to destroy the parameter structure, and
517 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as it's data argument.
518 The tension-specific parameter structures are contained in the domain
519 group field. See the Section \ref{sec.model_selection} for a
520 discussion on the command line framework. See the worm-like chain
521 implementation for an example (Section \ref{sec.wlc}).
523 The major limitation of our approach is that all of our tension models
524 are Markovian (memory-less), which is not really a problem for the
525 chains (see \citet{evans99} for WLC thermalization rate calculations).
527 \subsection{Transition rate}
528 \label{sec.transition_rate}
530 Each state transition is modeled as unimolecular, first order reaction
532 \text{State 1} \xrightarrow{k} \text{State 2}
534 With the rate constant $k$ defined by
536 \frac{d[\text{State 2}]}{dt} = -\frac{d[\text{State 1}]}{dt} = k [\text{State 1}].
538 So the probability of a given domain transitioning in a short time
544 <<does domain transition>>=
545 int domain_transitions(double F, double dt, environment_t *env, int num_domains, transition_t *transition)
546 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, number of 'reactant' domains, pointer to transition structure */
548 k = accel_k(transition->k, F, env, transition->k_params);
549 //(*transition->k)(F, env, domain->k_params);
550 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
551 return happens(k*dt*num_domains); /* N dice rolls for prob. k*dt event */
553 @ [[happens]] is a random decision making function defined in Appendix \ref{sec.utils}.
555 <<randomly transition domains>>=
556 int random_transitions(list_t *transition_list, double tension,
557 double dt, environment_t *env)
558 { /* returns 1 if there was a transition and 0 otherwise */
559 while (transition_list != NULL) {
560 if (T(transition_list)->initial_state->num_domains > 0
561 && domain_transitions(tension, dt, env, T(transition_list)->initial_state->num_domains, T(transition_list))) {
562 if (flags & FLAG_OUTPUT_TRANSITION_FORCES)
563 fprintf(stdout, "%g\t%s\t%s\n", tension, T(transition_list)->initial_state->name, T(transition_list)->final_state->name);
564 T(transition_list)->initial_state->num_domains--;
565 T(transition_list)->final_state->num_domains++;
566 return 1; /* our one domain has transitioned, stop looking for others */
568 transition_list = transition_list->next;
574 \subsection{Timescales}
575 \label{sec.timescales}
577 The simulation assumes that chain equilibration is instantaneous
578 relative to the simulation timestep $dt$, so that tension is uniform
579 along the chain. The quality of this assumption depends on your
580 particular chain. For example, a damped spring thermalizes on a
581 timescale of order $\tau = 1/\gamma$ where the damping ratio $\gamma
582 \equiv \omega_0(-\zeta \pm \sqrt{\zeta^2-1}$, the natural angular
583 frequency $omega_0 \equiv \sqrt{k/m}$, $\zeta \equiv b/2\sqrt{km}$,
584 and $b$ sets the drag $F_b = -bv$. For our cantilevers $\tau$ is
585 on the order of milliseconds, which is longer than a timestep.
586 However, the approximation is still reasonable, because there is
587 rarely unfolding during the cantilever return. You could, of course,
588 take cantilever, WLC, etc.\ response times into effect, but that
589 would significantly complicate a simulation that seems to work
590 well enough without bothering :p. For WLC and FJC relaxation times,
591 see the Appendix of \citet{evans99} and \citet{kroy07}.
593 Besides assuming our timestep is much greater than equilibration
594 times, we also force it to be short enough so that only one domain may
595 unfold in a given timestep. For an unfolding timescale of $dt_u =
596 1/k$ we adapt our timesteps to keep $dt \ll dt_u$, so the probability
597 of two domains unfolding in the same timestep is negligible. This
598 approach breaks down as the adaptive timestep scheme approaches $dt
599 \sim dt_u$, but $dt_u \sim 1\U{$\mu$s}$ for Ig27-like proteins
600 \citep{klimov00}, so this shouldn't be a problem. To reassure
601 yourself, you can ask the simulator to print the smallest timestep
602 that was used with TODO.
604 \subsection{Adaptive timesteps}
605 \label{sec.adaptive_dt}
607 We'd like to pick $dt$ so the probability of changing state
608 (e.g.\ unfolding) in the next timestep is small. If we don't adapt our
609 timestep, we also risk breaking our approximation that $F(x' \in [x,
610 x+v \cdot dt]) = F(x)$ or that only one domain changes state in a
611 given timestep. The simulation would have been accurate for
612 sufficiently small fixed timesteps, but adaptive timesteps will allow
613 us to move through low tension regions in fewer steps, leading to a
614 more efficient simulation.
616 Consider the two state folded $\rightarrow$ unfolded transition.
617 Because $F(x)$ is monotonically increasing between unfoldings,
618 excessively large timesteps will lead to erroneously small unfolding
619 rates (an thus, higher average unfolding force).
621 The actual adaptive timestep implementation is not particularly
622 interesting, since we are only required to reduce $dt$ to somewhere
623 below a set threshold, so I've removed it to Appendix
624 \ref{sec.adaptive_dt_implementation}.
630 The models provide the physics of the simulation, but the simulation
631 allows interchangeable models, and we are currently focusing on the
632 simulation itself, so we remove the actual model implementation to the
635 The tension models are defined in Appendix \ref{sec.tension_models}
636 and unfolding models are defined in Appendix \ref{sec.k_models}.
639 #define k_B 1.3806503e-23 /* J/K */
643 \section{Command line interface}
645 <<get option functions>>=
647 <<index tension model>>
652 \subsection{Model selection}
653 \label{sec.model_selection}
655 The main difficulty with the command line interface in \prog\ is
656 developing an intuitive method for accessing the possibly large number
657 of available models. We'll treat this generally by defining an array
658 of available models, containing their important parameters.
660 <<get options definitions>>=
661 <<model getopt definitions>>
664 <<create/destroy definitions>>=
665 typedef void *create_data_func_t(char **param_strings);
666 typedef void destroy_data_func_t(void *param_struct);
669 <<model getopt definitions>>=
670 <<tension model getopt definitions>>
671 <<k model getopt definitions>>
674 In order to choose models, we need a method of assembling a reaction
675 graph such as that in Figure \ref{fig.domain_states} and population
676 the graph with a set of domains. First we declare the domain states
679 [[-s <state-name>,<model>[,<group>] -S <model-params> [-N <initial-population]]] \\
683 [[-s unfolded,wlc,1 -S 1e-8,4e-10 -N 5]]
685 This creates the state named [[unfolded]], using the WLC model and one
686 as a group number (defaults to 0, see Section \ref{sec.find_tension}).
687 The tension parameters are then set to [[1e-8,4e-10]], which in the
688 case of the WLC are the contour and persistence lengths respectively.
689 Finally we populate the state by adding five domains to the state.
690 The name of the state is importent for identifying states later.
692 Once the domains have been created and populated, you can added
693 transition rates following
695 [[-k <initial-state>,<final-state>,<model> -K <model-params>]] \\
699 [[-k folded,unfolded,bell -K 3.3e-4,0.25e-9]]
701 This creates a reaction going from the [[folded]] state to the
702 [[unfolded]] state, with the rate constant given by the Bell model
703 (Appendix \ref{sec.bell}). The unfolding parameters are then set to
704 [[3.3e-4,0.25e-9]], which in the case of the Bell model are the
705 unforced rate and transition state distance respectively.
707 \subsubsection{Tension}
709 To access the various tension models from the command line, we define
710 a structure that contains all the neccessary information about the
712 <<tension model getopt definitions>>=
713 typedef struct tension_model_getopt_struct {
714 one_dim_fn_t *handler; /* fn implementing the model on a group */
715 one_dim_fn_t *inverse_handler; /* fn implementing inverse on group */
716 char *name; /* name string identifying the model */
717 char *description; /* a brief description of the model */
718 int num_params; /* number of extra params for the model */
719 char **param_descriptions; /* descriptions of extra params */
720 char *params; /* default values for the extra params */
721 create_data_func_t *creator; /* param-string -> model-data-struct fn */
722 destroy_data_func_t *destructor; /* fn to free the model data structure */
723 } tension_model_getopt_t;
724 @ Set [[inverse_handler]] to [[NULL]] if you wish to invert using
725 bisection. Obviously, this will be slower than computing an
726 analytical inverse and more prone to rounding errors, but it may be
727 easier if a closed-form inverse does not exist.
729 <<tension model getopt array definition>>=
730 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
731 <<null tension model getopt>>,
732 <<constant tension model getopt>>,
733 <<hooke tension model getopt>>,
734 <<worm-like chain tension model getopt>>,
735 <<freely-jointed chain tension model getopt>>
739 \subsubsection{Unfolding rate}
741 <<k model getopt definitions>>=
742 #define NUM_K_MODELS 6
744 typedef struct k_model_getopt_struct {
749 char **param_descriptions;
751 create_data_func_t *creator;
752 destroy_data_func_t *destructor;
756 <<k model getopt array definition>>=
757 k_model_getopt_t k_models[NUM_K_MODELS] = {
758 <<null k model getopt>>,
759 <<const k model getopt>>,
760 <<bell k model getopt>>,
761 <<kbell k model getopt>>,
762 <<kramers k model getopt>>,
763 <<kramers integ k model getopt>>
770 void help(char *prog_name, double P, double t_max, double dt_max, double v, double x_step,
771 state_t *stop_state, environment_t *env,
772 int n_tension_models, tension_model_getopt_t *tension_models,
773 int n_k_models, k_model_getopt_t *k_models,
774 int k_model, list_t *state_list)
779 if (state_list != NULL)
780 state = S(tail(state_list));
782 printf("usage: %s [options]\n", prog_name);
783 printf("Version %s\n\n", VERSION);
784 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
786 printf("Simulation options:\n");
788 printf("-q state\tStop simulation when this state empties (- for no limit, currently %s)\n", stop_state ? stop_state->name : "-" );
789 printf("-t tmax\tMaximum simulated time (-1 for no limit, currently %g s)\n", t_max);
790 printf("-d dt\tMaximum allowed timestep dt (currently %g s)\n", dt_max);
791 printf("-P P\tTarget probability for dt (currently %g)\n", P);
792 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
793 printf("-x dx\tPulling step size (currently %g m)\n", x_step);
794 printf("\tWhen dx = 0, the pulling is continuous.\n");
795 printf("\tWhen dx > 0, the pulling occurs in discrete steps.\n");
797 printf("Environmental options:\n");
798 printf("-T T\tTemperature T (currently %g K)\n", env->T);
799 printf("-C T\tYou can also set the temperature T in Celsius\n");
801 printf("State creation options:\n");
802 printf("-s name,model[[,group],params]\tCreate a new state.\n");
803 printf("Once you've set up the state, you may populate it with domains\n");
804 printf("-N n\tAdd n domains in the current state (%s)\n", state ? state->name : "-");
806 printf("Transition creation options:\n");
807 printf("-k initial_state,final_state,model[,params]\tCreate a new transition\n");
809 printf("To double check your reaction graph, you can produce graphviz dot output\n");
810 printf("describing the current states and transitions.\n");
811 printf("-D\tPrint dot description of states and transitions and exit.\n");
813 printf("Simulation output mode options:\n");
814 printf("There are two output modes. In standard mode, only the transition\n");
815 printf("events are printed. For example:\n");
816 printf(" #Force (N)\tInitial state\tFinal state\n");
817 printf(" 123.456e-12\tfolded\tunfolded\n");
819 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
820 printf(" #Distance (m)\tForce (N)\tStiffness (N/m)\n");
821 printf(" 0.001\t0.0023\t0.002\n");
823 printf("-V\tChange output to verbose mode.\n");
824 printf("-h\tPrint this help and exit.\n");
827 printf("Tension models:\n");
828 for (i=0; i<n_tension_models; i++) {
829 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
830 for (j=0; j < tension_models[i].num_params ; j++)
831 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
832 printf("\t\tdefaults: %s\n", tension_models[i].params);
834 printf("Unfolding rate models:\n");
835 for (i=0; i<n_k_models; i++) {
836 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
837 for (j=0; j < k_models[i].num_params ; j++)
838 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
839 printf("\t\tdefaults: %s\n", k_models[i].params);
842 printf("Examples:\n");
843 printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded. Stop once there are no more folded domains.\n");
844 printf(" $ %s -v1e-6 -s cantilever,hooke,0.05 -N1 -s folded,null -N8 -s 'unfolded,wlc,{0.39e-9,28e-9}' -k 'folded,unfolded,bell,{3.3e-4,0.25e-9}' -q folded\n", prog_name);
845 printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded. Stop after 0.25 seconds\n");
846 printf(" $ %s -v1e-6 -s cantilever,hooke,0.05 -N1 -s 'folded,wlc,{1e-8,4e-10}' -N8 -s 'unfolded,wlc,1,{0.4e-9,28.1e-9}' -k 'folded,unfolded,bell,{5e-5,0.225e-9}' -t0.25\n", prog_name);
850 \subsection{Get options}
853 void get_options(int argc, char **argv,
854 double *pP, double *pT_max, double *pDt_max, double *pV,
855 double *pX_step, state_t **pStop_state, environment_t *env,
856 int n_tension_models, tension_model_getopt_t *tension_models,
857 int n_k_models, k_model_getopt_t *k_models,
858 list_t **pState_list, list_t **pTransition_list)
860 char *prog_name = NULL;
861 char c, options[] = "q:t:d:P:v:x:T:C:s:N:k:DVh";
862 int tension_model=0, k_model=0, initial_state=-1, final_state=-1;
863 char *state_name=NULL;
864 int i, n, tension_group=0;
865 list_t *temp_list=NULL;
868 void *model_params; /* for INIT_MODEL() */
869 int num_param_args; /* for INIT_MODEL() */
870 char **param_args; /* for INIT_MODEL() */
872 extern int optind, optopt, opterr;
877 for (i=0; i<n_tension_models; i++) {
878 fprintf(stderr, "tension model %i %s:\t model's handler %p\t inverse handler %p\n",
879 i, tension_models[i].name, tension_models[i].handler, tension_models[i].inverse_handler);
880 assert(tension_models[i].handler == tension_handlers[i]);
885 flags = FLAG_OUTPUT_TRANSITION_FORCES;
888 *pT_max = -1; /* s, -1 removes this limit */
889 *pDt_max = 0.001; /* s */
890 *pP = 1e-3; /* % pop per s */
891 *pV = 1e-6; /* m/s */
892 *pX_step = 0.0; /* m */
893 env->T = 300.0; /* K */
895 *pTransition_list = NULL;
897 while ((c=getopt(argc, argv, options)) != -1) {
900 if (STRMATCH(optarg, "-")) {
903 temp_list = head(*pState_list);
904 while (temp_list != NULL) {
905 if (STRMATCH(S(temp_list)->name, optarg)) {
906 *pStop_state = S(temp_list);
909 temp_list = temp_list->next;
913 case 't': *pT_max = atof(optarg); break;
914 case 'd': *pDt_max = atof(optarg); break;
915 case 'P': *pP = atof(optarg); break;
916 case 'v': *pV = atof(optarg); break;
917 case 'x': *pX_step = atof(optarg); break;
918 case 'T': env->T = atof(optarg); break;
919 case 'C': env->T = atof(optarg)+273.15; break;
921 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
922 assert(num_strings >= 2 && num_strings <= 4);
924 state_name = strings[0];
925 if (state_index(*pState_list, state_name) != -1) {
926 fprintf(stderr, "Invalid option '%s', repeated state name '%s'\n", optarg, state_name);
929 tension_model = index_tension_model(n_tension_models, tension_models, strings[1]);
930 if (num_strings == 4)
931 tension_group = atoi(strings[2]);
934 if (num_strings >= 3)
935 tension_models[tension_model].params = strings[num_strings-1];
937 fprintf(stderr, "%s:\t%s -> state %s, model %s, params %s, group %d\n", __FUNCTION__, optarg, state_name, tension_models[tension_model].name, tension_models[tension_model].params, tension_group);
939 INIT_MODEL(state_name, &tension_models[tension_model], tension_models[tension_model].params, model_params);
940 push(pState_list, create_state(state_name,
941 tension_models[tension_model].handler,
942 tension_models[tension_model].inverse_handler,
945 tension_models[tension_model].destructor,
947 *pState_list = tail(*pState_list);
949 free_parsed_list(num_strings, strings);
954 fprintf(stderr, "%s:\tadding %d domains to %s\n", __FUNCTION__, n, S(*pState_list)->name);
956 S(*pState_list)->num_domains += n;
959 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
960 assert(num_strings >= 3 && num_strings <= 4);
961 initial_state = state_index(*pState_list, strings[0]);
962 final_state = state_index(*pState_list, strings[1]);
963 k_model = index_k_model(n_k_models, k_models, strings[2]);
965 fprintf(stderr, "%s:\t%s -> %s reaction from state %s (%d) -> state %s (%d)\n", __FUNCTION__, optarg, strings[2], strings[0], initial_state, strings[1], final_state);
967 assert(initial_state != final_state);
968 if (num_strings == 4)
969 k_models[k_model].params = strings[3];
970 INIT_MODEL(state_name, &k_models[k_model], k_models[k_model].params, model_params);
971 push(pTransition_list, create_transition(list_element(*pState_list, initial_state),
972 list_element(*pState_list, final_state),
975 k_models[k_model].destructor), 1);
976 free_parsed_list(num_strings, strings);
979 print_state_graph(stdout, *pState_list, *pTransition_list);
983 flags = FLAG_OUTPUT_FULL_CURVE;
986 fprintf(stderr, "unrecognized option '%c'\n", optopt);
987 /* fall through to default case */
989 help(prog_name, *pP, *pT_max, *pDt_max, *pV, *pX_step, *pStop_state, env, n_tension_models, tension_models, n_k_models, k_models, k_model, *pState_list);
993 /* check the arguments */
994 assert(*pState_list != NULL); /* no domains! */
995 assert(*pP > 0.0); assert(*pP < 1.0);
996 assert(*pDt_max > 0.0);
998 assert(env->T > 0.0);
1000 crosslink_groups(*pState_list, *pTransition_list);
1006 <<index tension model>>=
1007 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
1010 for (i=0; i<n_models; i++)
1011 if (STRMATCH(models[i].name, name))
1013 if (i == n_models) {
1014 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
1022 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
1025 for (i=0; i<n_models; i++)
1026 if (STRMATCH(models[i].name, name))
1028 if (i == n_models) {
1029 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
1036 <<initialize model definition>>=
1037 /* requires int num_param_args and char **param_args in the current scope
1039 * INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
1040 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
1042 #define INIT_MODEL(role, model, param_string, param_pointer) \
1044 parse_list_string((param_string), SEP, '{', '}', \
1045 &num_param_args, ¶m_args); \
1046 if (num_param_args != (model)->num_params) { \
1048 "%s model %s expected %d params,\n", \
1049 (role), (model)->name, (model)->num_params); \
1051 "not the %d params in '%s'\n", \
1052 num_param_args, (param_string)); \
1053 assert(num_param_args == (model)->num_params); \
1055 if ((model)->creator) \
1056 param_pointer = (*(model)->creator)(param_args); \
1058 param_pointer = NULL; \
1059 free_parsed_list(num_param_args, param_args); \
1065 \addcontentsline{toc}{section}{Appendicies}
1067 \section{sawsim.c details}
1071 The general layout of our simulation code is:
1082 We include [[math.h]], so don't forget to link to the libm with `-lm'.
1084 #include <assert.h> /* assert() */
1085 #include <stdlib.h> /* malloc(), free(), rand() */
1086 #include <stdio.h> /* fprintf(), stdout */
1087 #include <string.h> /* strlen, strtok() */
1088 #include <math.h> /* exp(), M_PI, sqrt() */
1089 #include <sys/types.h> /* pid_t (returned by getpid()) */
1090 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
1093 #include "tension_balance.h"
1094 #include "k_model.h"
1095 #include "tension_model.h"
1097 #include "accel_k.h"
1101 <<version definition>>
1102 <<flag definitions>>
1103 <<max/min definitions>>
1104 <<string comparison definition and globals>>
1105 <<initialize model definition>>
1106 <<get options definitions>>
1107 <<state definitions>>
1108 <<transition definitions>>
1117 <<create/destroy state>>
1118 <<destroy state list>>
1120 <<create/destroy transition>>
1121 <<destroy transition list>>
1122 <<print state graph>>
1124 <<simulation functions>>
1125 <<boilerplate functions>>
1128 <<boilerplate functions>>=
1130 <<get option functions>>
1136 srand(getpid()*time(NULL)); /* seed rand() */
1140 <<flag definitions>>=
1141 /* in octal b/c of prefixed '0' */
1142 #define FLAG_OUTPUT_FULL_CURVE 01
1143 #define FLAG_OUTPUT_TRANSITION_FORCES 02
1147 static unsigned long int flags = 0;
1150 \subsection{Utilities}
1153 <<max/min definitions>>=
1154 #define MAX(a,b) ((a)>(b) ? (a) : (b))
1155 #define MIN(a,b) ((a)<(b) ? (a) : (b))
1158 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
1159 <<string comparison definition and globals>>=
1160 // Check if two strings match, return 1 if they do
1161 static char *temp_string_A;
1162 static char *temp_string_B;
1163 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
1164 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
1165 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
1166 /* +1 to also compare the '\0' */
1169 We also define a macro for our [[check]] unit testing
1170 <<check relative error>>=
1171 #define CHECK_ERR(max_err, expected, received) \
1173 fail_unless( (received-expected)/expected < max_err, \
1174 "relative error %g >= %g in %s (Expected %g, received %g)", \
1175 (received-expected)/expected, max_err, #received, \
1176 expected, received); \
1177 fail_unless(-(received-expected)/expected < max_err, \
1178 "relative error %g >= %g in %s (Expected %g, received %g)", \
1179 -(received-expected)/expected, max_err, #received, \
1180 expected, received); \
1185 int happens(double probability)
1187 assert(probability >= 0.0); assert(probability <= 1.0);
1188 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*/
1193 \subsection{Adaptive timesteps}
1194 \label{sec.adaptive_dt_implementation}
1196 $F(x)$ increases with $x$, possibly exploding (e.g.\ in the worm-like
1197 chain model), so basing the timestep on the the unfolding probability
1198 at the current tension is dangerous, and we need to search for a $dt$
1199 for which $P(F(x+v*dt)) < P_\text{target}$. There are two cases to
1200 consider. In the most common, no domains have unfolded since the last
1201 step, and we expect the next step to be slightly shorter than the
1202 previous one. In the less common, domains did unfold in the last
1203 step, and we expect the next step to be considerably longer than the
1206 double search_dt(transition_t *transition,
1208 environment_t *env, double x,
1209 double target_prob, double max_dt, double v)
1210 { /* Returns a valid timestep dt in seconds for the given transition.
1211 * Takes a list of all states, a pointer env to the environmental data,
1212 * a starting separation x in m, a target_prob between 0 and 1,
1213 * max_dt in s, stretching velocity v in m/s.
1215 double F, k, dtCur, dtU, dtUCur, dtL, dt;
1217 /* get upper bound using the current position */
1218 F = find_tension(state_list, env, x, 0); /* BUG. repeated calculation */
1219 //printf("Start with x = %g (F = %g)\n", x, F);
1220 k = accel_k(transition->k, F, env, transition->k_params);
1221 //printf("x %g\tF %g\tk %g\n", x, F, k);
1222 dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */
1224 //printf("overshot max_dt\n");
1227 if (v == 0) /* discrete stepping, so force is temporarily constant */
1230 /* set a lower bound on dt too */
1233 /* The dt determined above may produce illegitimate forces or ks.
1234 * Reduce the upper bound until we have valid ks. */
1236 F = find_tension(state_list, env, x+v*dt, 0);
1237 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1240 F = find_tension(state_list, env, x+v*dt, 0);
1242 //printf("Try for dt = %g (F = %g)\n", dt, F);
1243 k = accel_k(transition->k, F, env, transition->k_params);
1244 /* returned k may be -1.0 */
1245 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1246 while (k == -1.0) { /* reduce step until we hit a valid k */
1248 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1249 F = find_tension(state_list, env, x+v*dt, 0);
1250 //printf("Try for dt = %g (F = %g)\n", dt, F);
1251 k = accel_k(transition->k, F, env, transition->k_params);
1252 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1254 assert(dtU > 1e-14); /* timestep to valid k too small */
1255 dtUCur = target_prob / k; /* safe timestep back from x+dtU */
1257 return dt; /* dtU is safe. */
1259 /* dtU wasn't safe, lets see what would be. */
1260 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1261 dt = (dtU + dtL) / 2.0;
1262 F = find_tension(state_list, env, x+v*dt, 0);
1263 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1264 k = accel_k(transition->k, F, env, transition->k_params);
1265 dtCur = target_prob / k;
1266 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1267 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1269 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1270 dtU = dt; /* ... stepping out only dtCur would be safe */
1273 break; /* dtCur = dt */
1275 return MAX(dtUCur, dtL);
1280 To determine $dt$ for an array of potentially different folded domains,
1281 we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
1284 double determine_dt(list_t *state_list, list_t *transition_list,
1285 environment_t *env, double x,
1286 double target_prob, double dt_max, double v)
1287 { /* Returns the timestep dt in seconds.
1288 * Takes the state and transition lists, a pointer to the environment,
1289 * an extension x in m, target_prob between 0 and 1, dt_max in s,
1291 double dt=dt_max, new_dt;
1292 assert(target_prob > 0.0); assert(target_prob < 1.0);
1293 assert(dt_max > 0.0);
1294 assert(state_list != NULL);
1295 assert(transition_list != NULL);
1296 transition_list = head(transition_list);
1297 while (transition_list != NULL) {
1298 if (T(transition_list)->initial_state->num_domains > 0) {
1299 new_dt = search_dt(T(transition_list), state_list, env, x, target_prob, dt, v);
1300 dt = MIN(dt, new_dt);
1302 transition_list = transition_list->next;
1309 \subsection{State data}
1310 \label{sec.state_data}
1312 The domains exist in one of an array of possible states. Each state
1313 has a tension model which determines it's force/extension curve
1314 (possibly [[null]] for rigid states, see Appendix
1315 \ref{sec.null_tension}).
1317 Groups are, for example, for two domain states with WLC tensions; one
1318 with $p=0.4\U{nm}$, $L=10\U{nm}$ and the other with $p=0.4\U{nm}$,
1319 $L=28\U{nm}$. Since the persistence lengths are the same, you would
1320 like to add the contour lengths of all the domains in \emph{both}
1321 states, and plug that total contour length into a single WLC
1322 evaluation (see Section \ref{sec.find_tension}).
1323 <<state definitions>>=
1324 typedef struct state_struct {
1325 char *name; /* name string */
1326 one_dim_fn_t *tension_handler; /* tension handler */
1327 one_dim_fn_t *inverse_tension_handler; /* inverse tension handler */
1328 int tension_group; /* for grouping similar tension states */
1329 void *tension_params; /* pointer to tension parameters */
1330 destroy_data_func_t *destroy_tension;
1331 int num_domains; /* number of domains currently in state */
1332 /* cached lookups generated from the above data */
1333 int num_out_trans; /* length of out_trans array */
1334 struct transition_struct **out_trans; /* transitions leaving this state */
1335 int num_grouped_states; /* length of grouped states array */
1336 struct state_struct **grouped_states; /* states grouped with this one (doesn't include self) */
1339 /* get the state data for the current list node */
1340 #define S(list) ((state_t *)(list)->d)
1342 @ [[tension_params]] is a pointer to the parameters used by the
1343 handler function when determining the tension. [[destroy_tension]]
1344 points to a function for freeing [[tension_params]]. We it in the
1345 state structure so that [[destroy_state]] and will have an easier job
1348 [[create_]] and [[destroy_state]] are simple wrappers around
1349 [[malloc]] and [[free]].
1350 <<create/destroy state>>=
1351 state_t *create_state(char *name,
1352 one_dim_fn_t *tension_handler,
1353 one_dim_fn_t *inverse_tension_handler,
1355 void *tension_params,
1356 destroy_data_func_t *destroy_tension,
1359 state_t *ret = (state_t *)malloc(sizeof(state_t));
1360 assert(ret != NULL);
1361 /* make a copy of the name, so we aren't dependent on the original */
1362 ret->name = (char *)malloc(sizeof(char)*(strlen(name)+1));
1363 assert(ret->name != NULL);
1364 strcpy(ret->name, name); /* we know ret->name is long enough */
1366 ret->tension_handler = tension_handler;
1367 ret->inverse_tension_handler = inverse_tension_handler;
1368 ret->tension_group = tension_group;
1369 ret->tension_params = tension_params;
1370 ret->destroy_tension = destroy_tension;
1371 ret->num_domains = num_domains;
1373 ret->num_out_trans = 0;
1374 ret->out_trans = NULL;
1375 ret->num_grouped_states = 0;
1376 ret->grouped_states = NULL;
1379 fprintf(stderr, "generated state %s (%p) with tension handler %p, inverse handler %p, params %p, group %d.\n", ret->name, ret, ret->tension_handler, ret->inverse_tension_handler, ret->tension_params, ret->tension_group);
1384 void destroy_state(state_t *state)
1388 fprintf(stderr, "freeing state %s (%p)\n", state->name, state);
1392 if (state->destroy_tension)
1393 (*state->destroy_tension)(state->tension_params);
1394 if (state->out_trans)
1395 free(state->out_trans);
1396 if (state->grouped_states)
1397 free(state->grouped_states);
1404 We'll be storing the states in a list (see Appendix \ref{sec.lists}),
1405 so we also define a convenience function for cleaning up.
1406 <<destroy state list>>=
1407 void destroy_state_list(list_t *state_list)
1409 state_list = head(state_list);
1410 while (state_list != NULL)
1411 destroy_state((state_t *) pop(&state_list));
1416 We can also define a convenient way to get a state index from it's name.
1418 int state_index(list_t *state_list, char *name)
1421 state_list = head(state_list);
1422 while (state_list != NULL) {
1423 if (STRMATCH(S(state_list)->name, name)) return i;
1424 state_list = state_list->next;
1432 \subsection{Transition data}
1433 \label{sec.transition_data}
1435 Once you've created states (Sections \ref{sec.main},
1436 \ref{sec.model_selection}, and \ref{sec.state_data}), you need to
1437 create transitions between them.
1438 <<transition definitions>>=
1439 typedef struct transition_struct {
1440 state_t *initial_state;
1441 state_t *final_state;
1442 k_func_t *k; /* transition rate model */
1443 void *k_params; /* pointer to k parameters */
1444 destroy_data_func_t *destroy_k;
1447 /* get the transition data for the current list node */
1448 #define T(list) ((transition_t *)(list)->d)
1449 @ [[k]] is a pointer to the function determining the transition rate
1450 for a given tension. [[k_params]] and [[destroy_k]] have similar
1451 roles with regards to [[k]] as [[tension_params]] and
1452 [[destroy_tension]] have with regards to [[tension_handler]] (see
1453 Section \ref{sec.state_data}).
1455 [[create_]] and [[destroy_transition]] are simple wrappers around
1456 [[malloc]] and [[free]].
1457 <<create/destroy transition>>=
1458 transition_t *create_transition(state_t *initial_state, state_t *final_state,
1461 destroy_data_func_t *destroy_k)
1463 transition_t *ret = (transition_t *)malloc(sizeof(transition_t));
1464 assert(ret != NULL);
1465 assert(initial_state != NULL);
1466 assert(final_state != NULL);
1467 ret->initial_state = initial_state;
1468 ret->final_state = final_state;
1470 ret->k_params = k_params;
1471 ret->destroy_k = destroy_k;
1473 fprintf(stderr, "generated transition %p from %s to %s, params %p\n", ret, ret->initial_state->name, ret->final_state->name, ret->k_params);
1478 void destroy_transition(transition_t *transition)
1482 fprintf(stderr, "freeing transition %p (%s to %s)\n", transition, transition->initial_state->name, transition->final_state->name);
1484 if (transition->destroy_k)
1485 (*transition->destroy_k)(transition->k_params);
1492 We'll be storing the transitions in a list (see Appendix
1493 \ref{sec.lists}), so we also define a convenience function for
1495 <<destroy transition list>>=
1496 void destroy_transition_list(list_t *transition_list)
1498 transition_list = head(transition_list);
1499 while (transition_list != NULL)
1500 destroy_transition((transition_t *) pop(&transition_list));
1505 \subsection{Graphviz output}
1506 \label{sec.graphviz_output}
1508 <<print state graph>>=
1509 void print_state_graph(FILE *file, list_t *state_list, list_t *transition_list)
1511 state_list = head(state_list);
1512 transition_list = head(transition_list);
1513 fprintf(file, "digraph states {\n node [shape=record,color=black,fillcolor=white,style=filled];\n\n");
1515 fprintf(file, " n%d [label=\"%s\"]\n", state_index(state_list, S(state_list)->name), S(state_list)->name);
1516 if (state_list->next == NULL) break;
1517 state_list = state_list->next;
1519 fprintf(file, "\n");
1520 while (transition_list != NULL) {
1521 fprintf(file, " n%d -> n%d;\n", state_index(state_list, T(transition_list)->initial_state->name), state_index(state_list, T(transition_list)->final_state->name));
1522 transition_list = transition_list->next;
1524 fprintf(file, "}\n");
1528 \subsection{Domain model and group handling}
1530 <<group functions>>=
1531 <<crosslink groups>>
1532 <<get active groups>>
1535 Scan through a list of states and construct an array of all the
1537 <<get active groups>>=
1538 void get_active_groups(list_t *state_list,
1539 int *pNum_active_groups,
1540 one_dim_fn_t ***pPer_group_handlers,
1541 one_dim_fn_t ***pPer_group_inverse_handlers,
1542 void ***pPer_group_data /* pointer to array of pointers to tension_handler_data_t */)
1544 void **active_groups=NULL;
1545 one_dim_fn_t *handler, *old_handler, *inverse_handler, *old_inverse_handler;
1546 one_dim_fn_t **per_group_handlers=NULL, **per_group_inverse_handlers=NULL;
1547 void **per_group_data=NULL; /* array of pointers to tension_handler_data_t */
1548 int i, j, num_domains, group, old_group, num_active_groups=0;
1549 list_t *active_states_list=NULL;
1550 tension_handler_data_t *tdata=NULL;
1553 /* get all the active groups in a list */
1554 state_list = head(state_list);
1556 fprintf(stderr, "%s:\t", __FUNCTION__);
1557 list_print(stderr, state_list, "states");
1559 while (state_list != NULL) {
1560 state = S(state_list);
1561 handler = state->tension_handler;
1562 inverse_handler = state->inverse_tension_handler;
1563 group = state->tension_group;
1564 num_domains = state->num_domains;
1565 if (list_index(active_states_list, state) == -1) {
1566 /* we haven't added this state (or it's associated group) yet */
1567 if (num_domains > 0 && handler != NULL) { /* add the state */
1568 num_active_groups++;
1569 push(&active_states_list, state, 1);
1571 fprintf(stderr, "%s:\tactive %s state (%p) with %d domain(s)\n", __FUNCTION__, state->name, state, state->num_domains);
1573 for (i=0; i < state->num_grouped_states; i++) {
1574 if (state->grouped_states[i]->num_domains > 0) {
1575 /* add this related (and active) state now */
1576 assert(state->grouped_states[i]->tension_handler == handler);
1577 assert(state->grouped_states[i]->tension_group == group);
1578 push(&active_states_list, state->grouped_states[i], 1);
1580 fprintf(stderr, "%s:\tactive grouped %s state (%p) with %d domain(s)\n", __FUNCTION__, state->grouped_states[i]->name, state->grouped_states[i], state->grouped_states[i]->num_domains);
1584 } /* else empty state or NULL handler */
1585 } /* else state already added */
1586 state_list = state_list->next;
1589 fprintf(stderr, "%s:\t", __FUNCTION__);
1590 list_print(stderr, active_states_list, "active states");
1593 assert(num_active_groups <= list_length(active_states_list));
1594 assert(num_active_groups > 0);
1596 /* allocate the arrays we'll be returning */
1597 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1598 assert(per_group_handlers != NULL);
1599 per_group_inverse_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1600 assert(per_group_inverse_handlers != NULL);
1601 per_group_data = (void **)malloc(sizeof(void *)*num_active_groups);
1602 assert(per_group_data != NULL);
1604 fprintf(stderr, "%s:\tallocating per group data %p:", __FUNCTION__, per_group_data);
1606 for (i=0; i<num_active_groups; i++) { /* we might want to recycle any group data structs allocated in the previous round */
1607 per_group_data[i] = malloc(sizeof(tension_handler_data_t));
1608 assert(per_group_data[i] != NULL);
1610 fprintf(stderr, " %d,%p", i, per_group_data[i]);
1614 fprintf(stderr, "\n");
1617 /* populate the arrays we'll be returning */
1618 active_states_list = head(active_states_list);
1620 old_inverse_handler = NULL;
1621 j = -1; /* j is the current active _group_ index */
1622 while (active_states_list != NULL) {
1623 state = (state_t *) pop(&active_states_list);
1624 handler = state->tension_handler;
1625 inverse_handler = state->inverse_tension_handler;
1626 group = state->tension_group;
1627 num_domains = state->num_domains;
1628 if (handler != old_handler || inverse_handler != old_inverse_handler || group != old_group) {
1629 j++; /* move to the next group */
1630 tdata = (tension_handler_data_t *) per_group_data[j];
1631 per_group_handlers[j] = handler;
1632 per_group_inverse_handlers[j] = inverse_handler;
1633 tdata->group_tension_params = NULL;
1635 tdata->persist = NULL; /* we might want to grab any appropriate data from the old persist here... */
1638 fprintf(stderr, "%s:\tadding active state %s to group %d (%p:%p:%d) with %d domain(s)\n", __FUNCTION__, state->name, j, handler, inverse_handler, group, num_domains);
1640 for (i=0; i<num_domains; i++) { /* add num_domains copies of the tension params */
1641 push(&tdata->group_tension_params, state->tension_params, 1);
1644 fprintf(stderr, "%s:\tpushed %d copies of %p onto data %p's param list %p\n", __FUNCTION__, num_domains, state->tension_params, tdata, tdata->group_tension_params);
1646 old_handler = handler;
1647 old_inverse_handler = inverse_handler;
1651 /* free old groups */
1652 if (*pPer_group_handlers != NULL)
1653 free(*pPer_group_handlers);
1654 if (*pPer_group_inverse_handlers != NULL)
1655 free(*pPer_group_inverse_handlers);
1656 if (*pPer_group_data != NULL) {
1657 for (i=0; i<*pNum_active_groups; i++)
1658 free((*pPer_group_data)[i]);
1659 free(*pPer_group_data);
1662 /* set new groups */
1664 fprintf(stderr, "%s:\t%d groups, tension handler array %p, inverse tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_inverse_handlers, per_group_data, per_group_data[0]);
1666 *pNum_active_groups = num_active_groups;
1667 *pPer_group_handlers = per_group_handlers;
1668 *pPer_group_inverse_handlers = per_group_inverse_handlers;
1669 *pPer_group_data = per_group_data;
1673 <<crosslink groups>>=
1674 void crosslink_groups(list_t *state_groups, list_t *transition_list)
1676 list_t *list, *out_trans=NULL, *associates=NULL;
1677 one_dim_fn_t *handler;
1680 state_groups = head(state_groups);
1681 while (state_groups != NULL) {
1682 /* find transitions leaving this state */
1684 fprintf(stderr, "finding transitions leaving %s\n", S(state_groups)->name);
1686 list = head(transition_list);
1687 while (list != NULL) {
1688 if (T(list)->initial_state == S(state_groups) && list_index(associates, T(list)->initial_state) != -1) {
1689 push(&out_trans, T(list), 1);
1693 n = list_length(out_trans);
1694 S(state_groups)->num_out_trans = n;
1695 S(state_groups)->out_trans = (transition_t **)malloc(sizeof(transition_t *)*n);
1696 assert(S(state_groups)->out_trans != NULL);
1697 for (i=0; i<n; i++) {
1698 S(state_groups)->out_trans[i] = (transition_t *)pop(&out_trans);
1701 /* find states grouped with this state */
1703 fprintf(stderr, "finding states grouped with %s\n", S(state_groups)->name);
1705 handler = S(state_groups)->tension_handler;
1706 group = S(state_groups)->tension_group;
1707 list = head(state_groups);
1708 while (list != NULL) {
1709 if (S(list)->tension_handler == handler && S(list)->tension_group == group && list != state_groups) {
1710 push(&associates, S(list), 1);
1714 n = list_length(associates);
1715 S(state_groups)->num_grouped_states = n;
1716 S(state_groups)->grouped_states = (state_t **)malloc(sizeof(state_t *)*n);
1717 assert(S(state_groups)->grouped_states != NULL);
1718 for (i=0; i<n; i++) {
1719 S(state_groups)->grouped_states[i] = (state_t *)pop(&associates);
1721 state_groups = state_groups->next;
1727 \section{String parsing}
1729 For handling command line arguments, we need parse delimited strings (e.g.\ [["param1,param2,param3..."]]).
1730 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1731 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1732 need to take care of parsing those parameters themselves.
1733 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1739 <<parse definitions>>
1740 <<parse declarations>>
1744 <<parse module makefile lines>>=
1745 NW_SAWSIM_MODS += parse
1746 CHECK_BINS += check_parse
1750 <<parse definitions>>=
1751 #define SEP ',' /* argument separator character */
1754 <<parse declarations>>=
1755 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1756 int *num, char ***string_array);
1757 extern void free_parsed_list(int num, char **string_array);
1760 [[parse_list_string]] allocates memory, don't forget to free it
1761 afterward with [[free_parsed_list]]. It does not alter the original.
1763 The string may start off with a [[deeper]] character
1764 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1765 [[x,y]], where the pointer is one character in on the copied string.
1766 However, when we go to free the memory, we need a pointer to the
1767 beginning of the string. In order to accommodate this for a string
1768 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1769 the first $N$ elements point to the separated arguments, and let the
1770 last element point to the start of the copied string regardless of
1772 <<parse delimited list string functions>>=
1773 /* TODO, split out into parse.hc */
1774 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1777 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1778 if (string[i] == deeper) {depth++;}
1779 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1785 void parse_list_string(char *string, char sep, char deeper, char shallower,
1786 int *num, char ***string_array)
1788 char *str=NULL, **ret=NULL;
1790 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
1792 *string_array = NULL;
1795 /* make a copy of the string, so we don't change the original */
1796 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1797 assert(str != NULL);
1798 strcpy(str, string); /* we know str is long enough */
1799 /* count the number of regions, so we can allocate pointers to them */
1802 n++; i++; /* move on to next argument */
1803 i += next_delim_index(str+i, sep, deeper, shallower);
1804 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1805 } while (str[i] != '\0');
1806 ret = (char **)malloc(sizeof(char *)*(n+1));
1807 assert(ret != NULL);
1808 /* replace the separators with '\0' & assign pointers */
1809 ret[n] = str; /* point to the front of the copied string */
1812 for(i=1; i<n; i++) {
1813 j += next_delim_index(str+j, sep, deeper, shallower);
1815 ret[i] = str+j; /* point to the separated arguments */
1817 /* strip off bounding braces, if any */
1818 for(i=0; i<n; i++) {
1819 if (ret[i][0] == deeper) {
1823 if (ret[i][strlen(ret[i])-1] == shallower)
1824 ret[i][strlen(ret[i])-1] = '\0';
1827 *string_array = ret;
1830 void free_parsed_list(int num, char **string_array)
1833 /* frees all items, since they were allocated as a single string*/
1834 free(string_array[num]);
1835 /* and free the array of pointers */
1841 \subsection{Parsing implementation}
1847 <<parse delimited list string functions>>
1851 #include <assert.h> /* assert() */
1852 #include <stdlib.h> /* NULL */
1853 #include <stdio.h> /* fprintf(), stdout *//*!!*/
1854 #include <string.h> /* strlen() */
1858 \subsection{Parsing unit tests}
1860 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1862 <<parse check includes>>
1863 <<string comparison definition and globals>>
1864 <<check relative error>>
1865 <<parse test suite>>
1866 <<main check program>>
1869 <<parse check includes>>=
1870 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
1871 #include <stdio.h> /* printf() */
1872 #include <assert.h> /* assert() */
1873 #include <string.h> /* strlen() */
1878 <<parse test suite>>=
1879 <<parse list string tests>>
1880 <<parse suite definition>>
1883 <<parse suite definition>>=
1884 Suite *test_suite (void)
1886 Suite *s = suite_create ("k model");
1887 <<parse list string test case defs>>
1889 <<parse list string test case add>>
1894 <<parse list string tests>>=
1897 START_TEST(test_next_delim_index)
1899 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
1900 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
1901 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
1902 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
1903 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
1907 START_TEST(test_parse_list_null)
1911 parse_list_string(NULL, SEP, '{', '}',
1912 &num_param_args, ¶m_args);
1913 fail_unless(num_param_args == 0, NULL);
1914 fail_unless(param_args == NULL, NULL);
1917 START_TEST(test_parse_list_single_simple)
1921 parse_list_string("arg", SEP, '{', '}',
1922 &num_param_args, ¶m_args);
1923 fail_unless(num_param_args == 1, NULL);
1924 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
1927 START_TEST(test_parse_list_single_compound)
1931 parse_list_string("{x,y,z}", SEP, '{', '}',
1932 &num_param_args, ¶m_args);
1933 fail_unless(num_param_args == 1, NULL);
1934 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
1937 START_TEST(test_parse_list_double_simple)
1941 parse_list_string("abc,def", SEP, '{', '}',
1942 &num_param_args, ¶m_args);
1943 fail_unless(num_param_args == 2, NULL);
1944 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
1945 fail_unless(STRMATCH(param_args[1],"def"), NULL);
1948 START_TEST(test_parse_list_compound)
1952 parse_list_string("abc,{def,ghi}", SEP, '{', '}',
1953 &num_param_args, ¶m_args);
1954 fail_unless(num_param_args == 2, NULL);
1955 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
1956 fail_unless(STRMATCH(param_args[1],"def,ghi"), NULL);
1960 <<parse list string test case defs>>=
1961 TCase *tc_parse_list_string = tcase_create("parse list string");
1963 <<parse list string test case add>>=
1964 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
1965 tcase_add_test(tc_parse_list_string, test_parse_list_null);
1966 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
1967 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
1968 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
1969 tcase_add_test(tc_parse_list_string, test_parse_list_compound);
1970 suite_add_tcase(s, tc_parse_list_string);
1974 \section{Unit tests}
1976 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1983 <<check relative error>>
1986 <<main check program>>
1998 <<determine dt tests>>
2000 <<does domain transition tests>>
2001 <<randomly unfold domains tests>>
2002 <<suite definition>>
2005 <<suite definition>>=
2006 Suite *test_suite (void)
2008 Suite *s = suite_create ("sawsim");
2009 <<F test case defs>>
2010 <<determine dt test case defs>>
2011 <<happens test case defs>>
2012 <<does domain transition test case defs>>
2013 <<randomly unfold domains test case defs>>
2016 <<determine dt test case add>>
2017 <<happens test case add>>
2018 <<does domain transition test case add>>
2019 <<randomly unfold domains test case add>>
2022 tcase_add_checked_fixture(tc_strip_address,
2023 setup_strip_address,
2024 teardown_strip_address);
2030 <<main check program>>=
2034 Suite *s = test_suite();
2035 SRunner *sr = srunner_create(s);
2036 srunner_run_all(sr, CK_ENV);
2037 number_failed = srunner_ntests_failed(sr);
2039 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
2043 \subsection{F tests}
2045 <<worm-like chain tests>>
2047 <<F test case defs>>=
2048 <<worm-like chain test case def>>
2050 <<F test case add>>=
2051 <<worm-like chain test case add>>
2054 <<worm-like chain tests>>=
2055 <<worm-like chain function>>
2057 START_TEST(test_wlc_at_zero)
2059 double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
2060 fail_unless(abs(wlc(x, T, p, L)) < lim, \
2061 "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
2062 x, T, p, L, abs(wlc(x, T, p, L)), lim);
2066 START_TEST(test_wlc_at_half)
2068 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
2069 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
2070 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
2072 * linear term = x/L = 0.5
2073 * nonlinear + linear = 0.75 + 0.5 = 1.25
2074 * wlc = 10e21*1.25 = 12.5e21
2076 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
2077 "wlc(%g, %g, %g, %g) = %g != %g",
2078 x, T, p, L, wlc(x, T, p, L), 12.5e21);
2082 <<worm-like chain test case def>>=
2083 TCase *tc_wlc = tcase_create("WLC");
2086 <<worm-like chain test case add>>=
2087 tcase_add_test(tc_wlc, test_wlc_at_zero);
2088 tcase_add_test(tc_wlc, test_wlc_at_half);
2089 suite_add_tcase(s, tc_wlc);
2092 \subsection{Model tests}
2094 Check the searching with [[linear_k]].
2095 Check overwhelming force treatment with the heavyside-esque step [[?]].
2097 Hmm, I'm not really sure what I was doing here...
2098 <<determine dt tests>>=
2099 double linear_k(double F, environment_t *env, void *params)
2101 double Fnot = *(double *)params;
2105 START_TEST(test_determine_dt_linear_k)
2108 transition_t unfold={0};
2109 environment_t env = {0};
2110 double F[]={0,1,2,3,4,5,6};
2111 double dt_max=3.0, Fnot=3.0, x=1.0, targ_p=0.1, v=1.0;
2112 list_t *state_list=NULL, *transition_list=NULL;
2115 folded.tension_handler = &hooke_handler;
2116 folded.num_domains = 1;
2117 unfold.initial_state = &folded;
2118 unfold.k = linear_k;
2119 unfold.k_params = &Fnot;
2120 push(&state_list, &folded, 1);
2121 push(&transition_list, &unfold, 1);
2122 for( i=0; i < sizeof(F)/sizeof(double); i++) {
2123 //fail_unless(determine_dt(state_list, transition_list, &env, x, targ_p, dt_max, v)==1, NULL);
2128 <<determine dt test case defs>>=
2129 TCase *tc_determine_dt = tcase_create("Determine dt");
2131 <<determine dt test case add>>=
2132 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
2133 suite_add_tcase(s, tc_determine_dt);
2138 <<happens test case defs>>=
2140 <<happens test case add>>=
2143 <<does domain transition tests>>=
2145 <<does domain transition test case defs>>=
2147 <<does domain transition test case add>>=
2150 <<randomly unfold domains tests>>=
2152 <<randomly unfold domains test case defs>>=
2154 <<randomly unfold domains test case add>>=
2158 \section{Balancing group extensions}
2159 \label{sec.tension_balance}
2161 For a given total extension $x$ (of the piezo), the various domain
2162 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
2163 amounts, and we need to tweak the portion that each extends to balance
2164 the tension amongst the active groups. Since the tension balancing is
2165 separable from the bulk of the simulation, we break it out into a
2166 separate module. The interface is defined in [[tension_balance.h]],
2167 the implementation in [[tension_balance.c]], and the unit testing in
2168 [[check_tension_balance.c]]
2170 <<tension-balance.h>>=
2171 #ifndef TENSION_BALANCE_H
2172 #define TENSION_BALANCE_H
2174 <<tension balance function declaration>>
2175 #endif /* TENSION_BALANCE_H */
2178 <<tension balance functions>>=
2179 <<one dimensional bracket>>
2180 <<one dimensional solve>>
2182 <<group stiffness function>>
2183 <<chain stiffness function>>
2184 <<tension balance function>>
2187 <<tension balance module makefile lines>>=
2188 NW_SAWSIM_MODS += tension_balance
2189 CHECK_BINS += check_tension_balance
2190 check_tension_balance_MODS =
2193 The entire force balancing problem reduces to a solving two nested
2194 one-dimensional root-finding problems. First we define the one
2195 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
2196 populated group). $x(x_0)$ is also strictly monotonically increasing,
2197 so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
2198 stores the last successful value of $x$, since we expect to be taking
2199 small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
2200 indicates that the groups have changed.
2201 <<tension balance function declaration>>=
2202 double tension_balance(int num_tension_groups,
2203 one_dim_fn_t **tension_handlers,
2204 one_dim_fn_t **inverse_tension_handlers,
2205 void **params, /* array of pointers to tension_handler_data_t */
2210 <<tension balance function>>=
2211 double tension_balance(int num_tension_groups,
2212 one_dim_fn_t **tension_handlers,
2213 one_dim_fn_t **inverse_tension_handlers,
2214 void **params, /* array of pointers to tension_handler_data_t */
2219 static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL,NULL};
2220 double F, xo_guess, xo, lb, ub;
2221 double min_relx=1e-6, min_rely=1e-6;
2222 int max_steps=200, i;
2224 assert(num_tension_groups > 0);
2225 assert(tension_handlers != NULL);
2226 assert(params != NULL);
2227 assert(num_tension_groups > 0);
2229 if (last_x == -1 || x_xo_data.xi == NULL) { /* new group setup, reset x_xo_data */
2230 /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
2231 if (x_xo_data.xi != NULL)
2233 /* construct new x_xo_data */
2234 x_xo_data.n_groups = num_tension_groups;
2235 x_xo_data.tension_handler = tension_handlers;
2236 x_xo_data.inverse_tension_handler = inverse_tension_handlers;
2237 x_xo_data.handler_data = params;
2239 fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p), x_of_xo_data %p\n", __FUNCTION__, x_xo_data.n_groups, x_xo_data.tension_handler, x_xo_data.handler_data, x_xo_data.handler_data[0], &x_xo_data);
2241 x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
2242 assert(x_xo_data.xi != NULL);
2243 for (i=0; i<num_tension_groups; i++)
2244 x_xo_data.xi[i] = last_x;
2247 if (num_tension_groups == 1) { /* only one group, no balancing required */
2249 x_xo_data.xi[0] = xo;
2251 fprintf(stderr, "%s:\tonly one active group, no balancing required\n", __FUNCTION__);
2255 fprintf(stderr, "%s:\tbalancing for x = %g with", __FUNCTION__, x);
2256 for(i=0; i<num_tension_groups; i++) fprintf(stderr, "\tx[%d]=%g", i, x_xo_data.xi[i]);
2257 fprintf(stderr, "\n");
2259 if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
2260 if (x == 0) xo_guess = length_scale;
2261 else xo_guess = x/num_tension_groups;
2263 fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
2265 oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
2266 } else { /* work off the last balanced point */
2268 fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
2271 lb = x_xo_data.xi[0];
2272 ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */
2273 } else if (x < last_x) {
2274 lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
2275 ub = x_xo_data.xi[0];
2276 } else { /* x == last_x */
2277 fprintf(stderr, "not moving\n");
2278 lb= x_xo_data.xi[0];
2279 ub= x_xo_data.xi[0];
2283 fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
2284 __FUNCTION__, x, lb, ub);
2286 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
2287 F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
2289 if (num_tension_groups > 1) {
2290 fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
2291 for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
2292 fprintf(stderr, "\n");
2297 F = (*tension_handlers[0])(xo, params[0]);
2299 if (stiffness != NULL)
2300 *stiffness = chain_stiffness(&x_xo_data, min_relx);
2307 <<tension balance internal definitions>>=
2308 <<x of x0 definitions>>
2311 <<x of x0 definitions>>=
2312 typedef struct x_of_xo_data_struct {
2314 one_dim_fn_t **tension_handler; /* array of fn pointers */
2315 one_dim_fn_t **inverse_tension_handler; /* inverse fn pointers */
2316 void **handler_data; /* array of pointers to tension_handler_data_t */
2317 double *xi; /* array of group extensions */
2322 double x_of_xo(double xo, void *pdata)
2324 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
2325 double F, x=0, xi, xi_guess, lb, ub, slope;
2327 double min_relx=1e-6, min_rely=1e-6;
2329 assert(data->n_groups > 0);
2332 fprintf(stderr, "%s:\tget initial tension at %g with %d groups, handler %p, data %p (x_xo_data %p)\n", __FUNCTION__, xo, data->n_groups, data->tension_handler[0], data->handler_data[0], data);
2335 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
2337 fprintf(stderr, "%s:\tfinding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
2341 if (data->xi) data->xi[0] = xo;
2342 for (i=1; i < data->n_groups; i++) {
2343 if (data->inverse_tension_handler[i] != NULL) {
2344 xi = (*data->inverse_tension_handler[i])(F, data->handler_data[i]);
2345 } else { /* invert numerically */
2346 if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
2347 else xi_guess = data->xi[i];
2349 fprintf(stderr, "%s:\tsearching for proper x[%d] for F[%d](x[%d]) = %g N (handler %p, data %p)\n", __FUNCTION__, i, i, i, F, data->tension_handler[i], data->handler_data[i]);
2351 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
2352 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
2357 fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
2361 fprintf(stderr, "%s:\tfound x(x0=%g) = %g\n", __FUNCTION__, xo, x);
2367 <<group stiffness function>>=
2368 double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double x, double relx)
2370 double dx, xl, F, Fl, stiffness;
2372 assert(relx > 0 && relx < 1);
2373 if (x == 0 || relx == 0) {
2374 dx = 1e-12; /* HACK, how to get length scale? */
2384 F = (*tension_handler)(x, handler_data);
2385 Fl = (*tension_handler)(xl, handler_data);
2386 stiffness = (F-Fl)/dx;
2387 assert(stiffness > 0);
2392 <<chain stiffness function>>=
2393 double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx)
2397 /* add group stiffnesses in series */
2398 for (i=0; i < x_xo_data->n_groups; i++) {
2400 fprintf(stderr, "%s:\tgetting stiffness of active state %d of %d for x[%d]=%g, relx=%g\n", __FUNCTION__, i, x_xo_data->n_groups, i, x_xo_data->xi[i], relx);
2403 stiffness += 1.0/group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx);
2405 stiffness = 1.0 / stiffness;
2411 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
2412 number of steps for monotonically increasing $f(x)$. Simple
2413 bisection, so it's robust and converges fairly quickly. You can set
2414 the maximum number of steps to take, but if you have enough processor
2415 power, it's better to limit your precision with the relative limits
2416 [[min_relx]] and [[y]]. These ensure that the width of the bracket is
2417 small on the length scale of it's larger side. Note that \emph{both}
2418 relative limits must be satisfied for the search to stop.
2419 <<one dimensional function definition>>=
2420 /* equivalent to gsl_func_t */
2421 typedef double one_dim_fn_t(double x, void *params);
2423 <<one dimensional solve declaration>>=
2424 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2425 double min_relx, double min_rely, int max_steps,
2429 <<one dimensional solve>>=
2430 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2431 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2432 double min_relx, double min_rely, int max_steps,
2435 double yx, ylb, yub, x;
2438 ylb = (*f)(lb, params);
2439 yub = (*f)(ub, params);
2441 /* check some simple cases */
2443 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bracketed */
2444 else return lb; /* any x on the range [lb, ub] would work */
2446 if (ylb == y) { x=lb; yx=ylb; goto end; }
2447 if (yub == y) { x=ub; yx=yub; goto end; }
2450 fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2452 assert(ylb < y); assert(yub > y);
2453 for (i=0; i<max_steps; i++) {
2455 yx = (*f)(x, params);
2457 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);
2459 if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2460 else if (yx > y) { ub=x; yub=yx; }
2461 else /* yx < y */ { lb=x; ylb=yx; }
2466 fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2472 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2473 Generate bracketing $x$ values through bisection or exponential growth.
2474 <<one dimensional bracket declaration>>=
2475 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2478 <<one dimensional bracket>>=
2479 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2481 double yx, step, x=xguess;
2484 fprintf(stderr, "%s:\tbracketing %g (params %p)\n", __FUNCTION__, y, params);
2486 yx = (*f)(x, params);
2488 fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2495 fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2499 if (x == 0) x = length_scale/2.0; /* HACK */
2502 yx = (*f)(x, params);
2504 fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2509 fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2512 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2516 \subsection{Balancing implementation}
2518 <<tension-balance.c>>=
2521 <<tension balance includes>>
2522 #include "tension_balance.h"
2524 double length_scale = 1e-10; /* HACK */
2526 <<tension balance internal definitions>>
2527 <<tension balance functions>>
2530 <<tension balance includes>>=
2531 #include <assert.h> /* assert() */
2532 #include <stdlib.h> /* NULL */
2533 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2534 #include <stdio.h> /* fprintf(), stdout */
2538 \subsection{Balancing unit tests}
2540 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2541 <<check-tension-balance.c>>=
2542 <<tension balance check includes>>
2543 <<tension balance test suite>>
2544 <<main check program>>
2547 <<tension balance check includes>>=
2548 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2549 #include <assert.h> /* assert() */
2552 #include "tension_balance.h"
2555 <<tension balance test suite>>=
2556 <<tension balance function tests>>
2557 <<tension balance suite definition>>
2560 <<tension balance suite definition>>=
2561 Suite *test_suite (void)
2563 Suite *s = suite_create ("tension balance");
2564 <<tension balance function test case defs>>
2566 <<tension balance function test case adds>>
2571 <<tension balance function tests>>=
2572 <<check relative error>>
2574 double hooke(double x, void *pK)
2577 return *((double*)pK) * x;
2580 START_TEST(test_single_function)
2582 double k=5, x=3, last_x=2.0, F;
2583 one_dim_fn_t *handlers[] = {&hooke};
2584 one_dim_fn_t *inverse_handlers[] = {NULL};
2585 void *data[] = {&k};
2586 F = tension_balance(1, handlers, inverse_handlers, data, last_x, x, NULL);
2587 fail_unless(F = k*x, NULL);
2592 We can also test balancing two springs with different spring constants.
2593 The analytic solution for a general number of springs is given in Appendix \ref{sec.math_hooke}.
2594 <<tension balance function tests>>=
2595 START_TEST(test_double_hooke)
2597 double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
2598 one_dim_fn_t *handlers[] = {&hooke, &hooke};
2599 one_dim_fn_t *inverse_handlers[] = {NULL, NULL};
2600 void *data[] = {&k1, &k2};
2601 F = tension_balance(2, handlers, inverse_handlers, data, last_x, x, NULL);
2605 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2606 //CHECK_ERR(1e-6, x1e, xi[0]);
2607 //CHECK_ERR(1e-6, x2e, xi[1]);
2608 CHECK_ERR(1e-6, Fe, F);
2612 TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
2614 <<tension balance function test case defs>>=
2615 TCase *tc_tbfunc = tcase_create("tension balance function");
2618 <<tension balance function test case adds>>=
2619 tcase_add_test(tc_tbfunc, test_single_function);
2620 tcase_add_test(tc_tbfunc, test_double_hooke);
2621 suite_add_tcase(s, tc_tbfunc);
2627 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2628 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2629 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2635 <<list definitions>>
2636 <<list declarations>>
2640 <<list declarations>>=
2641 <<head and tail declarations>>
2642 <<list length declaration>>
2643 <<push declaration>>
2645 <<list destroy declaration>>
2646 <<list to array declaration>>
2647 <<move declaration>>
2648 <<sort declaration>>
2649 <<list index declaration>>
2650 <<list element declaration>>
2651 <<unique declaration>>
2652 <<list print declaration>>
2656 <<create and destroy node>>
2671 <<list module makefile lines>>=
2672 NW_SAWSIM_MODS += list
2673 CHECK_BINS += check_list
2677 <<list definitions>>=
2678 typedef struct list_struct {
2679 struct list_struct *next;
2680 struct list_struct *prev;
2684 <<comparison function typedef>>
2687 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2688 <<head and tail declarations>>=
2689 list_t *head(list_t *list);
2690 list_t *tail(list_t *list);
2693 list_t *head(list_t *list)
2697 while (list->prev != NULL) {
2703 list_t *tail(list_t *list)
2707 while (list->next != NULL) {
2714 <<list length declaration>>=
2715 int list_length(list_t *list);
2718 int list_length(list_t *list)
2725 while (list->next != NULL) {
2733 [[push]] inserts a node relative to the current position in the list
2734 without changing the current position.
2735 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2736 so we set it to that of the pushed domain.
2737 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2738 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2739 <<push declaration>>=
2740 void push(list_t **pList, void *data, int below);
2743 void push(list_t **pList, void *data, int below)
2745 list_t *list, *new_node;
2746 assert(pList != NULL);
2748 new_node = create_node(data);
2751 else if (below==0) { /* insert above */
2752 if (list->prev != NULL)
2753 list->prev->next = new_node;
2754 new_node->prev = list->prev;
2755 list->prev = new_node;
2756 new_node->next = list;
2757 } else { /* insert below */
2758 if (list->next != NULL)
2759 list->next->prev = new_node;
2760 new_node->next = list->next;
2761 list->next = new_node;
2762 new_node->prev = list;
2767 [[pop]] removes the current domain node, moving the current position
2768 to the node after it, unless that node is [[NULL]], in which case move
2769 the current position to the node before it.
2770 <<pop declaration>>=
2771 void *pop(list_t **pList);
2774 void *pop(list_t **pList)
2776 list_t *list, *popped;
2778 assert(pList != NULL);
2780 assert(list != NULL); /* not an empty list */
2782 /* bypass the popped node */
2783 if (list->prev != NULL)
2784 list->prev->next = list->next;
2785 if (list->next != NULL) {
2786 list->next->prev = list->prev;
2787 *pList = list->next;
2789 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
2791 destroy_node(popped);
2796 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
2797 If the cleanup function is [[NULL]], no cleanup function is called.
2798 The nodes are not popped in any particular order.
2799 <<list destroy declaration>>=
2800 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
2803 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
2807 assert(pList != NULL);
2810 assert(list != NULL); /* not an empty list */
2811 while (list != NULL) {
2813 if (destroy != NULL)
2819 [[list_to_array]] converts a list to an array of pointers, preserving order.
2820 <<list to array declaration>>=
2821 void list_to_array(list_t *list, int *length, void ***array);
2824 void list_to_array(list_t *list, int *pLength, void ***pArray)
2828 assert(list != NULL);
2829 assert(pLength != NULL);
2830 assert(pArray != NULL);
2832 length = list_length(list);
2833 array = (void **)malloc(sizeof(void *)*length);
2834 assert(array != NULL);
2837 while (list != NULL) {
2838 array[i++] = list->d;
2846 [[move]] moves the current element along the list in either direction.
2847 <<move declaration>>=
2848 void move(list_t **pList, int downstream);
2851 void move(list_t **pList, int downstream)
2853 list_t *list, *flipper;
2855 assert(pList != NULL);
2856 list = *pList; /* a>B>c>d */
2857 assert(list != NULL); /* not an empty list */
2858 if (downstream == 0)
2859 flipper = list->prev; /* flipper is A */
2861 flipper = list->next; /* flipper is C */
2862 /* ensure there is somewhere to go in the direction we're heading */
2863 assert(flipper != NULL);
2864 /* remove the current element */
2865 data = pop(&list); /* data=B, a>C>d */
2866 /* and put it back in in the correct spot */
2868 if (downstream == 0) {
2869 push(&list, data, 0); /* b>A>c>d */
2870 list = list->prev; /* B>a>c>d */
2872 push(&list, data, 1); /* a>C>b>d */
2873 list = list->next; /* a>c>B>d */
2879 [[sort]] sorts a list from smallest to largest according to a user
2880 specified comparison.
2881 <<comparison function typedef>>=
2882 typedef int (comparison_fn_t)(void *A, void *B);
2885 <<int comparison function>>=
2886 int int_comparison(void *A, void *B)
2891 if (a > b) return 1;
2892 else if (a == b) return 0;
2897 <<sort declaration>>=
2898 void sort(list_t **pList, comparison_fn_t *comp);
2900 Here we use the bubble sort algorithm.
2902 void sort(list_t **pList, comparison_fn_t *comp)
2906 assert(pList != NULL);
2908 assert(list != NULL); /* not an empty list */
2909 while (stable == 0) {
2912 while (list->next != NULL) {
2913 if ((*comp)(list->d, list->next->d) > 0) {
2915 move(&list, 1 /* downstream */);
2925 [[list_index]] finds the location of [[data]] in [[list]]. Returns
2926 [[-1]] if [[data]] is not in [[list]].
2927 <<list index declaration>>=
2928 int list_index(list_t *list, void *data);
2932 int list_index(list_t *list, void *data)
2936 while (list != NULL) {
2937 if (list->d == data) return i;
2946 [[list_element]] returns the data bound to the $i^\text{th}$ node of the list.
2947 <<list element declaration>>=
2948 void *list_element(list_t *list, int i);
2951 void *list_element(list_t *list, int i)
2955 while (list != NULL) {
2956 if (i == j) return list->d;
2966 [[unique]] removes duplicates from a sorted list according to a user
2967 specified comparison. Don't do this unless you have other pointers
2968 any dynamically allocated data in your list, because [[unique]] just
2969 drops any non-unique data without freeing it.
2970 <<unique declaration>>=
2971 void unique(list_t **pList, comparison_fn_t *comp);
2974 void unique(list_t **pList, comparison_fn_t *comp)
2977 assert(pList != NULL);
2979 assert(list != NULL); /* not an empty list */
2981 while (list->next != NULL) {
2982 if ((*comp)(list->d, list->next->d) == 0) {
2991 [[list_print]] prints a list to a [[FILE *]] stream.
2992 <<list print declaration>>=
2993 void list_print(FILE *file, list_t *list, char *list_name);
2996 void list_print(FILE *file, list_t *list, char *list_name)
2998 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
3000 while (list != NULL) {
3001 fprintf(file, " %p", list->d);
3004 fprintf(file, "\n");
3008 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
3009 <<create and destroy node>>=
3010 list_t *create_node(void *data)
3012 list_t *ret = (list_t *)malloc(sizeof(list_t));
3013 assert(ret != NULL);
3020 void destroy_node(list_t *node)
3026 The user must free the data pointed to by the node on their own.
3028 \subsection{List implementation}
3038 #include <assert.h> /* assert() */
3039 #include <stdlib.h> /* malloc(), free(), rand() */
3040 #include <stdio.h> /* fprintf(), stdout, FILE */
3041 #include "global.h" /* destroy_data_func_t */
3044 \subsection{List unit tests}
3046 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3048 <<list check includes>>
3050 <<main check program>>
3053 <<list check includes>>=
3054 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3055 #include <stdio.h> /* FILE */
3061 <<list test suite>>=
3064 <<list suite definition>>
3067 <<list suite definition>>=
3068 Suite *test_suite (void)
3070 Suite *s = suite_create ("list");
3071 <<push test case defs>>
3072 <<pop test case defs>>
3074 <<push test case adds>>
3075 <<pop test case adds>>
3081 START_TEST(test_push)
3086 push(&list, (void *)i, 1);
3087 fail_unless(list == head(list), NULL);
3088 fail_unless( (int)list->d == 0 );
3089 for (i=0; i<n; i++) {
3090 /* we expect 0, n-1, n-2, ..., 1 */
3093 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
3099 <<push test case defs>>=
3100 TCase *tc_push = tcase_create("push");
3103 <<push test case adds>>=
3104 tcase_add_test(tc_push, test_push);
3105 suite_add_tcase(s, tc_push);
3110 <<pop test case defs>>=
3112 <<pop test case adds>>=
3115 \section{Function string evaluation}
3117 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).
3118 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
3119 The solution is to run a scripting language as a subprocess accessed via pipes.
3120 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
3122 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
3123 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
3124 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.
3125 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
3127 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.
3128 Then you can either statically or dynamically link to those hardcoded functions.
3129 While much less flexible, this approach would be a fairly simple-to-implement fallback.
3131 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
3132 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
3135 #ifndef STRING_EVAL_H
3136 #define STRING_EVAL_H
3138 <<string eval setup declaration>>
3139 <<string eval function declaration>>
3140 <<string eval teardown declaration>>
3141 #endif /* STRING_EVAL_H */
3144 <<string eval module makefile lines>>=
3145 NW_SAWSIM_MODS += string_eval
3146 CHECK_BINS += check_string_eval
3147 check_string_eval_MODS =
3150 For an introduction to POSIX process control, see\\
3151 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
3152 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
3153 and of course, the relavent [[man]] pages.
3155 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
3156 [[execvp]] replaces the calling process' program with a new program.
3157 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
3158 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
3159 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
3160 The new program needs command line arguments, just like it would if you were running it from a shell.
3161 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
3162 with the final array entry being a [[NULL]] pointer.
3164 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
3165 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
3166 <<string eval subprocess definitions>>=
3167 #define SUBPROCESS "python"
3168 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
3169 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
3170 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
3172 The [[i]] option lets Python know that it should run in interactive mode.
3173 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
3174 In interactive mode, python acts on every instruction as soon as it is recieved.
3175 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.
3176 %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.
3178 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
3179 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
3180 [[fork]] returns two copies of the same program, executing the original code.
3181 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
3182 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
3184 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.
3185 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
3186 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
3187 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
3188 <<string eval pipe definitions>>=
3189 #define PIPE_READ 0 /* the end you read from */
3190 #define PIPE_WRITE 1 /* the end you write to */
3192 #define STDIN 0 /* initial index of stdin pair */
3193 #define STDOUT 2 /* initial index of stdout pair */
3196 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
3198 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.
3200 <<string eval setup declaration>>=
3201 extern void string_eval_setup(FILE **pIn, FILE **pOut);
3203 <<string eval setup definition>>=
3204 void string_eval_setup(FILE **pIn, FILE **pOut)
3207 int pfd[4]; /* file descriptors for stdin and stdout */
3209 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
3210 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
3212 pid = fork(); /* split process into two copies */
3214 if (pid == -1) { /* fork error */
3215 perror("fork error");
3217 } else if (pid == 0) { /* child */
3218 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
3219 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
3220 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
3221 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
3222 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
3223 perror("exec error"); /* exec shouldn't return */
3225 } else { /* parent */
3226 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
3227 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
3228 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
3229 if ( *pIn == NULL ) {
3230 perror("fdopen (in)");
3233 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
3234 if ( *pOut == NULL ) {
3235 perror("fdopen (out)");
3242 To use the evaluating subprocess, we just pipe in our command, and read out the results.
3243 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
3244 <<string eval function declaration>>=
3245 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
3247 <<string eval function definition>>=
3248 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
3251 rc = fprintf(in, "%s", input);
3252 assert(rc == strlen(input));
3255 alarm(1); /* set a one second timeout on the read */
3256 assert( fgets(output, buflen, out) != NULL );
3257 alarm(0); /* cancel the timeout */
3258 //fprintf(stderr, "eval: %s ----> %s", input, output);
3261 The [[alarm]] calls set and clear a timeout on the returned output.
3262 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.
3263 This protects against invalid input for which a line of output is not printed to [[stdout]].
3264 Other invalid inputs (e.g.\ those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
3265 If you are getting strange results, check your python code seperately. TODO, better error handling.
3267 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
3268 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
3269 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
3270 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
3271 <<string eval teardown declaration>>=
3272 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
3275 <<string eval teardown definition>>=
3276 void string_eval_teardown(FILE **pIn, FILE **pOut)
3281 /* redirect Python's stderr */
3282 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
3286 assert( fclose(*pIn) == 0 );
3288 assert( fclose(*pOut) == 0 );
3291 /* wait for python to exit */
3293 pid = wait(&stat_loc);
3300 if (WIFEXITED(stat_loc)) {
3301 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
3302 } else if (WIFSIGNALED(stat_loc)) {
3303 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
3308 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
3310 \subsection{String evaluation implementation}
3314 <<string eval includes>>
3315 #include "string_eval.h"
3316 <<string eval internal definitions>>
3317 <<string eval functions>>
3320 <<string eval includes>>=
3321 #include <assert.h> /* assert() */
3322 #include <stdlib.h> /* NULL */
3323 #include <stdio.h> /* fprintf(), stdout, fdopen() */
3324 #include <string.h> /* strlen() */
3325 #include <sys/types.h> /* pid_t */
3326 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
3327 #include <wait.h> /* wait() */
3330 <<string eval internal definitions>>=
3331 <<string eval subprocess definitions>>
3332 <<string eval pipe definitions>>
3335 <<string eval functions>>=
3336 <<string eval setup definition>>
3337 <<string eval function definition>>
3338 <<string eval teardown definition>>
3341 \subsection{String evaluation unit tests}
3343 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3344 <<check-string-eval.c>>=
3345 <<string eval check includes>>
3346 <<string comparison definition and globals>>
3347 <<string eval test suite>>
3348 <<main check program>>
3351 <<string eval check includes>>=
3352 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
3353 #include <stdio.h> /* printf() */
3354 #include <assert.h> /* assert() */
3355 #include <string.h> /* strlen() */
3356 #include <signal.h> /* SIGKILL */
3358 #include "string_eval.h"
3361 <<string eval test suite>>=
3362 <<string eval tests>>
3363 <<string eval suite definition>>
3366 <<string eval suite definition>>=
3367 Suite *test_suite (void)
3369 Suite *s = suite_create ("string eval");
3370 <<string eval test case defs>>
3372 <<string eval test case add>>
3377 <<string eval tests>>=
3378 START_TEST(test_setup_teardown)
3381 string_eval_setup(&in, &out);
3382 string_eval_teardown(&in, &out);
3385 START_TEST(test_invalid_command)
3388 char input[80], output[80]={};
3389 string_eval_setup(&in, &out);
3390 sprintf(input, "print ABCDefg\n");
3391 string_eval(in, out, input, 80, output);
3392 string_eval_teardown(&in, &out);
3395 START_TEST(test_simple_eval)
3398 char input[80], output[80]={};
3399 string_eval_setup(&in, &out);
3400 sprintf(input, "print 3+4*5\n");
3401 string_eval(in, out, input, 80, output);
3402 fail_unless(STRMATCH(output,"23\n"), NULL);
3403 string_eval_teardown(&in, &out);
3406 START_TEST(test_multiple_evals)
3409 char input[80], output[80]={};
3410 string_eval_setup(&in, &out);
3411 sprintf(input, "print 3+4*5\n");
3412 string_eval(in, out, input, 80, output);
3413 fail_unless(STRMATCH(output,"23\n"), NULL);
3414 sprintf(input, "print (3**2 + 4**2)**0.5\n");
3415 string_eval(in, out, input, 80, output);
3416 fail_unless(STRMATCH(output,"5.0\n"), NULL);
3417 string_eval_teardown(&in, &out);
3420 START_TEST(test_eval_with_state)
3423 char input[80], output[80]={};
3424 string_eval_setup(&in, &out);
3425 sprintf(input, "print 3+4*5\n");
3426 fprintf(in, "A = 3\n");
3427 sprintf(input, "print A*3\n");
3428 string_eval(in, out, input, 80, output);
3429 fail_unless(STRMATCH(output,"9\n"), NULL);
3430 string_eval_teardown(&in, &out);
3434 <<string eval test case defs>>=
3435 TCase *tc_string_eval = tcase_create("string_eval");
3437 <<string eval test case add>>=
3438 tcase_add_test(tc_string_eval, test_setup_teardown);
3439 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
3440 tcase_add_test(tc_string_eval, test_simple_eval);
3441 tcase_add_test(tc_string_eval, test_multiple_evals);
3442 tcase_add_test(tc_string_eval, test_eval_with_state);
3443 suite_add_tcase(s, tc_string_eval);
3447 \section{Accelerating function evaluation}
3449 My first version-0.3 code was running very slowly.
3450 With the options suggested in the help
3451 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
3452 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
3453 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
3454 That's only 15 calls per solution, so the search algorithm seems reasonable.
3455 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
3460 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
3462 #endif /* ACCEL_K_H */
3465 <<accel k module makefile lines>>=
3466 NW_SAWSIM_MODS += accel_k
3467 #CHECK_BINS += check_accel_k
3468 check_accel_k_MODS =
3472 #include <assert.h> /* assert() */
3473 #include <stdlib.h> /* realloc(), free(), NULL */
3474 #include "global.h" /* environment_t */
3475 #include "k_model.h" /* k_func_t */
3476 #include "interp.h" /* interp_* */
3477 #include "accel_k.h"
3479 typedef struct accel_k_struct {
3480 interp_table_t *itable;
3486 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3487 static int num_accels = 0;
3488 static accel_k_t *accels=NULL;
3490 /* Wrap k in the standard f(x) acceleration form */
3491 static double k_wrap(double F, void *params)
3493 accel_k_t *p = (accel_k_t *)params;
3494 return (*p->k)(F, p->env, p->params);
3497 static int k_tol(double FA, double kA, double FB, double kB)
3500 if (FB-FA > 1e-12) {
3501 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3502 return 1; /* unnacceptable */
3504 //printf("acceptable tol\n");
3505 return 0; /* acceptable */
3509 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3512 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3513 assert(accels != NULL);
3514 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3516 accels[i].env = env;
3517 accels[i].params = params;
3524 for (i=0; i<num_accels; i++)
3525 interp_table_free(accels[i].itable);
3527 if (accels) free(accels);
3531 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3534 for (i=0; i<num_accels; i++) {
3535 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3536 /* We've already setup this function.
3537 * WARNING: we're assuming param and env are the same. */
3538 return interp_table_eval(accels[i].itable, accels+i, F);
3541 /* set up a new acceleration environment */
3542 i = add_accel_k(k, env, params);
3543 return interp_table_eval(accels[i].itable, accels+i, F);
3547 \section{Tension models}
3548 \label{sec.tension_models}
3550 TODO: tension model intro.
3551 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.
3553 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3554 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]].
3556 <<tension-model.h>>=
3557 #ifndef TENSION_MODEL_H
3558 #define TENSION_MODEL_H
3560 <<tension handler types>>
3561 <<constant tension model declarations>>
3562 <<hooke tension model declarations>>
3563 <<worm-like chain tension model declarations>>
3564 <<freely-jointed chain tension model declarations>>
3565 <<find tension definitions>>
3566 <<tension model global declarations>>
3567 #endif /* TENSION_MODEL_H */
3570 <<tension model module makefile lines>>=
3571 NW_SAWSIM_MODS += tension_model
3572 #CHECK_BINS += check_tension_model
3573 check_tension_model_MODS =
3575 <<tension model utils makefile lines>>=
3576 TENSION_MODEL_MODS = tension_model parse list tension_balance
3577 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3578 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3579 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3580 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
3581 notangle -Rtension-model-utils.c $< > $@
3582 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
3583 gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3584 clean_tension_model_utils :
3585 rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3586 @ The pipe symbol [[|]] marks the prerequisits following it (in this
3587 case, the directories) as ‘order-only’ prerequisites. The timestamp
3588 on these prerequisits does not effect whether the rules are executed.
3589 This is appropriate for directories, where we don't need to recompile
3590 after an unrelated has been added to the directory, but only when the
3591 source prerequisites change. See the [[make]] documentation for more
3593 (\url{http://www.gnu.org/software/make/manual/html_node/Prerequisite-Types.html}).
3596 \label{sec.null_tension}
3598 For unstretchable domains.
3600 <<null tension model getopt>>=
3601 {NULL, NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3604 \subsection{Constant}
3605 \label{sec.const_tension}
3607 An infinitely stretchable domain providing a constant tension.
3608 Obviously this cannot be inverted, so you can't put this domain in
3609 series with any other domains. We include it mostly for testing
3612 <<constant tension functions>>=
3613 <<constant tension function>>
3614 <<constant tension parameter create/destroy functions>>
3617 <<constant tension model declarations>>=
3618 <<constant tension function declaration>>
3619 <<constant tension parameter create/destroy function declarations>>
3620 <<constant tension model global declarations>>
3624 <<constant tension function declaration>>=
3625 extern double const_tension_handler(double x, void *pdata);
3627 <<constant tension function>>=
3628 double const_tension_handler(double x, void *pdata)
3630 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3631 list_t *list = data->group_tension_params;
3636 assert(list != NULL); /* empty active group?! */
3637 F = ((const_tension_param_t *)list->d)->F;
3639 fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3641 while (list != NULL) {
3642 assert(((const_tension_param_t *)list->d)->F == F);
3650 <<constant tension structure>>=
3651 typedef struct const_tension_param_struct {
3652 double F; /* tension (force) in N */
3653 } const_tension_param_t;
3657 <<constant tension parameter create/destroy function declarations>>=
3658 extern void *string_create_const_tension_param_t(char **param_strings);
3659 extern void destroy_const_tension_param_t(void *p);
3662 <<constant tension parameter create/destroy functions>>=
3663 const_tension_param_t *create_const_tension_param_t(double F)
3665 const_tension_param_t *ret
3666 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3667 assert(ret != NULL);
3672 void *string_create_const_tension_param_t(char **param_strings)
3674 return create_const_tension_param_t(atof(param_strings[0]));
3677 void destroy_const_tension_param_t(void *p)
3685 <<constant tension model global declarations>>=
3686 extern int num_const_tension_params;
3687 extern char *const_tension_param_descriptions[];
3688 extern char const_tension_param_string[];
3691 <<constant tension model globals>>=
3692 int num_const_tension_params = 1;
3693 char *const_tension_param_descriptions[] = {"tension F, N"};
3694 char const_tension_param_string[] = "0";
3697 <<constant tension model getopt>>=
3698 {&const_tension_handler, NULL, "const", "an infinitely stretchable domain with constant tension", 1, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t}
3704 <<hooke functions>>=
3705 <<internal hooke functions>>
3707 <<inverse hooke handler>>
3708 <<hooke parameter create/destroy functions>>
3711 <<hooke tension model declarations>>=
3712 <<hooke tension function declaration>>
3713 <<hooke tension parameter create/destroy function declarations>>
3714 <<hooke tension model global declarations>>
3718 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
3719 The behavior of a series of springs $k_i$ in series is given by
3721 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
3723 For a simple proof, see Appendix \ref{sec.math_hooke}.
3725 <<hooke tension function declaration>>=
3726 extern double hooke_handler(double x, void *pdata);
3727 extern double inverse_hooke_handler(double F, void *pdata);
3730 First we define a function that computes the effective spring constant
3731 (stored in a single [[hooke_param_t]]) for the entire group.
3732 <<internal hooke functions>>=
3733 static void hooke_param_grouper(tension_handler_data_t *data,
3734 hooke_param_t *param) {
3735 list_t *list = data->group_tension_params;
3739 assert(list != NULL); /* empty active group?! */
3740 while (list != NULL) {
3741 assert( ((hooke_param_t *)list->d)->k > 0 );
3742 k += 1.0/ ((hooke_param_t *)list->d)->k;
3744 fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
3745 __FUNCTION__, ((hooke_param_t *)list->d)->k);
3751 fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N (from data %p)\n",
3752 __FUNCTION__, k, x, k*x, data);
3759 Everyone knows Hooke's law: $F=kx$.
3761 double hooke_handler(double x, void *pdata)
3763 hooke_param_t param = {0};
3766 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
3772 Inverting Hooke's law gives $x=F/k$.
3773 <<inverse hooke handler>>=
3774 double inverse_hooke_handler(double F, void *pdata)
3776 hooke_param_t param = {0};
3779 hooke_param_grouper((tension_handler_data_t *)pdata, ¶m);
3785 Not much to keep track of for a single effective spring, just the
3786 spring constant $k$.
3787 <<hooke structure>>=
3788 typedef struct hooke_param_struct {
3789 double k; /* spring constant in N/m */
3794 We wrap up our Hooke implementation with some book-keeping functions.
3795 <<hooke tension parameter create/destroy function declarations>>=
3796 extern void *string_create_hooke_param_t(char **param_strings);
3797 extern void destroy_hooke_param_t(void *p);
3801 <<hooke parameter create/destroy functions>>=
3802 hooke_param_t *create_hooke_param_t(double k)
3804 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
3805 assert(ret != NULL);
3810 void *string_create_hooke_param_t(char **param_strings)
3812 return create_hooke_param_t(atof(param_strings[0]));
3815 void destroy_hooke_param_t(void *p)
3822 <<hooke tension model global declarations>>=
3823 extern int num_hooke_params;
3824 extern char *hooke_param_descriptions[];
3825 extern char hooke_param_string[];
3828 <<hooke tension model globals>>=
3829 int num_hooke_params = 1;
3830 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
3831 char hooke_param_string[]="0.05";
3834 <<hooke tension model getopt>>=
3835 {hooke_handler, inverse_hooke_handler, "hooke", "a Hookean spring (F=kx)", num_hooke_params, hooke_param_descriptions, (char *)hooke_param_string, &string_create_hooke_param_t, &destroy_hooke_param_t}
3838 \subsection{Worm-like chain}
3841 We can model several unfolded domains with a single worm-like chain.
3842 <<worm-like chain functions>>=
3843 <<internal worm-like chain functions>>
3844 <<worm-like chain handler>>
3845 <<inverse worm-like chain handler>>
3846 <<worm-like chain create/destroy functions>>
3849 <<worm-like chain tension model declarations>>=
3851 <<worm-like chain handler declaration>>
3852 <<inverse worm-like chain handler declaration>>
3853 <<worm-like chain create/destroy function declarations>>
3854 <<worm-like chain tension model global declarations>>
3858 <<internal worm-like chain functions>>=
3859 <<worm-like chain function>>
3860 <<inverse worm-like chain function>>
3861 <<worm-like chain parameter grouper>>
3864 The combination of all unfolded domains is modeled as a worm like chain,
3865 with the force $F_{\text{WLC}}$ approximately given by
3867 F_{\text{WLC}} \approx \frac{k_B T}{p}
3868 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
3870 where $x$ is the distance between the fixed ends,
3871 $k_B$ is Boltzmann's constant,
3872 $T$ is the absolute temperature,
3873 $p$ is the persistence length, and
3874 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
3875 <<worm-like chain function>>=
3876 static double wlc(double x, double T, double p, double L)
3878 double xL=0.0; /* uses global k_B */
3880 assert(T > 0); assert(p > 0); assert(L > 0);
3881 if (x >= L) return HUGE_VAL;
3882 xL = x/L; /* unitless */
3883 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
3884 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
3885 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
3890 Inverting the approximate WLC is fairly straightforward, if a bit tedious.
3891 Letting $x_L \equiv x/L$ and $F_T \equiv Fp/k_B T$ we have
3893 F_T &\approx \frac{1}{4} \p({\frac{1}{(1-x_L)^2} - 1}) + x_L \\
3894 F_T (1-x_L)^2 &\approx \frac{1}{4} - \frac{1}{4} (1-x_L)^2 + x_L(1-x_L)^2 \\
3895 0 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) (1-x_L)^2 + x_L(1-x_L)^2 \\
3896 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) \p({1-2x_L+x_L^2})
3897 + x_L - 2x_L^2 + x_L^3 \\
3898 &\approx \p[{\frac{1}{4} - \p({F_T + \frac{1}{4}})}]
3899 + x_L \p[{2F_T + \frac{1}{2} + 1}]
3900 + x_L^2 \p[{-\p({F_T + \frac{1}{4}}) - 2}]
3903 + x_L \p({2F_T + \frac{3}{2}})
3904 - x_L^2 \p({F_T + \frac{9}{4}})
3906 &\approx ax_L^3 + bx_L^2 + cx_L + d
3908 where $a \equiv 1$, $b \equiv -F_T - 9/4$, $c \equiv 2F_T + 3/2$, and
3910 % From \citet{wikipedia_cubic_function} the discriminant
3912 % \Delta &= -4b^3d + b^2c^2 - 4ac^3 + 18abcd - 27a^2d^2 \\
3913 % &= -4F_T\p({F_T + \frac{9}{4}})^3
3914 % + \p({F_T + \frac{9}{4}})^2 \cdot \p({2F_T + \frac{3}{2}})^2
3915 % - 4 \p({2F_T + \frac{3}{2}})^3
3916 % + 18F_T \p({F_T + \frac{9}{4}}) \cdot \p({2F_T + \frac{3}{2}})
3918 % &= -4F_T \p({F_T^3 + \frac{27}{4}F_T^2 +\frac{243}{16}F_T +\frac{729}{64}})
3919 %% a^3 + 3a^2b + 3ab^2 + b^3
3920 % + \p({F_T^2 + \frac{9}{2}F_T + \frac{81}{16}})
3921 % \cdot \p({4F_T^2 + 6F_T + \frac{9}{4}}) \\
3922 % &\qquad - 4 \p({8F_T^3 + 18F_T^2 + \frac{27}{2}F_T + \frac{27}{8}})
3923 %% a^3 + 3a^2b + 3ab^2 + b^3
3924 % + 18F_T \p({2F_T^2 + 6F_T + \frac{27}{8}})
3926 % &= -4F_T^4 - 27F_T^3 - \frac{243}{4}F_T^2 - \frac{729}{16}F_T
3927 % + 4F_T^4 + 24F_T^3 +\frac{99}{2}F_T^2 +\frac{81}{2}F_T +\frac{729}{64}\\
3928 % &\qquad - 32F_T^3 - 72F_T^2 - 54F_T - \frac{27}{2}
3929 % + 36F_T^3 + 108F_T^2 + \frac{243}{4}F_T
3931 % &= F_T^4 \p({-4 + 4}) + F_T^3 \p({-27 + 24 - 32 + 36})
3932 % + F_T^2 \p({-\frac{243}{4} + \frac{99}{2} - 72 + 108 -27}) \\
3933 % &\qquad + F_T \p({-\frac{729}{16} + \frac{81}{2} - 54 + \frac{243}{4}})
3934 % \p({\frac{729}{64} - \frac{27}{2}}) \\
3935 % &= F_T^3 - \frac{9}{4}F_T^2 + \frac{27}{16}F_T - \frac{135}{64}
3937 We can plug these coefficients into the GSL cubic solver to invert the
3938 WLC\citep{galassi05}. The applicable root is always the one which
3939 comes before the singularity, which will be the smallest real root.
3940 <<inverse worm-like chain function>>=
3941 static double inverse_wlc(double F, double T, double p, double L)
3943 double FT = F*p/(k_B*T); /* uses global k_B */
3944 double xL0, xL1, xL2;
3947 assert(T > 0); assert(p > 0); assert(L > 0);
3948 num_roots = gsl_poly_solve_cubic(-(FT+2.25),2*FT+1.5,-FT, &xL0, &xL1, &xL2);
3949 assert(num_roots > 0); assert(xL0 >= -DOUBLE_PRECISION); assert(xL0 < 1);
3950 if (xL0 < 0) xL0 = 0.0;
3956 First we define a function that computes the effective WLC parameters
3957 (stored in a single [[wlc_param_t]]) for the entire group.
3958 <<worm-like chain parameter grouper>>=
3959 static void wlc_param_grouper(tension_handler_data_t *data,
3960 wlc_param_t *param) {
3961 list_t *list = data->group_tension_params;
3962 double p=0.0, L=0.0;
3965 assert(list != NULL); /* empty active group?! */
3966 p = ((wlc_param_t *)list->d)->p;
3967 while (list != NULL) {
3968 /* enforce consistent p */
3969 assert( ((wlc_param_t *)list->d)->p == p);
3970 L += ((wlc_param_t *)list->d)->L;
3972 fprintf(stderr, "%s: WLC section %g m long in series\n",
3973 __FUNCTION__, ((wlc_param_t *)list->d)->L);
3978 fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
3979 __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
3987 <<worm-like chain handler declaration>>=
3988 extern double wlc_handler(double x, void *pdata);
3991 This model requires that each unfolded domain in the group have the
3992 same persistence length.
3993 <<worm-like chain handler>>=
3994 double wlc_handler(double x, void *pdata)
3996 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3997 wlc_param_t param = {0};
3999 wlc_param_grouper(data, ¶m);
4000 return wlc(x, data->env->T, param.p, param.L);
4005 <<inverse worm-like chain handler declaration>>=
4006 extern double inverse_wlc_handler(double F, void *pdata);
4009 <<inverse worm-like chain handler>>=
4010 double inverse_wlc_handler(double F, void *pdata)
4012 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4013 wlc_param_t param = {0};
4015 wlc_param_grouper(data, ¶m);
4016 return inverse_wlc(F, data->env->T, param.p, param.L);
4021 <<worm-like chain structure>>=
4022 typedef struct wlc_param_struct {
4023 double p; /* WLC persistence length (m) */
4024 double L; /* WLC contour length (m) */
4028 <<worm-like chain create/destroy function declarations>>=
4029 extern void *string_create_wlc_param_t(char **param_strings);
4030 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
4033 <<worm-like chain create/destroy functions>>=
4034 wlc_param_t *create_wlc_param_t(double p, double L)
4036 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
4037 assert(ret != NULL);
4043 void *string_create_wlc_param_t(char **param_strings)
4045 return create_wlc_param_t(atof(param_strings[0]), atof(param_strings[1]));
4048 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
4056 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.
4057 TODO (now we copy the string before parsing, but still do this for clarity).
4059 <<valid string write code>>=
4060 char string[] = "abc";
4063 <<valid string write code>>=
4064 char *string = "abc";
4068 <<worm-like chain tension model global declarations>>=
4069 extern int num_wlc_params;
4070 extern char *wlc_param_descriptions[];
4071 extern char wlc_param_string[];
4073 <<worm-like chain tension model globals>>=
4074 int num_wlc_params = 2;
4075 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
4076 char wlc_param_string[]="0.39e-9,28.4e-9";
4079 <<worm-like chain tension model getopt>>=
4080 {wlc_handler, inverse_wlc_handler, "wlc", "a worm-like chain (F=kbT/p{.25[(1-x/L)^-2 -1] + x/L})", num_wlc_params, wlc_param_descriptions, (char *)wlc_param_string, &string_create_wlc_param_t, &destroy_wlc_param_t}
4082 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
4084 \subsection{Freely-jointed chain}
4087 <<freely-jointed chain functions>>=
4088 <<freely-jointed chain function>>
4089 <<freely-jointed chain handler>>
4090 <<freely-jointed chain create/destroy functions>>
4093 <<freely-jointed chain tension model declarations>>=
4094 <<freely-jointed chain handler declaration>>
4095 <<freely-jointed chain create/destroy function declarations>>
4096 <<freely-jointed chain tension model global declarations>>
4100 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
4101 The tension of a chain of $N$ such links is given by
4103 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
4105 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}.
4106 <<freely-jointed chain function>>=
4107 double langevin(double x, void *params)
4110 return 1.0/tanh(x) - 1/x;
4112 <<one dimensional bracket declaration>>
4113 <<one dimensional solve declaration>>
4114 double inv_langevin(double x)
4116 double lb=0.0, ub=1.0, ret;
4117 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
4118 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
4122 double fjc(double x, double T, double l, int N)
4124 double L = l*(double)N;
4125 if (x == 0) return 0;
4126 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
4127 return k_B*T/l * inv_langevin(x/L);
4131 <<freely-jointed chain handler declaration>>=
4132 extern double fjc_handler(double x, void *pdata);
4134 <<freely-jointed chain handler>>=
4135 double fjc_handler(double x, void *pdata)
4137 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
4138 list_t *list = data->group_tension_params;
4143 assert(list != NULL); /* empty active group?! */
4144 l = ((fjc_param_t *)list->d)->l;
4145 while (list != NULL) {
4146 /* enforce consistent l */
4147 assert( ((fjc_param_t *)list->d)->l == l);
4148 N += ((fjc_param_t *)list->d)->N;
4150 fprintf(stderr, "%s: FJC section %d links long in series\n",
4151 __FUNCTION__, ((fjc_param_t *)list->d)->N);
4156 fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
4157 __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
4159 return fjc(x, data->env->T, l, N);
4163 <<freely-jointed chain structure>>=
4164 typedef struct fjc_param_struct {
4165 double l; /* FJC link length (m) */
4166 int N; /* FJC number of links */
4170 <<freely-jointed chain create/destroy function declarations>>=
4171 extern void *string_create_fjc_param_t(char **param_strings);
4172 extern void destroy_fjc_param_t(void *p);
4175 <<freely-jointed chain create/destroy functions>>=
4176 fjc_param_t *create_fjc_param_t(double l, double N)
4178 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
4179 assert(ret != NULL);
4185 void *string_create_fjc_param_t(char **param_strings)
4187 return create_fjc_param_t(atof(param_strings[0]), atof(param_strings[1]));
4190 void destroy_fjc_param_t(void *p)
4197 <<freely-jointed chain tension model global declarations>>=
4198 extern int num_fjc_params;
4199 extern char *fjc_param_descriptions[];
4200 extern char fjc_param_string[];
4203 <<freely-jointed chain tension model globals>>=
4204 int num_fjc_params = 2;
4205 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
4206 char fjc_param_string[]="0.5e-9,200";
4209 <<freely-jointed chain tension model getopt>>=
4210 {fjc_handler, NULL, "fjc", "a freely-jointed chain [F=kbT/l L^-1(x/L)]", num_fjc_params, fjc_param_descriptions, (char *)fjc_param_string, &string_create_fjc_param_t, &destroy_fjc_param_t}
4213 \subsection{Tension model implementation}
4215 <<tension-model.c>>=
4218 <<tension model includes>>
4219 #include "tension_model.h"
4220 <<tension model internal definitions>>
4221 <<tension model globals>>
4222 <<tension model functions>>
4225 <<tension model includes>>=
4226 #include <assert.h> /* assert() */
4227 #include <stdlib.h> /* NULL */
4228 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4229 #include <stdio.h> /* fprintf(), stdout */
4232 #include "tension_balance.h" /* oneD_*() */
4235 <<tension model internal definitions>>=
4236 <<constant tension structure>>
4238 <<worm-like chain structure>>
4239 <<freely-jointed chain structure>>
4242 <<tension model globals>>=
4243 <<tension handler array global>>
4244 <<constant tension model globals>>
4245 <<hooke tension model globals>>
4246 <<worm-like chain tension model globals>>
4247 <<freely-jointed chain tension model globals>>
4250 <<tension model functions>>=
4251 <<constant tension functions>>
4253 <<worm-like chain functions>>
4254 <<freely-jointed chain functions>>
4258 \subsection{Utilities}
4260 The tension models can get complicated, and users may want to reassure
4261 themselves that this portion of the input physics are functioning properly. The
4262 stand-alone program [[tension_model_utils]] demonstrates the tension model
4263 interface without the overhead of having to understand a full simulation.
4264 [[tension_model_utils]] takes command line model arguments like the full
4265 simulation, and prints $F(x)$ data to the screen for the selected model over a
4268 <<tension-model-utils.c>>=
4270 <<tension model utility includes>>
4271 <<tension model utility definitions>>
4272 <<tension model utility getopt functions>>
4274 <<tension model utility main>>
4277 <<tension model utility main>>=
4278 int main(int argc, char **argv)
4280 <<tension model getopt array definition>>
4281 tension_model_getopt_t *model = NULL;
4285 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
4286 tension_handler_data_t tdata;
4287 int num_param_args; /* for INIT_MODEL() */
4288 char **param_args; /* for INIT_MODEL() */
4290 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model,
4291 &Fmax, &Xmax, &flags);
4293 if (flags & VFLAG) {
4294 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4296 INIT_MODEL("utility", model, model->params, params);
4297 tdata.group_tension_params = NULL;
4298 push(&tdata.group_tension_params, params, 1);
4300 tdata.persist = NULL;
4301 if (model->handler == NULL) {
4302 printf("No tension function for model %s\n", model->name);
4308 printf("#Distance (m)\tForce (N)\n");
4309 for (i=0; i<=N; i++) {
4310 x = Xmax*i/(double)N;
4311 F = (*model->handler)(x, &tdata);
4312 if (F < 0 || F > Fmax) break;
4313 printf("%g\t%g\n", x, F);
4315 if (flags & VFLAG && i <= N) {
4316 /* explain exit condition */
4318 printf("Impossible force %g\n", F);
4320 printf("Reached large force limit %g > %g\n", F, Fmax);
4323 params = pop(&tdata.group_tension_params);
4324 if (model->destructor)
4325 (*model->destructor)(params);
4330 <<tension model utility includes>>=
4333 #include <string.h> /* strlen() */
4334 #include <assert.h> /* assert() */
4338 #include "tension_model.h"
4341 <<tension model utility definitions>>=
4342 <<version definition>>
4343 #define VFLAG 1 /* verbose */
4344 <<string comparison definition and globals>>
4345 <<tension model getopt definitions>>
4346 <<initialize model definition>>
4350 <<tension model utility getopt functions>>=
4351 <<index tension model>>
4352 <<tension model utility help>>
4353 <<tension model utility get options>>
4356 <<tension model utility help>>=
4357 void help(char *prog_name,
4359 int n_tension_models, tension_model_getopt_t *tension_models,
4360 int tension_model, double Fmax, double Xmax)
4363 printf("usage: %s [options]\n", prog_name);
4364 printf("Version %s\n\n", VERSION);
4365 printf("A utility for understanding the available tension models\n\n");
4366 printf("Environmental options:\n");
4367 printf("-T T\tTemperature T (currently %g K)\n", env->T);
4368 printf("-C T\tYou can also set the temperature T in Celsius\n");
4369 printf("Model options:\n");
4370 printf("-m model\tFolded tension model (currently %s)\n",
4371 tension_models[tension_model].name);
4372 printf("-a args\tFolded tension model argument string (currently %s)\n",
4373 tension_models[tension_model].params);
4374 printf("F(x) is calculated for a range of x and printed\n");
4375 printf("For example:\n");
4376 printf(" #Distance (m)\tForce (N)\n");
4377 printf(" 123.456\t7.89\n");
4379 printf("-F\tSet the maximum F value for the standard mode F(x) output (currently %g)\n", Fmax);
4380 printf("-X\tSet the maximum x value for the standard mode F(x) output (currently %g)\n", Xmax);
4381 printf("-V\tChange output to verbose mode\n");
4382 printf("-h\tPrint this help and exit\n");
4384 printf("Tension models:\n");
4385 for (i=0; i<n_tension_models; i++) {
4386 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
4387 for (j=0; j < tension_models[i].num_params ; j++)
4388 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
4389 printf("\t\tdefaults: %s\n", tension_models[i].params);
4394 <<tension model utility get options>>=
4395 void get_options(int argc, char **argv, environment_t *env,
4396 int n_tension_models, tension_model_getopt_t *tension_models,
4397 tension_model_getopt_t **model, double *Fmax, double *Xmax,
4398 unsigned int *flags)
4400 char *prog_name = NULL;
4401 char c, options[] = "T:C:m:a:F:X:Vh";
4402 int tension_model=0, num_strings;
4403 extern char *optarg;
4404 extern int optind, optopt, opterr;
4408 /* setup defaults */
4410 prog_name = argv[0];
4411 env->T = 300.0; /* K */
4415 *model = tension_models;
4417 while ((c=getopt(argc, argv, options)) != -1) {
4419 case 'T': env->T = atof(optarg); break;
4420 case 'C': env->T = atof(optarg)+273.15; break;
4422 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
4423 *model = tension_models+tension_model;
4426 tension_models[tension_model].params = optarg;
4428 case 'F': *Fmax = atof(optarg); break;
4429 case 'X': *Xmax = atof(optarg); break;
4430 case 'V': *flags |= VFLAG; break;
4432 fprintf(stderr, "unrecognized option '%c'\n", optopt);
4433 /* fall through to default case */
4435 help(prog_name, env, n_tension_models, tension_models, tension_model, *Fmax, *Xmax);
4444 \section{Unfolding rate models}
4445 \label{sec.k_models}
4447 We have two main choices for determining $k$: Bell's model, which assumes the
4448 folded domains exist in a single `folded' state and make instantaneous,
4449 stochastic transitions over a free energy barrier; and Kramers' model, which
4450 treats the unfolding as a thermalized diffusion process.
4451 We deal with these options like we dealt with the different tension models: by bundling all model-specific
4452 parameters into structures, with model specific functions for determining $k$.
4454 <<k func definition>>=
4455 typedef double k_func_t(double F, environment_t *env, void *params);
4458 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
4459 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]].
4461 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
4462 because \LaTeX\ doesn't like underscores outside math mode.
4463 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
4464 You could use spaces instead (`\verb+ +'),
4465 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
4466 which I don't like as much.
4472 <<k func definition>>
4473 <<null k declarations>>
4474 <<const k declarations>>
4475 <<bell k declarations>>
4476 <<kbell k declarations>>
4477 <<kramers k declarations>>
4478 <<kramers integ k declarations>>
4479 #endif /* K_MODEL_H */
4482 <<k model module makefile lines>>=
4483 NW_SAWSIM_MODS += k_model
4484 CHECK_BINS += check_k_model
4485 check_k_model_MODS = parse string_eval
4487 [[check_k_model]] also depends on the parse module.
4489 <<k model utils makefile lines>>=
4490 K_MODEL_MODS = k_model parse string_eval
4491 K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
4492 $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h)
4493 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
4494 notangle -Rk-model-utils.c $< > $@
4495 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
4496 gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
4497 clean_k_model_utils :
4498 rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
4504 For domains that are never folded.
4506 <<null k declarations>>=
4508 <<null k functions>>=
4513 <<null k model getopt>>=
4514 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
4517 \subsection{Constant rate model}
4520 This is a pretty boring model, but it is usefull for testing the engine.
4524 where $k_0$ is the constant unfolding rate.
4526 <<const k functions>>=
4527 <<const k function>>
4528 <<const k structure create/destroy functions>>
4531 <<const k declarations>>=
4532 <<const k function declaration>>
4533 <<const k structure create/destroy function declarations>>
4534 <<const k global declarations>>
4538 <<const k structure>>=
4539 typedef struct const_k_param_struct {
4540 double knot; /* transition rate k_0 (frac population per s) */
4544 <<const k function declaration>>=
4545 double const_k(double F, environment_t *env, void *const_k_params);
4547 <<const k function>>=
4548 double const_k(double F, environment_t *env, void *const_k_params)
4549 { /* Returns the rate constant k in frac pop per s. */
4550 const_k_param_t *p = (const_k_param_t *) const_k_params;
4552 assert(p->knot > 0);
4557 <<const k structure create/destroy function declarations>>=
4558 void *string_create_const_k_param_t(char **param_strings);
4559 void destroy_const_k_param_t(void *p);
4562 <<const k structure create/destroy functions>>=
4563 const_k_param_t *create_const_k_param_t(double knot)
4565 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
4566 assert(ret != NULL);
4571 void *string_create_const_k_param_t(char **param_strings)
4573 return create_const_k_param_t(atof(param_strings[0]));
4576 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
4583 <<const k global declarations>>=
4584 extern int num_const_k_params;
4585 extern char *const_k_param_descriptions[];
4586 extern char const_k_param_string[];
4589 <<const k globals>>=
4590 int num_const_k_params = 1;
4591 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
4592 char const_k_param_string[]="1";
4595 <<const k model getopt>>=
4596 {"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}
4599 \subsection{Bell's model}
4602 Quantitatively, Bell's model gives $k$ as
4604 k = k_0 \cdot e^{F dx / k_B T} \;,
4606 where $k_0$ is the unforced unfolding rate,
4607 $F$ is the applied unfolding force,
4608 $dx$ is the distance to the transition state, and
4609 $k_B T$ is the thermal energy\citep{schlierf06}.
4611 <<bell k functions>>=
4613 <<bell k structure create/destroy functions>>
4616 <<bell k declarations>>=
4617 <<bell k function declaration>>
4618 <<bell k structure create/destroy function declarations>>
4619 <<bell k global declarations>>
4623 <<bell k structure>>=
4624 typedef struct bell_param_struct {
4625 double knot; /* transition rate k_0 (frac population per s) */
4626 double dx; /* `distance to transition state' (nm) */
4630 <<bell k function declaration>>=
4631 double bell_k(double F, environment_t *env, void *bell_params);
4633 <<bell k function>>=
4634 double bell_k(double F, environment_t *env, void *bell_params)
4635 { /* Returns the rate constant k in frac pop per s.
4636 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4637 * uses global k_B in J/K */
4638 bell_param_t *p = (bell_param_t *) bell_params;
4639 assert(F >= 0); assert(env->T > 0);
4641 assert(p->knot > 0); //assert(p->dx >= 0); /* p < 0 for refolding reactions, etc. */
4643 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
4644 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
4645 p->knot * exp(F*p->dx / (k_B*env->T) ));
4646 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
4648 return p->knot * exp(F*p->dx / (k_B*env->T) );
4652 <<bell k structure create/destroy function declarations>>=
4653 void *string_create_bell_param_t(char **param_strings);
4654 void destroy_bell_param_t(void *p);
4657 <<bell k structure create/destroy functions>>=
4658 bell_param_t *create_bell_param_t(double knot, double dx)
4660 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
4661 assert(ret != NULL);
4667 void *string_create_bell_param_t(char **param_strings)
4669 return create_bell_param_t(atof(param_strings[0]), atof(param_strings[1]));
4672 void destroy_bell_param_t(void *p /* bell_param_t * */)
4679 <<bell k global declarations>>=
4680 extern int num_bell_params;
4681 extern char *bell_param_descriptions[];
4682 extern char bell_param_string[];
4686 int num_bell_params = 2;
4687 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4688 char bell_param_string[]="3.3e-4,0.25e-9";
4691 <<bell k model getopt>>=
4692 {"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}
4694 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4695 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4698 \subsection{Linker stiffness corrected Bell model}
4701 Walton et. al showed that the Bell model constant force approximation
4702 doesn't hold for biotin-streptavadin binding, and I extended their
4703 results to I27 for the 2009 Biophysical Society Annual
4704 Meeting\cite{walton08,king09}. More details to follow when I get done
4705 with the conference\ldots
4707 We adjust Bell's model to
4709 k = k_0 \cdot e^{\frac{F dx - \frac{1}{2} \kappa dx^2}{k_B T}} \;,
4711 where $k_0$ is the unforced unfolding rate,
4712 $F$ is the applied unfolding force,
4713 $\kappa$ is the effective spring constant,
4714 $dx$ is the distance to the transition state, and
4715 $k_B T$ is the thermal energy\citep{schlierf06}.
4717 <<kbell k functions>>=
4718 <<kbell k function>>
4719 <<kbell k structure create/destroy functions>>
4722 <<kbell k declarations>>=
4723 <<kbell k function declaration>>
4724 <<kbell k structure create/destroy function declarations>>
4725 <<kbell k global declarations>>
4729 <<kbell k structure>>=
4730 typedef struct kbell_param_struct {
4731 double knot; /* transition rate k_0 (frac population per s) */
4732 double dx; /* `distance to transition state' (nm) */
4736 <<kbell k function declaration>>=
4737 double kbell_k(double F, environment_t *env, void *kbell_params);
4739 <<kbell k function>>=
4740 double kbell_k(double F, environment_t *env, void *kbell_params)
4741 { /* Returns the rate constant k in frac pop per s.
4742 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4743 * uses global k_B in J/K */
4744 kbell_param_t *p = (kbell_param_t *) kbell_params;
4745 assert(F >= 0); assert(env->T > 0); assert(env->stiffness >= 0);
4747 assert(p->knot > 0); assert(p->dx >= 0);
4748 return p->knot * exp((F*p->dx - 0.5*env->stiffness*(p->dx)*(p->dx))/(k_B*env->T) );
4752 <<kbell k structure create/destroy function declarations>>=
4753 void *string_create_kbell_param_t(char **param_strings);
4754 void destroy_kbell_param_t(void *p);
4757 <<kbell k structure create/destroy functions>>=
4758 kbell_param_t *create_kbell_param_t(double knot, double dx)
4760 kbell_param_t *ret = (kbell_param_t *) malloc(sizeof(kbell_param_t));
4761 assert(ret != NULL);
4767 void *string_create_kbell_param_t(char **param_strings)
4769 return create_kbell_param_t(atof(param_strings[0]), atof(param_strings[1]));
4772 void destroy_kbell_param_t(void *p /* kbell_param_t * */)
4779 <<kbell k global declarations>>=
4780 extern int num_kbell_params;
4781 extern char *kbell_param_descriptions[];
4782 extern char kbell_param_string[];
4785 <<kbell k globals>>=
4786 int num_kbell_params = 2;
4787 char *kbell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4788 char kbell_param_string[]="3.3e-4,0.25e-9";
4791 <<kbell k model getopt>>=
4792 {"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}
4794 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4795 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4798 \subsection{Kramer's model}
4801 Kramer's model gives $k$ as
4803 % \frac{1}{k} = \frac{1}{D}
4804 % \int_{x_\text{min}}^{x_\text{max}}
4805 % e^{\frac{-U_F(x)}{k_B T}}
4806 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4809 %where $D$ is the diffusion constant,
4810 %$U_F(x)$ is the free energy along the unfolding coordinate, and
4811 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4812 % before the minimum and after the maximum, respectively \citep{schlierf06}.
4814 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
4816 where $D$ is the diffusion constant,
4818 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
4820 are length scales where
4821 $x_c(F)$ is the location of the bound state and
4822 $x_{ts}(F)$ is the location of the transition state,
4823 $E(F,x)$ is the free energy, and
4824 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
4825 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
4826 \citet{evans97} Eqn.~3).
4828 <<kramers k functions>>=
4829 <<kramers k function>>
4830 <<kramers k structure create/destroy functions>>
4833 <<kramers k declarations>>=
4834 <<kramers k function declaration>>
4835 <<kramers k structure create/destroy function declarations>>
4836 <<kramers k global declarations>>
4840 <<kramers k structure>>=
4841 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
4842 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
4844 typedef struct kramers_param_struct {
4845 double D; /* diffusion rate (um/s) */
4852 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
4853 //kramers_x_func_t *xts; /* function returning transition x (nm) */
4854 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
4855 //kramers_E_func_t *E; /* function returning E (J) */
4856 //double *E_params; /* parametrize the energy landscape */
4857 //int n_E_params; /* length of E_params array */
4861 <<kramers k function declaration>>=
4862 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
4863 extern double kramers_k(double F, environment_t *env, void *kramers_params);
4865 <<kramers k function>>=
4866 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
4868 static char input[160]={0};
4869 static char output[80]={0};
4870 /* setup the environment */
4871 fprintf(in, "F = %g; T = %g\n", F, T);
4872 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4873 string_eval(in, out, input, 80, output);
4874 return atof(output);
4877 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
4879 static char input[160]={0};
4880 static char output[80]={0};
4881 /* setup the environment */
4882 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
4883 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4884 string_eval(in, out, input, 80, output);
4885 return atof(output);
4888 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
4890 kramers_param_t *p = (kramers_param_t *) kramers_params;
4891 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
4894 double kramers_k(double F, environment_t *env, void *kramers_params)
4895 { /* Returns the rate constant k in frac pop per s.
4896 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4897 * uses global k_B in J/K */
4898 kramers_param_t *p = (kramers_param_t *) kramers_params;
4899 double kbT, xc, xts, lc, lts, Eb;
4900 assert(F >= 0); assert(env->T > 0);
4903 assert(p->xc != NULL); assert(p->xts != NULL);
4904 assert(p->ddEdxx != NULL); assert(p->E != NULL);
4905 assert(p->in != NULL); assert(p->out != NULL);
4907 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
4908 if (xc == -1.0) return -1.0; /* error (Too much force) */
4909 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
4910 if (xc == -1.0) return -1.0; /* error (Too much force) */
4911 lc = sqrt(2.0*M_PI*kbT /
4912 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
4913 lts = sqrt(-2.0*M_PI*kbT/
4914 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
4915 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
4916 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
4917 //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));
4918 return p->D/(lc*lts) * exp(-Eb/kbT);
4922 <<kramers k structure create/destroy function declarations>>=
4923 void *string_create_kramers_param_t(char **param_strings);
4924 void destroy_kramers_param_t(void *p);
4927 <<kramers k structure create/destroy functions>>=
4928 kramers_param_t *create_kramers_param_t(double D,
4929 char *xc_fn, char *xts_fn,
4930 char *E_fn, char *ddEdxx_fn,
4931 char *global_define_string)
4932 // kramers_x_func_t *xc, kramers_x_func_t *xts,
4933 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
4934 // double *E_params, int n_E_params)
4936 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
4937 assert(ret != NULL);
4938 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4939 assert(ret->xc != NULL);
4940 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4941 assert(ret->xts != NULL);
4942 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
4943 assert(ret->E != NULL);
4944 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
4945 assert(ret->ddEdxx != NULL);
4947 strcpy(ret->xc, xc_fn);
4948 strcpy(ret->xts, xts_fn);
4949 strcpy(ret->E, E_fn);
4950 strcpy(ret->ddEdxx, ddEdxx_fn);
4951 //ret->E_params = E_params;
4952 //ret->n_E_params = n_E_params;
4953 string_eval_setup(&ret->in, &ret->out);
4954 fprintf(ret->in, "kB = %g\n", k_B);
4955 fprintf(ret->in, "%s\n", global_define_string);
4959 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
4960 void *string_create_kramers_param_t(char **param_strings)
4962 return create_kramers_param_t(atof(param_strings[0]),
4970 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
4972 kramers_param_t *pk = (kramers_param_t *)p;
4974 string_eval_teardown(&pk->in, &pk->out);
4976 // free(pk->E_params);
4982 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
4983 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.
4984 We follow \citet{shillcock98} and use
4986 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
4987 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
4990 where TODO, variable meanings.%\citep{shillcock98}.
4991 The first and second derivatives of $E(F,E_b,x,x_s)$ are
4993 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\\
4994 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
4996 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
4997 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
4998 The roots are therefor at
5000 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
5001 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
5002 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
5005 <<kramers k global declarations>>=
5006 extern int num_kramers_params;
5007 extern char *kramers_param_descriptions[];
5008 extern char kramers_param_string[];
5011 <<kramers k globals>>=
5012 int num_kramers_params = 6;
5013 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"};
5014 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)}";
5016 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
5017 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
5018 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
5019 Returning an $x_c$ position of $-1\U{m}$ lets the simulation know that the rate is poorly defined for that force, and the timestepping algorithm in Section \ref{sec.adaptive_dt} reduces it's timestep appropriately.
5020 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
5021 It works so long as [[val_1]] is not [[false]].
5023 <<kramers k model getopt>>=
5024 {"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}
5027 \subsection{Kramer's model (natural cubic splines)}
5028 \label{sec.kramers_integ}
5030 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
5031 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
5032 \citet{schlierf06} Eqn.~2)
5034 \frac{1}{k} = \frac{1}{D}
5035 \int_{x_\text{min}}^{x_\text{max}}
5036 e^{\frac{U_F(x)}{k_B T}}
5037 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5040 where $D$ is the diffusion constant,
5041 $U_F(x)$ is the free energy along the unfolding coordinate, and
5042 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
5043 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
5045 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
5046 interpolating between them with cubic splines.
5047 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
5049 <<kramers integ k functions>>=
5050 <<spline functions>>
5051 <<kramers integ A k functions>>
5052 <<kramers integ B k functions>>
5053 <<kramers integ k function>>
5054 <<kramers integ k structure create/destroy functions>>
5057 <<kramers integ k declarations>>=
5058 <<kramers integ k function declaration>>
5059 <<kramers integ k structure create/destroy function declarations>>
5060 <<kramers integ k global declarations>>
5064 <<kramers integ k structure>>=
5065 typedef double func_t(double x, void *params);
5066 typedef struct kramers_integ_param_struct {
5067 double D; /* diffusion rate (um/s) */
5068 double x_min; /* integration bounds */
5070 func_t *E; /* function returning E (J) */
5071 void *E_params; /* parametrize the energy landscape */
5072 destroy_data_func_t *destroy_E_params;
5074 double F; /* for passing into GSL evaluations */
5076 } kramers_integ_param_t;
5079 <<spline param structure>>=
5080 typedef struct spline_param_struct {
5081 int n_knots; /* natural cubic spline knots for U(x) */
5084 gsl_spline *spline; /* GSL spline parameters */
5085 gsl_interp_accel *acc; /* GSL spline acceleration data */
5089 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
5091 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
5093 <<kramers integ A k functions>>=
5094 double kramers_integ_k_integrandA(double x, void *params)
5096 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5097 double U = (*p->E)(x, p->E_params);
5098 double Fx = p->F * x;
5099 double kbT = k_B * p->env->T;
5100 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
5101 return exp(-(U-Fx)/kbT);
5103 double kramers_integ_k_integralA(double x, void *params)
5105 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5107 double result, abserr;
5108 size_t neval; /* for qng */
5109 static gsl_integration_workspace *w = NULL;
5110 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5111 f.function = &kramers_integ_k_integrandA;
5113 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5114 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5115 //fprintf(stderr, "integralA = %g\n", result);
5121 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
5122 \int_{x_\text{min}}^{x_\text{max}}
5123 e^{\frac{U_F(x)}{k_B T}}
5124 \text{Integral}_A(x)
5127 <<kramers integ B k functions>>=
5128 double kramers_integ_k_integrandB(double x, void *params)
5130 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5131 double U = (*p->E)(x, p->E_params);
5132 double Fx = p->F * x;
5133 double kbT = k_B * p->env->T;
5134 double intA = kramers_integ_k_integralA(x, params);
5135 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
5136 return exp((U-Fx)/kbT)*intA;
5138 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
5140 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5142 double result, abserr;
5143 size_t neval; /* for qng */
5144 static gsl_integration_workspace *w = NULL;
5145 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
5146 f.function = &kramers_integ_k_integrandB;
5150 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
5151 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
5152 //fprintf(stderr, "integralB = %g\n", result);
5157 With the integrals taken care of, evaluating $k$ is simply
5159 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
5161 <<kramers integ k function declaration>>=
5162 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
5164 <<kramers integ k function>>=
5165 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
5166 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
5167 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
5168 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
5172 <<kramers integ k structure create/destroy function declarations>>=
5173 void *string_create_kramers_integ_param_t(char **param_strings);
5174 void destroy_kramers_integ_param_t(void *p);
5177 <<kramers integ k structure create/destroy functions>>=
5178 kramers_integ_param_t *create_kramers_integ_param_t(double D,
5179 double xmin, double xmax, func_t *E,
5181 destroy_data_func_t *destroy_E_params)
5183 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
5184 assert(ret != NULL);
5189 ret->E_params = E_params;
5190 ret->destroy_E_params = destroy_E_params;
5194 void *string_create_kramers_integ_param_t(char **param_strings)
5196 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
5197 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
5198 void *E_params = string_create_spline_param_t(param_strings+1);
5199 return create_kramers_integ_param_t(atof(param_strings[0]),
5200 atof(param_strings[2]),
5201 atof(param_strings[3]),
5202 &spline_eval, E_params,
5203 destroy_spline_param_t);
5206 void destroy_kramers_integ_param_t(void *params)
5208 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5211 (*p->destroy_E_params)(p->E_params);
5217 Finally we have the GSL spline wrappers:
5218 <<spline functions>>=
5220 <<create/destroy spline>>
5223 <<spline function>>=
5224 double spline_eval(double x, void *spline_params)
5226 spline_param_t *p = (spline_param_t *)(spline_params);
5227 return gsl_spline_eval(p->spline, x, p->acc);
5231 <<create/destroy spline>>=
5232 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
5234 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
5235 assert(ret != NULL);
5236 ret->n_knots = n_knots;
5239 ret->acc = gsl_interp_accel_alloc();
5240 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
5241 gsl_spline_init(ret->spline, x, y, n_knots);
5245 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
5246 void *string_create_spline_param_t(char **param_strings)
5250 double *x=NULL, *y=NULL;
5251 /* break into ordered pairs */
5252 parse_list_string(param_strings[0], SEP, '(', ')',
5253 &num_knots, &knot_coords);
5254 x = (double *)malloc(sizeof(double)*num_knots);
5256 y = (double *)malloc(sizeof(double)*num_knots);
5258 for (i=0; i<num_knots; i++) {
5259 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
5260 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
5263 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
5265 free_parsed_list(num_knots, knot_coords);
5266 return create_spline_param_t(num_knots, x, y);
5269 void destroy_spline_param_t(void *params)
5271 spline_param_t *p = (spline_param_t *)params;
5274 gsl_spline_free(p->spline);
5276 gsl_interp_accel_free(p->acc);
5286 <<kramers integ k global declarations>>=
5287 extern int num_kramers_integ_params;
5288 extern char *kramers_integ_param_descriptions[];
5289 extern char kramers_integ_param_string[];
5292 <<kramers integ k globals>>=
5293 int num_kramers_integ_params = 4;
5294 char *kramers_integ_param_descriptions[] = {
5295 "diffusion D, m^2 / s",
5296 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
5297 "starting integration bound x_min, m",
5298 "ending integration bound x_max, m"};
5299 //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";
5300 //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";
5301 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";
5302 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
5303 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
5304 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
5307 <<kramers integ k model getopt>>=
5308 {"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}
5310 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
5311 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5313 \subsection{Unfolding model implementation}
5317 <<k model includes>>
5318 #include "k_model.h"
5319 <<k model internal definitions>>
5321 <<k model functions>>
5324 <<k model includes>>=
5325 #include <assert.h> /* assert() */
5326 #include <stdlib.h> /* NULL, malloc() */
5327 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
5328 #include <stdio.h> /* fprintf(), stdout */
5329 #include <string.h> /* strlen(), strcpy() */
5330 #include <gsl/gsl_integration.h>
5331 #include <gsl/gsl_spline.h>
5336 <<k model internal definitions>>=
5337 <<const k structure>>
5338 <<bell k structure>>
5339 <<kbell k structure>>
5340 <<kramers k structure>>
5341 <<kramers integ k structure>>
5342 <<spline param structure>>
5345 <<k model globals>>=
5350 <<kramers k globals>>
5351 <<kramers integ k globals>>
5354 <<k model functions>>=
5355 <<null k functions>>
5356 <<const k functions>>
5357 <<bell k functions>>
5358 <<kbell k functions>>
5359 <<kramers k functions>>
5360 <<kramers integ k functions>>
5363 \subsection{Unfolding model unit tests}
5365 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
5366 <<check-k-model.c>>=
5367 <<k model check includes>>
5368 <<check relative error>>
5370 <<k model test suite>>
5371 <<main check program>>
5374 <<k model check includes>>=
5375 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
5376 #include <stdio.h> /* sprintf() */
5377 #include <assert.h> /* assert() */
5378 #include <math.h> /* exp() */
5381 #include "k_model.h"
5384 <<k model test suite>>=
5387 <<k model suite definition>>
5390 <<k model suite definition>>=
5391 Suite *test_suite (void)
5393 Suite *s = suite_create ("k model");
5394 <<const k test case defs>>
5395 <<bell k test case defs>>
5397 <<const k test case adds>>
5398 <<bell k test case adds>>
5403 \subsubsection{Constant}
5405 <<const k test case defs>>=
5406 TCase *tc_const_k = tcase_create("Constant k");
5409 <<const k test case adds>>=
5410 tcase_add_test(tc_const_k, test_const_k_create_destroy);
5411 tcase_add_test(tc_const_k, test_const_k_over_range);
5412 suite_add_tcase(s, tc_const_k);
5417 START_TEST(test_const_k_create_destroy)
5419 char *knot[] = {"1","2","3","4","5","6"};
5420 char *params[] = {knot[0]};
5423 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5424 params[0] = knot[i];
5425 p = string_create_const_k_param_t(params);
5426 destroy_const_k_param_t(p);
5431 START_TEST(test_const_k_over_range)
5433 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
5434 char *knot[] = {"1","2","3","4","5","6"};
5435 char *params[] = {knot[0]};
5438 char param_string[80];
5441 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5442 params[0] = knot[i];
5443 p = string_create_const_k_param_t(params);
5444 for ( F=Fm; F<FM; F+=dF ) {
5445 fail_unless(const_k(F, &env, p)==atof(knot[i]),
5446 "const_k(%g, %g, &{%s}) = %g != %s",
5447 F, T, knot[i], const_k(F, &env, p), knot[i]);
5449 destroy_const_k_param_t(p);
5455 \subsubsection{Bell}
5457 <<bell k test case defs>>=
5458 TCase *tc_bell = tcase_create("Bell's k");
5461 <<bell k test case adds>>=
5462 tcase_add_test(tc_bell, test_bell_k_create_destroy);
5463 tcase_add_test(tc_bell, test_bell_k_at_zero);
5464 tcase_add_test(tc_bell, test_bell_k_at_one);
5465 suite_add_tcase(s, tc_bell);
5469 START_TEST(test_bell_k_create_destroy)
5471 char *knot[] = {"1","2","3","4","5","6"};
5473 char *params[] = {knot[0], dx};
5476 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5477 params[0] = knot[i];
5478 p = string_create_bell_param_t(params);
5479 destroy_bell_param_t(p);
5484 START_TEST(test_bell_k_at_zero)
5486 double F=0.0, T=300.0;
5487 char *knot[] = {"1","2","3","4","5","6"};
5489 char *params[] = {knot[0], dx};
5492 char param_string[80];
5495 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5496 params[0] = knot[i];
5497 p = string_create_bell_param_t(params);
5498 fail_unless(bell_k(F, &env, p)==atof(knot[i]),
5499 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
5500 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
5501 destroy_bell_param_t(p);
5506 START_TEST(test_bell_k_at_one)
5509 char *knot[] = {"1","2","3","4","5","6"};
5511 char *params[] = {knot[0], dx};
5512 double F= k_B*T/atof(dx);
5515 char param_string[80];
5518 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5519 params[0] = knot[i];
5520 p = string_create_bell_param_t(params);
5521 CHECK_ERR(1e-6, atof(knot[i])*exp(1.0), bell_k(F, &env, p));
5522 //fail_unless(bell_k(F, &env, p)==atof(knot[i])*exp(1.0),
5523 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
5524 // F, T, knot[i], dx, bell_k(F, &env, p), atof(knot[i])*exp(1.0));
5525 destroy_bell_param_t(p);
5531 <<kramers k tests>>=
5534 <<kramers k test case def>>=
5537 <<kramers k test case add>>=
5540 <<k model function tests>>=
5541 <<check relative error>>
5543 START_TEST(test_something)
5545 double k=5, x=3, last_x=2.0, F;
5546 one_dim_fn_t *handlers[] = {&hooke};
5547 void *data[] = {&k};
5548 double xi[] = {0.0};
5550 int new_active_groups = 1;
5551 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
5552 fail_unless(F = k*x, NULL);
5557 \subsection{Utilities}
5559 The unfolding models can get complicated, and are sometimes hard to explain to others.
5560 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
5561 the overhead of having to understand a full simulation.
5562 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
5563 or special mode, where the operation depends on the specific model selected.
5565 <<k-model-utils.c>>=
5567 <<k model utility includes>>
5568 <<k model utility definitions>>
5569 <<k model utility getopt functions>>
5570 <<k model utility multi model E>>
5571 <<k model utility main>>
5574 <<k model utility main>>=
5575 int main(int argc, char **argv)
5577 <<k model getopt array definition>>
5578 k_model_getopt_t *model = NULL;
5583 int num_param_args; /* for INIT_MODEL() */
5584 char **param_args; /* for INIT_MODEL() */
5586 double special_xmin;
5587 double special_xmax;
5588 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
5589 &Fmax, &special_xmin, &special_xmax, &flags);
5590 if (flags & VFLAG) {
5591 printf("#initializing model %s with parameters %s\n", model->name, model->params);
5593 INIT_MODEL("utility", model, model->params, params);
5596 if (model->k == NULL) {
5597 printf("No k function for model %s\n", model->name);
5601 printf("Fmax = %g <= 0\n", Fmax);
5607 printf("#F (N)\tk (%% pop. per s)\n");
5608 for (i=0; i<=N; i++) {
5609 F = Fmax*i/(double)N;
5610 k = (*model->k)(F, &env, params);
5612 printf("%g\t%g\n", F, k);
5617 if (model->k == NULL || model->k == &bell_k) {
5618 printf("No special mode for model %s\n", model->name);
5621 if (special_xmax <= special_xmin) {
5622 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
5626 double Xrng=(special_xmax-special_xmin), x, E;
5628 printf("#x (m)\tE (J)\n");
5629 for (i=0; i<=N; i++) {
5630 x = special_xmin + Xrng*i/(double)N;
5631 E = multi_model_E(model->k, params, &env, x);
5632 printf("%g\t%g\n", x, E);
5639 if (model->destructor)
5640 (*model->destructor)(params);
5645 <<k model utility multi model E>>=
5646 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
5648 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
5650 if (k == kramers_integ_k)
5651 E = (*p->E)(x, p->E_params);
5653 E = kramers_E(0, env, model_params, x);
5659 <<k model utility includes>>=
5662 #include <string.h> /* strlen() */
5663 #include <assert.h> /* assert() */
5666 #include "string_eval.h"
5667 #include "k_model.h"
5670 <<k model utility definitions>>=
5671 <<version definition>>
5672 #define VFLAG 1 /* verbose */
5673 enum mode_t {M_K_OF_F, M_SPECIAL};
5674 <<string comparison definition and globals>>
5675 <<k model getopt definitions>>
5676 <<initialize model definition>>
5677 <<kramers k structure>>
5678 <<kramers integ k structure>>
5682 <<k model utility getopt functions>>=
5684 <<k model utility help>>
5685 <<k model utility get options>>
5688 <<k model utility help>>=
5689 void help(char *prog_name,
5691 int n_k_models, k_model_getopt_t *k_models,
5692 int k_model, double Fmax, double special_xmin, double special_xmax)
5695 printf("usage: %s [options]\n", prog_name);
5696 printf("Version %s\n\n", VERSION);
5697 printf("A utility for understanding the available unfolding models\n\n");
5698 printf("Environmental options:\n");
5699 printf("-T T\tTemperature T (currently %g K)\n", env->T);
5700 printf("-C T\tYou can also set the temperature T in Celsius\n");
5701 printf("Model options:\n");
5702 printf("-k model\tTransition rate model (currently %s)\n",
5703 k_models[k_model].name);
5704 printf("-K args\tTransition rate model argument string (currently %s)\n",
5705 k_models[k_model].params);
5706 printf("There are two output modes. In standard mode, k(F) is printed\n");
5707 printf("For example:\n");
5708 printf(" #Force (N)\tk (% pop. per s)\n");
5709 printf(" 123.456\t7.89\n");
5711 printf("In special mode, the output depends on the model.\n");
5712 printf("For models defining an energy function E(x), that function is printed\n");
5713 printf(" #Distance (m)\tE (J)\n");
5714 printf(" 0.0012\t0.0034\n");
5716 printf("-m\tChange output to standard mode\n");
5717 printf("-M\tChange output to special mode\n");
5718 printf("-F\tSet the maximum F value for the standard mode k(F) output (currently %g)\n", Fmax);
5719 printf("-x\tSet the minimum x value for the special mode E(x) output (currently %g)\n", special_xmin);
5720 printf("-X\tSet the maximum x value for the special mode E(x) output (currently %g)\n", special_xmax);
5721 printf("-V\tChange output to verbose mode\n");
5722 printf("-h\tPrint this help and exit\n");
5724 printf("Unfolding rate models:\n");
5725 for (i=0; i<n_k_models; i++) {
5726 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
5727 for (j=0; j < k_models[i].num_params ; j++)
5728 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
5729 printf("\t\tdefaults: %s\n", k_models[i].params);
5734 <<k model utility get options>>=
5735 void get_options(int argc, char **argv, environment_t *env,
5736 int n_k_models, k_model_getopt_t *k_models,
5737 enum mode_t *mode, k_model_getopt_t **model,
5738 double *Fmax, double *special_xmin, double *special_xmax,
5739 unsigned int *flags)
5741 char *prog_name = NULL;
5742 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
5744 extern char *optarg;
5745 extern int optind, optopt, opterr;
5749 /* setup defaults */
5751 prog_name = argv[0];
5752 env->T = 300.0; /* K */
5757 *special_xmax = 1e-8;
5758 *special_xmin = 0.0;
5760 while ((c=getopt(argc, argv, options)) != -1) {
5762 case 'T': env->T = atof(optarg); break;
5763 case 'C': env->T = atof(optarg)+273.15; break;
5765 k_model = index_k_model(n_k_models, k_models, optarg);
5766 *model = k_models+k_model;
5769 k_models[k_model].params = optarg;
5771 case 'm': *mode = M_K_OF_F; break;
5772 case 'M': *mode = M_SPECIAL; break;
5773 case 'F': *Fmax = atof(optarg); break;
5774 case 'x': *special_xmin = atof(optarg); break;
5775 case 'X': *special_xmax = atof(optarg); break;
5776 case 'V': *flags |= VFLAG; break;
5778 fprintf(stderr, "unrecognized option '%c'\n", optopt);
5779 /* fall through to default case */
5781 help(prog_name, env, n_k_models, k_models, k_model, *Fmax, *special_xmin, *special_xmax);
5790 \section{\LaTeX\ documentation}
5792 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
5793 The comment blocks are extracted (with nicely formatted code blocks), using
5794 <<latex makefile lines>>=
5795 $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5796 noweave -latex -index -delay $< > $@
5797 $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR)
5801 <<latex makefile lines>>=
5802 $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
5804 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5805 $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
5806 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5807 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5808 mv $(BUILD_DIR)/sawsim.pdf $@
5810 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
5811 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
5812 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
5814 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
5815 <<latex makefile lines>>=
5817 rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
5818 $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
5819 $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.toc \
5820 $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib
5821 clean_latex : semi-clean_latex
5822 rm -f $(DOC_DIR)/sawsim.pdf
5827 [[make]] is a common utility on *nix systems for managing dependencies while building software.
5828 In this case, we have several \emph{targets} that we'd like to build:
5829 the main simulation program \prog;
5830 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
5831 the documentation [[sawsim.pdf]];
5832 and the [[Makefile]] itself.
5833 Besides the generated files, there is the utility target
5834 [[clean]] for removing all generated files except [[Makefile]].
5836 % [[dist]] for generating a distributable tar file.
5838 Extract the makefile with
5839 `[[notangle -Rmakefile src/sawsim.nw | sed 's/ /\t/' > Makefile]]'.
5840 The sed is needed because notangle replaces tabs with eight spaces,
5841 but [[make]] needs tabs.
5843 # Decide what directories to put things in
5848 # Note: Cannot use, for example, `./build', because make eats the `./'
5849 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5851 # Modules (X.c, X.h) defined in the noweb file
5854 # Modules defined outside the noweb file
5855 FREE_SAWSIM_MODS = interp tavl
5857 <<list module makefile lines>>
5858 <<tension balance module makefile lines>>
5859 <<tension model module makefile lines>>
5860 <<k model module makefile lines>>
5861 <<parse module makefile lines>>
5862 <<string eval module makefile lines>>
5863 <<accel k module makefile lines>>
5865 SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS)
5867 # Everything needed for sawsim under one roof. sawsim.c must come first
5868 SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \
5869 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h)
5870 # Libraries needed to compile sawsim
5871 LIBS = gsl gslcblas m
5872 CHECK_LIBS = gsl gslcblas m check
5873 # The non-check binaries generated
5874 BINS = sawsim tension_model_utils k_model_utils sawsim_profile
5877 # Define the major targets
5878 all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
5880 view : $(DOC_DIR)/sawsim.pdf
5882 profile : $(BIN_DIR)/sawsim_profile | $(BIN_DIR)
5883 $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \
5884 -mnull -Mwlc -A0.39e-9,28e-9 -F8
5885 gprof $(BIN_DIR)/sawsim_profile gmon.out > $@
5886 check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim
5887 $(SHELL) -e -c 'for B in $^; do ./$$B; done'
5888 clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \
5889 clean_tension_model_utils \
5890 clean_k_model_utils clean_latex clean_check_sawsim
5891 rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \
5892 $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \
5893 $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
5894 $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
5895 $(BUILD_DIR)/global.h ./gmon.out
5896 $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
5898 # Various builds of sawsim
5899 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
5900 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5901 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
5902 gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5903 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
5904 gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5906 # Create the directories
5907 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
5910 # Copy over the source external to sawsim
5911 # Note: Cannot use, for example, `./build', because make eats the `./'
5912 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5914 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
5918 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
5922 ## Basic source generated with noweb
5923 # The central files sawsim.c and global.h...
5924 $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5926 $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5927 notangle -Rglobal.h $< > $@
5928 # ... and the modules
5929 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5930 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5931 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5932 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5933 $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5934 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5935 # Note: I use `_' as a space in my file names, but noweb doesn't like
5936 # that so we use `-' to name the noweb roots and substitute here.
5938 ## Building the unit-test programs
5940 $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
5941 notangle -Rchecks $< > $@
5942 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
5943 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
5944 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR)
5945 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
5946 clean_check_sawsim :
5947 rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
5948 # ... and the modules
5950 $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
5951 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \
5952 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \
5953 $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
5954 $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
5955 $$(BUILD_DIR)/global.h | $(BIN_DIR)
5956 gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
5957 $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
5959 # todo: clean up the required modules to
5960 $(CHECK_BINS:%=clean_%) :
5961 rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
5963 # Cleaning up the modules
5965 $(SAWSIM_MODS:%=clean_%) :
5966 rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h)
5968 <<tension model utils makefile lines>>
5969 <<k model utils makefile lines>>
5970 <<latex makefile lines>>
5972 Makefile : $(SRC_DIR)/sawsim.nw
5973 notangle -Rmakefile $< | sed 's/ /\t/' > $@
5975 This is fairly self-explanatory. We compile a dynamically linked
5976 version ([[sawsim]]) and a statically linked version
5977 ([[sawsim_static]]). The static version is about 9 times as big, but
5978 it runs without needing dynamic GSL libraries to link to, so it's a
5979 better format for distributing to a cluster for parallel evaluation.
5983 \subsection{Hookean springs in series}
5984 \label{sec.math_hooke}
5987 The effective spring constant for $n$ springs $k_i$ in series is given by
5989 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5995 F &= k_i x_i &&\forall i \le n \\
5996 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
5997 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
5998 F &= k_1 x_1 = k_\text{eff} x \\
5999 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
6000 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
6005 \addcontentsline{toc}{section}{References}
6006 \bibliography{sawsim}