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};
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, &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);
226 dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, v);
227 unfolding=random_unfoldings(domain_list, F, dt, &env, &num_folded);
228 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\n", x, F);
231 destroy_domain_list(domain_list);
235 @ The meat of the simulation is bundled into the three functions
236 [[find_tension]], [[determine_dt]], and [[random_unfoldings]].
237 [[find_tension]] is discussed in Section \ref{sec.find_tension},
238 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
239 [[random_unfoldings]] in Section \ref{sec.unfolding_rate}.
241 Environmental parameters important in determining reaction rates and
242 tensions (e.g. temperature, pH) are stored in a single structure to
243 facilitate extension to more complicated models in the future.
244 <<environment definition>>=
245 typedef struct environment_struct {
253 <<environment definition>>
254 <<one dimensional function definition>>
255 <<create/destroy definitions>>
257 #endif /* GLOBAL_H */
259 Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
260 included multiple times.
262 \section{Simulation functions}
264 <<simulation functions>>=
268 <<does domain unfold>>
269 <<randomly unfold domains>>
273 \label{sec.find_tension}
275 Because the stretched system may be made up of several parts (folded
276 domains, unfolded domains, spring-like cantilever, \ldots), we will
277 assign the domains to models and groups. The models are used to
278 determine a domain's tension (Hookean spring, WLC, \ldots). Groups
279 will represent collections of domains which the model should treat as
280 a single entity. A domain may belong to different groups or models
281 depending on its state. For example, a domain might be modeled as a
282 freely-jointed chain when it iss folded and as a worm-like chain class
283 when it is unfolded. The domains are assumed to be commutative, so
284 ordering is ignored. The interactions between the groups are assumed
285 to be linear, but the interactions between domains of the same group
286 need not be. This allows for non-linear group models such as th
287 worm-like or freely-jointed chains. Each model has a tension handler
288 function, which gives the tension $F$ needed to stretch a given group
289 of domains a distance $x$.
291 To understand the purpose of groups, consider a chain of proteins
292 where two folded proteins $A$ and $B$ are modeled as WLCs with
293 persistence length $p_f$ and two unfolded proteins $C$ and $D$ are
294 modeled as WLCs with persistence length $p_u$. The proteins are
295 attached to a cantilever $E$ qof spring constant $κ$. The string
296 would get split into three lists:
298 \begin{tabular}{llll}
299 Model & Group & List & Tension \\
300 WLC & 0 & $[A,B]$ & $F_\text{WLC}(x, p_f, L_A+L_B)$ \\
301 WLC & 1 & $[C,D]$ & $F_\text{WLC}(x, p_u, L_C+L_D)$ \\
302 Hooke & 0 & $[E]$ & $F_\text{Hooke}(x, κ)$ \\
305 Note that group numbers only matter \emph{within} model classes, since
306 grouping domains with seperate models doesn't make sense.
308 <<find tension definitions>>=
309 #define NUM_TENSION_MODELS 5
313 <<tension handler array global declaration>>=
314 extern one_dim_fn_t *tension_handlers[];
317 <<tension handler array global>>=
318 one_dim_fn_t *tension_handlers[] = {
320 &const_tension_handler,
328 <<tension model global declarations>>=
329 <<tension handler array global declaration>>
332 <<tension handler types>>=
333 typedef struct tension_handler_data_struct {
334 /* int num_domains; */
335 /* domain_t *domains; */
339 } tension_handler_data_t;
343 After sorting the chain into separate groups $G_i$, with tension
344 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
345 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
346 \forall i,j$. Note that there may be several groups within each model
347 class. [[find_tension]] is basically a wrapper to feed proper input
348 into the [[tension_balance]] function. [[unfolding]] should be set to
349 0 if there were no unfoldings since the previous call to
352 double find_tension(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list, environment_t *env, double x, int unfolding)
354 static int num_active_groups=0;
355 static one_dim_fn_t **per_group_handlers = NULL;
356 static void **per_group_params = NULL;
357 static double last_x;
358 tension_handler_data_t *tdata;
363 fprintf(stderr, "finding tension at distance %g%s\n", x, unfolding ? " after an unfolding" : "");
364 list_print(stderr, domain_list, "domain list");
367 if (unfolding != 0 || num_active_groups == 0) { /* setup statics */
368 /* free old statics */
369 if (per_group_handlers != NULL)
370 free(per_group_handlers);
371 if (per_group_params != NULL) {
372 for (i=0; i < num_active_groups; i++) {
373 assert(per_group_params[i] != NULL);
374 free(per_group_params[i]);
376 free(per_group_params);
379 /* get new statics */
380 get_active_groups(domain_list, num_tension_handlers, tension_handlers, &num_active_groups, &per_group_handlers, &per_group_params);
382 /* fill in the group handlers and params */
383 for (i=0; i<num_active_groups; i++) {
384 tdata = (tension_handler_data_t *) per_group_params[i];
386 fprintf(stderr, "active group %d of %d", i, num_active_groups);
387 list_print(stderr, tdata->group, " parameters");
390 /* tdata->persist continues unchanged */
393 } /* else roll with the current statics */
395 F = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_params, last_x, x);
399 @ For the moment, we will restrict our group boundaries to a single
400 dimension, so $\sum_i x_i = x$, our total extension (it is this
401 restriction that causes the groups to interact linearly). We'll also
402 restrict our extensions to all be positive. With these restrictions,
403 the problem of balancing the tensions reduces to solving a system of
404 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
405 the number of active groups. In general this can be fairly
406 complicated, but without much loss of practicality we can restrict
407 ourselves to strictly monotonically increasing, non-negative tension
408 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
409 simpler. For example, it guarantees the existence of a unique, real
410 solution for finite forces. See Appendix \ref{app.tension_balance}
411 for the tension balancing implementation.
413 Each group must define a parameter structure if the tension model is
415 a function to create the parameter structure,
416 a function to destroy the parameter structure, and
417 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as it's data argument.
418 The tension-specific parameter structures are contained in the domain
419 group field. See the Section \ref{app.model_selection} for a
420 discussion on the command line framework. See the worm-like chain
421 implementation for an example (Section \ref{sec.wlc}).
423 The major limitation of our approach is that all of our tension models
424 are Markovian (memory-less), which is not really a problem for the
425 chains (see \citet{evans99} for WLC thermalization rate calculations).
427 \subsection{Unfolding rate}
428 \label{sec.unfolding_rate}
430 Each folded domain is modeled as unimolecular, first order reaction
432 \text{Folded} \xrightarrow{k} % in tex: X atop Y
435 With the rate constant $k$ defined by
437 \frac{d[\text{Folded}]}{dt} = -k [\text{Folded}].
439 So the probability of a given protein unfolding in a short time $dt$
445 <<does domain unfold>>=
446 int domain_unfolds(double F, double dt, environment_t *env, domain_t *domain)
447 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to a folded domain */
449 k = accel_k(domain->k, F, env, domain->k_params);
450 //(*domain->k)(F, env, domain->k_params);
451 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
452 return happens(k*dt); /* dice roll for prob. k*dt event */
454 @ [[happens]] is a random decision making function defined in Appendix \ref{app.utils}.
456 Only one domain can unfold in each timestep, because the timescale of
457 a domain unfolding $dt_u$ is assumed to be less than the simulation
458 timestep $dt$, so a domain will completely unfold in a single
459 timestep. We adapt our timesteps to keep the probability of a single
460 domain unfolding low, so the probability of two domains unfolding in
461 the same timestep is negligible. This approach breaks down as the
462 adaptive timestep scheme approaches $dt \sim dt_u$, but $dt_u \sim
463 1\U{$\mu$s}$ for Ig27-like proteins \citep{klimov00}, so this
464 shouldn't be a problem. To reassure yourself, you can ask the
465 simulator to print the smallest timestep that was used with TODO.
466 <<randomly unfold domains>>=
467 int random_unfoldings(list_t *domain_list, double tension,
468 double dt, environment_t *env,
469 int *num_folded_domains)
470 { /* returns 1 if there was an unfolding and 0 otherwise */
471 while (domain_list != NULL) {
472 if (D(domain_list)->state == DS_FOLDED
473 && domain_unfolds(tension, dt, env, D(domain_list))) {
474 if (flags & FLAG_OUTPUT_UNFOLDING_FORCES)
475 fprintf(stdout, "%g\n", tension);
476 D(domain_list)->state = DS_UNFOLDED;
477 (*num_folded_domains)--;
478 return 1; /* our one domain has unfolded, stop looking for others */
480 domain_list = domain_list->next;
486 \subsection{Adaptive timesteps}
487 \label{sec.adaptive_dt}
489 We'd like to pick $dt$ so the probability of unfolding in the next
490 timestep is small. If we don't adapt our timestep, we risk breaking
491 our approximation that $F(x' \in [x, x+v \cdot dt]) = F(x)$ or that
492 only one domain unfolds in a given timestep. Because $F(x)$ is
493 monotonically increasing, excessively large timesteps will lead to
494 erroneously large unfolding forces. The simulation would have been
495 accurate for sufficiently small fixed timesteps, but adaptive
496 timesteps will allow us to move through low tension regions in fewer
497 steps, leading to a more efficient simulation.
499 The actual adaptive timestep implementation is not particularly
500 interesting, since we are only required to reduce $dt$ to somewhere
501 below a set threshold, so I've removed it to Appendix
502 \ref{app.adaptive_dt}.
508 The models provide the physics of the simulation, but the simulation
509 allows interchangeable models, and we are currently focusing on the
510 simulation itself, so we remove the actual model implementation to the
513 The tension models are defined in Appendix \ref{sec.tension_models}
514 and unfolding models are defined in Appendix \ref{sec.k_models}.
517 #define k_B 1.3806503e-23 /* J/K */
521 \section{Command line interface}
523 <<get option functions>>=
525 <<index tension model>>
531 \subsection{Model selection}
532 \label{app.model_selection}
534 The main difficulty with the command line interface in \prog\ is
535 developing an intuitive method for accessing the possibly large number
536 of available models. We'll treat this generally by defining an array
537 of available models, containing their important parameters.
539 <<get options definitions>>=
540 <<model getopt definitions>>
543 <<create/destroy definitions>>=
544 typedef void *create_data_func_t(char **param_strings);
545 typedef void destroy_data_func_t(void *param_struct);
548 <<model getopt definitions>>=
549 <<tension model getopt definitions>>
550 <<k model getopt definitions>>
554 \subsubsection{Tension}
556 To access the various tension models from the command line, we define
557 a structure that contains all the neccessary information about the
559 <<tension model getopt definitions>>=
560 typedef struct tension_model_getopt_struct {
561 one_dim_fn_t *handler; /* fn implementing the model on a group */
562 char *name; /* name string identifying the model */
563 char *description; /* a brief description of the model */
564 int num_params; /* number of extra params for the model */
565 char **param_descriptions; /* descriptions of extra params */
566 char *params; /* default values for the extra params */
567 create_data_func_t *creator; /* param-string -> model-data-struct fn */
568 destroy_data_func_t *destructor; /* fn to free the model data structure */
569 } tension_model_getopt_t;
572 <<tension model getopt array definition>>=
573 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
574 <<null tension model getopt>>,
575 <<constant tension model getopt>>,
576 <<hooke tension model getopt>>,
577 <<worm-like chain tension model getopt>>,
578 <<freely-jointed chain tension model getopt>>
582 \subsubsection{Unfolding rate}
584 <<k model getopt definitions>>=
585 #define NUM_K_MODELS 5
587 typedef struct k_model_getopt_struct {
592 char **param_descriptions;
594 create_data_func_t *creator;
595 destroy_data_func_t *destructor;
599 <<k model getopt array definition>>=
600 k_model_getopt_t k_models[NUM_K_MODELS] = {
601 <<null k model getopt>>,
602 <<const k model getopt>>,
603 <<bell k model getopt>>,
604 <<kramers k model getopt>>,
605 <<kramers integ k model getopt>>
612 void help(char *prog_name, double P, double dt_max, double v,
614 int n_tension_models, tension_model_getopt_t *tension_models,
615 int folded_tension_model, int unfolded_tension_model,
616 int n_k_models, k_model_getopt_t *k_models,
620 printf("usage: %s [options]\n", prog_name);
621 printf("Version %s\n\n", VERSION);
622 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
623 printf("Simulation options:\n");
624 printf("-P P\tTarget probability for dt (currently %g)\n", P);
625 printf("-t dt\tMaximum allowed timestep dt (currently %g)\n", dt_max);
626 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
627 printf("Environmental options:\n");
628 printf("-T T\tTemperature T (currently %g K)\n", env->T);
629 printf("-C T\tYou can also set the temperature T in Celsius\n");
630 printf("Model options:\n");
631 printf("The domains exist in either the folded or unfolded state\n");
632 printf("The following options change the default behavior in each state.\n");
633 printf("-m model[,group]\tFolded tension model (currently %s)\n",
634 tension_models[folded_tension_model].name);
635 printf("-a args\tFolded tension model argument string (currently %s)\n",
636 tension_models[folded_tension_model].params);
637 printf("-M model[,group]\tUnfolded tension model (currently %s)\n",
638 tension_models[unfolded_tension_model].name);
639 printf("-A args\tUnfolded tension model argument string (currently %s)\n",
640 tension_models[unfolded_tension_model].params);
641 printf("The following options change the unfolding rate.\n");
642 printf("-k model\tTransition rate model (currently %s)\n",
643 k_models[k_model].name);
644 printf("-K args\tTransition rate model argument string (currently %s)\n",
645 k_models[k_model].params);
646 printf("Domain creation options:\n");
647 printf("Once you've set up the appropriate default models, you need to add the domains\n");
648 printf("-F n\tAdd n folded domains with the current models\n");
649 printf("-U n\tAdd n unfolded domains with the current models\n");
650 printf("Output mode options:\n");
651 printf("There are two output modes. In standard mode, only the unfolding\n");
652 printf("events are printed. For example:\n");
653 printf(" #Unfolding Force (N)\n");
654 printf(" 123.456e-12\n");
656 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
657 printf(" #Position (m)\tForce (N)\n");
658 printf(" 0.001\t0.0023\n");
660 printf("-V\tChange output to verbose mode\n");
661 printf("-h\tPrint this help and exit\n");
663 printf("Tension models:\n");
664 for (i=0; i<n_tension_models; i++) {
665 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
666 for (j=0; j < tension_models[i].num_params ; j++)
667 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
668 printf("\t\tdefaults: %s\n", tension_models[i].params);
670 printf("Unfolding rate models:\n");
671 for (i=0; i<n_k_models; i++) {
672 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
673 for (j=0; j < k_models[i].num_params ; j++)
674 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
675 printf("\t\tdefaults: %s\n", k_models[i].params);
677 printf("Examples:\n");
678 printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded\n");
679 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);
680 printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded\n");
681 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);
685 \subsection{Get options}
688 void get_options(int argc, char **argv,
689 double *pP, double *pDt_max, double *pV,
691 int n_tension_models, tension_model_getopt_t *tension_models,
692 int n_k_models, k_model_getopt_t *k_models,
693 list_t **pDomain_list, int *pNum_folded)
695 char *prog_name = NULL;
696 char c, options[] = "P:t:v:T:C:m:a:M:A:k:K:F:U:Vh";
697 int ftension_model=0, utension_model=0, k_model=0;
698 char *ftension_params=NULL, *utension_params=NULL;
699 int i, n, ftension_group=0, utension_group=0;
703 extern int optind, optopt, opterr;
708 for (i=0; i<n_tension_models; i++) {
709 fprintf(stderr, "tension model %i %s:\t model's handler %p\n",
710 i, tension_models[i].name, tension_models[i].handler);
711 assert(tension_models[i].handler == tension_handlers[i]);
716 flags = FLAG_OUTPUT_UNFOLDING_FORCES;
718 *pP = 1e-3; /* % pop per s */
719 *pDt_max = 0.001; /* s */
720 *pV = 1e-6; /* m/s */
721 env->T = 300.0; /* K */
723 while ((c=getopt(argc, argv, options)) != -1) {
725 case 'P': *pP = atof(optarg); break;
726 case 't': *pDt_max = atof(optarg); break;
727 case 'v': *pV = atof(optarg); break;
728 case 'T': env->T = atof(optarg); break;
729 case 'C': env->T = atof(optarg)+273.15; break;
731 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
732 assert(num_strings == 1 || num_strings == 2);
733 if (num_strings == 1)
736 ftension_group = atoi(strings[1]);
737 ftension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
738 ftension_params = tension_models[ftension_model].params;
739 free_parsed_list(num_strings, strings);
742 ftension_params = optarg;
745 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
746 assert(num_strings == 1 || num_strings == 2);
747 if (num_strings == 1)
750 utension_group = atoi(strings[1]);
751 utension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
752 utension_params = tension_models[utension_model].params;
753 free_parsed_list(num_strings, strings);
756 utension_params = optarg;
759 k_model = index_k_model(n_k_models, k_models, optarg);
762 k_models[k_model].params = optarg;
767 for (i=0; i<n; i++) {
768 push(pDomain_list, generate_domain(DS_FOLDED,
769 tension_models+ftension_model,
772 tension_models+utension_model,
775 k_models+k_model), 1);
782 for (i=0; i<n; i++) {
783 push(pDomain_list, generate_domain(DS_UNFOLDED,
784 tension_models+ftension_model,
787 tension_models+utension_model,
790 k_models+k_model), 1);
794 flags = FLAG_OUTPUT_FULL_CURVE;
797 fprintf(stderr, "unrecognized option '%c'\n", optopt);
798 /* fall through to default case */
800 help(prog_name, *pP, *pDt_max, *pV, env, n_tension_models, tension_models, ftension_model, utension_model, n_k_models, k_models, k_model);
804 /* check the arguments */
805 assert(*pDomain_list != NULL); /* no domains! */
806 assert(*pP > 0.0); assert(*pP < 1.0);
807 assert(*pDt_max > 0.0);
809 assert(env->T > 0.0);
814 <<index tension model>>=
815 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
818 for (i=0; i<n_models; i++)
819 if (STRMATCH(models[i].name, name))
822 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
830 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
833 for (i=0; i<n_models; i++)
834 if (STRMATCH(models[i].name, name))
837 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
844 <<initialize model definition>>=
845 /* requires int num_param_args and char **param_args in the current scope
847 * INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
848 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
850 #define INIT_MODEL(role, model, param_string, param_pointer) \
852 parse_list_string(param_string, SEP, '{', '}', \
853 &num_param_args, ¶m_args); \
854 if (num_param_args != model->num_params) { \
856 "%s model %s expected %d params,\n", \
857 role, model->name, model->num_params); \
859 "not the %d params in '%s'\n", \
860 num_param_args, param_string); \
861 assert(num_param_args == model->num_params); \
863 if (model->creator) \
864 param_pointer = (*model->creator)(param_args); \
866 param_pointer = NULL; \
867 free_parsed_list(num_param_args, param_args); \
872 void *generate_domain(enum domain_state_t state,
873 tension_model_getopt_t *folded_model,
875 char *folded_param_string,
876 tension_model_getopt_t *unfolded_model,
878 char *unfolded_param_string,
879 k_model_getopt_t *k_model)
881 void *folded_params, *unfolded_params, *k_params;
882 int num_param_args; /* for INIT_MODEL() */
883 char **param_args; /* for INIT_MODEL() */
886 fprintf(stderr, "generating %s ", state==DS_FOLDED ? "folded" : "unfolded");
887 fprintf(stderr, "domain (k: \"%s\", \"%s\", f: \"%s\" (%p) group %d \"%s\", u: \"%s\" (%p) group %d \"%s\")\n",
888 k_model->name, k_model->params,
889 folded_model->name, folded_model->handler, folded_group, folded_param_string,
890 unfolded_model->name, unfolded_model->handler, unfolded_group, unfolded_param_string);
892 INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
893 INIT_MODEL("unfolded", unfolded_model, unfolded_param_string, unfolded_params);
894 INIT_MODEL("k", k_model, k_model->params, k_params);
895 return create_domain(state,
896 k_model->k, k_params, k_model->destructor,
897 folded_model->handler, folded_group, folded_params, folded_model->destructor,
898 unfolded_model->handler, unfolded_group, unfolded_params, unfolded_model->destructor
905 \addcontentsline{toc}{section}{Appendicies}
907 \section{sawsim.c details}
911 The general layout of our simulation code is:
921 We include [[math.h]], so don't forget to link to the libm with `-lm'.
923 #include <assert.h> /* assert() */
924 #include <stdlib.h> /* malloc(), free(), rand() */
925 #include <stdio.h> /* fprintf(), stdout */
926 #include <string.h> /* strlen, strtok() */
927 #include <math.h> /* exp(), M_PI, sqrt() */
928 #include <sys/types.h> /* pid_t (returned by getpid()) */
929 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
932 #include "tension_balance.h"
934 #include "tension_model.h"
940 <<version definition>>
942 <<max/min definitions>>
943 <<string comparison definition and globals>>
944 <<initialize model definition>>
945 <<get options definitions>>
946 <<domain definitions>>
955 <<create/destroy domain>>
956 <<destroy domain list>>
958 <<simulation functions>>
959 <<boilerplate functions>>
962 <<boilerplate functions>>=
964 <<get option functions>>
970 srand(getpid()*time(NULL)); /* seed rand() */
974 <<flag definitions>>=
975 /* in octal b/c of prefixed '0' */
976 #define FLAG_OUTPUT_FULL_CURVE 01
977 #define FLAG_OUTPUT_UNFOLDING_FORCES 02
981 static unsigned long int flags = 0;
984 \subsection{Utilities}
987 <<max/min definitions>>=
988 #define MAX(a,b) ((a)>(b) ? (a) : (b))
989 #define MIN(a,b) ((a)<(b) ? (a) : (b))
992 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
993 <<string comparison definition and globals>>=
994 // Check if two strings match, return 1 if they do
995 static char *temp_string_A;
996 static char *temp_string_B;
997 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
998 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
999 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
1000 /* +1 to also compare the '\0' */
1003 We also define a macro for our [[check]] unit testing
1004 <<check relative error>>=
1005 #define CHECK_ERR(max_err, expected, received) \
1007 fail_unless( (received-expected)/expected < max_err, \
1008 "relative error %g >= %g in %s (Expected %g, received %g)", \
1009 (received-expected)/expected, max_err, #received, \
1010 expected, received); \
1011 fail_unless(-(received-expected)/expected < max_err, \
1012 "relative error %g >= %g in %s (Expected %g, received %g)", \
1013 -(received-expected)/expected, max_err, #received, \
1014 expected, received); \
1019 int happens(double probability)
1021 assert(probability >= 0.0); assert(probability <= 1.0);
1022 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*/
1027 \subsection{Adaptive timesteps}
1028 \label{app.adaptive_dt}
1030 $F(x)$ increases with $x$, possibly exploding, as in the worm-like chain model,
1031 so basing the timestep on the the unfolding probability at the current tension
1032 is dangerous, and we need to search for a $dt$ for which
1033 $P(F(x+v*dt)) < P_\text{target}$.
1034 There are two cases to consider.
1035 In the most common, no domains have unfolded since the last step,
1036 and we expect the next step to be slightly shorter than the previous one.
1037 In the less common, domains did unfold in the last step,
1038 and we expect the next step to be considerably longer than the previous one.
1040 double search_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers,
1041 list_t *domain_list,
1042 environment_t *env, double x,
1043 double target_prob, double max_dt, double v)
1044 { /* Returns the timestep dt in seconds for the current folded domain.
1045 * Takes a list of tension handlers, the list of domains,
1046 * a pointer env to the environmental data, a starting separation x in m,
1047 * a target_prob between 0 and 1,
1048 * max_dt in s, stretching velocity v in m/s.
1050 double F, k, dtCur, dtU, dtUCur, dtL, dt;
1052 /* get upper bound using the current position */
1053 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x, 0); /* BUG. repeated calculation */
1054 //printf("Start with x = %g (F = %g)\n", x, F);
1055 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1056 //printf("x %g\tF %g\tk %g\n", x, F, k);
1057 dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */
1059 //printf("overshot max_dt\n");
1062 /* set a lower bound on dt too */
1065 /* The dt determined above may produce illegitimate forces or ks.
1066 * Reduce the upper bound until we have valid ks. */
1068 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1069 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1072 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1074 //printf("Try for dt = %g (F = %g)\n", dt, F);
1075 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1076 /* returned k may be -1.0 */
1077 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1078 while (k == -1.0) { /* reduce step until we hit a valid k */
1080 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1081 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1082 //printf("Try for dt = %g (F = %g)\n", dt, F);
1083 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1084 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1086 assert(dtU > 1e-14); /* timestep to valid k too small */
1087 dtUCur = target_prob / k; /* safe timestep back from x+dtU */
1089 return dt; /* dtU is safe. */
1091 /* dtU wasn't safe, lets see what would be. */
1092 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1093 dt = (dtU + dtL) / 2.0;
1094 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1095 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1096 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1097 dtCur = target_prob / k;
1098 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1099 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1101 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1102 dtU = dt; /* ... stepping out only dtCur would be safe */
1105 break; /* dtCur = dt */
1107 return MAX(dtUCur, dtL);
1111 To determine $dt$ for an array of potentially different folded domains,
1112 we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
1115 double determine_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list,
1116 environment_t *env, double x,
1117 double target_prob, double dt_max, double v)
1118 { /* Returns the timestep dt in seconds.
1119 * Takes the list of folded domains, target_prob between 0 and 1,
1120 * F in N, and T in K. */
1121 double dt=dt_max, new_dt;
1122 assert(target_prob > 0.0); assert(target_prob < 1.0);
1123 assert(dt_max > 0.0);
1125 /* .5 nm steps = v * dt */
1127 while (domain_list != NULL) {
1128 if (D(domain_list)->state == DS_FOLDED) {
1129 new_dt = search_dt(num_tension_handlers, tension_handlers, domain_list, env, x, target_prob, dt, v);
1130 dt = MIN(dt, new_dt);
1132 domain_list = domain_list->next;
1138 \subsection{Domain data}
1140 Currently domains exist in two states, folded and unfolded, and the
1141 only allowed transitions are folded $\rightarrow$ unfolded. Of
1142 course, it wouldn't be too complicated to extent this to a multi-state
1143 system, with an array containing the domains group for each possible
1144 state, and a matrix of transition-rate-calculating functions.
1145 However, at this point such generality seems unnecessary at this
1147 <<domain definitions>>=
1148 enum domain_state_t {DS_FOLDED,
1152 typedef struct domain_struct {
1153 enum domain_state_t state;
1154 one_dim_fn_t *folded_handler;
1156 one_dim_fn_t *unfolded_handler;
1158 k_func_t *k; /* function returning unfolding rate */
1159 void *folded_params; /* pointer to folded parameters */
1160 void *unfolded_params; /* pointer to unfolded parameters */
1161 void *k_params; /* pointer to k parameters */
1162 destroy_data_func_t *destroy_folded;
1163 destroy_data_func_t *destroy_unfolded;
1164 destroy_data_func_t *destroy_k;
1167 /* get the domain data for the current list node */
1168 #define D(list) ((domain_t *)(list)->d)
1169 /* get the tension params for the current list node */
1170 #define D_TP(list) (((domain_t *)(list)->d)->state == DS_FOLDED \
1171 ? ((domain_t *)(list)->d)->folded_params \
1172 : ((domain_t *)(list)->d)->unfolded_params)
1174 [[k]] is a pointer to the function determining the unfolding rate for a given tension.
1175 [[folded_params]] is a pointer to the parameters used by the function pointed to by [[k]].
1176 [[unfolded_params]] is a pointer to the parameters used by the group-appropriate handler function when determining the tension.
1177 The [[destroy_*]] pointers point to functions for freeing the memory [[*_params]].
1178 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.
1180 [[create_]] and [[destroy_domain]] are simple wrappers around [[malloc]] and [[free]].
1181 <<create/destroy domain>>=
1182 domain_t *create_domain(enum domain_state_t state,
1185 destroy_data_func_t *destroy_k,
1186 one_dim_fn_t *folded_handler,
1188 void *folded_params,
1189 destroy_data_func_t *destroy_folded,
1190 one_dim_fn_t *unfolded_handler,
1192 void *unfolded_params,
1193 destroy_data_func_t *destroy_unfolded)
1195 domain_t *ret = (domain_t *)malloc(sizeof(domain_t));
1196 assert(ret != NULL);
1197 if (state == DS_FOLDED) {
1198 assert(k != NULL); /* the pointer points somewhere valid */
1199 assert(*k != NULL); /* and there is something useful there */
1201 assert(state == DS_UNFOLDED);
1203 ret->folded_handler = folded_handler;
1204 ret->folded_group = folded_group;
1205 ret->unfolded_handler = unfolded_handler;
1206 ret->unfolded_group = unfolded_group;
1208 ret->folded_params = folded_params;
1209 ret->unfolded_params = unfolded_params;
1210 ret->k_params = k_params;
1211 ret->destroy_folded = destroy_folded;
1212 ret->destroy_unfolded = destroy_unfolded;
1213 ret->destroy_k = destroy_k;
1215 fprintf(stderr, "generated domain %p. state %d\n", ret, ret->state);
1220 void destroy_domain(domain_t *domain)
1223 //printf("domain %p & %p\n", *domain, domain);
1224 if (domain->destroy_folded)
1225 (*domain->destroy_folded)(domain->folded_params);
1226 if (domain->destroy_unfolded)
1227 (*domain->destroy_unfolded)(domain->unfolded_params);
1228 if (domain->destroy_k)
1229 (*domain->destroy_k)(domain->k_params);
1235 <<destroy domain list>>=
1236 void destroy_domain_list(list_t *domain_list)
1238 domain_list = head(domain_list);
1239 while (domain_list != NULL)
1240 destroy_domain((domain_t *) pop(&domain_list));
1244 \subsection{Domain model and group handling}
1246 <<group functions>>=
1247 <<get tension handler>>
1249 <<int comparison function>>
1250 <<find possible groups>>
1253 <<get active groups>>
1256 <<get tension handler>>=
1257 one_dim_fn_t *get_tension_handler(domain_t *domain)
1259 if (domain->state == DS_FOLDED)
1260 return domain->folded_handler;
1262 if (domain->state != DS_UNFOLDED)
1263 fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
1264 assert(domain->state == DS_UNFOLDED);
1265 return domain->unfolded_handler;
1271 int get_group(domain_t *domain)
1273 if (domain->state == DS_FOLDED)
1274 return domain->folded_group;
1276 if (domain->state != DS_UNFOLDED)
1277 fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
1278 assert(domain->state == DS_UNFOLDED);
1279 return domain->unfolded_group;
1284 We already know all possible tension classes, since they are all
1285 hardcoded into \prog. However, there may be any number of possible
1286 groups. We define a function [[find_possible_groups]] to search for
1287 possible (not neccessarily active) groups. It can search for
1288 subgroups of a single [[handler]], or by setting [[handler]] to
1289 [[NULL]], subgroups of \emph{any} handler. It returns a list of group
1291 <<find possible groups>>=
1292 list_t *find_possible_groups(list_t *list, one_dim_fn_t *handler) {
1295 while (list != NULL) {
1296 if (handler == NULL || get_tension_handler(D(list)) == handler) {
1297 push(&ret, &D(list)->folded_group, 1);
1298 push(&ret, &D(list)->unfolded_group, 1);
1304 return ret; /* no groups with this handler, no processing needed */
1306 /* sort the ret list, and remove duplicates */
1307 sort(&ret, &int_comparison);
1308 unique(&ret, &int_comparison);
1310 fprintf(stderr, "handler %p possible groups:", handler);
1312 while (list != NULL) {
1313 fprintf(stderr, " %d", *((int *)list->d));
1316 fprintf(stderr, "\n");
1322 Search a [[list]] of domains, and determine whether a particular model
1323 class and group number combination is active.
1324 <<is group active>>=
1325 int is_group_active(list_t *list, one_dim_fn_t *handler, int group)
1328 while (list != NULL) {
1329 if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group)
1337 Search a [[list]] of domains, and return all domains matching a
1338 particular model class and group number.
1340 list_t *get_group_list(list_t *list, one_dim_fn_t *handler, int group)
1344 while (list != NULL) {
1345 if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group) {
1346 push(&ret, D_TP(list), 1); /* add a pointer to the appropriate tension parameters to our new list. */
1348 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);
1356 Because all the node data in lists returned by [[get_group_list]] is
1357 also in the main domain list, you shouldn't destroy the node data
1358 popped off when destroying the group lists. It will all get cleaned
1359 up when the main domain list is destroyed.
1361 Put all this together to scan through a list of domains and construct
1362 an array of all the active groups.
1363 <<get active groups>>=
1364 void get_active_groups(list_t *list,
1365 int num_tension_handlers, one_dim_fn_t **tension_handlers,
1366 int *pNum_active_groups, one_dim_fn_t ***pPer_group_handlers, void ***pActive_groups)
1368 void **active_groups=NULL;
1369 one_dim_fn_t *handler, **per_group_handlers=NULL;
1370 int i, num_active_groups, *group;
1371 list_t *possible_group_numbers, *active_group_list, *active_groups_list=NULL, *per_group_handlers_list=NULL;
1373 /* get all the active groups in a list */
1374 for (i=0; i<num_tension_handlers; i++) {
1375 handler = tension_handlers[i];
1376 if (handler == NULL) continue; /* tensionless handler */
1377 possible_group_numbers = head(find_possible_groups(list, handler));
1378 while (possible_group_numbers != NULL) {
1379 group = (int *)pop(&possible_group_numbers);
1380 if (is_group_active(list, handler, *group) == 1) {
1381 active_group_list = get_group_list(list, handler, *group);
1382 push(&active_groups_list, active_group_list, 1);
1383 push(&per_group_handlers_list, handler, 1);
1385 fprintf(stderr, "active group %p for %p %d headed by tension parameters %p\n", active_group_list, handler, *group, head(active_group_list)->d);
1386 list_print(stderr, active_group_list, "active group");
1392 /* allocate the array we'll be returning */
1393 num_active_groups = list_length(active_groups_list);
1394 active_groups = (void **) malloc(sizeof(void *)*num_active_groups);
1395 assert(active_groups != NULL);
1396 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1397 assert(per_group_handlers != NULL);
1399 /* populate the array we'll be returning */
1400 active_groups_list = head(active_groups_list);
1401 for (i=0; i<num_active_groups; i++) {
1402 handler = pop(&per_group_handlers_list);
1403 per_group_handlers[i] = handler;
1405 active_groups[i] = (void *) malloc(sizeof(tension_handler_data_t));
1406 assert(active_groups[i] != NULL);
1407 active_group_list = pop(&active_groups_list);
1408 ((tension_handler_data_t *)active_groups[i])->group = active_group_list;
1409 ((tension_handler_data_t *)active_groups[i])->env = NULL;
1410 ((tension_handler_data_t *)active_groups[i])->persist = NULL;
1412 assert(active_groups_list == NULL);
1413 assert(per_group_handlers_list == NULL);
1415 *pNum_active_groups = num_active_groups;
1416 *pPer_group_handlers = per_group_handlers;
1417 *pActive_groups = active_groups;
1422 \section{String parsing}
1424 For handling command line arguments, we need parse delimited strings (e.g. [["param1,param2,param3..."]]).
1425 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1426 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1427 need to take care of parsing those parameters themselves.
1428 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1434 <<parse definitions>>
1435 <<parse declarations>>
1439 <<parse module makefile lines>>=
1440 NW_SAWSIM_MODS += parse
1441 CHECK_BINS += check_parse
1445 <<parse definitions>>=
1446 #define SEP ',' /* argument separator character */
1449 <<parse declarations>>=
1450 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1451 int *num, char ***string_array);
1452 extern void free_parsed_list(int num, char **string_array);
1455 [[parse_list_string]] allocates memory, don't forget to free it
1456 afterward with [[free_parsed_list]]. It does not alter the original.
1458 The string may start off with a [[deeper]] character
1459 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1460 [[x,y]], where the pointer is one character in on the copied string.
1461 However, when we go to free the memory, we need a pointer to the
1462 beginning of the string. In order to accommodate this for a string
1463 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1464 the first $N$ elements point to the separated arguments, and let the
1465 last element point to the start of the copied string regardless of
1467 <<parse delimited list string functions>>=
1468 /* TODO, split out into parse.hc */
1469 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1472 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1473 if (string[i] == deeper) {depth++;}
1474 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1480 void parse_list_string(char *string, char sep, char deeper, char shallower,
1481 int *num, char ***string_array)
1483 char *str=NULL, **ret=NULL;
1485 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
1487 *string_array = NULL;
1490 /* make a copy of the string, so we don't change the original */
1491 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1492 assert(str != NULL);
1493 strcpy(str, string); /* we know str is long enough */
1494 /* count the number of regions, so we can allocate pointers to them */
1497 n++; i++; /* move on to next argument */
1498 i += next_delim_index(str+i, sep, deeper, shallower);
1499 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1500 } while (str[i] != '\0');
1501 ret = (char **)malloc(sizeof(char *)*(n+1));
1502 assert(ret != NULL);
1503 /* replace the separators with '\0' & assign pointers */
1504 ret[n] = str; /* point to the front of the copied string */
1507 for(i=1; i<n; i++) {
1508 j += next_delim_index(str+j, sep, deeper, shallower);
1510 ret[i] = str+j; /* point to the separated arguments */
1512 /* strip off bounding braces, if any */
1513 for(i=0; i<n; i++) {
1514 if (ret[i][0] == deeper) {
1518 if (ret[i][strlen(ret[i])-1] == shallower)
1519 ret[i][strlen(ret[i])-1] = '\0';
1522 *string_array = ret;
1525 void free_parsed_list(int num, char **string_array)
1528 /* frees all items, since they were allocated as a single string*/
1529 free(string_array[num]);
1530 /* and free the array of pointers */
1536 \subsection{Parsing implementation}
1542 <<parse delimited list string functions>>
1546 #include <assert.h> /* assert() */
1547 #include <stdlib.h> /* NULL */
1548 #include <stdio.h> /* fprintf(), stdout *//*!!*/
1549 #include <string.h> /* strlen() */
1553 \subsection{Parsing unit tests}
1555 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1557 <<parse check includes>>
1558 <<string comparison definition and globals>>
1559 <<check relative error>>
1560 <<parse test suite>>
1561 <<main check program>>
1564 <<parse check includes>>=
1565 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
1566 #include <stdio.h> /* printf() */
1567 #include <assert.h> /* assert() */
1568 #include <string.h> /* strlen() */
1573 <<parse test suite>>=
1574 <<parse list string tests>>
1575 <<parse suite definition>>
1578 <<parse suite definition>>=
1579 Suite *test_suite (void)
1581 Suite *s = suite_create ("k model");
1582 <<parse list string test case defs>>
1584 <<parse list string test case add>>
1589 <<parse list string tests>>=
1592 START_TEST(test_next_delim_index)
1594 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
1595 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
1596 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
1597 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
1598 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
1602 START_TEST(test_parse_list_null)
1606 parse_list_string(NULL, SEP, '{', '}',
1607 &num_param_args, ¶m_args);
1608 fail_unless(num_param_args == 0, NULL);
1609 fail_unless(param_args == NULL, NULL);
1612 START_TEST(test_parse_list_single_simple)
1616 parse_list_string("arg", SEP, '{', '}',
1617 &num_param_args, ¶m_args);
1618 fail_unless(num_param_args == 1, NULL);
1619 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
1622 START_TEST(test_parse_list_single_compound)
1626 parse_list_string("{x,y,z}", SEP, '{', '}',
1627 &num_param_args, ¶m_args);
1628 fail_unless(num_param_args == 1, NULL);
1629 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
1632 START_TEST(test_parse_list_double_simple)
1636 parse_list_string("abc,def", SEP, '{', '}',
1637 &num_param_args, ¶m_args);
1638 fail_unless(num_param_args == 2, NULL);
1639 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
1640 fail_unless(STRMATCH(param_args[1],"def"), NULL);
1644 <<parse list string test case defs>>=
1645 TCase *tc_parse_list_string = tcase_create("parse list string");
1647 <<parse list string test case add>>=
1648 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
1649 tcase_add_test(tc_parse_list_string, test_parse_list_null);
1650 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
1651 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
1652 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
1653 suite_add_tcase(s, tc_parse_list_string);
1657 \section{Unit tests}
1659 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1666 <<check relative error>>
1669 <<main check program>>
1681 <<determine dt tests>>
1683 <<does domain unfold tests>>
1684 <<randomly unfold domains tests>>
1685 <<suite definition>>
1688 <<suite definition>>=
1689 Suite *test_suite (void)
1691 Suite *s = suite_create ("sawsim");
1692 <<F test case defs>>
1693 <<determine dt test case defs>>
1694 <<happens test case defs>>
1695 <<does domain unfold test case defs>>
1696 <<randomly unfold domains test case defs>>
1699 <<determine dt test case add>>
1700 <<happens test case add>>
1701 <<does domain unfold test case add>>
1702 <<randomly unfold domains test case add>>
1705 tcase_add_checked_fixture(tc_strip_address,
1706 setup_strip_address,
1707 teardown_strip_address);
1713 <<main check program>>=
1717 Suite *s = test_suite();
1718 SRunner *sr = srunner_create(s);
1719 srunner_run_all(sr, CK_ENV);
1720 number_failed = srunner_ntests_failed(sr);
1722 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
1726 \subsection{F tests}
1728 <<worm-like chain tests>>
1730 <<F test case defs>>=
1731 <<worm-like chain test case def>>
1733 <<F test case add>>=
1734 <<worm-like chain test case add>>
1737 <<worm-like chain tests>>=
1738 extern double wlc(double x, double T, double p, double L);
1739 START_TEST(test_wlc_at_zero)
1741 double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
1742 fail_unless(abs(wlc(x, T, p, L)) < lim, \
1743 "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
1744 x, T, p, L, abs(wlc(x, T, p, L)), lim);
1748 START_TEST(test_wlc_at_half)
1750 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
1751 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
1752 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
1754 * linear term = x/L = 0.5
1755 * nonlinear + linear = 0.75 + 0.5 = 1.25
1756 * wlc = 10e21*1.25 = 12.5e21
1758 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
1759 "wlc(%g, %g, %g, %g) = %g != %g",
1760 x, T, p, L, wlc(x, T, p, L), 12.5e21);
1764 <<worm-like chain test case def>>=
1765 TCase *tc_wlc = tcase_create("WLC");
1768 <<worm-like chain test case add>>=
1769 tcase_add_test(tc_wlc, test_wlc_at_zero);
1770 tcase_add_test(tc_wlc, test_wlc_at_half);
1771 suite_add_tcase(s, tc_wlc);
1774 \subsection{Model tests}
1776 Check the searching with [[linear_k]].
1777 Check overwhelming force treatment with the heavyside-esque step [[?]].
1778 <<determine dt tests>>=
1779 double linear_k(double F, environment_t *env, void *params)
1781 double Fnot = *(double *)params;
1785 START_TEST(test_determine_dt_linear_k)
1788 double dt_max=3.0, Fnot=3.0;
1789 double F[]={0,1,2,3,4,5,6};
1790 domain_t dom; /* use both parts at once for folded/unfolded */
1794 dom->next = dom->prev = NULL;
1795 dom->k_func_t = linear_k;
1796 dom->folded_params = &Fnot;
1797 dom->unfolded_params = !!!!!!!!!
1798 dom->destroy_folded = dom->destroy_unfolded = NULL;
1799 for( i=0; i < sizeof(F)/sizeof(double); i++) {
1800 fail_unless(determine_dt(folded, unfolded, targ_p, dt_max, v, x, T)==1, NULL);
1806 <<determine dt test case defs>>=
1807 TCase *tc_determine_dt = tcase_create("Determine dt");
1809 <<determine dt test case add>>=
1810 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
1811 suite_add_tcase(s, tc_determine_dt);
1816 <<happens test case defs>>=
1818 <<happens test case add>>=
1821 <<does domain unfold tests>>=
1823 <<does domain unfold test case defs>>=
1825 <<does domain unfold test case add>>=
1828 <<randomly unfold domains tests>>=
1830 <<randomly unfold domains test case defs>>=
1832 <<randomly unfold domains test case add>>=
1836 \section{Balancing group extensions}
1837 \label{app.tension_balance}
1839 For a given total extension $x$ (of the piezo), the various domain
1840 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
1841 amounts, and we need to tweak the portion that each extends to balance
1842 the tension amongst the active groups. Since the tension balancing is
1843 separable from the bulk of the simulation, we break it out into a
1844 separate module. The interface is defined in [[tension_balance.h]],
1845 the implementation in [[tension_balance.c]], and the unit testing in
1846 [[check_tension_balance.c]]
1848 <<tension-balance.h>>=
1849 #ifndef TENSION_BALANCE_H
1850 #define TENSION_BALANCE_H
1852 <<tension balance function declaration>>
1853 #endif /* TENSION_BALANCE_H */
1856 <<tension balance functions>>=
1857 <<one dimensional bracket>>
1858 <<one dimensional solve>>
1860 <<tension balance function>>
1863 <<tension balance module makefile lines>>=
1864 NW_SAWSIM_MODS += tension_balance
1865 CHECK_BINS += check_tension_balance
1866 check_tension_balance_MODS =
1869 The entire force balancing problem reduces to a solving two nested
1870 one-dimensional root-finding problems. First we define the one
1871 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
1872 non-empty group). $x(x_0)$ is also strictly monotonically increasing,
1873 so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
1874 stores the last successful value of $x$, since we expect to be taking
1875 small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
1876 indicates that the groups have changed.
1877 <<tension balance function declaration>>=
1878 double tension_balance(int num_tension_groups,
1879 one_dim_fn_t **tension_handlers,
1884 <<tension balance function>>=
1885 double tension_balance(int num_tension_groups,
1886 one_dim_fn_t **tension_handlers,
1891 static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL};
1892 double F, xo_guess, xo, lb, ub;
1893 double min_relx=1e-6, min_rely=1e-6;
1894 int max_steps=200, i;
1896 assert(num_tension_groups > 0);
1897 assert(tension_handlers != NULL);
1898 assert(params != NULL);
1899 assert(num_tension_groups > 0);
1901 if (num_tension_groups == 1) { /* only one group, no balancing required */
1904 //fprintf(stderr, "balancing for x = %g with ", x);
1905 //for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
1906 //fprintf(stderr, "\n");
1907 if (last_x == -1) { /* new group setup, reset x_xo_data */
1908 /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
1909 if (x_xo_data.xi != NULL)
1911 /* construct new x_xo_data */
1912 x_xo_data.n_groups = num_tension_groups;
1913 x_xo_data.tension_handler = tension_handlers;
1914 x_xo_data.handler_data = params;
1915 x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
1916 for (i=0; i<num_tension_groups; i++)
1917 x_xo_data.xi[i] = -1.0;
1919 if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
1920 if (x == 0) xo_guess = length_scale;
1921 else xo_guess = x/num_tension_groups;
1923 fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
1925 oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
1926 } else { /* work off the last balanced point */
1928 fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
1931 lb = x_xo_data.xi[0];
1932 ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */
1933 } else if (x < last_x) {
1934 lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
1935 ub = x_xo_data.xi[0];
1936 } else { /* x == last_x */
1937 fprintf(stderr, "not moving\n");
1938 lb= x_xo_data.xi[0];
1939 ub= x_xo_data.xi[0];
1943 fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
1944 __FUNCTION__, x, lb, ub);
1946 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
1947 F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
1949 if (num_tension_groups > 1) {
1950 fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
1951 for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
1952 fprintf(stderr, "\n");
1956 F = (*tension_handlers[0])(xo, params[0]);
1961 <<tension balance internal definitions>>=
1962 <<x of x0 definitions>>
1965 <<x of x0 definitions>>=
1966 typedef struct x_of_xo_data_struct {
1968 one_dim_fn_t **tension_handler; /* array of fn pointers */
1969 void **handler_data; /* array of void* pointers */
1975 double x_of_xo(double xo, void *pdata)
1977 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
1978 double F, x=0, xi, xi_guess, lb, ub;
1980 double min_relx=1e-6, min_rely=1e-6;
1982 assert(data->n_groups > 0);
1984 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
1986 fprintf(stderr, "%s: finding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
1989 if (data->xi) data->xi[0] = xo;
1990 for (i=1; i < data->n_groups; i++) {
1991 if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
1992 else xi_guess = data->xi[i];
1994 fprintf(stderr, "%s: searching for proper x[%d] for F[%d](x[%d]) = %g N\n", __FUNCTION__, i, i, i, F);
1996 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
1997 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
1999 fprintf(stderr, "%s: found F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
2003 if (data->xi) data->xi[i] = xi;
2006 fprintf(stderr, "%s: found x(x0=%g) = %g\n", __FUNCTION__, xo, x);
2012 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
2013 number of steps for monotonically increasing $f(x)$. Simple
2014 bisection, so it's robust and converges fairly quickly. You can set
2015 the maximum number of steps to take, but if you have enough processor
2016 power, it's better to limit your precision with the relative limits
2017 [[min_relx]] and [[y]]. These ensure that the width of the bracket is
2018 small on the length scale of it's larger side. Note that \emph{both}
2019 relative limits must be satisfied for the search to stop.
2020 <<one dimensional function definition>>=
2021 /* equivalent to gsl_func_t */
2022 typedef double one_dim_fn_t(double x, void *params);
2024 <<one dimensional solve declaration>>=
2025 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2026 double min_relx, double min_rely, int max_steps,
2030 <<one dimensional solve>>=
2031 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2032 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2033 double min_relx, double min_rely, int max_steps,
2036 double yx, ylb, yub, x;
2039 ylb = (*f)(lb, params);
2040 yub = (*f)(ub, params);
2042 /* check some simple cases */
2044 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bounded */
2045 else return lb; /* any x on the range [lb, ub] would work */
2047 if (ylb == y) { x=lb; yx=ylb; goto end; }
2048 if (yub == y) { x=ub; yx=yub; goto end; }
2051 fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2053 assert(ylb < y); assert(yub > y);
2054 for (i=0; i<max_steps; i++) {
2056 yx = (*f)(x, params);
2058 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);
2060 if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2061 else if (yx > y) { ub=x; yub=yx; }
2062 else /* yx < y */ { lb=x; ylb=yx; }
2067 fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2073 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2074 Generate bracketing $x$ values through bisection or exponential growth.
2075 <<one dimensional bracket declaration>>=
2076 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2079 <<one dimensional bracket>>=
2080 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2082 double yx, step, x=xguess;
2084 yx = (*f)(x, params);
2086 fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2093 fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2097 if (x == 0) x = length_scale/2.0; /* HACK */
2100 yx = (*f)(x, params);
2102 fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2107 fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2110 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2114 \subsection{Balancing implementation}
2116 <<tension-balance.c>>=
2119 <<tension balance includes>>
2120 #include "tension_balance.h"
2122 double length_scale = 1e-10; /* HACK */
2124 <<tension balance internal definitions>>
2125 <<tension balance functions>>
2128 <<tension balance includes>>=
2129 #include <assert.h> /* assert() */
2130 #include <stdlib.h> /* NULL */
2131 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2132 #include <stdio.h> /* fprintf(), stdout */
2136 \subsection{Balancing unit tests}
2138 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2139 <<check-tension-balance.c>>=
2140 <<tension balance check includes>>
2141 <<tension balance test suite>>
2142 <<main check program>>
2145 <<tension balance check includes>>=
2146 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2147 #include <assert.h> /* assert() */
2150 #include "tension_balance.h"
2153 <<tension balance test suite>>=
2154 <<tension balance function tests>>
2155 <<tension balance suite definition>>
2158 <<tension balance suite definition>>=
2159 Suite *test_suite (void)
2161 Suite *s = suite_create ("tension balance");
2162 <<tension balance function test case defs>>
2164 <<tension balance function test case adds>>
2169 <<tension balance function tests>>=
2170 <<check relative error>>
2172 double hooke(double x, void *pK)
2175 return *((double*)pK) * x;
2178 START_TEST(test_single_function)
2180 double k=5, x=3, last_x=2.0, F;
2181 one_dim_fn_t *handlers[] = {&hooke};
2182 void *data[] = {&k};
2183 F = tension_balance(1, handlers, data, last_x, x);
2184 fail_unless(F = k*x, NULL);
2189 We can also test balancing two springs with different spring constants.
2190 The analytic solution for a general number of springs is given in Appendix \ref{app.math_hooke}.
2191 <<tension balance function tests>>=
2192 START_TEST(test_double_hooke)
2194 double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
2195 one_dim_fn_t *handlers[] = {&hooke, &hooke};
2196 void *data[] = {&k1, &k2};
2197 F = tension_balance(2, handlers, data, last_x, x);
2201 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2202 //CHECK_ERR(1e-6, x1e, xi[0]);
2203 //CHECK_ERR(1e-6, x2e, xi[1]);
2204 CHECK_ERR(1e-6, Fe, F);
2208 TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
2210 <<tension balance function test case defs>>=
2211 TCase *tc_tbfunc = tcase_create("tension balance function");
2214 <<tension balance function test case adds>>=
2215 tcase_add_test(tc_tbfunc, test_single_function);
2216 tcase_add_test(tc_tbfunc, test_double_hooke);
2217 suite_add_tcase(s, tc_tbfunc);
2222 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2223 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2224 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2230 <<list definitions>>
2231 <<list declarations>>
2235 <<list declarations>>=
2236 <<head and tail declarations>>
2237 <<list length declaration>>
2238 <<push declaration>>
2240 <<list destroy declaration>>
2241 <<list to array declaration>>
2242 <<move declaration>>
2243 <<sort declaration>>
2244 <<unique declaration>>
2245 <<list print declaration>>
2249 <<create and destroy node>>
2262 <<list module makefile lines>>=
2263 NW_SAWSIM_MODS += list
2264 CHECK_BINS += check_list
2268 <<list definitions>>=
2269 typedef struct list_struct {
2270 struct list_struct *next;
2271 struct list_struct *prev;
2275 <<comparison function typedef>>
2278 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2279 <<head and tail declarations>>=
2280 list_t *head(list_t *list);
2281 list_t *tail(list_t *list);
2284 list_t *head(list_t *list)
2288 while (list->prev != NULL) {
2294 list_t *tail(list_t *list)
2298 while (list->next != NULL) {
2305 <<list length declaration>>=
2306 int list_length(list_t *list);
2309 int list_length(list_t *list)
2316 while (list->next != NULL) {
2324 [[push]] inserts a node relative to the current position in the list
2325 without changing the current position.
2326 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2327 so we set it to that of the pushed domain.
2328 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2329 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2330 <<push declaration>>=
2331 void push(list_t **pList, void *data, int below);
2334 void push(list_t **pList, void *data, int below)
2336 list_t *list, *new_node;
2337 assert(pList != NULL);
2339 new_node = create_node(data);
2342 else if (below==0) { /* insert above */
2343 if (list->prev != NULL)
2344 list->prev->next = new_node;
2345 new_node->prev = list->prev;
2346 list->prev = new_node;
2347 new_node->next = list;
2348 } else { /* insert below */
2349 if (list->next != NULL)
2350 list->next->prev = new_node;
2351 new_node->next = list->next;
2352 list->next = new_node;
2353 new_node->prev = list;
2358 [[pop]] removes the current domain node, moving the current position
2359 to the node after it, unless that node is [[NULL]], in which case move
2360 the current position to the node before it.
2361 <<pop declaration>>=
2362 void *pop(list_t **pList);
2365 void *pop(list_t **pList)
2367 list_t *list, *popped;
2369 assert(pList != NULL);
2371 assert(list != NULL); /* not an empty list */
2373 /* bypass the popped node */
2374 if (list->prev != NULL)
2375 list->prev->next = list->next;
2376 if (list->next != NULL) {
2377 list->next->prev = list->prev;
2378 *pList = list->next;
2380 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
2382 destroy_node(popped);
2387 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
2388 If the cleanup function is [[NULL]], no cleanup function is called.
2389 The nodes are not popped in any particular order.
2390 <<list destroy declaration>>=
2391 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
2394 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
2398 assert(pList != NULL);
2401 assert(list != NULL); /* not an empty list */
2402 while (list != NULL) {
2404 if (destroy != NULL)
2410 [[list_to_array]] converts a list to an array of pointers, preserving order.
2411 <<list to array declaration>>=
2412 void list_to_array(list_t *list, int *length, void ***array);
2415 void list_to_array(list_t *list, int *pLength, void ***pArray)
2419 assert(list != NULL);
2420 assert(pLength != NULL);
2421 assert(pArray != NULL);
2423 length = list_length(list);
2424 array = (void **)malloc(sizeof(void *)*length);
2425 assert(array != NULL);
2428 while (list != NULL) {
2429 array[i++] = list->d;
2437 [[move]] moves the current element along the list in either direction.
2438 <<move declaration>>=
2439 void move(list_t **pList, int downstream);
2442 void move(list_t **pList, int downstream)
2444 list_t *list, *flipper;
2446 assert(pList != NULL);
2447 list = *pList; /* a>B>c>d */
2448 assert(list != NULL); /* not an empty list */
2449 if (downstream == 0)
2450 flipper = list->prev; /* flipper is A */
2452 flipper = list->next; /* flipper is C */
2453 /* ensure there is somewhere to go in the direction we're heading */
2454 assert(flipper != NULL);
2455 /* remove the current element */
2456 data = pop(&list); /* data=B, a>C>d */
2457 /* and put it back in in the correct spot */
2459 if (downstream == 0) {
2460 push(&list, data, 0); /* b>A>c>d */
2461 list = list->prev; /* B>a>c>d */
2463 push(&list, data, 1); /* a>C>b>d */
2464 list = list->next; /* a>c>B>d */
2470 [[sort]] sorts a list from smallest to largest according to a user
2471 specified comparison.
2472 <<comparison function typedef>>=
2473 typedef int (comparison_fn_t)(void *A, void *B);
2476 <<int comparison function>>=
2477 int int_comparison(void *A, void *B)
2482 if (a > b) return 1;
2483 else if (a == b) return 0;
2487 <<sort declaration>>=
2488 void sort(list_t **pList, comparison_fn_t *comp);
2490 Here we use the bubble sort algorithm.
2492 void sort(list_t **pList, comparison_fn_t *comp)
2496 assert(pList != NULL);
2498 assert(list != NULL); /* not an empty list */
2499 while (stable == 0) {
2502 while (list->next != NULL) {
2503 if ((*comp)(list->d, list->next->d) > 0) {
2505 move(&list, 1 /* downstream */);
2514 [[unique]] removes duplicates from a sorted list according to a user
2515 specified comparison. Don't do this unless you have other pointers
2516 any dynamically allocated data in your list, because [[unique]] just
2517 drops any non-unique data without freeing it.
2518 <<unique declaration>>=
2519 void unique(list_t **pList, comparison_fn_t *comp);
2522 void unique(list_t **pList, comparison_fn_t *comp)
2525 assert(pList != NULL);
2527 assert(list != NULL); /* not an empty list */
2529 while (list->next != NULL) {
2530 if ((*comp)(list->d, list->next->d) == 0) {
2539 [[list_print]] prints a list to a [[FILE *]] stream.
2540 <<list print declaration>>=
2541 void list_print(FILE *file, list_t *list, char *list_name);
2544 void list_print(FILE *file, list_t *list, char *list_name)
2546 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
2548 while (list != NULL) {
2549 fprintf(file, " %p", list->d);
2552 fprintf(file, "\n");
2556 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
2557 <<create and destroy node>>=
2558 list_t *create_node(void *data)
2560 list_t *ret = (list_t *)malloc(sizeof(list_t));
2561 assert(ret != NULL);
2568 void destroy_node(list_t *node)
2574 The user must free the data pointed to by the node on their own.
2576 \subsection{List implementation}
2586 #include <assert.h> /* assert() */
2587 #include <stdlib.h> /* malloc(), free(), rand() */
2588 #include <stdio.h> /* fprintf(), stdout, FILE */
2589 #include "global.h" /* destroy_data_func_t */
2592 \subsection{List unit tests}
2594 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2596 <<list check includes>>
2598 <<main check program>>
2601 <<list check includes>>=
2602 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2603 #include <stdio.h> /* FILE */
2609 <<list test suite>>=
2612 <<list suite definition>>
2615 <<list suite definition>>=
2616 Suite *test_suite (void)
2618 Suite *s = suite_create ("list");
2619 <<push test case defs>>
2620 <<pop test case defs>>
2622 <<push test case adds>>
2623 <<pop test case adds>>
2629 START_TEST(test_push)
2634 push(&list, (void *)i, 1);
2635 fail_unless(list == head(list), NULL);
2636 fail_unless( (int)list->d == 0 );
2637 for (i=0; i<n; i++) {
2638 /* we expect 0, n-1, n-2, ..., 1 */
2641 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
2647 <<push test case defs>>=
2648 TCase *tc_push = tcase_create("push");
2651 <<push test case adds>>=
2652 tcase_add_test(tc_push, test_push);
2653 suite_add_tcase(s, tc_push);
2658 <<pop test case defs>>=
2660 <<pop test case adds>>=
2663 \section{Function string evaluation}
2665 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).
2666 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
2667 The solution is to run a scripting language as a subprocess accessed via pipes.
2668 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
2670 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
2671 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
2672 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.
2673 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
2675 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.
2676 Then you can either statically or dynamically link to those hardcoded functions.
2677 While much less flexible, this approach would be a fairly simple-to-implement fallback.
2679 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
2680 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
2683 #ifndef STRING_EVAL_H
2684 #define STRING_EVAL_H
2686 <<string eval setup declaration>>
2687 <<string eval function declaration>>
2688 <<string eval teardown declaration>>
2689 #endif /* STRING_EVAL_H */
2692 <<string eval module makefile lines>>=
2693 NW_SAWSIM_MODS += string_eval
2694 CHECK_BINS += check_string_eval
2695 check_string_eval_MODS =
2698 For an introduction to POSIX process control, see\\
2699 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
2700 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
2701 and of course, the relavent [[man]] pages.
2703 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
2704 [[execvp]] replaces the calling process' program with a new program.
2705 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
2706 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
2707 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
2708 The new program needs command line arguments, just like it would if you were running it from a shell.
2709 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
2710 with the final array entry being a [[NULL]] pointer.
2712 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
2713 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
2714 <<string eval subprocess definitions>>=
2715 #define SUBPROCESS "python"
2716 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
2717 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
2718 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
2720 The [[i]] option lets Python know that it should run in interactive mode.
2721 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
2722 In interactive mode, python acts on every instruction as soon as it is recieved.
2723 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.
2724 %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.
2726 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
2727 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
2728 [[fork]] returns two copies of the same program, executing the original code.
2729 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
2730 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
2732 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.
2733 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
2734 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
2735 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
2736 <<string eval pipe definitions>>=
2737 #define PIPE_READ 0 /* the end you read from */
2738 #define PIPE_WRITE 1 /* the end you write to */
2740 #define STDIN 0 /* initial index of stdin pair */
2741 #define STDOUT 2 /* initial index of stdout pair */
2744 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
2746 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.
2748 <<string eval setup declaration>>=
2749 extern void string_eval_setup(FILE **pIn, FILE **pOut);
2751 <<string eval setup definition>>=
2752 void string_eval_setup(FILE **pIn, FILE **pOut)
2755 int pfd[4]; /* file descriptors for stdin and stdout */
2757 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
2758 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
2760 pid = fork(); /* split process into two copies */
2762 if (pid == -1) { /* fork error */
2763 perror("fork error");
2765 } else if (pid == 0) { /* child */
2766 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
2767 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
2768 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
2769 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
2770 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
2771 perror("exec error"); /* exec shouldn't return */
2773 } else { /* parent */
2774 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
2775 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
2776 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
2777 if ( *pIn == NULL ) {
2778 perror("fdopen (in)");
2781 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
2782 if ( *pOut == NULL ) {
2783 perror("fdopen (out)");
2790 To use the evaluating subprocess, we just pipe in our command, and read out the results.
2791 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
2792 <<string eval function declaration>>=
2793 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
2795 <<string eval function definition>>=
2796 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
2799 rc = fprintf(in, "%s", input);
2800 assert(rc == strlen(input));
2803 alarm(1); /* set a one second timeout on the read */
2804 assert( fgets(output, buflen, out) != NULL );
2805 alarm(0); /* cancel the timeout */
2806 //fprintf(stderr, "eval: %s ----> %s", input, output);
2809 The [[alarm]] calls set and clear a timeout on the returned output.
2810 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.
2811 This protects against invalid input for which a line of output is not printed to [[stdout]].
2812 Other invalid inputs (e.g. those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
2813 If you are getting strange results, check your python code seperately. TODO, better error handling.
2815 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
2816 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
2817 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
2818 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
2819 <<string eval teardown declaration>>=
2820 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
2823 <<string eval teardown definition>>=
2824 void string_eval_teardown(FILE **pIn, FILE **pOut)
2829 /* redirect Python's stderr */
2830 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
2834 assert( fclose(*pIn) == 0 );
2836 assert( fclose(*pOut) == 0 );
2839 /* wait for python to exit */
2841 pid = wait(&stat_loc);
2848 if (WIFEXITED(stat_loc)) {
2849 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
2850 } else if (WIFSIGNALED(stat_loc)) {
2851 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
2856 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
2858 \subsection{String evaluation implementation}
2862 <<string eval includes>>
2863 #include "string_eval.h"
2864 <<string eval internal definitions>>
2865 <<string eval functions>>
2868 <<string eval includes>>=
2869 #include <assert.h> /* assert() */
2870 #include <stdlib.h> /* NULL */
2871 #include <stdio.h> /* fprintf(), stdout, fdopen() */
2872 #include <string.h> /* strlen() */
2873 #include <sys/types.h> /* pid_t */
2874 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
2875 #include <wait.h> /* wait() */
2878 <<string eval internal definitions>>=
2879 <<string eval subprocess definitions>>
2880 <<string eval pipe definitions>>
2883 <<string eval functions>>=
2884 <<string eval setup definition>>
2885 <<string eval function definition>>
2886 <<string eval teardown definition>>
2889 \subsection{String evaluation unit tests}
2891 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2892 <<check-string-eval.c>>=
2893 <<string eval check includes>>
2894 <<string comparison definition and globals>>
2895 <<string eval test suite>>
2896 <<main check program>>
2899 <<string eval check includes>>=
2900 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
2901 #include <stdio.h> /* printf() */
2902 #include <assert.h> /* assert() */
2903 #include <string.h> /* strlen() */
2904 #include <signal.h> /* SIGKILL */
2906 #include "string_eval.h"
2909 <<string eval test suite>>=
2910 <<string eval tests>>
2911 <<string eval suite definition>>
2914 <<string eval suite definition>>=
2915 Suite *test_suite (void)
2917 Suite *s = suite_create ("string eval");
2918 <<string eval test case defs>>
2920 <<string eval test case add>>
2925 <<string eval tests>>=
2926 START_TEST(test_setup_teardown)
2929 string_eval_setup(&in, &out);
2930 string_eval_teardown(&in, &out);
2933 START_TEST(test_invalid_command)
2936 char input[80], output[80]={};
2937 string_eval_setup(&in, &out);
2938 sprintf(input, "print ABCDefg\n");
2939 string_eval(in, out, input, 80, output);
2940 string_eval_teardown(&in, &out);
2943 START_TEST(test_simple_eval)
2946 char input[80], output[80]={};
2947 string_eval_setup(&in, &out);
2948 sprintf(input, "print 3+4*5\n");
2949 string_eval(in, out, input, 80, output);
2950 fail_unless(STRMATCH(output,"23\n"), NULL);
2951 string_eval_teardown(&in, &out);
2954 START_TEST(test_multiple_evals)
2957 char input[80], output[80]={};
2958 string_eval_setup(&in, &out);
2959 sprintf(input, "print 3+4*5\n");
2960 string_eval(in, out, input, 80, output);
2961 fail_unless(STRMATCH(output,"23\n"), NULL);
2962 sprintf(input, "print (3**2 + 4**2)**0.5\n");
2963 string_eval(in, out, input, 80, output);
2964 fail_unless(STRMATCH(output,"5.0\n"), NULL);
2965 string_eval_teardown(&in, &out);
2968 START_TEST(test_eval_with_state)
2971 char input[80], output[80]={};
2972 string_eval_setup(&in, &out);
2973 sprintf(input, "print 3+4*5\n");
2974 fprintf(in, "A = 3\n");
2975 sprintf(input, "print A*3\n");
2976 string_eval(in, out, input, 80, output);
2977 fail_unless(STRMATCH(output,"9\n"), NULL);
2978 string_eval_teardown(&in, &out);
2982 <<string eval test case defs>>=
2983 TCase *tc_string_eval = tcase_create("string_eval");
2985 <<string eval test case add>>=
2986 tcase_add_test(tc_string_eval, test_setup_teardown);
2987 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
2988 tcase_add_test(tc_string_eval, test_simple_eval);
2989 tcase_add_test(tc_string_eval, test_multiple_evals);
2990 tcase_add_test(tc_string_eval, test_eval_with_state);
2991 suite_add_tcase(s, tc_string_eval);
2995 \section{Accelerating function evaluation}
2997 My first version-0.3 code was running very slowly.
2998 With the options suggested in the help
2999 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
3000 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
3001 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
3002 That's only 15 calls per solution, so the search algorithm seems reasonable.
3003 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
3008 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
3010 #endif /* ACCEL_K_H */
3013 <<accel k module makefile lines>>=
3014 NW_SAWSIM_MODS += accel_k
3015 #CHECK_BINS += check_accel_k
3016 check_accel_k_MODS =
3020 #include <assert.h> /* assert() */
3021 #include <stdlib.h> /* realloc(), free(), NULL */
3022 #include "global.h" /* environment_t */
3023 #include "k_model.h" /* k_func_t */
3024 #include "interp.h" /* interp_* */
3025 #include "accel_k.h"
3027 typedef struct accel_k_struct {
3028 interp_table_t *itable;
3034 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3035 static int num_accels = 0;
3036 static accel_k_t *accels=NULL;
3038 /* Wrap k in the standard f(x) acceleration form */
3039 static double k_wrap(double F, void *params)
3041 accel_k_t *p = (accel_k_t *)params;
3042 return (*p->k)(F, p->env, p->params);
3045 static int k_tol(double FA, double kA, double FB, double kB)
3048 if (FB-FA > 1e-12) {
3049 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3050 return 1; /* unnacceptable */
3052 //printf("acceptable tol\n");
3053 return 0; /* acceptable */
3057 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3060 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3061 assert(accels != NULL);
3062 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3064 accels[i].env = env;
3065 accels[i].params = params;
3072 for (i=0; i<num_accels; i++)
3073 interp_table_free(accels[i].itable);
3075 if (accels) free(accels);
3079 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3082 for (i=0; i<num_accels; i++) {
3083 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3084 /* We've already setup this function.
3085 * WARNING: we're assuming param and env are the same. */
3086 return interp_table_eval(accels[i].itable, accels+i, F);
3089 /* set up a new acceleration environment */
3090 i = add_accel_k(k, env, params);
3091 return interp_table_eval(accels[i].itable, accels+i, F);
3095 \section{Tension models}
3096 \label{sec.tension_models}
3098 TODO: tension model intro.
3099 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.
3101 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3102 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]].
3104 <<tension-model.h>>=
3105 #ifndef TENSION_MODEL_H
3106 #define TENSION_MODEL_H
3108 <<tension handler types>>
3109 <<constant tension model declarations>>
3110 <<hooke tension model declarations>>
3111 <<worm-like chain tension model declarations>>
3112 <<freely-jointed chain tension model declarations>>
3113 <<find tension definitions>>
3114 <<tension model global declarations>>
3115 #endif /* TENSION_MODEL_H */
3118 <<tension model module makefile lines>>=
3119 NW_SAWSIM_MODS += tension_model
3120 #CHECK_BINS += check_tension_model
3121 check_tension_model_MODS =
3123 <<tension model utils makefile lines>>=
3124 TENSION_MODEL_MODS = tension_model parse list tension_balance
3125 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3126 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3127 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3128 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
3129 notangle -Rtension-model-utils.c $< > $@
3130 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
3131 gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3132 clean_tension_model_utils :
3133 rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3134 @ The pipe symbol [[|]] marks the prerequisits following it (in this
3135 case, the directories) as ‘order-only’ prerequisites. The timestamp
3136 on these prerequisits does not effect whether the rules are executed.
3137 This is appropriate for directories, where we don't need to recompile
3138 after an unrelated has been added to the directory, but only when the
3139 source prerequisites change. See the [[make]] documentation for more
3141 (\url{http://www.gnu.org/software/make/manual/html_node/Prerequisite-Types.html}).
3144 \label{sec.null_tension}
3146 For unstretchable domains.
3148 <<null tension model getopt>>=
3149 {NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3152 \subsection{Constant}
3153 \label{sec.const_tension}
3155 <<constant tension functions>>=
3156 <<constant tension function>>
3157 <<constant tension parameter create/destroy functions>>
3160 <<constant tension model declarations>>=
3161 <<constant tension function declaration>>
3162 <<constant tension parameter create/destroy function declarations>>
3163 <<constant tension model global declarations>>
3167 An infinitely stretchable domain providing a constant tension.
3169 <<constant tension function declaration>>=
3170 extern double const_tension_handler(double x, void *pdata);
3172 <<constant tension function>>=
3173 double const_tension_handler(double x, void *pdata)
3175 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3176 list_t *list = data->group;
3181 assert(list != NULL); /* empty active group?! */
3182 F = ((const_tension_param_t *)list->d)->F;
3184 fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3186 while (list != NULL) {
3187 assert(((const_tension_param_t *)list->d)->F == F);
3195 <<constant tension structure>>=
3196 typedef struct const_tension_param_struct {
3197 double F; /* tension (force) in N */
3198 } const_tension_param_t;
3202 <<constant tension parameter create/destroy function declarations>>=
3203 extern void *string_create_const_tension_param_t(char **param_strings);
3204 extern void destroy_const_tension_param_t(void *p);
3207 <<constant tension parameter create/destroy functions>>=
3208 const_tension_param_t *create_const_tension_param_t(double F)
3210 const_tension_param_t *ret
3211 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3212 assert(ret != NULL);
3217 void *string_create_const_tension_param_t(char **param_strings)
3219 return create_const_tension_param_t(atof(param_strings[0]));
3222 void destroy_const_tension_param_t(void *p)
3230 <<constant tension model global declarations>>=
3231 extern int num_const_tension_params;
3232 extern char *const_tension_param_descriptions[];
3233 extern char const_tension_param_string[];
3236 <<constant tension model globals>>=
3237 int num_const_tension_params = 1;
3238 char *const_tension_param_descriptions[] = {"tension F, N"};
3239 char const_tension_param_string[] = "0";
3242 <<constant tension model getopt>>=
3243 {&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}
3249 <<hooke functions>>=
3251 <<hooke parameter create/destroy functions>>
3254 <<hooke tension model declarations>>=
3255 <<hooke tension function declaration>>
3256 <<hooke tension parameter create/destroy function declarations>>
3257 <<hooke tension model global declarations>>
3261 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
3262 The behavior of a series of springs $k_i$ in series is given by
3264 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
3266 For a simple proof, see Appendix \ref{app.math_hooke}.
3268 <<hooke tension function declaration>>=
3269 extern double hooke_handler(double x, void *pdata);
3273 double hooke_handler(double x, void *pdata)
3275 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3276 list_t *list = data->group;
3281 assert(list != NULL); /* empty active group?! */
3282 while (list != NULL) {
3283 assert( ((hooke_param_t *)list->d)->k > 0 );
3284 k += 1.0/ ((hooke_param_t *)list->d)->k;
3286 fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
3287 __FUNCTION__, ((hooke_param_t *)list->d)->k);
3293 fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N\n",
3294 __FUNCTION__, k, x, k*x);
3300 <<hooke structure>>=
3301 typedef struct hooke_param_struct {
3302 double k; /* spring constant in N/m */
3306 <<hooke tension parameter create/destroy function declarations>>=
3307 extern void *string_create_hooke_param_t(char **param_strings);
3308 extern void destroy_hooke_param_t(void *p);
3311 <<hooke parameter create/destroy functions>>=
3312 hooke_param_t *create_hooke_param_t(double k)
3314 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
3315 assert(ret != NULL);
3320 void *string_create_hooke_param_t(char **param_strings)
3322 return create_hooke_param_t(atof(param_strings[0]));
3325 void destroy_hooke_param_t(void *p)
3332 <<hooke tension model global declarations>>=
3333 extern int num_hooke_params;
3334 extern char *hooke_param_descriptions[];
3335 extern char hooke_param_string[];
3338 <<hooke tension model globals>>=
3339 int num_hooke_params = 1;
3340 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
3341 char hooke_param_string[]="0.05";
3344 <<hooke tension model getopt>>=
3345 {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}
3348 \subsection{Worm-like chain}
3351 We can model several unfolded domains with a single worm-like chain.
3352 <<worm-like chain functions>>=
3353 <<worm-like chain function>>
3354 <<worm-like chain handler>>
3355 <<worm-like chain create/destroy functions>>
3358 <<worm-like chain tension model declarations>>=
3359 <<worm-like chain handler declaration>>
3360 <<worm-like chain create/destroy function declarations>>
3361 <<worm-like chain tension model global declarations>>
3365 The combination of all unfolded domains is modeled as a worm like chain,
3366 with the force $F_{\text{WLC}}$ approximately given by
3368 F_{\text{WLC}} \approx \frac{k_B T}{p}
3369 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
3371 where $x$ is the distance between the fixed ends,
3372 $k_B$ is Boltzmann's constant,
3373 $T$ is the absolute temperature,
3374 $p$ is the persistence length, and
3375 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
3376 <<worm-like chain function>>=
3377 double wlc(double x, double T, double p, double L)
3379 double xL=0.0; /* uses global k_B */
3381 assert(T > 0); assert(p > 0); assert(L > 0);
3382 if (x >= L) return HUGE_VAL;
3383 xL = x/L; /* unitless */
3384 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
3385 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
3386 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
3389 This model assumes that each unfolded domain has the same persistence length.
3391 <<worm-like chain handler declaration>>=
3392 extern double wlc_handler(double x, void *pdata);
3395 <<worm-like chain handler>>=
3396 double wlc_handler(double x, void *pdata)
3398 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3399 list_t *list = data->group;
3403 assert(list != NULL); /* empty active group?! */
3404 p = ((wlc_param_t *)list->d)->p;
3405 while (list != NULL) {
3406 /* enforce consistent p */
3407 assert( ((wlc_param_t *)list->d)->p == p);
3408 L += ((wlc_param_t *)list->d)->L;
3410 fprintf(stderr, "%s: WLC section %g m long in series\n",
3411 __FUNCTION__, ((wlc_param_t *)list->d)->L);
3416 fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
3417 __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
3419 return wlc(x, data->env->T, p, L);
3423 <<worm-like chain structure>>=
3424 typedef struct wlc_param_struct {
3425 double p; /* WLC persistence length (m) */
3426 double L; /* WLC contour length (m) */
3430 <<worm-like chain create/destroy function declarations>>=
3431 extern void *string_create_wlc_param_t(char **param_strings);
3432 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
3435 <<worm-like chain create/destroy functions>>=
3436 wlc_param_t *create_wlc_param_t(double p, double L)
3438 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
3439 assert(ret != NULL);
3445 void *string_create_wlc_param_t(char **param_strings)
3447 return create_wlc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3450 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
3457 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.
3458 TODO (now we copy the string before parsing, but still do this for clarity).
3460 <<valid string write code>>=
3461 char string[] = "abc";
3464 <<valid string write code>>=
3465 char *string = "abc";
3469 <<worm-like chain tension model global declarations>>=
3470 extern int num_wlc_params;
3471 extern char *wlc_param_descriptions[];
3472 extern char wlc_param_string[];
3474 <<worm-like chain tension model globals>>=
3475 int num_wlc_params = 2;
3476 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
3477 char wlc_param_string[]="0.39e-9,28.4e-9";
3480 <<worm-like chain tension model getopt>>=
3481 {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}
3483 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
3485 \subsection{Freely-jointed chain}
3488 <<freely-jointed chain functions>>=
3489 <<freely-jointed chain function>>
3490 <<freely-jointed chain handler>>
3491 <<freely-jointed chain create/destroy functions>>
3494 <<freely-jointed chain tension model declarations>>=
3495 <<freely-jointed chain handler declaration>>
3496 <<freely-jointed chain create/destroy function declarations>>
3497 <<freely-jointed chain tension model global declarations>>
3501 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
3502 The tension of a chain of $N$ such links is given by
3504 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
3506 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}.
3507 <<freely-jointed chain function>>=
3508 double langevin(double x, void *params)
3511 return 1.0/tanh(x) - 1/x;
3513 <<one dimensional bracket declaration>>
3514 <<one dimensional solve declaration>>
3515 double inv_langevin(double x)
3517 double lb=0.0, ub=1.0, ret;
3518 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
3519 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
3523 double fjc(double x, double T, double l, int N)
3525 double L = l*(double)N;
3526 if (x == 0) return 0;
3527 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
3528 return k_B*T/l * inv_langevin(x/L);
3532 <<freely-jointed chain handler declaration>>=
3533 extern double fjc_handler(double x, void *pdata);
3535 <<freely-jointed chain handler>>=
3536 double fjc_handler(double x, void *pdata)
3538 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3539 list_t *list = data->group;
3544 assert(list != NULL); /* empty active group?! */
3545 l = ((fjc_param_t *)list->d)->l;
3546 while (list != NULL) {
3547 /* enforce consistent l */
3548 assert( ((fjc_param_t *)list->d)->l == l);
3549 N += ((fjc_param_t *)list->d)->N;
3551 fprintf(stderr, "%s: FJC section %d links long in series\n",
3552 __FUNCTION__, ((fjc_param_t *)list->d)->N);
3557 fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
3558 __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
3560 return fjc(x, data->env->T, l, N);
3564 <<freely-jointed chain structure>>=
3565 typedef struct fjc_param_struct {
3566 double l; /* FJC link length (m) */
3567 int N; /* FJC number of links */
3571 <<freely-jointed chain create/destroy function declarations>>=
3572 extern void *string_create_fjc_param_t(char **param_strings);
3573 extern void destroy_fjc_param_t(void *p);
3576 <<freely-jointed chain create/destroy functions>>=
3577 fjc_param_t *create_fjc_param_t(double l, double N)
3579 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
3580 assert(ret != NULL);
3586 void *string_create_fjc_param_t(char **param_strings)
3588 return create_fjc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3591 void destroy_fjc_param_t(void *p)
3598 <<freely-jointed chain tension model global declarations>>=
3599 extern int num_fjc_params;
3600 extern char *fjc_param_descriptions[];
3601 extern char fjc_param_string[];
3604 <<freely-jointed chain tension model globals>>=
3605 int num_fjc_params = 2;
3606 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
3607 char fjc_param_string[]="0.5e-9,200";
3610 <<freely-jointed chain tension model getopt>>=
3611 {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}
3614 \subsection{Tension model implementation}
3616 <<tension-model.c>>=
3619 <<tension model includes>>
3620 #include "tension_model.h"
3621 <<tension model internal definitions>>
3622 <<tension model globals>>
3623 <<tension model functions>>
3626 <<tension model includes>>=
3627 #include <assert.h> /* assert() */
3628 #include <stdlib.h> /* NULL */
3629 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
3630 #include <stdio.h> /* fprintf(), stdout */
3633 #include "tension_balance.h" /* oneD_*() */
3636 <<tension model internal definitions>>=
3637 <<constant tension structure>>
3639 <<worm-like chain structure>>
3640 <<freely-jointed chain structure>>
3643 <<tension model globals>>=
3644 <<tension handler array global>>
3645 <<constant tension model globals>>
3646 <<hooke tension model globals>>
3647 <<worm-like chain tension model globals>>
3648 <<freely-jointed chain tension model globals>>
3651 <<tension model functions>>=
3652 <<constant tension functions>>
3654 <<worm-like chain functions>>
3655 <<freely-jointed chain functions>>
3659 \subsection{Utilities}
3661 The tension models can get complicated, and users may want to reassure
3662 themselves that this portion of the input physics are functioning properly. The
3663 stand-alone program [[tension_model_utils]] demonstrates the tension model
3664 interface without the overhead of having to understand a full simulation.
3665 [[tension_model_utils]] takes command line model arguments like the full
3666 simulation, and prints $F(x)$ data to the screen for the selected model over a
3669 <<tension-model-utils.c>>=
3671 <<tension model utility includes>>
3672 <<tension model utility definitions>>
3673 <<tension model utility getopt functions>>
3675 <<tension model utility main>>
3678 <<tension model utility main>>=
3679 int main(int argc, char **argv)
3681 <<tension model getopt array definition>>
3682 tension_model_getopt_t *model = NULL;
3686 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
3687 tension_handler_data_t tdata;
3688 int num_param_args; /* for INIT_MODEL() */
3689 char **param_args; /* for INIT_MODEL() */
3690 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model, &flags);
3692 if (flags & VFLAG) {
3693 printf("#initializing model %s with parameters %s\n", model->name, model->params);
3695 INIT_MODEL("utility", model, model->params, params);
3697 push(&tdata.group, params, 1);
3699 tdata.persist = NULL;
3700 if (model->handler == NULL) {
3701 printf("No tension function for model %s\n", model->name);
3705 double dx=1e-10, x=0, F=0;
3706 printf("#F (N)\tk (%% pop. per s)\n");
3707 while (F >= 0 && F < 1e5 && x < 1e-6) {
3708 F = (*model->handler)(x, &tdata);
3709 printf("%g\t%g\n", x, F);
3713 params = pop(&tdata.group);
3714 if (model->destructor)
3715 (*model->destructor)(params);
3720 <<tension model utility includes>>=
3723 #include <string.h> /* strlen() */
3724 #include <assert.h> /* assert() */
3728 #include "tension_model.h"
3731 <<tension model utility definitions>>=
3732 <<version definition>>
3733 #define VFLAG 1 /* verbose */
3734 <<string comparison definition and globals>>
3735 <<tension model getopt definitions>>
3736 <<initialize model definition>>
3740 <<tension model utility getopt functions>>=
3741 <<index tension model>>
3742 <<tension model utility help>>
3743 <<tension model utility get options>>
3746 <<tension model utility help>>=
3747 void help(char *prog_name,
3749 int n_tension_models, tension_model_getopt_t *tension_models,
3753 printf("usage: %s [options]\n", prog_name);
3754 printf("Version %s\n\n", VERSION);
3755 printf("A utility for understanding the available tension models\n\n");
3756 printf("Environmental options:\n");
3757 printf("-T T\tTemperature T (currently %g K)\n", env->T);
3758 printf("-C T\tYou can also set the temperature T in Celsius\n");
3759 printf("Model options:\n");
3760 printf("-m model\tFolded tension model (currently %s)\n",
3761 tension_models[tension_model].name);
3762 printf("-a args\tFolded tension model argument string (currently %s)\n",
3763 tension_models[tension_model].params);
3764 printf("F(x) is calculated for a range of x and printed\n");
3765 printf("For example:\n");
3766 printf(" #Distance (x)\tForce (N)\n");
3767 printf(" 123.456\t7.89\n");
3769 printf("-V\tChange output to verbose mode\n");
3770 printf("-h\tPrint this help and exit\n");
3772 printf("Tension models:\n");
3773 for (i=0; i<n_tension_models; i++) {
3774 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
3775 for (j=0; j < tension_models[i].num_params ; j++)
3776 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
3777 printf("\t\tdefaults: %s\n", tension_models[i].params);
3782 <<tension model utility get options>>=
3783 void get_options(int argc, char **argv, environment_t *env,
3784 int n_tension_models, tension_model_getopt_t *tension_models,
3785 tension_model_getopt_t **model,
3786 unsigned int *flags)
3788 char *prog_name = NULL;
3789 char c, options[] = "T:C:m:a:Vh";
3790 int tension_model=0, num_strings;
3791 extern char *optarg;
3792 extern int optind, optopt, opterr;
3796 /* setup defaults */
3798 prog_name = argv[0];
3799 env->T = 300.0; /* K */
3801 *model = tension_models;
3803 while ((c=getopt(argc, argv, options)) != -1) {
3805 case 'T': env->T = atof(optarg); break;
3806 case 'C': env->T = atof(optarg)+273.15; break;
3808 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
3809 *model = tension_models+tension_model;
3812 tension_models[tension_model].params = optarg;
3814 case 'V': *flags |= VFLAG; break;
3816 fprintf(stderr, "unrecognized option '%c'\n", optopt);
3817 /* fall through to default case */
3819 help(prog_name, env, n_tension_models, tension_models, tension_model);
3828 \section{Unfolding rate models}
3829 \label{sec.k_models}
3831 We have two main choices for determining $k$: Bell's model, which assumes the
3832 folded domains exist in a single `folded' state and make instantaneous,
3833 stochastic transitions over a free energy barrier; and Kramers' model, which
3834 treats the unfolding as a thermalized diffusion process.
3835 We deal with these options like we dealt with the different tension models: by bundling all model-specific
3836 parameters into structures, with model specific functions for determining $k$.
3838 <<k func definition>>=
3839 typedef double k_func_t(double F, environment_t *env, void *params);
3842 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3843 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]].
3845 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
3846 because \LaTeX\ doesn't like underscores outside math mode.
3847 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
3848 You could use spaces instead (`\verb+ +'),
3849 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
3850 which I don't like as much.
3856 <<k func definition>>
3857 <<null k declarations>>
3858 <<const k declarations>>
3859 <<bell k declarations>>
3860 <<kramers k declarations>>
3861 <<kramers integ k declarations>>
3862 #endif /* K_MODEL_H */
3865 <<k model module makefile lines>>=
3866 NW_SAWSIM_MODS += k_model
3867 CHECK_BINS += check_k_model
3868 check_k_model_MODS = parse string_eval
3870 [[check_k_model]] also depends on the parse module.
3872 <<k model utils makefile lines>>=
3873 K_MODEL_MODS = k_model parse string_eval
3874 K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
3875 $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3876 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
3877 notangle -Rk-model-utils.c $< > $@
3878 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
3879 gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3880 clean_k_model_utils :
3881 rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
3887 For domains that are never folded.
3889 <<null k declarations>>=
3891 <<null k functions>>=
3896 <<null k model getopt>>=
3897 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
3900 \subsection{Constant rate model}
3903 This is a pretty boring model, but it is usefull for testing the engine.
3907 where $k_0$ is the constant unfolding rate.
3909 <<const k functions>>=
3910 <<const k function>>
3911 <<const k structure create/destroy functions>>
3914 <<const k declarations>>=
3915 <<const k function declaration>>
3916 <<const k structure create/destroy function declarations>>
3917 <<const k global declarations>>
3921 <<const k structure>>=
3922 typedef struct const_k_param_struct {
3923 double knot; /* transition rate k_0 (frac population per s) */
3927 <<const k function declaration>>=
3928 double const_k(double F, environment_t *env, void *const_k_params);
3930 <<const k function>>=
3931 double const_k(double F, environment_t *env, void *const_k_params)
3932 { /* Returns the rate constant k in frac pop per s. */
3933 const_k_param_t *p = (const_k_param_t *) const_k_params;
3935 assert(p->knot > 0);
3940 <<const k structure create/destroy function declarations>>=
3941 void *string_create_const_k_param_t(char **param_strings);
3942 void destroy_const_k_param_t(void *p);
3945 <<const k structure create/destroy functions>>=
3946 const_k_param_t *create_const_k_param_t(double knot)
3948 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
3949 assert(ret != NULL);
3954 void *string_create_const_k_param_t(char **param_strings)
3956 return create_const_k_param_t(atof(param_strings[0]));
3959 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
3966 <<const k global declarations>>=
3967 extern int num_const_k_params;
3968 extern char *const_k_param_descriptions[];
3969 extern char const_k_param_string[];
3972 <<const k globals>>=
3973 int num_const_k_params = 1;
3974 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
3975 char const_k_param_string[]="1";
3978 <<const k model getopt>>=
3979 {"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}
3982 \subsection{Bell's model}
3985 Quantitatively, Bell's model gives $k$ as
3987 k = k_0 \cdot e^{F dx / k_B T} \;,
3989 where $k_0$ is the unforced unfolding rate,
3990 $F$ is the applied unfolding force,
3991 $dx$ is the distance to the transition state, and
3992 $k_B T$ is the thermal energy\citep{schlierf06}.
3994 <<bell k functions>>=
3996 <<bell k structure create/destroy functions>>
3999 <<bell k declarations>>=
4000 <<bell k function declaration>>
4001 <<bell k structure create/destroy function declarations>>
4002 <<bell k global declarations>>
4006 <<bell k structure>>=
4007 typedef struct bell_param_struct {
4008 double knot; /* transition rate k_0 (frac population per s) */
4009 double dx; /* `distance to transition state' (nm) */
4013 <<bell k function declaration>>=
4014 double bell_k(double F, environment_t *env, void *bell_params);
4016 <<bell k function>>=
4017 double bell_k(double F, environment_t *env, void *bell_params)
4018 { /* Returns the rate constant k in frac pop per s.
4019 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4020 * uses global k_B in J/K */
4021 bell_param_t *p = (bell_param_t *) bell_params;
4022 assert(F >= 0); assert(env->T > 0);
4024 assert(p->knot > 0); assert(p->dx >= 0);
4026 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
4027 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
4028 p->knot * exp(F*p->dx / (k_B*env->T) ));
4029 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
4031 return p->knot * exp(F*p->dx / (k_B*env->T) );
4035 <<bell k structure create/destroy function declarations>>=
4036 void *string_create_bell_param_t(char **param_strings);
4037 void destroy_bell_param_t(void *p);
4040 <<bell k structure create/destroy functions>>=
4041 bell_param_t *create_bell_param_t(double knot, double dx)
4043 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
4044 assert(ret != NULL);
4050 void *string_create_bell_param_t(char **param_strings)
4052 return create_bell_param_t(atof(param_strings[0]), atof(param_strings[1]));
4055 void destroy_bell_param_t(void *p /* bell_param_t * */)
4062 <<bell k global declarations>>=
4063 extern int num_bell_params;
4064 extern char *bell_param_descriptions[];
4065 extern char bell_param_string[];
4069 int num_bell_params = 2;
4070 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4071 char bell_param_string[]="3.3e-4,0.25e-9";
4074 <<bell k model getopt>>=
4075 {"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}
4077 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4078 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4080 \subsection{Kramer's model}
4083 Kramer's model gives $k$ as
4085 % \frac{1}{k} = \frac{1}{D}
4086 % \int_{x_\text{min}}^{x_\text{max}}
4087 % e^{\frac{-U_F(x)}{k_B T}}
4088 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4091 %where $D$ is the diffusion constant,
4092 %$U_F(x)$ is the free energy along the unfolding coordinate, and
4093 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4094 % before the minimum and after the maximum, respectively \citep{schlierf06}.
4096 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
4098 where $D$ is the diffusion constant,
4100 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
4102 are length scales where
4103 $x_c(F)$ is the location of the bound state and
4104 $x_{ts}(F)$ is the location of the transition state,
4105 $E(F,x)$ is the free energy, and
4106 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
4107 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
4108 \citet{evans97} Eqn.~3).
4110 <<kramers k functions>>=
4111 <<kramers k function>>
4112 <<kramers k structure create/destroy functions>>
4115 <<kramers k declarations>>=
4116 <<kramers k function declaration>>
4117 <<kramers k structure create/destroy function declarations>>
4118 <<kramers k global declarations>>
4122 <<kramers k structure>>=
4123 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
4124 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
4126 typedef struct kramers_param_struct {
4127 double D; /* diffusion rate (um/s) */
4134 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
4135 //kramers_x_func_t *xts; /* function returning transition x (nm) */
4136 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
4137 //kramers_E_func_t *E; /* function returning E (J) */
4138 //double *E_params; /* parametrize the energy landscape */
4139 //int n_E_params; /* length of E_params array */
4143 <<kramers k function declaration>>=
4144 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
4145 extern double kramers_k(double F, environment_t *env, void *kramers_params);
4147 <<kramers k function>>=
4148 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
4150 static char input[160]={0};
4151 static char output[80]={0};
4152 /* setup the environment */
4153 fprintf(in, "F = %g; T = %g\n", F, T);
4154 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4155 string_eval(in, out, input, 80, output);
4156 return atof(output);
4159 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
4161 static char input[160]={0};
4162 static char output[80]={0};
4163 /* setup the environment */
4164 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
4165 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4166 string_eval(in, out, input, 80, output);
4167 return atof(output);
4170 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
4172 kramers_param_t *p = (kramers_param_t *) kramers_params;
4173 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
4176 double kramers_k(double F, environment_t *env, void *kramers_params)
4177 { /* Returns the rate constant k in frac pop per s.
4178 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4179 * uses global k_B in J/K */
4180 kramers_param_t *p = (kramers_param_t *) kramers_params;
4181 double kbT, xc, xts, lc, lts, Eb;
4182 assert(F >= 0); assert(env->T > 0);
4185 assert(p->xc != NULL); assert(p->xts != NULL);
4186 assert(p->ddEdxx != NULL); assert(p->E != NULL);
4187 assert(p->in != NULL); assert(p->out != NULL);
4189 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
4190 if (xc == -1.0) return -1.0; /* error (Too much force) */
4191 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
4192 if (xc == -1.0) return -1.0; /* error (Too much force) */
4193 lc = sqrt(2.0*M_PI*kbT /
4194 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
4195 lts = sqrt(-2.0*M_PI*kbT/
4196 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
4197 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
4198 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
4199 //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));
4200 return p->D/(lc*lts) * exp(-Eb/kbT);
4204 <<kramers k structure create/destroy function declarations>>=
4205 void *string_create_kramers_param_t(char **param_strings);
4206 void destroy_kramers_param_t(void *p);
4209 <<kramers k structure create/destroy functions>>=
4210 kramers_param_t *create_kramers_param_t(double D,
4211 char *xc_fn, char *xts_fn,
4212 char *E_fn, char *ddEdxx_fn,
4213 char *global_define_string)
4214 // kramers_x_func_t *xc, kramers_x_func_t *xts,
4215 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
4216 // double *E_params, int n_E_params)
4218 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
4219 assert(ret != NULL);
4220 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4221 assert(ret->xc != NULL);
4222 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4223 assert(ret->xts != NULL);
4224 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
4225 assert(ret->E != NULL);
4226 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
4227 assert(ret->ddEdxx != NULL);
4229 strcpy(ret->xc, xc_fn);
4230 strcpy(ret->xts, xts_fn);
4231 strcpy(ret->E, E_fn);
4232 strcpy(ret->ddEdxx, ddEdxx_fn);
4233 //ret->E_params = E_params;
4234 //ret->n_E_params = n_E_params;
4235 string_eval_setup(&ret->in, &ret->out);
4236 fprintf(ret->in, "kB = %g\n", k_B);
4237 fprintf(ret->in, "%s\n", global_define_string);
4241 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
4242 void *string_create_kramers_param_t(char **param_strings)
4244 return create_kramers_param_t(atof(param_strings[0]),
4252 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
4254 kramers_param_t *pk = (kramers_param_t *)p;
4256 string_eval_teardown(&pk->in, &pk->out);
4258 // free(pk->E_params);
4264 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
4265 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.
4266 We follow \citet{shillcock98} and use
4268 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
4269 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
4272 where TODO, variable meanings.%\citep{shillcock98}.
4273 The first and second derivatives of $E(F,E_b,x,x_s)$ are
4275 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\\
4276 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
4278 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
4279 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
4280 The roots are therefor at
4282 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
4283 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
4284 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
4287 <<kramers k global declarations>>=
4288 extern int num_kramers_params;
4289 extern char *kramers_param_descriptions[];
4290 extern char kramers_param_string[];
4293 <<kramers k globals>>=
4294 int num_kramers_params = 6;
4295 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"};
4296 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)}";
4298 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
4299 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
4300 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
4301 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.
4302 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
4303 It works so long as [[val_1]] is not [[false]].
4305 <<kramers k model getopt>>=
4306 {"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}
4309 \subsection{Kramer's model (natural cubic splines)}
4310 \label{sec.kramers_integ}
4312 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
4313 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
4314 \citet{schlierf06} Eqn.~2)
4316 \frac{1}{k} = \frac{1}{D}
4317 \int_{x_\text{min}}^{x_\text{max}}
4318 e^{\frac{U_F(x)}{k_B T}}
4319 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4322 where $D$ is the diffusion constant,
4323 $U_F(x)$ is the free energy along the unfolding coordinate, and
4324 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4325 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
4327 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
4328 interpolating between them with cubic splines.
4329 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
4331 <<kramers integ k functions>>=
4332 <<spline functions>>
4333 <<kramers integ A k functions>>
4334 <<kramers integ B k functions>>
4335 <<kramers integ k function>>
4336 <<kramers integ k structure create/destroy functions>>
4339 <<kramers integ k declarations>>=
4340 <<kramers integ k function declaration>>
4341 <<kramers integ k structure create/destroy function declarations>>
4342 <<kramers integ k global declarations>>
4346 <<kramers integ k structure>>=
4347 typedef double func_t(double x, void *params);
4348 typedef struct kramers_integ_param_struct {
4349 double D; /* diffusion rate (um/s) */
4350 double x_min; /* integration bounds */
4352 func_t *E; /* function returning E (J) */
4353 void *E_params; /* parametrize the energy landscape */
4354 destroy_data_func_t *destroy_E_params;
4356 double F; /* for passing into GSL evaluations */
4358 } kramers_integ_param_t;
4361 <<spline param structure>>=
4362 typedef struct spline_param_struct {
4363 int n_knots; /* natural cubic spline knots for U(x) */
4366 gsl_spline *spline; /* GSL spline parameters */
4367 gsl_interp_accel *acc; /* GSL spline acceleration data */
4371 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
4373 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4375 <<kramers integ A k functions>>=
4376 double kramers_integ_k_integrandA(double x, void *params)
4378 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4379 double U = (*p->E)(x, p->E_params);
4380 double Fx = p->F * x;
4381 double kbT = k_B * p->env->T;
4382 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
4383 return exp(-(U-Fx)/kbT);
4385 double kramers_integ_k_integralA(double x, void *params)
4387 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4389 double result, abserr;
4390 size_t neval; /* for qng */
4391 static gsl_integration_workspace *w = NULL;
4392 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4393 f.function = &kramers_integ_k_integrandA;
4395 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4396 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4397 //fprintf(stderr, "integralA = %g\n", result);
4403 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
4404 \int_{x_\text{min}}^{x_\text{max}}
4405 e^{\frac{U_F(x)}{k_B T}}
4406 \text{Integral}_A(x)
4409 <<kramers integ B k functions>>=
4410 double kramers_integ_k_integrandB(double x, void *params)
4412 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4413 double U = (*p->E)(x, p->E_params);
4414 double Fx = p->F * x;
4415 double kbT = k_B * p->env->T;
4416 double intA = kramers_integ_k_integralA(x, params);
4417 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
4418 return exp((U-Fx)/kbT)*intA;
4420 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
4422 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4424 double result, abserr;
4425 size_t neval; /* for qng */
4426 static gsl_integration_workspace *w = NULL;
4427 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4428 f.function = &kramers_integ_k_integrandB;
4432 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4433 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4434 //fprintf(stderr, "integralB = %g\n", result);
4439 With the integrals taken care of, evaluating $k$ is simply
4441 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
4443 <<kramers integ k function declaration>>=
4444 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
4446 <<kramers integ k function>>=
4447 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
4448 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
4449 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
4450 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
4454 <<kramers integ k structure create/destroy function declarations>>=
4455 void *string_create_kramers_integ_param_t(char **param_strings);
4456 void destroy_kramers_integ_param_t(void *p);
4459 <<kramers integ k structure create/destroy functions>>=
4460 kramers_integ_param_t *create_kramers_integ_param_t(double D,
4461 double xmin, double xmax, func_t *E,
4463 destroy_data_func_t *destroy_E_params)
4465 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
4466 assert(ret != NULL);
4471 ret->E_params = E_params;
4472 ret->destroy_E_params = destroy_E_params;
4476 void *string_create_kramers_integ_param_t(char **param_strings)
4478 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
4479 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
4480 void *E_params = string_create_spline_param_t(param_strings+1);
4481 return create_kramers_integ_param_t(atof(param_strings[0]),
4482 atof(param_strings[2]),
4483 atof(param_strings[3]),
4484 &spline_eval, E_params,
4485 destroy_spline_param_t);
4488 void destroy_kramers_integ_param_t(void *params)
4490 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4493 (*p->destroy_E_params)(p->E_params);
4499 Finally we have the GSL spline wrappers:
4500 <<spline functions>>=
4502 <<create/destroy spline>>
4505 <<spline function>>=
4506 double spline_eval(double x, void *spline_params)
4508 spline_param_t *p = (spline_param_t *)(spline_params);
4509 return gsl_spline_eval(p->spline, x, p->acc);
4513 <<create/destroy spline>>=
4514 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
4516 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
4517 assert(ret != NULL);
4518 ret->n_knots = n_knots;
4521 ret->acc = gsl_interp_accel_alloc();
4522 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
4523 gsl_spline_init(ret->spline, x, y, n_knots);
4527 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
4528 void *string_create_spline_param_t(char **param_strings)
4532 double *x=NULL, *y=NULL;
4533 /* break into ordered pairs */
4534 parse_list_string(param_strings[0], SEP, '(', ')',
4535 &num_knots, &knot_coords);
4536 x = (double *)malloc(sizeof(double)*num_knots);
4538 y = (double *)malloc(sizeof(double)*num_knots);
4540 for (i=0; i<num_knots; i++) {
4541 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
4542 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
4545 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
4547 free_parsed_list(num_knots, knot_coords);
4548 return create_spline_param_t(num_knots, x, y);
4551 void destroy_spline_param_t(void *params)
4553 spline_param_t *p = (spline_param_t *)params;
4556 gsl_spline_free(p->spline);
4558 gsl_interp_accel_free(p->acc);
4568 <<kramers integ k global declarations>>=
4569 extern int num_kramers_integ_params;
4570 extern char *kramers_integ_param_descriptions[];
4571 extern char kramers_integ_param_string[];
4574 <<kramers integ k globals>>=
4575 int num_kramers_integ_params = 4;
4576 char *kramers_integ_param_descriptions[] = {
4577 "diffusion D, m^2 / s",
4578 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
4579 "starting integration bound x_min, m",
4580 "ending integration bound x_max, m"};
4581 //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";
4582 //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";
4583 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";
4584 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
4585 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
4586 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
4589 <<kramers integ k model getopt>>=
4590 {"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}
4592 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4593 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4595 \subsection{Unfolding model implementation}
4599 <<k model includes>>
4600 #include "k_model.h"
4601 <<k model internal definitions>>
4603 <<k model functions>>
4606 <<k model includes>>=
4607 #include <assert.h> /* assert() */
4608 #include <stdlib.h> /* NULL, malloc() */
4609 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4610 #include <stdio.h> /* fprintf(), stdout */
4611 #include <string.h> /* strlen(), strcpy() */
4612 #include <gsl/gsl_integration.h>
4613 #include <gsl/gsl_spline.h>
4618 <<k model internal definitions>>=
4619 <<const k structure>>
4620 <<bell k structure>>
4621 <<kramers k structure>>
4622 <<kramers integ k structure>>
4623 <<spline param structure>>
4626 <<k model globals>>=
4630 <<kramers k globals>>
4631 <<kramers integ k globals>>
4634 <<k model functions>>=
4635 <<null k functions>>
4636 <<const k functions>>
4637 <<bell k functions>>
4638 <<kramers k functions>>
4639 <<kramers integ k functions>>
4642 \subsection{Unfolding model unit tests}
4644 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
4645 <<check-k-model.c>>=
4646 <<k model check includes>>
4647 <<check relative error>>
4649 <<k model test suite>>
4650 <<main check program>>
4653 <<k model check includes>>=
4654 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
4655 #include <stdio.h> /* sprintf() */
4656 #include <assert.h> /* assert() */
4657 #include <math.h> /* exp() */
4660 #include "k_model.h"
4663 <<k model test suite>>=
4666 <<k model suite definition>>
4669 <<k model suite definition>>=
4670 Suite *test_suite (void)
4672 Suite *s = suite_create ("k model");
4673 <<const k test case defs>>
4674 <<bell k test case defs>>
4676 <<const k test case adds>>
4677 <<bell k test case adds>>
4682 \subsubsection{Constant}
4684 <<const k test case defs>>=
4685 TCase *tc_const_k = tcase_create("Constant k");
4688 <<const k test case adds>>=
4689 tcase_add_test(tc_const_k, test_const_k_create_destroy);
4690 tcase_add_test(tc_const_k, test_const_k_over_range);
4691 suite_add_tcase(s, tc_const_k);
4696 START_TEST(test_const_k_create_destroy)
4698 char *knot[] = {"1","2","3","4","5","6"};
4699 char *params[] = {knot[0]};
4702 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4703 params[0] = knot[i];
4704 p = string_create_const_k_param_t(params);
4705 destroy_const_k_param_t(p);
4710 START_TEST(test_const_k_over_range)
4712 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
4713 char *knot[] = {"1","2","3","4","5","6"};
4714 char *params[] = {knot[0]};
4717 char param_string[80];
4720 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4721 params[0] = knot[i];
4722 p = string_create_const_k_param_t(params);
4723 for ( F=Fm; F<FM; F+=dF ) {
4724 fail_unless(const_k(F, &env, p)==atof(knot[i]),
4725 "const_k(%g, %g, &{%s}) = %g != %s",
4726 F, T, knot[i], const_k(F, &env, p), knot[i]);
4728 destroy_const_k_param_t(p);
4734 \subsubsection{Bell}
4736 <<bell k test case defs>>=
4737 TCase *tc_bell = tcase_create("Bell's k");
4740 <<bell k test case adds>>=
4741 tcase_add_test(tc_bell, test_bell_k_create_destroy);
4742 tcase_add_test(tc_bell, test_bell_k_at_zero);
4743 tcase_add_test(tc_bell, test_bell_k_at_one);
4744 suite_add_tcase(s, tc_bell);
4748 START_TEST(test_bell_k_create_destroy)
4750 char *knot[] = {"1","2","3","4","5","6"};
4752 char *params[] = {knot[0], dx};
4755 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4756 params[0] = knot[i];
4757 p = string_create_bell_param_t(params);
4758 destroy_bell_param_t(p);
4763 START_TEST(test_bell_k_at_zero)
4765 double F=0.0, T=300.0;
4766 char *knot[] = {"1","2","3","4","5","6"};
4768 char *params[] = {knot[0], dx};
4771 char param_string[80];
4774 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4775 params[0] = knot[i];
4776 p = string_create_bell_param_t(params);
4777 fail_unless(bell_k(F, &env, p)==atof(knot[i]),
4778 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
4779 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
4780 destroy_bell_param_t(p);
4785 START_TEST(test_bell_k_at_one)
4788 char *knot[] = {"1","2","3","4","5","6"};
4790 char *params[] = {knot[0], dx};
4791 double F= k_B*T/atof(dx);
4794 char param_string[80];
4797 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4798 params[0] = knot[i];
4799 p = string_create_bell_param_t(params);
4800 CHECK_ERR(1e-6, atof(knot[i])*exp(1.0), bell_k(F, &env, p));
4801 //fail_unless(bell_k(F, &env, p)==atof(knot[i])*exp(1.0),
4802 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
4803 // F, T, knot[i], dx, bell_k(F, &env, p), atof(knot[i])*exp(1.0));
4804 destroy_bell_param_t(p);
4810 <<kramers k tests>>=
4813 <<kramers k test case def>>=
4816 <<kramers k test case add>>=
4819 <<k model function tests>>=
4820 <<check relative error>>
4822 START_TEST(test_something)
4824 double k=5, x=3, last_x=2.0, F;
4825 one_dim_fn_t *handlers[] = {&hooke};
4826 void *data[] = {&k};
4827 double xi[] = {0.0};
4829 int new_active_groups = 1;
4830 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
4831 fail_unless(F = k*x, NULL);
4836 \subsection{Utilities}
4838 The unfolding models can get complicated, and are sometimes hard to explain to others.
4839 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
4840 the overhead of having to understand a full simulation.
4841 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
4842 or special mode, where the operation depends on the specific model selected.
4844 <<k-model-utils.c>>=
4846 <<k model utility includes>>
4847 <<k model utility definitions>>
4848 <<k model utility getopt functions>>
4849 <<k model utility multi model E>>
4850 <<k model utility main>>
4853 <<k model utility main>>=
4854 int main(int argc, char **argv)
4856 <<k model getopt array definition>>
4857 k_model_getopt_t *model = NULL;
4862 int num_param_args; /* for INIT_MODEL() */
4863 char **param_args; /* for INIT_MODEL() */
4865 double special_xmin;
4866 double special_xmax;
4867 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
4868 &Fmax, &special_xmin, &special_xmax, &flags);
4869 if (flags & VFLAG) {
4870 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4872 INIT_MODEL("utility", model, model->params, params);
4875 if (model->k == NULL) {
4876 printf("No k function for model %s\n", model->name);
4880 printf("Fmax = %g <= 0\n", Fmax);
4886 printf("#F (N)\tk (%% pop. per s)\n");
4887 for (i=0; i<=N; i++) {
4888 F = Fmax*i/(double)N;
4889 k = (*model->k)(F, &env, params);
4891 printf("%g\t%g\n", F, k);
4896 if (model->k == NULL || model->k == &bell_k) {
4897 printf("No special mode for model %s\n", model->name);
4900 if (special_xmax <= special_xmin) {
4901 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
4905 double Xrng=(special_xmax-special_xmin), x, E;
4907 printf("#x (m)\tE (J)\n");
4908 for (i=0; i<=N; i++) {
4909 x = special_xmin + Xrng*i/(double)N;
4910 E = multi_model_E(model->k, params, &env, x);
4911 printf("%g\t%g\n", x, E);
4918 if (model->destructor)
4919 (*model->destructor)(params);
4924 <<k model utility multi model E>>=
4925 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
4927 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
4929 if (k == kramers_integ_k)
4930 E = (*p->E)(x, p->E_params);
4932 E = kramers_E(0, env, model_params, x);
4938 <<k model utility includes>>=
4941 #include <string.h> /* strlen() */
4942 #include <assert.h> /* assert() */
4945 #include "string_eval.h"
4946 #include "k_model.h"
4949 <<k model utility definitions>>=
4950 <<version definition>>
4951 #define VFLAG 1 /* verbose */
4952 enum mode_t {M_K_OF_F, M_SPECIAL};
4953 <<string comparison definition and globals>>
4954 <<k model getopt definitions>>
4955 <<initialize model definition>>
4956 <<kramers k structure>>
4957 <<kramers integ k structure>>
4961 <<k model utility getopt functions>>=
4963 <<k model utility help>>
4964 <<k model utility get options>>
4967 <<k model utility help>>=
4968 void help(char *prog_name,
4970 int n_k_models, k_model_getopt_t *k_models,
4974 printf("usage: %s [options]\n", prog_name);
4975 printf("Version %s\n\n", VERSION);
4976 printf("A utility for understanding the available unfolding models\n\n");
4977 printf("Environmental options:\n");
4978 printf("-T T\tTemperature T (currently %g K)\n", env->T);
4979 printf("-C T\tYou can also set the temperature T in Celsius\n");
4980 printf("Model options:\n");
4981 printf("-k model\tTransition rate model (currently %s)\n",
4982 k_models[k_model].name);
4983 printf("-K args\tTransition rate model argument string (currently %s)\n",
4984 k_models[k_model].params);
4985 printf("There are two output modes. In standard mode, k(F) is printed\n");
4986 printf("For example:\n");
4987 printf(" #Force (N)\tk (% pop. per s)\n");
4988 printf(" 123.456\t7.89\n");
4990 printf("In special mode, the output depends on the model.\n");
4991 printf("For models defining an energy function E(x), that function is printed\n");
4992 printf(" #Position (m)\tE (J)\n");
4993 printf(" 0.0012\t0.0034\n");
4995 printf("-m\tChange output to standard mode\n");
4996 printf("-M\tChange output to special mode\n");
4997 printf("-F\tSet the maximum F value for the standard mode k(F) output\n");
4998 printf("-x\tSet the minimum x value for the special mode E(x) output\n");
4999 printf("-X\tSet the maximum x value for the special mode E(x) output\n");
5000 printf("-V\tChange output to verbose mode\n");
5001 printf("-h\tPrint this help and exit\n");
5003 printf("Unfolding rate models:\n");
5004 for (i=0; i<n_k_models; i++) {
5005 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
5006 for (j=0; j < k_models[i].num_params ; j++)
5007 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
5008 printf("\t\tdefaults: %s\n", k_models[i].params);
5013 <<k model utility get options>>=
5014 void get_options(int argc, char **argv, environment_t *env,
5015 int n_k_models, k_model_getopt_t *k_models,
5016 enum mode_t *mode, k_model_getopt_t **model,
5017 double *Fmax, double *special_xmin, double *special_xmax,
5018 unsigned int *flags)
5020 char *prog_name = NULL;
5021 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
5023 extern char *optarg;
5024 extern int optind, optopt, opterr;
5028 /* setup defaults */
5030 prog_name = argv[0];
5031 env->T = 300.0; /* K */
5036 *special_xmax = 1e-8;
5037 *special_xmin = 0.0;
5039 while ((c=getopt(argc, argv, options)) != -1) {
5041 case 'T': env->T = atof(optarg); break;
5042 case 'C': env->T = atof(optarg)+273.15; break;
5044 k_model = index_k_model(n_k_models, k_models, optarg);
5045 *model = k_models+k_model;
5048 k_models[k_model].params = optarg;
5050 case 'm': *mode = M_K_OF_F; break;
5051 case 'M': *mode = M_SPECIAL; break;
5052 case 'F': *Fmax = atof(optarg); break;
5053 case 'x': *special_xmin = atof(optarg); break;
5054 case 'X': *special_xmax = atof(optarg); break;
5055 case 'V': *flags |= VFLAG; break;
5057 fprintf(stderr, "unrecognized option '%c'\n", optopt);
5058 /* fall through to default case */
5060 help(prog_name, env, n_k_models, k_models, k_model);
5069 \section{\LaTeX\ documentation}
5071 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
5072 The comment blocks are extracted (with nicely formatted code blocks), using
5073 <<latex makefile lines>>=
5074 $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5075 noweave -latex -index -delay $< > $@
5076 $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR)
5080 <<latex makefile lines>>=
5081 $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
5083 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5084 $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
5085 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5086 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5087 mv $(BUILD_DIR)/sawsim.pdf $@
5089 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
5090 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
5091 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
5093 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
5094 <<latex makefile lines>>=
5096 rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
5097 $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
5098 $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.tex \
5099 $(BUILD_DIR)/sawsim.bib
5100 clean_latex : semi-clean_latex
5101 rm -f $(DOC_DIR)/sawsim.pdf
5106 [[make]] is a common utility on *nix systems for managing dependencies while building software.
5107 In this case, we have several \emph{targets} that we'd like to build:
5108 the main simulation program \prog;
5109 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
5110 the documentation [[sawsim.pdf]];
5111 and the [[Makefile]] itself.
5112 Besides the generated files, there is the utility target
5113 [[clean]] for removing all generated files except [[Makefile]].
5115 % [[dist]] for generating a distributable tar file.
5117 Extract the makefile with
5118 `[[notangle -Rmakefile src/sawsim.nw | sed 's/ /\t/' > Makefile]]'.
5119 The sed is needed because notangle replaces tabs with eight spaces,
5120 but [[make]] needs tabs.
5122 # Decide what directories to put things in
5127 # Note: Cannot use, for example, `./build', because make eats the `./'
5128 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5130 # Modules (X.c, X.h) defined in the noweb file
5133 # Modules defined outside the noweb file
5134 FREE_SAWSIM_MODS = interp tavl
5136 <<list module makefile lines>>
5137 <<tension balance module makefile lines>>
5138 <<tension model module makefile lines>>
5139 <<k model module makefile lines>>
5140 <<parse module makefile lines>>
5141 <<string eval module makefile lines>>
5142 <<accel k module makefile lines>>
5144 SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS)
5146 # Everything needed for sawsim under one roof. sawsim.c must come first
5147 SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \
5148 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h)
5149 # Libraries needed to compile sawsim
5150 LIBS = gsl gslcblas m
5151 CHECK_LIBS = gsl gslcblas m check
5152 # The non-check binaries generated
5153 BINS = sawsim tension_model_utils k_model_utils sawsim_profile
5156 # Define the major targets
5157 all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
5159 view : $(DOC_DIR)/sawsim.pdf
5161 profile : $(BIN_DIR)/sawsim_profile | $(BIN_DIR)
5162 $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \
5163 -mnull -Mwlc -A0.39e-9,28e-9 -F8
5164 gprof $(BIN_DIR)/sawsim_profile gmon.out > $@
5165 check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim
5166 $(SHELL) -e -c 'for B in $^; do ./$$B; done'
5167 clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \
5168 clean_tension_model_utils \
5169 clean_k_model_utils clean_latex clean_check_sawsim
5170 rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \
5171 $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \
5172 $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
5173 $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
5174 $(BUILD_DIR)/global.h ./gmon.out
5175 $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
5177 # Various builds of sawsim
5178 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
5179 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5180 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
5181 gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5182 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
5183 gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5185 # Create the directories
5186 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
5189 # Copy over the source external to sawsim
5190 # Note: Cannot use, for example, `./build', because make eats the `./'
5191 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5193 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
5197 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
5201 ## Basic source generated with noweb
5202 # The central files sawsim.c and global.h...
5203 $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5205 $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5206 notangle -Rglobal.h $< > $@
5207 # ... and the modules
5208 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5209 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5210 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5211 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5212 $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
5213 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5214 # Note: I use `_' as a space in my file names, but noweb doesn't like
5215 # that so we use `-' to name the noweb roots and substitute here.
5217 ## Building the unit-test programs
5219 $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
5220 notangle -Rchecks $< > $@
5221 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
5222 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
5223 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR)
5224 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
5225 clean_check_sawsim :
5226 rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
5227 # ... and the modules
5229 $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
5230 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \
5231 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \
5232 $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
5233 $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
5234 $$(BUILD_DIR)/global.h | $(BIN_DIR)
5235 gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
5236 $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
5238 # todo: clean up the required modules to
5239 $(CHECK_BINS:%=clean_%) :
5240 rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
5242 # Cleaning up the modules
5244 $(SAWSIM_MODS:%=clean_%) :
5245 rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h)
5247 <<tension model utils makefile lines>>
5248 <<k model utils makefile lines>>
5249 <<latex makefile lines>>
5251 Makefile : $(SRC_DIR)/sawsim.nw
5252 notangle -Rmakefile $< | sed 's/ /\t/' > $@
5254 This is fairly self-explanatory. We compile a dynamically linked
5255 version ([[sawsim]]) and a statically linked version
5256 ([[sawsim_static]]). The static version is about 9 times as big, but
5257 it runs without needing dynamic GSL libraries to link to, so it's a
5258 better format for distributing to a cluster for parallel evaluation.
5262 \subsection{Hookean springs in series}
5263 \label{app.math_hooke}
5266 The effective spring constant for $n$ springs $k_i$ in series is given by
5268 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5274 F &= k_i x_i &&\forall i \le n \\
5275 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
5276 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
5277 F &= k_1 x_1 = k_\text{eff} x \\
5278 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
5279 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5284 \addcontentsline{toc}{section}{References}
5285 \bibliography{sawsim}