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 %\bibliographystyle{ieeetr} % pick the bibliography style, short and sweet
24 %\bibliographystyle{plain} % pick the bibliography style, includes dates
25 \bibliographystyle{plainnat}
26 \defcitealias{sw:noweb}{{\tt noweb}}
27 \defcitealias{galassi05}{GSL}
28 \defcitealias{sw:check}{{\tt check}}
29 \defcitealias{sw:python}{Python}
34 \textheight 9.5in % leave a bit of extra room for page numbers (at top of page)
37 \setlength{\parindent}{0pt}
38 \setlength{\parskip}{5pt}
40 % For some reason, the \maketitle wants to go on it's own page
41 % so we'll just hardcode the following in our first page.
42 %\title{Sawsim: a sawtooth protein unfolding simulator}
43 %\author{W.~Trevor King}
46 \newcommand{\prog}{[[sawsim]]}
47 \newcommand{\lang}{[[C]]}
48 \newcommand{\p}[3]{\left#1 #2 \right#3} % e.g. \p({complicated expr})
49 \newcommand{\Th}{\ensuremath{\sup{\text{th}}}} % ordinal th
50 \newcommand{\U}[1]{\ensuremath{\text{ #1}}} % units: $1\U{m/s}$
52 % single spaced lists, from Alvin Alexander
53 % http://www.devdaily.com/blog/post/latex/control-line-spacing-in-itemize-enumerate-tags/
54 \newenvironment{packed_item}{
56 \setlength{\itemsep}{1pt}
57 \setlength{\parskip}{0pt}
58 \setlength{\parsep}{0pt}
62 \nwfilename{sawsim.nw}
67 Sawsim: a sawtooth protein unfolding simulator \\
72 %%%%%%%%%%%%%%%%%%%%%%%%%%%% End LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%%
74 \section{Introduction}
76 The unfolding of multi-globular protein strings is a stochastic
77 process. In the \prog\ program, we use Monte Carlo methods to
78 simulate this unfolding through an explicit time-stepping approach.
79 We develop a program in \lang\ to simulate probability of unfolding
80 under a constant extension velocity of a chain of globular domains.
82 \subsection{Related work}
86 \subsection{About this document}
88 This document was written in \citetalias{sw:noweb} with code blocks in
89 \lang\ and comment blocks in \LaTeX.
91 \subsection{Change Log}
93 Version 0.0 used only the unfolded domain WLC to determine the tension
94 and had a fixed timestep $dt$.
96 Version 0.1 added adaptive timesteps, adjusting so that the unfolding
97 probability for a given domain was constant.
99 Version 0.2 added Kramers' model unfolding rate calculations, and a
100 switch to select Bell's or Kramers' model. This induced a major
101 rewrite, introducing the tension group abstraction, and a split to
102 multiple \lang\ source files. Also added a unit-test suites for
103 sanity checking, and switched to SI units throughout.
105 Version 0.3 added integrated Kramers' model unfolding rate
106 calculations. Also replaced some of my hand-coded routines with GNU
107 Scientific Library \citepalias{galassi05} calls.
109 Version 0.4 added Python evaluation for the saddle-point Kramers'
110 model unfolding rate calculations, so that the model functions could
111 be controlled from the command line. Also added interpolating lookup
112 tables to accelerate slow unfolding rate model evaluation, and command
113 line control of tension groups.
115 <<version definition>>=
116 #define VERSION "0.4"
122 sawsim - program for simulating single molecule mechanical unfolding.
123 Copyright (C) 2008, William Trevor King
125 This program is free software; you can redistribute it and/or
126 modify it under the terms of the GNU General Public License as
127 published by the Free Software Foundation; either version 3 of the
128 License, or (at your option) any later version.
130 This program is distributed in the hope that it will be useful, but
131 WITHOUT ANY WARRANTY; without even the implied warranty of
132 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
133 See the GNU General Public License for more details.
135 You should have received a copy of the GNU General Public License
136 along with this program; if not, write to the Free Software
137 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
140 The author may be contacted at <wking@drexel.edu> on the Internet, or
141 write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
142 Philadelphia PA 19104, USA.
154 The unfolding system is stored as a chain of \emph{domains}. Domains
155 can be folded, globular protein domains, unfolded protein linkers, AFM
156 cantilevers, or other stretchable system components.
158 Each domain contributes to the total tension, which depends on the
159 chain extension. The particular model for the tension as a function
160 of extension is handled generally with function pointers. So far the
161 following models have been implemented:
163 \item Null (Appendix \ref{sec.null_tension}),
164 \item Constant (Appendix \ref{sec.const_tension}),
165 \item Hookean spring (Appendix \ref{sec.hooke}),
166 \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
167 \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
170 The tension model and parameters used for a particular domain depend
171 on whether the domain is folded or unfolded. The transition rate from
172 the folded state to the unfolded state is also handled generally with
173 function pointers, with implementations for
175 \item Null (Appendix \ref{sec.null_k}),
176 \item Bell's model (Appendix \ref{sec.bell}),
177 \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
178 \item Kramers' saddle point model (Appendix \ref{sec.kramers}).
181 The unfolding of the chain domains is modeled in two stages. First
182 the tension of the chain is calculated. Then the tension (assumed to
183 be constant over the length of the chain), is applied to each folded
184 domain, and used to calculate the probability of that subdomain
185 unfolding in the next timestep $dt$. Then the time is stepped
186 forward, and the process is repeated.
189 int main(int argc, char **argv)
191 <<tension model getopt array definition>>
192 <<k model getopt array definition>>
193 list_t *domain_list=NULL; /* variables initialized in get_options() */
194 environment_t env={0};
196 int num_folded=0, unfolding;
197 double x=0, dt, F; /* variables used in the simulation loop */
198 get_options(argc, argv, &P, &dt_max, &v, &env, NUM_TENSION_MODELS,
199 tension_models, NUM_K_MODELS, k_models, &domain_list, &num_folded);
201 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Position (m)\tForce (N)\n");
202 if (flags & FLAG_OUTPUT_UNFOLDING_FORCES) printf("#Unfolding Force (N)\n");
203 while (num_folded > 0) {
204 F = find_tension(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, unfolding);
205 dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, v);
206 unfolding=random_unfoldings(domain_list, F, dt, &env, &num_folded);
207 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\n", x, F);
210 destroy_domain_list(domain_list);
214 @ The meat of the simulation is bundled into the three functions
215 [[find_tension]], [[determine_dt]], and [[random_unfoldings]].
216 [[find_tension]] is discussed in Section \ref{sec.find_tension},
217 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
218 [[random_unfoldings]] in Section \ref{sec.unfolding_rate}.
220 Environmental parameters important in determining reaction rates and
221 tensions (e.g. temperature, pH) are stored in a single structure to
222 facilitate extension to more complicated models in the future.
223 <<environment definition>>=
224 typedef struct environment_struct {
230 <<environment definition>>
231 <<one dimensional function definition>>
232 <<create/destroy definitions>>
236 \section{Simulation functions}
238 <<simulation functions>>=
242 <<does domain unfold>>
243 <<randomly unfold domains>>
247 \label{sec.find_tension}
249 Because the stretched system may be made up of several parts (folded
250 domains, unfolded domains, spring-like cantilever, \ldots), we will
251 assign the domains to models and groups. The models are used to
252 determine a domain's tension (Hookean spring, WLC, \ldots). Groups
253 will represent collections of domains which the model should treat as
254 a single entity. A domain may belong to different groups or models
255 depending on its state. For example, a domain might be modeled as a
256 freely-jointed chain when it iss folded and as a worm-like chain class
257 when it is unfolded. The domains are assumed to be commutative, so
258 ordering is ignored. The interactions between the groups are assumed
259 to be linear, but the interactions between domains of the same group
260 need not be. This allows for non-linear group models such as th
261 worm-like or freely-jointed chains. Each model has a tension handler
262 function, which gives the tension $F$ needed to stretch a given group
263 of domains a distance $x$.
265 To understand the purpose of groups, consider a chain of proteins
266 where two folded proteins $A$ and $B$ are modeled as WLCs with
267 persistence length $p_f$ and two unfolded proteins $C$ and $D$ are
268 modeled as WLCs with persistence length $p_u$. The proteins are
269 attached to a cantilever $E$ qof spring constant $κ$. The string
270 would get split into three lists:
272 \begin{tabular}{llll}
273 Model & Group & List & Tension \\
274 WLC & 0 & $[A,B]$ & $F_\text{WLC}(x, p_f, L_A+L_B)$ \\
275 WLC & 1 & $[C,D]$ & $F_\text{WLC}(x, p_u, L_C+L_D)$ \\
276 Hooke & 0 & $[E]$ & $F_\text{Hooke}(x, κ)$ \\
279 Note that group numbers only matter \emph{within} model classes, since
280 grouping domains with seperate models doesn't make sense.
282 <<find tension definitions>>=
283 #define NUM_TENSION_MODELS 5
287 <<tension handler array global declaration>>=
288 extern one_dim_fn_t *tension_handlers[];
291 <<tension handler array global>>=
292 one_dim_fn_t *tension_handlers[] = {
294 &const_tension_handler,
302 <<tension model global declarations>>=
303 <<tension handler array global declaration>>
306 <<tension handler types>>=
307 typedef struct tension_handler_data_struct {
308 /* int num_domains; */
309 /* domain_t *domains; */
313 } tension_handler_data_t;
317 After sorting the chain into separate groups $G_i$, with tension
318 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
319 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
320 \forall i,j$. Note that there may be several groups within each model
321 class. [[find_tension]] is basically a wrapper to feed proper input
322 into the [[tension_balance]] function. [[unfolding]] should be set to
323 0 if there were no unfoldings since the previous call to
326 double find_tension(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list, environment_t *env, double x, int unfolding)
328 static int num_active_groups=0;
329 static one_dim_fn_t **per_group_handlers = NULL;
330 static void **per_group_params = NULL;
331 static double last_x;
332 tension_handler_data_t *tdata;
337 fprintf(stderr, "finding tension at distance %g%s\n", x, unfolding ? " after an unfolding" : "");
338 list_print(stderr, domain_list, "domain list");
341 if (unfolding != 0 || num_active_groups == 0) { /* setup statics */
342 /* free old statics */
343 if (per_group_handlers != NULL)
344 free(per_group_handlers);
345 if (per_group_params != NULL) {
346 for (i=0; i < num_active_groups; i++) {
347 assert(per_group_params[i] != NULL);
348 free(per_group_params[i]);
350 free(per_group_params);
353 /* get new statics */
354 get_active_groups(domain_list, num_tension_handlers, tension_handlers, &num_active_groups, &per_group_handlers, &per_group_params);
356 /* fill in the group handlers and params */
357 for (i=0; i<num_active_groups; i++) {
358 tdata = (tension_handler_data_t *) per_group_params[i];
360 fprintf(stderr, "active group %d of %d", i, num_active_groups);
361 list_print(stderr, tdata->group, " parameters");
364 /* tdata->persist continues unchanged */
367 } /* else roll with the current statics */
369 F = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_params, last_x, x);
373 @ For the moment, we will restrict our group boundaries to a single
374 dimension, so $\sum_i x_i = x$, our total extension (it is this
375 restriction that causes the groups to interact linearly). We'll also
376 restrict our extensions to all be positive. With these restrictions,
377 the problem of balancing the tensions reduces to solving a system of
378 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
379 the number of active groups. In general this can be fairly
380 complicated, but without much loss of practicality we can restrict
381 ourselves to strictly monotonically increasing, non-negative tension
382 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
383 simpler. For example, it guarantees the existence of a unique, real
384 solution for finite forces. See Appendix \ref{app.tension_balance}
385 for the tension balancing implementation.
387 Each group must define a parameter structure if the tension model is
389 a function to create the parameter structure,
390 a function to destroy the parameter structure, and
391 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as it's data argument.
392 The tension-specific parameter structures are contained in the domain
393 group field. See the Section \ref{app.model_selection} for a
394 discussion on the command line framework. See the worm-like chain
395 implementation for an example (Section \ref{sec.wlc}).
397 The major limitation of our approach is that all of our tension models
398 are Markovian (memory-less), which is not really a problem for the
399 chains (see \citet{evans99} for WLC thermalization rate calculations).
401 \subsection{Unfolding rate}
402 \label{sec.unfolding_rate}
404 Each folded domain is modeled as unimolecular, first order reaction
406 \text{Folded} \xrightarrow{k} % in tex: X atop Y
409 With the rate constant $k$ defined by
411 \frac{d[\text{Folded}]}{dt} = -k [\text{Folded}].
413 So the probability of a given protein unfolding in a short time $dt$
419 <<does domain unfold>>=
420 int domain_unfolds(double F, double dt, environment_t *env, domain_t *domain)
421 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to a folded domain */
423 k = accel_k(domain->k, F, env, domain->k_params);
424 //(*domain->k)(F, env, domain->k_params);
425 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
426 return happens(k*dt); /* dice roll for prob. k*dt event */
428 @ [[happens]] is a random decision making function defined in Appendix \ref{app.utils}.
430 Only one domain can unfold in each timestep, because the timescale of
431 a domain unfolding $dt_u$ is assumed to be less than the simulation
432 timestep $dt$, so a domain will completely unfold in a single
433 timestep. We adapt our timesteps to keep the probability of a single
434 domain unfolding low, so the probability of two domains unfolding in
435 the same timestep is negligible. This approach breaks down as the
436 adaptive timestep scheme approaches $dt \sim dt_u$, but $dt_u \sim
437 1\U{$\mu$s}$ for Ig27-like proteins \citep{klimov00}, so this
438 shouldn't be a problem. To reassure yourself, you can ask the
439 simulator to print the smallest timestep that was used with TODO.
440 <<randomly unfold domains>>=
441 int random_unfoldings(list_t *domain_list, double tension,
442 double dt, environment_t *env,
443 int *num_folded_domains)
444 { /* returns 1 if there was an unfolding and 0 otherwise */
445 while (domain_list != NULL) {
446 if (D(domain_list)->state == DS_FOLDED
447 && domain_unfolds(tension, dt, env, D(domain_list))) {
448 if (flags & FLAG_OUTPUT_UNFOLDING_FORCES)
449 fprintf(stdout, "%g\n", tension);
450 D(domain_list)->state = DS_UNFOLDED;
451 (*num_folded_domains)--;
452 return 1; /* our one domain has unfolded, stop looking for others */
454 domain_list = domain_list->next;
460 \subsection{Adaptive timesteps}
461 \label{sec.adaptive_dt}
463 We'd like to pick $dt$ so the probability of unfolding in the next
464 timestep is small. If we don't adapt our timestep, we risk breaking
465 our approximation that $F(x' \in [x, x+v \cdot dt]) = F(x)$ or that
466 only one domain unfolds in a given timestep. Because $F(x)$ is
467 monotonically increasing, excessively large timesteps will lead to
468 erroneously large unfolding forces. The simulation would have been
469 accurate for sufficiently small fixed timesteps, but adaptive
470 timesteps will allow us to move through low tension regions in fewer
471 steps, leading to a more efficient simulation.
473 The actual adaptive timestep implementation is not particularly
474 interesting, since we are only required to reduce $dt$ to somewhere
475 below a set threshold, so I've removed it to Appendix
476 \ref{app.adaptive_dt}.
482 The models provide the physics of the simulation, but the simulation
483 allows interchangeable models, and we are currently focusing on the
484 simulation itself, so we remove the actual model implementation to the
487 The tension models are defined in Appendix \ref{sec.tension_models}
488 and unfolding models are defined in Appendix \ref{sec.k_models}.
491 #define k_B 1.3806503e-23 /* J/K */
495 \section{Command line interface}
497 <<get option functions>>=
499 <<index tension model>>
505 \subsection{Model selection}
506 \label{app.model_selection}
508 The main difficulty with the command line interface in \prog\ is
509 developing an intuitive method for accessing the possibly large number
510 of available models. We'll treat this generally by defining an array
511 of available models, containing their important parameters.
513 <<get options definitions>>=
514 <<model getopt definitions>>
517 <<create/destroy definitions>>=
518 typedef void *create_data_func_t(char **param_strings);
519 typedef void destroy_data_func_t(void *param_struct);
522 <<model getopt definitions>>=
523 <<tension model getopt definitions>>
524 <<k model getopt definitions>>
528 \subsubsection{Tension}
530 To access the various tension models from the command line, we define
531 a structure that contains all the neccessary information about the
533 <<tension model getopt definitions>>=
534 typedef struct tension_model_getopt_struct {
535 one_dim_fn_t *handler; /* fn implementing the model on a group */
536 char *name; /* name string identifying the model */
537 char *description; /* a brief description of the model */
538 int num_params; /* number of extra params for the model */
539 char **param_descriptions; /* descriptions of extra params */
540 char *params; /* default values for the extra params */
541 create_data_func_t *creator; /* param-string -> model-data-struct fn */
542 destroy_data_func_t *destructor; /* fn to free the model data structure */
543 } tension_model_getopt_t;
546 <<tension model getopt array definition>>=
547 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
548 <<null tension model getopt>>,
549 <<constant tension model getopt>>,
550 <<hooke tension model getopt>>,
551 <<worm-like chain tension model getopt>>,
552 <<freely-jointed chain tension model getopt>>
556 \subsubsection{Unfolding rate}
558 <<k model getopt definitions>>=
559 #define NUM_K_MODELS 5
561 typedef struct k_model_getopt_struct {
566 char **param_descriptions;
568 create_data_func_t *creator;
569 destroy_data_func_t *destructor;
573 <<k model getopt array definition>>=
574 k_model_getopt_t k_models[NUM_K_MODELS] = {
575 <<null k model getopt>>,
576 <<const k model getopt>>,
577 <<bell k model getopt>>,
578 <<kramers k model getopt>>,
579 <<kramers integ k model getopt>>
586 void help(char *prog_name, double P, double dt_max, double v,
588 int n_tension_models, tension_model_getopt_t *tension_models,
589 int folded_tension_model, int unfolded_tension_model,
590 int n_k_models, k_model_getopt_t *k_models,
594 printf("usage: %s [options]\n", prog_name);
595 printf("Version %s\n\n", VERSION);
596 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
597 printf("Simulation options:\n");
598 printf("-P P\tTarget probability for dt (currently %g)\n", P);
599 printf("-t dt\tMaximum allowed timestep dt (currently %g)\n", dt_max);
600 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
601 printf("Environmental options:\n");
602 printf("-T T\tTemperature T (currently %g K)\n", env->T);
603 printf("-C T\tYou can also set the temperature T in Celsius\n");
604 printf("Model options:\n");
605 printf("The domains exist in either the folded or unfolded state\n");
606 printf("The following options change the default behavior in each state.\n");
607 printf("-m model[,group]\tFolded tension model (currently %s)\n",
608 tension_models[folded_tension_model].name);
609 printf("-a args\tFolded tension model argument string (currently %s)\n",
610 tension_models[folded_tension_model].params);
611 printf("-M model[,group]\tUnfolded tension model (currently %s)\n",
612 tension_models[unfolded_tension_model].name);
613 printf("-A args\tUnfolded tension model argument string (currently %s)\n",
614 tension_models[unfolded_tension_model].params);
615 printf("The following options change the unfolding rate.\n");
616 printf("-k model\tTransition rate model (currently %s)\n",
617 k_models[k_model].name);
618 printf("-K args\tTransition rate model argument string (currently %s)\n",
619 k_models[k_model].params);
620 printf("Domain creation options:\n");
621 printf("Once you've set up the appropriate default models, you need to add the domains\n");
622 printf("-F n\tAdd n folded domains with the current models\n");
623 printf("-U n\tAdd n unfolded domains with the current models\n");
624 printf("Output mode options:\n");
625 printf("There are two output modes. In standard mode, only the unfolding\n");
626 printf("events are printed. For example:\n");
627 printf(" #Unfolding Force (N)\n");
628 printf(" 123.456e-12\n");
630 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
631 printf(" #Position (m)\tForce (N)\n");
632 printf(" 0.001\t0.0023\n");
634 printf("-V\tChange output to verbose mode\n");
635 printf("-h\tPrint this help and exit\n");
637 printf("Tension models:\n");
638 for (i=0; i<n_tension_models; i++) {
639 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
640 for (j=0; j < tension_models[i].num_params ; j++)
641 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
642 printf("\t\tdefaults: %s\n", tension_models[i].params);
644 printf("Unfolding rate models:\n");
645 for (i=0; i<n_k_models; i++) {
646 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
647 for (j=0; j < k_models[i].num_params ; j++)
648 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
649 printf("\t\tdefaults: %s\n", k_models[i].params);
651 printf("Examples:\n");
652 printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded\n");
653 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);
654 printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded\n");
655 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);
659 \subsection{Get options}
662 void get_options(int argc, char **argv,
663 double *pP, double *pDt_max, double *pV,
665 int n_tension_models, tension_model_getopt_t *tension_models,
666 int n_k_models, k_model_getopt_t *k_models,
667 list_t **pDomain_list, int *pNum_folded)
669 char *prog_name = NULL;
670 char c, options[] = "P:t:v:T:C:m:a:M:A:k:K:F:U:Vh";
671 int ftension_model=0, utension_model=0, k_model=0;
672 char *ftension_params=NULL, *utension_params=NULL;
673 int i, n, ftension_group=0, utension_group=0;
677 extern int optind, optopt, opterr;
682 for (i=0; i<n_tension_models; i++) {
683 fprintf(stderr, "tension model %i %s:\t model's handler %p\n",
684 i, tension_models[i].name, tension_models[i].handler);
685 assert(tension_models[i].handler == tension_handlers[i]);
690 flags = FLAG_OUTPUT_UNFOLDING_FORCES;
692 *pP = 1e-3; /* % pop per s */
693 *pDt_max = 0.001; /* s */
694 *pV = 1e-6; /* m/s */
695 env->T = 300.0; /* K */
697 while ((c=getopt(argc, argv, options)) != -1) {
699 case 'P': *pP = atof(optarg); break;
700 case 't': *pDt_max = atof(optarg); break;
701 case 'v': *pV = atof(optarg); break;
702 case 'T': env->T = atof(optarg); break;
703 case 'C': env->T = atof(optarg)+273.15; break;
705 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
706 assert(num_strings == 1 || num_strings == 2);
707 if (num_strings == 1)
710 ftension_group = atoi(strings[1]);
711 ftension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
712 ftension_params = tension_models[ftension_model].params;
713 free_parsed_list(num_strings, strings);
716 ftension_params = optarg;
719 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
720 assert(num_strings == 1 || num_strings == 2);
721 if (num_strings == 1)
724 utension_group = atoi(strings[1]);
725 utension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
726 utension_params = tension_models[utension_model].params;
727 free_parsed_list(num_strings, strings);
730 utension_params = optarg;
733 k_model = index_k_model(n_k_models, k_models, optarg);
736 k_models[k_model].params = optarg;
741 for (i=0; i<n; i++) {
742 push(pDomain_list, generate_domain(DS_FOLDED,
743 tension_models+ftension_model,
746 tension_models+utension_model,
749 k_models+k_model), 1);
756 for (i=0; i<n; i++) {
757 push(pDomain_list, generate_domain(DS_UNFOLDED,
758 tension_models+ftension_model,
761 tension_models+utension_model,
764 k_models+k_model), 1);
768 flags = FLAG_OUTPUT_FULL_CURVE;
771 fprintf(stderr, "unrecognized option '%c'\n", optopt);
772 /* fall through to default case */
774 help(prog_name, *pP, *pDt_max, *pV, env, n_tension_models, tension_models, ftension_model, utension_model, n_k_models, k_models, k_model);
778 /* check the arguments */
779 assert(*pDomain_list != NULL); /* no domains! */
780 assert(*pP > 0.0); assert(*pP < 1.0);
781 assert(*pDt_max > 0.0);
783 assert(env->T > 0.0);
788 <<index tension model>>=
789 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
792 for (i=0; i<n_models; i++)
793 if (STRMATCH(models[i].name, name))
796 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
804 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
807 for (i=0; i<n_models; i++)
808 if (STRMATCH(models[i].name, name))
811 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
818 <<initialize model definition>>=
819 /* requires int num_param_args and char **param_args in the current scope
821 * INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
822 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
824 #define INIT_MODEL(role, model, param_string, param_pointer) \
826 parse_list_string(param_string, SEP, '{', '}', \
827 &num_param_args, ¶m_args); \
828 if (num_param_args != model->num_params) { \
830 "%s model %s expected %d params,\n", \
831 role, model->name, model->num_params); \
833 "not the %d params in '%s'\n", \
834 num_param_args, param_string); \
835 assert(num_param_args == model->num_params); \
837 if (model->creator) \
838 param_pointer = (*model->creator)(param_args); \
840 param_pointer = NULL; \
841 free_parsed_list(num_param_args, param_args); \
846 void *generate_domain(enum domain_state_t state,
847 tension_model_getopt_t *folded_model,
849 char *folded_param_string,
850 tension_model_getopt_t *unfolded_model,
852 char *unfolded_param_string,
853 k_model_getopt_t *k_model)
855 void *folded_params, *unfolded_params, *k_params;
856 int num_param_args; /* for INIT_MODEL() */
857 char **param_args; /* for INIT_MODEL() */
860 fprintf(stderr, "generating %s ", state==DS_FOLDED ? "folded" : "unfolded");
861 fprintf(stderr, "domain (k: \"%s\", \"%s\", f: \"%s\" (%p) group %d \"%s\", u: \"%s\" (%p) group %d \"%s\")\n",
862 k_model->name, k_model->params,
863 folded_model->name, folded_model->handler, folded_group, folded_param_string,
864 unfolded_model->name, unfolded_model->handler, unfolded_group, unfolded_param_string);
866 INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
867 INIT_MODEL("unfolded", unfolded_model, unfolded_param_string, unfolded_params);
868 INIT_MODEL("k", k_model, k_model->params, k_params);
869 return create_domain(state,
870 k_model->k, k_params, k_model->destructor,
871 folded_model->handler, folded_group, folded_params, folded_model->destructor,
872 unfolded_model->handler, unfolded_group, unfolded_params, unfolded_model->destructor
879 \addcontentsline{toc}{section}{Appendicies}
881 \section{sawsim.c details}
885 The general layout of our simulation code is:
895 We include [[math.h]], so don't forget to link to the libm with `-lm'.
897 #include <assert.h> /* assert() */
898 #include <stdlib.h> /* malloc(), free(), rand() */
899 #include <stdio.h> /* fprintf(), stdout */
900 #include <string.h> /* strlen, strtok() */
901 #include <math.h> /* exp(), M_PI, sqrt() */
902 #include <sys/types.h> /* pid_t (returned by getpid()) */
903 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
906 #include "tension_balance.h"
908 #include "tension_model.h"
914 <<version definition>>
916 <<max/min definitions>>
917 <<string comparison definition and globals>>
918 <<initialize model definition>>
919 <<get options definitions>>
920 <<domain definitions>>
929 <<create/destroy domain>>
930 <<destroy domain list>>
932 <<simulation functions>>
933 <<boilerplate functions>>
936 <<boilerplate functions>>=
938 <<get option functions>>
944 srand(getpid()*time(NULL)); /* seed rand() */
948 <<flag definitions>>=
949 /* in octal b/c of prefixed '0' */
950 #define FLAG_OUTPUT_FULL_CURVE 01
951 #define FLAG_OUTPUT_UNFOLDING_FORCES 02
955 static unsigned long int flags = 0;
958 \subsection{Utilities}
961 <<max/min definitions>>=
962 #define MAX(a,b) ((a)>(b) ? (a) : (b))
963 #define MIN(a,b) ((a)<(b) ? (a) : (b))
966 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
967 <<string comparison definition and globals>>=
968 // Check if two strings match, return 1 if they do
969 static char *temp_string_A;
970 static char *temp_string_B;
971 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
972 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
973 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
974 /* +1 to also compare the '\0' */
977 We also define a macro for our [[check]] unit testing
978 <<check relative error>>=
979 #define CHECK_ERR(max_err, expected, received) \
981 fail_unless( (received-expected)/expected < max_err, \
982 "relative error %g >= %g in %s (Expected %g, received %g)", \
983 (received-expected)/expected, max_err, #received, \
984 expected, received); \
985 fail_unless(-(received-expected)/expected < max_err, \
986 "relative error %g >= %g in %s (Expected %g, received %g)", \
987 -(received-expected)/expected, max_err, #received, \
988 expected, received); \
993 int happens(double probability)
995 assert(probability >= 0.0); assert(probability <= 1.0);
996 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*/
1001 \subsection{Adaptive timesteps}
1002 \label{app.adaptive_dt}
1004 $F(x)$ increases with $x$, possibly exploding, as in the worm-like chain model,
1005 so basing the timestep on the the unfolding probability at the current tension
1006 is dangerous, and we need to search for a $dt$ for which
1007 $P(F(x+v*dt)) < P_\text{target}$.
1008 There are two cases to consider.
1009 In the most common, no domains have unfolded since the last step,
1010 and we expect the next step to be slightly shorter than the previous one.
1011 In the less common, domains did unfold in the last step,
1012 and we expect the next step to be considerably longer than the previous one.
1014 double search_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers,
1015 list_t *domain_list,
1016 environment_t *env, double x,
1017 double target_prob, double max_dt, double v)
1018 { /* Returns the timestep dt in seconds for the current folded domain.
1019 * Takes a list of tension handlers, the list of domains,
1020 * a pointer env to the environmental data, a starting separation x in m,
1021 * a target_prob between 0 and 1,
1022 * max_dt in s, stretching velocity v in m/s.
1024 double F, k, dtCur, dtU, dtUCur, dtL, dt;
1026 /* get upper bound using the current position */
1027 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x, 0); /* BUG. repeated calculation */
1028 //printf("Start with x = %g (F = %g)\n", x, F);
1029 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1030 //printf("x %g\tF %g\tk %g\n", x, F, k);
1031 dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */
1033 //printf("overshot max_dt\n");
1036 /* set a lower bound on dt too */
1039 /* The dt determined above may produce illegitimate forces or ks.
1040 * Reduce the upper bound until we have valid ks. */
1042 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1043 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1046 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1048 //printf("Try for dt = %g (F = %g)\n", dt, F);
1049 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1050 /* returned k may be -1.0 */
1051 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1052 while (k == -1.0) { /* reduce step until we hit a valid k */
1054 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1055 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1056 //printf("Try for dt = %g (F = %g)\n", dt, F);
1057 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1058 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1060 assert(dtU > 1e-14); /* timestep to valid k too small */
1061 dtUCur = target_prob / k; /* safe timestep back from x+dtU */
1063 return dt; /* dtU is safe. */
1065 /* dtU wasn't safe, lets see what would be. */
1066 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1067 dt = (dtU + dtL) / 2.0;
1068 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1069 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1070 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1071 dtCur = target_prob / k;
1072 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1073 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1075 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1076 dtU = dt; /* ... stepping out only dtCur would be safe */
1079 break; /* dtCur = dt */
1081 return MAX(dtUCur, dtL);
1085 To determine $dt$ for an array of potentially different folded domains,
1086 we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
1089 double determine_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list,
1090 environment_t *env, double x,
1091 double target_prob, double dt_max, double v)
1092 { /* Returns the timestep dt in seconds.
1093 * Takes the list of folded domains, target_prob between 0 and 1,
1094 * F in N, and T in K. */
1095 double dt=dt_max, new_dt;
1096 assert(target_prob > 0.0); assert(target_prob < 1.0);
1097 assert(dt_max > 0.0);
1099 /* .5 nm steps = v * dt */
1101 while (domain_list != NULL) {
1102 if (D(domain_list)->state == DS_FOLDED) {
1103 new_dt = search_dt(num_tension_handlers, tension_handlers, domain_list, env, x, target_prob, dt, v);
1104 dt = MIN(dt, new_dt);
1106 domain_list = domain_list->next;
1112 \subsection{Domain data}
1114 Currently domains exist in two states, folded and unfolded, and the
1115 only allowed transitions are folded $\rightarrow$ unfolded. Of
1116 course, it wouldn't be too complicated to extent this to a multi-state
1117 system, with an array containing the domains group for each possible
1118 state, and a matrix of transition-rate-calculating functions.
1119 However, at this point such generality seems unnecessary at this
1121 <<domain definitions>>=
1122 enum domain_state_t {DS_FOLDED,
1126 typedef struct domain_struct {
1127 enum domain_state_t state;
1128 one_dim_fn_t *folded_handler;
1130 one_dim_fn_t *unfolded_handler;
1132 k_func_t *k; /* function returning unfolding rate */
1133 void *folded_params; /* pointer to folded parameters */
1134 void *unfolded_params; /* pointer to unfolded parameters */
1135 void *k_params; /* pointer to k parameters */
1136 destroy_data_func_t *destroy_folded;
1137 destroy_data_func_t *destroy_unfolded;
1138 destroy_data_func_t *destroy_k;
1141 /* get the domain data for the current list node */
1142 #define D(list) ((domain_t *)(list)->d)
1143 /* get the tension params for the current list node */
1144 #define D_TP(list) (((domain_t *)(list)->d)->state == DS_FOLDED \
1145 ? ((domain_t *)(list)->d)->folded_params \
1146 : ((domain_t *)(list)->d)->unfolded_params)
1148 [[k]] is a pointer to the function determining the unfolding rate for a given tension.
1149 [[folded_params]] is a pointer to the parameters used by the function pointed to by [[k]].
1150 [[unfolded_params]] is a pointer to the parameters used by the group-appropriate handler function when determining the tension.
1151 The [[destroy_*]] pointers point to functions for freeing the memory [[*_params]].
1152 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.
1154 [[create_]] and [[destroy_domain]] are simple wrappers around [[malloc]] and [[free]].
1155 <<create/destroy domain>>=
1156 domain_t *create_domain(enum domain_state_t state,
1159 destroy_data_func_t *destroy_k,
1160 one_dim_fn_t *folded_handler,
1162 void *folded_params,
1163 destroy_data_func_t *destroy_folded,
1164 one_dim_fn_t *unfolded_handler,
1166 void *unfolded_params,
1167 destroy_data_func_t *destroy_unfolded)
1169 domain_t *ret = (domain_t *)malloc(sizeof(domain_t));
1170 assert(ret != NULL);
1171 if (state == DS_FOLDED) {
1172 assert(k != NULL); /* the pointer points somewhere valid */
1173 assert(*k != NULL); /* and there is something useful there */
1175 assert(state == DS_UNFOLDED);
1177 ret->folded_handler = folded_handler;
1178 ret->folded_group = folded_group;
1179 ret->unfolded_handler = unfolded_handler;
1180 ret->unfolded_group = unfolded_group;
1182 ret->folded_params = folded_params;
1183 ret->unfolded_params = unfolded_params;
1184 ret->k_params = k_params;
1185 ret->destroy_folded = destroy_folded;
1186 ret->destroy_unfolded = destroy_unfolded;
1187 ret->destroy_k = destroy_k;
1189 fprintf(stderr, "generated domain %p. state %d\n", ret, ret->state);
1194 void destroy_domain(domain_t *domain)
1197 //printf("domain %p & %p\n", *domain, domain);
1198 if (domain->destroy_folded)
1199 (*domain->destroy_folded)(domain->folded_params);
1200 if (domain->destroy_unfolded)
1201 (*domain->destroy_unfolded)(domain->unfolded_params);
1202 if (domain->destroy_k)
1203 (*domain->destroy_k)(domain->k_params);
1209 <<destroy domain list>>=
1210 void destroy_domain_list(list_t *domain_list)
1212 domain_list = head(domain_list);
1213 while (domain_list != NULL)
1214 destroy_domain((domain_t *) pop(&domain_list));
1218 \subsection{Domain model and group handling}
1220 <<group functions>>=
1221 <<get tension handler>>
1223 <<int comparison function>>
1224 <<find possible groups>>
1227 <<get active groups>>
1230 <<get tension handler>>=
1231 one_dim_fn_t *get_tension_handler(domain_t *domain)
1233 if (domain->state == DS_FOLDED)
1234 return domain->folded_handler;
1236 if (domain->state != DS_UNFOLDED)
1237 fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
1238 assert(domain->state == DS_UNFOLDED);
1239 return domain->unfolded_handler;
1245 int get_group(domain_t *domain)
1247 if (domain->state == DS_FOLDED)
1248 return domain->folded_group;
1250 if (domain->state != DS_UNFOLDED)
1251 fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
1252 assert(domain->state == DS_UNFOLDED);
1253 return domain->unfolded_group;
1258 We already know all possible tension classes, since they are all
1259 hardcoded into \prog. However, there may be any number of possible
1260 groups. We define a function [[find_possible_groups]] to search for
1261 possible (not neccessarily active) groups. It can search for
1262 subgroups of a single [[handler]], or by setting [[handler]] to
1263 [[NULL]], subgroups of \emph{any} handler. It returns a list of group
1265 <<find possible groups>>=
1266 list_t *find_possible_groups(list_t *list, one_dim_fn_t *handler) {
1269 while (list != NULL) {
1270 if (handler == NULL || get_tension_handler(D(list)) == handler) {
1271 push(&ret, &D(list)->folded_group, 1);
1272 push(&ret, &D(list)->unfolded_group, 1);
1278 return ret; /* no groups with this handler, no processing needed */
1280 /* sort the ret list, and remove duplicates */
1281 sort(&ret, &int_comparison);
1282 unique(&ret, &int_comparison);
1284 fprintf(stderr, "handler %p possible groups:", handler);
1286 while (list != NULL) {
1287 fprintf(stderr, " %d", *((int *)list->d));
1290 fprintf(stderr, "\n");
1296 Search a [[list]] of domains, and determine whether a particular model
1297 class and group number combination is active.
1298 <<is group active>>=
1299 int is_group_active(list_t *list, one_dim_fn_t *handler, int group)
1302 while (list != NULL) {
1303 if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group)
1311 Search a [[list]] of domains, and return all domains matching a
1312 particular model class and group number.
1314 list_t *get_group_list(list_t *list, one_dim_fn_t *handler, int group)
1318 while (list != NULL) {
1319 if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group) {
1320 push(&ret, D_TP(list), 1); /* add a pointer to the appropriate tension parameters to our new list. */
1322 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);
1330 Because all the node data in lists returned by [[get_group_list]] is
1331 also in the main domain list, you shouldn't destroy the node data
1332 popped off when destroying the group lists. It will all get cleaned
1333 up when the main domain list is destroyed.
1335 Put all this together to scan through a list of domains and construct
1336 an array of all the active groups.
1337 <<get active groups>>=
1338 void get_active_groups(list_t *list,
1339 int num_tension_handlers, one_dim_fn_t **tension_handlers,
1340 int *pNum_active_groups, one_dim_fn_t ***pPer_group_handlers, void ***pActive_groups)
1342 void **active_groups=NULL;
1343 one_dim_fn_t *handler, **per_group_handlers=NULL;
1344 int i, num_active_groups, *group;
1345 list_t *possible_group_numbers, *active_group_list, *active_groups_list=NULL, *per_group_handlers_list=NULL;
1347 /* get all the active groups in a list */
1348 for (i=0; i<num_tension_handlers; i++) {
1349 handler = tension_handlers[i];
1350 if (handler == NULL) continue; /* tensionless handler */
1351 possible_group_numbers = head(find_possible_groups(list, handler));
1352 while (possible_group_numbers != NULL) {
1353 group = (int *)pop(&possible_group_numbers);
1354 if (is_group_active(list, handler, *group) == 1) {
1355 active_group_list = get_group_list(list, handler, *group);
1356 push(&active_groups_list, active_group_list, 1);
1357 push(&per_group_handlers_list, handler, 1);
1359 fprintf(stderr, "active group %p for %p %d headed by tension parameters %p\n", active_group_list, handler, *group, head(active_group_list)->d);
1360 list_print(stderr, active_group_list, "active group");
1366 /* allocate the array we'll be returning */
1367 num_active_groups = list_length(active_groups_list);
1368 active_groups = (void **) malloc(sizeof(void *)*num_active_groups);
1369 assert(active_groups != NULL);
1370 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1371 assert(per_group_handlers != NULL);
1373 /* populate the array we'll be returning */
1374 active_groups_list = head(active_groups_list);
1375 for (i=0; i<num_active_groups; i++) {
1376 handler = pop(&per_group_handlers_list);
1377 per_group_handlers[i] = handler;
1379 active_groups[i] = (void *) malloc(sizeof(tension_handler_data_t));
1380 assert(active_groups[i] != NULL);
1381 active_group_list = pop(&active_groups_list);
1382 ((tension_handler_data_t *)active_groups[i])->group = active_group_list;
1383 ((tension_handler_data_t *)active_groups[i])->env = NULL;
1384 ((tension_handler_data_t *)active_groups[i])->persist = NULL;
1386 assert(active_groups_list == NULL);
1387 assert(per_group_handlers_list == NULL);
1389 *pNum_active_groups = num_active_groups;
1390 *pPer_group_handlers = per_group_handlers;
1391 *pActive_groups = active_groups;
1396 \section{String parsing}
1398 For handling command line arguments, we need parse delimited strings (e.g. [["param1,param2,param3..."]]).
1399 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1400 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1401 need to take care of parsing those parameters themselves.
1402 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1406 <<parse definitions>>
1407 <<parse declarations>>
1410 <<parse module makefile lines>>=
1411 NW_SAWSIM_MODS += parse
1412 CHECK_BINS += check_parse
1416 <<parse definitions>>=
1417 #define SEP ',' /* argument separator character */
1420 <<parse declarations>>=
1421 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1422 int *num, char ***string_array);
1423 extern void free_parsed_list(int num, char **string_array);
1426 [[parse_list_string]] allocates memory, don't forget to free it
1427 afterward with [[free_parsed_list]]. It does not alter the original.
1429 The string may start off with a [[deeper]] character
1430 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1431 [[x,y]], where the pointer is one character in on the copied string.
1432 However, when we go to free the memory, we need a pointer to the
1433 beginning of the string. In order to accommodate this for a string
1434 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1435 the first $N$ elements point to the separated arguments, and let the
1436 last element point to the start of the copied string regardless of
1438 <<parse delimited list string functions>>=
1439 /* TODO, split out into parse.hc */
1440 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1443 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1444 if (string[i] == deeper) {depth++;}
1445 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1451 void parse_list_string(char *string, char sep, char deeper, char shallower,
1452 int *num, char ***string_array)
1454 char *str=NULL, **ret=NULL;
1456 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
1458 *string_array = NULL;
1461 /* make a copy of the string, so we don't change the original */
1462 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1463 assert(str != NULL);
1464 strcpy(str, string); /* we know str is long enough */
1465 /* count the number of regions, so we can allocate pointers to them */
1468 n++; i++; /* move on to next argument */
1469 i += next_delim_index(str+i, sep, deeper, shallower);
1470 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1471 } while (str[i] != '\0');
1472 ret = (char **)malloc(sizeof(char *)*(n+1));
1473 assert(ret != NULL);
1474 /* replace the separators with '\0' & assign pointers */
1475 ret[n] = str; /* point to the front of the copied string */
1478 for(i=1; i<n; i++) {
1479 j += next_delim_index(str+j, sep, deeper, shallower);
1481 ret[i] = str+j; /* point to the separated arguments */
1483 /* strip off bounding braces, if any */
1484 for(i=0; i<n; i++) {
1485 if (ret[i][0] == deeper) {
1489 if (ret[i][strlen(ret[i])-1] == shallower)
1490 ret[i][strlen(ret[i])-1] = '\0';
1493 *string_array = ret;
1496 void free_parsed_list(int num, char **string_array)
1499 /* frees all items, since they were allocated as a single string*/
1500 free(string_array[num]);
1501 /* and free the array of pointers */
1507 \subsection{Parsing implementation}
1513 <<parse delimited list string functions>>
1517 #include <assert.h> /* assert() */
1518 #include <stdlib.h> /* NULL */
1519 #include <stdio.h> /* fprintf(), stdout *//*!!*/
1520 #include <string.h> /* strlen() */
1524 \subsection{Parsing unit tests}
1526 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1528 <<parse check includes>>
1529 <<string comparison definition and globals>>
1530 <<check relative error>>
1531 <<parse test suite>>
1532 <<main check program>>
1535 <<parse check includes>>=
1536 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
1537 #include <stdio.h> /* printf() */
1538 #include <assert.h> /* assert() */
1539 #include <string.h> /* strlen() */
1544 <<parse test suite>>=
1545 <<parse list string tests>>
1546 <<parse suite definition>>
1549 <<parse suite definition>>=
1550 Suite *test_suite (void)
1552 Suite *s = suite_create ("k model");
1553 <<parse list string test case defs>>
1555 <<parse list string test case add>>
1560 <<parse list string tests>>=
1563 START_TEST(test_next_delim_index)
1565 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
1566 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
1567 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
1568 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
1569 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
1573 START_TEST(test_parse_list_null)
1577 parse_list_string(NULL, SEP, '{', '}',
1578 &num_param_args, ¶m_args);
1579 fail_unless(num_param_args == 0, NULL);
1580 fail_unless(param_args == NULL, NULL);
1583 START_TEST(test_parse_list_single_simple)
1587 parse_list_string("arg", SEP, '{', '}',
1588 &num_param_args, ¶m_args);
1589 fail_unless(num_param_args == 1, NULL);
1590 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
1593 START_TEST(test_parse_list_single_compound)
1597 parse_list_string("{x,y,z}", SEP, '{', '}',
1598 &num_param_args, ¶m_args);
1599 fail_unless(num_param_args == 1, NULL);
1600 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
1603 START_TEST(test_parse_list_double_simple)
1607 parse_list_string("abc,def", SEP, '{', '}',
1608 &num_param_args, ¶m_args);
1609 fail_unless(num_param_args == 2, NULL);
1610 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
1611 fail_unless(STRMATCH(param_args[1],"def"), NULL);
1615 <<parse list string test case defs>>=
1616 TCase *tc_parse_list_string = tcase_create("parse list string");
1618 <<parse list string test case add>>=
1619 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
1620 tcase_add_test(tc_parse_list_string, test_parse_list_null);
1621 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
1622 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
1623 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
1624 suite_add_tcase(s, tc_parse_list_string);
1628 \section{Unit tests}
1630 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1639 <<main check program>>
1651 <<determine dt tests>>
1653 <<does domain unfold tests>>
1654 <<randomly unfold domains tests>>
1655 <<suite definition>>
1658 <<suite definition>>=
1659 Suite *test_suite (void)
1661 Suite *s = suite_create ("sawsim");
1662 <<F test case defs>>
1663 <<determine dt test case defs>>
1664 <<happens test case defs>>
1665 <<does domain unfold test case defs>>
1666 <<randomly unfold domains test case defs>>
1669 <<determine dt test case add>>
1670 <<happens test case add>>
1671 <<does domain unfold test case add>>
1672 <<randomly unfold domains test case add>>
1675 tcase_add_checked_fixture(tc_strip_address,
1676 setup_strip_address,
1677 teardown_strip_address);
1683 <<main check program>>=
1687 Suite *s = test_suite();
1688 SRunner *sr = srunner_create(s);
1689 srunner_run_all(sr, CK_ENV);
1690 number_failed = srunner_ntests_failed(sr);
1692 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
1696 \subsection{F tests}
1698 <<worm-like chain tests>>
1700 <<F test case defs>>=
1701 <<worm-like chain test case def>>
1703 <<F test case add>>=
1704 <<worm-like chain test case add>>
1707 <<worm-like chain tests>>=
1708 START_TEST(test_wlc_at_zero)
1710 double T=1.0, L=1.0, p=0.1, x=0.0;
1711 fail_unless(wlc(x, T, p, L)==0, NULL);
1715 START_TEST(test_wlc_at_half)
1717 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
1718 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
1719 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
1721 * linear term = x/L = 0.5
1722 * nonlinear + linear = 0.75 + 0.5 = 1.25
1723 * wlc = 10e21*1.25 = 12.5e21
1725 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
1726 "wlc(%g, %g, %g, %g) = %g != %g",
1727 x, T, p, L, wlc(x, T, p, L), 12.5e21);
1731 <<worm-like chain test case def>>=
1732 TCase *tc_wlc = tcase_create("WLC");
1735 <<worm-like chain test case add>>=
1736 tcase_add_test(tc_wlc, test_wlc_at_zero);
1737 tcase_add_test(tc_wlc, test_wlc_at_half);
1738 suite_add_tcase(s, tc_wlc);
1741 \subsection{Model tests}
1743 Check the searching with [[linear_k]].
1744 Check overwhelming force treatment with the heavyside-esque step [[?]].
1745 <<determine dt tests>>=
1746 double linear_k(double F, environment_t *env, void *params)
1748 double Fnot = *(double *)params;
1752 START_TEST(test_determine_dt_linear_k)
1755 double dt_max=3.0, Fnot=3.0;
1756 double F[]={0,1,2,3,4,5,6};
1757 domain_t dom; /* use both parts at once for folded/unfolded */
1761 dom->next = dom->prev = NULL;
1762 dom->k_func_t = linear_k;
1763 dom->folded_params = &Fnot;
1764 dom->unfolded_params = !!!!!!!!!
1765 dom->destroy_folded = dom->destroy_unfolded = NULL;
1766 for( i=0; i < sizeof(F)/sizeof(double); i++) {
1767 fail_unless(determine_dt(folded, unfolded, targ_p, dt_max, v, x, T)==1, NULL);
1773 <<determine dt test case defs>>=
1774 TCase *tc_determine_dt = tcase_create("Determine dt");
1776 <<determine dt test case add>>=
1777 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
1778 suite_add_tcase(s, tc_determine_dt);
1783 <<happens test case defs>>=
1785 <<happens test case add>>=
1788 <<does domain unfold tests>>=
1790 <<does domain unfold test case defs>>=
1792 <<does domain unfold test case add>>=
1795 <<randomly unfold domains tests>>=
1797 <<randomly unfold domains test case defs>>=
1799 <<randomly unfold domains test case add>>=
1803 \section{Balancing group extensions}
1804 \label{app.tension_balance}
1806 For a given total extension $x$ (of the piezo), the various domain
1807 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
1808 amounts, and we need to tweak the portion that each extends to balance
1809 the tension amongst the active groups. Since the tension balancing is
1810 separable from the bulk of the simulation, we break it out into a
1811 separate module. The interface is defined in [[tension_balance.h]],
1812 the implementation in [[tension_balance.c]], and the unit testing in
1813 [[check_tension_balance.c]]
1815 <<tension-balance.h>>=
1817 <<tension balance function declaration>>
1820 <<tension balance functions>>=
1821 <<one dimensional bracket>>
1822 <<one dimensional solve>>
1824 <<tension balance function>>
1827 <<tension balance module makefile lines>>=
1828 NW_SAWSIM_MODS += tension_balance
1829 CHECK_BINS += check_tension_balance
1830 check_tension_balance_MODS =
1833 The entire force balancing problem reduces to a solving two nested
1834 one-dimensional root-finding problems. First we define the one
1835 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
1836 non-empty group). $x(x_0)$ is also strictly monotonically increasing,
1837 so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
1838 stores the last successful value of $x$, since we expect to be taking
1839 small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
1840 indicates that the groups have changed.
1841 <<tension balance function declaration>>=
1842 double tension_balance(int num_tension_groups,
1843 one_dim_fn_t **tension_handlers,
1848 <<tension balance function>>=
1849 double tension_balance(int num_tension_groups,
1850 one_dim_fn_t **tension_handlers,
1855 static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL};
1856 double F, xo_guess, xo, lb, ub;
1857 double min_relx=1e-6, min_rely=1e-6;
1858 int max_steps=200, i;
1860 assert(num_tension_groups > 0);
1861 assert(tension_handlers != NULL);
1862 assert(params != NULL);
1863 assert(num_tension_groups > 0);
1865 if (num_tension_groups == 1) { /* only one group, no balancing required */
1868 //fprintf(stderr, "balancing for x = %g with ", x);
1869 //for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
1870 //fprintf(stderr, "\n");
1871 if (last_x == -1) { /* new group setup, reset x_xo_data */
1872 /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
1873 if (x_xo_data.xi != NULL)
1875 /* construct new x_xo_data */
1876 x_xo_data.n_groups = num_tension_groups;
1877 x_xo_data.tension_handler = tension_handlers;
1878 x_xo_data.handler_data = params;
1879 x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
1880 for (i=0; i<num_tension_groups; i++)
1881 x_xo_data.xi[i] = -1.0;
1883 if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
1884 if (x == 0) xo_guess = length_scale;
1885 else xo_guess = x/num_tension_groups;
1887 fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
1889 oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
1890 } else { /* work off the last balanced point */
1892 fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
1895 lb = x_xo_data.xi[0];
1896 ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */
1897 } else if (x < last_x) {
1898 lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
1899 ub = x_xo_data.xi[0];
1900 } else { /* x == last_x */
1901 fprintf(stderr, "not moving\n");
1902 lb= x_xo_data.xi[0];
1903 ub= x_xo_data.xi[0];
1907 fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
1908 __FUNCTION__, x, lb, ub);
1910 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
1911 F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
1913 if (num_tension_groups > 1) {
1914 fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
1915 for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
1916 fprintf(stderr, "\n");
1920 F = (*tension_handlers[0])(xo, params[0]);
1925 <<tension balance internal definitions>>=
1926 <<x of x0 definitions>>
1929 <<x of x0 definitions>>=
1930 typedef struct x_of_xo_data_struct {
1932 one_dim_fn_t **tension_handler; /* array of fn pointers */
1933 void **handler_data; /* array of void* pointers */
1939 double x_of_xo(double xo, void *pdata)
1941 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
1942 double F, x=0, xi, xi_guess, lb, ub;
1944 double min_relx=1e-6, min_rely=1e-6;
1946 assert(data->n_groups > 0);
1948 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
1950 fprintf(stderr, "%s: finding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
1953 if (data->xi) data->xi[0] = xo;
1954 for (i=1; i < data->n_groups; i++) {
1955 if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
1956 else xi_guess = data->xi[i];
1958 fprintf(stderr, "%s: searching for proper x[%d] for F[%d](x[%d]) = %g N\n", __FUNCTION__, i, i, i, F);
1960 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
1961 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
1963 fprintf(stderr, "%s: found F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
1967 if (data->xi) data->xi[i] = xi;
1970 fprintf(stderr, "%s: found x(x0=%g) = %g\n", __FUNCTION__, xo, x);
1976 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
1977 number of steps for monotonically increasing $f(x)$. Simple
1978 bisection, so it's robust and converges fairly quickly. You can set
1979 the maximum number of steps to take, but if you have enough processor
1980 power, it's better to limit your precision with the relative limits
1981 [[min_relx]] and [[y]]. These ensure that the width of the bracket is
1982 small on the length scale of it's larger side. Note that \emph{both}
1983 relative limits must be satisfied for the search to stop.
1984 <<one dimensional function definition>>=
1985 /* equivalent to gsl_func_t */
1986 typedef double one_dim_fn_t(double x, void *params);
1988 <<one dimensional solve declaration>>=
1989 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
1990 double min_relx, double min_rely, int max_steps,
1994 <<one dimensional solve>>=
1995 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
1996 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
1997 double min_relx, double min_rely, int max_steps,
2000 double yx, ylb, yub, x;
2003 ylb = (*f)(lb, params);
2004 yub = (*f)(ub, params);
2006 /* check some simple cases */
2008 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bounded */
2009 else return lb; /* any x on the range [lb, ub] would work */
2011 if (ylb == y) { x=lb; yx=ylb; goto end; }
2012 if (yub == y) { x=ub; yx=yub; goto end; }
2015 fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2017 assert(ylb < y); assert(yub > y);
2018 for (i=0; i<max_steps; i++) {
2020 yx = (*f)(x, params);
2022 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);
2024 if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2025 else if (yx > y) { ub=x; yub=yx; }
2026 else /* yx < y */ { lb=x; ylb=yx; }
2031 fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2037 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2038 Generate bracketing $x$ values through bisection or exponential growth.
2039 <<one dimensional bracket declaration>>=
2040 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2043 <<one dimensional bracket>>=
2044 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2046 double yx, step, x=xguess;
2048 yx = (*f)(x, params);
2050 fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2057 fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2061 if (x == 0) x = length_scale/2.0; /* HACK */
2064 yx = (*f)(x, params);
2066 fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2071 fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2074 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2078 \subsection{Balancing implementation}
2080 <<tension-balance.c>>=
2083 <<tension balance includes>>
2084 #include "tension_balance.h"
2086 double length_scale = 1e-10; /* HACK */
2088 <<tension balance internal definitions>>
2089 <<tension balance functions>>
2092 <<tension balance includes>>=
2093 #include <assert.h> /* assert() */
2094 #include <stdlib.h> /* NULL */
2095 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2096 #include <stdio.h> /* fprintf(), stdout */
2100 \subsection{Balancing unit tests}
2102 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2103 <<check-tension-balance.c>>=
2104 <<tension balance check includes>>
2105 <<tension balance test suite>>
2106 <<main check program>>
2109 <<tension balance check includes>>=
2110 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2111 #include <assert.h> /* assert() */
2114 #include "tension_balance.h"
2117 <<tension balance test suite>>=
2118 <<tension balance function tests>>
2119 <<tension balance suite definition>>
2122 <<tension balance suite definition>>=
2123 Suite *test_suite (void)
2125 Suite *s = suite_create ("tension balance");
2126 <<tension balance function test case defs>>
2128 <<tension balance function test case adds>>
2133 <<tension balance function tests>>=
2134 <<check relative error>>
2136 double hooke(double x, void *pK)
2139 return *((double*)pK) * x;
2142 START_TEST(test_single_function)
2144 double k=5, x=3, last_x=2.0, F;
2145 one_dim_fn_t *handlers[] = {&hooke};
2146 void *data[] = {&k};
2147 F = tension_balance(1, handlers, data, last_x, x);
2148 fail_unless(F = k*x, NULL);
2153 We can also test balancing two springs with different spring constants.
2154 The analytic solution for a general number of springs is given in Appendix \ref{app.math_hooke}.
2155 <<tension balance function tests>>=
2156 START_TEST(test_double_hooke)
2158 double k1=5, k2=4, x=3, last_x=2.0, F, Fe, x1e, x2e;
2159 one_dim_fn_t *handlers[] = {&hooke, &hooke};
2160 void *data[] = {&k1, &k2};
2161 F = tension_balance(2, handlers, data, last_x, x);
2165 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2166 //CHECK_ERR(1e-6, x1e, xi[0]);
2167 //CHECK_ERR(1e-6, x2e, xi[1]);
2168 CHECK_ERR(1e-6, Fe, F);
2173 <<tension balance function test case defs>>=
2174 TCase *tc_tbfunc = tcase_create("tension balance function");
2177 <<tension balance function test case adds>>=
2178 tcase_add_test(tc_tbfunc, test_single_function);
2179 tcase_add_test(tc_tbfunc, test_double_hooke);
2180 suite_add_tcase(s, tc_tbfunc);
2185 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2186 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2187 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2191 <<list definitions>>
2192 <<list declarations>>
2195 <<list declarations>>=
2196 <<head and tail declarations>>
2197 <<list length declaration>>
2198 <<push declaration>>
2200 <<list destroy declaration>>
2201 <<list to array declaration>>
2202 <<move declaration>>
2203 <<sort declaration>>
2204 <<unique declaration>>
2205 <<list print declaration>>
2209 <<create and destroy node>>
2222 <<list module makefile lines>>=
2223 NW_SAWSIM_MODS += list
2224 CHECK_BINS += check_list
2228 <<list definitions>>=
2229 typedef struct list_struct {
2230 struct list_struct *next;
2231 struct list_struct *prev;
2235 <<comparison function typedef>>
2238 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2239 <<head and tail declarations>>=
2240 list_t *head(list_t *list);
2241 list_t *tail(list_t *list);
2244 list_t *head(list_t *list)
2248 while (list->prev != NULL) {
2254 list_t *tail(list_t *list)
2258 while (list->next != NULL) {
2265 <<list length declaration>>=
2266 int list_length(list_t *list);
2269 int list_length(list_t *list)
2276 while (list->next != NULL) {
2284 [[push]] inserts a node relative to the current position in the list
2285 without changing the current position.
2286 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2287 so we set it to that of the pushed domain.
2288 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2289 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2290 <<push declaration>>=
2291 void push(list_t **pList, void *data, int below);
2294 void push(list_t **pList, void *data, int below)
2296 list_t *list, *new_node;
2297 assert(pList != NULL);
2299 new_node = create_node(data);
2302 else if (below==0) { /* insert above */
2303 if (list->prev != NULL)
2304 list->prev->next = new_node;
2305 new_node->prev = list->prev;
2306 list->prev = new_node;
2307 new_node->next = list;
2308 } else { /* insert below */
2309 if (list->next != NULL)
2310 list->next->prev = new_node;
2311 new_node->next = list->next;
2312 list->next = new_node;
2313 new_node->prev = list;
2318 [[pop]] removes the current domain node, moving the current position
2319 to the node after it, unless that node is [[NULL]], in which case move
2320 the current position to the node before it.
2321 <<pop declaration>>=
2322 void *pop(list_t **pList);
2325 void *pop(list_t **pList)
2327 list_t *list, *popped;
2329 assert(pList != NULL);
2331 assert(list != NULL); /* not an empty list */
2333 /* bypass the popped node */
2334 if (list->prev != NULL)
2335 list->prev->next = list->next;
2336 if (list->next != NULL) {
2337 list->next->prev = list->prev;
2338 *pList = list->next;
2340 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
2342 destroy_node(popped);
2347 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
2348 If the cleanup function is [[NULL]], no cleanup function is called.
2349 The nodes are not popped in any particular order.
2350 <<list destroy declaration>>=
2351 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
2354 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
2358 assert(pList != NULL);
2361 assert(list != NULL); /* not an empty list */
2362 while (list != NULL) {
2364 if (destroy != NULL)
2370 [[list_to_array]] converts a list to an array of pointers, preserving order.
2371 <<list to array declaration>>=
2372 void list_to_array(list_t *list, int *length, void ***array);
2375 void list_to_array(list_t *list, int *pLength, void ***pArray)
2379 assert(list != NULL);
2380 assert(pLength != NULL);
2381 assert(pArray != NULL);
2383 length = list_length(list);
2384 array = (void **)malloc(sizeof(void *)*length);
2385 assert(array != NULL);
2388 while (list != NULL) {
2389 array[i++] = list->d;
2397 [[move]] moves the current element along the list in either direction.
2398 <<move declaration>>=
2399 void move(list_t **pList, int downstream);
2402 void move(list_t **pList, int downstream)
2404 list_t *list, *flipper;
2406 assert(pList != NULL);
2407 list = *pList; /* a>B>c>d */
2408 assert(list != NULL); /* not an empty list */
2409 if (downstream == 0)
2410 flipper = list->prev; /* flipper is A */
2412 flipper = list->next; /* flipper is C */
2413 /* ensure there is somewhere to go in the direction we're heading */
2414 assert(flipper != NULL);
2415 /* remove the current element */
2416 data = pop(&list); /* data=B, a>C>d */
2417 /* and put it back in in the correct spot */
2419 if (downstream == 0) {
2420 push(&list, data, 0); /* b>A>c>d */
2421 list = list->prev; /* B>a>c>d */
2423 push(&list, data, 1); /* a>C>b>d */
2424 list = list->next; /* a>c>B>d */
2430 [[sort]] sorts a list from smallest to largest according to a user
2431 specified comparison.
2432 <<comparison function typedef>>=
2433 typedef int (comparison_fn_t)(void *A, void *B);
2436 <<int comparison function>>=
2437 int int_comparison(void *A, void *B)
2442 if (a > b) return 1;
2443 else if (a == b) return 0;
2447 <<sort declaration>>=
2448 void sort(list_t **pList, comparison_fn_t *comp);
2450 Here we use the bubble sort algorithm.
2452 void sort(list_t **pList, comparison_fn_t *comp)
2456 assert(pList != NULL);
2458 assert(list != NULL); /* not an empty list */
2459 while (stable == 0) {
2462 while (list->next != NULL) {
2463 if ((*comp)(list->d, list->next->d) > 0) {
2465 move(&list, 1 /* downstream */);
2474 [[unique]] removes duplicates from a sorted list according to a user
2475 specified comparison. Don't do this unless you have other pointers
2476 any dynamically allocated data in your list, because [[unique]] just
2477 drops any non-unique data without freeing it.
2478 <<unique declaration>>=
2479 void unique(list_t **pList, comparison_fn_t *comp);
2482 void unique(list_t **pList, comparison_fn_t *comp)
2485 assert(pList != NULL);
2487 assert(list != NULL); /* not an empty list */
2489 while (list->next != NULL) {
2490 if ((*comp)(list->d, list->next->d) == 0) {
2499 [[list_print]] prints a list to a [[FILE *]] stream.
2500 <<list print declaration>>=
2501 void list_print(FILE *file, list_t *list, char *list_name);
2504 void list_print(FILE *file, list_t *list, char *list_name)
2506 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
2508 while (list != NULL) {
2509 fprintf(file, " %p", list->d);
2512 fprintf(file, "\n");
2516 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
2517 <<create and destroy node>>=
2518 list_t *create_node(void *data)
2520 list_t *ret = (list_t *)malloc(sizeof(list_t));
2521 assert(ret != NULL);
2528 void destroy_node(list_t *node)
2534 The user must free the data pointed to by the node on their own.
2536 \subsection{List implementation}
2546 #include <assert.h> /* assert() */
2547 #include <stdlib.h> /* malloc(), free(), rand() */
2548 #include <stdio.h> /* fprintf(), stdout, FILE */
2549 #include "global.h" /* destroy_data_func_t */
2552 \subsection{List unit tests}
2554 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2556 <<list check includes>>
2558 <<main check program>>
2561 <<list check includes>>=
2562 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2563 #include <stdio.h> /* FILE */
2569 <<list test suite>>=
2572 <<list suite definition>>
2575 <<list suite definition>>=
2576 Suite *test_suite (void)
2578 Suite *s = suite_create ("list");
2579 <<push test case defs>>
2580 <<pop test case defs>>
2582 <<push test case adds>>
2583 <<pop test case adds>>
2589 START_TEST(test_push)
2594 push(&list, (void *)i, 1);
2595 fail_unless(list == head(list), NULL);
2596 fail_unless( (int)list->d == 0 );
2597 for (i=0; i<n; i++) {
2598 /* we expect 0, n-1, n-2, ..., 1 */
2601 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
2607 <<push test case defs>>=
2608 TCase *tc_push = tcase_create("push");
2611 <<push test case adds>>=
2612 tcase_add_test(tc_push, test_push);
2613 suite_add_tcase(s, tc_push);
2618 <<pop test case defs>>=
2620 <<pop test case adds>>=
2623 \section{Function string evaluation}
2625 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).
2626 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
2627 The solution is to run a scripting language as a subprocess accessed via pipes.
2628 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
2630 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
2631 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
2632 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.
2633 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
2635 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.
2636 Then you can either statically or dynamically link to those hardcoded functions.
2637 While much less flexible, this approach would be a fairly simple-to-implement fallback.
2639 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
2640 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
2644 <<string eval setup declaration>>
2645 <<string eval function declaration>>
2646 <<string eval teardown declaration>>
2649 <<string eval module makefile lines>>=
2650 NW_SAWSIM_MODS += string_eval
2651 CHECK_BINS += check_string_eval
2652 check_string_eval_MODS =
2655 For an introduction to POSIX process control, see\\
2656 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
2657 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
2658 and of course, the relavent [[man]] pages.
2660 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
2661 [[execvp]] replaces the calling process' program with a new program.
2662 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
2663 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
2664 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
2665 The new program needs command line arguments, just like it would if you were running it from a shell.
2666 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
2667 with the final array entry being a [[NULL]] pointer.
2669 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
2670 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
2671 <<string eval subprocess definitions>>=
2672 #define SUBPROCESS "python"
2673 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
2674 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
2675 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
2677 The [[i]] option lets Python know that it should run in interactive mode.
2678 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
2679 In interactive mode, python acts on every instruction as soon as it is recieved.
2680 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.
2681 %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.
2683 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
2684 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
2685 [[fork]] returns two copies of the same program, executing the original code.
2686 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
2687 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
2689 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.
2690 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
2691 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
2692 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
2693 <<string eval pipe definitions>>=
2694 #define PIPE_READ 0 /* the end you read from */
2695 #define PIPE_WRITE 1 /* the end you write to */
2697 #define STDIN 0 /* initial index of stdin pair */
2698 #define STDOUT 2 /* initial index of stdout pair */
2701 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
2703 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.
2705 <<string eval setup declaration>>=
2706 extern void string_eval_setup(FILE **pIn, FILE **pOut);
2708 <<string eval setup definition>>=
2709 void string_eval_setup(FILE **pIn, FILE **pOut)
2712 int pfd[4]; /* file descriptors for stdin and stdout */
2714 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
2715 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
2717 pid = fork(); /* split process into two copies */
2719 if (pid == -1) { /* fork error */
2720 perror("fork error");
2722 } else if (pid == 0) { /* child */
2723 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
2724 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
2725 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
2726 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
2727 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
2728 perror("exec error"); /* exec shouldn't return */
2730 } else { /* parent */
2731 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
2732 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
2733 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
2734 if ( *pIn == NULL ) {
2735 perror("fdopen (in)");
2738 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
2739 if ( *pOut == NULL ) {
2740 perror("fdopen (out)");
2747 To use the evaluating subprocess, we just pipe in our command, and read out the results.
2748 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
2749 <<string eval function declaration>>=
2750 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
2752 <<string eval function definition>>=
2753 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
2756 rc = fprintf(in, "%s", input);
2757 assert(rc == strlen(input));
2760 alarm(1); /* set a one second timeout on the read */
2761 assert( fgets(output, buflen, out) != NULL );
2762 alarm(0); /* cancel the timeout */
2763 //fprintf(stderr, "eval: %s ----> %s", input, output);
2766 The [[alarm]] calls set and clear a timeout on the returned output.
2767 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.
2768 This protects against invalid input for which a line of output is not printed to [[stdout]].
2769 Other invalid inputs (e.g. those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
2770 If you are getting strange results, check your python code seperately. TODO, better error handling.
2772 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
2773 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
2774 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
2775 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
2776 <<string eval teardown declaration>>=
2777 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
2780 <<string eval teardown definition>>=
2781 void string_eval_teardown(FILE **pIn, FILE **pOut)
2786 /* redirect Python's stderr */
2787 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
2791 assert( fclose(*pIn) == 0 );
2793 assert( fclose(*pOut) == 0 );
2796 /* wait for python to exit */
2798 pid = wait(&stat_loc);
2805 if (WIFEXITED(stat_loc)) {
2806 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
2807 } else if (WIFSIGNALED(stat_loc)) {
2808 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
2813 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
2815 \subsection{String evaluation implementation}
2819 <<string eval includes>>
2820 #include "string_eval.h"
2821 <<string eval internal definitions>>
2822 <<string eval functions>>
2825 <<string eval includes>>=
2826 #include <assert.h> /* assert() */
2827 #include <stdlib.h> /* NULL */
2828 #include <stdio.h> /* fprintf(), stdout, fdopen() */
2829 #include <string.h> /* strlen() */
2830 #include <sys/types.h> /* pid_t */
2831 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
2832 #include <wait.h> /* wait() */
2835 <<string eval internal definitions>>=
2836 <<string eval subprocess definitions>>
2837 <<string eval pipe definitions>>
2840 <<string eval functions>>=
2841 <<string eval setup definition>>
2842 <<string eval function definition>>
2843 <<string eval teardown definition>>
2846 \subsection{String evaluation unit tests}
2848 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2849 <<check-string-eval.c>>=
2850 <<string eval check includes>>
2851 <<string comparison definition and globals>>
2852 <<string eval test suite>>
2853 <<main check program>>
2856 <<string eval check includes>>=
2857 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
2858 #include <stdio.h> /* printf() */
2859 #include <assert.h> /* assert() */
2860 #include <string.h> /* strlen() */
2861 #include <signal.h> /* SIGKILL */
2863 #include "string_eval.h"
2866 <<string eval test suite>>=
2867 <<string eval tests>>
2868 <<string eval suite definition>>
2871 <<string eval suite definition>>=
2872 Suite *test_suite (void)
2874 Suite *s = suite_create ("string eval");
2875 <<string eval test case defs>>
2877 <<string eval test case add>>
2882 <<string eval tests>>=
2883 START_TEST(test_setup_teardown)
2886 string_eval_setup(&in, &out);
2887 string_eval_teardown(&in, &out);
2890 START_TEST(test_invalid_command)
2893 char input[80], output[80]={};
2894 string_eval_setup(&in, &out);
2895 sprintf(input, "print ABCDefg\n");
2896 string_eval(in, out, input, 80, output);
2897 string_eval_teardown(&in, &out);
2900 START_TEST(test_simple_eval)
2903 char input[80], output[80]={};
2904 string_eval_setup(&in, &out);
2905 sprintf(input, "print 3+4*5\n");
2906 string_eval(in, out, input, 80, output);
2907 fail_unless(STRMATCH(output,"23\n"), NULL);
2908 string_eval_teardown(&in, &out);
2911 START_TEST(test_multiple_evals)
2914 char input[80], output[80]={};
2915 string_eval_setup(&in, &out);
2916 sprintf(input, "print 3+4*5\n");
2917 string_eval(in, out, input, 80, output);
2918 fail_unless(STRMATCH(output,"23\n"), NULL);
2919 sprintf(input, "print (3**2 + 4**2)**0.5\n");
2920 string_eval(in, out, input, 80, output);
2921 fail_unless(STRMATCH(output,"5.0\n"), NULL);
2922 string_eval_teardown(&in, &out);
2925 START_TEST(test_eval_with_state)
2928 char input[80], output[80]={};
2929 string_eval_setup(&in, &out);
2930 sprintf(input, "print 3+4*5\n");
2931 fprintf(in, "A = 3\n");
2932 sprintf(input, "print A*3\n");
2933 string_eval(in, out, input, 80, output);
2934 fail_unless(STRMATCH(output,"9\n"), NULL);
2935 string_eval_teardown(&in, &out);
2939 <<string eval test case defs>>=
2940 TCase *tc_string_eval = tcase_create("string_eval");
2942 <<string eval test case add>>=
2943 tcase_add_test(tc_string_eval, test_setup_teardown);
2944 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
2945 tcase_add_test(tc_string_eval, test_simple_eval);
2946 tcase_add_test(tc_string_eval, test_multiple_evals);
2947 tcase_add_test(tc_string_eval, test_eval_with_state);
2948 suite_add_tcase(s, tc_string_eval);
2952 \section{Accelerating function evaluation}
2954 My first version-0.3 code was running very slowly.
2955 With the options suggested in the help
2956 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
2957 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
2958 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
2959 That's only 15 calls per solution, so the search algorithm seems reasonable.
2960 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
2963 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
2967 <<accel k module makefile lines>>=
2968 NW_SAWSIM_MODS += accel_k
2969 #CHECK_BINS += check_accel_k
2970 check_accel_k_MODS =
2974 #include <assert.h> /* assert() */
2975 #include <stdlib.h> /* realloc(), free(), NULL */
2976 #include "global.h" /* environment_t */
2977 #include "k_model.h" /* k_func_t */
2978 #include "interp.h" /* interp_* */
2979 #include "accel_k.h"
2981 typedef struct accel_k_struct {
2982 interp_table_t *itable;
2988 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
2989 static int num_accels = 0;
2990 static accel_k_t *accels=NULL;
2992 /* Wrap k in the standard f(x) acceleration form */
2993 static double k_wrap(double F, void *params)
2995 accel_k_t *p = (accel_k_t *)params;
2996 return (*p->k)(F, p->env, p->params);
2999 static int k_tol(double FA, double kA, double FB, double kB)
3002 if (FB-FA > 1e-12) {
3003 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3004 return 1; /* unnacceptable */
3006 //printf("acceptable tol\n");
3007 return 0; /* acceptable */
3011 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3014 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3015 assert(accels != NULL);
3016 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3018 accels[i].env = env;
3019 accels[i].params = params;
3026 for (i=0; i<num_accels; i++)
3027 interp_table_free(accels[i].itable);
3029 if (accels) free(accels);
3033 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3036 for (i=0; i<num_accels; i++) {
3037 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3038 /* We've already setup this function.
3039 * WARNING: we're assuming param and env are the same. */
3040 return interp_table_eval(accels[i].itable, accels+i, F);
3043 /* set up a new acceleration environment */
3044 i = add_accel_k(k, env, params);
3045 return interp_table_eval(accels[i].itable, accels+i, F);
3049 \section{Tension models}
3050 \label{sec.tension_models}
3052 TODO: tension model intro.
3053 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.
3055 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3056 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]].
3058 <<tension-model.h>>=
3060 <<tension handler types>>
3061 <<constant tension model declarations>>
3062 <<hooke tension model declarations>>
3063 <<worm-like chain tension model declarations>>
3064 <<freely-jointed chain tension model declarations>>
3065 <<find tension definitions>>
3066 <<tension model global declarations>>
3069 <<tension model module makefile lines>>=
3070 NW_SAWSIM_MODS += tension_model
3071 #CHECK_BINS += check_tension_model
3072 check_tension_model_MODS =
3074 <<tension model utils makefile lines>>=
3075 TENSION_MODEL_MODS = tension_model parse list tension_balance
3076 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3077 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3078 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3079 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw $(BUILD_DIR)
3080 notangle -Rtension-model-utils.c $< > $@
3081 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) $(BIN_DIR)
3082 gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3083 clean_tension_model_utils :
3084 rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3088 \label{sec.null_tension}
3090 For unstretchable domains.
3092 <<null tension model getopt>>=
3093 {NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3096 \subsection{Constant}
3097 \label{sec.const_tension}
3099 <<constant tension functions>>=
3100 <<constant tension function>>
3101 <<constant tension parameter create/destroy functions>>
3104 <<constant tension model declarations>>=
3105 <<constant tension function declaration>>
3106 <<constant tension parameter create/destroy function declarations>>
3107 <<constant tension model global declarations>>
3111 An infinitely stretchable domain providing a constant tension.
3113 <<constant tension function declaration>>=
3114 extern double const_tension_handler(double x, void *pdata);
3116 <<constant tension function>>=
3117 double const_tension_handler(double x, void *pdata)
3119 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3120 list_t *list = data->group;
3125 assert(list != NULL); /* empty active group?! */
3126 F = ((const_tension_param_t *)list->d)->F;
3128 fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3130 while (list != NULL) {
3131 assert(((const_tension_param_t *)list->d)->F == F);
3139 <<constant tension structure>>=
3140 typedef struct const_tension_param_struct {
3141 double F; /* tension (force) in N */
3142 } const_tension_param_t;
3146 <<constant tension parameter create/destroy function declarations>>=
3147 extern void *string_create_const_tension_param_t(char **param_strings);
3148 extern void destroy_const_tension_param_t(void *p);
3151 <<constant tension parameter create/destroy functions>>=
3152 const_tension_param_t *create_const_tension_param_t(double F)
3154 const_tension_param_t *ret
3155 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3156 assert(ret != NULL);
3161 void *string_create_const_tension_param_t(char **param_strings)
3163 return create_const_tension_param_t(atof(param_strings[0]));
3166 void destroy_const_tension_param_t(void *p)
3174 <<constant tension model global declarations>>=
3175 extern int num_const_tension_params;
3176 extern char *const_tension_param_descriptions[];
3177 extern char const_tension_param_string[];
3180 <<constant tension model globals>>=
3181 int num_const_tension_params = 1;
3182 char *const_tension_param_descriptions[] = {"tension F, N"};
3183 char const_tension_param_string[] = "0";
3186 <<constant tension model getopt>>=
3187 {&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}
3193 <<hooke functions>>=
3195 <<hooke parameter create/destroy functions>>
3198 <<hooke tension model declarations>>=
3199 <<hooke tension function declaration>>
3200 <<hooke tension parameter create/destroy function declarations>>
3201 <<hooke tension model global declarations>>
3205 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
3206 The behavior of a series of springs $k_i$ in series is given by
3208 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
3210 For a simple proof, see Appendix \ref{app.math_hooke}.
3212 <<hooke tension function declaration>>=
3213 extern double hooke_handler(double x, void *pdata);
3217 double hooke_handler(double x, void *pdata)
3219 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3220 list_t *list = data->group;
3225 assert(list != NULL); /* empty active group?! */
3226 while (list != NULL) {
3227 assert( ((hooke_param_t *)list->d)->k > 0 );
3228 k += 1.0/ ((hooke_param_t *)list->d)->k;
3230 fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
3231 __FUNCTION__, ((hooke_param_t *)list->d)->k);
3237 fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N\n",
3238 __FUNCTION__, k, x, k*x);
3244 <<hooke structure>>=
3245 typedef struct hooke_param_struct {
3246 double k; /* spring constant in N/m */
3250 <<hooke tension parameter create/destroy function declarations>>=
3251 extern void *string_create_hooke_param_t(char **param_strings);
3252 extern void destroy_hooke_param_t(void *p);
3255 <<hooke parameter create/destroy functions>>=
3256 hooke_param_t *create_hooke_param_t(double k)
3258 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
3259 assert(ret != NULL);
3264 void *string_create_hooke_param_t(char **param_strings)
3266 return create_hooke_param_t(atof(param_strings[0]));
3269 void destroy_hooke_param_t(void *p)
3276 <<hooke tension model global declarations>>=
3277 extern int num_hooke_params;
3278 extern char *hooke_param_descriptions[];
3279 extern char hooke_param_string[];
3282 <<hooke tension model globals>>=
3283 int num_hooke_params = 1;
3284 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
3285 char hooke_param_string[]="0.05";
3288 <<hooke tension model getopt>>=
3289 {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}
3292 \subsection{Worm-like chain}
3295 We can model several unfolded domains with a single worm-like chain.
3296 <<worm-like chain functions>>=
3297 <<worm-like chain function>>
3298 <<worm-like chain handler>>
3299 <<worm-like chain create/destroy functions>>
3302 <<worm-like chain tension model declarations>>=
3303 <<worm-like chain handler declaration>>
3304 <<worm-like chain create/destroy function declarations>>
3305 <<worm-like chain tension model global declarations>>
3309 The combination of all unfolded domains is modeled as a worm like chain,
3310 with the force $F_{\text{WLC}}$ approximately given by
3312 F_{\text{WLC}} \approx \frac{k_B T}{p}
3313 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
3315 where $x$ is the distance between the fixed ends,
3316 $k_B$ is Boltzmann's constant,
3317 $T$ is the absolute temperature,
3318 $p$ is the persistence length, and
3319 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
3320 <<worm-like chain function>>=
3321 double wlc(double x, double T, double p, double L)
3323 double xL=0.0; /* uses global k_B */
3325 assert(T > 0); assert(p > 0); assert(L > 0);
3326 if (x >= L) return HUGE_VAL;
3327 xL = x/L; /* unitless */
3328 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
3329 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
3330 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
3333 This model assumes that each unfolded domain has the same persistence length.
3335 <<worm-like chain handler declaration>>=
3336 extern double wlc_handler(double x, void *pdata);
3339 <<worm-like chain handler>>=
3340 double wlc_handler(double x, void *pdata)
3342 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3343 list_t *list = data->group;
3347 assert(list != NULL); /* empty active group?! */
3348 p = ((wlc_param_t *)list->d)->p;
3349 while (list != NULL) {
3350 /* enforce consistent p */
3351 assert( ((wlc_param_t *)list->d)->p == p);
3352 L += ((wlc_param_t *)list->d)->L;
3354 fprintf(stderr, "%s: WLC section %g m long in series\n",
3355 __FUNCTION__, ((wlc_param_t *)list->d)->L);
3360 fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
3361 __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
3363 return wlc(x, data->env->T, p, L);
3367 <<worm-like chain structure>>=
3368 typedef struct wlc_param_struct {
3369 double p; /* WLC persistence length (m) */
3370 double L; /* WLC contour length (m) */
3374 <<worm-like chain create/destroy function declarations>>=
3375 extern void *string_create_wlc_param_t(char **param_strings);
3376 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
3379 <<worm-like chain create/destroy functions>>=
3380 wlc_param_t *create_wlc_param_t(double p, double L)
3382 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
3383 assert(ret != NULL);
3389 void *string_create_wlc_param_t(char **param_strings)
3391 return create_wlc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3394 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
3401 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.
3402 TODO (now we copy the string before parsing, but still do this for clarity).
3404 <<valid string write code>>=
3405 char string[] = "abc";
3408 <<valid string write code>>=
3409 char *string = "abc";
3413 <<worm-like chain tension model global declarations>>=
3414 extern int num_wlc_params;
3415 extern char *wlc_param_descriptions[];
3416 extern char wlc_param_string[];
3418 <<worm-like chain tension model globals>>=
3419 int num_wlc_params = 2;
3420 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
3421 char wlc_param_string[]="0.39e-9,28.4e-9";
3424 <<worm-like chain tension model getopt>>=
3425 {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}
3427 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
3429 \subsection{Freely-jointed chain}
3432 <<freely-jointed chain functions>>=
3433 <<freely-jointed chain function>>
3434 <<freely-jointed chain handler>>
3435 <<freely-jointed chain create/destroy functions>>
3438 <<freely-jointed chain tension model declarations>>=
3439 <<freely-jointed chain handler declaration>>
3440 <<freely-jointed chain create/destroy function declarations>>
3441 <<freely-jointed chain tension model global declarations>>
3445 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
3446 The tension of a chain of $N$ such links is given by
3448 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
3450 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}.
3451 <<freely-jointed chain function>>=
3452 double langevin(double x, void *params)
3455 return 1.0/tanh(x) - 1/x;
3457 <<one dimensional bracket declaration>>
3458 <<one dimensional solve declaration>>
3459 double inv_langevin(double x)
3461 double lb=0.0, ub=1.0, ret;
3462 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
3463 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
3467 double fjc(double x, double T, double l, int N)
3469 double L = l*(double)N;
3470 if (x == 0) return 0;
3471 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
3472 return k_B*T/l * inv_langevin(x/L);
3476 <<freely-jointed chain handler declaration>>=
3477 extern double fjc_handler(double x, void *pdata);
3479 <<freely-jointed chain handler>>=
3480 double fjc_handler(double x, void *pdata)
3482 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3483 list_t *list = data->group;
3488 assert(list != NULL); /* empty active group?! */
3489 l = ((fjc_param_t *)list->d)->l;
3490 while (list != NULL) {
3491 /* enforce consistent l */
3492 assert( ((fjc_param_t *)list->d)->l == l);
3493 N += ((fjc_param_t *)list->d)->N;
3495 fprintf(stderr, "%s: FJC section %d links long in series\n",
3496 __FUNCTION__, ((fjc_param_t *)list->d)->N);
3501 fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
3502 __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
3504 return fjc(x, data->env->T, l, N);
3508 <<freely-jointed chain structure>>=
3509 typedef struct fjc_param_struct {
3510 double l; /* FJC link length (m) */
3511 int N; /* FJC number of links */
3515 <<freely-jointed chain create/destroy function declarations>>=
3516 extern void *string_create_fjc_param_t(char **param_strings);
3517 extern void destroy_fjc_param_t(void *p);
3520 <<freely-jointed chain create/destroy functions>>=
3521 fjc_param_t *create_fjc_param_t(double l, double N)
3523 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
3524 assert(ret != NULL);
3530 void *string_create_fjc_param_t(char **param_strings)
3532 return create_fjc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3535 void destroy_fjc_param_t(void *p)
3542 <<freely-jointed chain tension model global declarations>>=
3543 extern int num_fjc_params;
3544 extern char *fjc_param_descriptions[];
3545 extern char fjc_param_string[];
3548 <<freely-jointed chain tension model globals>>=
3549 int num_fjc_params = 2;
3550 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
3551 char fjc_param_string[]="0.5e-9,200";
3554 <<freely-jointed chain tension model getopt>>=
3555 {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}
3558 \subsection{Tension model implementation}
3560 <<tension-model.c>>=
3563 <<tension model includes>>
3564 #include "tension_model.h"
3565 <<tension model internal definitions>>
3566 <<tension model globals>>
3567 <<tension model functions>>
3570 <<tension model includes>>=
3571 #include <assert.h> /* assert() */
3572 #include <stdlib.h> /* NULL */
3573 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
3574 #include <stdio.h> /* fprintf(), stdout */
3577 #include "tension_balance.h" /* oneD_*() */
3580 <<tension model internal definitions>>=
3581 <<constant tension structure>>
3583 <<worm-like chain structure>>
3584 <<freely-jointed chain structure>>
3587 <<tension model globals>>=
3588 <<tension handler array global>>
3589 <<constant tension model globals>>
3590 <<hooke tension model globals>>
3591 <<worm-like chain tension model globals>>
3592 <<freely-jointed chain tension model globals>>
3595 <<tension model functions>>=
3596 <<constant tension functions>>
3598 <<worm-like chain functions>>
3599 <<freely-jointed chain functions>>
3603 \subsection{Utilities}
3605 The tension models can get complicated, and users may want to reassure
3606 themselves that this portion of the input physics are functioning properly. The
3607 stand-alone program [[tension_model_utils]] demonstrates the tension model
3608 interface without the overhead of having to understand a full simulation.
3609 [[tension_model_utils]] takes command line model arguments like the full
3610 simulation, and prints $F(x)$ data to the screen for the selected model over a
3613 <<tension-model-utils.c>>=
3615 <<tension model utility includes>>
3616 <<tension model utility definitions>>
3617 <<tension model utility getopt functions>>
3619 <<tension model utility main>>
3622 <<tension model utility main>>=
3623 int main(int argc, char **argv)
3625 <<tension model getopt array definition>>
3626 tension_model_getopt_t *model = NULL;
3630 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
3631 tension_handler_data_t tdata;
3632 int num_param_args; /* for INIT_MODEL() */
3633 char **param_args; /* for INIT_MODEL() */
3634 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model, &flags);
3636 if (flags & VFLAG) {
3637 printf("#initializing model %s with parameters %s\n", model->name, model->params);
3639 INIT_MODEL("utility", model, model->params, params);
3641 push(&tdata.group, params, 1);
3643 tdata.persist = NULL;
3644 if (model->handler == NULL) {
3645 printf("No tension function for model %s\n", model->name);
3649 double dx=1e-10, x=0, F=0;
3650 printf("#F (N)\tk (%% pop. per s)\n");
3651 while (F >= 0 && F < 1e5 && x < 1e-6) {
3652 F = (*model->handler)(x, &tdata);
3653 printf("%g\t%g\n", x, F);
3657 params = pop(&tdata.group);
3658 if (model->destructor)
3659 (*model->destructor)(params);
3664 <<tension model utility includes>>=
3667 #include <string.h> /* strlen() */
3668 #include <assert.h> /* assert() */
3672 #include "tension_model.h"
3675 <<tension model utility definitions>>=
3676 <<version definition>>
3677 #define VFLAG 1 /* verbose */
3678 <<string comparison definition and globals>>
3679 <<tension model getopt definitions>>
3680 <<initialize model definition>>
3684 <<tension model utility getopt functions>>=
3685 <<index tension model>>
3686 <<tension model utility help>>
3687 <<tension model utility get options>>
3690 <<tension model utility help>>=
3691 void help(char *prog_name,
3693 int n_tension_models, tension_model_getopt_t *tension_models,
3697 printf("usage: %s [options]\n", prog_name);
3698 printf("Version %s\n\n", VERSION);
3699 printf("A utility for understanding the available tension models\n\n");
3700 printf("Environmental options:\n");
3701 printf("-T T\tTemperature T (currently %g K)\n", env->T);
3702 printf("-C T\tYou can also set the temperature T in Celsius\n");
3703 printf("Model options:\n");
3704 printf("-m model\tFolded tension model (currently %s)\n",
3705 tension_models[tension_model].name);
3706 printf("-a args\tFolded tension model argument string (currently %s)\n",
3707 tension_models[tension_model].params);
3708 printf("F(x) is calculated for a range of x and printed\n");
3709 printf("For example:\n");
3710 printf(" #Distance (x)\tForce (N)\n");
3711 printf(" 123.456\t7.89\n");
3713 printf("-V\tChange output to verbose mode\n");
3714 printf("-h\tPrint this help and exit\n");
3716 printf("Tension models:\n");
3717 for (i=0; i<n_tension_models; i++) {
3718 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
3719 for (j=0; j < tension_models[i].num_params ; j++)
3720 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
3721 printf("\t\tdefaults: %s\n", tension_models[i].params);
3726 <<tension model utility get options>>=
3727 void get_options(int argc, char **argv, environment_t *env,
3728 int n_tension_models, tension_model_getopt_t *tension_models,
3729 tension_model_getopt_t **model,
3730 unsigned int *flags)
3732 char *prog_name = NULL;
3733 char c, options[] = "T:C:m:a:Vh";
3734 int tension_model=0, num_strings;
3735 extern char *optarg;
3736 extern int optind, optopt, opterr;
3740 /* setup defaults */
3742 prog_name = argv[0];
3743 env->T = 300.0; /* K */
3745 *model = tension_models;
3747 while ((c=getopt(argc, argv, options)) != -1) {
3749 case 'T': env->T = atof(optarg); break;
3750 case 'C': env->T = atof(optarg)+273.15; break;
3752 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
3753 *model = tension_models+tension_model;
3756 tension_models[tension_model].params = optarg;
3758 case 'V': *flags |= VFLAG; break;
3760 fprintf(stderr, "unrecognized option '%c'\n", optopt);
3761 /* fall through to default case */
3763 help(prog_name, env, n_tension_models, tension_models, tension_model);
3772 \section{Unfolding rate models}
3773 \label{sec.k_models}
3775 We have two main choices for determining $k$: Bell's model, which assumes the
3776 folded domains exist in a single `folded' state and make instantaneous,
3777 stochastic transitions over a free energy barrier; and Kramers' model, which
3778 treats the unfolding as a thermalized diffusion process.
3779 We deal with these options like we dealt with the different tension models: by bundling all model-specific
3780 parameters into structures, with model specific functions for determining $k$.
3782 <<k func definition>>=
3783 typedef double k_func_t(double F, environment_t *env, void *params);
3786 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3787 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]].
3789 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
3790 because \LaTeX\ doesn't like underscores outside math mode.
3791 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
3792 You could use spaces instead (`\verb+ +'),
3793 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
3794 which I don't like as much.
3798 <<k func definition>>
3799 <<null k declarations>>
3800 <<const k declarations>>
3801 <<bell k declarations>>
3802 <<kramers k declarations>>
3803 <<kramers integ k declarations>>
3806 <<k model module makefile lines>>=
3807 NW_SAWSIM_MODS += k_model
3808 CHECK_BINS += check_k_model
3809 check_k_model_MODS = parse string_eval
3811 [[check_k_model]] also depends on the parse module.
3813 <<k model utils makefile lines>>=
3814 K_MODEL_MODS = k_model parse string_eval
3815 K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
3816 $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3817 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw $(BUILD)
3818 notangle -Rk-model-utils.c $< > $@
3819 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) $(BIN_DIR)
3820 gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3821 clean_k_model_utils :
3822 rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
3828 For domains that are never folded.
3830 <<null k declarations>>=
3832 <<null k functions>>=
3837 <<null k model getopt>>=
3838 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
3841 \subsection{Constant rate model}
3844 This is a pretty boring model, but it is usefull for testing the engine.
3848 where $k_0$ is the constant unfolding rate.
3850 <<const k functions>>=
3851 <<const k function>>
3852 <<const k structure create/destroy functions>>
3855 <<const k declarations>>=
3856 <<const k function declaration>>
3857 <<const k structure create/destroy function declarations>>
3858 <<const k global declarations>>
3862 <<const k structure>>=
3863 typedef struct const_k_param_struct {
3864 double knot; /* transition rate k_0 (frac population per s) */
3868 <<const k function declaration>>=
3869 double const_k(double F, environment_t *env, void *const_k_params);
3871 <<const k function>>=
3872 double const_k(double F, environment_t *env, void *const_k_params)
3873 { /* Returns the rate constant k in frac pop per s. */
3874 const_k_param_t *p = (const_k_param_t *) const_k_params;
3876 assert(p->knot > 0);
3881 <<const k structure create/destroy function declarations>>=
3882 void *string_create_const_k_param_t(char **param_strings);
3883 void destroy_const_k_param_t(void *p);
3886 <<const k structure create/destroy functions>>=
3887 const_k_param_t *create_const_k_param_t(double knot)
3889 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
3890 assert(ret != NULL);
3895 void *string_create_const_k_param_t(char **param_strings)
3897 return create_const_k_param_t(atof(param_strings[0]));
3900 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
3907 <<const k global declarations>>=
3908 extern int num_const_k_params;
3909 extern char *const_k_param_descriptions[];
3910 extern char const_k_param_string[];
3913 <<const k globals>>=
3914 int num_const_k_params = 1;
3915 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
3916 char const_k_param_string[]="1";
3919 <<const k model getopt>>=
3920 {"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}
3923 \subsection{Bell's model}
3926 Quantitatively, Bell's model gives $k$ as
3928 k = k_0 \cdot e^{F dx / k_B T} \;,
3930 where $k_0$ is the unforced unfolding rate,
3931 $F$ is the applied unfolding force,
3932 $dx$ is the distance to the transition state, and
3933 $k_B T$ is the thermal energy\citep{schlierf06}.
3935 <<bell k functions>>=
3937 <<bell k structure create/destroy functions>>
3940 <<bell k declarations>>=
3941 <<bell k function declaration>>
3942 <<bell k structure create/destroy function declarations>>
3943 <<bell k global declarations>>
3947 <<bell k structure>>=
3948 typedef struct bell_param_struct {
3949 double knot; /* transition rate k_0 (frac population per s) */
3950 double dx; /* `distance to transition state' (nm) */
3954 <<bell k function declaration>>=
3955 double bell_k(double F, environment_t *env, void *bell_params);
3957 <<bell k function>>=
3958 double bell_k(double F, environment_t *env, void *bell_params)
3959 { /* Returns the rate constant k in frac pop per s.
3960 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
3961 * uses global k_B in J/K */
3962 bell_param_t *p = (bell_param_t *) bell_params;
3963 assert(F >= 0); assert(env->T > 0);
3965 assert(p->knot > 0); assert(p->dx >= 0);
3967 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
3968 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
3969 p->knot * exp(F*p->dx / (k_B*env->T) ));
3970 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
3972 return p->knot * exp(F*p->dx / (k_B*env->T) );
3976 <<bell k structure create/destroy function declarations>>=
3977 void *string_create_bell_param_t(char **param_strings);
3978 void destroy_bell_param_t(void *p);
3981 <<bell k structure create/destroy functions>>=
3982 bell_param_t *create_bell_param_t(double knot, double dx)
3984 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
3985 assert(ret != NULL);
3991 void *string_create_bell_param_t(char **param_strings)
3993 return create_bell_param_t(atof(param_strings[0]), atof(param_strings[1]));
3996 void destroy_bell_param_t(void *p /* bell_param_t * */)
4003 <<bell k global declarations>>=
4004 extern int num_bell_params;
4005 extern char *bell_param_descriptions[];
4006 extern char bell_param_string[];
4010 int num_bell_params = 2;
4011 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4012 char bell_param_string[]="3.3e-4,0.25e-9";
4015 <<bell k model getopt>>=
4016 {"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}
4018 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4019 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4021 \subsection{Kramer's model}
4024 Kramer's model gives $k$ as
4026 % \frac{1}{k} = \frac{1}{D}
4027 % \int_{x_\text{min}}^{x_\text{max}}
4028 % e^{\frac{-U_F(x)}{k_B T}}
4029 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4032 %where $D$ is the diffusion constant,
4033 %$U_F(x)$ is the free energy along the unfolding coordinate, and
4034 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4035 % before the minimum and after the maximum, respectively \citep{schlierf06}.
4037 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
4039 where $D$ is the diffusion constant,
4041 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
4043 are length scales where
4044 $x_c(F)$ is the location of the bound state and
4045 $x_{ts}(F)$ is the location of the transition state,
4046 $E(F,x)$ is the free energy, and
4047 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
4048 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
4049 \citet{evans97} Eqn.~3).
4051 <<kramers k functions>>=
4052 <<kramers k function>>
4053 <<kramers k structure create/destroy functions>>
4056 <<kramers k declarations>>=
4057 <<kramers k function declaration>>
4058 <<kramers k structure create/destroy function declarations>>
4059 <<kramers k global declarations>>
4063 <<kramers k structure>>=
4064 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
4065 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
4067 typedef struct kramers_param_struct {
4068 double D; /* diffusion rate (um/s) */
4075 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
4076 //kramers_x_func_t *xts; /* function returning transition x (nm) */
4077 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
4078 //kramers_E_func_t *E; /* function returning E (J) */
4079 //double *E_params; /* parametrize the energy landscape */
4080 //int n_E_params; /* length of E_params array */
4084 <<kramers k function declaration>>=
4085 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
4086 extern double kramers_k(double F, environment_t *env, void *kramers_params);
4088 <<kramers k function>>=
4089 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
4091 static char input[160]={0};
4092 static char output[80]={0};
4093 /* setup the environment */
4094 fprintf(in, "F = %g; T = %g\n", F, T);
4095 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4096 string_eval(in, out, input, 80, output);
4097 return atof(output);
4100 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
4102 static char input[160]={0};
4103 static char output[80]={0};
4104 /* setup the environment */
4105 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
4106 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4107 string_eval(in, out, input, 80, output);
4108 return atof(output);
4111 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
4113 kramers_param_t *p = (kramers_param_t *) kramers_params;
4114 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
4117 double kramers_k(double F, environment_t *env, void *kramers_params)
4118 { /* Returns the rate constant k in frac pop per s.
4119 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4120 * uses global k_B in J/K */
4121 kramers_param_t *p = (kramers_param_t *) kramers_params;
4122 double kbT, xc, xts, lc, lts, Eb;
4123 assert(F >= 0); assert(env->T > 0);
4126 assert(p->xc != NULL); assert(p->xts != NULL);
4127 assert(p->ddEdxx != NULL); assert(p->E != NULL);
4128 assert(p->in != NULL); assert(p->out != NULL);
4130 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
4131 if (xc == -1.0) return -1.0; /* error (Too much force) */
4132 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
4133 if (xc == -1.0) return -1.0; /* error (Too much force) */
4134 lc = sqrt(2.0*M_PI*kbT /
4135 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
4136 lts = sqrt(-2.0*M_PI*kbT/
4137 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
4138 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
4139 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
4140 //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));
4141 return p->D/(lc*lts) * exp(-Eb/kbT);
4145 <<kramers k structure create/destroy function declarations>>=
4146 void *string_create_kramers_param_t(char **param_strings);
4147 void destroy_kramers_param_t(void *p);
4150 <<kramers k structure create/destroy functions>>=
4151 kramers_param_t *create_kramers_param_t(double D,
4152 char *xc_fn, char *xts_fn,
4153 char *E_fn, char *ddEdxx_fn,
4154 char *global_define_string)
4155 // kramers_x_func_t *xc, kramers_x_func_t *xts,
4156 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
4157 // double *E_params, int n_E_params)
4159 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
4160 assert(ret != NULL);
4161 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4162 assert(ret->xc != NULL);
4163 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4164 assert(ret->xts != NULL);
4165 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
4166 assert(ret->E != NULL);
4167 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
4168 assert(ret->ddEdxx != NULL);
4170 strcpy(ret->xc, xc_fn);
4171 strcpy(ret->xts, xts_fn);
4172 strcpy(ret->E, E_fn);
4173 strcpy(ret->ddEdxx, ddEdxx_fn);
4174 //ret->E_params = E_params;
4175 //ret->n_E_params = n_E_params;
4176 string_eval_setup(&ret->in, &ret->out);
4177 fprintf(ret->in, "kB = %g\n", k_B);
4178 fprintf(ret->in, "%s\n", global_define_string);
4182 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
4183 void *string_create_kramers_param_t(char **param_strings)
4185 return create_kramers_param_t(atof(param_strings[0]),
4193 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
4195 kramers_param_t *pk = (kramers_param_t *)p;
4197 string_eval_teardown(&pk->in, &pk->out);
4199 // free(pk->E_params);
4205 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
4206 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.
4207 We follow \citet{shillcock98} and use
4209 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
4210 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
4213 where TODO, variable meanings.%\citep{shillcock98}.
4214 The first and second derivatives of $E(F,E_b,x,x_s)$ are
4216 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\\
4217 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
4219 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
4220 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
4221 The roots are therefor at
4223 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
4224 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
4225 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
4228 <<kramers k global declarations>>=
4229 extern int num_kramers_params;
4230 extern char *kramers_param_descriptions[];
4231 extern char kramers_param_string[];
4234 <<kramers k globals>>=
4235 int num_kramers_params = 6;
4236 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"};
4237 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)}";
4239 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
4240 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
4241 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
4242 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.
4243 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
4244 It works so long as [[val_1]] is not [[false]].
4246 <<kramers k model getopt>>=
4247 {"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}
4250 \subsection{Kramer's model (natural cubic splines)}
4251 \label{sec.kramers_integ}
4253 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
4254 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
4255 \citet{schlierf06} Eqn.~2)
4257 \frac{1}{k} = \frac{1}{D}
4258 \int_{x_\text{min}}^{x_\text{max}}
4259 e^{\frac{U_F(x)}{k_B T}}
4260 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4263 where $D$ is the diffusion constant,
4264 $U_F(x)$ is the free energy along the unfolding coordinate, and
4265 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4266 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
4268 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
4269 interpolating between them with cubic splines.
4270 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
4272 <<kramers integ k functions>>=
4273 <<spline functions>>
4274 <<kramers integ A k functions>>
4275 <<kramers integ B k functions>>
4276 <<kramers integ k function>>
4277 <<kramers integ k structure create/destroy functions>>
4280 <<kramers integ k declarations>>=
4281 <<kramers integ k function declaration>>
4282 <<kramers integ k structure create/destroy function declarations>>
4283 <<kramers integ k global declarations>>
4287 <<kramers integ k structure>>=
4288 typedef double func_t(double x, void *params);
4289 typedef struct kramers_integ_param_struct {
4290 double D; /* diffusion rate (um/s) */
4291 double x_min; /* integration bounds */
4293 func_t *E; /* function returning E (J) */
4294 void *E_params; /* parametrize the energy landscape */
4295 destroy_data_func_t *destroy_E_params;
4297 double F; /* for passing into GSL evaluations */
4299 } kramers_integ_param_t;
4302 <<spline param structure>>=
4303 typedef struct spline_param_struct {
4304 int n_knots; /* natural cubic spline knots for U(x) */
4307 gsl_spline *spline; /* GSL spline parameters */
4308 gsl_interp_accel *acc; /* GSL spline acceleration data */
4312 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
4314 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4316 <<kramers integ A k functions>>=
4317 double kramers_integ_k_integrandA(double x, void *params)
4319 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4320 double U = (*p->E)(x, p->E_params);
4321 double Fx = p->F * x;
4322 double kbT = k_B * p->env->T;
4323 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
4324 return exp(-(U-Fx)/kbT);
4326 double kramers_integ_k_integralA(double x, void *params)
4328 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4330 double result, abserr;
4331 size_t neval; /* for qng */
4332 static gsl_integration_workspace *w = NULL;
4333 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4334 f.function = &kramers_integ_k_integrandA;
4336 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4337 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4338 //fprintf(stderr, "integralA = %g\n", result);
4344 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
4345 \int_{x_\text{min}}^{x_\text{max}}
4346 e^{\frac{U_F(x)}{k_B T}}
4347 \text{Integral}_A(x)
4350 <<kramers integ B k functions>>=
4351 double kramers_integ_k_integrandB(double x, void *params)
4353 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4354 double U = (*p->E)(x, p->E_params);
4355 double Fx = p->F * x;
4356 double kbT = k_B * p->env->T;
4357 double intA = kramers_integ_k_integralA(x, params);
4358 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
4359 return exp((U-Fx)/kbT)*intA;
4361 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
4363 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4365 double result, abserr;
4366 size_t neval; /* for qng */
4367 static gsl_integration_workspace *w = NULL;
4368 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4369 f.function = &kramers_integ_k_integrandB;
4373 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4374 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4375 //fprintf(stderr, "integralB = %g\n", result);
4380 With the integrals taken care of, evaluating $k$ is simply
4382 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
4384 <<kramers integ k function declaration>>=
4385 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
4387 <<kramers integ k function>>=
4388 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
4389 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
4390 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
4391 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
4395 <<kramers integ k structure create/destroy function declarations>>=
4396 void *string_create_kramers_integ_param_t(char **param_strings);
4397 void destroy_kramers_integ_param_t(void *p);
4400 <<kramers integ k structure create/destroy functions>>=
4401 kramers_integ_param_t *create_kramers_integ_param_t(double D,
4402 double xmin, double xmax, func_t *E,
4404 destroy_data_func_t *destroy_E_params)
4406 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
4407 assert(ret != NULL);
4412 ret->E_params = E_params;
4413 ret->destroy_E_params = destroy_E_params;
4417 void *string_create_kramers_integ_param_t(char **param_strings)
4419 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
4420 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
4421 void *E_params = string_create_spline_param_t(param_strings+1);
4422 return create_kramers_integ_param_t(atof(param_strings[0]),
4423 atof(param_strings[2]),
4424 atof(param_strings[3]),
4425 &spline_eval, E_params,
4426 destroy_spline_param_t);
4429 void destroy_kramers_integ_param_t(void *params)
4431 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4434 (*p->destroy_E_params)(p->E_params);
4440 Finally we have the GSL spline wrappers:
4441 <<spline functions>>=
4443 <<create/destroy spline>>
4446 <<spline function>>=
4447 double spline_eval(double x, void *spline_params)
4449 spline_param_t *p = (spline_param_t *)(spline_params);
4450 return gsl_spline_eval(p->spline, x, p->acc);
4454 <<create/destroy spline>>=
4455 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
4457 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
4458 assert(ret != NULL);
4459 ret->n_knots = n_knots;
4462 ret->acc = gsl_interp_accel_alloc();
4463 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
4464 gsl_spline_init(ret->spline, x, y, n_knots);
4468 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
4469 void *string_create_spline_param_t(char **param_strings)
4473 double *x=NULL, *y=NULL;
4474 /* break into ordered pairs */
4475 parse_list_string(param_strings[0], SEP, '(', ')',
4476 &num_knots, &knot_coords);
4477 x = (double *)malloc(sizeof(double)*num_knots);
4479 y = (double *)malloc(sizeof(double)*num_knots);
4481 for (i=0; i<num_knots; i++) {
4482 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
4483 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
4486 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
4488 free_parsed_list(num_knots, knot_coords);
4489 return create_spline_param_t(num_knots, x, y);
4492 void destroy_spline_param_t(void *params)
4494 spline_param_t *p = (spline_param_t *)params;
4497 gsl_spline_free(p->spline);
4499 gsl_interp_accel_free(p->acc);
4509 <<kramers integ k global declarations>>=
4510 extern int num_kramers_integ_params;
4511 extern char *kramers_integ_param_descriptions[];
4512 extern char kramers_integ_param_string[];
4515 <<kramers integ k globals>>=
4516 int num_kramers_integ_params = 4;
4517 char *kramers_integ_param_descriptions[] = {
4518 "diffusion D, m^2 / s",
4519 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
4520 "starting integration bound x_min, m",
4521 "ending integration bound x_max, m"};
4522 //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";
4523 //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";
4524 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";
4525 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
4526 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
4527 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
4530 <<kramers integ k model getopt>>=
4531 {"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}
4533 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4534 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4536 \subsection{Unfolding model implementation}
4540 <<k model includes>>
4541 #include "k_model.h"
4542 <<k model internal definitions>>
4544 <<k model functions>>
4547 <<k model includes>>=
4548 #include <assert.h> /* assert() */
4549 #include <stdlib.h> /* NULL, malloc() */
4550 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4551 #include <stdio.h> /* fprintf(), stdout */
4552 #include <string.h> /* strlen(), strcpy() */
4553 #include <gsl/gsl_integration.h>
4554 #include <gsl/gsl_spline.h>
4559 <<k model internal definitions>>=
4560 <<const k structure>>
4561 <<bell k structure>>
4562 <<kramers k structure>>
4563 <<kramers integ k structure>>
4564 <<spline param structure>>
4567 <<k model globals>>=
4571 <<kramers k globals>>
4572 <<kramers integ k globals>>
4575 <<k model functions>>=
4576 <<null k functions>>
4577 <<const k functions>>
4578 <<bell k functions>>
4579 <<kramers k functions>>
4580 <<kramers integ k functions>>
4583 \subsection{Unfolding model unit tests}
4585 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
4586 <<check-k-model.c>>=
4587 <<k model check includes>>
4588 <<check relative error>>
4590 <<k model test suite>>
4591 <<main check program>>
4594 <<k model check includes>>=
4595 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
4596 #include <stdio.h> /* sprintf() */
4597 #include <assert.h> /* assert() */
4598 #include <math.h> /* exp() */
4601 #include "k_model.h"
4604 <<k model test suite>>=
4607 <<k model suite definition>>
4610 <<k model suite definition>>=
4611 Suite *test_suite (void)
4613 Suite *s = suite_create ("k model");
4614 <<const k test case defs>>
4615 <<bell k test case defs>>
4617 <<const k test case adds>>
4618 <<bell k test case adds>>
4623 \subsubsection{Constant}
4625 <<const k test case defs>>=
4626 TCase *tc_const_k = tcase_create("Constant k");
4629 <<const k test case adds>>=
4630 tcase_add_test(tc_const_k, test_const_k_create_destroy);
4631 tcase_add_test(tc_const_k, test_const_k_over_range);
4632 suite_add_tcase(s, tc_const_k);
4637 START_TEST(test_const_k_create_destroy)
4639 char *knot[] = {"1","2","3","4","5","6"};
4640 char *params[] = {knot[0]};
4643 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4644 params[0] = knot[i];
4645 p = string_create_const_k_param_t(params);
4646 destroy_const_k_param_t(p);
4651 START_TEST(test_const_k_over_range)
4653 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
4654 char *knot[] = {"1","2","3","4","5","6"};
4655 char *params[] = {knot[0]};
4658 char param_string[80];
4661 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4662 params[0] = knot[i];
4663 p = string_create_const_k_param_t(params);
4664 for ( F=Fm; F<FM; F+=dF ) {
4665 fail_unless(const_k(F, &env, p)==atof(knot[i]),
4666 "const_k(%g, %g, &{%s}) = %g != %s",
4667 F, T, knot[i], const_k(F, &env, p), knot[i]);
4669 destroy_const_k_param_t(p);
4675 \subsubsection{Bell}
4677 <<bell k test case defs>>=
4678 TCase *tc_bell = tcase_create("Bell's k");
4681 <<bell k test case adds>>=
4682 tcase_add_test(tc_bell, test_bell_k_create_destroy);
4683 tcase_add_test(tc_bell, test_bell_k_at_zero);
4684 tcase_add_test(tc_bell, test_bell_k_at_one);
4685 suite_add_tcase(s, tc_bell);
4689 START_TEST(test_bell_k_create_destroy)
4691 char *knot[] = {"1","2","3","4","5","6"};
4693 char *params[] = {knot[0], dx};
4696 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4697 params[0] = knot[i];
4698 p = string_create_bell_param_t(params);
4699 destroy_bell_param_t(p);
4704 START_TEST(test_bell_k_at_zero)
4706 double F=0.0, T=300.0;
4707 char *knot[] = {"1","2","3","4","5","6"};
4709 char *params[] = {knot[0], dx};
4712 char param_string[80];
4715 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4716 params[0] = knot[i];
4717 p = string_create_bell_param_t(params);
4718 fail_unless(bell_k(F, &env, p)==atof(knot[i]),
4719 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
4720 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
4721 destroy_bell_param_t(p);
4726 START_TEST(test_bell_k_at_one)
4729 char *knot[] = {"1","2","3","4","5","6"};
4731 char *params[] = {knot[0], dx};
4732 double F= k_B*T/atof(dx);
4735 char param_string[80];
4738 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4739 params[0] = knot[i];
4740 p = string_create_bell_param_t(params);
4741 CHECK_ERR(1e-6, atof(knot[i])*exp(1.0), bell_k(F, &env, p));
4742 //fail_unless(bell_k(F, &env, p)==atof(knot[i])*exp(1.0),
4743 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
4744 // F, T, knot[i], dx, bell_k(F, &env, p), atof(knot[i])*exp(1.0));
4745 destroy_bell_param_t(p);
4751 <<kramers k tests>>=
4754 <<kramers k test case def>>=
4757 <<kramers k test case add>>=
4760 <<k model function tests>>=
4761 <<check relative error>>
4763 START_TEST(test_something)
4765 double k=5, x=3, last_x=2.0, F;
4766 one_dim_fn_t *handlers[] = {&hooke};
4767 void *data[] = {&k};
4768 double xi[] = {0.0};
4770 int new_active_groups = 1;
4771 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
4772 fail_unless(F = k*x, NULL);
4777 \subsection{Utilities}
4779 The unfolding models can get complicated, and are sometimes hard to explain to others.
4780 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
4781 the overhead of having to understand a full simulation.
4782 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
4783 or special mode, where the operation depends on the specific model selected.
4785 <<k-model-utils.c>>=
4787 <<k model utility includes>>
4788 <<k model utility definitions>>
4789 <<k model utility getopt functions>>
4790 <<k model utility multi model E>>
4791 <<k model utility main>>
4794 <<k model utility main>>=
4795 int main(int argc, char **argv)
4797 <<k model getopt array definition>>
4798 k_model_getopt_t *model = NULL;
4803 int num_param_args; /* for INIT_MODEL() */
4804 char **param_args; /* for INIT_MODEL() */
4806 double special_xmin;
4807 double special_xmax;
4808 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
4809 &Fmax, &special_xmin, &special_xmax, &flags);
4810 if (flags & VFLAG) {
4811 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4813 INIT_MODEL("utility", model, model->params, params);
4816 if (model->k == NULL) {
4817 printf("No k function for model %s\n", model->name);
4821 printf("Fmax = %g <= 0\n", Fmax);
4827 printf("#F (N)\tk (%% pop. per s)\n");
4828 for (i=0; i<=N; i++) {
4829 F = Fmax*i/(double)N;
4830 k = (*model->k)(F, &env, params);
4832 printf("%g\t%g\n", F, k);
4837 if (model->k == NULL || model->k == &bell_k) {
4838 printf("No special mode for model %s\n", model->name);
4841 if (special_xmax <= special_xmin) {
4842 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
4846 double Xrng=(special_xmax-special_xmin), x, E;
4848 printf("#x (m)\tE (J)\n");
4849 for (i=0; i<=N; i++) {
4850 x = special_xmin + Xrng*i/(double)N;
4851 E = multi_model_E(model->k, params, &env, x);
4852 printf("%g\t%g\n", x, E);
4859 if (model->destructor)
4860 (*model->destructor)(params);
4865 <<k model utility multi model E>>=
4866 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
4868 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
4870 if (k == kramers_integ_k)
4871 E = (*p->E)(x, p->E_params);
4873 E = kramers_E(0, env, model_params, x);
4879 <<k model utility includes>>=
4882 #include <string.h> /* strlen() */
4883 #include <assert.h> /* assert() */
4886 #include "string_eval.h"
4887 #include "k_model.h"
4890 <<k model utility definitions>>=
4891 <<version definition>>
4892 #define VFLAG 1 /* verbose */
4893 enum mode_t {M_K_OF_F, M_SPECIAL};
4894 <<string comparison definition and globals>>
4895 <<k model getopt definitions>>
4896 <<initialize model definition>>
4897 <<kramers k structure>>
4898 <<kramers integ k structure>>
4902 <<k model utility getopt functions>>=
4904 <<k model utility help>>
4905 <<k model utility get options>>
4908 <<k model utility help>>=
4909 void help(char *prog_name,
4911 int n_k_models, k_model_getopt_t *k_models,
4915 printf("usage: %s [options]\n", prog_name);
4916 printf("Version %s\n\n", VERSION);
4917 printf("A utility for understanding the available unfolding models\n\n");
4918 printf("Environmental options:\n");
4919 printf("-T T\tTemperature T (currently %g K)\n", env->T);
4920 printf("-C T\tYou can also set the temperature T in Celsius\n");
4921 printf("Model options:\n");
4922 printf("-k model\tTransition rate model (currently %s)\n",
4923 k_models[k_model].name);
4924 printf("-K args\tTransition rate model argument string (currently %s)\n",
4925 k_models[k_model].params);
4926 printf("There are two output modes. In standard mode, k(F) is printed\n");
4927 printf("For example:\n");
4928 printf(" #Force (N)\tk (% pop. per s)\n");
4929 printf(" 123.456\t7.89\n");
4931 printf("In special mode, the output depends on the model.\n");
4932 printf("For models defining an energy function E(x), that function is printed\n");
4933 printf(" #Position (m)\tE (J)\n");
4934 printf(" 0.0012\t0.0034\n");
4936 printf("-m\tChange output to standard mode\n");
4937 printf("-M\tChange output to special mode\n");
4938 printf("-F\tSet the maximum F value for the standard mode k(F) output\n");
4939 printf("-x\tSet the minimum x value for the special mode E(x) output\n");
4940 printf("-X\tSet the maximum x value for the special mode E(x) output\n");
4941 printf("-V\tChange output to verbose mode\n");
4942 printf("-h\tPrint this help and exit\n");
4944 printf("Unfolding rate models:\n");
4945 for (i=0; i<n_k_models; i++) {
4946 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
4947 for (j=0; j < k_models[i].num_params ; j++)
4948 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
4949 printf("\t\tdefaults: %s\n", k_models[i].params);
4954 <<k model utility get options>>=
4955 void get_options(int argc, char **argv, environment_t *env,
4956 int n_k_models, k_model_getopt_t *k_models,
4957 enum mode_t *mode, k_model_getopt_t **model,
4958 double *Fmax, double *special_xmin, double *special_xmax,
4959 unsigned int *flags)
4961 char *prog_name = NULL;
4962 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
4964 extern char *optarg;
4965 extern int optind, optopt, opterr;
4969 /* setup defaults */
4971 prog_name = argv[0];
4972 env->T = 300.0; /* K */
4977 *special_xmax = 1e-8;
4978 *special_xmin = 0.0;
4980 while ((c=getopt(argc, argv, options)) != -1) {
4982 case 'T': env->T = atof(optarg); break;
4983 case 'C': env->T = atof(optarg)+273.15; break;
4985 k_model = index_k_model(n_k_models, k_models, optarg);
4986 *model = k_models+k_model;
4989 k_models[k_model].params = optarg;
4991 case 'm': *mode = M_K_OF_F; break;
4992 case 'M': *mode = M_SPECIAL; break;
4993 case 'F': *Fmax = atof(optarg); break;
4994 case 'x': *special_xmin = atof(optarg); break;
4995 case 'X': *special_xmax = atof(optarg); break;
4996 case 'V': *flags |= VFLAG; break;
4998 fprintf(stderr, "unrecognized option '%c'\n", optopt);
4999 /* fall through to default case */
5001 help(prog_name, env, n_k_models, k_models, k_model);
5010 \section{\LaTeX\ documentation}
5012 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
5013 The comment blocks are extracted (with nicely formatted code blocks), using
5014 <<latex makefile lines>>=
5015 $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw $(BUILD_DIR)
5016 noweave -latex -index -delay $< > $@
5017 $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib $(BUILD_DIR)
5021 <<latex makefile lines>>=
5022 $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
5024 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5025 $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
5026 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5027 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5028 mv $(BUILD_DIR)/sawsim.pdf $@
5030 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
5031 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
5032 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
5034 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
5035 <<latex makefile lines>>=
5037 rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
5038 $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
5039 $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.tex \
5040 $(BUILD_DIR)/sawsim.bib
5041 clean_latex : semi-clean_latex
5042 rm -f $(DOC_DIR)/sawsim.pdf
5047 [[make]] is a common utility on *nix systems for managing dependencies while building software.
5048 In this case, we have several \emph{targets} that we'd like to build:
5049 the main simulation program \prog;
5050 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
5051 the documentation [[sawsim.pdf]];
5052 and the [[Makefile]] itself.
5053 Besides the generated files, there is the utility target
5054 [[clean]] for removing all generated files except [[Makefile]].
5056 % [[dist]] for generating a distributable tar file.
5058 Extract the makefile with
5059 `[[notangle -Rmakefile src/sawsim.nw | sed 's/ /\t/' > Makefile]]'.
5060 The sed is needed because notangle replaces tabs with eight spaces,
5061 but [[make]] needs tabs.
5063 # Decide what directories to put things in
5068 # Note: Cannot use, for example, `./build', because make eats the `./'
5069 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5071 # Modules (X.c, X.h) defined in the noweb file
5074 # Modules defined outside the noweb file
5075 FREE_SAWSIM_MODS = interp tavl
5077 <<list module makefile lines>>
5078 <<tension balance module makefile lines>>
5079 <<tension model module makefile lines>>
5080 <<k model module makefile lines>>
5081 <<parse module makefile lines>>
5082 <<string eval module makefile lines>>
5083 <<accel k module makefile lines>>
5085 SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS)
5087 # Everything needed for sawsim under one roof. sawsim.c must come first
5088 SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \
5089 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h)
5090 # Libraries needed to compile sawsim
5091 LIBS = gsl gslcblas m
5092 CHECK_LIBS = gsl gslcblas m check
5093 # The non-check binaries generated
5094 BINS = sawsim tension_model_utils k_model_utils sawsim_profile
5097 # Define the major targets
5098 all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
5100 view : $(DOC_DIR)/sawsim.pdf
5102 profile : $(BIN_DIR)/sawsim_profile $(BIN_DIR)
5103 $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \
5104 -mnull -Mwlc -A0.39e-9,28e-9 -F8
5105 gprof $(BIN_DIR)/sawsim_profile gmon.out > $@
5106 check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim
5107 $(SHELL) -e -c 'for B in $^; do ./$$B; done'
5108 clean : $(CHECK_BINS:%=clean_%) clean_tension_model_utils \
5109 clean_k_model_utils clean_latex
5110 rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \
5111 $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \
5112 $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
5113 $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
5114 $(BUILD_DIR)/global.h ./gmon.out
5115 $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
5117 # Various builds of sawsim
5118 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) $(BIN_DIR)
5119 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5120 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) $(BIN_DIR)
5121 gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5122 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) $(BIN_DIR)
5123 gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5125 # Create the directories
5126 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
5129 # Copy over the source external to sawsim
5130 # Note: Cannot use, for example, `./build', because make eats the `./'
5131 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5133 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
5137 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
5141 ## Basic source generated with noweb
5142 # The central files sawsim.c and global.h...
5143 $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw $(BUILD_DIR)
5145 $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw $(BUILD_DIR)
5146 notangle -Rglobal.h $< > $@
5147 # ... and the modules
5148 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw $(BUILD_DIR)
5149 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5150 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw $(BUILD_DIR)
5151 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5152 $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw $(BUILD_DIR)
5153 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5154 # Note: I use `_' as a space in my file names, but noweb doesn't like
5155 # that so we use `-' to name the noweb roots and substitute here.
5157 ## Building the unit-test programs
5159 $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
5160 notangle -Rchecks $< > $@
5161 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
5162 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
5163 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) $(BIN_DIR)
5164 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
5165 clean_check_sawsim :
5166 rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
5167 # ... and the modules
5169 $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
5170 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \
5171 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \
5172 $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
5173 $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
5174 $$(BUILD_DIR)/global.h $(BIN_DIR)
5175 gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
5176 $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
5178 # todo: clean up the required modules to
5179 $(CHECK_BINS:%=clean_%) :
5180 rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
5182 # Cleaning up the modules
5184 $(SAWSIM_MODS:%=clean_%) :
5185 rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h)
5187 <<tension model utils makefile lines>>
5188 <<k model utils makefile lines>>
5189 <<latex makefile lines>>
5191 Makefile : $(SRC_DIR)/sawsim.nw
5192 notangle -Rmakefile $< | sed 's/ /\t/' > $@
5194 This is fairly self-explanatory. We compile a dynamically linked
5195 version ([[sawsim]]) and a statically linked version
5196 ([[sawsim_static]]). The static version is about 9 times as big, but
5197 it runs without needing dynamic GSL libraries to link to, so it's a
5198 better format for distributing to a cluster for parallel evaluation.
5202 \subsection{Hookean springs in series}
5203 \label{app.math_hooke}
5206 The effective spring constant for $n$ springs $k_i$ in series is given by
5208 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5214 F &= k_i x_i &&\forall i \le n \\
5215 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
5216 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
5217 F &= k_1 x_1 = k_\text{eff} x \\
5218 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
5219 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5224 \addcontentsline{toc}{section}{References}
5225 \bibliography{sawsim}