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 definition>>=
153 #define VERSION "0.6"
159 sawsim - program for simulating single molecule mechanical unfolding.
160 Copyright (C) 2008-2009, William Trevor King
162 This program is free software; you can redistribute it and/or
163 modify it under the terms of the GNU General Public License as
164 published by the Free Software Foundation; either version 3 of the
165 License, or (at your option) any later version.
167 This program is distributed in the hope that it will be useful, but
168 WITHOUT ANY WARRANTY; without even the implied warranty of
169 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
170 See the GNU General Public License for more details.
172 You should have received a copy of the GNU General Public License
173 along with this program; if not, write to the Free Software
174 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
177 The author may be contacted at <wking@drexel.edu> on the Internet, or
178 write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
179 Philadelphia PA 19104, USA.
192 The unfolding system is stored as a chain of \emph{domains} (Figure
193 \ref{fig.domain_chain}). Domains can be folded, globular protein
194 domains, unfolded protein linkers, AFM cantilevers, or other
195 stretchable system components. The system is modeled as graph of
196 possible domain states with transitions between them (Figure
197 \ref{fig.domain_states}).
201 \subfloat[][]{\label{fig.domain_chain}
202 \begin{tikzpicture}[->,node distance=2cm,font=\footnotesize]
203 \tikzstyle{every state}=[fill,draw=red!50,very thick,fill=red!20]
204 \node[state] (A) {domain 1};
205 \node[state] (B) [below of=A] {domain 2};
206 \node[state] (C) [below of=B] {.~.~.};
207 \node[state] (D) [below of=C] {domain $N$};
208 \node (S) [below of=D] {Surface};
209 \node (E) [above of=A] {};
211 \path[-] (A) edge (B)
212 (B) edge node [right] {Tension} (C)
215 \path[->,green] (A) edge node [right,black] {Extension} (E);
219 \subfloat[][]{\label{fig.domain_states}
220 \begin{tikzpicture}[->,node distance=2.5cm,shorten <=1pt,shorten >=1pt,font=\footnotesize]
221 \tikzstyle{every state}=[fill,draw=blue!50,very thick,fill=blue!20]
222 \node[state] (A) {cantilever};
223 \node[state] (C) [below of=A] {transition};
224 \node[state] (B) [left of=C] {folded};
225 \node[state] (D) [right of=C] {unfolded};
227 \path (B) edge [bend left] node [above] {$k_1$} (C)
228 (C) edge [bend left] node [below] {$k_1'$} (B)
229 edge [bend left] node [above] {$k_2$} (D)
230 (D) edge [bend left] node [below] {$k_2'$} (C);
233 \caption{\subref{fig.domain_chain} Extending a chain of domains. One end
234 of the chain is fixed, while the other is extended at a constant
235 speed. The domains are coupled with rigid linkers, so the domains
236 themselves must stretch to accomodate the extension.
237 \subref{fig.domain_states} Each domain exists in a discrete state. At
238 each timestep, it may transition into another state following a
239 user-defined state matrix such as this one, showing a metastable
240 transition state and an explicit ``cantilever'' domain.}
244 Each domain contributes to the total tension, which depends on the
245 chain extension. The particular model for the tension as a function
246 of extension is handled generally with function pointers. So far the
247 following models have been implemented:
249 \item Null (Appendix \ref{sec.null_tension}),
250 \item Constant (Appendix \ref{sec.const_tension}),
251 \item Hookean spring (Appendix \ref{sec.hooke}),
252 \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
253 \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
256 The tension model and parameters used for a particular domain depend
257 on whether the domain's current state. The transition rates between
258 states are also handled generally with function pointers, with
261 \item Null (Appendix \ref{sec.null_k}),
262 \item Bell's model (Appendix \ref{sec.bell}),
263 \item Spring constant corrected Bell model (Appendix \ref{sec.kbell}),
264 \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
265 \item Kramers' saddle point model (Appendix \ref{sec.kramers}).
268 The unfolding of the chain domains is modeled in two stages. First
269 the tension of the chain is calculated. Then the tension (assumed to
270 be constant over the length of the chain, see Section
271 \ref{sec.timescales}), is applied to each domain, and used to
272 calculate the probability of that domain changing state in the next
273 timestep $dt$. Then the time is stepped forward, and the process is
274 repeated. The simulation is complete when a pre-set time limit
275 [[t_max]] is reached or a pre-selected state [[stop_state]] is empty.
277 int main(int argc, char **argv)
279 <<tension model getopt array definition>>
280 <<k model getopt array definition>>
282 /* variables initialized in get_options() */
283 list_t *state_list=NULL, *transition_list=NULL;
284 environment_t env={0};
285 double P, t_max, dt_max, v, x_step;
286 state_t *stop_state=NULL;
288 /* variables used in the simulation loop */
289 double t=0, x=0, dt, F; /* time, distance, timestep, force */
290 int transition=1; /* boolean marking a transition for tension optimization */
292 get_options(argc, argv, &P, &t_max, &dt_max, &v, &x_step, &stop_state, &env,
293 NUM_TENSION_MODELS, tension_models, NUM_K_MODELS, k_models,
294 &state_list, &transition_list);
297 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Distance (m)\tForce (N)\tStiffness (N/m)\n");
298 if (flags & FLAG_OUTPUT_TRANSITION_FORCES) printf("#Force (N)\tInitial state\tFinal state\n");
299 while ((t_max > 0 && t <= t_max) || (stop_state != NULL && stop_state->num_domains > 0)) {
300 F = find_tension(state_list, &env, x, transition);
302 dt = determine_dt(state_list, transition_list, &env, x, P, dt_max, v);
304 dt = determine_dt(state_list, transition_list, &env, x, P, dt_max, 0);
305 transition=random_transitions(transition_list, F, dt, &env);
306 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness);
311 if (dt == dt_max) { /* step completed */
314 } else { /* still working on this step */
319 destroy_state_list(state_list);
320 destroy_transition_list(transition_list);
324 @ The meat of the simulation is bundled into the three functions
325 [[find_tension]], [[determine_dt]], and [[random_transitions]].
326 [[find_tension]] is discussed in Section \ref{sec.find_tension},
327 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
328 [[random_transitions]] in Section \ref{sec.transition_rate}.
330 The stretched distance is extended in one of two modes depending on
331 [[xstep]]. If [[xstep == 0]], then the pulling is continuous, at
332 least within the limits of the inherent discretization of a time
333 stepped approach. At any rate, the timestep is limited by the maximum
334 allowed timestep [[dt_max]] and the maximum allowed unfolding
335 probability in a given step [[P]]. In the second mode [[xstep > 0]],
336 and the end to end distance increases in discrete steps of that
337 length. The time between steps is chosen to maintain an average
338 unfolding velocity [[v]]. This approach is less work to simulate than
339 smooth pulling and also closer to the reality of many experiments, but
340 it is practically impossible to treat analytically. With both pulling
341 styles implemented, the effects of the discretization can be easily
342 calculated, bridging the gap between theory and experiment.
344 Environmental parameters important in determining reaction rates and
345 tensions (e.g.\ temperature, pH) are stored in a single structure to
346 facilitate extension to more complicated models in the future.
347 <<environment definition>>=
348 typedef struct environment_struct {
357 <<environment definition>>
358 <<one dimensional function definition>>
359 <<create/destroy definitions>>
361 #endif /* GLOBAL_H */
363 Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
364 included multiple times.
366 \section{Simulation functions}
368 <<simulation functions>>=
372 <<does domain transition>>
373 <<randomly transition domains>>
377 \label{sec.find_tension}
379 Because the stretched system may be made up of several parts (folded
380 domains, unfolded domains, spring-like cantilever, \ldots), we will
381 assign the domains to states with tension models and groups. The
382 models are used to determine the tension of all the domain in a given
383 state following some model (Hookean spring, WLC, \ldots). The domains
384 are assumed to be commutative, so ordering is ignored. The
385 interactions between the groups are assumed to be linear, but the
386 interactions between domains of the same group need not be. Each
387 model has a tension handler function, which gives the tension $F$
388 needed to stretch a given group of domains a distance $x$.
390 GROUPS ARE CURRENTLY DISABLED.
392 Groups represent collections of states which the model should treat as
393 a single entity. To understand the purpose of groups, consider a
394 system with two unfolded domain states (e.g.\ for two different
395 proteins) that were both modeled as WLCs. If both unfolded states had
396 the same persistence length, it would be wise to assign them to the
397 same group. This would reduce both the computational cost of
398 balancing the tension and the error associated with the inter-group
399 interaction linearity. Note that group numbers only matter
400 \emph{within} model classes, since grouping domains with seperate
401 models doesn't make sense.
403 <<find tension definitions>>=
404 #define NUM_TENSION_MODELS 5
408 <<tension handler array global declaration>>=
409 extern one_dim_fn_t *tension_handlers[];
412 <<tension handler array global>>=
413 one_dim_fn_t *tension_handlers[] = {
415 &const_tension_handler,
423 <<tension model global declarations>>=
424 <<tension handler array global declaration>>
427 <<tension handler types>>=
428 typedef struct tension_handler_data_struct {
429 list_t *group_tension_params; /* one item for each domain in the group */
432 } tension_handler_data_t;
436 After sorting the chain into separate groups $G_i$, with tension
437 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
438 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
439 \forall i,j$. Note that there may be several states within each
440 group. [[find_tension]] is basically a wrapper to feed proper input
441 into the [[tension_balance]] function. [[transition]] should be set
442 to 0 if there were no transitions since the previous call to
443 [[find_tension]] to support some optimizations that assume a small
444 increase in tension between steps (only true if no transition has
445 occured). While were messing with the tension handlers, we also grab
446 a measure of the current linker stiffness for the environmental
447 variable, which is needed by some of the unfolding rate models
448 (e.g.\ the linker-stiffness corrected Bell model, Appendix
451 double find_tension(list_t *state_list, environment_t *env, double x, int transition)
453 static int num_active_groups=0;
454 static one_dim_fn_t **per_group_handlers = NULL;
455 static void **per_group_data = NULL; /* array of pointers to tension_handler_data_t */
456 static double last_x;
457 tension_handler_data_t *tdata;
462 fprintf(stderr, "%s:\tfinding tension at distance %g%s\n", __FUNCTION__, x, transition ? " after a transition" : "");
465 if (transition != 0 || num_active_groups == 0) { /* setup statics */
466 /* get new statics, freeing the old ones if they aren't NULL */
467 get_active_groups(state_list, &num_active_groups, &per_group_handlers, &per_group_data);
469 /* fill in the group handlers and params */
471 fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_data, per_group_data[0]);
473 for (i=0; i<num_active_groups; i++) {
474 tdata = (tension_handler_data_t *) per_group_data[i];
476 fprintf(stderr, "%s:\tactive group %d of %d (data %p, param list %p)", __FUNCTION__, i, num_active_groups, tdata, tdata->group_tension_params);
477 list_print(stderr, tdata->group_tension_params, " parameters");
480 /* tdata->persist continues unchanged */
483 } /* else continue with the current statics */
485 F = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_data, last_x, x, &env->stiffness);
491 @ For the moment, we will restrict our group boundaries to a single
492 dimension, so $\sum_i x_i = x$, our total extension (it is this
493 restriction that causes the groups to interact linearly). We'll also
494 restrict our extensions to all be positive. With these restrictions,
495 the problem of balancing the tensions reduces to solving a system of
496 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
497 the number of active groups. In general this can be fairly
498 complicated, but without much loss of practicality we can restrict
499 ourselves to strictly monotonically increasing, non-negative tension
500 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
501 simpler. For example, it guarantees the existence of a unique, real
502 solution for finite forces. See Appendix \ref{sec.tension_balance}
503 for the tension balancing implementation.
505 Each group must define a parameter structure if the tension model is
507 a function to create the parameter structure,
508 a function to destroy the parameter structure, and
509 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as it's data argument.
510 The tension-specific parameter structures are contained in the domain
511 group field. See the Section \ref{sec.model_selection} for a
512 discussion on the command line framework. See the worm-like chain
513 implementation for an example (Section \ref{sec.wlc}).
515 The major limitation of our approach is that all of our tension models
516 are Markovian (memory-less), which is not really a problem for the
517 chains (see \citet{evans99} for WLC thermalization rate calculations).
519 \subsection{Transition rate}
520 \label{sec.transition_rate}
522 Each state transition is modeled as unimolecular, first order reaction
524 \text{State 1} \xrightarrow{k} \text{State 2}
526 With the rate constant $k$ defined by
528 \frac{d[\text{State 2}]}{dt} = -\frac{d[\text{State 1}]}{dt} = k [\text{State 1}].
530 So the probability of a given domain transitioning in a short time
536 <<does domain transition>>=
537 int domain_transitions(double F, double dt, environment_t *env, int num_domains, transition_t *transition)
538 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, number of 'reactant' domains, pointer to transition structure */
540 k = accel_k(transition->k, F, env, transition->k_params);
541 //(*transition->k)(F, env, domain->k_params);
542 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
543 return happens(k*dt*num_domains); /* N dice rolls for prob. k*dt event */
545 @ [[happens]] is a random decision making function defined in Appendix \ref{sec.utils}.
547 <<randomly transition domains>>=
548 int random_transitions(list_t *transition_list, double tension,
549 double dt, environment_t *env)
550 { /* returns 1 if there was a transition and 0 otherwise */
551 while (transition_list != NULL) {
552 if (T(transition_list)->initial_state->num_domains > 0
553 && domain_transitions(tension, dt, env, T(transition_list)->initial_state->num_domains, T(transition_list))) {
554 if (flags & FLAG_OUTPUT_TRANSITION_FORCES)
555 fprintf(stdout, "%g\t%s\t%s\n", tension, T(transition_list)->initial_state->name, T(transition_list)->final_state->name);
556 T(transition_list)->initial_state->num_domains--;
557 T(transition_list)->final_state->num_domains++;
558 return 1; /* our one domain has transitioned, stop looking for others */
560 transition_list = transition_list->next;
566 \subsection{Timescales}
567 \label{sec.timescales}
569 The simulation assumes that chain equilibration is instantaneous
570 relative to the simulation timestep $dt$, so that tension is uniform
571 along the chain. The quality of this assumption depends on your
572 particular chain. For example, a damped spring thermalizes on a
573 timescale of order $\tau = 1/\gamma$ where the damping ratio $\gamma
574 \equiv \omega_0(-\zeta \pm \sqrt{\zeta^2-1}$, the natural angular
575 frequency $omega_0 \equiv \sqrt{k/m}$, $\zeta \equiv b/2\sqrt{km}$,
576 and $b$ sets the drag $F_b = -bv$. For our cantilevers $\tau$ is
577 on the order of milliseconds, which is longer than a timestep.
578 However, the approximation is still reasonable, because there is
579 rarely unfolding during the cantilever return. You could, of course,
580 take cantilever, WLC, etc.\ response times into effect, but that
581 would significantly complicate a simulation that seems to work
582 well enough without bothering :p. For WLC and FJC relaxation times,
583 see the Appendix of \citet{evans99} and \citet{kroy07}.
585 Besides assuming our timestep is much greater than equilibration
586 times, we also force it to be short enough so that only one domain may
587 unfold in a given timestep. For an unfolding timescale of $dt_u =
588 1/k$ we adapt our timesteps to keep $dt \ll dt_u$, so the probability
589 of two domains unfolding in the same timestep is negligible. This
590 approach breaks down as the adaptive timestep scheme approaches $dt
591 \sim dt_u$, but $dt_u \sim 1\U{$\mu$s}$ for Ig27-like proteins
592 \citep{klimov00}, so this shouldn't be a problem. To reassure
593 yourself, you can ask the simulator to print the smallest timestep
594 that was used with TODO.
596 \subsection{Adaptive timesteps}
597 \label{sec.adaptive_dt}
599 We'd like to pick $dt$ so the probability of changing state
600 (e.g.\ unfolding) in the next timestep is small. If we don't adapt our
601 timestep, we also risk breaking our approximation that $F(x' \in [x,
602 x+v \cdot dt]) = F(x)$ or that only one domain changes state in a
603 given timestep. The simulation would have been accurate for
604 sufficiently small fixed timesteps, but adaptive timesteps will allow
605 us to move through low tension regions in fewer steps, leading to a
606 more efficient simulation.
608 Consider the two state folded $\rightarrow$ unfolded transition.
609 Because $F(x)$ is monotonically increasing between unfoldings,
610 excessively large timesteps will lead to erroneously small unfolding
611 rates (an thus, higher average unfolding force).
613 The actual adaptive timestep implementation is not particularly
614 interesting, since we are only required to reduce $dt$ to somewhere
615 below a set threshold, so I've removed it to Appendix
616 \ref{sec.adaptive_dt_implementation}.
622 The models provide the physics of the simulation, but the simulation
623 allows interchangeable models, and we are currently focusing on the
624 simulation itself, so we remove the actual model implementation to the
627 The tension models are defined in Appendix \ref{sec.tension_models}
628 and unfolding models are defined in Appendix \ref{sec.k_models}.
631 #define k_B 1.3806503e-23 /* J/K */
635 \section{Command line interface}
637 <<get option functions>>=
639 <<index tension model>>
644 \subsection{Model selection}
645 \label{sec.model_selection}
647 The main difficulty with the command line interface in \prog\ is
648 developing an intuitive method for accessing the possibly large number
649 of available models. We'll treat this generally by defining an array
650 of available models, containing their important parameters.
652 <<get options definitions>>=
653 <<model getopt definitions>>
656 <<create/destroy definitions>>=
657 typedef void *create_data_func_t(char **param_strings);
658 typedef void destroy_data_func_t(void *param_struct);
661 <<model getopt definitions>>=
662 <<tension model getopt definitions>>
663 <<k model getopt definitions>>
666 In order to choose models, we need a method of assembling a reaction
667 graph such as that in Figure \ref{fig.domain_states} and population
668 the graph with a set of domains. First we declare the domain states
671 [[-s <state-name>,<model>[,<group>] -S <model-params> [-N <initial-population]]] \\
675 [[-s unfolded,wlc,1 -S 1e-8,4e-10 -N 5]]
677 This creates the state named [[unfolded]], using the WLC model and one
678 as a group number (defaults to 0, see Section \ref{sec.find_tension}).
679 The tension parameters are then set to [[1e-8,4e-10]], which in the
680 case of the WLC are the contour and persistence lengths respectively.
681 Finally we populate the state by adding five domains to the state.
682 The name of the state is importent for identifying states later.
684 Once the domains have been created and populated, you can added
685 transition rates following
687 [[-k <initial-state>,<final-state>,<model> -K <model-params>]] \\
691 [[-k folded,unfolded,bell -K 3.3e-4,0.25e-9]]
693 This creates a reaction going from the [[folded]] state to the
694 [[unfolded]] state, with the rate constant given by the Bell model
695 (Appendix \ref{sec.bell}). The unfolding parameters are then set to
696 [[3.3e-4,0.25e-9]], which in the case of the Bell model are the
697 unforced rate and transition state distance respectively.
699 \subsubsection{Tension}
701 To access the various tension models from the command line, we define
702 a structure that contains all the neccessary information about the
704 <<tension model getopt definitions>>=
705 typedef struct tension_model_getopt_struct {
706 one_dim_fn_t *handler; /* fn implementing the model on a group */
707 char *name; /* name string identifying the model */
708 char *description; /* a brief description of the model */
709 int num_params; /* number of extra params for the model */
710 char **param_descriptions; /* descriptions of extra params */
711 char *params; /* default values for the extra params */
712 create_data_func_t *creator; /* param-string -> model-data-struct fn */
713 destroy_data_func_t *destructor; /* fn to free the model data structure */
714 } tension_model_getopt_t;
717 <<tension model getopt array definition>>=
718 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
719 <<null tension model getopt>>,
720 <<constant tension model getopt>>,
721 <<hooke tension model getopt>>,
722 <<worm-like chain tension model getopt>>,
723 <<freely-jointed chain tension model getopt>>
727 \subsubsection{Unfolding rate}
729 <<k model getopt definitions>>=
730 #define NUM_K_MODELS 6
732 typedef struct k_model_getopt_struct {
737 char **param_descriptions;
739 create_data_func_t *creator;
740 destroy_data_func_t *destructor;
744 <<k model getopt array definition>>=
745 k_model_getopt_t k_models[NUM_K_MODELS] = {
746 <<null k model getopt>>,
747 <<const k model getopt>>,
748 <<bell k model getopt>>,
749 <<kbell k model getopt>>,
750 <<kramers k model getopt>>,
751 <<kramers integ k model getopt>>
758 void help(char *prog_name, double P, double t_max, double dt_max, double v, double x_step,
759 state_t *stop_state, environment_t *env,
760 int n_tension_models, tension_model_getopt_t *tension_models,
761 int n_k_models, k_model_getopt_t *k_models,
762 int k_model, list_t *state_list)
767 if (state_list != NULL)
768 state = S(tail(state_list));
770 printf("usage: %s [options]\n", prog_name);
771 printf("Version %s\n\n", VERSION);
772 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
774 printf("Simulation options:\n");
776 printf("-q state\tStop simulation when this state empties (- for no limit, currently %s)\n", stop_state ? stop_state->name : "-" );
777 printf("-t tmax\tMaximum simulated time (-1 for no limit, currently %g s)\n", t_max);
778 printf("-d dt\tMaximum allowed timestep dt (currently %g s)\n", dt_max);
779 printf("-P P\tTarget probability for dt (currently %g)\n", P);
780 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
781 printf("-x dx\tPulling step size (currently %g m)\n", x_step);
782 printf("\tWhen dx = 0, the pulling is continuous.\n");
783 printf("\tWhen dx > 0, the pulling occurs in discrete steps.\n");
785 printf("Environmental options:\n");
786 printf("-T T\tTemperature T (currently %g K)\n", env->T);
787 printf("-C T\tYou can also set the temperature T in Celsius\n");
789 printf("State creation options:\n");
790 printf("-s name,model[[,group],params]\tCreate a new state.\n");
791 printf("Once you've set up the state, you may populate it with domains\n");
792 printf("-N n\tAdd n domains in the current state (%s)\n", state ? state->name : "-");
794 printf("Transition creation options:\n");
795 printf("-k initial_state,final_state,model[,params]\tCreate a new transition\n");
797 printf("To double check your reaction graph, you can produce graphviz dot output\n");
798 printf("describing the current states and transitions.\n");
799 printf("-D\tPrint dot description of states and transitions and exit.\n");
801 printf("Simulation output mode options:\n");
802 printf("There are two output modes. In standard mode, only the transition\n");
803 printf("events are printed. For example:\n");
804 printf(" #Force (N)\tInitial state\tFinal state\n");
805 printf(" 123.456e-12\tfolded\tunfolded\n");
807 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
808 printf(" #Distance (m)\tForce (N)\tStiffness (N/m)\n");
809 printf(" 0.001\t0.0023\t0.002\n");
811 printf("-V\tChange output to verbose mode.\n");
812 printf("-h\tPrint this help and exit.\n");
815 printf("Tension models:\n");
816 for (i=0; i<n_tension_models; i++) {
817 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
818 for (j=0; j < tension_models[i].num_params ; j++)
819 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
820 printf("\t\tdefaults: %s\n", tension_models[i].params);
822 printf("Unfolding rate models:\n");
823 for (i=0; i<n_k_models; i++) {
824 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
825 for (j=0; j < k_models[i].num_params ; j++)
826 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
827 printf("\t\tdefaults: %s\n", k_models[i].params);
830 printf("Examples:\n");
831 printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded. Stop once there are no more folded domains.\n");
832 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);
833 printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded. Stop after 0.25 seconds\n");
834 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);
838 \subsection{Get options}
841 void get_options(int argc, char **argv,
842 double *pP, double *pT_max, double *pDt_max, double *pV,
843 double *pX_step, state_t **pStop_state, environment_t *env,
844 int n_tension_models, tension_model_getopt_t *tension_models,
845 int n_k_models, k_model_getopt_t *k_models,
846 list_t **pState_list, list_t **pTransition_list)
848 char *prog_name = NULL;
849 char c, options[] = "q:t:d:P:v:x:T:C:s:N:k:DVh";
850 int tension_model=0, k_model=0, initial_state=-1, final_state=-1;
851 char *state_name=NULL;
852 int i, n, tension_group=0;
853 list_t *temp_list=NULL;
856 void *model_params; /* for INIT_MODEL() */
857 int num_param_args; /* for INIT_MODEL() */
858 char **param_args; /* for INIT_MODEL() */
860 extern int optind, optopt, opterr;
865 for (i=0; i<n_tension_models; i++) {
866 fprintf(stderr, "tension model %i %s:\t model's handler %p\n",
867 i, tension_models[i].name, tension_models[i].handler);
868 assert(tension_models[i].handler == tension_handlers[i]);
873 flags = FLAG_OUTPUT_TRANSITION_FORCES;
876 *pT_max = -1; /* s, -1 removes this limit */
877 *pDt_max = 0.001; /* s */
878 *pP = 1e-3; /* % pop per s */
879 *pV = 1e-6; /* m/s */
880 *pX_step = 0.0; /* m */
881 env->T = 300.0; /* K */
883 *pTransition_list = NULL;
885 while ((c=getopt(argc, argv, options)) != -1) {
888 if (STRMATCH(optarg, "-")) {
891 temp_list = head(*pState_list);
892 while (temp_list != NULL) {
893 if (STRMATCH(S(temp_list)->name, optarg)) {
894 *pStop_state = S(temp_list);
897 temp_list = temp_list->next;
901 case 't': *pT_max = atof(optarg); break;
902 case 'd': *pDt_max = atof(optarg); break;
903 case 'P': *pP = atof(optarg); break;
904 case 'v': *pV = atof(optarg); break;
905 case 'x': *pX_step = atof(optarg); break;
906 case 'T': env->T = atof(optarg); break;
907 case 'C': env->T = atof(optarg)+273.15; break;
909 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
910 assert(num_strings >= 2 && num_strings <= 4);
912 state_name = strings[0];
913 if (state_index(*pState_list, state_name) != -1) {
914 fprintf(stderr, "Invalid option '%s', repeated state name '%s'\n", optarg, state_name);
917 tension_model = index_tension_model(n_tension_models, tension_models, strings[1]);
918 if (num_strings == 4)
919 tension_group = atoi(strings[2]);
922 if (num_strings >= 3)
923 tension_models[tension_model].params = strings[num_strings-1];
925 fprintf(stderr, "%s:\t%s -> state %s, model %s, params %s (%s), group %g\n", __FUNCTION__, optarg, state_name, tension_models[tension_model].name, tension_models[tension_model].params, tension_group);
927 INIT_MODEL(state_name, &tension_models[tension_model], tension_models[tension_model].params, model_params);
928 push(pState_list, create_state(state_name,
929 tension_models[tension_model].handler,
932 tension_models[tension_model].destructor,
934 *pState_list = tail(*pState_list);
936 free_parsed_list(num_strings, strings);
941 fprintf(stderr, "%s:\tadding %d domains to %s\n", __FUNCTION__, n, S(*pState_list)->name);
943 S(*pState_list)->num_domains += n;
946 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
947 assert(num_strings >= 3 && num_strings <= 4);
948 initial_state = state_index(*pState_list, strings[0]);
949 final_state = state_index(*pState_list, strings[1]);
950 k_model = index_k_model(n_k_models, k_models, strings[2]);
952 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);
954 assert(initial_state != final_state);
955 if (num_strings == 4)
956 k_models[k_model].params = strings[3];
957 INIT_MODEL(state_name, &k_models[k_model], k_models[k_model].params, model_params);
958 push(pTransition_list, create_transition(list_element(*pState_list, initial_state),
959 list_element(*pState_list, final_state),
962 k_models[k_model].destructor), 1);
963 free_parsed_list(num_strings, strings);
966 print_state_graph(stdout, *pState_list, *pTransition_list);
970 flags = FLAG_OUTPUT_FULL_CURVE;
973 fprintf(stderr, "unrecognized option '%c'\n", optopt);
974 /* fall through to default case */
976 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);
980 /* check the arguments */
981 assert(*pState_list != NULL); /* no domains! */
982 assert(*pP > 0.0); assert(*pP < 1.0);
983 assert(*pDt_max > 0.0);
985 assert(env->T > 0.0);
987 crosslink_groups(*pState_list, *pTransition_list);
993 <<index tension model>>=
994 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
997 for (i=0; i<n_models; i++)
998 if (STRMATCH(models[i].name, name))
1000 if (i == n_models) {
1001 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
1009 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
1012 for (i=0; i<n_models; i++)
1013 if (STRMATCH(models[i].name, name))
1015 if (i == n_models) {
1016 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
1023 <<initialize model definition>>=
1024 /* requires int num_param_args and char **param_args in the current scope
1026 * INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
1027 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
1029 #define INIT_MODEL(role, model, param_string, param_pointer) \
1031 parse_list_string((param_string), SEP, '{', '}', \
1032 &num_param_args, ¶m_args); \
1033 if (num_param_args != (model)->num_params) { \
1035 "%s model %s expected %d params,\n", \
1036 (role), (model)->name, (model)->num_params); \
1038 "not the %d params in '%s'\n", \
1039 num_param_args, (param_string)); \
1040 assert(num_param_args == (model)->num_params); \
1042 if ((model)->creator) \
1043 param_pointer = (*(model)->creator)(param_args); \
1045 param_pointer = NULL; \
1046 free_parsed_list(num_param_args, param_args); \
1052 \addcontentsline{toc}{section}{Appendicies}
1054 \section{sawsim.c details}
1058 The general layout of our simulation code is:
1069 We include [[math.h]], so don't forget to link to the libm with `-lm'.
1071 #include <assert.h> /* assert() */
1072 #include <stdlib.h> /* malloc(), free(), rand() */
1073 #include <stdio.h> /* fprintf(), stdout */
1074 #include <string.h> /* strlen, strtok() */
1075 #include <math.h> /* exp(), M_PI, sqrt() */
1076 #include <sys/types.h> /* pid_t (returned by getpid()) */
1077 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
1080 #include "tension_balance.h"
1081 #include "k_model.h"
1082 #include "tension_model.h"
1084 #include "accel_k.h"
1088 <<version definition>>
1089 <<flag definitions>>
1090 <<max/min definitions>>
1091 <<string comparison definition and globals>>
1092 <<initialize model definition>>
1093 <<get options definitions>>
1094 <<state definitions>>
1095 <<transition definitions>>
1104 <<create/destroy state>>
1105 <<destroy state list>>
1107 <<create/destroy transition>>
1108 <<destroy transition list>>
1109 <<print state graph>>
1111 <<simulation functions>>
1112 <<boilerplate functions>>
1115 <<boilerplate functions>>=
1117 <<get option functions>>
1123 srand(getpid()*time(NULL)); /* seed rand() */
1127 <<flag definitions>>=
1128 /* in octal b/c of prefixed '0' */
1129 #define FLAG_OUTPUT_FULL_CURVE 01
1130 #define FLAG_OUTPUT_TRANSITION_FORCES 02
1134 static unsigned long int flags = 0;
1137 \subsection{Utilities}
1140 <<max/min definitions>>=
1141 #define MAX(a,b) ((a)>(b) ? (a) : (b))
1142 #define MIN(a,b) ((a)<(b) ? (a) : (b))
1145 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
1146 <<string comparison definition and globals>>=
1147 // Check if two strings match, return 1 if they do
1148 static char *temp_string_A;
1149 static char *temp_string_B;
1150 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
1151 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
1152 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
1153 /* +1 to also compare the '\0' */
1156 We also define a macro for our [[check]] unit testing
1157 <<check relative error>>=
1158 #define CHECK_ERR(max_err, expected, received) \
1160 fail_unless( (received-expected)/expected < max_err, \
1161 "relative error %g >= %g in %s (Expected %g, received %g)", \
1162 (received-expected)/expected, max_err, #received, \
1163 expected, received); \
1164 fail_unless(-(received-expected)/expected < max_err, \
1165 "relative error %g >= %g in %s (Expected %g, received %g)", \
1166 -(received-expected)/expected, max_err, #received, \
1167 expected, received); \
1172 int happens(double probability)
1174 assert(probability >= 0.0); assert(probability <= 1.0);
1175 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*/
1180 \subsection{Adaptive timesteps}
1181 \label{sec.adaptive_dt_implementation}
1183 $F(x)$ increases with $x$, possibly exploding (e.g.\ in the worm-like
1184 chain model), so basing the timestep on the the unfolding probability
1185 at the current tension is dangerous, and we need to search for a $dt$
1186 for which $P(F(x+v*dt)) < P_\text{target}$. There are two cases to
1187 consider. In the most common, no domains have unfolded since the last
1188 step, and we expect the next step to be slightly shorter than the
1189 previous one. In the less common, domains did unfold in the last
1190 step, and we expect the next step to be considerably longer than the
1193 double search_dt(transition_t *transition,
1195 environment_t *env, double x,
1196 double target_prob, double max_dt, double v)
1197 { /* Returns a valid timestep dt in seconds for the given transition.
1198 * Takes a list of all states, a pointer env to the environmental data,
1199 * a starting separation x in m, a target_prob between 0 and 1,
1200 * max_dt in s, stretching velocity v in m/s.
1202 double F, k, dtCur, dtU, dtUCur, dtL, dt;
1204 /* get upper bound using the current position */
1205 F = find_tension(state_list, env, x, 0); /* BUG. repeated calculation */
1206 //printf("Start with x = %g (F = %g)\n", x, F);
1207 k = accel_k(transition->k, F, env, transition->k_params);
1208 //printf("x %g\tF %g\tk %g\n", x, F, k);
1209 dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */
1211 //printf("overshot max_dt\n");
1214 if (v == 0) /* discrete stepping, so force is temporarily constant */
1217 /* set a lower bound on dt too */
1220 /* The dt determined above may produce illegitimate forces or ks.
1221 * Reduce the upper bound until we have valid ks. */
1223 F = find_tension(state_list, env, x+v*dt, 0);
1224 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1227 F = find_tension(state_list, env, x+v*dt, 0);
1229 //printf("Try for dt = %g (F = %g)\n", dt, F);
1230 k = accel_k(transition->k, F, env, transition->k_params);
1231 /* returned k may be -1.0 */
1232 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1233 while (k == -1.0) { /* reduce step until we hit a valid k */
1235 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1236 F = find_tension(state_list, env, x+v*dt, 0);
1237 //printf("Try for dt = %g (F = %g)\n", dt, F);
1238 k = accel_k(transition->k, F, env, transition->k_params);
1239 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1241 assert(dtU > 1e-14); /* timestep to valid k too small */
1242 dtUCur = target_prob / k; /* safe timestep back from x+dtU */
1244 return dt; /* dtU is safe. */
1246 /* dtU wasn't safe, lets see what would be. */
1247 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1248 dt = (dtU + dtL) / 2.0;
1249 F = find_tension(state_list, env, x+v*dt, 0);
1250 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1251 k = accel_k(transition->k, F, env, transition->k_params);
1252 dtCur = target_prob / k;
1253 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1254 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1256 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1257 dtU = dt; /* ... stepping out only dtCur would be safe */
1260 break; /* dtCur = dt */
1262 return MAX(dtUCur, dtL);
1267 To determine $dt$ for an array of potentially different folded domains,
1268 we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
1271 double determine_dt(list_t *state_list, list_t *transition_list,
1272 environment_t *env, double x,
1273 double target_prob, double dt_max, double v)
1274 { /* Returns the timestep dt in seconds.
1275 * Takes the state and transition lists, a pointer to the environment,
1276 * an extension x in m, target_prob between 0 and 1, dt_max in s,
1278 double dt=dt_max, new_dt;
1279 assert(target_prob > 0.0); assert(target_prob < 1.0);
1280 assert(dt_max > 0.0);
1281 assert(state_list != NULL);
1282 assert(transition_list != NULL);
1283 transition_list = head(transition_list);
1284 while (transition_list != NULL) {
1285 if (T(transition_list)->initial_state->num_domains > 0) {
1286 new_dt = search_dt(T(transition_list), state_list, env, x, target_prob, dt, v);
1287 dt = MIN(dt, new_dt);
1289 transition_list = transition_list->next;
1296 \subsection{State data}
1297 \label{sec.state_data}
1299 The domains exist in one of an array of possible states. Each state
1300 has a tension model which determines it's force/extension curve
1301 (possibly [[null]] for rigid states, see Appendix
1302 \ref{sec.null_tension}).
1304 Groups are, for example, for two domain states with WLC tensions; one
1305 with $p=0.4\U{nm}$, $L=10\U{nm}$ and the other with $p=0.4\U{nm}$,
1306 $L=28\U{nm}$. Since the persistence lengths are the same, you would
1307 like to add the contour lengths of all the domains in \emph{both}
1308 states, and plug that total contour length into a single WLC
1309 evaluation (see Section \ref{sec.find_tension}).
1310 <<state definitions>>=
1311 typedef struct state_struct {
1312 char *name; /* name string */
1313 one_dim_fn_t *tension_handler; /* tension handler */
1314 int tension_group; /* for grouping similar tension states */
1315 void *tension_params; /* pointer to tension parameters */
1316 destroy_data_func_t *destroy_tension;
1317 int num_domains; /* number of domains currently in state */
1318 /* cached lookups generated from the above data */
1319 int num_out_trans; /* length of out_trans array */
1320 struct transition_struct **out_trans; /* transitions leaving this state */
1321 int num_grouped_states; /* length of grouped states array */
1322 struct state_struct **grouped_states; /* states grouped with this one (doesn't include self) */
1325 /* get the state data for the current list node */
1326 #define S(list) ((state_t *)(list)->d)
1328 @ [[tension_params]] is a pointer to the parameters used by the
1329 handler function when determining the tension. [[destroy_tension]]
1330 points to a function for freeing [[tension_params]]. We it in the
1331 state structure so that [[destroy_state]] and will have an easier job
1334 [[create_]] and [[destroy_state]] are simple wrappers around
1335 [[malloc]] and [[free]].
1336 <<create/destroy state>>=
1337 state_t *create_state(char *name,
1338 one_dim_fn_t *tension_handler,
1340 void *tension_params,
1341 destroy_data_func_t *destroy_tension,
1344 state_t *ret = (state_t *)malloc(sizeof(state_t));
1345 assert(ret != NULL);
1346 /* make a copy of the name, so we aren't dependent on the original */
1347 ret->name = (char *)malloc(sizeof(char)*(strlen(name)+1));
1348 assert(ret->name != NULL);
1349 strcpy(ret->name, name); /* we know ret->name is long enough */
1351 ret->tension_handler = tension_handler;
1352 ret->tension_group = tension_group;
1353 ret->tension_params = tension_params;
1354 ret->destroy_tension = destroy_tension;
1355 ret->num_domains = num_domains;
1357 ret->num_out_trans = 0;
1358 ret->out_trans = NULL;
1359 ret->num_grouped_states = 0;
1360 ret->grouped_states = NULL;
1363 fprintf(stderr, "generated state %s (%p) with tension handler %p, params %p, group %d.\n", ret->name, ret, ret->tension_handler, ret->tension_params, ret->tension_group);
1368 void destroy_state(state_t *state)
1372 fprintf(stderr, "freeing state %s (%p)\n", state->name, state);
1376 if (state->destroy_tension)
1377 (*state->destroy_tension)(state->tension_params);
1378 if (state->out_trans)
1379 free(state->out_trans);
1380 if (state->grouped_states)
1381 free(state->grouped_states);
1388 We'll be storing the states in a list (see Appendix \ref{sec.lists}),
1389 so we also define a convenience function for cleaning up.
1390 <<destroy state list>>=
1391 void destroy_state_list(list_t *state_list)
1393 state_list = head(state_list);
1394 while (state_list != NULL)
1395 destroy_state((state_t *) pop(&state_list));
1400 We can also define a convenient way to get a state index from it's name.
1402 int state_index(list_t *state_list, char *name)
1405 state_list = head(state_list);
1406 while (state_list != NULL) {
1407 if (STRMATCH(S(state_list)->name, name)) return i;
1408 state_list = state_list->next;
1416 \subsection{Transition data}
1417 \label{sec.transition_data}
1419 Once you've created states (Sections \ref{sec.main},
1420 \ref{sec.model_selection}, and \ref{sec.state_data}), you need to
1421 create transitions between them.
1422 <<transition definitions>>=
1423 typedef struct transition_struct {
1424 state_t *initial_state;
1425 state_t *final_state;
1426 k_func_t *k; /* transition rate model */
1427 void *k_params; /* pointer to k parameters */
1428 destroy_data_func_t *destroy_k;
1431 /* get the transition data for the current list node */
1432 #define T(list) ((transition_t *)(list)->d)
1433 @ [[k]] is a pointer to the function determining the transition rate
1434 for a given tension. [[k_params]] and [[destroy_k]] have similar
1435 roles with regards to [[k]] as [[tension_params]] and
1436 [[destroy_tension]] have with regards to [[tension_handler]] (see
1437 Section \ref{sec.state_data}).
1439 [[create_]] and [[destroy_transition]] are simple wrappers around
1440 [[malloc]] and [[free]].
1441 <<create/destroy transition>>=
1442 transition_t *create_transition(state_t *initial_state, state_t *final_state,
1445 destroy_data_func_t *destroy_k)
1447 transition_t *ret = (transition_t *)malloc(sizeof(transition_t));
1448 assert(ret != NULL);
1449 assert(initial_state != NULL);
1450 assert(final_state != NULL);
1451 ret->initial_state = initial_state;
1452 ret->final_state = final_state;
1454 ret->k_params = k_params;
1455 ret->destroy_k = destroy_k;
1457 fprintf(stderr, "generated transition %p from %s to %s, params %p\n", ret, ret->initial_state->name, ret->final_state->name, ret->k_params);
1462 void destroy_transition(transition_t *transition)
1466 fprintf(stderr, "freeing transition %p (%s to %s)\n", transition, transition->initial_state->name, transition->final_state->name);
1468 if (transition->destroy_k)
1469 (*transition->destroy_k)(transition->k_params);
1476 We'll be storing the transitions in a list (see Appendix
1477 \ref{sec.lists}), so we also define a convenience function for
1479 <<destroy transition list>>=
1480 void destroy_transition_list(list_t *transition_list)
1482 transition_list = head(transition_list);
1483 while (transition_list != NULL)
1484 destroy_transition((transition_t *) pop(&transition_list));
1489 \subsection{Graphviz output}
1490 \label{sec.graphviz_output}
1492 <<print state graph>>=
1493 void print_state_graph(FILE *file, list_t *state_list, list_t *transition_list)
1495 state_list = head(state_list);
1496 transition_list = head(transition_list);
1497 fprintf(file, "digraph states {\n node [shape=record,color=black,fillcolor=white,style=filled];\n\n");
1499 fprintf(file, " n%d [label=\"%s\"]\n", state_index(state_list, S(state_list)->name), S(state_list)->name);
1500 if (state_list->next == NULL) break;
1501 state_list = state_list->next;
1503 fprintf(file, "\n");
1504 while (transition_list != NULL) {
1505 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));
1506 transition_list = transition_list->next;
1508 fprintf(file, "}\n");
1512 \subsection{Domain model and group handling}
1514 <<group functions>>=
1515 <<crosslink groups>>
1516 <<get active groups>>
1519 Scan through a list of states and construct an array of all the
1521 <<get active groups>>=
1522 void get_active_groups(list_t *state_list,
1523 int *pNum_active_groups,
1524 one_dim_fn_t ***pPer_group_handlers,
1525 void ***pPer_group_data /* pointer to array of pointers to tension_handler_data_t */)
1527 void **active_groups=NULL;
1528 one_dim_fn_t *handler, *old_handler, **per_group_handlers=NULL;
1529 void **per_group_data=NULL; /* array of pointers to tension_handler_data_t */
1530 int i, j, num_domains, group, old_group, num_active_groups=0;
1531 list_t *active_states_list=NULL;
1532 tension_handler_data_t *tdata=NULL;
1535 /* get all the active groups in a list */
1536 state_list = head(state_list);
1538 fprintf(stderr, "%s:\t", __FUNCTION__);
1539 list_print(stderr, state_list, "states");
1541 while (state_list != NULL) {
1542 state = S(state_list);
1543 handler = state->tension_handler;
1544 group = state->tension_group;
1545 num_domains = state->num_domains;
1546 if (list_index(active_states_list, state) == -1) {
1547 /* we haven't added this state (or it's associated group) yet */
1548 if (num_domains > 0 && handler != NULL) { /* add the state */
1549 num_active_groups++;
1550 push(&active_states_list, state, 1);
1552 fprintf(stderr, "%s:\tactive %s state (%p) with %d domain(s)\n", __FUNCTION__, state->name, state, state->num_domains);
1554 for (i=0; i < state->num_grouped_states; i++) {
1555 if (state->grouped_states[i]->num_domains > 0) {
1556 /* add this related (and active) state now */
1557 assert(state->grouped_states[i]->tension_handler == handler);
1558 assert(state->grouped_states[i]->tension_group == group);
1559 push(&active_states_list, state->grouped_states[i], 1);
1561 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);
1565 } /* else empty state or NULL handler */
1566 } /* else state already added */
1567 state_list = state_list->next;
1570 fprintf(stderr, "%s:\t", __FUNCTION__);
1571 list_print(stderr, active_states_list, "active states");
1574 assert(num_active_groups <= list_length(active_states_list));
1575 assert(num_active_groups > 0);
1577 /* allocate the arrays we'll be returning */
1578 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1579 assert(per_group_handlers != NULL);
1580 per_group_data = (void **)malloc(sizeof(void *)*num_active_groups);
1581 assert(per_group_data != NULL);
1583 fprintf(stderr, "%s:\tallocating per group data %p:", __FUNCTION__, per_group_data);
1585 for (i=0; i<num_active_groups; i++) { /* we might want to recycle any group data structs allocated in the previous round */
1586 per_group_data[i] = malloc(sizeof(tension_handler_data_t));
1587 assert(per_group_data[i] != NULL);
1589 fprintf(stderr, " %d,%p", i, per_group_data[i]);
1593 fprintf(stderr, "\n");
1596 /* populate the arrays we'll be returning */
1597 active_states_list = head(active_states_list);
1599 j = -1; /* j is the current active _group_ index */
1600 while (active_states_list != NULL) {
1601 state = (state_t *) pop(&active_states_list);
1602 handler = state->tension_handler;
1603 group = state->tension_group;
1604 num_domains = state->num_domains;
1605 if (handler != old_handler || group != old_group) {
1606 j++; /* move to the next group */
1607 tdata = (tension_handler_data_t *) per_group_data[j];
1608 per_group_handlers[j] = handler;
1609 tdata->group_tension_params = NULL;
1611 tdata->persist = NULL; /* we might want to grab any appropriate data from the old persist here... */
1614 fprintf(stderr, "%s:\tadding active state %s to group %d (%p:%d) with %d domain(s)\n", __FUNCTION__, state->name, j, handler, group, num_domains);
1616 for (i=0; i<num_domains; i++) { /* add num_domains copies of the tension params */
1617 push(&tdata->group_tension_params, state->tension_params, 1);
1620 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);
1622 old_handler = handler;
1626 /* free old groups */
1627 if (*pPer_group_handlers != NULL)
1628 free(*pPer_group_handlers);
1629 if (*pPer_group_data != NULL) {
1630 for (i=0; i<*pNum_active_groups; i++)
1631 free((*pPer_group_data)[i]);
1632 free(*pPer_group_data);
1635 /* set new groups */
1637 fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_data, per_group_data[0]);
1639 *pNum_active_groups = num_active_groups;
1640 *pPer_group_handlers = per_group_handlers;
1641 *pPer_group_data = per_group_data;
1645 <<crosslink groups>>=
1646 void crosslink_groups(list_t *state_groups, list_t *transition_list)
1648 list_t *list, *out_trans=NULL, *associates=NULL;
1649 one_dim_fn_t *handler;
1652 state_groups = head(state_groups);
1653 while (state_groups != NULL) {
1654 /* find transitions leaving this state */
1656 fprintf(stderr, "finding transitions leaving %s\n", S(state_groups)->name);
1658 list = head(transition_list);
1659 while (list != NULL) {
1660 if (T(list)->initial_state == S(state_groups) && list_index(associates, T(list)->initial_state) != -1) {
1661 push(&out_trans, T(list), 1);
1665 n = list_length(out_trans);
1666 S(state_groups)->num_out_trans = n;
1667 S(state_groups)->out_trans = (transition_t **)malloc(sizeof(transition_t *)*n);
1668 assert(S(state_groups)->out_trans != NULL);
1669 for (i=0; i<n; i++) {
1670 S(state_groups)->out_trans[i] = (transition_t *)pop(&out_trans);
1673 /* find states grouped with this state */
1675 fprintf(stderr, "finding states grouped with %s\n", S(state_groups)->name);
1677 handler = S(state_groups)->tension_handler;
1678 group = S(state_groups)->tension_group;
1679 list = head(state_groups);
1680 while (list != NULL) {
1681 if (S(list)->tension_handler == handler && S(list)->tension_group == group && list != state_groups) {
1682 push(&associates, S(list), 1);
1686 n = list_length(associates);
1687 S(state_groups)->num_grouped_states = n;
1688 S(state_groups)->grouped_states = (state_t **)malloc(sizeof(state_t *)*n);
1689 assert(S(state_groups)->grouped_states != NULL);
1690 for (i=0; i<n; i++) {
1691 S(state_groups)->grouped_states[i] = (state_t *)pop(&associates);
1693 state_groups = state_groups->next;
1699 \section{String parsing}
1701 For handling command line arguments, we need parse delimited strings (e.g.\ [["param1,param2,param3..."]]).
1702 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1703 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1704 need to take care of parsing those parameters themselves.
1705 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1711 <<parse definitions>>
1712 <<parse declarations>>
1716 <<parse module makefile lines>>=
1717 NW_SAWSIM_MODS += parse
1718 CHECK_BINS += check_parse
1722 <<parse definitions>>=
1723 #define SEP ',' /* argument separator character */
1726 <<parse declarations>>=
1727 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1728 int *num, char ***string_array);
1729 extern void free_parsed_list(int num, char **string_array);
1732 [[parse_list_string]] allocates memory, don't forget to free it
1733 afterward with [[free_parsed_list]]. It does not alter the original.
1735 The string may start off with a [[deeper]] character
1736 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1737 [[x,y]], where the pointer is one character in on the copied string.
1738 However, when we go to free the memory, we need a pointer to the
1739 beginning of the string. In order to accommodate this for a string
1740 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1741 the first $N$ elements point to the separated arguments, and let the
1742 last element point to the start of the copied string regardless of
1744 <<parse delimited list string functions>>=
1745 /* TODO, split out into parse.hc */
1746 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1749 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1750 if (string[i] == deeper) {depth++;}
1751 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1757 void parse_list_string(char *string, char sep, char deeper, char shallower,
1758 int *num, char ***string_array)
1760 char *str=NULL, **ret=NULL;
1762 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
1764 *string_array = NULL;
1767 /* make a copy of the string, so we don't change the original */
1768 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1769 assert(str != NULL);
1770 strcpy(str, string); /* we know str is long enough */
1771 /* count the number of regions, so we can allocate pointers to them */
1774 n++; i++; /* move on to next argument */
1775 i += next_delim_index(str+i, sep, deeper, shallower);
1776 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1777 } while (str[i] != '\0');
1778 ret = (char **)malloc(sizeof(char *)*(n+1));
1779 assert(ret != NULL);
1780 /* replace the separators with '\0' & assign pointers */
1781 ret[n] = str; /* point to the front of the copied string */
1784 for(i=1; i<n; i++) {
1785 j += next_delim_index(str+j, sep, deeper, shallower);
1787 ret[i] = str+j; /* point to the separated arguments */
1789 /* strip off bounding braces, if any */
1790 for(i=0; i<n; i++) {
1791 if (ret[i][0] == deeper) {
1795 if (ret[i][strlen(ret[i])-1] == shallower)
1796 ret[i][strlen(ret[i])-1] = '\0';
1799 *string_array = ret;
1802 void free_parsed_list(int num, char **string_array)
1805 /* frees all items, since they were allocated as a single string*/
1806 free(string_array[num]);
1807 /* and free the array of pointers */
1813 \subsection{Parsing implementation}
1819 <<parse delimited list string functions>>
1823 #include <assert.h> /* assert() */
1824 #include <stdlib.h> /* NULL */
1825 #include <stdio.h> /* fprintf(), stdout *//*!!*/
1826 #include <string.h> /* strlen() */
1830 \subsection{Parsing unit tests}
1832 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1834 <<parse check includes>>
1835 <<string comparison definition and globals>>
1836 <<check relative error>>
1837 <<parse test suite>>
1838 <<main check program>>
1841 <<parse check includes>>=
1842 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
1843 #include <stdio.h> /* printf() */
1844 #include <assert.h> /* assert() */
1845 #include <string.h> /* strlen() */
1850 <<parse test suite>>=
1851 <<parse list string tests>>
1852 <<parse suite definition>>
1855 <<parse suite definition>>=
1856 Suite *test_suite (void)
1858 Suite *s = suite_create ("k model");
1859 <<parse list string test case defs>>
1861 <<parse list string test case add>>
1866 <<parse list string tests>>=
1869 START_TEST(test_next_delim_index)
1871 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
1872 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
1873 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
1874 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
1875 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
1879 START_TEST(test_parse_list_null)
1883 parse_list_string(NULL, SEP, '{', '}',
1884 &num_param_args, ¶m_args);
1885 fail_unless(num_param_args == 0, NULL);
1886 fail_unless(param_args == NULL, NULL);
1889 START_TEST(test_parse_list_single_simple)
1893 parse_list_string("arg", SEP, '{', '}',
1894 &num_param_args, ¶m_args);
1895 fail_unless(num_param_args == 1, NULL);
1896 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
1899 START_TEST(test_parse_list_single_compound)
1903 parse_list_string("{x,y,z}", SEP, '{', '}',
1904 &num_param_args, ¶m_args);
1905 fail_unless(num_param_args == 1, NULL);
1906 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
1909 START_TEST(test_parse_list_double_simple)
1913 parse_list_string("abc,def", SEP, '{', '}',
1914 &num_param_args, ¶m_args);
1915 fail_unless(num_param_args == 2, NULL);
1916 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
1917 fail_unless(STRMATCH(param_args[1],"def"), NULL);
1920 START_TEST(test_parse_list_compound)
1924 parse_list_string("abc,{def,ghi}", SEP, '{', '}',
1925 &num_param_args, ¶m_args);
1926 fail_unless(num_param_args == 2, NULL);
1927 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
1928 fail_unless(STRMATCH(param_args[1],"def,ghi"), NULL);
1932 <<parse list string test case defs>>=
1933 TCase *tc_parse_list_string = tcase_create("parse list string");
1935 <<parse list string test case add>>=
1936 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
1937 tcase_add_test(tc_parse_list_string, test_parse_list_null);
1938 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
1939 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
1940 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
1941 tcase_add_test(tc_parse_list_string, test_parse_list_compound);
1942 suite_add_tcase(s, tc_parse_list_string);
1946 \section{Unit tests}
1948 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1955 <<check relative error>>
1958 <<main check program>>
1970 <<determine dt tests>>
1972 <<does domain transition tests>>
1973 <<randomly unfold domains tests>>
1974 <<suite definition>>
1977 <<suite definition>>=
1978 Suite *test_suite (void)
1980 Suite *s = suite_create ("sawsim");
1981 <<F test case defs>>
1982 <<determine dt test case defs>>
1983 <<happens test case defs>>
1984 <<does domain transition test case defs>>
1985 <<randomly unfold domains test case defs>>
1988 <<determine dt test case add>>
1989 <<happens test case add>>
1990 <<does domain transition test case add>>
1991 <<randomly unfold domains test case add>>
1994 tcase_add_checked_fixture(tc_strip_address,
1995 setup_strip_address,
1996 teardown_strip_address);
2002 <<main check program>>=
2006 Suite *s = test_suite();
2007 SRunner *sr = srunner_create(s);
2008 srunner_run_all(sr, CK_ENV);
2009 number_failed = srunner_ntests_failed(sr);
2011 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
2015 \subsection{F tests}
2017 <<worm-like chain tests>>
2019 <<F test case defs>>=
2020 <<worm-like chain test case def>>
2022 <<F test case add>>=
2023 <<worm-like chain test case add>>
2026 <<worm-like chain tests>>=
2027 extern double wlc(double x, double T, double p, double L);
2028 START_TEST(test_wlc_at_zero)
2030 double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
2031 fail_unless(abs(wlc(x, T, p, L)) < lim, \
2032 "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
2033 x, T, p, L, abs(wlc(x, T, p, L)), lim);
2037 START_TEST(test_wlc_at_half)
2039 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
2040 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
2041 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
2043 * linear term = x/L = 0.5
2044 * nonlinear + linear = 0.75 + 0.5 = 1.25
2045 * wlc = 10e21*1.25 = 12.5e21
2047 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
2048 "wlc(%g, %g, %g, %g) = %g != %g",
2049 x, T, p, L, wlc(x, T, p, L), 12.5e21);
2053 <<worm-like chain test case def>>=
2054 TCase *tc_wlc = tcase_create("WLC");
2057 <<worm-like chain test case add>>=
2058 tcase_add_test(tc_wlc, test_wlc_at_zero);
2059 tcase_add_test(tc_wlc, test_wlc_at_half);
2060 suite_add_tcase(s, tc_wlc);
2063 \subsection{Model tests}
2065 Check the searching with [[linear_k]].
2066 Check overwhelming force treatment with the heavyside-esque step [[?]].
2068 Hmm, I'm not really sure what I was doing here...
2069 <<determine dt tests>>=
2070 double linear_k(double F, environment_t *env, void *params)
2072 double Fnot = *(double *)params;
2076 START_TEST(test_determine_dt_linear_k)
2079 transition_t unfold={0};
2080 environment_t env = {0};
2081 double F[]={0,1,2,3,4,5,6};
2082 double dt_max=3.0, Fnot=3.0, x=1.0, targ_p=0.1, v=1.0;
2083 list_t *state_list=NULL, *transition_list=NULL;
2086 folded.tension_handler = &hooke_handler;
2087 folded.num_domains = 1;
2088 unfold.initial_state = &folded;
2089 unfold.k = linear_k;
2090 unfold.k_params = &Fnot;
2091 push(&state_list, &folded, 1);
2092 push(&transition_list, &unfold, 1);
2093 for( i=0; i < sizeof(F)/sizeof(double); i++) {
2094 //fail_unless(determine_dt(state_list, transition_list, &env, x, targ_p, dt_max, v)==1, NULL);
2099 <<determine dt test case defs>>=
2100 TCase *tc_determine_dt = tcase_create("Determine dt");
2102 <<determine dt test case add>>=
2103 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
2104 suite_add_tcase(s, tc_determine_dt);
2109 <<happens test case defs>>=
2111 <<happens test case add>>=
2114 <<does domain transition tests>>=
2116 <<does domain transition test case defs>>=
2118 <<does domain transition test case add>>=
2121 <<randomly unfold domains tests>>=
2123 <<randomly unfold domains test case defs>>=
2125 <<randomly unfold domains test case add>>=
2129 \section{Balancing group extensions}
2130 \label{sec.tension_balance}
2132 For a given total extension $x$ (of the piezo), the various domain
2133 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
2134 amounts, and we need to tweak the portion that each extends to balance
2135 the tension amongst the active groups. Since the tension balancing is
2136 separable from the bulk of the simulation, we break it out into a
2137 separate module. The interface is defined in [[tension_balance.h]],
2138 the implementation in [[tension_balance.c]], and the unit testing in
2139 [[check_tension_balance.c]]
2141 <<tension-balance.h>>=
2142 #ifndef TENSION_BALANCE_H
2143 #define TENSION_BALANCE_H
2145 <<tension balance function declaration>>
2146 #endif /* TENSION_BALANCE_H */
2149 <<tension balance functions>>=
2150 <<one dimensional bracket>>
2151 <<one dimensional solve>>
2153 <<group stiffness function>>
2154 <<chain stiffness function>>
2155 <<tension balance function>>
2158 <<tension balance module makefile lines>>=
2159 NW_SAWSIM_MODS += tension_balance
2160 CHECK_BINS += check_tension_balance
2161 check_tension_balance_MODS =
2164 The entire force balancing problem reduces to a solving two nested
2165 one-dimensional root-finding problems. First we define the one
2166 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
2167 populated group). $x(x_0)$ is also strictly monotonically increasing,
2168 so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
2169 stores the last successful value of $x$, since we expect to be taking
2170 small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
2171 indicates that the groups have changed.
2172 <<tension balance function declaration>>=
2173 double tension_balance(int num_tension_groups,
2174 one_dim_fn_t **tension_handlers,
2175 void **params, /* array of pointers to tension_handler_data_t */
2180 <<tension balance function>>=
2181 double tension_balance(int num_tension_groups,
2182 one_dim_fn_t **tension_handlers,
2183 void **params, /* array of pointers to tension_handler_data_t */
2188 static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL};
2189 double F, xo_guess, xo, lb, ub;
2190 double min_relx=1e-6, min_rely=1e-6;
2191 int max_steps=200, i;
2193 assert(num_tension_groups > 0);
2194 assert(tension_handlers != NULL);
2195 assert(params != NULL);
2196 assert(num_tension_groups > 0);
2198 if (last_x == -1 || x_xo_data.xi == NULL) { /* new group setup, reset x_xo_data */
2199 /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
2200 if (x_xo_data.xi != NULL)
2202 /* construct new x_xo_data */
2203 x_xo_data.n_groups = num_tension_groups;
2204 x_xo_data.tension_handler = tension_handlers;
2205 x_xo_data.handler_data = params;
2207 fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, x_xo_data.n_groups, x_xo_data.tension_handler, x_xo_data.handler_data, x_xo_data.handler_data[0]);
2209 x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
2210 assert(x_xo_data.xi != NULL);
2211 for (i=0; i<num_tension_groups; i++)
2212 x_xo_data.xi[i] = last_x;
2215 if (num_tension_groups == 1) { /* only one group, no balancing required */
2217 x_xo_data.xi[0] = xo;
2219 fprintf(stderr, "%s:\tonly one active group, no balancing required\n", __FUNCTION__);
2223 fprintf(stderr, "%s:\tbalancing for x = %g with", __FUNCTION__, x);
2224 for(i=0; i<num_tension_groups; i++) fprintf(stderr, "\tx[%d]=%g", i, x_xo_data.xi[i]);
2225 fprintf(stderr, "\n");
2227 if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
2228 if (x == 0) xo_guess = length_scale;
2229 else xo_guess = x/num_tension_groups;
2231 fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
2233 oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
2234 } else { /* work off the last balanced point */
2236 fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
2239 lb = x_xo_data.xi[0];
2240 ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */
2241 } else if (x < last_x) {
2242 lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
2243 ub = x_xo_data.xi[0];
2244 } else { /* x == last_x */
2245 fprintf(stderr, "not moving\n");
2246 lb= x_xo_data.xi[0];
2247 ub= x_xo_data.xi[0];
2251 fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
2252 __FUNCTION__, x, lb, ub);
2254 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
2255 F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
2257 if (num_tension_groups > 1) {
2258 fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
2259 for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
2260 fprintf(stderr, "\n");
2265 F = (*tension_handlers[0])(xo, params[0]);
2267 if (stiffness != NULL)
2268 *stiffness = chain_stiffness(&x_xo_data, min_relx);
2275 <<tension balance internal definitions>>=
2276 <<x of x0 definitions>>
2279 <<x of x0 definitions>>=
2280 typedef struct x_of_xo_data_struct {
2282 one_dim_fn_t **tension_handler; /* array of fn pointers */
2283 void **handler_data; /* array of pointers to tension_handler_data_t */
2284 double *xi; /* array of group extensions */
2289 double x_of_xo(double xo, void *pdata)
2291 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
2292 double F, x=0, xi, xi_guess, lb, ub, slope;
2294 double min_relx=1e-6, min_rely=1e-6;
2296 assert(data->n_groups > 0);
2299 fprintf(stderr, "%s:\tget initial tension at %g with %d groups, handler %p, data %p\n", __FUNCTION__, xo, data->n_groups, data->tension_handler[0], data->handler_data[0]);
2302 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
2304 fprintf(stderr, "%s:\tfinding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
2308 if (data->xi) data->xi[0] = xo;
2309 for (i=1; i < data->n_groups; i++) {
2310 if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
2311 else xi_guess = data->xi[i];
2313 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]);
2315 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
2316 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
2318 fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
2322 if (data->xi) data->xi[i] = xi;
2325 fprintf(stderr, "%s:\tfound x(x0=%g) = %g\n", __FUNCTION__, xo, x);
2331 <<group stiffness function>>=
2332 double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double x, double relx)
2334 double dx, xl, F, Fl, stiffness;
2336 assert(relx > 0 && relx < 1);
2337 if (x == 0 || relx == 0) {
2338 dx = 1e-12; /* HACK, how to get length scale? */
2348 F = (*tension_handler)(x, handler_data);
2349 Fl = (*tension_handler)(xl, handler_data);
2350 stiffness = (F-Fl)/dx;
2351 assert(stiffness > 0);
2356 <<chain stiffness function>>=
2357 double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx)
2361 /* add group stiffnesses in series */
2362 for (i=0; i < x_xo_data->n_groups; i++) {
2364 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);
2367 stiffness += 1.0/group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx);
2369 stiffness = 1.0 / stiffness;
2375 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
2376 number of steps for monotonically increasing $f(x)$. Simple
2377 bisection, so it's robust and converges fairly quickly. You can set
2378 the maximum number of steps to take, but if you have enough processor
2379 power, it's better to limit your precision with the relative limits
2380 [[min_relx]] and [[y]]. These ensure that the width of the bracket is
2381 small on the length scale of it's larger side. Note that \emph{both}
2382 relative limits must be satisfied for the search to stop.
2383 <<one dimensional function definition>>=
2384 /* equivalent to gsl_func_t */
2385 typedef double one_dim_fn_t(double x, void *params);
2387 <<one dimensional solve declaration>>=
2388 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2389 double min_relx, double min_rely, int max_steps,
2393 <<one dimensional solve>>=
2394 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2395 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2396 double min_relx, double min_rely, int max_steps,
2399 double yx, ylb, yub, x;
2402 ylb = (*f)(lb, params);
2403 yub = (*f)(ub, params);
2405 /* check some simple cases */
2407 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bracketed */
2408 else return lb; /* any x on the range [lb, ub] would work */
2410 if (ylb == y) { x=lb; yx=ylb; goto end; }
2411 if (yub == y) { x=ub; yx=yub; goto end; }
2414 fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2416 assert(ylb < y); assert(yub > y);
2417 for (i=0; i<max_steps; i++) {
2419 yx = (*f)(x, params);
2421 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);
2423 if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2424 else if (yx > y) { ub=x; yub=yx; }
2425 else /* yx < y */ { lb=x; ylb=yx; }
2430 fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2436 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2437 Generate bracketing $x$ values through bisection or exponential growth.
2438 <<one dimensional bracket declaration>>=
2439 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2442 <<one dimensional bracket>>=
2443 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2445 double yx, step, x=xguess;
2447 yx = (*f)(x, params);
2449 fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2456 fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2460 if (x == 0) x = length_scale/2.0; /* HACK */
2463 yx = (*f)(x, params);
2465 fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2470 fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2473 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2477 \subsection{Balancing implementation}
2479 <<tension-balance.c>>=
2482 <<tension balance includes>>
2483 #include "tension_balance.h"
2485 double length_scale = 1e-10; /* HACK */
2487 <<tension balance internal definitions>>
2488 <<tension balance functions>>
2491 <<tension balance includes>>=
2492 #include <assert.h> /* assert() */
2493 #include <stdlib.h> /* NULL */
2494 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2495 #include <stdio.h> /* fprintf(), stdout */
2499 \subsection{Balancing unit tests}
2501 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2502 <<check-tension-balance.c>>=
2503 <<tension balance check includes>>
2504 <<tension balance test suite>>
2505 <<main check program>>
2508 <<tension balance check includes>>=
2509 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2510 #include <assert.h> /* assert() */
2513 #include "tension_balance.h"
2516 <<tension balance test suite>>=
2517 <<tension balance function tests>>
2518 <<tension balance suite definition>>
2521 <<tension balance suite definition>>=
2522 Suite *test_suite (void)
2524 Suite *s = suite_create ("tension balance");
2525 <<tension balance function test case defs>>
2527 <<tension balance function test case adds>>
2532 <<tension balance function tests>>=
2533 <<check relative error>>
2535 double hooke(double x, void *pK)
2538 return *((double*)pK) * x;
2541 START_TEST(test_single_function)
2543 double k=5, x=3, last_x=2.0, F;
2544 one_dim_fn_t *handlers[] = {&hooke};
2545 void *data[] = {&k};
2546 F = tension_balance(1, handlers, data, last_x, x, NULL);
2547 fail_unless(F = k*x, NULL);
2552 We can also test balancing two springs with different spring constants.
2553 The analytic solution for a general number of springs is given in Appendix \ref{sec.math_hooke}.
2554 <<tension balance function tests>>=
2555 START_TEST(test_double_hooke)
2557 double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
2558 one_dim_fn_t *handlers[] = {&hooke, &hooke};
2559 void *data[] = {&k1, &k2};
2560 F = tension_balance(2, handlers, data, last_x, x, NULL);
2564 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2565 //CHECK_ERR(1e-6, x1e, xi[0]);
2566 //CHECK_ERR(1e-6, x2e, xi[1]);
2567 CHECK_ERR(1e-6, Fe, F);
2571 TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
2573 <<tension balance function test case defs>>=
2574 TCase *tc_tbfunc = tcase_create("tension balance function");
2577 <<tension balance function test case adds>>=
2578 tcase_add_test(tc_tbfunc, test_single_function);
2579 tcase_add_test(tc_tbfunc, test_double_hooke);
2580 suite_add_tcase(s, tc_tbfunc);
2586 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2587 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2588 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2594 <<list definitions>>
2595 <<list declarations>>
2599 <<list declarations>>=
2600 <<head and tail declarations>>
2601 <<list length declaration>>
2602 <<push declaration>>
2604 <<list destroy declaration>>
2605 <<list to array declaration>>
2606 <<move declaration>>
2607 <<sort declaration>>
2608 <<list index declaration>>
2609 <<list element declaration>>
2610 <<unique declaration>>
2611 <<list print declaration>>
2615 <<create and destroy node>>
2630 <<list module makefile lines>>=
2631 NW_SAWSIM_MODS += list
2632 CHECK_BINS += check_list
2636 <<list definitions>>=
2637 typedef struct list_struct {
2638 struct list_struct *next;
2639 struct list_struct *prev;
2643 <<comparison function typedef>>
2646 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2647 <<head and tail declarations>>=
2648 list_t *head(list_t *list);
2649 list_t *tail(list_t *list);
2652 list_t *head(list_t *list)
2656 while (list->prev != NULL) {
2662 list_t *tail(list_t *list)
2666 while (list->next != NULL) {
2673 <<list length declaration>>=
2674 int list_length(list_t *list);
2677 int list_length(list_t *list)
2684 while (list->next != NULL) {
2692 [[push]] inserts a node relative to the current position in the list
2693 without changing the current position.
2694 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2695 so we set it to that of the pushed domain.
2696 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2697 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2698 <<push declaration>>=
2699 void push(list_t **pList, void *data, int below);
2702 void push(list_t **pList, void *data, int below)
2704 list_t *list, *new_node;
2705 assert(pList != NULL);
2707 new_node = create_node(data);
2710 else if (below==0) { /* insert above */
2711 if (list->prev != NULL)
2712 list->prev->next = new_node;
2713 new_node->prev = list->prev;
2714 list->prev = new_node;
2715 new_node->next = list;
2716 } else { /* insert below */
2717 if (list->next != NULL)
2718 list->next->prev = new_node;
2719 new_node->next = list->next;
2720 list->next = new_node;
2721 new_node->prev = list;
2726 [[pop]] removes the current domain node, moving the current position
2727 to the node after it, unless that node is [[NULL]], in which case move
2728 the current position to the node before it.
2729 <<pop declaration>>=
2730 void *pop(list_t **pList);
2733 void *pop(list_t **pList)
2735 list_t *list, *popped;
2737 assert(pList != NULL);
2739 assert(list != NULL); /* not an empty list */
2741 /* bypass the popped node */
2742 if (list->prev != NULL)
2743 list->prev->next = list->next;
2744 if (list->next != NULL) {
2745 list->next->prev = list->prev;
2746 *pList = list->next;
2748 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
2750 destroy_node(popped);
2755 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
2756 If the cleanup function is [[NULL]], no cleanup function is called.
2757 The nodes are not popped in any particular order.
2758 <<list destroy declaration>>=
2759 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
2762 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
2766 assert(pList != NULL);
2769 assert(list != NULL); /* not an empty list */
2770 while (list != NULL) {
2772 if (destroy != NULL)
2778 [[list_to_array]] converts a list to an array of pointers, preserving order.
2779 <<list to array declaration>>=
2780 void list_to_array(list_t *list, int *length, void ***array);
2783 void list_to_array(list_t *list, int *pLength, void ***pArray)
2787 assert(list != NULL);
2788 assert(pLength != NULL);
2789 assert(pArray != NULL);
2791 length = list_length(list);
2792 array = (void **)malloc(sizeof(void *)*length);
2793 assert(array != NULL);
2796 while (list != NULL) {
2797 array[i++] = list->d;
2805 [[move]] moves the current element along the list in either direction.
2806 <<move declaration>>=
2807 void move(list_t **pList, int downstream);
2810 void move(list_t **pList, int downstream)
2812 list_t *list, *flipper;
2814 assert(pList != NULL);
2815 list = *pList; /* a>B>c>d */
2816 assert(list != NULL); /* not an empty list */
2817 if (downstream == 0)
2818 flipper = list->prev; /* flipper is A */
2820 flipper = list->next; /* flipper is C */
2821 /* ensure there is somewhere to go in the direction we're heading */
2822 assert(flipper != NULL);
2823 /* remove the current element */
2824 data = pop(&list); /* data=B, a>C>d */
2825 /* and put it back in in the correct spot */
2827 if (downstream == 0) {
2828 push(&list, data, 0); /* b>A>c>d */
2829 list = list->prev; /* B>a>c>d */
2831 push(&list, data, 1); /* a>C>b>d */
2832 list = list->next; /* a>c>B>d */
2838 [[sort]] sorts a list from smallest to largest according to a user
2839 specified comparison.
2840 <<comparison function typedef>>=
2841 typedef int (comparison_fn_t)(void *A, void *B);
2844 <<int comparison function>>=
2845 int int_comparison(void *A, void *B)
2850 if (a > b) return 1;
2851 else if (a == b) return 0;
2856 <<sort declaration>>=
2857 void sort(list_t **pList, comparison_fn_t *comp);
2859 Here we use the bubble sort algorithm.
2861 void sort(list_t **pList, comparison_fn_t *comp)
2865 assert(pList != NULL);
2867 assert(list != NULL); /* not an empty list */
2868 while (stable == 0) {
2871 while (list->next != NULL) {
2872 if ((*comp)(list->d, list->next->d) > 0) {
2874 move(&list, 1 /* downstream */);
2884 [[list_index]] finds the location of [[data]] in [[list]]. Returns
2885 [[-1]] if [[data]] is not in [[list]].
2886 <<list index declaration>>=
2887 int list_index(list_t *list, void *data);
2891 int list_index(list_t *list, void *data)
2895 while (list != NULL) {
2896 if (list->d == data) return i;
2905 [[list_element]] returns the data bound to the $i^\text{th}$ node of the list.
2906 <<list element declaration>>=
2907 void *list_element(list_t *list, int i);
2910 void *list_element(list_t *list, int i)
2914 while (list != NULL) {
2915 if (i == j) return list->d;
2925 [[unique]] removes duplicates from a sorted list according to a user
2926 specified comparison. Don't do this unless you have other pointers
2927 any dynamically allocated data in your list, because [[unique]] just
2928 drops any non-unique data without freeing it.
2929 <<unique declaration>>=
2930 void unique(list_t **pList, comparison_fn_t *comp);
2933 void unique(list_t **pList, comparison_fn_t *comp)
2936 assert(pList != NULL);
2938 assert(list != NULL); /* not an empty list */
2940 while (list->next != NULL) {
2941 if ((*comp)(list->d, list->next->d) == 0) {
2950 [[list_print]] prints a list to a [[FILE *]] stream.
2951 <<list print declaration>>=
2952 void list_print(FILE *file, list_t *list, char *list_name);
2955 void list_print(FILE *file, list_t *list, char *list_name)
2957 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
2959 while (list != NULL) {
2960 fprintf(file, " %p", list->d);
2963 fprintf(file, "\n");
2967 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
2968 <<create and destroy node>>=
2969 list_t *create_node(void *data)
2971 list_t *ret = (list_t *)malloc(sizeof(list_t));
2972 assert(ret != NULL);
2979 void destroy_node(list_t *node)
2985 The user must free the data pointed to by the node on their own.
2987 \subsection{List implementation}
2997 #include <assert.h> /* assert() */
2998 #include <stdlib.h> /* malloc(), free(), rand() */
2999 #include <stdio.h> /* fprintf(), stdout, FILE */
3000 #include "global.h" /* destroy_data_func_t */
3003 \subsection{List unit tests}
3005 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3007 <<list check includes>>
3009 <<main check program>>
3012 <<list check includes>>=
3013 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
3014 #include <stdio.h> /* FILE */
3020 <<list test suite>>=
3023 <<list suite definition>>
3026 <<list suite definition>>=
3027 Suite *test_suite (void)
3029 Suite *s = suite_create ("list");
3030 <<push test case defs>>
3031 <<pop test case defs>>
3033 <<push test case adds>>
3034 <<pop test case adds>>
3040 START_TEST(test_push)
3045 push(&list, (void *)i, 1);
3046 fail_unless(list == head(list), NULL);
3047 fail_unless( (int)list->d == 0 );
3048 for (i=0; i<n; i++) {
3049 /* we expect 0, n-1, n-2, ..., 1 */
3052 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
3058 <<push test case defs>>=
3059 TCase *tc_push = tcase_create("push");
3062 <<push test case adds>>=
3063 tcase_add_test(tc_push, test_push);
3064 suite_add_tcase(s, tc_push);
3069 <<pop test case defs>>=
3071 <<pop test case adds>>=
3074 \section{Function string evaluation}
3076 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).
3077 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
3078 The solution is to run a scripting language as a subprocess accessed via pipes.
3079 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
3081 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
3082 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
3083 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.
3084 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
3086 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.
3087 Then you can either statically or dynamically link to those hardcoded functions.
3088 While much less flexible, this approach would be a fairly simple-to-implement fallback.
3090 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
3091 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
3094 #ifndef STRING_EVAL_H
3095 #define STRING_EVAL_H
3097 <<string eval setup declaration>>
3098 <<string eval function declaration>>
3099 <<string eval teardown declaration>>
3100 #endif /* STRING_EVAL_H */
3103 <<string eval module makefile lines>>=
3104 NW_SAWSIM_MODS += string_eval
3105 CHECK_BINS += check_string_eval
3106 check_string_eval_MODS =
3109 For an introduction to POSIX process control, see\\
3110 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
3111 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
3112 and of course, the relavent [[man]] pages.
3114 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
3115 [[execvp]] replaces the calling process' program with a new program.
3116 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
3117 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
3118 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
3119 The new program needs command line arguments, just like it would if you were running it from a shell.
3120 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
3121 with the final array entry being a [[NULL]] pointer.
3123 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
3124 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
3125 <<string eval subprocess definitions>>=
3126 #define SUBPROCESS "python"
3127 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
3128 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
3129 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
3131 The [[i]] option lets Python know that it should run in interactive mode.
3132 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
3133 In interactive mode, python acts on every instruction as soon as it is recieved.
3134 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.
3135 %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.
3137 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
3138 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
3139 [[fork]] returns two copies of the same program, executing the original code.
3140 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
3141 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
3143 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.
3144 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
3145 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
3146 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
3147 <<string eval pipe definitions>>=
3148 #define PIPE_READ 0 /* the end you read from */
3149 #define PIPE_WRITE 1 /* the end you write to */
3151 #define STDIN 0 /* initial index of stdin pair */
3152 #define STDOUT 2 /* initial index of stdout pair */
3155 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
3157 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.
3159 <<string eval setup declaration>>=
3160 extern void string_eval_setup(FILE **pIn, FILE **pOut);
3162 <<string eval setup definition>>=
3163 void string_eval_setup(FILE **pIn, FILE **pOut)
3166 int pfd[4]; /* file descriptors for stdin and stdout */
3168 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
3169 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
3171 pid = fork(); /* split process into two copies */
3173 if (pid == -1) { /* fork error */
3174 perror("fork error");
3176 } else if (pid == 0) { /* child */
3177 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
3178 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
3179 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
3180 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
3181 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
3182 perror("exec error"); /* exec shouldn't return */
3184 } else { /* parent */
3185 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
3186 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
3187 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
3188 if ( *pIn == NULL ) {
3189 perror("fdopen (in)");
3192 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
3193 if ( *pOut == NULL ) {
3194 perror("fdopen (out)");
3201 To use the evaluating subprocess, we just pipe in our command, and read out the results.
3202 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
3203 <<string eval function declaration>>=
3204 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
3206 <<string eval function definition>>=
3207 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
3210 rc = fprintf(in, "%s", input);
3211 assert(rc == strlen(input));
3214 alarm(1); /* set a one second timeout on the read */
3215 assert( fgets(output, buflen, out) != NULL );
3216 alarm(0); /* cancel the timeout */
3217 //fprintf(stderr, "eval: %s ----> %s", input, output);
3220 The [[alarm]] calls set and clear a timeout on the returned output.
3221 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.
3222 This protects against invalid input for which a line of output is not printed to [[stdout]].
3223 Other invalid inputs (e.g.\ those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
3224 If you are getting strange results, check your python code seperately. TODO, better error handling.
3226 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
3227 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
3228 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
3229 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
3230 <<string eval teardown declaration>>=
3231 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
3234 <<string eval teardown definition>>=
3235 void string_eval_teardown(FILE **pIn, FILE **pOut)
3240 /* redirect Python's stderr */
3241 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
3245 assert( fclose(*pIn) == 0 );
3247 assert( fclose(*pOut) == 0 );
3250 /* wait for python to exit */
3252 pid = wait(&stat_loc);
3259 if (WIFEXITED(stat_loc)) {
3260 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
3261 } else if (WIFSIGNALED(stat_loc)) {
3262 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
3267 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
3269 \subsection{String evaluation implementation}
3273 <<string eval includes>>
3274 #include "string_eval.h"
3275 <<string eval internal definitions>>
3276 <<string eval functions>>
3279 <<string eval includes>>=
3280 #include <assert.h> /* assert() */
3281 #include <stdlib.h> /* NULL */
3282 #include <stdio.h> /* fprintf(), stdout, fdopen() */
3283 #include <string.h> /* strlen() */
3284 #include <sys/types.h> /* pid_t */
3285 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
3286 #include <wait.h> /* wait() */
3289 <<string eval internal definitions>>=
3290 <<string eval subprocess definitions>>
3291 <<string eval pipe definitions>>
3294 <<string eval functions>>=
3295 <<string eval setup definition>>
3296 <<string eval function definition>>
3297 <<string eval teardown definition>>
3300 \subsection{String evaluation unit tests}
3302 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
3303 <<check-string-eval.c>>=
3304 <<string eval check includes>>
3305 <<string comparison definition and globals>>
3306 <<string eval test suite>>
3307 <<main check program>>
3310 <<string eval check includes>>=
3311 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
3312 #include <stdio.h> /* printf() */
3313 #include <assert.h> /* assert() */
3314 #include <string.h> /* strlen() */
3315 #include <signal.h> /* SIGKILL */
3317 #include "string_eval.h"
3320 <<string eval test suite>>=
3321 <<string eval tests>>
3322 <<string eval suite definition>>
3325 <<string eval suite definition>>=
3326 Suite *test_suite (void)
3328 Suite *s = suite_create ("string eval");
3329 <<string eval test case defs>>
3331 <<string eval test case add>>
3336 <<string eval tests>>=
3337 START_TEST(test_setup_teardown)
3340 string_eval_setup(&in, &out);
3341 string_eval_teardown(&in, &out);
3344 START_TEST(test_invalid_command)
3347 char input[80], output[80]={};
3348 string_eval_setup(&in, &out);
3349 sprintf(input, "print ABCDefg\n");
3350 string_eval(in, out, input, 80, output);
3351 string_eval_teardown(&in, &out);
3354 START_TEST(test_simple_eval)
3357 char input[80], output[80]={};
3358 string_eval_setup(&in, &out);
3359 sprintf(input, "print 3+4*5\n");
3360 string_eval(in, out, input, 80, output);
3361 fail_unless(STRMATCH(output,"23\n"), NULL);
3362 string_eval_teardown(&in, &out);
3365 START_TEST(test_multiple_evals)
3368 char input[80], output[80]={};
3369 string_eval_setup(&in, &out);
3370 sprintf(input, "print 3+4*5\n");
3371 string_eval(in, out, input, 80, output);
3372 fail_unless(STRMATCH(output,"23\n"), NULL);
3373 sprintf(input, "print (3**2 + 4**2)**0.5\n");
3374 string_eval(in, out, input, 80, output);
3375 fail_unless(STRMATCH(output,"5.0\n"), NULL);
3376 string_eval_teardown(&in, &out);
3379 START_TEST(test_eval_with_state)
3382 char input[80], output[80]={};
3383 string_eval_setup(&in, &out);
3384 sprintf(input, "print 3+4*5\n");
3385 fprintf(in, "A = 3\n");
3386 sprintf(input, "print A*3\n");
3387 string_eval(in, out, input, 80, output);
3388 fail_unless(STRMATCH(output,"9\n"), NULL);
3389 string_eval_teardown(&in, &out);
3393 <<string eval test case defs>>=
3394 TCase *tc_string_eval = tcase_create("string_eval");
3396 <<string eval test case add>>=
3397 tcase_add_test(tc_string_eval, test_setup_teardown);
3398 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
3399 tcase_add_test(tc_string_eval, test_simple_eval);
3400 tcase_add_test(tc_string_eval, test_multiple_evals);
3401 tcase_add_test(tc_string_eval, test_eval_with_state);
3402 suite_add_tcase(s, tc_string_eval);
3406 \section{Accelerating function evaluation}
3408 My first version-0.3 code was running very slowly.
3409 With the options suggested in the help
3410 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
3411 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
3412 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
3413 That's only 15 calls per solution, so the search algorithm seems reasonable.
3414 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
3419 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
3421 #endif /* ACCEL_K_H */
3424 <<accel k module makefile lines>>=
3425 NW_SAWSIM_MODS += accel_k
3426 #CHECK_BINS += check_accel_k
3427 check_accel_k_MODS =
3431 #include <assert.h> /* assert() */
3432 #include <stdlib.h> /* realloc(), free(), NULL */
3433 #include "global.h" /* environment_t */
3434 #include "k_model.h" /* k_func_t */
3435 #include "interp.h" /* interp_* */
3436 #include "accel_k.h"
3438 typedef struct accel_k_struct {
3439 interp_table_t *itable;
3445 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3446 static int num_accels = 0;
3447 static accel_k_t *accels=NULL;
3449 /* Wrap k in the standard f(x) acceleration form */
3450 static double k_wrap(double F, void *params)
3452 accel_k_t *p = (accel_k_t *)params;
3453 return (*p->k)(F, p->env, p->params);
3456 static int k_tol(double FA, double kA, double FB, double kB)
3459 if (FB-FA > 1e-12) {
3460 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3461 return 1; /* unnacceptable */
3463 //printf("acceptable tol\n");
3464 return 0; /* acceptable */
3468 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3471 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3472 assert(accels != NULL);
3473 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3475 accels[i].env = env;
3476 accels[i].params = params;
3483 for (i=0; i<num_accels; i++)
3484 interp_table_free(accels[i].itable);
3486 if (accels) free(accels);
3490 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3493 for (i=0; i<num_accels; i++) {
3494 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3495 /* We've already setup this function.
3496 * WARNING: we're assuming param and env are the same. */
3497 return interp_table_eval(accels[i].itable, accels+i, F);
3500 /* set up a new acceleration environment */
3501 i = add_accel_k(k, env, params);
3502 return interp_table_eval(accels[i].itable, accels+i, F);
3506 \section{Tension models}
3507 \label{sec.tension_models}
3509 TODO: tension model intro.
3510 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.
3512 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3513 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]].
3515 <<tension-model.h>>=
3516 #ifndef TENSION_MODEL_H
3517 #define TENSION_MODEL_H
3519 <<tension handler types>>
3520 <<constant tension model declarations>>
3521 <<hooke tension model declarations>>
3522 <<worm-like chain tension model declarations>>
3523 <<freely-jointed chain tension model declarations>>
3524 <<find tension definitions>>
3525 <<tension model global declarations>>
3526 #endif /* TENSION_MODEL_H */
3529 <<tension model module makefile lines>>=
3530 NW_SAWSIM_MODS += tension_model
3531 #CHECK_BINS += check_tension_model
3532 check_tension_model_MODS =
3534 <<tension model utils makefile lines>>=
3535 TENSION_MODEL_MODS = tension_model parse list tension_balance
3536 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3537 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3538 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3539 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
3540 notangle -Rtension-model-utils.c $< > $@
3541 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
3542 gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3543 clean_tension_model_utils :
3544 rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3545 @ The pipe symbol [[|]] marks the prerequisits following it (in this
3546 case, the directories) as ‘order-only’ prerequisites. The timestamp
3547 on these prerequisits does not effect whether the rules are executed.
3548 This is appropriate for directories, where we don't need to recompile
3549 after an unrelated has been added to the directory, but only when the
3550 source prerequisites change. See the [[make]] documentation for more
3552 (\url{http://www.gnu.org/software/make/manual/html_node/Prerequisite-Types.html}).
3555 \label{sec.null_tension}
3557 For unstretchable domains.
3559 <<null tension model getopt>>=
3560 {NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3563 \subsection{Constant}
3564 \label{sec.const_tension}
3566 <<constant tension functions>>=
3567 <<constant tension function>>
3568 <<constant tension parameter create/destroy functions>>
3571 <<constant tension model declarations>>=
3572 <<constant tension function declaration>>
3573 <<constant tension parameter create/destroy function declarations>>
3574 <<constant tension model global declarations>>
3578 An infinitely stretchable domain providing a constant tension.
3580 <<constant tension function declaration>>=
3581 extern double const_tension_handler(double x, void *pdata);
3583 <<constant tension function>>=
3584 double const_tension_handler(double x, void *pdata)
3586 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3587 list_t *list = data->group_tension_params;
3592 assert(list != NULL); /* empty active group?! */
3593 F = ((const_tension_param_t *)list->d)->F;
3595 fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3597 while (list != NULL) {
3598 assert(((const_tension_param_t *)list->d)->F == F);
3606 <<constant tension structure>>=
3607 typedef struct const_tension_param_struct {
3608 double F; /* tension (force) in N */
3609 } const_tension_param_t;
3613 <<constant tension parameter create/destroy function declarations>>=
3614 extern void *string_create_const_tension_param_t(char **param_strings);
3615 extern void destroy_const_tension_param_t(void *p);
3618 <<constant tension parameter create/destroy functions>>=
3619 const_tension_param_t *create_const_tension_param_t(double F)
3621 const_tension_param_t *ret
3622 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3623 assert(ret != NULL);
3628 void *string_create_const_tension_param_t(char **param_strings)
3630 return create_const_tension_param_t(atof(param_strings[0]));
3633 void destroy_const_tension_param_t(void *p)
3641 <<constant tension model global declarations>>=
3642 extern int num_const_tension_params;
3643 extern char *const_tension_param_descriptions[];
3644 extern char const_tension_param_string[];
3647 <<constant tension model globals>>=
3648 int num_const_tension_params = 1;
3649 char *const_tension_param_descriptions[] = {"tension F, N"};
3650 char const_tension_param_string[] = "0";
3653 <<constant tension model getopt>>=
3654 {&const_tension_handler, "const", "an infinitely stretchable domain with constant tension", 1, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t}
3660 <<hooke functions>>=
3662 <<hooke parameter create/destroy functions>>
3665 <<hooke tension model declarations>>=
3666 <<hooke tension function declaration>>
3667 <<hooke tension parameter create/destroy function declarations>>
3668 <<hooke tension model global declarations>>
3672 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
3673 The behavior of a series of springs $k_i$ in series is given by
3675 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
3677 For a simple proof, see Appendix \ref{sec.math_hooke}.
3679 <<hooke tension function declaration>>=
3680 extern double hooke_handler(double x, void *pdata);
3684 double hooke_handler(double x, void *pdata)
3686 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3687 list_t *list = data->group_tension_params;
3692 assert(list != NULL); /* empty active group?! */
3693 while (list != NULL) {
3694 assert( ((hooke_param_t *)list->d)->k > 0 );
3695 k += 1.0/ ((hooke_param_t *)list->d)->k;
3697 fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
3698 __FUNCTION__, ((hooke_param_t *)list->d)->k);
3704 fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N (from data %p)\n",
3705 __FUNCTION__, k, x, k*x, data);
3711 <<hooke structure>>=
3712 typedef struct hooke_param_struct {
3713 double k; /* spring constant in N/m */
3717 <<hooke tension parameter create/destroy function declarations>>=
3718 extern void *string_create_hooke_param_t(char **param_strings);
3719 extern void destroy_hooke_param_t(void *p);
3722 <<hooke parameter create/destroy functions>>=
3723 hooke_param_t *create_hooke_param_t(double k)
3725 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
3726 assert(ret != NULL);
3731 void *string_create_hooke_param_t(char **param_strings)
3733 return create_hooke_param_t(atof(param_strings[0]));
3736 void destroy_hooke_param_t(void *p)
3743 <<hooke tension model global declarations>>=
3744 extern int num_hooke_params;
3745 extern char *hooke_param_descriptions[];
3746 extern char hooke_param_string[];
3749 <<hooke tension model globals>>=
3750 int num_hooke_params = 1;
3751 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
3752 char hooke_param_string[]="0.05";
3755 <<hooke tension model getopt>>=
3756 {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}
3759 \subsection{Worm-like chain}
3762 We can model several unfolded domains with a single worm-like chain.
3763 <<worm-like chain functions>>=
3764 <<worm-like chain function>>
3765 <<worm-like chain handler>>
3766 <<worm-like chain create/destroy functions>>
3769 <<worm-like chain tension model declarations>>=
3770 <<worm-like chain handler declaration>>
3771 <<worm-like chain create/destroy function declarations>>
3772 <<worm-like chain tension model global declarations>>
3776 The combination of all unfolded domains is modeled as a worm like chain,
3777 with the force $F_{\text{WLC}}$ approximately given by
3779 F_{\text{WLC}} \approx \frac{k_B T}{p}
3780 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
3782 where $x$ is the distance between the fixed ends,
3783 $k_B$ is Boltzmann's constant,
3784 $T$ is the absolute temperature,
3785 $p$ is the persistence length, and
3786 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
3787 <<worm-like chain function>>=
3788 double wlc(double x, double T, double p, double L)
3790 double xL=0.0; /* uses global k_B */
3792 assert(T > 0); assert(p > 0); assert(L > 0);
3793 if (x >= L) return HUGE_VAL;
3794 xL = x/L; /* unitless */
3795 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
3796 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
3797 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
3800 This model assumes that each unfolded domain has the same persistence length.
3802 <<worm-like chain handler declaration>>=
3803 extern double wlc_handler(double x, void *pdata);
3806 <<worm-like chain handler>>=
3807 double wlc_handler(double x, void *pdata)
3809 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3810 list_t *list = data->group_tension_params;
3814 assert(list != NULL); /* empty active group?! */
3815 p = ((wlc_param_t *)list->d)->p;
3816 while (list != NULL) {
3817 /* enforce consistent p */
3818 assert( ((wlc_param_t *)list->d)->p == p);
3819 L += ((wlc_param_t *)list->d)->L;
3821 fprintf(stderr, "%s: WLC section %g m long in series\n",
3822 __FUNCTION__, ((wlc_param_t *)list->d)->L);
3827 fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
3828 __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
3830 return wlc(x, data->env->T, p, L);
3834 <<worm-like chain structure>>=
3835 typedef struct wlc_param_struct {
3836 double p; /* WLC persistence length (m) */
3837 double L; /* WLC contour length (m) */
3841 <<worm-like chain create/destroy function declarations>>=
3842 extern void *string_create_wlc_param_t(char **param_strings);
3843 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
3846 <<worm-like chain create/destroy functions>>=
3847 wlc_param_t *create_wlc_param_t(double p, double L)
3849 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
3850 assert(ret != NULL);
3856 void *string_create_wlc_param_t(char **param_strings)
3858 return create_wlc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3861 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
3868 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.
3869 TODO (now we copy the string before parsing, but still do this for clarity).
3871 <<valid string write code>>=
3872 char string[] = "abc";
3875 <<valid string write code>>=
3876 char *string = "abc";
3880 <<worm-like chain tension model global declarations>>=
3881 extern int num_wlc_params;
3882 extern char *wlc_param_descriptions[];
3883 extern char wlc_param_string[];
3885 <<worm-like chain tension model globals>>=
3886 int num_wlc_params = 2;
3887 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
3888 char wlc_param_string[]="0.39e-9,28.4e-9";
3891 <<worm-like chain tension model getopt>>=
3892 {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}
3894 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
3896 \subsection{Freely-jointed chain}
3899 <<freely-jointed chain functions>>=
3900 <<freely-jointed chain function>>
3901 <<freely-jointed chain handler>>
3902 <<freely-jointed chain create/destroy functions>>
3905 <<freely-jointed chain tension model declarations>>=
3906 <<freely-jointed chain handler declaration>>
3907 <<freely-jointed chain create/destroy function declarations>>
3908 <<freely-jointed chain tension model global declarations>>
3912 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
3913 The tension of a chain of $N$ such links is given by
3915 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
3917 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}.
3918 <<freely-jointed chain function>>=
3919 double langevin(double x, void *params)
3922 return 1.0/tanh(x) - 1/x;
3924 <<one dimensional bracket declaration>>
3925 <<one dimensional solve declaration>>
3926 double inv_langevin(double x)
3928 double lb=0.0, ub=1.0, ret;
3929 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
3930 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
3934 double fjc(double x, double T, double l, int N)
3936 double L = l*(double)N;
3937 if (x == 0) return 0;
3938 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
3939 return k_B*T/l * inv_langevin(x/L);
3943 <<freely-jointed chain handler declaration>>=
3944 extern double fjc_handler(double x, void *pdata);
3946 <<freely-jointed chain handler>>=
3947 double fjc_handler(double x, void *pdata)
3949 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3950 list_t *list = data->group_tension_params;
3955 assert(list != NULL); /* empty active group?! */
3956 l = ((fjc_param_t *)list->d)->l;
3957 while (list != NULL) {
3958 /* enforce consistent l */
3959 assert( ((fjc_param_t *)list->d)->l == l);
3960 N += ((fjc_param_t *)list->d)->N;
3962 fprintf(stderr, "%s: FJC section %d links long in series\n",
3963 __FUNCTION__, ((fjc_param_t *)list->d)->N);
3968 fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
3969 __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
3971 return fjc(x, data->env->T, l, N);
3975 <<freely-jointed chain structure>>=
3976 typedef struct fjc_param_struct {
3977 double l; /* FJC link length (m) */
3978 int N; /* FJC number of links */
3982 <<freely-jointed chain create/destroy function declarations>>=
3983 extern void *string_create_fjc_param_t(char **param_strings);
3984 extern void destroy_fjc_param_t(void *p);
3987 <<freely-jointed chain create/destroy functions>>=
3988 fjc_param_t *create_fjc_param_t(double l, double N)
3990 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
3991 assert(ret != NULL);
3997 void *string_create_fjc_param_t(char **param_strings)
3999 return create_fjc_param_t(atof(param_strings[0]), atof(param_strings[1]));
4002 void destroy_fjc_param_t(void *p)
4009 <<freely-jointed chain tension model global declarations>>=
4010 extern int num_fjc_params;
4011 extern char *fjc_param_descriptions[];
4012 extern char fjc_param_string[];
4015 <<freely-jointed chain tension model globals>>=
4016 int num_fjc_params = 2;
4017 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
4018 char fjc_param_string[]="0.5e-9,200";
4021 <<freely-jointed chain tension model getopt>>=
4022 {fjc_handler, "fjc", "a freely-jointed chain [F=kbT/l L^-1(x/L)]", num_fjc_params, fjc_param_descriptions, (char *)fjc_param_string, &string_create_fjc_param_t, &destroy_fjc_param_t}
4025 \subsection{Tension model implementation}
4027 <<tension-model.c>>=
4030 <<tension model includes>>
4031 #include "tension_model.h"
4032 <<tension model internal definitions>>
4033 <<tension model globals>>
4034 <<tension model functions>>
4037 <<tension model includes>>=
4038 #include <assert.h> /* assert() */
4039 #include <stdlib.h> /* NULL */
4040 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4041 #include <stdio.h> /* fprintf(), stdout */
4044 #include "tension_balance.h" /* oneD_*() */
4047 <<tension model internal definitions>>=
4048 <<constant tension structure>>
4050 <<worm-like chain structure>>
4051 <<freely-jointed chain structure>>
4054 <<tension model globals>>=
4055 <<tension handler array global>>
4056 <<constant tension model globals>>
4057 <<hooke tension model globals>>
4058 <<worm-like chain tension model globals>>
4059 <<freely-jointed chain tension model globals>>
4062 <<tension model functions>>=
4063 <<constant tension functions>>
4065 <<worm-like chain functions>>
4066 <<freely-jointed chain functions>>
4070 \subsection{Utilities}
4072 The tension models can get complicated, and users may want to reassure
4073 themselves that this portion of the input physics are functioning properly. The
4074 stand-alone program [[tension_model_utils]] demonstrates the tension model
4075 interface without the overhead of having to understand a full simulation.
4076 [[tension_model_utils]] takes command line model arguments like the full
4077 simulation, and prints $F(x)$ data to the screen for the selected model over a
4080 <<tension-model-utils.c>>=
4082 <<tension model utility includes>>
4083 <<tension model utility definitions>>
4084 <<tension model utility getopt functions>>
4086 <<tension model utility main>>
4089 <<tension model utility main>>=
4090 int main(int argc, char **argv)
4092 <<tension model getopt array definition>>
4093 tension_model_getopt_t *model = NULL;
4097 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
4098 tension_handler_data_t tdata;
4099 int num_param_args; /* for INIT_MODEL() */
4100 char **param_args; /* for INIT_MODEL() */
4102 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model,
4103 &Fmax, &Xmax, &flags);
4105 if (flags & VFLAG) {
4106 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4108 INIT_MODEL("utility", model, model->params, params);
4109 tdata.group_tension_params = NULL;
4110 push(&tdata.group_tension_params, params, 1);
4112 tdata.persist = NULL;
4113 if (model->handler == NULL) {
4114 printf("No tension function for model %s\n", model->name);
4120 printf("#Distance (m)\tForce (N)\n");
4121 for (i=0; i<=N; i++) {
4122 x = Xmax*i/(double)N;
4123 F = (*model->handler)(x, &tdata);
4124 if (F < 0 || F > Fmax) break;
4125 printf("%g\t%g\n", x, F);
4127 if (flags & VFLAG && i <= N) {
4128 /* explain exit condition */
4130 printf("Impossible force %g\n", F);
4132 printf("Reached large force limit %g > %g\n", F, Fmax);
4135 params = pop(&tdata.group_tension_params);
4136 if (model->destructor)
4137 (*model->destructor)(params);
4142 <<tension model utility includes>>=
4145 #include <string.h> /* strlen() */
4146 #include <assert.h> /* assert() */
4150 #include "tension_model.h"
4153 <<tension model utility definitions>>=
4154 <<version definition>>
4155 #define VFLAG 1 /* verbose */
4156 <<string comparison definition and globals>>
4157 <<tension model getopt definitions>>
4158 <<initialize model definition>>
4162 <<tension model utility getopt functions>>=
4163 <<index tension model>>
4164 <<tension model utility help>>
4165 <<tension model utility get options>>
4168 <<tension model utility help>>=
4169 void help(char *prog_name,
4171 int n_tension_models, tension_model_getopt_t *tension_models,
4172 int tension_model, double Fmax, double Xmax)
4175 printf("usage: %s [options]\n", prog_name);
4176 printf("Version %s\n\n", VERSION);
4177 printf("A utility for understanding the available tension models\n\n");
4178 printf("Environmental options:\n");
4179 printf("-T T\tTemperature T (currently %g K)\n", env->T);
4180 printf("-C T\tYou can also set the temperature T in Celsius\n");
4181 printf("Model options:\n");
4182 printf("-m model\tFolded tension model (currently %s)\n",
4183 tension_models[tension_model].name);
4184 printf("-a args\tFolded tension model argument string (currently %s)\n",
4185 tension_models[tension_model].params);
4186 printf("F(x) is calculated for a range of x and printed\n");
4187 printf("For example:\n");
4188 printf(" #Distance (m)\tForce (N)\n");
4189 printf(" 123.456\t7.89\n");
4191 printf("-F\tSet the maximum F value for the standard mode F(x) output (currently %g)\n", Fmax);
4192 printf("-X\tSet the maximum x value for the standard mode F(x) output (currently %g)\n", Xmax);
4193 printf("-V\tChange output to verbose mode\n");
4194 printf("-h\tPrint this help and exit\n");
4196 printf("Tension models:\n");
4197 for (i=0; i<n_tension_models; i++) {
4198 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
4199 for (j=0; j < tension_models[i].num_params ; j++)
4200 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
4201 printf("\t\tdefaults: %s\n", tension_models[i].params);
4206 <<tension model utility get options>>=
4207 void get_options(int argc, char **argv, environment_t *env,
4208 int n_tension_models, tension_model_getopt_t *tension_models,
4209 tension_model_getopt_t **model, double *Fmax, double *Xmax,
4210 unsigned int *flags)
4212 char *prog_name = NULL;
4213 char c, options[] = "T:C:m:a:F:X:Vh";
4214 int tension_model=0, num_strings;
4215 extern char *optarg;
4216 extern int optind, optopt, opterr;
4220 /* setup defaults */
4222 prog_name = argv[0];
4223 env->T = 300.0; /* K */
4227 *model = tension_models;
4229 while ((c=getopt(argc, argv, options)) != -1) {
4231 case 'T': env->T = atof(optarg); break;
4232 case 'C': env->T = atof(optarg)+273.15; break;
4234 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
4235 *model = tension_models+tension_model;
4238 tension_models[tension_model].params = optarg;
4240 case 'F': *Fmax = atof(optarg); break;
4241 case 'X': *Xmax = atof(optarg); break;
4242 case 'V': *flags |= VFLAG; break;
4244 fprintf(stderr, "unrecognized option '%c'\n", optopt);
4245 /* fall through to default case */
4247 help(prog_name, env, n_tension_models, tension_models, tension_model, *Fmax, *Xmax);
4256 \section{Unfolding rate models}
4257 \label{sec.k_models}
4259 We have two main choices for determining $k$: Bell's model, which assumes the
4260 folded domains exist in a single `folded' state and make instantaneous,
4261 stochastic transitions over a free energy barrier; and Kramers' model, which
4262 treats the unfolding as a thermalized diffusion process.
4263 We deal with these options like we dealt with the different tension models: by bundling all model-specific
4264 parameters into structures, with model specific functions for determining $k$.
4266 <<k func definition>>=
4267 typedef double k_func_t(double F, environment_t *env, void *params);
4270 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
4271 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]].
4273 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
4274 because \LaTeX\ doesn't like underscores outside math mode.
4275 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
4276 You could use spaces instead (`\verb+ +'),
4277 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
4278 which I don't like as much.
4284 <<k func definition>>
4285 <<null k declarations>>
4286 <<const k declarations>>
4287 <<bell k declarations>>
4288 <<kbell k declarations>>
4289 <<kramers k declarations>>
4290 <<kramers integ k declarations>>
4291 #endif /* K_MODEL_H */
4294 <<k model module makefile lines>>=
4295 NW_SAWSIM_MODS += k_model
4296 CHECK_BINS += check_k_model
4297 check_k_model_MODS = parse string_eval
4299 [[check_k_model]] also depends on the parse module.
4301 <<k model utils makefile lines>>=
4302 K_MODEL_MODS = k_model parse string_eval
4303 K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
4304 $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h)
4305 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
4306 notangle -Rk-model-utils.c $< > $@
4307 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
4308 gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
4309 clean_k_model_utils :
4310 rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
4316 For domains that are never folded.
4318 <<null k declarations>>=
4320 <<null k functions>>=
4325 <<null k model getopt>>=
4326 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
4329 \subsection{Constant rate model}
4332 This is a pretty boring model, but it is usefull for testing the engine.
4336 where $k_0$ is the constant unfolding rate.
4338 <<const k functions>>=
4339 <<const k function>>
4340 <<const k structure create/destroy functions>>
4343 <<const k declarations>>=
4344 <<const k function declaration>>
4345 <<const k structure create/destroy function declarations>>
4346 <<const k global declarations>>
4350 <<const k structure>>=
4351 typedef struct const_k_param_struct {
4352 double knot; /* transition rate k_0 (frac population per s) */
4356 <<const k function declaration>>=
4357 double const_k(double F, environment_t *env, void *const_k_params);
4359 <<const k function>>=
4360 double const_k(double F, environment_t *env, void *const_k_params)
4361 { /* Returns the rate constant k in frac pop per s. */
4362 const_k_param_t *p = (const_k_param_t *) const_k_params;
4364 assert(p->knot > 0);
4369 <<const k structure create/destroy function declarations>>=
4370 void *string_create_const_k_param_t(char **param_strings);
4371 void destroy_const_k_param_t(void *p);
4374 <<const k structure create/destroy functions>>=
4375 const_k_param_t *create_const_k_param_t(double knot)
4377 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
4378 assert(ret != NULL);
4383 void *string_create_const_k_param_t(char **param_strings)
4385 return create_const_k_param_t(atof(param_strings[0]));
4388 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
4395 <<const k global declarations>>=
4396 extern int num_const_k_params;
4397 extern char *const_k_param_descriptions[];
4398 extern char const_k_param_string[];
4401 <<const k globals>>=
4402 int num_const_k_params = 1;
4403 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
4404 char const_k_param_string[]="1";
4407 <<const k model getopt>>=
4408 {"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}
4411 \subsection{Bell's model}
4414 Quantitatively, Bell's model gives $k$ as
4416 k = k_0 \cdot e^{F dx / k_B T} \;,
4418 where $k_0$ is the unforced unfolding rate,
4419 $F$ is the applied unfolding force,
4420 $dx$ is the distance to the transition state, and
4421 $k_B T$ is the thermal energy\citep{schlierf06}.
4423 <<bell k functions>>=
4425 <<bell k structure create/destroy functions>>
4428 <<bell k declarations>>=
4429 <<bell k function declaration>>
4430 <<bell k structure create/destroy function declarations>>
4431 <<bell k global declarations>>
4435 <<bell k structure>>=
4436 typedef struct bell_param_struct {
4437 double knot; /* transition rate k_0 (frac population per s) */
4438 double dx; /* `distance to transition state' (nm) */
4442 <<bell k function declaration>>=
4443 double bell_k(double F, environment_t *env, void *bell_params);
4445 <<bell k function>>=
4446 double bell_k(double F, environment_t *env, void *bell_params)
4447 { /* Returns the rate constant k in frac pop per s.
4448 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4449 * uses global k_B in J/K */
4450 bell_param_t *p = (bell_param_t *) bell_params;
4451 assert(F >= 0); assert(env->T > 0);
4453 assert(p->knot > 0); //assert(p->dx >= 0); /* p < 0 for refolding reactions, etc. */
4455 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
4456 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
4457 p->knot * exp(F*p->dx / (k_B*env->T) ));
4458 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
4460 return p->knot * exp(F*p->dx / (k_B*env->T) );
4464 <<bell k structure create/destroy function declarations>>=
4465 void *string_create_bell_param_t(char **param_strings);
4466 void destroy_bell_param_t(void *p);
4469 <<bell k structure create/destroy functions>>=
4470 bell_param_t *create_bell_param_t(double knot, double dx)
4472 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
4473 assert(ret != NULL);
4479 void *string_create_bell_param_t(char **param_strings)
4481 return create_bell_param_t(atof(param_strings[0]), atof(param_strings[1]));
4484 void destroy_bell_param_t(void *p /* bell_param_t * */)
4491 <<bell k global declarations>>=
4492 extern int num_bell_params;
4493 extern char *bell_param_descriptions[];
4494 extern char bell_param_string[];
4498 int num_bell_params = 2;
4499 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4500 char bell_param_string[]="3.3e-4,0.25e-9";
4503 <<bell k model getopt>>=
4504 {"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}
4506 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4507 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4510 \subsection{Linker stiffness corrected Bell model}
4513 Walton et. al showed that the Bell model constant force approximation
4514 doesn't hold for biotin-streptavadin binding, and I extended their
4515 results to I27 for the 2009 Biophysical Society Annual
4516 Meeting\cite{walton08,king09}. More details to follow when I get done
4517 with the conference\ldots
4519 We adjust Bell's model to
4521 k = k_0 \cdot e^{\frac{F dx - \frac{1}{2} \kappa dx^2}{k_B T}} \;,
4523 where $k_0$ is the unforced unfolding rate,
4524 $F$ is the applied unfolding force,
4525 $\kappa$ is the effective spring constant,
4526 $dx$ is the distance to the transition state, and
4527 $k_B T$ is the thermal energy\citep{schlierf06}.
4529 <<kbell k functions>>=
4530 <<kbell k function>>
4531 <<kbell k structure create/destroy functions>>
4534 <<kbell k declarations>>=
4535 <<kbell k function declaration>>
4536 <<kbell k structure create/destroy function declarations>>
4537 <<kbell k global declarations>>
4541 <<kbell k structure>>=
4542 typedef struct kbell_param_struct {
4543 double knot; /* transition rate k_0 (frac population per s) */
4544 double dx; /* `distance to transition state' (nm) */
4548 <<kbell k function declaration>>=
4549 double kbell_k(double F, environment_t *env, void *kbell_params);
4551 <<kbell k function>>=
4552 double kbell_k(double F, environment_t *env, void *kbell_params)
4553 { /* Returns the rate constant k in frac pop per s.
4554 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4555 * uses global k_B in J/K */
4556 kbell_param_t *p = (kbell_param_t *) kbell_params;
4557 assert(F >= 0); assert(env->T > 0); assert(env->stiffness >= 0);
4559 assert(p->knot > 0); assert(p->dx >= 0);
4560 return p->knot * exp((F*p->dx - 0.5*env->stiffness*(p->dx)*(p->dx))/(k_B*env->T) );
4564 <<kbell k structure create/destroy function declarations>>=
4565 void *string_create_kbell_param_t(char **param_strings);
4566 void destroy_kbell_param_t(void *p);
4569 <<kbell k structure create/destroy functions>>=
4570 kbell_param_t *create_kbell_param_t(double knot, double dx)
4572 kbell_param_t *ret = (kbell_param_t *) malloc(sizeof(kbell_param_t));
4573 assert(ret != NULL);
4579 void *string_create_kbell_param_t(char **param_strings)
4581 return create_kbell_param_t(atof(param_strings[0]), atof(param_strings[1]));
4584 void destroy_kbell_param_t(void *p /* kbell_param_t * */)
4591 <<kbell k global declarations>>=
4592 extern int num_kbell_params;
4593 extern char *kbell_param_descriptions[];
4594 extern char kbell_param_string[];
4597 <<kbell k globals>>=
4598 int num_kbell_params = 2;
4599 char *kbell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4600 char kbell_param_string[]="3.3e-4,0.25e-9";
4603 <<kbell k model getopt>>=
4604 {"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}
4606 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4607 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4610 \subsection{Kramer's model}
4613 Kramer's model gives $k$ as
4615 % \frac{1}{k} = \frac{1}{D}
4616 % \int_{x_\text{min}}^{x_\text{max}}
4617 % e^{\frac{-U_F(x)}{k_B T}}
4618 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4621 %where $D$ is the diffusion constant,
4622 %$U_F(x)$ is the free energy along the unfolding coordinate, and
4623 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4624 % before the minimum and after the maximum, respectively \citep{schlierf06}.
4626 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
4628 where $D$ is the diffusion constant,
4630 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
4632 are length scales where
4633 $x_c(F)$ is the location of the bound state and
4634 $x_{ts}(F)$ is the location of the transition state,
4635 $E(F,x)$ is the free energy, and
4636 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
4637 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
4638 \citet{evans97} Eqn.~3).
4640 <<kramers k functions>>=
4641 <<kramers k function>>
4642 <<kramers k structure create/destroy functions>>
4645 <<kramers k declarations>>=
4646 <<kramers k function declaration>>
4647 <<kramers k structure create/destroy function declarations>>
4648 <<kramers k global declarations>>
4652 <<kramers k structure>>=
4653 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
4654 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
4656 typedef struct kramers_param_struct {
4657 double D; /* diffusion rate (um/s) */
4664 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
4665 //kramers_x_func_t *xts; /* function returning transition x (nm) */
4666 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
4667 //kramers_E_func_t *E; /* function returning E (J) */
4668 //double *E_params; /* parametrize the energy landscape */
4669 //int n_E_params; /* length of E_params array */
4673 <<kramers k function declaration>>=
4674 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
4675 extern double kramers_k(double F, environment_t *env, void *kramers_params);
4677 <<kramers k function>>=
4678 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
4680 static char input[160]={0};
4681 static char output[80]={0};
4682 /* setup the environment */
4683 fprintf(in, "F = %g; T = %g\n", F, T);
4684 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4685 string_eval(in, out, input, 80, output);
4686 return atof(output);
4689 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
4691 static char input[160]={0};
4692 static char output[80]={0};
4693 /* setup the environment */
4694 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
4695 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4696 string_eval(in, out, input, 80, output);
4697 return atof(output);
4700 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
4702 kramers_param_t *p = (kramers_param_t *) kramers_params;
4703 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
4706 double kramers_k(double F, environment_t *env, void *kramers_params)
4707 { /* Returns the rate constant k in frac pop per s.
4708 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4709 * uses global k_B in J/K */
4710 kramers_param_t *p = (kramers_param_t *) kramers_params;
4711 double kbT, xc, xts, lc, lts, Eb;
4712 assert(F >= 0); assert(env->T > 0);
4715 assert(p->xc != NULL); assert(p->xts != NULL);
4716 assert(p->ddEdxx != NULL); assert(p->E != NULL);
4717 assert(p->in != NULL); assert(p->out != NULL);
4719 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
4720 if (xc == -1.0) return -1.0; /* error (Too much force) */
4721 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
4722 if (xc == -1.0) return -1.0; /* error (Too much force) */
4723 lc = sqrt(2.0*M_PI*kbT /
4724 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
4725 lts = sqrt(-2.0*M_PI*kbT/
4726 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
4727 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
4728 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
4729 //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));
4730 return p->D/(lc*lts) * exp(-Eb/kbT);
4734 <<kramers k structure create/destroy function declarations>>=
4735 void *string_create_kramers_param_t(char **param_strings);
4736 void destroy_kramers_param_t(void *p);
4739 <<kramers k structure create/destroy functions>>=
4740 kramers_param_t *create_kramers_param_t(double D,
4741 char *xc_fn, char *xts_fn,
4742 char *E_fn, char *ddEdxx_fn,
4743 char *global_define_string)
4744 // kramers_x_func_t *xc, kramers_x_func_t *xts,
4745 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
4746 // double *E_params, int n_E_params)
4748 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
4749 assert(ret != NULL);
4750 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4751 assert(ret->xc != NULL);
4752 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4753 assert(ret->xts != NULL);
4754 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
4755 assert(ret->E != NULL);
4756 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
4757 assert(ret->ddEdxx != NULL);
4759 strcpy(ret->xc, xc_fn);
4760 strcpy(ret->xts, xts_fn);
4761 strcpy(ret->E, E_fn);
4762 strcpy(ret->ddEdxx, ddEdxx_fn);
4763 //ret->E_params = E_params;
4764 //ret->n_E_params = n_E_params;
4765 string_eval_setup(&ret->in, &ret->out);
4766 fprintf(ret->in, "kB = %g\n", k_B);
4767 fprintf(ret->in, "%s\n", global_define_string);
4771 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
4772 void *string_create_kramers_param_t(char **param_strings)
4774 return create_kramers_param_t(atof(param_strings[0]),
4782 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
4784 kramers_param_t *pk = (kramers_param_t *)p;
4786 string_eval_teardown(&pk->in, &pk->out);
4788 // free(pk->E_params);
4794 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
4795 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.
4796 We follow \citet{shillcock98} and use
4798 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
4799 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
4802 where TODO, variable meanings.%\citep{shillcock98}.
4803 The first and second derivatives of $E(F,E_b,x,x_s)$ are
4805 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\\
4806 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
4808 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
4809 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
4810 The roots are therefor at
4812 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
4813 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
4814 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
4817 <<kramers k global declarations>>=
4818 extern int num_kramers_params;
4819 extern char *kramers_param_descriptions[];
4820 extern char kramers_param_string[];
4823 <<kramers k globals>>=
4824 int num_kramers_params = 6;
4825 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"};
4826 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)}";
4828 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
4829 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
4830 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
4831 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.
4832 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
4833 It works so long as [[val_1]] is not [[false]].
4835 <<kramers k model getopt>>=
4836 {"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}
4839 \subsection{Kramer's model (natural cubic splines)}
4840 \label{sec.kramers_integ}
4842 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
4843 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
4844 \citet{schlierf06} Eqn.~2)
4846 \frac{1}{k} = \frac{1}{D}
4847 \int_{x_\text{min}}^{x_\text{max}}
4848 e^{\frac{U_F(x)}{k_B T}}
4849 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4852 where $D$ is the diffusion constant,
4853 $U_F(x)$ is the free energy along the unfolding coordinate, and
4854 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4855 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
4857 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
4858 interpolating between them with cubic splines.
4859 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
4861 <<kramers integ k functions>>=
4862 <<spline functions>>
4863 <<kramers integ A k functions>>
4864 <<kramers integ B k functions>>
4865 <<kramers integ k function>>
4866 <<kramers integ k structure create/destroy functions>>
4869 <<kramers integ k declarations>>=
4870 <<kramers integ k function declaration>>
4871 <<kramers integ k structure create/destroy function declarations>>
4872 <<kramers integ k global declarations>>
4876 <<kramers integ k structure>>=
4877 typedef double func_t(double x, void *params);
4878 typedef struct kramers_integ_param_struct {
4879 double D; /* diffusion rate (um/s) */
4880 double x_min; /* integration bounds */
4882 func_t *E; /* function returning E (J) */
4883 void *E_params; /* parametrize the energy landscape */
4884 destroy_data_func_t *destroy_E_params;
4886 double F; /* for passing into GSL evaluations */
4888 } kramers_integ_param_t;
4891 <<spline param structure>>=
4892 typedef struct spline_param_struct {
4893 int n_knots; /* natural cubic spline knots for U(x) */
4896 gsl_spline *spline; /* GSL spline parameters */
4897 gsl_interp_accel *acc; /* GSL spline acceleration data */
4901 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
4903 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4905 <<kramers integ A k functions>>=
4906 double kramers_integ_k_integrandA(double x, void *params)
4908 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4909 double U = (*p->E)(x, p->E_params);
4910 double Fx = p->F * x;
4911 double kbT = k_B * p->env->T;
4912 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
4913 return exp(-(U-Fx)/kbT);
4915 double kramers_integ_k_integralA(double x, void *params)
4917 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4919 double result, abserr;
4920 size_t neval; /* for qng */
4921 static gsl_integration_workspace *w = NULL;
4922 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4923 f.function = &kramers_integ_k_integrandA;
4925 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4926 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4927 //fprintf(stderr, "integralA = %g\n", result);
4933 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
4934 \int_{x_\text{min}}^{x_\text{max}}
4935 e^{\frac{U_F(x)}{k_B T}}
4936 \text{Integral}_A(x)
4939 <<kramers integ B k functions>>=
4940 double kramers_integ_k_integrandB(double x, void *params)
4942 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4943 double U = (*p->E)(x, p->E_params);
4944 double Fx = p->F * x;
4945 double kbT = k_B * p->env->T;
4946 double intA = kramers_integ_k_integralA(x, params);
4947 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
4948 return exp((U-Fx)/kbT)*intA;
4950 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
4952 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4954 double result, abserr;
4955 size_t neval; /* for qng */
4956 static gsl_integration_workspace *w = NULL;
4957 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4958 f.function = &kramers_integ_k_integrandB;
4962 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4963 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4964 //fprintf(stderr, "integralB = %g\n", result);
4969 With the integrals taken care of, evaluating $k$ is simply
4971 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
4973 <<kramers integ k function declaration>>=
4974 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
4976 <<kramers integ k function>>=
4977 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
4978 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
4979 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
4980 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
4984 <<kramers integ k structure create/destroy function declarations>>=
4985 void *string_create_kramers_integ_param_t(char **param_strings);
4986 void destroy_kramers_integ_param_t(void *p);
4989 <<kramers integ k structure create/destroy functions>>=
4990 kramers_integ_param_t *create_kramers_integ_param_t(double D,
4991 double xmin, double xmax, func_t *E,
4993 destroy_data_func_t *destroy_E_params)
4995 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
4996 assert(ret != NULL);
5001 ret->E_params = E_params;
5002 ret->destroy_E_params = destroy_E_params;
5006 void *string_create_kramers_integ_param_t(char **param_strings)
5008 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
5009 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
5010 void *E_params = string_create_spline_param_t(param_strings+1);
5011 return create_kramers_integ_param_t(atof(param_strings[0]),
5012 atof(param_strings[2]),
5013 atof(param_strings[3]),
5014 &spline_eval, E_params,
5015 destroy_spline_param_t);
5018 void destroy_kramers_integ_param_t(void *params)
5020 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
5023 (*p->destroy_E_params)(p->E_params);
5029 Finally we have the GSL spline wrappers:
5030 <<spline functions>>=
5032 <<create/destroy spline>>
5035 <<spline function>>=
5036 double spline_eval(double x, void *spline_params)
5038 spline_param_t *p = (spline_param_t *)(spline_params);
5039 return gsl_spline_eval(p->spline, x, p->acc);
5043 <<create/destroy spline>>=
5044 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
5046 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
5047 assert(ret != NULL);
5048 ret->n_knots = n_knots;
5051 ret->acc = gsl_interp_accel_alloc();
5052 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
5053 gsl_spline_init(ret->spline, x, y, n_knots);
5057 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
5058 void *string_create_spline_param_t(char **param_strings)
5062 double *x=NULL, *y=NULL;
5063 /* break into ordered pairs */
5064 parse_list_string(param_strings[0], SEP, '(', ')',
5065 &num_knots, &knot_coords);
5066 x = (double *)malloc(sizeof(double)*num_knots);
5068 y = (double *)malloc(sizeof(double)*num_knots);
5070 for (i=0; i<num_knots; i++) {
5071 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
5072 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
5075 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
5077 free_parsed_list(num_knots, knot_coords);
5078 return create_spline_param_t(num_knots, x, y);
5081 void destroy_spline_param_t(void *params)
5083 spline_param_t *p = (spline_param_t *)params;
5086 gsl_spline_free(p->spline);
5088 gsl_interp_accel_free(p->acc);
5098 <<kramers integ k global declarations>>=
5099 extern int num_kramers_integ_params;
5100 extern char *kramers_integ_param_descriptions[];
5101 extern char kramers_integ_param_string[];
5104 <<kramers integ k globals>>=
5105 int num_kramers_integ_params = 4;
5106 char *kramers_integ_param_descriptions[] = {
5107 "diffusion D, m^2 / s",
5108 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
5109 "starting integration bound x_min, m",
5110 "ending integration bound x_max, m"};
5111 //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";
5112 //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";
5113 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";
5114 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
5115 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
5116 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
5119 <<kramers integ k model getopt>>=
5120 {"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}
5122 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
5123 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
5125 \subsection{Unfolding model implementation}
5129 <<k model includes>>
5130 #include "k_model.h"
5131 <<k model internal definitions>>
5133 <<k model functions>>
5136 <<k model includes>>=
5137 #include <assert.h> /* assert() */
5138 #include <stdlib.h> /* NULL, malloc() */
5139 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
5140 #include <stdio.h> /* fprintf(), stdout */
5141 #include <string.h> /* strlen(), strcpy() */
5142 #include <gsl/gsl_integration.h>
5143 #include <gsl/gsl_spline.h>
5148 <<k model internal definitions>>=
5149 <<const k structure>>
5150 <<bell k structure>>
5151 <<kbell k structure>>
5152 <<kramers k structure>>
5153 <<kramers integ k structure>>
5154 <<spline param structure>>
5157 <<k model globals>>=
5162 <<kramers k globals>>
5163 <<kramers integ k globals>>
5166 <<k model functions>>=
5167 <<null k functions>>
5168 <<const k functions>>
5169 <<bell k functions>>
5170 <<kbell k functions>>
5171 <<kramers k functions>>
5172 <<kramers integ k functions>>
5175 \subsection{Unfolding model unit tests}
5177 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
5178 <<check-k-model.c>>=
5179 <<k model check includes>>
5180 <<check relative error>>
5182 <<k model test suite>>
5183 <<main check program>>
5186 <<k model check includes>>=
5187 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
5188 #include <stdio.h> /* sprintf() */
5189 #include <assert.h> /* assert() */
5190 #include <math.h> /* exp() */
5193 #include "k_model.h"
5196 <<k model test suite>>=
5199 <<k model suite definition>>
5202 <<k model suite definition>>=
5203 Suite *test_suite (void)
5205 Suite *s = suite_create ("k model");
5206 <<const k test case defs>>
5207 <<bell k test case defs>>
5209 <<const k test case adds>>
5210 <<bell k test case adds>>
5215 \subsubsection{Constant}
5217 <<const k test case defs>>=
5218 TCase *tc_const_k = tcase_create("Constant k");
5221 <<const k test case adds>>=
5222 tcase_add_test(tc_const_k, test_const_k_create_destroy);
5223 tcase_add_test(tc_const_k, test_const_k_over_range);
5224 suite_add_tcase(s, tc_const_k);
5229 START_TEST(test_const_k_create_destroy)
5231 char *knot[] = {"1","2","3","4","5","6"};
5232 char *params[] = {knot[0]};
5235 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5236 params[0] = knot[i];
5237 p = string_create_const_k_param_t(params);
5238 destroy_const_k_param_t(p);
5243 START_TEST(test_const_k_over_range)
5245 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
5246 char *knot[] = {"1","2","3","4","5","6"};
5247 char *params[] = {knot[0]};
5250 char param_string[80];
5253 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5254 params[0] = knot[i];
5255 p = string_create_const_k_param_t(params);
5256 for ( F=Fm; F<FM; F+=dF ) {
5257 fail_unless(const_k(F, &env, p)==atof(knot[i]),
5258 "const_k(%g, %g, &{%s}) = %g != %s",
5259 F, T, knot[i], const_k(F, &env, p), knot[i]);
5261 destroy_const_k_param_t(p);
5267 \subsubsection{Bell}
5269 <<bell k test case defs>>=
5270 TCase *tc_bell = tcase_create("Bell's k");
5273 <<bell k test case adds>>=
5274 tcase_add_test(tc_bell, test_bell_k_create_destroy);
5275 tcase_add_test(tc_bell, test_bell_k_at_zero);
5276 tcase_add_test(tc_bell, test_bell_k_at_one);
5277 suite_add_tcase(s, tc_bell);
5281 START_TEST(test_bell_k_create_destroy)
5283 char *knot[] = {"1","2","3","4","5","6"};
5285 char *params[] = {knot[0], dx};
5288 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5289 params[0] = knot[i];
5290 p = string_create_bell_param_t(params);
5291 destroy_bell_param_t(p);
5296 START_TEST(test_bell_k_at_zero)
5298 double F=0.0, T=300.0;
5299 char *knot[] = {"1","2","3","4","5","6"};
5301 char *params[] = {knot[0], dx};
5304 char param_string[80];
5307 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5308 params[0] = knot[i];
5309 p = string_create_bell_param_t(params);
5310 fail_unless(bell_k(F, &env, p)==atof(knot[i]),
5311 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
5312 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
5313 destroy_bell_param_t(p);
5318 START_TEST(test_bell_k_at_one)
5321 char *knot[] = {"1","2","3","4","5","6"};
5323 char *params[] = {knot[0], dx};
5324 double F= k_B*T/atof(dx);
5327 char param_string[80];
5330 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
5331 params[0] = knot[i];
5332 p = string_create_bell_param_t(params);
5333 CHECK_ERR(1e-6, atof(knot[i])*exp(1.0), bell_k(F, &env, p));
5334 //fail_unless(bell_k(F, &env, p)==atof(knot[i])*exp(1.0),
5335 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
5336 // F, T, knot[i], dx, bell_k(F, &env, p), atof(knot[i])*exp(1.0));
5337 destroy_bell_param_t(p);
5343 <<kramers k tests>>=
5346 <<kramers k test case def>>=
5349 <<kramers k test case add>>=
5352 <<k model function tests>>=
5353 <<check relative error>>
5355 START_TEST(test_something)
5357 double k=5, x=3, last_x=2.0, F;
5358 one_dim_fn_t *handlers[] = {&hooke};
5359 void *data[] = {&k};
5360 double xi[] = {0.0};
5362 int new_active_groups = 1;
5363 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
5364 fail_unless(F = k*x, NULL);
5369 \subsection{Utilities}
5371 The unfolding models can get complicated, and are sometimes hard to explain to others.
5372 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
5373 the overhead of having to understand a full simulation.
5374 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
5375 or special mode, where the operation depends on the specific model selected.
5377 <<k-model-utils.c>>=
5379 <<k model utility includes>>
5380 <<k model utility definitions>>
5381 <<k model utility getopt functions>>
5382 <<k model utility multi model E>>
5383 <<k model utility main>>
5386 <<k model utility main>>=
5387 int main(int argc, char **argv)
5389 <<k model getopt array definition>>
5390 k_model_getopt_t *model = NULL;
5395 int num_param_args; /* for INIT_MODEL() */
5396 char **param_args; /* for INIT_MODEL() */
5398 double special_xmin;
5399 double special_xmax;
5400 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
5401 &Fmax, &special_xmin, &special_xmax, &flags);
5402 if (flags & VFLAG) {
5403 printf("#initializing model %s with parameters %s\n", model->name, model->params);
5405 INIT_MODEL("utility", model, model->params, params);
5408 if (model->k == NULL) {
5409 printf("No k function for model %s\n", model->name);
5413 printf("Fmax = %g <= 0\n", Fmax);
5419 printf("#F (N)\tk (%% pop. per s)\n");
5420 for (i=0; i<=N; i++) {
5421 F = Fmax*i/(double)N;
5422 k = (*model->k)(F, &env, params);
5424 printf("%g\t%g\n", F, k);
5429 if (model->k == NULL || model->k == &bell_k) {
5430 printf("No special mode for model %s\n", model->name);
5433 if (special_xmax <= special_xmin) {
5434 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
5438 double Xrng=(special_xmax-special_xmin), x, E;
5440 printf("#x (m)\tE (J)\n");
5441 for (i=0; i<=N; i++) {
5442 x = special_xmin + Xrng*i/(double)N;
5443 E = multi_model_E(model->k, params, &env, x);
5444 printf("%g\t%g\n", x, E);
5451 if (model->destructor)
5452 (*model->destructor)(params);
5457 <<k model utility multi model E>>=
5458 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
5460 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
5462 if (k == kramers_integ_k)
5463 E = (*p->E)(x, p->E_params);
5465 E = kramers_E(0, env, model_params, x);
5471 <<k model utility includes>>=
5474 #include <string.h> /* strlen() */
5475 #include <assert.h> /* assert() */
5478 #include "string_eval.h"
5479 #include "k_model.h"
5482 <<k model utility definitions>>=
5483 <<version definition>>
5484 #define VFLAG 1 /* verbose */
5485 enum mode_t {M_K_OF_F, M_SPECIAL};
5486 <<string comparison definition and globals>>
5487 <<k model getopt definitions>>
5488 <<initialize model definition>>
5489 <<kramers k structure>>
5490 <<kramers integ k structure>>
5494 <<k model utility getopt functions>>=
5496 <<k model utility help>>
5497 <<k model utility get options>>
5500 <<k model utility help>>=
5501 void help(char *prog_name,
5503 int n_k_models, k_model_getopt_t *k_models,
5504 int k_model, double Fmax, double special_xmin, double special_xmax)
5507 printf("usage: %s [options]\n", prog_name);
5508 printf("Version %s\n\n", VERSION);
5509 printf("A utility for understanding the available unfolding models\n\n");
5510 printf("Environmental options:\n");
5511 printf("-T T\tTemperature T (currently %g K)\n", env->T);
5512 printf("-C T\tYou can also set the temperature T in Celsius\n");
5513 printf("Model options:\n");
5514 printf("-k model\tTransition rate model (currently %s)\n",
5515 k_models[k_model].name);
5516 printf("-K args\tTransition rate model argument string (currently %s)\n",
5517 k_models[k_model].params);
5518 printf("There are two output modes. In standard mode, k(F) is printed\n");
5519 printf("For example:\n");
5520 printf(" #Force (N)\tk (% pop. per s)\n");
5521 printf(" 123.456\t7.89\n");
5523 printf("In special mode, the output depends on the model.\n");
5524 printf("For models defining an energy function E(x), that function is printed\n");
5525 printf(" #Distance (m)\tE (J)\n");
5526 printf(" 0.0012\t0.0034\n");
5528 printf("-m\tChange output to standard mode\n");
5529 printf("-M\tChange output to special mode\n");
5530 printf("-F\tSet the maximum F value for the standard mode k(F) output (currently %g)\n", Fmax);
5531 printf("-x\tSet the minimum x value for the special mode E(x) output (currently %g)\n", special_xmin);
5532 printf("-X\tSet the maximum x value for the special mode E(x) output (currently %g)\n", special_xmax);
5533 printf("-V\tChange output to verbose mode\n");
5534 printf("-h\tPrint this help and exit\n");
5536 printf("Unfolding rate models:\n");
5537 for (i=0; i<n_k_models; i++) {
5538 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
5539 for (j=0; j < k_models[i].num_params ; j++)
5540 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
5541 printf("\t\tdefaults: %s\n", k_models[i].params);
5546 <<k model utility get options>>=
5547 void get_options(int argc, char **argv, environment_t *env,
5548 int n_k_models, k_model_getopt_t *k_models,
5549 enum mode_t *mode, k_model_getopt_t **model,
5550 double *Fmax, double *special_xmin, double *special_xmax,
5551 unsigned int *flags)
5553 char *prog_name = NULL;
5554 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
5556 extern char *optarg;
5557 extern int optind, optopt, opterr;
5561 /* setup defaults */
5563 prog_name = argv[0];
5564 env->T = 300.0; /* K */
5569 *special_xmax = 1e-8;
5570 *special_xmin = 0.0;
5572 while ((c=getopt(argc, argv, options)) != -1) {
5574 case 'T': env->T = atof(optarg); break;
5575 case 'C': env->T = atof(optarg)+273.15; break;
5577 k_model = index_k_model(n_k_models, k_models, optarg);
5578 *model = k_models+k_model;
5581 k_models[k_model].params = optarg;
5583 case 'm': *mode = M_K_OF_F; break;
5584 case 'M': *mode = M_SPECIAL; break;
5585 case 'F': *Fmax = atof(optarg); break;
5586 case 'x': *special_xmin = atof(optarg); break;
5587 case 'X': *special_xmax = atof(optarg); break;
5588 case 'V': *flags |= VFLAG; break;
5590 fprintf(stderr, "unrecognized option '%c'\n", optopt);
5591 /* fall through to default case */
5593 help(prog_name, env, n_k_models, k_models, k_model, *Fmax, *special_xmin, *special_xmax);
5602 \section{\LaTeX\ documentation}
5604 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
5605 The comment blocks are extracted (with nicely formatted code blocks), using
5606 <<latex makefile lines>>=
5607 $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5608 noweave -latex -index -delay $< > $@
5609 $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR)
5613 <<latex makefile lines>>=
5614 $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
5616 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5617 $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
5618 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5619 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5620 mv $(BUILD_DIR)/sawsim.pdf $@
5622 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
5623 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
5624 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
5626 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
5627 <<latex makefile lines>>=
5629 rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
5630 $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
5631 $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.tex \
5632 $(BUILD_DIR)/sawsim.bib
5633 clean_latex : semi-clean_latex
5634 rm -f $(DOC_DIR)/sawsim.pdf
5639 [[make]] is a common utility on *nix systems for managing dependencies while building software.
5640 In this case, we have several \emph{targets} that we'd like to build:
5641 the main simulation program \prog;
5642 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
5643 the documentation [[sawsim.pdf]];
5644 and the [[Makefile]] itself.
5645 Besides the generated files, there is the utility target
5646 [[clean]] for removing all generated files except [[Makefile]].
5648 % [[dist]] for generating a distributable tar file.
5650 Extract the makefile with
5651 `[[notangle -Rmakefile src/sawsim.nw | sed 's/ /\t/' > Makefile]]'.
5652 The sed is needed because notangle replaces tabs with eight spaces,
5653 but [[make]] needs tabs.
5655 # Decide what directories to put things in
5660 # Note: Cannot use, for example, `./build', because make eats the `./'
5661 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5663 # Modules (X.c, X.h) defined in the noweb file
5666 # Modules defined outside the noweb file
5667 FREE_SAWSIM_MODS = interp tavl
5669 <<list module makefile lines>>
5670 <<tension balance module makefile lines>>
5671 <<tension model module makefile lines>>
5672 <<k model module makefile lines>>
5673 <<parse module makefile lines>>
5674 <<string eval module makefile lines>>
5675 <<accel k module makefile lines>>
5677 SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS)
5679 # Everything needed for sawsim under one roof. sawsim.c must come first
5680 SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \
5681 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h)
5682 # Libraries needed to compile sawsim
5683 LIBS = gsl gslcblas m
5684 CHECK_LIBS = gsl gslcblas m check
5685 # The non-check binaries generated
5686 BINS = sawsim tension_model_utils k_model_utils sawsim_profile
5689 # Define the major targets
5690 all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
5692 view : $(DOC_DIR)/sawsim.pdf
5694 profile : $(BIN_DIR)/sawsim_profile | $(BIN_DIR)
5695 $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \
5696 -mnull -Mwlc -A0.39e-9,28e-9 -F8
5697 gprof $(BIN_DIR)/sawsim_profile gmon.out > $@
5698 check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim
5699 $(SHELL) -e -c 'for B in $^; do ./$$B; done'
5700 clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \
5701 clean_tension_model_utils \
5702 clean_k_model_utils clean_latex clean_check_sawsim
5703 rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \
5704 $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \
5705 $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
5706 $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
5707 $(BUILD_DIR)/global.h ./gmon.out
5708 $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
5710 # Various builds of sawsim
5711 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
5712 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5713 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
5714 gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5715 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
5716 gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5718 # Create the directories
5719 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
5722 # Copy over the source external to sawsim
5723 # Note: Cannot use, for example, `./build', because make eats the `./'
5724 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5726 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
5730 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
5734 ## Basic source generated with noweb
5735 # The central files sawsim.c and global.h...
5736 $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5738 $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5739 notangle -Rglobal.h $< > $@
5740 # ... and the modules
5741 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5742 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5743 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5744 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5745 $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5746 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5747 # Note: I use `_' as a space in my file names, but noweb doesn't like
5748 # that so we use `-' to name the noweb roots and substitute here.
5750 ## Building the unit-test programs
5752 $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
5753 notangle -Rchecks $< > $@
5754 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
5755 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
5756 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR)
5757 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
5758 clean_check_sawsim :
5759 rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
5760 # ... and the modules
5762 $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
5763 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \
5764 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \
5765 $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
5766 $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
5767 $$(BUILD_DIR)/global.h | $(BIN_DIR)
5768 gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
5769 $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
5771 # todo: clean up the required modules to
5772 $(CHECK_BINS:%=clean_%) :
5773 rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
5775 # Cleaning up the modules
5777 $(SAWSIM_MODS:%=clean_%) :
5778 rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h)
5780 <<tension model utils makefile lines>>
5781 <<k model utils makefile lines>>
5782 <<latex makefile lines>>
5784 Makefile : $(SRC_DIR)/sawsim.nw
5785 notangle -Rmakefile $< | sed 's/ /\t/' > $@
5787 This is fairly self-explanatory. We compile a dynamically linked
5788 version ([[sawsim]]) and a statically linked version
5789 ([[sawsim_static]]). The static version is about 9 times as big, but
5790 it runs without needing dynamic GSL libraries to link to, so it's a
5791 better format for distributing to a cluster for parallel evaluation.
5795 \subsection{Hookean springs in series}
5796 \label{sec.math_hooke}
5799 The effective spring constant for $n$ springs $k_i$ in series is given by
5801 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5807 F &= k_i x_i &&\forall i \le n \\
5808 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
5809 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
5810 F &= k_1 x_1 = k_\text{eff} x \\
5811 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
5812 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5817 \addcontentsline{toc}{section}{References}
5818 \bibliography{sawsim}