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{doc} % BibTeX
19 \usepackage[super,sort&compress]{natbib} % fancy citation extensions
20 % super selects citations in superscript mode
21 % sort&compress automatically sorts and compresses compound citations (\cite{a,b,...})
23 \usepackage[mathletters]{ucs} % use math fonts to allow Greek UTF-8 in the text
24 \usepackage[utf8x]{inputenc} % allow UTF-8 characters
26 %\bibliographystyle{ieeetr} % pick the bibliography style, short and sweet
27 %\bibliographystyle{plain} % pick the bibliography style, includes dates
28 \bibliographystyle{plainnat}
29 \defcitealias{sw:noweb}{{\tt noweb}}
30 \defcitealias{galassi05}{GSL}
31 \defcitealias{sw:check}{{\tt check}}
32 \defcitealias{sw:python}{Python}
37 \textheight 9.5in % leave a bit of extra room for page numbers (at top of page)
40 \setlength{\parindent}{0pt}
41 \setlength{\parskip}{5pt}
43 % For some reason, the \maketitle wants to go on it's own page
44 % so we'll just hardcode the following in our first page.
45 %\title{Sawsim: a sawtooth protein unfolding simulator}
46 %\author{W.~Trevor King}
49 \newcommand{\prog}{[[sawsim]]}
50 \newcommand{\lang}{[[C]]}
51 \newcommand{\p}[3]{\left#1 #2 \right#3} % e.g. \p({complicated expr})
52 \newcommand{\Th}{\ensuremath{\sup{\text{th}}}} % ordinal th
53 \newcommand{\U}[1]{\ensuremath{\text{ #1}}} % units: $1\U{m/s}$
55 % single spaced lists, from Alvin Alexander
56 % http://www.devdaily.com/blog/post/latex/control-line-spacing-in-itemize-enumerate-tags/
57 \newenvironment{packed_item}{
59 \setlength{\itemsep}{1pt}
60 \setlength{\parskip}{0pt}
61 \setlength{\parsep}{0pt}
65 \nwfilename{sawsim.nw}
70 Sawsim: a sawtooth protein unfolding simulator \\
75 %%%%%%%%%%%%%%%%%%%%%%%%%%%% End LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%%
77 \section{Introduction}
79 The unfolding of multi-globular protein strings is a stochastic
80 process. In the \prog\ program, we use Monte Carlo methods to
81 simulate this unfolding through an explicit time-stepping approach.
82 We develop a program in \lang\ to simulate probability of unfolding
83 under a constant extension velocity of a chain of globular domains.
85 In order to extract physical meaning from many mechanical unfolding
86 simulations, it is necessary to compare the experimental results
87 against the results of simulations using different models and
88 parameters. The model and parameters that best match the experimental
89 data provide the ‘best’ model for that experiment.
91 Previous mechanical unfolding studies have rolled their own
92 simulations for modelling their experiment, and there has been little
93 exchange and standardization of these simulations between groups.
94 The \prog\ program attempts to fill this gap, providing a flexible,
95 extensible, \emph{community} program, to make it easier for groups
96 to share the burden and increase the transparency of data analysis.
98 ‘Sawsim’ is blend of ‘sawtooth simulation’.
100 \subsection{Related work}
104 \subsection{About this document}
106 This document was written in \citetalias{sw:noweb} with code blocks in
107 \lang\ and comment blocks in \LaTeX.
109 \subsection{System Requirements and Dependencies}
112 \subsection{Change Log}
114 Version 0.0 used only the unfolded domain WLC to determine the tension
115 and had a fixed timestep $dt$.
117 Version 0.1 added adaptive timesteps, adjusting so that the unfolding
118 probability for a given domain was constant.
120 Version 0.2 added Kramers' model unfolding rate calculations, and a
121 switch to select Bell's or Kramers' model. This induced a major
122 rewrite, introducing the tension group abstraction, and a split to
123 multiple \lang\ source files. Also added a unit-test suites for
124 sanity checking, and switched to SI units throughout.
126 Version 0.3 added integrated Kramers' model unfolding rate
127 calculations. Also replaced some of my hand-coded routines with GNU
128 Scientific Library \citepalias{galassi05} calls.
130 Version 0.4 added Python evaluation for the saddle-point Kramers'
131 model unfolding rate calculations, so that the model functions could
132 be controlled from the command line. Also added interpolating lookup
133 tables to accelerate slow unfolding rate model evaluation, and command
134 line control of tension groups.
136 <<version definition>>=
137 #define VERSION "0.4"
143 sawsim - program for simulating single molecule mechanical unfolding.
144 Copyright (C) 2008, William Trevor King
146 This program is free software; you can redistribute it and/or
147 modify it under the terms of the GNU General Public License as
148 published by the Free Software Foundation; either version 3 of the
149 License, or (at your option) any later version.
151 This program is distributed in the hope that it will be useful, but
152 WITHOUT ANY WARRANTY; without even the implied warranty of
153 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
154 See the GNU General Public License for more details.
156 You should have received a copy of the GNU General Public License
157 along with this program; if not, write to the Free Software
158 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
161 The author may be contacted at <wking@drexel.edu> on the Internet, or
162 write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
163 Philadelphia PA 19104, USA.
175 The unfolding system is stored as a chain of \emph{domains}. Domains
176 can be folded, globular protein domains, unfolded protein linkers, AFM
177 cantilevers, or other stretchable system components.
179 Each domain contributes to the total tension, which depends on the
180 chain extension. The particular model for the tension as a function
181 of extension is handled generally with function pointers. So far the
182 following models have been implemented:
184 \item Null (Appendix \ref{sec.null_tension}),
185 \item Constant (Appendix \ref{sec.const_tension}),
186 \item Hookean spring (Appendix \ref{sec.hooke}),
187 \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
188 \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
191 The tension model and parameters used for a particular domain depend
192 on whether the domain is folded or unfolded. The transition rate from
193 the folded state to the unfolded state is also handled generally with
194 function pointers, with implementations for
196 \item Null (Appendix \ref{sec.null_k}),
197 \item Bell's model (Appendix \ref{sec.bell}),
198 \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
199 \item Kramers' saddle point model (Appendix \ref{sec.kramers}).
202 The unfolding of the chain domains is modeled in two stages. First
203 the tension of the chain is calculated. Then the tension (assumed to
204 be constant over the length of the chain), is applied to each folded
205 domain, and used to calculate the probability of that subdomain
206 unfolding in the next timestep $dt$. Then the time is stepped
207 forward, and the process is repeated.
210 int main(int argc, char **argv)
212 <<tension model getopt array definition>>
213 <<k model getopt array definition>>
214 list_t *domain_list=NULL; /* variables initialized in get_options() */
215 environment_t env={0};
216 double P, dt_max, v, xstep;
217 int num_folded=0, unfolding;
218 double x=0, dt, F; /* variables used in the simulation loop */
219 get_options(argc, argv, &P, &dt_max, &v, &xstep, &env, NUM_TENSION_MODELS,
220 tension_models, NUM_K_MODELS, k_models, &domain_list, &num_folded);
222 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Position (m)\tForce (N)\n");
223 if (flags & FLAG_OUTPUT_UNFOLDING_FORCES) printf("#Unfolding Force (N)\n");
224 while (num_folded > 0) {
225 F = find_tension(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, unfolding);
227 dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, v);
229 dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, 0);
230 unfolding=random_unfoldings(domain_list, F, dt, &env, &num_folded);
231 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\n", x, F);
235 if (dt == dt_max) { /* step completed */
238 } else { /* still working on this step */
243 destroy_domain_list(domain_list);
247 @ The meat of the simulation is bundled into the three functions
248 [[find_tension]], [[determine_dt]], and [[random_unfoldings]].
249 [[find_tension]] is discussed in Section \ref{sec.find_tension},
250 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
251 [[random_unfoldings]] in Section \ref{sec.unfolding_rate}.
253 The stretched distance is extended in one of two modes depending on
254 [[xstep]]. If [[xstep == 0]], then the pulling is continuous, at
255 least within the limits of the inherent discretization of a time
256 stepped approach. At any rate, the timestep is limited by the maximum
257 allowed timestep [[dt_max]] and the maximum allowed unfolding
258 probability in a given step [[P]]. In the second mode [[xstep > 0]],
259 and the end to end distance increases in discrete steps of that
260 length. The time between steps is chosen to maintain an average
261 unfolding velocity [[v]]. This approach is less work to simulate than
262 smooth pulling and also closer to the reality of many experiments, but
263 it is practically impossible to treat analytically. With both pulling
264 styles implemented, the effects of the discretization can be easily
265 calculated, bridging the gap between theory and experiment.
267 Environmental parameters important in determining reaction rates and
268 tensions (e.g. temperature, pH) are stored in a single structure to
269 facilitate extension to more complicated models in the future.
270 <<environment definition>>=
271 typedef struct environment_struct {
279 <<environment definition>>
280 <<one dimensional function definition>>
281 <<create/destroy definitions>>
283 #endif /* GLOBAL_H */
285 Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
286 included multiple times.
288 \section{Simulation functions}
290 <<simulation functions>>=
294 <<does domain unfold>>
295 <<randomly unfold domains>>
299 \label{sec.find_tension}
301 Because the stretched system may be made up of several parts (folded
302 domains, unfolded domains, spring-like cantilever, \ldots), we will
303 assign the domains to models and groups. The models are used to
304 determine a domain's tension (Hookean spring, WLC, \ldots). Groups
305 will represent collections of domains which the model should treat as
306 a single entity. A domain may belong to different groups or models
307 depending on its state. For example, a domain might be modeled as a
308 freely-jointed chain when it is folded and as a worm-like chain when
309 it is unfolded. The domains are assumed to be commutative, so
310 ordering is ignored. The interactions between the groups are assumed
311 to be linear, but the interactions between domains of the same group
312 need not be. This allows for non-linear group models such as th
313 worm-like or freely-jointed chains. Each model has a tension handler
314 function, which gives the tension $F$ needed to stretch a given group
315 of domains a distance $x$.
317 To understand the purpose of groups, consider a chain of proteins
318 where two folded proteins $A$ and $B$ are modeled as WLCs with
319 persistence length $p_f$ and two unfolded proteins $C$ and $D$ are
320 modeled as WLCs with persistence length $p_u$. The proteins are
321 attached to a cantilever $E$ qof spring constant $κ$. The string
322 would get split into three lists:
324 \begin{tabular}{llll}
325 Model & Group & List & Tension \\
326 WLC & 0 & $[A,B]$ & $F_\text{WLC}(x, p_f, L_A+L_B)$ \\
327 WLC & 1 & $[C,D]$ & $F_\text{WLC}(x, p_u, L_C+L_D)$ \\
328 Hooke & 0 & $[E]$ & $F_\text{Hooke}(x, κ)$ \\
331 Note that group numbers only matter \emph{within} model classes, since
332 grouping domains with seperate models doesn't make sense.
334 <<find tension definitions>>=
335 #define NUM_TENSION_MODELS 5
339 <<tension handler array global declaration>>=
340 extern one_dim_fn_t *tension_handlers[];
343 <<tension handler array global>>=
344 one_dim_fn_t *tension_handlers[] = {
346 &const_tension_handler,
354 <<tension model global declarations>>=
355 <<tension handler array global declaration>>
358 <<tension handler types>>=
359 typedef struct tension_handler_data_struct {
360 /* int num_domains; */
361 /* domain_t *domains; */
365 } tension_handler_data_t;
369 After sorting the chain into separate groups $G_i$, with tension
370 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
371 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
372 \forall i,j$. Note that there may be several groups within each model
373 class. [[find_tension]] is basically a wrapper to feed proper input
374 into the [[tension_balance]] function. [[unfolding]] should be set to
375 0 if there were no unfoldings since the previous call to
378 double find_tension(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list, environment_t *env, double x, int unfolding)
380 static int num_active_groups=0;
381 static one_dim_fn_t **per_group_handlers = NULL;
382 static void **per_group_params = NULL;
383 static double last_x;
384 tension_handler_data_t *tdata;
389 fprintf(stderr, "finding tension at distance %g%s\n", x, unfolding ? " after an unfolding" : "");
390 list_print(stderr, domain_list, "domain list");
393 if (unfolding != 0 || num_active_groups == 0) { /* setup statics */
394 /* free old statics */
395 if (per_group_handlers != NULL)
396 free(per_group_handlers);
397 if (per_group_params != NULL) {
398 for (i=0; i < num_active_groups; i++) {
399 assert(per_group_params[i] != NULL);
400 free(per_group_params[i]);
402 free(per_group_params);
405 /* get new statics */
406 get_active_groups(domain_list, num_tension_handlers, tension_handlers, &num_active_groups, &per_group_handlers, &per_group_params);
408 /* fill in the group handlers and params */
409 for (i=0; i<num_active_groups; i++) {
410 tdata = (tension_handler_data_t *) per_group_params[i];
412 fprintf(stderr, "active group %d of %d", i, num_active_groups);
413 list_print(stderr, tdata->group, " parameters");
416 /* tdata->persist continues unchanged */
419 } /* else roll with the current statics */
421 F = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_params, last_x, x);
425 @ For the moment, we will restrict our group boundaries to a single
426 dimension, so $\sum_i x_i = x$, our total extension (it is this
427 restriction that causes the groups to interact linearly). We'll also
428 restrict our extensions to all be positive. With these restrictions,
429 the problem of balancing the tensions reduces to solving a system of
430 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
431 the number of active groups. In general this can be fairly
432 complicated, but without much loss of practicality we can restrict
433 ourselves to strictly monotonically increasing, non-negative tension
434 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
435 simpler. For example, it guarantees the existence of a unique, real
436 solution for finite forces. See Appendix \ref{app.tension_balance}
437 for the tension balancing implementation.
439 Each group must define a parameter structure if the tension model is
441 a function to create the parameter structure,
442 a function to destroy the parameter structure, and
443 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as it's data argument.
444 The tension-specific parameter structures are contained in the domain
445 group field. See the Section \ref{app.model_selection} for a
446 discussion on the command line framework. See the worm-like chain
447 implementation for an example (Section \ref{sec.wlc}).
449 The major limitation of our approach is that all of our tension models
450 are Markovian (memory-less), which is not really a problem for the
451 chains (see \citet{evans99} for WLC thermalization rate calculations).
453 \subsection{Unfolding rate}
454 \label{sec.unfolding_rate}
456 Each folded domain is modeled as unimolecular, first order reaction
458 \text{Folded} \xrightarrow{k} % in tex: X atop Y
461 With the rate constant $k$ defined by
463 \frac{d[\text{Folded}]}{dt} = -k [\text{Folded}].
465 So the probability of a given protein unfolding in a short time $dt$
471 <<does domain unfold>>=
472 int domain_unfolds(double F, double dt, environment_t *env, domain_t *domain)
473 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to a folded domain */
475 k = accel_k(domain->k, F, env, domain->k_params);
476 //(*domain->k)(F, env, domain->k_params);
477 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
478 return happens(k*dt); /* dice roll for prob. k*dt event */
480 @ [[happens]] is a random decision making function defined in Appendix \ref{app.utils}.
482 Only one domain can unfold in each timestep, because the timescale of
483 a domain unfolding $dt_u$ is assumed to be less than the simulation
484 timestep $dt$, so a domain will completely unfold in a single
485 timestep. We adapt our timesteps to keep the probability of a single
486 domain unfolding low, so the probability of two domains unfolding in
487 the same timestep is negligible. This approach breaks down as the
488 adaptive timestep scheme approaches $dt \sim dt_u$, but $dt_u \sim
489 1\U{$\mu$s}$ for Ig27-like proteins \citep{klimov00}, so this
490 shouldn't be a problem. To reassure yourself, you can ask the
491 simulator to print the smallest timestep that was used with TODO.
492 <<randomly unfold domains>>=
493 int random_unfoldings(list_t *domain_list, double tension,
494 double dt, environment_t *env,
495 int *num_folded_domains)
496 { /* returns 1 if there was an unfolding and 0 otherwise */
497 while (domain_list != NULL) {
498 if (D(domain_list)->state == DS_FOLDED
499 && domain_unfolds(tension, dt, env, D(domain_list))) {
500 if (flags & FLAG_OUTPUT_UNFOLDING_FORCES)
501 fprintf(stdout, "%g\n", tension);
502 D(domain_list)->state = DS_UNFOLDED;
503 (*num_folded_domains)--;
504 return 1; /* our one domain has unfolded, stop looking for others */
506 domain_list = domain_list->next;
512 \subsection{Adaptive timesteps}
513 \label{sec.adaptive_dt}
515 We'd like to pick $dt$ so the probability of unfolding in the next
516 timestep is small. If we don't adapt our timestep, we risk breaking
517 our approximation that $F(x' \in [x, x+v \cdot dt]) = F(x)$ or that
518 only one domain unfolds in a given timestep. Because $F(x)$ is
519 monotonically increasing, excessively large timesteps will lead to
520 erroneously large unfolding forces. The simulation would have been
521 accurate for sufficiently small fixed timesteps, but adaptive
522 timesteps will allow us to move through low tension regions in fewer
523 steps, leading to a more efficient simulation.
525 The actual adaptive timestep implementation is not particularly
526 interesting, since we are only required to reduce $dt$ to somewhere
527 below a set threshold, so I've removed it to Appendix
528 \ref{app.adaptive_dt}.
534 The models provide the physics of the simulation, but the simulation
535 allows interchangeable models, and we are currently focusing on the
536 simulation itself, so we remove the actual model implementation to the
539 The tension models are defined in Appendix \ref{sec.tension_models}
540 and unfolding models are defined in Appendix \ref{sec.k_models}.
543 #define k_B 1.3806503e-23 /* J/K */
547 \section{Command line interface}
549 <<get option functions>>=
551 <<index tension model>>
557 \subsection{Model selection}
558 \label{app.model_selection}
560 The main difficulty with the command line interface in \prog\ is
561 developing an intuitive method for accessing the possibly large number
562 of available models. We'll treat this generally by defining an array
563 of available models, containing their important parameters.
565 <<get options definitions>>=
566 <<model getopt definitions>>
569 <<create/destroy definitions>>=
570 typedef void *create_data_func_t(char **param_strings);
571 typedef void destroy_data_func_t(void *param_struct);
574 <<model getopt definitions>>=
575 <<tension model getopt definitions>>
576 <<k model getopt definitions>>
580 \subsubsection{Tension}
582 To access the various tension models from the command line, we define
583 a structure that contains all the neccessary information about the
585 <<tension model getopt definitions>>=
586 typedef struct tension_model_getopt_struct {
587 one_dim_fn_t *handler; /* fn implementing the model on a group */
588 char *name; /* name string identifying the model */
589 char *description; /* a brief description of the model */
590 int num_params; /* number of extra params for the model */
591 char **param_descriptions; /* descriptions of extra params */
592 char *params; /* default values for the extra params */
593 create_data_func_t *creator; /* param-string -> model-data-struct fn */
594 destroy_data_func_t *destructor; /* fn to free the model data structure */
595 } tension_model_getopt_t;
598 <<tension model getopt array definition>>=
599 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
600 <<null tension model getopt>>,
601 <<constant tension model getopt>>,
602 <<hooke tension model getopt>>,
603 <<worm-like chain tension model getopt>>,
604 <<freely-jointed chain tension model getopt>>
608 \subsubsection{Unfolding rate}
610 <<k model getopt definitions>>=
611 #define NUM_K_MODELS 5
613 typedef struct k_model_getopt_struct {
618 char **param_descriptions;
620 create_data_func_t *creator;
621 destroy_data_func_t *destructor;
625 <<k model getopt array definition>>=
626 k_model_getopt_t k_models[NUM_K_MODELS] = {
627 <<null k model getopt>>,
628 <<const k model getopt>>,
629 <<bell k model getopt>>,
630 <<kramers k model getopt>>,
631 <<kramers integ k model getopt>>
638 void help(char *prog_name, double P, double dt_max, double v, double xstep,
640 int n_tension_models, tension_model_getopt_t *tension_models,
641 int folded_tension_model, int unfolded_tension_model,
642 int n_k_models, k_model_getopt_t *k_models,
646 printf("usage: %s [options]\n", prog_name);
647 printf("Version %s\n\n", VERSION);
648 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
649 printf("Simulation options:\n");
650 printf("-P P\tTarget probability for dt (currently %g)\n", P);
651 printf("-t dt\tMaximum allowed timestep dt (currently %g)\n", dt_max);
652 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
653 printf("-x dx\tPulling step size (currently %g m)\n", xstep);
654 printf("\tWhen dx = 0, the pulling is continuous\n");
655 printf("\tWhen dx > 0, the pulling occurs in discrete steps\n");
656 printf("Environmental options:\n");
657 printf("-T T\tTemperature T (currently %g K)\n", env->T);
658 printf("-C T\tYou can also set the temperature T in Celsius\n");
659 printf("Model options:\n");
660 printf("The domains exist in either the folded or unfolded state\n");
661 printf("The following options change the default behavior in each state.\n");
662 printf("-m model[,group]\tFolded tension model (currently %s)\n",
663 tension_models[folded_tension_model].name);
664 printf("-a args\tFolded tension model argument string (currently %s)\n",
665 tension_models[folded_tension_model].params);
666 printf("-M model[,group]\tUnfolded tension model (currently %s)\n",
667 tension_models[unfolded_tension_model].name);
668 printf("-A args\tUnfolded tension model argument string (currently %s)\n",
669 tension_models[unfolded_tension_model].params);
670 printf("The following options change the unfolding rate.\n");
671 printf("-k model\tTransition rate model (currently %s)\n",
672 k_models[k_model].name);
673 printf("-K args\tTransition rate model argument string (currently %s)\n",
674 k_models[k_model].params);
675 printf("Domain creation options:\n");
676 printf("Once you've set up the appropriate default models, you need to add the domains\n");
677 printf("-F n\tAdd n folded domains with the current models\n");
678 printf("-U n\tAdd n unfolded domains with the current models\n");
679 printf("Output mode options:\n");
680 printf("There are two output modes. In standard mode, only the unfolding\n");
681 printf("events are printed. For example:\n");
682 printf(" #Unfolding Force (N)\n");
683 printf(" 123.456e-12\n");
685 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
686 printf(" #Position (m)\tForce (N)\n");
687 printf(" 0.001\t0.0023\n");
689 printf("-V\tChange output to verbose mode\n");
690 printf("-h\tPrint this help and exit\n");
692 printf("Tension models:\n");
693 for (i=0; i<n_tension_models; i++) {
694 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
695 for (j=0; j < tension_models[i].num_params ; j++)
696 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
697 printf("\t\tdefaults: %s\n", tension_models[i].params);
699 printf("Unfolding rate models:\n");
700 for (i=0; i<n_k_models; i++) {
701 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
702 for (j=0; j < k_models[i].num_params ; j++)
703 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
704 printf("\t\tdefaults: %s\n", k_models[i].params);
706 printf("Examples:\n");
707 printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded\n");
708 printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8\n", prog_name);
709 printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded\n");
710 printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K5e-5,.225e-9 -mwlc -a1e-8,4e-10 -Mwlc,1 -A0.4e-9,28.1e-9 -F8\n", prog_name);
714 \subsection{Get options}
717 void get_options(int argc, char **argv,
718 double *pP, double *pDt_max, double *pV, double *pXstep,
720 int n_tension_models, tension_model_getopt_t *tension_models,
721 int n_k_models, k_model_getopt_t *k_models,
722 list_t **pDomain_list, int *pNum_folded)
724 char *prog_name = NULL;
725 char c, options[] = "P:t:v:x:T:C:m:a:M:A:k:K:F:U:Vh";
726 int ftension_model=0, utension_model=0, k_model=0;
727 char *ftension_params=NULL, *utension_params=NULL;
728 int i, n, ftension_group=0, utension_group=0;
732 extern int optind, optopt, opterr;
737 for (i=0; i<n_tension_models; i++) {
738 fprintf(stderr, "tension model %i %s:\t model's handler %p\n",
739 i, tension_models[i].name, tension_models[i].handler);
740 assert(tension_models[i].handler == tension_handlers[i]);
745 flags = FLAG_OUTPUT_UNFOLDING_FORCES;
747 *pP = 1e-3; /* % pop per s */
748 *pDt_max = 0.001; /* s */
749 *pV = 1e-6; /* m/s */
750 *pXstep = 0.0; /* m */
751 env->T = 300.0; /* K */
753 while ((c=getopt(argc, argv, options)) != -1) {
755 case 'P': *pP = atof(optarg); break;
756 case 't': *pDt_max = atof(optarg); break;
757 case 'v': *pV = atof(optarg); break;
758 case 'x': *pXstep = atof(optarg); break;
759 case 'T': env->T = atof(optarg); break;
760 case 'C': env->T = atof(optarg)+273.15; break;
762 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
763 assert(num_strings == 1 || num_strings == 2);
764 if (num_strings == 1)
767 ftension_group = atoi(strings[1]);
768 ftension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
769 ftension_params = tension_models[ftension_model].params;
770 free_parsed_list(num_strings, strings);
773 ftension_params = optarg;
776 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
777 assert(num_strings == 1 || num_strings == 2);
778 if (num_strings == 1)
781 utension_group = atoi(strings[1]);
782 utension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
783 utension_params = tension_models[utension_model].params;
784 free_parsed_list(num_strings, strings);
787 utension_params = optarg;
790 k_model = index_k_model(n_k_models, k_models, optarg);
793 k_models[k_model].params = optarg;
798 for (i=0; i<n; i++) {
799 push(pDomain_list, generate_domain(DS_FOLDED,
800 tension_models+ftension_model,
803 tension_models+utension_model,
806 k_models+k_model), 1);
813 for (i=0; i<n; i++) {
814 push(pDomain_list, generate_domain(DS_UNFOLDED,
815 tension_models+ftension_model,
818 tension_models+utension_model,
821 k_models+k_model), 1);
825 flags = FLAG_OUTPUT_FULL_CURVE;
828 fprintf(stderr, "unrecognized option '%c'\n", optopt);
829 /* fall through to default case */
831 help(prog_name, *pP, *pDt_max, *pV, *pXstep, env, n_tension_models, tension_models, ftension_model, utension_model, n_k_models, k_models, k_model);
835 /* check the arguments */
836 assert(*pDomain_list != NULL); /* no domains! */
837 assert(*pP > 0.0); assert(*pP < 1.0);
838 assert(*pDt_max > 0.0);
840 assert(env->T > 0.0);
845 <<index tension model>>=
846 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
849 for (i=0; i<n_models; i++)
850 if (STRMATCH(models[i].name, name))
853 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
861 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
864 for (i=0; i<n_models; i++)
865 if (STRMATCH(models[i].name, name))
868 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
875 <<initialize model definition>>=
876 /* requires int num_param_args and char **param_args in the current scope
878 * INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
879 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
881 #define INIT_MODEL(role, model, param_string, param_pointer) \
883 parse_list_string(param_string, SEP, '{', '}', \
884 &num_param_args, ¶m_args); \
885 if (num_param_args != model->num_params) { \
887 "%s model %s expected %d params,\n", \
888 role, model->name, model->num_params); \
890 "not the %d params in '%s'\n", \
891 num_param_args, param_string); \
892 assert(num_param_args == model->num_params); \
894 if (model->creator) \
895 param_pointer = (*model->creator)(param_args); \
897 param_pointer = NULL; \
898 free_parsed_list(num_param_args, param_args); \
903 void *generate_domain(enum domain_state_t state,
904 tension_model_getopt_t *folded_model,
906 char *folded_param_string,
907 tension_model_getopt_t *unfolded_model,
909 char *unfolded_param_string,
910 k_model_getopt_t *k_model)
912 void *folded_params, *unfolded_params, *k_params;
913 int num_param_args; /* for INIT_MODEL() */
914 char **param_args; /* for INIT_MODEL() */
917 fprintf(stderr, "generating %s ", state==DS_FOLDED ? "folded" : "unfolded");
918 fprintf(stderr, "domain (k: \"%s\", \"%s\", f: \"%s\" (%p) group %d \"%s\", u: \"%s\" (%p) group %d \"%s\")\n",
919 k_model->name, k_model->params,
920 folded_model->name, folded_model->handler, folded_group, folded_param_string,
921 unfolded_model->name, unfolded_model->handler, unfolded_group, unfolded_param_string);
923 INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
924 INIT_MODEL("unfolded", unfolded_model, unfolded_param_string, unfolded_params);
925 INIT_MODEL("k", k_model, k_model->params, k_params);
926 return create_domain(state,
927 k_model->k, k_params, k_model->destructor,
928 folded_model->handler, folded_group, folded_params, folded_model->destructor,
929 unfolded_model->handler, unfolded_group, unfolded_params, unfolded_model->destructor
936 \addcontentsline{toc}{section}{Appendicies}
938 \section{sawsim.c details}
942 The general layout of our simulation code is:
952 We include [[math.h]], so don't forget to link to the libm with `-lm'.
954 #include <assert.h> /* assert() */
955 #include <stdlib.h> /* malloc(), free(), rand() */
956 #include <stdio.h> /* fprintf(), stdout */
957 #include <string.h> /* strlen, strtok() */
958 #include <math.h> /* exp(), M_PI, sqrt() */
959 #include <sys/types.h> /* pid_t (returned by getpid()) */
960 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
963 #include "tension_balance.h"
965 #include "tension_model.h"
971 <<version definition>>
973 <<max/min definitions>>
974 <<string comparison definition and globals>>
975 <<initialize model definition>>
976 <<get options definitions>>
977 <<domain definitions>>
986 <<create/destroy domain>>
987 <<destroy domain list>>
989 <<simulation functions>>
990 <<boilerplate functions>>
993 <<boilerplate functions>>=
995 <<get option functions>>
1001 srand(getpid()*time(NULL)); /* seed rand() */
1005 <<flag definitions>>=
1006 /* in octal b/c of prefixed '0' */
1007 #define FLAG_OUTPUT_FULL_CURVE 01
1008 #define FLAG_OUTPUT_UNFOLDING_FORCES 02
1012 static unsigned long int flags = 0;
1015 \subsection{Utilities}
1018 <<max/min definitions>>=
1019 #define MAX(a,b) ((a)>(b) ? (a) : (b))
1020 #define MIN(a,b) ((a)<(b) ? (a) : (b))
1023 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
1024 <<string comparison definition and globals>>=
1025 // Check if two strings match, return 1 if they do
1026 static char *temp_string_A;
1027 static char *temp_string_B;
1028 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
1029 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
1030 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
1031 /* +1 to also compare the '\0' */
1034 We also define a macro for our [[check]] unit testing
1035 <<check relative error>>=
1036 #define CHECK_ERR(max_err, expected, received) \
1038 fail_unless( (received-expected)/expected < max_err, \
1039 "relative error %g >= %g in %s (Expected %g, received %g)", \
1040 (received-expected)/expected, max_err, #received, \
1041 expected, received); \
1042 fail_unless(-(received-expected)/expected < max_err, \
1043 "relative error %g >= %g in %s (Expected %g, received %g)", \
1044 -(received-expected)/expected, max_err, #received, \
1045 expected, received); \
1050 int happens(double probability)
1052 assert(probability >= 0.0); assert(probability <= 1.0);
1053 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*/
1058 \subsection{Adaptive timesteps}
1059 \label{app.adaptive_dt}
1061 $F(x)$ increases with $x$, possibly exploding, as in the worm-like chain model,
1062 so basing the timestep on the the unfolding probability at the current tension
1063 is dangerous, and we need to search for a $dt$ for which
1064 $P(F(x+v*dt)) < P_\text{target}$.
1065 There are two cases to consider.
1066 In the most common, no domains have unfolded since the last step,
1067 and we expect the next step to be slightly shorter than the previous one.
1068 In the less common, domains did unfold in the last step,
1069 and we expect the next step to be considerably longer than the previous one.
1071 double search_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers,
1072 list_t *domain_list,
1073 environment_t *env, double x,
1074 double target_prob, double max_dt, double v)
1075 { /* Returns the timestep dt in seconds for the current folded domain.
1076 * Takes a list of tension handlers, the list of domains,
1077 * a pointer env to the environmental data, a starting separation x in m,
1078 * a target_prob between 0 and 1,
1079 * max_dt in s, stretching velocity v in m/s.
1081 double F, k, dtCur, dtU, dtUCur, dtL, dt;
1083 /* get upper bound using the current position */
1084 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x, 0); /* BUG. repeated calculation */
1085 //printf("Start with x = %g (F = %g)\n", x, F);
1086 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1087 //printf("x %g\tF %g\tk %g\n", x, F, k);
1088 dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */
1090 //printf("overshot max_dt\n");
1093 if (v == 0) /* discrete stepping, so force is temporatily constant */
1096 /* set a lower bound on dt too */
1099 /* The dt determined above may produce illegitimate forces or ks.
1100 * Reduce the upper bound until we have valid ks. */
1102 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1103 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1106 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1108 //printf("Try for dt = %g (F = %g)\n", dt, F);
1109 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1110 /* returned k may be -1.0 */
1111 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1112 while (k == -1.0) { /* reduce step until we hit a valid k */
1114 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1115 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1116 //printf("Try for dt = %g (F = %g)\n", dt, F);
1117 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1118 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1120 assert(dtU > 1e-14); /* timestep to valid k too small */
1121 dtUCur = target_prob / k; /* safe timestep back from x+dtU */
1123 return dt; /* dtU is safe. */
1125 /* dtU wasn't safe, lets see what would be. */
1126 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1127 dt = (dtU + dtL) / 2.0;
1128 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1129 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1130 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1131 dtCur = target_prob / k;
1132 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1133 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1135 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1136 dtU = dt; /* ... stepping out only dtCur would be safe */
1139 break; /* dtCur = dt */
1141 return MAX(dtUCur, dtL);
1145 To determine $dt$ for an array of potentially different folded domains,
1146 we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
1149 double determine_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list,
1150 environment_t *env, double x,
1151 double target_prob, double dt_max, double v)
1152 { /* Returns the timestep dt in seconds.
1153 * Takes the list of folded domains, target_prob between 0 and 1,
1154 * F in N, and T in K. */
1155 double dt=dt_max, new_dt;
1156 assert(target_prob > 0.0); assert(target_prob < 1.0);
1157 assert(dt_max > 0.0);
1159 while (domain_list != NULL) {
1160 if (D(domain_list)->state == DS_FOLDED) {
1161 new_dt = search_dt(num_tension_handlers, tension_handlers, domain_list, env, x, target_prob, dt, v);
1162 dt = MIN(dt, new_dt);
1164 domain_list = domain_list->next;
1170 \subsection{Domain data}
1172 Currently domains exist in two states, folded and unfolded, and the
1173 only allowed transitions are folded $\rightarrow$ unfolded. Of
1174 course, it wouldn't be too complicated to extent this to a multi-state
1175 system, with an array containing the domains group for each possible
1176 state, and a matrix of transition-rate-calculating functions.
1177 However, at this point such generality seems unnecessary at this
1179 <<domain definitions>>=
1180 enum domain_state_t {DS_FOLDED,
1184 typedef struct domain_struct {
1185 enum domain_state_t state;
1186 one_dim_fn_t *folded_handler;
1188 one_dim_fn_t *unfolded_handler;
1190 k_func_t *k; /* function returning unfolding rate */
1191 void *folded_params; /* pointer to folded parameters */
1192 void *unfolded_params; /* pointer to unfolded parameters */
1193 void *k_params; /* pointer to k parameters */
1194 destroy_data_func_t *destroy_folded;
1195 destroy_data_func_t *destroy_unfolded;
1196 destroy_data_func_t *destroy_k;
1199 /* get the domain data for the current list node */
1200 #define D(list) ((domain_t *)(list)->d)
1201 /* get the tension params for the current list node */
1202 #define D_TP(list) (((domain_t *)(list)->d)->state == DS_FOLDED \
1203 ? ((domain_t *)(list)->d)->folded_params \
1204 : ((domain_t *)(list)->d)->unfolded_params)
1206 [[k]] is a pointer to the function determining the unfolding rate for a given tension.
1207 [[folded_params]] is a pointer to the parameters used by the function pointed to by [[k]].
1208 [[unfolded_params]] is a pointer to the parameters used by the group-appropriate handler function when determining the tension.
1209 The [[destroy_*]] pointers point to functions for freeing the memory [[*_params]].
1210 We store them with the domain data so that [[destroy_domain]] doesn't have to know which type of domain it's cleaning up after.
1212 [[create_]] and [[destroy_domain]] are simple wrappers around [[malloc]] and [[free]].
1213 <<create/destroy domain>>=
1214 domain_t *create_domain(enum domain_state_t state,
1217 destroy_data_func_t *destroy_k,
1218 one_dim_fn_t *folded_handler,
1220 void *folded_params,
1221 destroy_data_func_t *destroy_folded,
1222 one_dim_fn_t *unfolded_handler,
1224 void *unfolded_params,
1225 destroy_data_func_t *destroy_unfolded)
1227 domain_t *ret = (domain_t *)malloc(sizeof(domain_t));
1228 assert(ret != NULL);
1229 if (state == DS_FOLDED) {
1230 assert(k != NULL); /* the pointer points somewhere valid */
1231 assert(*k != NULL); /* and there is something useful there */
1233 assert(state == DS_UNFOLDED);
1235 ret->folded_handler = folded_handler;
1236 ret->folded_group = folded_group;
1237 ret->unfolded_handler = unfolded_handler;
1238 ret->unfolded_group = unfolded_group;
1240 ret->folded_params = folded_params;
1241 ret->unfolded_params = unfolded_params;
1242 ret->k_params = k_params;
1243 ret->destroy_folded = destroy_folded;
1244 ret->destroy_unfolded = destroy_unfolded;
1245 ret->destroy_k = destroy_k;
1247 fprintf(stderr, "generated domain %p. state %d\n", ret, ret->state);
1252 void destroy_domain(domain_t *domain)
1255 //printf("domain %p & %p\n", *domain, domain);
1256 if (domain->destroy_folded)
1257 (*domain->destroy_folded)(domain->folded_params);
1258 if (domain->destroy_unfolded)
1259 (*domain->destroy_unfolded)(domain->unfolded_params);
1260 if (domain->destroy_k)
1261 (*domain->destroy_k)(domain->k_params);
1267 <<destroy domain list>>=
1268 void destroy_domain_list(list_t *domain_list)
1270 domain_list = head(domain_list);
1271 while (domain_list != NULL)
1272 destroy_domain((domain_t *) pop(&domain_list));
1276 \subsection{Domain model and group handling}
1278 <<group functions>>=
1279 <<get tension handler>>
1281 <<int comparison function>>
1282 <<find possible groups>>
1285 <<get active groups>>
1288 <<get tension handler>>=
1289 one_dim_fn_t *get_tension_handler(domain_t *domain)
1291 if (domain->state == DS_FOLDED)
1292 return domain->folded_handler;
1294 if (domain->state != DS_UNFOLDED)
1295 fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
1296 assert(domain->state == DS_UNFOLDED);
1297 return domain->unfolded_handler;
1303 int get_group(domain_t *domain)
1305 if (domain->state == DS_FOLDED)
1306 return domain->folded_group;
1308 if (domain->state != DS_UNFOLDED)
1309 fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
1310 assert(domain->state == DS_UNFOLDED);
1311 return domain->unfolded_group;
1316 We already know all possible tension classes, since they are all
1317 hardcoded into \prog. However, there may be any number of possible
1318 groups. We define a function [[find_possible_groups]] to search for
1319 possible (not neccessarily active) groups. It can search for
1320 subgroups of a single [[handler]], or by setting [[handler]] to
1321 [[NULL]], subgroups of \emph{any} handler. It returns a list of group
1323 <<find possible groups>>=
1324 list_t *find_possible_groups(list_t *list, one_dim_fn_t *handler) {
1327 while (list != NULL) {
1328 if (handler == NULL || get_tension_handler(D(list)) == handler) {
1329 push(&ret, &D(list)->folded_group, 1);
1330 push(&ret, &D(list)->unfolded_group, 1);
1336 return ret; /* no groups with this handler, no processing needed */
1338 /* sort the ret list, and remove duplicates */
1339 sort(&ret, &int_comparison);
1340 unique(&ret, &int_comparison);
1342 fprintf(stderr, "handler %p possible groups:", handler);
1344 while (list != NULL) {
1345 fprintf(stderr, " %d", *((int *)list->d));
1348 fprintf(stderr, "\n");
1354 Search a [[list]] of domains, and determine whether a particular model
1355 class and group number combination is active.
1356 <<is group active>>=
1357 int is_group_active(list_t *list, one_dim_fn_t *handler, int group)
1360 while (list != NULL) {
1361 if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group)
1369 Search a [[list]] of domains, and return all domains matching a
1370 particular model class and group number.
1372 list_t *get_group_list(list_t *list, one_dim_fn_t *handler, int group)
1376 while (list != NULL) {
1377 if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group) {
1378 push(&ret, D_TP(list), 1); /* add a pointer to the appropriate tension parameters to our new list. */
1380 fprintf(stderr, "push domain %p %s tension parameters %p onto active group %p %d list %p\n", D(list), D(list)->state == DS_FOLDED ? "folded" : "unfolded", D_TP(list), handler, group, ret);
1388 Because all the node data in lists returned by [[get_group_list]] is
1389 also in the main domain list, you shouldn't destroy the node data
1390 popped off when destroying the group lists. It will all get cleaned
1391 up when the main domain list is destroyed.
1393 Put all this together to scan through a list of domains and construct
1394 an array of all the active groups.
1395 <<get active groups>>=
1396 void get_active_groups(list_t *list,
1397 int num_tension_handlers, one_dim_fn_t **tension_handlers,
1398 int *pNum_active_groups, one_dim_fn_t ***pPer_group_handlers, void ***pActive_groups)
1400 void **active_groups=NULL;
1401 one_dim_fn_t *handler, **per_group_handlers=NULL;
1402 int i, num_active_groups, *group;
1403 list_t *possible_group_numbers, *active_group_list, *active_groups_list=NULL, *per_group_handlers_list=NULL;
1405 /* get all the active groups in a list */
1406 for (i=0; i<num_tension_handlers; i++) {
1407 handler = tension_handlers[i];
1408 if (handler == NULL) continue; /* tensionless handler */
1409 possible_group_numbers = head(find_possible_groups(list, handler));
1410 while (possible_group_numbers != NULL) {
1411 group = (int *)pop(&possible_group_numbers);
1412 if (is_group_active(list, handler, *group) == 1) {
1413 active_group_list = get_group_list(list, handler, *group);
1414 push(&active_groups_list, active_group_list, 1);
1415 push(&per_group_handlers_list, handler, 1);
1417 fprintf(stderr, "active group %p for %p %d headed by tension parameters %p\n", active_group_list, handler, *group, head(active_group_list)->d);
1418 list_print(stderr, active_group_list, "active group");
1424 /* allocate the array we'll be returning */
1425 num_active_groups = list_length(active_groups_list);
1426 active_groups = (void **) malloc(sizeof(void *)*num_active_groups);
1427 assert(active_groups != NULL);
1428 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1429 assert(per_group_handlers != NULL);
1431 /* populate the array we'll be returning */
1432 active_groups_list = head(active_groups_list);
1433 for (i=0; i<num_active_groups; i++) {
1434 handler = pop(&per_group_handlers_list);
1435 per_group_handlers[i] = handler;
1437 active_groups[i] = (void *) malloc(sizeof(tension_handler_data_t));
1438 assert(active_groups[i] != NULL);
1439 active_group_list = pop(&active_groups_list);
1440 ((tension_handler_data_t *)active_groups[i])->group = active_group_list;
1441 ((tension_handler_data_t *)active_groups[i])->env = NULL;
1442 ((tension_handler_data_t *)active_groups[i])->persist = NULL;
1444 assert(active_groups_list == NULL);
1445 assert(per_group_handlers_list == NULL);
1447 *pNum_active_groups = num_active_groups;
1448 *pPer_group_handlers = per_group_handlers;
1449 *pActive_groups = active_groups;
1454 \section{String parsing}
1456 For handling command line arguments, we need parse delimited strings (e.g. [["param1,param2,param3..."]]).
1457 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1458 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1459 need to take care of parsing those parameters themselves.
1460 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1466 <<parse definitions>>
1467 <<parse declarations>>
1471 <<parse module makefile lines>>=
1472 NW_SAWSIM_MODS += parse
1473 CHECK_BINS += check_parse
1477 <<parse definitions>>=
1478 #define SEP ',' /* argument separator character */
1481 <<parse declarations>>=
1482 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1483 int *num, char ***string_array);
1484 extern void free_parsed_list(int num, char **string_array);
1487 [[parse_list_string]] allocates memory, don't forget to free it
1488 afterward with [[free_parsed_list]]. It does not alter the original.
1490 The string may start off with a [[deeper]] character
1491 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1492 [[x,y]], where the pointer is one character in on the copied string.
1493 However, when we go to free the memory, we need a pointer to the
1494 beginning of the string. In order to accommodate this for a string
1495 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1496 the first $N$ elements point to the separated arguments, and let the
1497 last element point to the start of the copied string regardless of
1499 <<parse delimited list string functions>>=
1500 /* TODO, split out into parse.hc */
1501 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1504 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1505 if (string[i] == deeper) {depth++;}
1506 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1512 void parse_list_string(char *string, char sep, char deeper, char shallower,
1513 int *num, char ***string_array)
1515 char *str=NULL, **ret=NULL;
1517 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
1519 *string_array = NULL;
1522 /* make a copy of the string, so we don't change the original */
1523 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1524 assert(str != NULL);
1525 strcpy(str, string); /* we know str is long enough */
1526 /* count the number of regions, so we can allocate pointers to them */
1529 n++; i++; /* move on to next argument */
1530 i += next_delim_index(str+i, sep, deeper, shallower);
1531 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1532 } while (str[i] != '\0');
1533 ret = (char **)malloc(sizeof(char *)*(n+1));
1534 assert(ret != NULL);
1535 /* replace the separators with '\0' & assign pointers */
1536 ret[n] = str; /* point to the front of the copied string */
1539 for(i=1; i<n; i++) {
1540 j += next_delim_index(str+j, sep, deeper, shallower);
1542 ret[i] = str+j; /* point to the separated arguments */
1544 /* strip off bounding braces, if any */
1545 for(i=0; i<n; i++) {
1546 if (ret[i][0] == deeper) {
1550 if (ret[i][strlen(ret[i])-1] == shallower)
1551 ret[i][strlen(ret[i])-1] = '\0';
1554 *string_array = ret;
1557 void free_parsed_list(int num, char **string_array)
1560 /* frees all items, since they were allocated as a single string*/
1561 free(string_array[num]);
1562 /* and free the array of pointers */
1568 \subsection{Parsing implementation}
1574 <<parse delimited list string functions>>
1578 #include <assert.h> /* assert() */
1579 #include <stdlib.h> /* NULL */
1580 #include <stdio.h> /* fprintf(), stdout *//*!!*/
1581 #include <string.h> /* strlen() */
1585 \subsection{Parsing unit tests}
1587 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1589 <<parse check includes>>
1590 <<string comparison definition and globals>>
1591 <<check relative error>>
1592 <<parse test suite>>
1593 <<main check program>>
1596 <<parse check includes>>=
1597 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
1598 #include <stdio.h> /* printf() */
1599 #include <assert.h> /* assert() */
1600 #include <string.h> /* strlen() */
1605 <<parse test suite>>=
1606 <<parse list string tests>>
1607 <<parse suite definition>>
1610 <<parse suite definition>>=
1611 Suite *test_suite (void)
1613 Suite *s = suite_create ("k model");
1614 <<parse list string test case defs>>
1616 <<parse list string test case add>>
1621 <<parse list string tests>>=
1624 START_TEST(test_next_delim_index)
1626 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
1627 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
1628 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
1629 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
1630 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
1634 START_TEST(test_parse_list_null)
1638 parse_list_string(NULL, SEP, '{', '}',
1639 &num_param_args, ¶m_args);
1640 fail_unless(num_param_args == 0, NULL);
1641 fail_unless(param_args == NULL, NULL);
1644 START_TEST(test_parse_list_single_simple)
1648 parse_list_string("arg", SEP, '{', '}',
1649 &num_param_args, ¶m_args);
1650 fail_unless(num_param_args == 1, NULL);
1651 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
1654 START_TEST(test_parse_list_single_compound)
1658 parse_list_string("{x,y,z}", SEP, '{', '}',
1659 &num_param_args, ¶m_args);
1660 fail_unless(num_param_args == 1, NULL);
1661 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
1664 START_TEST(test_parse_list_double_simple)
1668 parse_list_string("abc,def", SEP, '{', '}',
1669 &num_param_args, ¶m_args);
1670 fail_unless(num_param_args == 2, NULL);
1671 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
1672 fail_unless(STRMATCH(param_args[1],"def"), NULL);
1676 <<parse list string test case defs>>=
1677 TCase *tc_parse_list_string = tcase_create("parse list string");
1679 <<parse list string test case add>>=
1680 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
1681 tcase_add_test(tc_parse_list_string, test_parse_list_null);
1682 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
1683 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
1684 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
1685 suite_add_tcase(s, tc_parse_list_string);
1689 \section{Unit tests}
1691 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1698 <<check relative error>>
1701 <<main check program>>
1713 <<determine dt tests>>
1715 <<does domain unfold tests>>
1716 <<randomly unfold domains tests>>
1717 <<suite definition>>
1720 <<suite definition>>=
1721 Suite *test_suite (void)
1723 Suite *s = suite_create ("sawsim");
1724 <<F test case defs>>
1725 <<determine dt test case defs>>
1726 <<happens test case defs>>
1727 <<does domain unfold test case defs>>
1728 <<randomly unfold domains test case defs>>
1731 <<determine dt test case add>>
1732 <<happens test case add>>
1733 <<does domain unfold test case add>>
1734 <<randomly unfold domains test case add>>
1737 tcase_add_checked_fixture(tc_strip_address,
1738 setup_strip_address,
1739 teardown_strip_address);
1745 <<main check program>>=
1749 Suite *s = test_suite();
1750 SRunner *sr = srunner_create(s);
1751 srunner_run_all(sr, CK_ENV);
1752 number_failed = srunner_ntests_failed(sr);
1754 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
1758 \subsection{F tests}
1760 <<worm-like chain tests>>
1762 <<F test case defs>>=
1763 <<worm-like chain test case def>>
1765 <<F test case add>>=
1766 <<worm-like chain test case add>>
1769 <<worm-like chain tests>>=
1770 extern double wlc(double x, double T, double p, double L);
1771 START_TEST(test_wlc_at_zero)
1773 double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
1774 fail_unless(abs(wlc(x, T, p, L)) < lim, \
1775 "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
1776 x, T, p, L, abs(wlc(x, T, p, L)), lim);
1780 START_TEST(test_wlc_at_half)
1782 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
1783 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
1784 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
1786 * linear term = x/L = 0.5
1787 * nonlinear + linear = 0.75 + 0.5 = 1.25
1788 * wlc = 10e21*1.25 = 12.5e21
1790 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
1791 "wlc(%g, %g, %g, %g) = %g != %g",
1792 x, T, p, L, wlc(x, T, p, L), 12.5e21);
1796 <<worm-like chain test case def>>=
1797 TCase *tc_wlc = tcase_create("WLC");
1800 <<worm-like chain test case add>>=
1801 tcase_add_test(tc_wlc, test_wlc_at_zero);
1802 tcase_add_test(tc_wlc, test_wlc_at_half);
1803 suite_add_tcase(s, tc_wlc);
1806 \subsection{Model tests}
1808 Check the searching with [[linear_k]].
1809 Check overwhelming force treatment with the heavyside-esque step [[?]].
1810 <<determine dt tests>>=
1811 double linear_k(double F, environment_t *env, void *params)
1813 double Fnot = *(double *)params;
1817 START_TEST(test_determine_dt_linear_k)
1820 double dt_max=3.0, Fnot=3.0;
1821 double F[]={0,1,2,3,4,5,6};
1822 domain_t dom; /* use both parts at once for folded/unfolded */
1826 dom->next = dom->prev = NULL;
1827 dom->k_func_t = linear_k;
1828 dom->folded_params = &Fnot;
1829 dom->unfolded_params = !!!!!!!!!
1830 dom->destroy_folded = dom->destroy_unfolded = NULL;
1831 for( i=0; i < sizeof(F)/sizeof(double); i++) {
1832 fail_unless(determine_dt(folded, unfolded, targ_p, dt_max, v, x, T)==1, NULL);
1838 <<determine dt test case defs>>=
1839 TCase *tc_determine_dt = tcase_create("Determine dt");
1841 <<determine dt test case add>>=
1842 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
1843 suite_add_tcase(s, tc_determine_dt);
1848 <<happens test case defs>>=
1850 <<happens test case add>>=
1853 <<does domain unfold tests>>=
1855 <<does domain unfold test case defs>>=
1857 <<does domain unfold test case add>>=
1860 <<randomly unfold domains tests>>=
1862 <<randomly unfold domains test case defs>>=
1864 <<randomly unfold domains test case add>>=
1868 \section{Balancing group extensions}
1869 \label{app.tension_balance}
1871 For a given total extension $x$ (of the piezo), the various domain
1872 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
1873 amounts, and we need to tweak the portion that each extends to balance
1874 the tension amongst the active groups. Since the tension balancing is
1875 separable from the bulk of the simulation, we break it out into a
1876 separate module. The interface is defined in [[tension_balance.h]],
1877 the implementation in [[tension_balance.c]], and the unit testing in
1878 [[check_tension_balance.c]]
1880 <<tension-balance.h>>=
1881 #ifndef TENSION_BALANCE_H
1882 #define TENSION_BALANCE_H
1884 <<tension balance function declaration>>
1885 #endif /* TENSION_BALANCE_H */
1888 <<tension balance functions>>=
1889 <<one dimensional bracket>>
1890 <<one dimensional solve>>
1892 <<tension balance function>>
1895 <<tension balance module makefile lines>>=
1896 NW_SAWSIM_MODS += tension_balance
1897 CHECK_BINS += check_tension_balance
1898 check_tension_balance_MODS =
1901 The entire force balancing problem reduces to a solving two nested
1902 one-dimensional root-finding problems. First we define the one
1903 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
1904 non-empty group). $x(x_0)$ is also strictly monotonically increasing,
1905 so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
1906 stores the last successful value of $x$, since we expect to be taking
1907 small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
1908 indicates that the groups have changed.
1909 <<tension balance function declaration>>=
1910 double tension_balance(int num_tension_groups,
1911 one_dim_fn_t **tension_handlers,
1916 <<tension balance function>>=
1917 double tension_balance(int num_tension_groups,
1918 one_dim_fn_t **tension_handlers,
1923 static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL};
1924 double F, xo_guess, xo, lb, ub;
1925 double min_relx=1e-6, min_rely=1e-6;
1926 int max_steps=200, i;
1928 assert(num_tension_groups > 0);
1929 assert(tension_handlers != NULL);
1930 assert(params != NULL);
1931 assert(num_tension_groups > 0);
1933 if (num_tension_groups == 1) { /* only one group, no balancing required */
1936 //fprintf(stderr, "balancing for x = %g with ", x);
1937 //for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
1938 //fprintf(stderr, "\n");
1939 if (last_x == -1) { /* new group setup, reset x_xo_data */
1940 /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
1941 if (x_xo_data.xi != NULL)
1943 /* construct new x_xo_data */
1944 x_xo_data.n_groups = num_tension_groups;
1945 x_xo_data.tension_handler = tension_handlers;
1946 x_xo_data.handler_data = params;
1947 x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
1948 for (i=0; i<num_tension_groups; i++)
1949 x_xo_data.xi[i] = -1.0;
1951 if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
1952 if (x == 0) xo_guess = length_scale;
1953 else xo_guess = x/num_tension_groups;
1955 fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
1957 oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
1958 } else { /* work off the last balanced point */
1960 fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
1963 lb = x_xo_data.xi[0];
1964 ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */
1965 } else if (x < last_x) {
1966 lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
1967 ub = x_xo_data.xi[0];
1968 } else { /* x == last_x */
1969 fprintf(stderr, "not moving\n");
1970 lb= x_xo_data.xi[0];
1971 ub= x_xo_data.xi[0];
1975 fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
1976 __FUNCTION__, x, lb, ub);
1978 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
1979 F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
1981 if (num_tension_groups > 1) {
1982 fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
1983 for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
1984 fprintf(stderr, "\n");
1988 F = (*tension_handlers[0])(xo, params[0]);
1993 <<tension balance internal definitions>>=
1994 <<x of x0 definitions>>
1997 <<x of x0 definitions>>=
1998 typedef struct x_of_xo_data_struct {
2000 one_dim_fn_t **tension_handler; /* array of fn pointers */
2001 void **handler_data; /* array of void* pointers */
2007 double x_of_xo(double xo, void *pdata)
2009 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
2010 double F, x=0, xi, xi_guess, lb, ub;
2012 double min_relx=1e-6, min_rely=1e-6;
2014 assert(data->n_groups > 0);
2016 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
2018 fprintf(stderr, "%s: finding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
2021 if (data->xi) data->xi[0] = xo;
2022 for (i=1; i < data->n_groups; i++) {
2023 if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
2024 else xi_guess = data->xi[i];
2026 fprintf(stderr, "%s: searching for proper x[%d] for F[%d](x[%d]) = %g N\n", __FUNCTION__, i, i, i, F);
2028 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
2029 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
2031 fprintf(stderr, "%s: found F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
2035 if (data->xi) data->xi[i] = xi;
2038 fprintf(stderr, "%s: found x(x0=%g) = %g\n", __FUNCTION__, xo, x);
2044 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
2045 number of steps for monotonically increasing $f(x)$. Simple
2046 bisection, so it's robust and converges fairly quickly. You can set
2047 the maximum number of steps to take, but if you have enough processor
2048 power, it's better to limit your precision with the relative limits
2049 [[min_relx]] and [[y]]. These ensure that the width of the bracket is
2050 small on the length scale of it's larger side. Note that \emph{both}
2051 relative limits must be satisfied for the search to stop.
2052 <<one dimensional function definition>>=
2053 /* equivalent to gsl_func_t */
2054 typedef double one_dim_fn_t(double x, void *params);
2056 <<one dimensional solve declaration>>=
2057 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2058 double min_relx, double min_rely, int max_steps,
2062 <<one dimensional solve>>=
2063 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2064 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2065 double min_relx, double min_rely, int max_steps,
2068 double yx, ylb, yub, x;
2071 ylb = (*f)(lb, params);
2072 yub = (*f)(ub, params);
2074 /* check some simple cases */
2076 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bounded */
2077 else return lb; /* any x on the range [lb, ub] would work */
2079 if (ylb == y) { x=lb; yx=ylb; goto end; }
2080 if (yub == y) { x=ub; yx=yub; goto end; }
2083 fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2085 assert(ylb < y); assert(yub > y);
2086 for (i=0; i<max_steps; i++) {
2088 yx = (*f)(x, params);
2090 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);
2092 if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2093 else if (yx > y) { ub=x; yub=yx; }
2094 else /* yx < y */ { lb=x; ylb=yx; }
2099 fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2105 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2106 Generate bracketing $x$ values through bisection or exponential growth.
2107 <<one dimensional bracket declaration>>=
2108 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2111 <<one dimensional bracket>>=
2112 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2114 double yx, step, x=xguess;
2116 yx = (*f)(x, params);
2118 fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2125 fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2129 if (x == 0) x = length_scale/2.0; /* HACK */
2132 yx = (*f)(x, params);
2134 fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2139 fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2142 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2146 \subsection{Balancing implementation}
2148 <<tension-balance.c>>=
2151 <<tension balance includes>>
2152 #include "tension_balance.h"
2154 double length_scale = 1e-10; /* HACK */
2156 <<tension balance internal definitions>>
2157 <<tension balance functions>>
2160 <<tension balance includes>>=
2161 #include <assert.h> /* assert() */
2162 #include <stdlib.h> /* NULL */
2163 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2164 #include <stdio.h> /* fprintf(), stdout */
2168 \subsection{Balancing unit tests}
2170 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2171 <<check-tension-balance.c>>=
2172 <<tension balance check includes>>
2173 <<tension balance test suite>>
2174 <<main check program>>
2177 <<tension balance check includes>>=
2178 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2179 #include <assert.h> /* assert() */
2182 #include "tension_balance.h"
2185 <<tension balance test suite>>=
2186 <<tension balance function tests>>
2187 <<tension balance suite definition>>
2190 <<tension balance suite definition>>=
2191 Suite *test_suite (void)
2193 Suite *s = suite_create ("tension balance");
2194 <<tension balance function test case defs>>
2196 <<tension balance function test case adds>>
2201 <<tension balance function tests>>=
2202 <<check relative error>>
2204 double hooke(double x, void *pK)
2207 return *((double*)pK) * x;
2210 START_TEST(test_single_function)
2212 double k=5, x=3, last_x=2.0, F;
2213 one_dim_fn_t *handlers[] = {&hooke};
2214 void *data[] = {&k};
2215 F = tension_balance(1, handlers, data, last_x, x);
2216 fail_unless(F = k*x, NULL);
2221 We can also test balancing two springs with different spring constants.
2222 The analytic solution for a general number of springs is given in Appendix \ref{app.math_hooke}.
2223 <<tension balance function tests>>=
2224 START_TEST(test_double_hooke)
2226 double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
2227 one_dim_fn_t *handlers[] = {&hooke, &hooke};
2228 void *data[] = {&k1, &k2};
2229 F = tension_balance(2, handlers, data, last_x, x);
2233 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2234 //CHECK_ERR(1e-6, x1e, xi[0]);
2235 //CHECK_ERR(1e-6, x2e, xi[1]);
2236 CHECK_ERR(1e-6, Fe, F);
2240 TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
2242 <<tension balance function test case defs>>=
2243 TCase *tc_tbfunc = tcase_create("tension balance function");
2246 <<tension balance function test case adds>>=
2247 tcase_add_test(tc_tbfunc, test_single_function);
2248 tcase_add_test(tc_tbfunc, test_double_hooke);
2249 suite_add_tcase(s, tc_tbfunc);
2254 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2255 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2256 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2262 <<list definitions>>
2263 <<list declarations>>
2267 <<list declarations>>=
2268 <<head and tail declarations>>
2269 <<list length declaration>>
2270 <<push declaration>>
2272 <<list destroy declaration>>
2273 <<list to array declaration>>
2274 <<move declaration>>
2275 <<sort declaration>>
2276 <<unique declaration>>
2277 <<list print declaration>>
2281 <<create and destroy node>>
2294 <<list module makefile lines>>=
2295 NW_SAWSIM_MODS += list
2296 CHECK_BINS += check_list
2300 <<list definitions>>=
2301 typedef struct list_struct {
2302 struct list_struct *next;
2303 struct list_struct *prev;
2307 <<comparison function typedef>>
2310 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2311 <<head and tail declarations>>=
2312 list_t *head(list_t *list);
2313 list_t *tail(list_t *list);
2316 list_t *head(list_t *list)
2320 while (list->prev != NULL) {
2326 list_t *tail(list_t *list)
2330 while (list->next != NULL) {
2337 <<list length declaration>>=
2338 int list_length(list_t *list);
2341 int list_length(list_t *list)
2348 while (list->next != NULL) {
2356 [[push]] inserts a node relative to the current position in the list
2357 without changing the current position.
2358 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2359 so we set it to that of the pushed domain.
2360 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2361 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2362 <<push declaration>>=
2363 void push(list_t **pList, void *data, int below);
2366 void push(list_t **pList, void *data, int below)
2368 list_t *list, *new_node;
2369 assert(pList != NULL);
2371 new_node = create_node(data);
2374 else if (below==0) { /* insert above */
2375 if (list->prev != NULL)
2376 list->prev->next = new_node;
2377 new_node->prev = list->prev;
2378 list->prev = new_node;
2379 new_node->next = list;
2380 } else { /* insert below */
2381 if (list->next != NULL)
2382 list->next->prev = new_node;
2383 new_node->next = list->next;
2384 list->next = new_node;
2385 new_node->prev = list;
2390 [[pop]] removes the current domain node, moving the current position
2391 to the node after it, unless that node is [[NULL]], in which case move
2392 the current position to the node before it.
2393 <<pop declaration>>=
2394 void *pop(list_t **pList);
2397 void *pop(list_t **pList)
2399 list_t *list, *popped;
2401 assert(pList != NULL);
2403 assert(list != NULL); /* not an empty list */
2405 /* bypass the popped node */
2406 if (list->prev != NULL)
2407 list->prev->next = list->next;
2408 if (list->next != NULL) {
2409 list->next->prev = list->prev;
2410 *pList = list->next;
2412 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
2414 destroy_node(popped);
2419 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
2420 If the cleanup function is [[NULL]], no cleanup function is called.
2421 The nodes are not popped in any particular order.
2422 <<list destroy declaration>>=
2423 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
2426 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
2430 assert(pList != NULL);
2433 assert(list != NULL); /* not an empty list */
2434 while (list != NULL) {
2436 if (destroy != NULL)
2442 [[list_to_array]] converts a list to an array of pointers, preserving order.
2443 <<list to array declaration>>=
2444 void list_to_array(list_t *list, int *length, void ***array);
2447 void list_to_array(list_t *list, int *pLength, void ***pArray)
2451 assert(list != NULL);
2452 assert(pLength != NULL);
2453 assert(pArray != NULL);
2455 length = list_length(list);
2456 array = (void **)malloc(sizeof(void *)*length);
2457 assert(array != NULL);
2460 while (list != NULL) {
2461 array[i++] = list->d;
2469 [[move]] moves the current element along the list in either direction.
2470 <<move declaration>>=
2471 void move(list_t **pList, int downstream);
2474 void move(list_t **pList, int downstream)
2476 list_t *list, *flipper;
2478 assert(pList != NULL);
2479 list = *pList; /* a>B>c>d */
2480 assert(list != NULL); /* not an empty list */
2481 if (downstream == 0)
2482 flipper = list->prev; /* flipper is A */
2484 flipper = list->next; /* flipper is C */
2485 /* ensure there is somewhere to go in the direction we're heading */
2486 assert(flipper != NULL);
2487 /* remove the current element */
2488 data = pop(&list); /* data=B, a>C>d */
2489 /* and put it back in in the correct spot */
2491 if (downstream == 0) {
2492 push(&list, data, 0); /* b>A>c>d */
2493 list = list->prev; /* B>a>c>d */
2495 push(&list, data, 1); /* a>C>b>d */
2496 list = list->next; /* a>c>B>d */
2502 [[sort]] sorts a list from smallest to largest according to a user
2503 specified comparison.
2504 <<comparison function typedef>>=
2505 typedef int (comparison_fn_t)(void *A, void *B);
2508 <<int comparison function>>=
2509 int int_comparison(void *A, void *B)
2514 if (a > b) return 1;
2515 else if (a == b) return 0;
2519 <<sort declaration>>=
2520 void sort(list_t **pList, comparison_fn_t *comp);
2522 Here we use the bubble sort algorithm.
2524 void sort(list_t **pList, comparison_fn_t *comp)
2528 assert(pList != NULL);
2530 assert(list != NULL); /* not an empty list */
2531 while (stable == 0) {
2534 while (list->next != NULL) {
2535 if ((*comp)(list->d, list->next->d) > 0) {
2537 move(&list, 1 /* downstream */);
2546 [[unique]] removes duplicates from a sorted list according to a user
2547 specified comparison. Don't do this unless you have other pointers
2548 any dynamically allocated data in your list, because [[unique]] just
2549 drops any non-unique data without freeing it.
2550 <<unique declaration>>=
2551 void unique(list_t **pList, comparison_fn_t *comp);
2554 void unique(list_t **pList, comparison_fn_t *comp)
2557 assert(pList != NULL);
2559 assert(list != NULL); /* not an empty list */
2561 while (list->next != NULL) {
2562 if ((*comp)(list->d, list->next->d) == 0) {
2571 [[list_print]] prints a list to a [[FILE *]] stream.
2572 <<list print declaration>>=
2573 void list_print(FILE *file, list_t *list, char *list_name);
2576 void list_print(FILE *file, list_t *list, char *list_name)
2578 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
2580 while (list != NULL) {
2581 fprintf(file, " %p", list->d);
2584 fprintf(file, "\n");
2588 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
2589 <<create and destroy node>>=
2590 list_t *create_node(void *data)
2592 list_t *ret = (list_t *)malloc(sizeof(list_t));
2593 assert(ret != NULL);
2600 void destroy_node(list_t *node)
2606 The user must free the data pointed to by the node on their own.
2608 \subsection{List implementation}
2618 #include <assert.h> /* assert() */
2619 #include <stdlib.h> /* malloc(), free(), rand() */
2620 #include <stdio.h> /* fprintf(), stdout, FILE */
2621 #include "global.h" /* destroy_data_func_t */
2624 \subsection{List unit tests}
2626 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2628 <<list check includes>>
2630 <<main check program>>
2633 <<list check includes>>=
2634 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2635 #include <stdio.h> /* FILE */
2641 <<list test suite>>=
2644 <<list suite definition>>
2647 <<list suite definition>>=
2648 Suite *test_suite (void)
2650 Suite *s = suite_create ("list");
2651 <<push test case defs>>
2652 <<pop test case defs>>
2654 <<push test case adds>>
2655 <<pop test case adds>>
2661 START_TEST(test_push)
2666 push(&list, (void *)i, 1);
2667 fail_unless(list == head(list), NULL);
2668 fail_unless( (int)list->d == 0 );
2669 for (i=0; i<n; i++) {
2670 /* we expect 0, n-1, n-2, ..., 1 */
2673 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
2679 <<push test case defs>>=
2680 TCase *tc_push = tcase_create("push");
2683 <<push test case adds>>=
2684 tcase_add_test(tc_push, test_push);
2685 suite_add_tcase(s, tc_push);
2690 <<pop test case defs>>=
2692 <<pop test case adds>>=
2695 \section{Function string evaluation}
2697 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).
2698 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
2699 The solution is to run a scripting language as a subprocess accessed via pipes.
2700 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
2702 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
2703 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
2704 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.
2705 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
2707 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.
2708 Then you can either statically or dynamically link to those hardcoded functions.
2709 While much less flexible, this approach would be a fairly simple-to-implement fallback.
2711 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
2712 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
2715 #ifndef STRING_EVAL_H
2716 #define STRING_EVAL_H
2718 <<string eval setup declaration>>
2719 <<string eval function declaration>>
2720 <<string eval teardown declaration>>
2721 #endif /* STRING_EVAL_H */
2724 <<string eval module makefile lines>>=
2725 NW_SAWSIM_MODS += string_eval
2726 CHECK_BINS += check_string_eval
2727 check_string_eval_MODS =
2730 For an introduction to POSIX process control, see\\
2731 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
2732 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
2733 and of course, the relavent [[man]] pages.
2735 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
2736 [[execvp]] replaces the calling process' program with a new program.
2737 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
2738 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
2739 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
2740 The new program needs command line arguments, just like it would if you were running it from a shell.
2741 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
2742 with the final array entry being a [[NULL]] pointer.
2744 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
2745 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
2746 <<string eval subprocess definitions>>=
2747 #define SUBPROCESS "python"
2748 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
2749 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
2750 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
2752 The [[i]] option lets Python know that it should run in interactive mode.
2753 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
2754 In interactive mode, python acts on every instruction as soon as it is recieved.
2755 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.
2756 %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.
2758 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
2759 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
2760 [[fork]] returns two copies of the same program, executing the original code.
2761 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
2762 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
2764 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.
2765 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
2766 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
2767 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
2768 <<string eval pipe definitions>>=
2769 #define PIPE_READ 0 /* the end you read from */
2770 #define PIPE_WRITE 1 /* the end you write to */
2772 #define STDIN 0 /* initial index of stdin pair */
2773 #define STDOUT 2 /* initial index of stdout pair */
2776 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
2778 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.
2780 <<string eval setup declaration>>=
2781 extern void string_eval_setup(FILE **pIn, FILE **pOut);
2783 <<string eval setup definition>>=
2784 void string_eval_setup(FILE **pIn, FILE **pOut)
2787 int pfd[4]; /* file descriptors for stdin and stdout */
2789 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
2790 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
2792 pid = fork(); /* split process into two copies */
2794 if (pid == -1) { /* fork error */
2795 perror("fork error");
2797 } else if (pid == 0) { /* child */
2798 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
2799 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
2800 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
2801 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
2802 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
2803 perror("exec error"); /* exec shouldn't return */
2805 } else { /* parent */
2806 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
2807 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
2808 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
2809 if ( *pIn == NULL ) {
2810 perror("fdopen (in)");
2813 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
2814 if ( *pOut == NULL ) {
2815 perror("fdopen (out)");
2822 To use the evaluating subprocess, we just pipe in our command, and read out the results.
2823 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
2824 <<string eval function declaration>>=
2825 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
2827 <<string eval function definition>>=
2828 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
2831 rc = fprintf(in, "%s", input);
2832 assert(rc == strlen(input));
2835 alarm(1); /* set a one second timeout on the read */
2836 assert( fgets(output, buflen, out) != NULL );
2837 alarm(0); /* cancel the timeout */
2838 //fprintf(stderr, "eval: %s ----> %s", input, output);
2841 The [[alarm]] calls set and clear a timeout on the returned output.
2842 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.
2843 This protects against invalid input for which a line of output is not printed to [[stdout]].
2844 Other invalid inputs (e.g. those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
2845 If you are getting strange results, check your python code seperately. TODO, better error handling.
2847 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
2848 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
2849 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
2850 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
2851 <<string eval teardown declaration>>=
2852 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
2855 <<string eval teardown definition>>=
2856 void string_eval_teardown(FILE **pIn, FILE **pOut)
2861 /* redirect Python's stderr */
2862 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
2866 assert( fclose(*pIn) == 0 );
2868 assert( fclose(*pOut) == 0 );
2871 /* wait for python to exit */
2873 pid = wait(&stat_loc);
2880 if (WIFEXITED(stat_loc)) {
2881 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
2882 } else if (WIFSIGNALED(stat_loc)) {
2883 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
2888 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
2890 \subsection{String evaluation implementation}
2894 <<string eval includes>>
2895 #include "string_eval.h"
2896 <<string eval internal definitions>>
2897 <<string eval functions>>
2900 <<string eval includes>>=
2901 #include <assert.h> /* assert() */
2902 #include <stdlib.h> /* NULL */
2903 #include <stdio.h> /* fprintf(), stdout, fdopen() */
2904 #include <string.h> /* strlen() */
2905 #include <sys/types.h> /* pid_t */
2906 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
2907 #include <wait.h> /* wait() */
2910 <<string eval internal definitions>>=
2911 <<string eval subprocess definitions>>
2912 <<string eval pipe definitions>>
2915 <<string eval functions>>=
2916 <<string eval setup definition>>
2917 <<string eval function definition>>
2918 <<string eval teardown definition>>
2921 \subsection{String evaluation unit tests}
2923 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2924 <<check-string-eval.c>>=
2925 <<string eval check includes>>
2926 <<string comparison definition and globals>>
2927 <<string eval test suite>>
2928 <<main check program>>
2931 <<string eval check includes>>=
2932 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
2933 #include <stdio.h> /* printf() */
2934 #include <assert.h> /* assert() */
2935 #include <string.h> /* strlen() */
2936 #include <signal.h> /* SIGKILL */
2938 #include "string_eval.h"
2941 <<string eval test suite>>=
2942 <<string eval tests>>
2943 <<string eval suite definition>>
2946 <<string eval suite definition>>=
2947 Suite *test_suite (void)
2949 Suite *s = suite_create ("string eval");
2950 <<string eval test case defs>>
2952 <<string eval test case add>>
2957 <<string eval tests>>=
2958 START_TEST(test_setup_teardown)
2961 string_eval_setup(&in, &out);
2962 string_eval_teardown(&in, &out);
2965 START_TEST(test_invalid_command)
2968 char input[80], output[80]={};
2969 string_eval_setup(&in, &out);
2970 sprintf(input, "print ABCDefg\n");
2971 string_eval(in, out, input, 80, output);
2972 string_eval_teardown(&in, &out);
2975 START_TEST(test_simple_eval)
2978 char input[80], output[80]={};
2979 string_eval_setup(&in, &out);
2980 sprintf(input, "print 3+4*5\n");
2981 string_eval(in, out, input, 80, output);
2982 fail_unless(STRMATCH(output,"23\n"), NULL);
2983 string_eval_teardown(&in, &out);
2986 START_TEST(test_multiple_evals)
2989 char input[80], output[80]={};
2990 string_eval_setup(&in, &out);
2991 sprintf(input, "print 3+4*5\n");
2992 string_eval(in, out, input, 80, output);
2993 fail_unless(STRMATCH(output,"23\n"), NULL);
2994 sprintf(input, "print (3**2 + 4**2)**0.5\n");
2995 string_eval(in, out, input, 80, output);
2996 fail_unless(STRMATCH(output,"5.0\n"), NULL);
2997 string_eval_teardown(&in, &out);
3000 START_TEST(test_eval_with_state)
3003 char input[80], output[80]={};
3004 string_eval_setup(&in, &out);
3005 sprintf(input, "print 3+4*5\n");
3006 fprintf(in, "A = 3\n");
3007 sprintf(input, "print A*3\n");
3008 string_eval(in, out, input, 80, output);
3009 fail_unless(STRMATCH(output,"9\n"), NULL);
3010 string_eval_teardown(&in, &out);
3014 <<string eval test case defs>>=
3015 TCase *tc_string_eval = tcase_create("string_eval");
3017 <<string eval test case add>>=
3018 tcase_add_test(tc_string_eval, test_setup_teardown);
3019 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
3020 tcase_add_test(tc_string_eval, test_simple_eval);
3021 tcase_add_test(tc_string_eval, test_multiple_evals);
3022 tcase_add_test(tc_string_eval, test_eval_with_state);
3023 suite_add_tcase(s, tc_string_eval);
3027 \section{Accelerating function evaluation}
3029 My first version-0.3 code was running very slowly.
3030 With the options suggested in the help
3031 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
3032 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
3033 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
3034 That's only 15 calls per solution, so the search algorithm seems reasonable.
3035 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
3040 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
3042 #endif /* ACCEL_K_H */
3045 <<accel k module makefile lines>>=
3046 NW_SAWSIM_MODS += accel_k
3047 #CHECK_BINS += check_accel_k
3048 check_accel_k_MODS =
3052 #include <assert.h> /* assert() */
3053 #include <stdlib.h> /* realloc(), free(), NULL */
3054 #include "global.h" /* environment_t */
3055 #include "k_model.h" /* k_func_t */
3056 #include "interp.h" /* interp_* */
3057 #include "accel_k.h"
3059 typedef struct accel_k_struct {
3060 interp_table_t *itable;
3066 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3067 static int num_accels = 0;
3068 static accel_k_t *accels=NULL;
3070 /* Wrap k in the standard f(x) acceleration form */
3071 static double k_wrap(double F, void *params)
3073 accel_k_t *p = (accel_k_t *)params;
3074 return (*p->k)(F, p->env, p->params);
3077 static int k_tol(double FA, double kA, double FB, double kB)
3080 if (FB-FA > 1e-12) {
3081 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3082 return 1; /* unnacceptable */
3084 //printf("acceptable tol\n");
3085 return 0; /* acceptable */
3089 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3092 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3093 assert(accels != NULL);
3094 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3096 accels[i].env = env;
3097 accels[i].params = params;
3104 for (i=0; i<num_accels; i++)
3105 interp_table_free(accels[i].itable);
3107 if (accels) free(accels);
3111 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3114 for (i=0; i<num_accels; i++) {
3115 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3116 /* We've already setup this function.
3117 * WARNING: we're assuming param and env are the same. */
3118 return interp_table_eval(accels[i].itable, accels+i, F);
3121 /* set up a new acceleration environment */
3122 i = add_accel_k(k, env, params);
3123 return interp_table_eval(accels[i].itable, accels+i, F);
3127 \section{Tension models}
3128 \label{sec.tension_models}
3130 TODO: tension model intro.
3131 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.
3133 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3134 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]].
3136 <<tension-model.h>>=
3137 #ifndef TENSION_MODEL_H
3138 #define TENSION_MODEL_H
3140 <<tension handler types>>
3141 <<constant tension model declarations>>
3142 <<hooke tension model declarations>>
3143 <<worm-like chain tension model declarations>>
3144 <<freely-jointed chain tension model declarations>>
3145 <<find tension definitions>>
3146 <<tension model global declarations>>
3147 #endif /* TENSION_MODEL_H */
3150 <<tension model module makefile lines>>=
3151 NW_SAWSIM_MODS += tension_model
3152 #CHECK_BINS += check_tension_model
3153 check_tension_model_MODS =
3155 <<tension model utils makefile lines>>=
3156 TENSION_MODEL_MODS = tension_model parse list tension_balance
3157 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3158 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3159 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3160 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
3161 notangle -Rtension-model-utils.c $< > $@
3162 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
3163 gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3164 clean_tension_model_utils :
3165 rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3166 @ The pipe symbol [[|]] marks the prerequisits following it (in this
3167 case, the directories) as ‘order-only’ prerequisites. The timestamp
3168 on these prerequisits does not effect whether the rules are executed.
3169 This is appropriate for directories, where we don't need to recompile
3170 after an unrelated has been added to the directory, but only when the
3171 source prerequisites change. See the [[make]] documentation for more
3173 (\url{http://www.gnu.org/software/make/manual/html_node/Prerequisite-Types.html}).
3176 \label{sec.null_tension}
3178 For unstretchable domains.
3180 <<null tension model getopt>>=
3181 {NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3184 \subsection{Constant}
3185 \label{sec.const_tension}
3187 <<constant tension functions>>=
3188 <<constant tension function>>
3189 <<constant tension parameter create/destroy functions>>
3192 <<constant tension model declarations>>=
3193 <<constant tension function declaration>>
3194 <<constant tension parameter create/destroy function declarations>>
3195 <<constant tension model global declarations>>
3199 An infinitely stretchable domain providing a constant tension.
3201 <<constant tension function declaration>>=
3202 extern double const_tension_handler(double x, void *pdata);
3204 <<constant tension function>>=
3205 double const_tension_handler(double x, void *pdata)
3207 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3208 list_t *list = data->group;
3213 assert(list != NULL); /* empty active group?! */
3214 F = ((const_tension_param_t *)list->d)->F;
3216 fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3218 while (list != NULL) {
3219 assert(((const_tension_param_t *)list->d)->F == F);
3227 <<constant tension structure>>=
3228 typedef struct const_tension_param_struct {
3229 double F; /* tension (force) in N */
3230 } const_tension_param_t;
3234 <<constant tension parameter create/destroy function declarations>>=
3235 extern void *string_create_const_tension_param_t(char **param_strings);
3236 extern void destroy_const_tension_param_t(void *p);
3239 <<constant tension parameter create/destroy functions>>=
3240 const_tension_param_t *create_const_tension_param_t(double F)
3242 const_tension_param_t *ret
3243 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3244 assert(ret != NULL);
3249 void *string_create_const_tension_param_t(char **param_strings)
3251 return create_const_tension_param_t(atof(param_strings[0]));
3254 void destroy_const_tension_param_t(void *p)
3262 <<constant tension model global declarations>>=
3263 extern int num_const_tension_params;
3264 extern char *const_tension_param_descriptions[];
3265 extern char const_tension_param_string[];
3268 <<constant tension model globals>>=
3269 int num_const_tension_params = 1;
3270 char *const_tension_param_descriptions[] = {"tension F, N"};
3271 char const_tension_param_string[] = "0";
3274 <<constant tension model getopt>>=
3275 {&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}
3281 <<hooke functions>>=
3283 <<hooke parameter create/destroy functions>>
3286 <<hooke tension model declarations>>=
3287 <<hooke tension function declaration>>
3288 <<hooke tension parameter create/destroy function declarations>>
3289 <<hooke tension model global declarations>>
3293 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
3294 The behavior of a series of springs $k_i$ in series is given by
3296 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
3298 For a simple proof, see Appendix \ref{app.math_hooke}.
3300 <<hooke tension function declaration>>=
3301 extern double hooke_handler(double x, void *pdata);
3305 double hooke_handler(double x, void *pdata)
3307 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3308 list_t *list = data->group;
3313 assert(list != NULL); /* empty active group?! */
3314 while (list != NULL) {
3315 assert( ((hooke_param_t *)list->d)->k > 0 );
3316 k += 1.0/ ((hooke_param_t *)list->d)->k;
3318 fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
3319 __FUNCTION__, ((hooke_param_t *)list->d)->k);
3325 fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N\n",
3326 __FUNCTION__, k, x, k*x);
3332 <<hooke structure>>=
3333 typedef struct hooke_param_struct {
3334 double k; /* spring constant in N/m */
3338 <<hooke tension parameter create/destroy function declarations>>=
3339 extern void *string_create_hooke_param_t(char **param_strings);
3340 extern void destroy_hooke_param_t(void *p);
3343 <<hooke parameter create/destroy functions>>=
3344 hooke_param_t *create_hooke_param_t(double k)
3346 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
3347 assert(ret != NULL);
3352 void *string_create_hooke_param_t(char **param_strings)
3354 return create_hooke_param_t(atof(param_strings[0]));
3357 void destroy_hooke_param_t(void *p)
3364 <<hooke tension model global declarations>>=
3365 extern int num_hooke_params;
3366 extern char *hooke_param_descriptions[];
3367 extern char hooke_param_string[];
3370 <<hooke tension model globals>>=
3371 int num_hooke_params = 1;
3372 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
3373 char hooke_param_string[]="0.05";
3376 <<hooke tension model getopt>>=
3377 {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}
3380 \subsection{Worm-like chain}
3383 We can model several unfolded domains with a single worm-like chain.
3384 <<worm-like chain functions>>=
3385 <<worm-like chain function>>
3386 <<worm-like chain handler>>
3387 <<worm-like chain create/destroy functions>>
3390 <<worm-like chain tension model declarations>>=
3391 <<worm-like chain handler declaration>>
3392 <<worm-like chain create/destroy function declarations>>
3393 <<worm-like chain tension model global declarations>>
3397 The combination of all unfolded domains is modeled as a worm like chain,
3398 with the force $F_{\text{WLC}}$ approximately given by
3400 F_{\text{WLC}} \approx \frac{k_B T}{p}
3401 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
3403 where $x$ is the distance between the fixed ends,
3404 $k_B$ is Boltzmann's constant,
3405 $T$ is the absolute temperature,
3406 $p$ is the persistence length, and
3407 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
3408 <<worm-like chain function>>=
3409 double wlc(double x, double T, double p, double L)
3411 double xL=0.0; /* uses global k_B */
3413 assert(T > 0); assert(p > 0); assert(L > 0);
3414 if (x >= L) return HUGE_VAL;
3415 xL = x/L; /* unitless */
3416 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
3417 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
3418 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
3421 This model assumes that each unfolded domain has the same persistence length.
3423 <<worm-like chain handler declaration>>=
3424 extern double wlc_handler(double x, void *pdata);
3427 <<worm-like chain handler>>=
3428 double wlc_handler(double x, void *pdata)
3430 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3431 list_t *list = data->group;
3435 assert(list != NULL); /* empty active group?! */
3436 p = ((wlc_param_t *)list->d)->p;
3437 while (list != NULL) {
3438 /* enforce consistent p */
3439 assert( ((wlc_param_t *)list->d)->p == p);
3440 L += ((wlc_param_t *)list->d)->L;
3442 fprintf(stderr, "%s: WLC section %g m long in series\n",
3443 __FUNCTION__, ((wlc_param_t *)list->d)->L);
3448 fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
3449 __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
3451 return wlc(x, data->env->T, p, L);
3455 <<worm-like chain structure>>=
3456 typedef struct wlc_param_struct {
3457 double p; /* WLC persistence length (m) */
3458 double L; /* WLC contour length (m) */
3462 <<worm-like chain create/destroy function declarations>>=
3463 extern void *string_create_wlc_param_t(char **param_strings);
3464 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
3467 <<worm-like chain create/destroy functions>>=
3468 wlc_param_t *create_wlc_param_t(double p, double L)
3470 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
3471 assert(ret != NULL);
3477 void *string_create_wlc_param_t(char **param_strings)
3479 return create_wlc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3482 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
3489 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.
3490 TODO (now we copy the string before parsing, but still do this for clarity).
3492 <<valid string write code>>=
3493 char string[] = "abc";
3496 <<valid string write code>>=
3497 char *string = "abc";
3501 <<worm-like chain tension model global declarations>>=
3502 extern int num_wlc_params;
3503 extern char *wlc_param_descriptions[];
3504 extern char wlc_param_string[];
3506 <<worm-like chain tension model globals>>=
3507 int num_wlc_params = 2;
3508 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
3509 char wlc_param_string[]="0.39e-9,28.4e-9";
3512 <<worm-like chain tension model getopt>>=
3513 {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}
3515 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
3517 \subsection{Freely-jointed chain}
3520 <<freely-jointed chain functions>>=
3521 <<freely-jointed chain function>>
3522 <<freely-jointed chain handler>>
3523 <<freely-jointed chain create/destroy functions>>
3526 <<freely-jointed chain tension model declarations>>=
3527 <<freely-jointed chain handler declaration>>
3528 <<freely-jointed chain create/destroy function declarations>>
3529 <<freely-jointed chain tension model global declarations>>
3533 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
3534 The tension of a chain of $N$ such links is given by
3536 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
3538 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}.
3539 <<freely-jointed chain function>>=
3540 double langevin(double x, void *params)
3543 return 1.0/tanh(x) - 1/x;
3545 <<one dimensional bracket declaration>>
3546 <<one dimensional solve declaration>>
3547 double inv_langevin(double x)
3549 double lb=0.0, ub=1.0, ret;
3550 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
3551 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
3555 double fjc(double x, double T, double l, int N)
3557 double L = l*(double)N;
3558 if (x == 0) return 0;
3559 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
3560 return k_B*T/l * inv_langevin(x/L);
3564 <<freely-jointed chain handler declaration>>=
3565 extern double fjc_handler(double x, void *pdata);
3567 <<freely-jointed chain handler>>=
3568 double fjc_handler(double x, void *pdata)
3570 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3571 list_t *list = data->group;
3576 assert(list != NULL); /* empty active group?! */
3577 l = ((fjc_param_t *)list->d)->l;
3578 while (list != NULL) {
3579 /* enforce consistent l */
3580 assert( ((fjc_param_t *)list->d)->l == l);
3581 N += ((fjc_param_t *)list->d)->N;
3583 fprintf(stderr, "%s: FJC section %d links long in series\n",
3584 __FUNCTION__, ((fjc_param_t *)list->d)->N);
3589 fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
3590 __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
3592 return fjc(x, data->env->T, l, N);
3596 <<freely-jointed chain structure>>=
3597 typedef struct fjc_param_struct {
3598 double l; /* FJC link length (m) */
3599 int N; /* FJC number of links */
3603 <<freely-jointed chain create/destroy function declarations>>=
3604 extern void *string_create_fjc_param_t(char **param_strings);
3605 extern void destroy_fjc_param_t(void *p);
3608 <<freely-jointed chain create/destroy functions>>=
3609 fjc_param_t *create_fjc_param_t(double l, double N)
3611 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
3612 assert(ret != NULL);
3618 void *string_create_fjc_param_t(char **param_strings)
3620 return create_fjc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3623 void destroy_fjc_param_t(void *p)
3630 <<freely-jointed chain tension model global declarations>>=
3631 extern int num_fjc_params;
3632 extern char *fjc_param_descriptions[];
3633 extern char fjc_param_string[];
3636 <<freely-jointed chain tension model globals>>=
3637 int num_fjc_params = 2;
3638 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
3639 char fjc_param_string[]="0.5e-9,200";
3642 <<freely-jointed chain tension model getopt>>=
3643 {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}
3646 \subsection{Tension model implementation}
3648 <<tension-model.c>>=
3651 <<tension model includes>>
3652 #include "tension_model.h"
3653 <<tension model internal definitions>>
3654 <<tension model globals>>
3655 <<tension model functions>>
3658 <<tension model includes>>=
3659 #include <assert.h> /* assert() */
3660 #include <stdlib.h> /* NULL */
3661 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
3662 #include <stdio.h> /* fprintf(), stdout */
3665 #include "tension_balance.h" /* oneD_*() */
3668 <<tension model internal definitions>>=
3669 <<constant tension structure>>
3671 <<worm-like chain structure>>
3672 <<freely-jointed chain structure>>
3675 <<tension model globals>>=
3676 <<tension handler array global>>
3677 <<constant tension model globals>>
3678 <<hooke tension model globals>>
3679 <<worm-like chain tension model globals>>
3680 <<freely-jointed chain tension model globals>>
3683 <<tension model functions>>=
3684 <<constant tension functions>>
3686 <<worm-like chain functions>>
3687 <<freely-jointed chain functions>>
3691 \subsection{Utilities}
3693 The tension models can get complicated, and users may want to reassure
3694 themselves that this portion of the input physics are functioning properly. The
3695 stand-alone program [[tension_model_utils]] demonstrates the tension model
3696 interface without the overhead of having to understand a full simulation.
3697 [[tension_model_utils]] takes command line model arguments like the full
3698 simulation, and prints $F(x)$ data to the screen for the selected model over a
3701 <<tension-model-utils.c>>=
3703 <<tension model utility includes>>
3704 <<tension model utility definitions>>
3705 <<tension model utility getopt functions>>
3707 <<tension model utility main>>
3710 <<tension model utility main>>=
3711 int main(int argc, char **argv)
3713 <<tension model getopt array definition>>
3714 tension_model_getopt_t *model = NULL;
3718 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
3719 tension_handler_data_t tdata;
3720 int num_param_args; /* for INIT_MODEL() */
3721 char **param_args; /* for INIT_MODEL() */
3723 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model,
3724 &Fmax, &Xmax, &flags);
3726 if (flags & VFLAG) {
3727 printf("#initializing model %s with parameters %s\n", model->name, model->params);
3729 INIT_MODEL("utility", model, model->params, params);
3731 push(&tdata.group, params, 1);
3733 tdata.persist = NULL;
3734 if (model->handler == NULL) {
3735 printf("No tension function for model %s\n", model->name);
3741 printf("#Distance (m)\tForce (N)\n");
3742 for (i=0; i<=N; i++) {
3743 x = Xmax*i/(double)N;
3744 F = (*model->handler)(x, &tdata);
3745 if (F < 0 || F > Fmax) break;
3746 printf("%g\t%g\n", x, F);
3748 if (flags & VFLAG && i <= N) {
3749 /* explain exit condition */
3751 printf("Impossible force %g\n", F);
3753 printf("Reached large force limit %g > %g\n", F, Fmax);
3756 params = pop(&tdata.group);
3757 if (model->destructor)
3758 (*model->destructor)(params);
3763 <<tension model utility includes>>=
3766 #include <string.h> /* strlen() */
3767 #include <assert.h> /* assert() */
3771 #include "tension_model.h"
3774 <<tension model utility definitions>>=
3775 <<version definition>>
3776 #define VFLAG 1 /* verbose */
3777 <<string comparison definition and globals>>
3778 <<tension model getopt definitions>>
3779 <<initialize model definition>>
3783 <<tension model utility getopt functions>>=
3784 <<index tension model>>
3785 <<tension model utility help>>
3786 <<tension model utility get options>>
3789 <<tension model utility help>>=
3790 void help(char *prog_name,
3792 int n_tension_models, tension_model_getopt_t *tension_models,
3793 int tension_model, double Fmax, double Xmax)
3796 printf("usage: %s [options]\n", prog_name);
3797 printf("Version %s\n\n", VERSION);
3798 printf("A utility for understanding the available tension models\n\n");
3799 printf("Environmental options:\n");
3800 printf("-T T\tTemperature T (currently %g K)\n", env->T);
3801 printf("-C T\tYou can also set the temperature T in Celsius\n");
3802 printf("Model options:\n");
3803 printf("-m model\tFolded tension model (currently %s)\n",
3804 tension_models[tension_model].name);
3805 printf("-a args\tFolded tension model argument string (currently %s)\n",
3806 tension_models[tension_model].params);
3807 printf("F(x) is calculated for a range of x and printed\n");
3808 printf("For example:\n");
3809 printf(" #Distance (m)\tForce (N)\n");
3810 printf(" 123.456\t7.89\n");
3812 printf("-F\tSet the maximum F value for the standard mode F(x) output (currently %g)\n", Fmax);
3813 printf("-X\tSet the maximum x value for the standard mode F(x) output (currently %g)\n", Xmax);
3814 printf("-V\tChange output to verbose mode\n");
3815 printf("-h\tPrint this help and exit\n");
3817 printf("Tension models:\n");
3818 for (i=0; i<n_tension_models; i++) {
3819 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
3820 for (j=0; j < tension_models[i].num_params ; j++)
3821 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
3822 printf("\t\tdefaults: %s\n", tension_models[i].params);
3827 <<tension model utility get options>>=
3828 void get_options(int argc, char **argv, environment_t *env,
3829 int n_tension_models, tension_model_getopt_t *tension_models,
3830 tension_model_getopt_t **model, double *Fmax, double *Xmax,
3831 unsigned int *flags)
3833 char *prog_name = NULL;
3834 char c, options[] = "T:C:m:a:F:X:Vh";
3835 int tension_model=0, num_strings;
3836 extern char *optarg;
3837 extern int optind, optopt, opterr;
3841 /* setup defaults */
3843 prog_name = argv[0];
3844 env->T = 300.0; /* K */
3848 *model = tension_models;
3850 while ((c=getopt(argc, argv, options)) != -1) {
3852 case 'T': env->T = atof(optarg); break;
3853 case 'C': env->T = atof(optarg)+273.15; break;
3855 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
3856 *model = tension_models+tension_model;
3859 tension_models[tension_model].params = optarg;
3861 case 'F': *Fmax = atof(optarg); break;
3862 case 'X': *Xmax = atof(optarg); break;
3863 case 'V': *flags |= VFLAG; break;
3865 fprintf(stderr, "unrecognized option '%c'\n", optopt);
3866 /* fall through to default case */
3868 help(prog_name, env, n_tension_models, tension_models, tension_model, *Fmax, *Xmax);
3877 \section{Unfolding rate models}
3878 \label{sec.k_models}
3880 We have two main choices for determining $k$: Bell's model, which assumes the
3881 folded domains exist in a single `folded' state and make instantaneous,
3882 stochastic transitions over a free energy barrier; and Kramers' model, which
3883 treats the unfolding as a thermalized diffusion process.
3884 We deal with these options like we dealt with the different tension models: by bundling all model-specific
3885 parameters into structures, with model specific functions for determining $k$.
3887 <<k func definition>>=
3888 typedef double k_func_t(double F, environment_t *env, void *params);
3891 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3892 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]].
3894 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
3895 because \LaTeX\ doesn't like underscores outside math mode.
3896 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
3897 You could use spaces instead (`\verb+ +'),
3898 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
3899 which I don't like as much.
3905 <<k func definition>>
3906 <<null k declarations>>
3907 <<const k declarations>>
3908 <<bell k declarations>>
3909 <<kramers k declarations>>
3910 <<kramers integ k declarations>>
3911 #endif /* K_MODEL_H */
3914 <<k model module makefile lines>>=
3915 NW_SAWSIM_MODS += k_model
3916 CHECK_BINS += check_k_model
3917 check_k_model_MODS = parse string_eval
3919 [[check_k_model]] also depends on the parse module.
3921 <<k model utils makefile lines>>=
3922 K_MODEL_MODS = k_model parse string_eval
3923 K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
3924 $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3925 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
3926 notangle -Rk-model-utils.c $< > $@
3927 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
3928 gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3929 clean_k_model_utils :
3930 rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
3936 For domains that are never folded.
3938 <<null k declarations>>=
3940 <<null k functions>>=
3945 <<null k model getopt>>=
3946 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
3949 \subsection{Constant rate model}
3952 This is a pretty boring model, but it is usefull for testing the engine.
3956 where $k_0$ is the constant unfolding rate.
3958 <<const k functions>>=
3959 <<const k function>>
3960 <<const k structure create/destroy functions>>
3963 <<const k declarations>>=
3964 <<const k function declaration>>
3965 <<const k structure create/destroy function declarations>>
3966 <<const k global declarations>>
3970 <<const k structure>>=
3971 typedef struct const_k_param_struct {
3972 double knot; /* transition rate k_0 (frac population per s) */
3976 <<const k function declaration>>=
3977 double const_k(double F, environment_t *env, void *const_k_params);
3979 <<const k function>>=
3980 double const_k(double F, environment_t *env, void *const_k_params)
3981 { /* Returns the rate constant k in frac pop per s. */
3982 const_k_param_t *p = (const_k_param_t *) const_k_params;
3984 assert(p->knot > 0);
3989 <<const k structure create/destroy function declarations>>=
3990 void *string_create_const_k_param_t(char **param_strings);
3991 void destroy_const_k_param_t(void *p);
3994 <<const k structure create/destroy functions>>=
3995 const_k_param_t *create_const_k_param_t(double knot)
3997 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
3998 assert(ret != NULL);
4003 void *string_create_const_k_param_t(char **param_strings)
4005 return create_const_k_param_t(atof(param_strings[0]));
4008 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
4015 <<const k global declarations>>=
4016 extern int num_const_k_params;
4017 extern char *const_k_param_descriptions[];
4018 extern char const_k_param_string[];
4021 <<const k globals>>=
4022 int num_const_k_params = 1;
4023 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
4024 char const_k_param_string[]="1";
4027 <<const k model getopt>>=
4028 {"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}
4031 \subsection{Bell's model}
4034 Quantitatively, Bell's model gives $k$ as
4036 k = k_0 \cdot e^{F dx / k_B T} \;,
4038 where $k_0$ is the unforced unfolding rate,
4039 $F$ is the applied unfolding force,
4040 $dx$ is the distance to the transition state, and
4041 $k_B T$ is the thermal energy\citep{schlierf06}.
4043 <<bell k functions>>=
4045 <<bell k structure create/destroy functions>>
4048 <<bell k declarations>>=
4049 <<bell k function declaration>>
4050 <<bell k structure create/destroy function declarations>>
4051 <<bell k global declarations>>
4055 <<bell k structure>>=
4056 typedef struct bell_param_struct {
4057 double knot; /* transition rate k_0 (frac population per s) */
4058 double dx; /* `distance to transition state' (nm) */
4062 <<bell k function declaration>>=
4063 double bell_k(double F, environment_t *env, void *bell_params);
4065 <<bell k function>>=
4066 double bell_k(double F, environment_t *env, void *bell_params)
4067 { /* Returns the rate constant k in frac pop per s.
4068 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4069 * uses global k_B in J/K */
4070 bell_param_t *p = (bell_param_t *) bell_params;
4071 assert(F >= 0); assert(env->T > 0);
4073 assert(p->knot > 0); assert(p->dx >= 0);
4075 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
4076 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
4077 p->knot * exp(F*p->dx / (k_B*env->T) ));
4078 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
4080 return p->knot * exp(F*p->dx / (k_B*env->T) );
4084 <<bell k structure create/destroy function declarations>>=
4085 void *string_create_bell_param_t(char **param_strings);
4086 void destroy_bell_param_t(void *p);
4089 <<bell k structure create/destroy functions>>=
4090 bell_param_t *create_bell_param_t(double knot, double dx)
4092 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
4093 assert(ret != NULL);
4099 void *string_create_bell_param_t(char **param_strings)
4101 return create_bell_param_t(atof(param_strings[0]), atof(param_strings[1]));
4104 void destroy_bell_param_t(void *p /* bell_param_t * */)
4111 <<bell k global declarations>>=
4112 extern int num_bell_params;
4113 extern char *bell_param_descriptions[];
4114 extern char bell_param_string[];
4118 int num_bell_params = 2;
4119 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4120 char bell_param_string[]="3.3e-4,0.25e-9";
4123 <<bell k model getopt>>=
4124 {"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}
4126 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4127 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4129 \subsection{Kramer's model}
4132 Kramer's model gives $k$ as
4134 % \frac{1}{k} = \frac{1}{D}
4135 % \int_{x_\text{min}}^{x_\text{max}}
4136 % e^{\frac{-U_F(x)}{k_B T}}
4137 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4140 %where $D$ is the diffusion constant,
4141 %$U_F(x)$ is the free energy along the unfolding coordinate, and
4142 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4143 % before the minimum and after the maximum, respectively \citep{schlierf06}.
4145 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
4147 where $D$ is the diffusion constant,
4149 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
4151 are length scales where
4152 $x_c(F)$ is the location of the bound state and
4153 $x_{ts}(F)$ is the location of the transition state,
4154 $E(F,x)$ is the free energy, and
4155 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
4156 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
4157 \citet{evans97} Eqn.~3).
4159 <<kramers k functions>>=
4160 <<kramers k function>>
4161 <<kramers k structure create/destroy functions>>
4164 <<kramers k declarations>>=
4165 <<kramers k function declaration>>
4166 <<kramers k structure create/destroy function declarations>>
4167 <<kramers k global declarations>>
4171 <<kramers k structure>>=
4172 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
4173 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
4175 typedef struct kramers_param_struct {
4176 double D; /* diffusion rate (um/s) */
4183 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
4184 //kramers_x_func_t *xts; /* function returning transition x (nm) */
4185 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
4186 //kramers_E_func_t *E; /* function returning E (J) */
4187 //double *E_params; /* parametrize the energy landscape */
4188 //int n_E_params; /* length of E_params array */
4192 <<kramers k function declaration>>=
4193 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
4194 extern double kramers_k(double F, environment_t *env, void *kramers_params);
4196 <<kramers k function>>=
4197 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
4199 static char input[160]={0};
4200 static char output[80]={0};
4201 /* setup the environment */
4202 fprintf(in, "F = %g; T = %g\n", F, T);
4203 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4204 string_eval(in, out, input, 80, output);
4205 return atof(output);
4208 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
4210 static char input[160]={0};
4211 static char output[80]={0};
4212 /* setup the environment */
4213 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
4214 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4215 string_eval(in, out, input, 80, output);
4216 return atof(output);
4219 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
4221 kramers_param_t *p = (kramers_param_t *) kramers_params;
4222 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
4225 double kramers_k(double F, environment_t *env, void *kramers_params)
4226 { /* Returns the rate constant k in frac pop per s.
4227 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4228 * uses global k_B in J/K */
4229 kramers_param_t *p = (kramers_param_t *) kramers_params;
4230 double kbT, xc, xts, lc, lts, Eb;
4231 assert(F >= 0); assert(env->T > 0);
4234 assert(p->xc != NULL); assert(p->xts != NULL);
4235 assert(p->ddEdxx != NULL); assert(p->E != NULL);
4236 assert(p->in != NULL); assert(p->out != NULL);
4238 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
4239 if (xc == -1.0) return -1.0; /* error (Too much force) */
4240 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
4241 if (xc == -1.0) return -1.0; /* error (Too much force) */
4242 lc = sqrt(2.0*M_PI*kbT /
4243 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
4244 lts = sqrt(-2.0*M_PI*kbT/
4245 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
4246 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
4247 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
4248 //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));
4249 return p->D/(lc*lts) * exp(-Eb/kbT);
4253 <<kramers k structure create/destroy function declarations>>=
4254 void *string_create_kramers_param_t(char **param_strings);
4255 void destroy_kramers_param_t(void *p);
4258 <<kramers k structure create/destroy functions>>=
4259 kramers_param_t *create_kramers_param_t(double D,
4260 char *xc_fn, char *xts_fn,
4261 char *E_fn, char *ddEdxx_fn,
4262 char *global_define_string)
4263 // kramers_x_func_t *xc, kramers_x_func_t *xts,
4264 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
4265 // double *E_params, int n_E_params)
4267 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
4268 assert(ret != NULL);
4269 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4270 assert(ret->xc != NULL);
4271 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4272 assert(ret->xts != NULL);
4273 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
4274 assert(ret->E != NULL);
4275 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
4276 assert(ret->ddEdxx != NULL);
4278 strcpy(ret->xc, xc_fn);
4279 strcpy(ret->xts, xts_fn);
4280 strcpy(ret->E, E_fn);
4281 strcpy(ret->ddEdxx, ddEdxx_fn);
4282 //ret->E_params = E_params;
4283 //ret->n_E_params = n_E_params;
4284 string_eval_setup(&ret->in, &ret->out);
4285 fprintf(ret->in, "kB = %g\n", k_B);
4286 fprintf(ret->in, "%s\n", global_define_string);
4290 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
4291 void *string_create_kramers_param_t(char **param_strings)
4293 return create_kramers_param_t(atof(param_strings[0]),
4301 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
4303 kramers_param_t *pk = (kramers_param_t *)p;
4305 string_eval_teardown(&pk->in, &pk->out);
4307 // free(pk->E_params);
4313 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
4314 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.
4315 We follow \citet{shillcock98} and use
4317 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
4318 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
4321 where TODO, variable meanings.%\citep{shillcock98}.
4322 The first and second derivatives of $E(F,E_b,x,x_s)$ are
4324 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\\
4325 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
4327 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
4328 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
4329 The roots are therefor at
4331 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
4332 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
4333 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
4336 <<kramers k global declarations>>=
4337 extern int num_kramers_params;
4338 extern char *kramers_param_descriptions[];
4339 extern char kramers_param_string[];
4342 <<kramers k globals>>=
4343 int num_kramers_params = 6;
4344 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"};
4345 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)}";
4347 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
4348 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
4349 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
4350 Returning an $x_c$ position of $-1\U{m}$ lets the simulation know that the rate is poorly defined for that force, and the timestepping algorithm in Section \ref{app.adaptive_dt} reduces it's timestep appropriately.
4351 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
4352 It works so long as [[val_1]] is not [[false]].
4354 <<kramers k model getopt>>=
4355 {"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}
4358 \subsection{Kramer's model (natural cubic splines)}
4359 \label{sec.kramers_integ}
4361 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
4362 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
4363 \citet{schlierf06} Eqn.~2)
4365 \frac{1}{k} = \frac{1}{D}
4366 \int_{x_\text{min}}^{x_\text{max}}
4367 e^{\frac{U_F(x)}{k_B T}}
4368 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4371 where $D$ is the diffusion constant,
4372 $U_F(x)$ is the free energy along the unfolding coordinate, and
4373 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4374 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
4376 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
4377 interpolating between them with cubic splines.
4378 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
4380 <<kramers integ k functions>>=
4381 <<spline functions>>
4382 <<kramers integ A k functions>>
4383 <<kramers integ B k functions>>
4384 <<kramers integ k function>>
4385 <<kramers integ k structure create/destroy functions>>
4388 <<kramers integ k declarations>>=
4389 <<kramers integ k function declaration>>
4390 <<kramers integ k structure create/destroy function declarations>>
4391 <<kramers integ k global declarations>>
4395 <<kramers integ k structure>>=
4396 typedef double func_t(double x, void *params);
4397 typedef struct kramers_integ_param_struct {
4398 double D; /* diffusion rate (um/s) */
4399 double x_min; /* integration bounds */
4401 func_t *E; /* function returning E (J) */
4402 void *E_params; /* parametrize the energy landscape */
4403 destroy_data_func_t *destroy_E_params;
4405 double F; /* for passing into GSL evaluations */
4407 } kramers_integ_param_t;
4410 <<spline param structure>>=
4411 typedef struct spline_param_struct {
4412 int n_knots; /* natural cubic spline knots for U(x) */
4415 gsl_spline *spline; /* GSL spline parameters */
4416 gsl_interp_accel *acc; /* GSL spline acceleration data */
4420 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
4422 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4424 <<kramers integ A k functions>>=
4425 double kramers_integ_k_integrandA(double x, void *params)
4427 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4428 double U = (*p->E)(x, p->E_params);
4429 double Fx = p->F * x;
4430 double kbT = k_B * p->env->T;
4431 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
4432 return exp(-(U-Fx)/kbT);
4434 double kramers_integ_k_integralA(double x, void *params)
4436 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4438 double result, abserr;
4439 size_t neval; /* for qng */
4440 static gsl_integration_workspace *w = NULL;
4441 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4442 f.function = &kramers_integ_k_integrandA;
4444 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4445 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4446 //fprintf(stderr, "integralA = %g\n", result);
4452 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
4453 \int_{x_\text{min}}^{x_\text{max}}
4454 e^{\frac{U_F(x)}{k_B T}}
4455 \text{Integral}_A(x)
4458 <<kramers integ B k functions>>=
4459 double kramers_integ_k_integrandB(double x, void *params)
4461 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4462 double U = (*p->E)(x, p->E_params);
4463 double Fx = p->F * x;
4464 double kbT = k_B * p->env->T;
4465 double intA = kramers_integ_k_integralA(x, params);
4466 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
4467 return exp((U-Fx)/kbT)*intA;
4469 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
4471 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4473 double result, abserr;
4474 size_t neval; /* for qng */
4475 static gsl_integration_workspace *w = NULL;
4476 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4477 f.function = &kramers_integ_k_integrandB;
4481 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4482 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4483 //fprintf(stderr, "integralB = %g\n", result);
4488 With the integrals taken care of, evaluating $k$ is simply
4490 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
4492 <<kramers integ k function declaration>>=
4493 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
4495 <<kramers integ k function>>=
4496 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
4497 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
4498 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
4499 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
4503 <<kramers integ k structure create/destroy function declarations>>=
4504 void *string_create_kramers_integ_param_t(char **param_strings);
4505 void destroy_kramers_integ_param_t(void *p);
4508 <<kramers integ k structure create/destroy functions>>=
4509 kramers_integ_param_t *create_kramers_integ_param_t(double D,
4510 double xmin, double xmax, func_t *E,
4512 destroy_data_func_t *destroy_E_params)
4514 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
4515 assert(ret != NULL);
4520 ret->E_params = E_params;
4521 ret->destroy_E_params = destroy_E_params;
4525 void *string_create_kramers_integ_param_t(char **param_strings)
4527 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
4528 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
4529 void *E_params = string_create_spline_param_t(param_strings+1);
4530 return create_kramers_integ_param_t(atof(param_strings[0]),
4531 atof(param_strings[2]),
4532 atof(param_strings[3]),
4533 &spline_eval, E_params,
4534 destroy_spline_param_t);
4537 void destroy_kramers_integ_param_t(void *params)
4539 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4542 (*p->destroy_E_params)(p->E_params);
4548 Finally we have the GSL spline wrappers:
4549 <<spline functions>>=
4551 <<create/destroy spline>>
4554 <<spline function>>=
4555 double spline_eval(double x, void *spline_params)
4557 spline_param_t *p = (spline_param_t *)(spline_params);
4558 return gsl_spline_eval(p->spline, x, p->acc);
4562 <<create/destroy spline>>=
4563 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
4565 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
4566 assert(ret != NULL);
4567 ret->n_knots = n_knots;
4570 ret->acc = gsl_interp_accel_alloc();
4571 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
4572 gsl_spline_init(ret->spline, x, y, n_knots);
4576 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
4577 void *string_create_spline_param_t(char **param_strings)
4581 double *x=NULL, *y=NULL;
4582 /* break into ordered pairs */
4583 parse_list_string(param_strings[0], SEP, '(', ')',
4584 &num_knots, &knot_coords);
4585 x = (double *)malloc(sizeof(double)*num_knots);
4587 y = (double *)malloc(sizeof(double)*num_knots);
4589 for (i=0; i<num_knots; i++) {
4590 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
4591 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
4594 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
4596 free_parsed_list(num_knots, knot_coords);
4597 return create_spline_param_t(num_knots, x, y);
4600 void destroy_spline_param_t(void *params)
4602 spline_param_t *p = (spline_param_t *)params;
4605 gsl_spline_free(p->spline);
4607 gsl_interp_accel_free(p->acc);
4617 <<kramers integ k global declarations>>=
4618 extern int num_kramers_integ_params;
4619 extern char *kramers_integ_param_descriptions[];
4620 extern char kramers_integ_param_string[];
4623 <<kramers integ k globals>>=
4624 int num_kramers_integ_params = 4;
4625 char *kramers_integ_param_descriptions[] = {
4626 "diffusion D, m^2 / s",
4627 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
4628 "starting integration bound x_min, m",
4629 "ending integration bound x_max, m"};
4630 //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";
4631 //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";
4632 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";
4633 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
4634 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
4635 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
4638 <<kramers integ k model getopt>>=
4639 {"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}
4641 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4642 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4644 \subsection{Unfolding model implementation}
4648 <<k model includes>>
4649 #include "k_model.h"
4650 <<k model internal definitions>>
4652 <<k model functions>>
4655 <<k model includes>>=
4656 #include <assert.h> /* assert() */
4657 #include <stdlib.h> /* NULL, malloc() */
4658 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4659 #include <stdio.h> /* fprintf(), stdout */
4660 #include <string.h> /* strlen(), strcpy() */
4661 #include <gsl/gsl_integration.h>
4662 #include <gsl/gsl_spline.h>
4667 <<k model internal definitions>>=
4668 <<const k structure>>
4669 <<bell k structure>>
4670 <<kramers k structure>>
4671 <<kramers integ k structure>>
4672 <<spline param structure>>
4675 <<k model globals>>=
4679 <<kramers k globals>>
4680 <<kramers integ k globals>>
4683 <<k model functions>>=
4684 <<null k functions>>
4685 <<const k functions>>
4686 <<bell k functions>>
4687 <<kramers k functions>>
4688 <<kramers integ k functions>>
4691 \subsection{Unfolding model unit tests}
4693 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
4694 <<check-k-model.c>>=
4695 <<k model check includes>>
4696 <<check relative error>>
4698 <<k model test suite>>
4699 <<main check program>>
4702 <<k model check includes>>=
4703 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
4704 #include <stdio.h> /* sprintf() */
4705 #include <assert.h> /* assert() */
4706 #include <math.h> /* exp() */
4709 #include "k_model.h"
4712 <<k model test suite>>=
4715 <<k model suite definition>>
4718 <<k model suite definition>>=
4719 Suite *test_suite (void)
4721 Suite *s = suite_create ("k model");
4722 <<const k test case defs>>
4723 <<bell k test case defs>>
4725 <<const k test case adds>>
4726 <<bell k test case adds>>
4731 \subsubsection{Constant}
4733 <<const k test case defs>>=
4734 TCase *tc_const_k = tcase_create("Constant k");
4737 <<const k test case adds>>=
4738 tcase_add_test(tc_const_k, test_const_k_create_destroy);
4739 tcase_add_test(tc_const_k, test_const_k_over_range);
4740 suite_add_tcase(s, tc_const_k);
4745 START_TEST(test_const_k_create_destroy)
4747 char *knot[] = {"1","2","3","4","5","6"};
4748 char *params[] = {knot[0]};
4751 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4752 params[0] = knot[i];
4753 p = string_create_const_k_param_t(params);
4754 destroy_const_k_param_t(p);
4759 START_TEST(test_const_k_over_range)
4761 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
4762 char *knot[] = {"1","2","3","4","5","6"};
4763 char *params[] = {knot[0]};
4766 char param_string[80];
4769 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4770 params[0] = knot[i];
4771 p = string_create_const_k_param_t(params);
4772 for ( F=Fm; F<FM; F+=dF ) {
4773 fail_unless(const_k(F, &env, p)==atof(knot[i]),
4774 "const_k(%g, %g, &{%s}) = %g != %s",
4775 F, T, knot[i], const_k(F, &env, p), knot[i]);
4777 destroy_const_k_param_t(p);
4783 \subsubsection{Bell}
4785 <<bell k test case defs>>=
4786 TCase *tc_bell = tcase_create("Bell's k");
4789 <<bell k test case adds>>=
4790 tcase_add_test(tc_bell, test_bell_k_create_destroy);
4791 tcase_add_test(tc_bell, test_bell_k_at_zero);
4792 tcase_add_test(tc_bell, test_bell_k_at_one);
4793 suite_add_tcase(s, tc_bell);
4797 START_TEST(test_bell_k_create_destroy)
4799 char *knot[] = {"1","2","3","4","5","6"};
4801 char *params[] = {knot[0], dx};
4804 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4805 params[0] = knot[i];
4806 p = string_create_bell_param_t(params);
4807 destroy_bell_param_t(p);
4812 START_TEST(test_bell_k_at_zero)
4814 double F=0.0, T=300.0;
4815 char *knot[] = {"1","2","3","4","5","6"};
4817 char *params[] = {knot[0], dx};
4820 char param_string[80];
4823 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4824 params[0] = knot[i];
4825 p = string_create_bell_param_t(params);
4826 fail_unless(bell_k(F, &env, p)==atof(knot[i]),
4827 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
4828 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
4829 destroy_bell_param_t(p);
4834 START_TEST(test_bell_k_at_one)
4837 char *knot[] = {"1","2","3","4","5","6"};
4839 char *params[] = {knot[0], dx};
4840 double F= k_B*T/atof(dx);
4843 char param_string[80];
4846 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4847 params[0] = knot[i];
4848 p = string_create_bell_param_t(params);
4849 CHECK_ERR(1e-6, atof(knot[i])*exp(1.0), bell_k(F, &env, p));
4850 //fail_unless(bell_k(F, &env, p)==atof(knot[i])*exp(1.0),
4851 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
4852 // F, T, knot[i], dx, bell_k(F, &env, p), atof(knot[i])*exp(1.0));
4853 destroy_bell_param_t(p);
4859 <<kramers k tests>>=
4862 <<kramers k test case def>>=
4865 <<kramers k test case add>>=
4868 <<k model function tests>>=
4869 <<check relative error>>
4871 START_TEST(test_something)
4873 double k=5, x=3, last_x=2.0, F;
4874 one_dim_fn_t *handlers[] = {&hooke};
4875 void *data[] = {&k};
4876 double xi[] = {0.0};
4878 int new_active_groups = 1;
4879 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
4880 fail_unless(F = k*x, NULL);
4885 \subsection{Utilities}
4887 The unfolding models can get complicated, and are sometimes hard to explain to others.
4888 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
4889 the overhead of having to understand a full simulation.
4890 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
4891 or special mode, where the operation depends on the specific model selected.
4893 <<k-model-utils.c>>=
4895 <<k model utility includes>>
4896 <<k model utility definitions>>
4897 <<k model utility getopt functions>>
4898 <<k model utility multi model E>>
4899 <<k model utility main>>
4902 <<k model utility main>>=
4903 int main(int argc, char **argv)
4905 <<k model getopt array definition>>
4906 k_model_getopt_t *model = NULL;
4911 int num_param_args; /* for INIT_MODEL() */
4912 char **param_args; /* for INIT_MODEL() */
4914 double special_xmin;
4915 double special_xmax;
4916 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
4917 &Fmax, &special_xmin, &special_xmax, &flags);
4918 if (flags & VFLAG) {
4919 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4921 INIT_MODEL("utility", model, model->params, params);
4924 if (model->k == NULL) {
4925 printf("No k function for model %s\n", model->name);
4929 printf("Fmax = %g <= 0\n", Fmax);
4935 printf("#F (N)\tk (%% pop. per s)\n");
4936 for (i=0; i<=N; i++) {
4937 F = Fmax*i/(double)N;
4938 k = (*model->k)(F, &env, params);
4940 printf("%g\t%g\n", F, k);
4945 if (model->k == NULL || model->k == &bell_k) {
4946 printf("No special mode for model %s\n", model->name);
4949 if (special_xmax <= special_xmin) {
4950 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
4954 double Xrng=(special_xmax-special_xmin), x, E;
4956 printf("#x (m)\tE (J)\n");
4957 for (i=0; i<=N; i++) {
4958 x = special_xmin + Xrng*i/(double)N;
4959 E = multi_model_E(model->k, params, &env, x);
4960 printf("%g\t%g\n", x, E);
4967 if (model->destructor)
4968 (*model->destructor)(params);
4973 <<k model utility multi model E>>=
4974 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
4976 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
4978 if (k == kramers_integ_k)
4979 E = (*p->E)(x, p->E_params);
4981 E = kramers_E(0, env, model_params, x);
4987 <<k model utility includes>>=
4990 #include <string.h> /* strlen() */
4991 #include <assert.h> /* assert() */
4994 #include "string_eval.h"
4995 #include "k_model.h"
4998 <<k model utility definitions>>=
4999 <<version definition>>
5000 #define VFLAG 1 /* verbose */
5001 enum mode_t {M_K_OF_F, M_SPECIAL};
5002 <<string comparison definition and globals>>
5003 <<k model getopt definitions>>
5004 <<initialize model definition>>
5005 <<kramers k structure>>
5006 <<kramers integ k structure>>
5010 <<k model utility getopt functions>>=
5012 <<k model utility help>>
5013 <<k model utility get options>>
5016 <<k model utility help>>=
5017 void help(char *prog_name,
5019 int n_k_models, k_model_getopt_t *k_models,
5020 int k_model, double Fmax, double special_xmin, double special_xmax)
5023 printf("usage: %s [options]\n", prog_name);
5024 printf("Version %s\n\n", VERSION);
5025 printf("A utility for understanding the available unfolding models\n\n");
5026 printf("Environmental options:\n");
5027 printf("-T T\tTemperature T (currently %g K)\n", env->T);
5028 printf("-C T\tYou can also set the temperature T in Celsius\n");
5029 printf("Model options:\n");
5030 printf("-k model\tTransition rate model (currently %s)\n",
5031 k_models[k_model].name);
5032 printf("-K args\tTransition rate model argument string (currently %s)\n",
5033 k_models[k_model].params);
5034 printf("There are two output modes. In standard mode, k(F) is printed\n");
5035 printf("For example:\n");
5036 printf(" #Force (N)\tk (% pop. per s)\n");
5037 printf(" 123.456\t7.89\n");
5039 printf("In special mode, the output depends on the model.\n");
5040 printf("For models defining an energy function E(x), that function is printed\n");
5041 printf(" #Position (m)\tE (J)\n");
5042 printf(" 0.0012\t0.0034\n");
5044 printf("-m\tChange output to standard mode\n");
5045 printf("-M\tChange output to special mode\n");
5046 printf("-F\tSet the maximum F value for the standard mode k(F) output (currently %g)\n", Fmax);
5047 printf("-x\tSet the minimum x value for the special mode E(x) output (currently %g)\n", special_xmin);
5048 printf("-X\tSet the maximum x value for the special mode E(x) output (currently %g)\n", special_xmax);
5049 printf("-V\tChange output to verbose mode\n");
5050 printf("-h\tPrint this help and exit\n");
5052 printf("Unfolding rate models:\n");
5053 for (i=0; i<n_k_models; i++) {
5054 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
5055 for (j=0; j < k_models[i].num_params ; j++)
5056 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
5057 printf("\t\tdefaults: %s\n", k_models[i].params);
5062 <<k model utility get options>>=
5063 void get_options(int argc, char **argv, environment_t *env,
5064 int n_k_models, k_model_getopt_t *k_models,
5065 enum mode_t *mode, k_model_getopt_t **model,
5066 double *Fmax, double *special_xmin, double *special_xmax,
5067 unsigned int *flags)
5069 char *prog_name = NULL;
5070 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
5072 extern char *optarg;
5073 extern int optind, optopt, opterr;
5077 /* setup defaults */
5079 prog_name = argv[0];
5080 env->T = 300.0; /* K */
5085 *special_xmax = 1e-8;
5086 *special_xmin = 0.0;
5088 while ((c=getopt(argc, argv, options)) != -1) {
5090 case 'T': env->T = atof(optarg); break;
5091 case 'C': env->T = atof(optarg)+273.15; break;
5093 k_model = index_k_model(n_k_models, k_models, optarg);
5094 *model = k_models+k_model;
5097 k_models[k_model].params = optarg;
5099 case 'm': *mode = M_K_OF_F; break;
5100 case 'M': *mode = M_SPECIAL; break;
5101 case 'F': *Fmax = atof(optarg); break;
5102 case 'x': *special_xmin = atof(optarg); break;
5103 case 'X': *special_xmax = atof(optarg); break;
5104 case 'V': *flags |= VFLAG; break;
5106 fprintf(stderr, "unrecognized option '%c'\n", optopt);
5107 /* fall through to default case */
5109 help(prog_name, env, n_k_models, k_models, k_model, *Fmax, *special_xmin, *special_xmax);
5118 \section{\LaTeX\ documentation}
5120 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
5121 The comment blocks are extracted (with nicely formatted code blocks), using
5122 <<latex makefile lines>>=
5123 $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5124 noweave -latex -index -delay $< > $@
5125 $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR)
5129 <<latex makefile lines>>=
5130 $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
5132 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5133 $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
5134 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5135 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5136 mv $(BUILD_DIR)/sawsim.pdf $@
5138 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
5139 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
5140 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
5142 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
5143 <<latex makefile lines>>=
5145 rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
5146 $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
5147 $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.tex \
5148 $(BUILD_DIR)/sawsim.bib
5149 clean_latex : semi-clean_latex
5150 rm -f $(DOC_DIR)/sawsim.pdf
5155 [[make]] is a common utility on *nix systems for managing dependencies while building software.
5156 In this case, we have several \emph{targets} that we'd like to build:
5157 the main simulation program \prog;
5158 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
5159 the documentation [[sawsim.pdf]];
5160 and the [[Makefile]] itself.
5161 Besides the generated files, there is the utility target
5162 [[clean]] for removing all generated files except [[Makefile]].
5164 % [[dist]] for generating a distributable tar file.
5166 Extract the makefile with
5167 `[[notangle -Rmakefile src/sawsim.nw | sed 's/ /\t/' > Makefile]]'.
5168 The sed is needed because notangle replaces tabs with eight spaces,
5169 but [[make]] needs tabs.
5171 # Decide what directories to put things in
5176 # Note: Cannot use, for example, `./build', because make eats the `./'
5177 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5179 # Modules (X.c, X.h) defined in the noweb file
5182 # Modules defined outside the noweb file
5183 FREE_SAWSIM_MODS = interp tavl
5185 <<list module makefile lines>>
5186 <<tension balance module makefile lines>>
5187 <<tension model module makefile lines>>
5188 <<k model module makefile lines>>
5189 <<parse module makefile lines>>
5190 <<string eval module makefile lines>>
5191 <<accel k module makefile lines>>
5193 SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS)
5195 # Everything needed for sawsim under one roof. sawsim.c must come first
5196 SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \
5197 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h)
5198 # Libraries needed to compile sawsim
5199 LIBS = gsl gslcblas m
5200 CHECK_LIBS = gsl gslcblas m check
5201 # The non-check binaries generated
5202 BINS = sawsim tension_model_utils k_model_utils sawsim_profile
5205 # Define the major targets
5206 all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
5208 view : $(DOC_DIR)/sawsim.pdf
5210 profile : $(BIN_DIR)/sawsim_profile | $(BIN_DIR)
5211 $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \
5212 -mnull -Mwlc -A0.39e-9,28e-9 -F8
5213 gprof $(BIN_DIR)/sawsim_profile gmon.out > $@
5214 check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim
5215 $(SHELL) -e -c 'for B in $^; do ./$$B; done'
5216 clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \
5217 clean_tension_model_utils \
5218 clean_k_model_utils clean_latex clean_check_sawsim
5219 rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \
5220 $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \
5221 $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
5222 $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
5223 $(BUILD_DIR)/global.h ./gmon.out
5224 $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
5226 # Various builds of sawsim
5227 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
5228 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5229 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
5230 gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5231 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
5232 gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5234 # Create the directories
5235 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
5238 # Copy over the source external to sawsim
5239 # Note: Cannot use, for example, `./build', because make eats the `./'
5240 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5242 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
5246 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
5250 ## Basic source generated with noweb
5251 # The central files sawsim.c and global.h...
5252 $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5254 $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5255 notangle -Rglobal.h $< > $@
5256 # ... and the modules
5257 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5258 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5259 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5260 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5261 $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5262 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5263 # Note: I use `_' as a space in my file names, but noweb doesn't like
5264 # that so we use `-' to name the noweb roots and substitute here.
5266 ## Building the unit-test programs
5268 $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
5269 notangle -Rchecks $< > $@
5270 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
5271 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
5272 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR)
5273 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
5274 clean_check_sawsim :
5275 rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
5276 # ... and the modules
5278 $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
5279 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \
5280 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \
5281 $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
5282 $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
5283 $$(BUILD_DIR)/global.h | $(BIN_DIR)
5284 gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
5285 $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
5287 # todo: clean up the required modules to
5288 $(CHECK_BINS:%=clean_%) :
5289 rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
5291 # Cleaning up the modules
5293 $(SAWSIM_MODS:%=clean_%) :
5294 rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h)
5296 <<tension model utils makefile lines>>
5297 <<k model utils makefile lines>>
5298 <<latex makefile lines>>
5300 Makefile : $(SRC_DIR)/sawsim.nw
5301 notangle -Rmakefile $< | sed 's/ /\t/' > $@
5303 This is fairly self-explanatory. We compile a dynamically linked
5304 version ([[sawsim]]) and a statically linked version
5305 ([[sawsim_static]]). The static version is about 9 times as big, but
5306 it runs without needing dynamic GSL libraries to link to, so it's a
5307 better format for distributing to a cluster for parallel evaluation.
5311 \subsection{Hookean springs in series}
5312 \label{app.math_hooke}
5315 The effective spring constant for $n$ springs $k_i$ in series is given by
5317 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5323 F &= k_i x_i &&\forall i \le n \\
5324 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
5325 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
5326 F &= k_1 x_1 = k_\text{eff} x \\
5327 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
5328 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5333 \addcontentsline{toc}{section}{References}
5334 \bibliography{sawsim}