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 {
232 <<environment definition>>
233 <<one dimensional function definition>>
234 <<create/destroy definitions>>
236 #endif /* GLOBAL_H */
238 Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
239 included multiple times.
241 \section{Simulation functions}
243 <<simulation functions>>=
247 <<does domain unfold>>
248 <<randomly unfold domains>>
252 \label{sec.find_tension}
254 Because the stretched system may be made up of several parts (folded
255 domains, unfolded domains, spring-like cantilever, \ldots), we will
256 assign the domains to models and groups. The models are used to
257 determine a domain's tension (Hookean spring, WLC, \ldots). Groups
258 will represent collections of domains which the model should treat as
259 a single entity. A domain may belong to different groups or models
260 depending on its state. For example, a domain might be modeled as a
261 freely-jointed chain when it iss folded and as a worm-like chain class
262 when it is unfolded. The domains are assumed to be commutative, so
263 ordering is ignored. The interactions between the groups are assumed
264 to be linear, but the interactions between domains of the same group
265 need not be. This allows for non-linear group models such as th
266 worm-like or freely-jointed chains. Each model has a tension handler
267 function, which gives the tension $F$ needed to stretch a given group
268 of domains a distance $x$.
270 To understand the purpose of groups, consider a chain of proteins
271 where two folded proteins $A$ and $B$ are modeled as WLCs with
272 persistence length $p_f$ and two unfolded proteins $C$ and $D$ are
273 modeled as WLCs with persistence length $p_u$. The proteins are
274 attached to a cantilever $E$ qof spring constant $κ$. The string
275 would get split into three lists:
277 \begin{tabular}{llll}
278 Model & Group & List & Tension \\
279 WLC & 0 & $[A,B]$ & $F_\text{WLC}(x, p_f, L_A+L_B)$ \\
280 WLC & 1 & $[C,D]$ & $F_\text{WLC}(x, p_u, L_C+L_D)$ \\
281 Hooke & 0 & $[E]$ & $F_\text{Hooke}(x, κ)$ \\
284 Note that group numbers only matter \emph{within} model classes, since
285 grouping domains with seperate models doesn't make sense.
287 <<find tension definitions>>=
288 #define NUM_TENSION_MODELS 5
292 <<tension handler array global declaration>>=
293 extern one_dim_fn_t *tension_handlers[];
296 <<tension handler array global>>=
297 one_dim_fn_t *tension_handlers[] = {
299 &const_tension_handler,
307 <<tension model global declarations>>=
308 <<tension handler array global declaration>>
311 <<tension handler types>>=
312 typedef struct tension_handler_data_struct {
313 /* int num_domains; */
314 /* domain_t *domains; */
318 } tension_handler_data_t;
322 After sorting the chain into separate groups $G_i$, with tension
323 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
324 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
325 \forall i,j$. Note that there may be several groups within each model
326 class. [[find_tension]] is basically a wrapper to feed proper input
327 into the [[tension_balance]] function. [[unfolding]] should be set to
328 0 if there were no unfoldings since the previous call to
331 double find_tension(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list, environment_t *env, double x, int unfolding)
333 static int num_active_groups=0;
334 static one_dim_fn_t **per_group_handlers = NULL;
335 static void **per_group_params = NULL;
336 static double last_x;
337 tension_handler_data_t *tdata;
342 fprintf(stderr, "finding tension at distance %g%s\n", x, unfolding ? " after an unfolding" : "");
343 list_print(stderr, domain_list, "domain list");
346 if (unfolding != 0 || num_active_groups == 0) { /* setup statics */
347 /* free old statics */
348 if (per_group_handlers != NULL)
349 free(per_group_handlers);
350 if (per_group_params != NULL) {
351 for (i=0; i < num_active_groups; i++) {
352 assert(per_group_params[i] != NULL);
353 free(per_group_params[i]);
355 free(per_group_params);
358 /* get new statics */
359 get_active_groups(domain_list, num_tension_handlers, tension_handlers, &num_active_groups, &per_group_handlers, &per_group_params);
361 /* fill in the group handlers and params */
362 for (i=0; i<num_active_groups; i++) {
363 tdata = (tension_handler_data_t *) per_group_params[i];
365 fprintf(stderr, "active group %d of %d", i, num_active_groups);
366 list_print(stderr, tdata->group, " parameters");
369 /* tdata->persist continues unchanged */
372 } /* else roll with the current statics */
374 F = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_params, last_x, x);
378 @ For the moment, we will restrict our group boundaries to a single
379 dimension, so $\sum_i x_i = x$, our total extension (it is this
380 restriction that causes the groups to interact linearly). We'll also
381 restrict our extensions to all be positive. With these restrictions,
382 the problem of balancing the tensions reduces to solving a system of
383 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
384 the number of active groups. In general this can be fairly
385 complicated, but without much loss of practicality we can restrict
386 ourselves to strictly monotonically increasing, non-negative tension
387 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
388 simpler. For example, it guarantees the existence of a unique, real
389 solution for finite forces. See Appendix \ref{app.tension_balance}
390 for the tension balancing implementation.
392 Each group must define a parameter structure if the tension model is
394 a function to create the parameter structure,
395 a function to destroy the parameter structure, and
396 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as it's data argument.
397 The tension-specific parameter structures are contained in the domain
398 group field. See the Section \ref{app.model_selection} for a
399 discussion on the command line framework. See the worm-like chain
400 implementation for an example (Section \ref{sec.wlc}).
402 The major limitation of our approach is that all of our tension models
403 are Markovian (memory-less), which is not really a problem for the
404 chains (see \citet{evans99} for WLC thermalization rate calculations).
406 \subsection{Unfolding rate}
407 \label{sec.unfolding_rate}
409 Each folded domain is modeled as unimolecular, first order reaction
411 \text{Folded} \xrightarrow{k} % in tex: X atop Y
414 With the rate constant $k$ defined by
416 \frac{d[\text{Folded}]}{dt} = -k [\text{Folded}].
418 So the probability of a given protein unfolding in a short time $dt$
424 <<does domain unfold>>=
425 int domain_unfolds(double F, double dt, environment_t *env, domain_t *domain)
426 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to a folded domain */
428 k = accel_k(domain->k, F, env, domain->k_params);
429 //(*domain->k)(F, env, domain->k_params);
430 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
431 return happens(k*dt); /* dice roll for prob. k*dt event */
433 @ [[happens]] is a random decision making function defined in Appendix \ref{app.utils}.
435 Only one domain can unfold in each timestep, because the timescale of
436 a domain unfolding $dt_u$ is assumed to be less than the simulation
437 timestep $dt$, so a domain will completely unfold in a single
438 timestep. We adapt our timesteps to keep the probability of a single
439 domain unfolding low, so the probability of two domains unfolding in
440 the same timestep is negligible. This approach breaks down as the
441 adaptive timestep scheme approaches $dt \sim dt_u$, but $dt_u \sim
442 1\U{$\mu$s}$ for Ig27-like proteins \citep{klimov00}, so this
443 shouldn't be a problem. To reassure yourself, you can ask the
444 simulator to print the smallest timestep that was used with TODO.
445 <<randomly unfold domains>>=
446 int random_unfoldings(list_t *domain_list, double tension,
447 double dt, environment_t *env,
448 int *num_folded_domains)
449 { /* returns 1 if there was an unfolding and 0 otherwise */
450 while (domain_list != NULL) {
451 if (D(domain_list)->state == DS_FOLDED
452 && domain_unfolds(tension, dt, env, D(domain_list))) {
453 if (flags & FLAG_OUTPUT_UNFOLDING_FORCES)
454 fprintf(stdout, "%g\n", tension);
455 D(domain_list)->state = DS_UNFOLDED;
456 (*num_folded_domains)--;
457 return 1; /* our one domain has unfolded, stop looking for others */
459 domain_list = domain_list->next;
465 \subsection{Adaptive timesteps}
466 \label{sec.adaptive_dt}
468 We'd like to pick $dt$ so the probability of unfolding in the next
469 timestep is small. If we don't adapt our timestep, we risk breaking
470 our approximation that $F(x' \in [x, x+v \cdot dt]) = F(x)$ or that
471 only one domain unfolds in a given timestep. Because $F(x)$ is
472 monotonically increasing, excessively large timesteps will lead to
473 erroneously large unfolding forces. The simulation would have been
474 accurate for sufficiently small fixed timesteps, but adaptive
475 timesteps will allow us to move through low tension regions in fewer
476 steps, leading to a more efficient simulation.
478 The actual adaptive timestep implementation is not particularly
479 interesting, since we are only required to reduce $dt$ to somewhere
480 below a set threshold, so I've removed it to Appendix
481 \ref{app.adaptive_dt}.
487 The models provide the physics of the simulation, but the simulation
488 allows interchangeable models, and we are currently focusing on the
489 simulation itself, so we remove the actual model implementation to the
492 The tension models are defined in Appendix \ref{sec.tension_models}
493 and unfolding models are defined in Appendix \ref{sec.k_models}.
496 #define k_B 1.3806503e-23 /* J/K */
500 \section{Command line interface}
502 <<get option functions>>=
504 <<index tension model>>
510 \subsection{Model selection}
511 \label{app.model_selection}
513 The main difficulty with the command line interface in \prog\ is
514 developing an intuitive method for accessing the possibly large number
515 of available models. We'll treat this generally by defining an array
516 of available models, containing their important parameters.
518 <<get options definitions>>=
519 <<model getopt definitions>>
522 <<create/destroy definitions>>=
523 typedef void *create_data_func_t(char **param_strings);
524 typedef void destroy_data_func_t(void *param_struct);
527 <<model getopt definitions>>=
528 <<tension model getopt definitions>>
529 <<k model getopt definitions>>
533 \subsubsection{Tension}
535 To access the various tension models from the command line, we define
536 a structure that contains all the neccessary information about the
538 <<tension model getopt definitions>>=
539 typedef struct tension_model_getopt_struct {
540 one_dim_fn_t *handler; /* fn implementing the model on a group */
541 char *name; /* name string identifying the model */
542 char *description; /* a brief description of the model */
543 int num_params; /* number of extra params for the model */
544 char **param_descriptions; /* descriptions of extra params */
545 char *params; /* default values for the extra params */
546 create_data_func_t *creator; /* param-string -> model-data-struct fn */
547 destroy_data_func_t *destructor; /* fn to free the model data structure */
548 } tension_model_getopt_t;
551 <<tension model getopt array definition>>=
552 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
553 <<null tension model getopt>>,
554 <<constant tension model getopt>>,
555 <<hooke tension model getopt>>,
556 <<worm-like chain tension model getopt>>,
557 <<freely-jointed chain tension model getopt>>
561 \subsubsection{Unfolding rate}
563 <<k model getopt definitions>>=
564 #define NUM_K_MODELS 5
566 typedef struct k_model_getopt_struct {
571 char **param_descriptions;
573 create_data_func_t *creator;
574 destroy_data_func_t *destructor;
578 <<k model getopt array definition>>=
579 k_model_getopt_t k_models[NUM_K_MODELS] = {
580 <<null k model getopt>>,
581 <<const k model getopt>>,
582 <<bell k model getopt>>,
583 <<kramers k model getopt>>,
584 <<kramers integ k model getopt>>
591 void help(char *prog_name, double P, double dt_max, double v,
593 int n_tension_models, tension_model_getopt_t *tension_models,
594 int folded_tension_model, int unfolded_tension_model,
595 int n_k_models, k_model_getopt_t *k_models,
599 printf("usage: %s [options]\n", prog_name);
600 printf("Version %s\n\n", VERSION);
601 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
602 printf("Simulation options:\n");
603 printf("-P P\tTarget probability for dt (currently %g)\n", P);
604 printf("-t dt\tMaximum allowed timestep dt (currently %g)\n", dt_max);
605 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
606 printf("Environmental options:\n");
607 printf("-T T\tTemperature T (currently %g K)\n", env->T);
608 printf("-C T\tYou can also set the temperature T in Celsius\n");
609 printf("Model options:\n");
610 printf("The domains exist in either the folded or unfolded state\n");
611 printf("The following options change the default behavior in each state.\n");
612 printf("-m model[,group]\tFolded tension model (currently %s)\n",
613 tension_models[folded_tension_model].name);
614 printf("-a args\tFolded tension model argument string (currently %s)\n",
615 tension_models[folded_tension_model].params);
616 printf("-M model[,group]\tUnfolded tension model (currently %s)\n",
617 tension_models[unfolded_tension_model].name);
618 printf("-A args\tUnfolded tension model argument string (currently %s)\n",
619 tension_models[unfolded_tension_model].params);
620 printf("The following options change the unfolding rate.\n");
621 printf("-k model\tTransition rate model (currently %s)\n",
622 k_models[k_model].name);
623 printf("-K args\tTransition rate model argument string (currently %s)\n",
624 k_models[k_model].params);
625 printf("Domain creation options:\n");
626 printf("Once you've set up the appropriate default models, you need to add the domains\n");
627 printf("-F n\tAdd n folded domains with the current models\n");
628 printf("-U n\tAdd n unfolded domains with the current models\n");
629 printf("Output mode options:\n");
630 printf("There are two output modes. In standard mode, only the unfolding\n");
631 printf("events are printed. For example:\n");
632 printf(" #Unfolding Force (N)\n");
633 printf(" 123.456e-12\n");
635 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
636 printf(" #Position (m)\tForce (N)\n");
637 printf(" 0.001\t0.0023\n");
639 printf("-V\tChange output to verbose mode\n");
640 printf("-h\tPrint this help and exit\n");
642 printf("Tension models:\n");
643 for (i=0; i<n_tension_models; i++) {
644 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
645 for (j=0; j < tension_models[i].num_params ; j++)
646 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
647 printf("\t\tdefaults: %s\n", tension_models[i].params);
649 printf("Unfolding rate models:\n");
650 for (i=0; i<n_k_models; i++) {
651 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
652 for (j=0; j < k_models[i].num_params ; j++)
653 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
654 printf("\t\tdefaults: %s\n", k_models[i].params);
656 printf("Examples:\n");
657 printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded\n");
658 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);
659 printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded\n");
660 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);
664 \subsection{Get options}
667 void get_options(int argc, char **argv,
668 double *pP, double *pDt_max, double *pV,
670 int n_tension_models, tension_model_getopt_t *tension_models,
671 int n_k_models, k_model_getopt_t *k_models,
672 list_t **pDomain_list, int *pNum_folded)
674 char *prog_name = NULL;
675 char c, options[] = "P:t:v:T:C:m:a:M:A:k:K:F:U:Vh";
676 int ftension_model=0, utension_model=0, k_model=0;
677 char *ftension_params=NULL, *utension_params=NULL;
678 int i, n, ftension_group=0, utension_group=0;
682 extern int optind, optopt, opterr;
687 for (i=0; i<n_tension_models; i++) {
688 fprintf(stderr, "tension model %i %s:\t model's handler %p\n",
689 i, tension_models[i].name, tension_models[i].handler);
690 assert(tension_models[i].handler == tension_handlers[i]);
695 flags = FLAG_OUTPUT_UNFOLDING_FORCES;
697 *pP = 1e-3; /* % pop per s */
698 *pDt_max = 0.001; /* s */
699 *pV = 1e-6; /* m/s */
700 env->T = 300.0; /* K */
702 while ((c=getopt(argc, argv, options)) != -1) {
704 case 'P': *pP = atof(optarg); break;
705 case 't': *pDt_max = atof(optarg); break;
706 case 'v': *pV = atof(optarg); break;
707 case 'T': env->T = atof(optarg); break;
708 case 'C': env->T = atof(optarg)+273.15; break;
710 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
711 assert(num_strings == 1 || num_strings == 2);
712 if (num_strings == 1)
715 ftension_group = atoi(strings[1]);
716 ftension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
717 ftension_params = tension_models[ftension_model].params;
718 free_parsed_list(num_strings, strings);
721 ftension_params = optarg;
724 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
725 assert(num_strings == 1 || num_strings == 2);
726 if (num_strings == 1)
729 utension_group = atoi(strings[1]);
730 utension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
731 utension_params = tension_models[utension_model].params;
732 free_parsed_list(num_strings, strings);
735 utension_params = optarg;
738 k_model = index_k_model(n_k_models, k_models, optarg);
741 k_models[k_model].params = optarg;
746 for (i=0; i<n; i++) {
747 push(pDomain_list, generate_domain(DS_FOLDED,
748 tension_models+ftension_model,
751 tension_models+utension_model,
754 k_models+k_model), 1);
761 for (i=0; i<n; i++) {
762 push(pDomain_list, generate_domain(DS_UNFOLDED,
763 tension_models+ftension_model,
766 tension_models+utension_model,
769 k_models+k_model), 1);
773 flags = FLAG_OUTPUT_FULL_CURVE;
776 fprintf(stderr, "unrecognized option '%c'\n", optopt);
777 /* fall through to default case */
779 help(prog_name, *pP, *pDt_max, *pV, env, n_tension_models, tension_models, ftension_model, utension_model, n_k_models, k_models, k_model);
783 /* check the arguments */
784 assert(*pDomain_list != NULL); /* no domains! */
785 assert(*pP > 0.0); assert(*pP < 1.0);
786 assert(*pDt_max > 0.0);
788 assert(env->T > 0.0);
793 <<index tension model>>=
794 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
797 for (i=0; i<n_models; i++)
798 if (STRMATCH(models[i].name, name))
801 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
809 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
812 for (i=0; i<n_models; i++)
813 if (STRMATCH(models[i].name, name))
816 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
823 <<initialize model definition>>=
824 /* requires int num_param_args and char **param_args in the current scope
826 * INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
827 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
829 #define INIT_MODEL(role, model, param_string, param_pointer) \
831 parse_list_string(param_string, SEP, '{', '}', \
832 &num_param_args, ¶m_args); \
833 if (num_param_args != model->num_params) { \
835 "%s model %s expected %d params,\n", \
836 role, model->name, model->num_params); \
838 "not the %d params in '%s'\n", \
839 num_param_args, param_string); \
840 assert(num_param_args == model->num_params); \
842 if (model->creator) \
843 param_pointer = (*model->creator)(param_args); \
845 param_pointer = NULL; \
846 free_parsed_list(num_param_args, param_args); \
851 void *generate_domain(enum domain_state_t state,
852 tension_model_getopt_t *folded_model,
854 char *folded_param_string,
855 tension_model_getopt_t *unfolded_model,
857 char *unfolded_param_string,
858 k_model_getopt_t *k_model)
860 void *folded_params, *unfolded_params, *k_params;
861 int num_param_args; /* for INIT_MODEL() */
862 char **param_args; /* for INIT_MODEL() */
865 fprintf(stderr, "generating %s ", state==DS_FOLDED ? "folded" : "unfolded");
866 fprintf(stderr, "domain (k: \"%s\", \"%s\", f: \"%s\" (%p) group %d \"%s\", u: \"%s\" (%p) group %d \"%s\")\n",
867 k_model->name, k_model->params,
868 folded_model->name, folded_model->handler, folded_group, folded_param_string,
869 unfolded_model->name, unfolded_model->handler, unfolded_group, unfolded_param_string);
871 INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
872 INIT_MODEL("unfolded", unfolded_model, unfolded_param_string, unfolded_params);
873 INIT_MODEL("k", k_model, k_model->params, k_params);
874 return create_domain(state,
875 k_model->k, k_params, k_model->destructor,
876 folded_model->handler, folded_group, folded_params, folded_model->destructor,
877 unfolded_model->handler, unfolded_group, unfolded_params, unfolded_model->destructor
884 \addcontentsline{toc}{section}{Appendicies}
886 \section{sawsim.c details}
890 The general layout of our simulation code is:
900 We include [[math.h]], so don't forget to link to the libm with `-lm'.
902 #include <assert.h> /* assert() */
903 #include <stdlib.h> /* malloc(), free(), rand() */
904 #include <stdio.h> /* fprintf(), stdout */
905 #include <string.h> /* strlen, strtok() */
906 #include <math.h> /* exp(), M_PI, sqrt() */
907 #include <sys/types.h> /* pid_t (returned by getpid()) */
908 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
911 #include "tension_balance.h"
913 #include "tension_model.h"
919 <<version definition>>
921 <<max/min definitions>>
922 <<string comparison definition and globals>>
923 <<initialize model definition>>
924 <<get options definitions>>
925 <<domain definitions>>
934 <<create/destroy domain>>
935 <<destroy domain list>>
937 <<simulation functions>>
938 <<boilerplate functions>>
941 <<boilerplate functions>>=
943 <<get option functions>>
949 srand(getpid()*time(NULL)); /* seed rand() */
953 <<flag definitions>>=
954 /* in octal b/c of prefixed '0' */
955 #define FLAG_OUTPUT_FULL_CURVE 01
956 #define FLAG_OUTPUT_UNFOLDING_FORCES 02
960 static unsigned long int flags = 0;
963 \subsection{Utilities}
966 <<max/min definitions>>=
967 #define MAX(a,b) ((a)>(b) ? (a) : (b))
968 #define MIN(a,b) ((a)<(b) ? (a) : (b))
971 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
972 <<string comparison definition and globals>>=
973 // Check if two strings match, return 1 if they do
974 static char *temp_string_A;
975 static char *temp_string_B;
976 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
977 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
978 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
979 /* +1 to also compare the '\0' */
982 We also define a macro for our [[check]] unit testing
983 <<check relative error>>=
984 #define CHECK_ERR(max_err, expected, received) \
986 fail_unless( (received-expected)/expected < max_err, \
987 "relative error %g >= %g in %s (Expected %g, received %g)", \
988 (received-expected)/expected, max_err, #received, \
989 expected, received); \
990 fail_unless(-(received-expected)/expected < max_err, \
991 "relative error %g >= %g in %s (Expected %g, received %g)", \
992 -(received-expected)/expected, max_err, #received, \
993 expected, received); \
998 int happens(double probability)
1000 assert(probability >= 0.0); assert(probability <= 1.0);
1001 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*/
1006 \subsection{Adaptive timesteps}
1007 \label{app.adaptive_dt}
1009 $F(x)$ increases with $x$, possibly exploding, as in the worm-like chain model,
1010 so basing the timestep on the the unfolding probability at the current tension
1011 is dangerous, and we need to search for a $dt$ for which
1012 $P(F(x+v*dt)) < P_\text{target}$.
1013 There are two cases to consider.
1014 In the most common, no domains have unfolded since the last step,
1015 and we expect the next step to be slightly shorter than the previous one.
1016 In the less common, domains did unfold in the last step,
1017 and we expect the next step to be considerably longer than the previous one.
1019 double search_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers,
1020 list_t *domain_list,
1021 environment_t *env, double x,
1022 double target_prob, double max_dt, double v)
1023 { /* Returns the timestep dt in seconds for the current folded domain.
1024 * Takes a list of tension handlers, the list of domains,
1025 * a pointer env to the environmental data, a starting separation x in m,
1026 * a target_prob between 0 and 1,
1027 * max_dt in s, stretching velocity v in m/s.
1029 double F, k, dtCur, dtU, dtUCur, dtL, dt;
1031 /* get upper bound using the current position */
1032 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x, 0); /* BUG. repeated calculation */
1033 //printf("Start with x = %g (F = %g)\n", x, F);
1034 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1035 //printf("x %g\tF %g\tk %g\n", x, F, k);
1036 dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */
1038 //printf("overshot max_dt\n");
1041 /* set a lower bound on dt too */
1044 /* The dt determined above may produce illegitimate forces or ks.
1045 * Reduce the upper bound until we have valid ks. */
1047 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1048 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1051 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1053 //printf("Try for dt = %g (F = %g)\n", dt, F);
1054 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1055 /* returned k may be -1.0 */
1056 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1057 while (k == -1.0) { /* reduce step until we hit a valid k */
1059 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1060 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1061 //printf("Try for dt = %g (F = %g)\n", dt, F);
1062 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1063 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1065 assert(dtU > 1e-14); /* timestep to valid k too small */
1066 dtUCur = target_prob / k; /* safe timestep back from x+dtU */
1068 return dt; /* dtU is safe. */
1070 /* dtU wasn't safe, lets see what would be. */
1071 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1072 dt = (dtU + dtL) / 2.0;
1073 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1074 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1075 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1076 dtCur = target_prob / k;
1077 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1078 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1080 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1081 dtU = dt; /* ... stepping out only dtCur would be safe */
1084 break; /* dtCur = dt */
1086 return MAX(dtUCur, dtL);
1090 To determine $dt$ for an array of potentially different folded domains,
1091 we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
1094 double determine_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list,
1095 environment_t *env, double x,
1096 double target_prob, double dt_max, double v)
1097 { /* Returns the timestep dt in seconds.
1098 * Takes the list of folded domains, target_prob between 0 and 1,
1099 * F in N, and T in K. */
1100 double dt=dt_max, new_dt;
1101 assert(target_prob > 0.0); assert(target_prob < 1.0);
1102 assert(dt_max > 0.0);
1104 /* .5 nm steps = v * dt */
1106 while (domain_list != NULL) {
1107 if (D(domain_list)->state == DS_FOLDED) {
1108 new_dt = search_dt(num_tension_handlers, tension_handlers, domain_list, env, x, target_prob, dt, v);
1109 dt = MIN(dt, new_dt);
1111 domain_list = domain_list->next;
1117 \subsection{Domain data}
1119 Currently domains exist in two states, folded and unfolded, and the
1120 only allowed transitions are folded $\rightarrow$ unfolded. Of
1121 course, it wouldn't be too complicated to extent this to a multi-state
1122 system, with an array containing the domains group for each possible
1123 state, and a matrix of transition-rate-calculating functions.
1124 However, at this point such generality seems unnecessary at this
1126 <<domain definitions>>=
1127 enum domain_state_t {DS_FOLDED,
1131 typedef struct domain_struct {
1132 enum domain_state_t state;
1133 one_dim_fn_t *folded_handler;
1135 one_dim_fn_t *unfolded_handler;
1137 k_func_t *k; /* function returning unfolding rate */
1138 void *folded_params; /* pointer to folded parameters */
1139 void *unfolded_params; /* pointer to unfolded parameters */
1140 void *k_params; /* pointer to k parameters */
1141 destroy_data_func_t *destroy_folded;
1142 destroy_data_func_t *destroy_unfolded;
1143 destroy_data_func_t *destroy_k;
1146 /* get the domain data for the current list node */
1147 #define D(list) ((domain_t *)(list)->d)
1148 /* get the tension params for the current list node */
1149 #define D_TP(list) (((domain_t *)(list)->d)->state == DS_FOLDED \
1150 ? ((domain_t *)(list)->d)->folded_params \
1151 : ((domain_t *)(list)->d)->unfolded_params)
1153 [[k]] is a pointer to the function determining the unfolding rate for a given tension.
1154 [[folded_params]] is a pointer to the parameters used by the function pointed to by [[k]].
1155 [[unfolded_params]] is a pointer to the parameters used by the group-appropriate handler function when determining the tension.
1156 The [[destroy_*]] pointers point to functions for freeing the memory [[*_params]].
1157 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.
1159 [[create_]] and [[destroy_domain]] are simple wrappers around [[malloc]] and [[free]].
1160 <<create/destroy domain>>=
1161 domain_t *create_domain(enum domain_state_t state,
1164 destroy_data_func_t *destroy_k,
1165 one_dim_fn_t *folded_handler,
1167 void *folded_params,
1168 destroy_data_func_t *destroy_folded,
1169 one_dim_fn_t *unfolded_handler,
1171 void *unfolded_params,
1172 destroy_data_func_t *destroy_unfolded)
1174 domain_t *ret = (domain_t *)malloc(sizeof(domain_t));
1175 assert(ret != NULL);
1176 if (state == DS_FOLDED) {
1177 assert(k != NULL); /* the pointer points somewhere valid */
1178 assert(*k != NULL); /* and there is something useful there */
1180 assert(state == DS_UNFOLDED);
1182 ret->folded_handler = folded_handler;
1183 ret->folded_group = folded_group;
1184 ret->unfolded_handler = unfolded_handler;
1185 ret->unfolded_group = unfolded_group;
1187 ret->folded_params = folded_params;
1188 ret->unfolded_params = unfolded_params;
1189 ret->k_params = k_params;
1190 ret->destroy_folded = destroy_folded;
1191 ret->destroy_unfolded = destroy_unfolded;
1192 ret->destroy_k = destroy_k;
1194 fprintf(stderr, "generated domain %p. state %d\n", ret, ret->state);
1199 void destroy_domain(domain_t *domain)
1202 //printf("domain %p & %p\n", *domain, domain);
1203 if (domain->destroy_folded)
1204 (*domain->destroy_folded)(domain->folded_params);
1205 if (domain->destroy_unfolded)
1206 (*domain->destroy_unfolded)(domain->unfolded_params);
1207 if (domain->destroy_k)
1208 (*domain->destroy_k)(domain->k_params);
1214 <<destroy domain list>>=
1215 void destroy_domain_list(list_t *domain_list)
1217 domain_list = head(domain_list);
1218 while (domain_list != NULL)
1219 destroy_domain((domain_t *) pop(&domain_list));
1223 \subsection{Domain model and group handling}
1225 <<group functions>>=
1226 <<get tension handler>>
1228 <<int comparison function>>
1229 <<find possible groups>>
1232 <<get active groups>>
1235 <<get tension handler>>=
1236 one_dim_fn_t *get_tension_handler(domain_t *domain)
1238 if (domain->state == DS_FOLDED)
1239 return domain->folded_handler;
1241 if (domain->state != DS_UNFOLDED)
1242 fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
1243 assert(domain->state == DS_UNFOLDED);
1244 return domain->unfolded_handler;
1250 int get_group(domain_t *domain)
1252 if (domain->state == DS_FOLDED)
1253 return domain->folded_group;
1255 if (domain->state != DS_UNFOLDED)
1256 fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
1257 assert(domain->state == DS_UNFOLDED);
1258 return domain->unfolded_group;
1263 We already know all possible tension classes, since they are all
1264 hardcoded into \prog. However, there may be any number of possible
1265 groups. We define a function [[find_possible_groups]] to search for
1266 possible (not neccessarily active) groups. It can search for
1267 subgroups of a single [[handler]], or by setting [[handler]] to
1268 [[NULL]], subgroups of \emph{any} handler. It returns a list of group
1270 <<find possible groups>>=
1271 list_t *find_possible_groups(list_t *list, one_dim_fn_t *handler) {
1274 while (list != NULL) {
1275 if (handler == NULL || get_tension_handler(D(list)) == handler) {
1276 push(&ret, &D(list)->folded_group, 1);
1277 push(&ret, &D(list)->unfolded_group, 1);
1283 return ret; /* no groups with this handler, no processing needed */
1285 /* sort the ret list, and remove duplicates */
1286 sort(&ret, &int_comparison);
1287 unique(&ret, &int_comparison);
1289 fprintf(stderr, "handler %p possible groups:", handler);
1291 while (list != NULL) {
1292 fprintf(stderr, " %d", *((int *)list->d));
1295 fprintf(stderr, "\n");
1301 Search a [[list]] of domains, and determine whether a particular model
1302 class and group number combination is active.
1303 <<is group active>>=
1304 int is_group_active(list_t *list, one_dim_fn_t *handler, int group)
1307 while (list != NULL) {
1308 if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group)
1316 Search a [[list]] of domains, and return all domains matching a
1317 particular model class and group number.
1319 list_t *get_group_list(list_t *list, one_dim_fn_t *handler, int group)
1323 while (list != NULL) {
1324 if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group) {
1325 push(&ret, D_TP(list), 1); /* add a pointer to the appropriate tension parameters to our new list. */
1327 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);
1335 Because all the node data in lists returned by [[get_group_list]] is
1336 also in the main domain list, you shouldn't destroy the node data
1337 popped off when destroying the group lists. It will all get cleaned
1338 up when the main domain list is destroyed.
1340 Put all this together to scan through a list of domains and construct
1341 an array of all the active groups.
1342 <<get active groups>>=
1343 void get_active_groups(list_t *list,
1344 int num_tension_handlers, one_dim_fn_t **tension_handlers,
1345 int *pNum_active_groups, one_dim_fn_t ***pPer_group_handlers, void ***pActive_groups)
1347 void **active_groups=NULL;
1348 one_dim_fn_t *handler, **per_group_handlers=NULL;
1349 int i, num_active_groups, *group;
1350 list_t *possible_group_numbers, *active_group_list, *active_groups_list=NULL, *per_group_handlers_list=NULL;
1352 /* get all the active groups in a list */
1353 for (i=0; i<num_tension_handlers; i++) {
1354 handler = tension_handlers[i];
1355 if (handler == NULL) continue; /* tensionless handler */
1356 possible_group_numbers = head(find_possible_groups(list, handler));
1357 while (possible_group_numbers != NULL) {
1358 group = (int *)pop(&possible_group_numbers);
1359 if (is_group_active(list, handler, *group) == 1) {
1360 active_group_list = get_group_list(list, handler, *group);
1361 push(&active_groups_list, active_group_list, 1);
1362 push(&per_group_handlers_list, handler, 1);
1364 fprintf(stderr, "active group %p for %p %d headed by tension parameters %p\n", active_group_list, handler, *group, head(active_group_list)->d);
1365 list_print(stderr, active_group_list, "active group");
1371 /* allocate the array we'll be returning */
1372 num_active_groups = list_length(active_groups_list);
1373 active_groups = (void **) malloc(sizeof(void *)*num_active_groups);
1374 assert(active_groups != NULL);
1375 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1376 assert(per_group_handlers != NULL);
1378 /* populate the array we'll be returning */
1379 active_groups_list = head(active_groups_list);
1380 for (i=0; i<num_active_groups; i++) {
1381 handler = pop(&per_group_handlers_list);
1382 per_group_handlers[i] = handler;
1384 active_groups[i] = (void *) malloc(sizeof(tension_handler_data_t));
1385 assert(active_groups[i] != NULL);
1386 active_group_list = pop(&active_groups_list);
1387 ((tension_handler_data_t *)active_groups[i])->group = active_group_list;
1388 ((tension_handler_data_t *)active_groups[i])->env = NULL;
1389 ((tension_handler_data_t *)active_groups[i])->persist = NULL;
1391 assert(active_groups_list == NULL);
1392 assert(per_group_handlers_list == NULL);
1394 *pNum_active_groups = num_active_groups;
1395 *pPer_group_handlers = per_group_handlers;
1396 *pActive_groups = active_groups;
1401 \section{String parsing}
1403 For handling command line arguments, we need parse delimited strings (e.g. [["param1,param2,param3..."]]).
1404 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1405 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1406 need to take care of parsing those parameters themselves.
1407 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1413 <<parse definitions>>
1414 <<parse declarations>>
1418 <<parse module makefile lines>>=
1419 NW_SAWSIM_MODS += parse
1420 CHECK_BINS += check_parse
1424 <<parse definitions>>=
1425 #define SEP ',' /* argument separator character */
1428 <<parse declarations>>=
1429 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1430 int *num, char ***string_array);
1431 extern void free_parsed_list(int num, char **string_array);
1434 [[parse_list_string]] allocates memory, don't forget to free it
1435 afterward with [[free_parsed_list]]. It does not alter the original.
1437 The string may start off with a [[deeper]] character
1438 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1439 [[x,y]], where the pointer is one character in on the copied string.
1440 However, when we go to free the memory, we need a pointer to the
1441 beginning of the string. In order to accommodate this for a string
1442 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1443 the first $N$ elements point to the separated arguments, and let the
1444 last element point to the start of the copied string regardless of
1446 <<parse delimited list string functions>>=
1447 /* TODO, split out into parse.hc */
1448 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1451 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1452 if (string[i] == deeper) {depth++;}
1453 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1459 void parse_list_string(char *string, char sep, char deeper, char shallower,
1460 int *num, char ***string_array)
1462 char *str=NULL, **ret=NULL;
1464 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
1466 *string_array = NULL;
1469 /* make a copy of the string, so we don't change the original */
1470 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1471 assert(str != NULL);
1472 strcpy(str, string); /* we know str is long enough */
1473 /* count the number of regions, so we can allocate pointers to them */
1476 n++; i++; /* move on to next argument */
1477 i += next_delim_index(str+i, sep, deeper, shallower);
1478 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1479 } while (str[i] != '\0');
1480 ret = (char **)malloc(sizeof(char *)*(n+1));
1481 assert(ret != NULL);
1482 /* replace the separators with '\0' & assign pointers */
1483 ret[n] = str; /* point to the front of the copied string */
1486 for(i=1; i<n; i++) {
1487 j += next_delim_index(str+j, sep, deeper, shallower);
1489 ret[i] = str+j; /* point to the separated arguments */
1491 /* strip off bounding braces, if any */
1492 for(i=0; i<n; i++) {
1493 if (ret[i][0] == deeper) {
1497 if (ret[i][strlen(ret[i])-1] == shallower)
1498 ret[i][strlen(ret[i])-1] = '\0';
1501 *string_array = ret;
1504 void free_parsed_list(int num, char **string_array)
1507 /* frees all items, since they were allocated as a single string*/
1508 free(string_array[num]);
1509 /* and free the array of pointers */
1515 \subsection{Parsing implementation}
1521 <<parse delimited list string functions>>
1525 #include <assert.h> /* assert() */
1526 #include <stdlib.h> /* NULL */
1527 #include <stdio.h> /* fprintf(), stdout *//*!!*/
1528 #include <string.h> /* strlen() */
1532 \subsection{Parsing unit tests}
1534 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1536 <<parse check includes>>
1537 <<string comparison definition and globals>>
1538 <<check relative error>>
1539 <<parse test suite>>
1540 <<main check program>>
1543 <<parse check includes>>=
1544 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
1545 #include <stdio.h> /* printf() */
1546 #include <assert.h> /* assert() */
1547 #include <string.h> /* strlen() */
1552 <<parse test suite>>=
1553 <<parse list string tests>>
1554 <<parse suite definition>>
1557 <<parse suite definition>>=
1558 Suite *test_suite (void)
1560 Suite *s = suite_create ("k model");
1561 <<parse list string test case defs>>
1563 <<parse list string test case add>>
1568 <<parse list string tests>>=
1571 START_TEST(test_next_delim_index)
1573 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
1574 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
1575 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
1576 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
1577 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
1581 START_TEST(test_parse_list_null)
1585 parse_list_string(NULL, SEP, '{', '}',
1586 &num_param_args, ¶m_args);
1587 fail_unless(num_param_args == 0, NULL);
1588 fail_unless(param_args == NULL, NULL);
1591 START_TEST(test_parse_list_single_simple)
1595 parse_list_string("arg", SEP, '{', '}',
1596 &num_param_args, ¶m_args);
1597 fail_unless(num_param_args == 1, NULL);
1598 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
1601 START_TEST(test_parse_list_single_compound)
1605 parse_list_string("{x,y,z}", SEP, '{', '}',
1606 &num_param_args, ¶m_args);
1607 fail_unless(num_param_args == 1, NULL);
1608 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
1611 START_TEST(test_parse_list_double_simple)
1615 parse_list_string("abc,def", SEP, '{', '}',
1616 &num_param_args, ¶m_args);
1617 fail_unless(num_param_args == 2, NULL);
1618 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
1619 fail_unless(STRMATCH(param_args[1],"def"), NULL);
1623 <<parse list string test case defs>>=
1624 TCase *tc_parse_list_string = tcase_create("parse list string");
1626 <<parse list string test case add>>=
1627 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
1628 tcase_add_test(tc_parse_list_string, test_parse_list_null);
1629 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
1630 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
1631 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
1632 suite_add_tcase(s, tc_parse_list_string);
1636 \section{Unit tests}
1638 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1645 <<check relative error>>
1648 <<main check program>>
1660 <<determine dt tests>>
1662 <<does domain unfold tests>>
1663 <<randomly unfold domains tests>>
1664 <<suite definition>>
1667 <<suite definition>>=
1668 Suite *test_suite (void)
1670 Suite *s = suite_create ("sawsim");
1671 <<F test case defs>>
1672 <<determine dt test case defs>>
1673 <<happens test case defs>>
1674 <<does domain unfold test case defs>>
1675 <<randomly unfold domains test case defs>>
1678 <<determine dt test case add>>
1679 <<happens test case add>>
1680 <<does domain unfold test case add>>
1681 <<randomly unfold domains test case add>>
1684 tcase_add_checked_fixture(tc_strip_address,
1685 setup_strip_address,
1686 teardown_strip_address);
1692 <<main check program>>=
1696 Suite *s = test_suite();
1697 SRunner *sr = srunner_create(s);
1698 srunner_run_all(sr, CK_ENV);
1699 number_failed = srunner_ntests_failed(sr);
1701 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
1705 \subsection{F tests}
1707 <<worm-like chain tests>>
1709 <<F test case defs>>=
1710 <<worm-like chain test case def>>
1712 <<F test case add>>=
1713 <<worm-like chain test case add>>
1716 <<worm-like chain tests>>=
1717 extern double wlc(double x, double T, double p, double L);
1718 START_TEST(test_wlc_at_zero)
1720 double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
1721 fail_unless(abs(wlc(x, T, p, L)) < lim, \
1722 "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
1723 x, T, p, L, abs(wlc(x, T, p, L)), lim);
1727 START_TEST(test_wlc_at_half)
1729 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
1730 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
1731 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
1733 * linear term = x/L = 0.5
1734 * nonlinear + linear = 0.75 + 0.5 = 1.25
1735 * wlc = 10e21*1.25 = 12.5e21
1737 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
1738 "wlc(%g, %g, %g, %g) = %g != %g",
1739 x, T, p, L, wlc(x, T, p, L), 12.5e21);
1743 <<worm-like chain test case def>>=
1744 TCase *tc_wlc = tcase_create("WLC");
1747 <<worm-like chain test case add>>=
1748 tcase_add_test(tc_wlc, test_wlc_at_zero);
1749 tcase_add_test(tc_wlc, test_wlc_at_half);
1750 suite_add_tcase(s, tc_wlc);
1753 \subsection{Model tests}
1755 Check the searching with [[linear_k]].
1756 Check overwhelming force treatment with the heavyside-esque step [[?]].
1757 <<determine dt tests>>=
1758 double linear_k(double F, environment_t *env, void *params)
1760 double Fnot = *(double *)params;
1764 START_TEST(test_determine_dt_linear_k)
1767 double dt_max=3.0, Fnot=3.0;
1768 double F[]={0,1,2,3,4,5,6};
1769 domain_t dom; /* use both parts at once for folded/unfolded */
1773 dom->next = dom->prev = NULL;
1774 dom->k_func_t = linear_k;
1775 dom->folded_params = &Fnot;
1776 dom->unfolded_params = !!!!!!!!!
1777 dom->destroy_folded = dom->destroy_unfolded = NULL;
1778 for( i=0; i < sizeof(F)/sizeof(double); i++) {
1779 fail_unless(determine_dt(folded, unfolded, targ_p, dt_max, v, x, T)==1, NULL);
1785 <<determine dt test case defs>>=
1786 TCase *tc_determine_dt = tcase_create("Determine dt");
1788 <<determine dt test case add>>=
1789 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
1790 suite_add_tcase(s, tc_determine_dt);
1795 <<happens test case defs>>=
1797 <<happens test case add>>=
1800 <<does domain unfold tests>>=
1802 <<does domain unfold test case defs>>=
1804 <<does domain unfold test case add>>=
1807 <<randomly unfold domains tests>>=
1809 <<randomly unfold domains test case defs>>=
1811 <<randomly unfold domains test case add>>=
1815 \section{Balancing group extensions}
1816 \label{app.tension_balance}
1818 For a given total extension $x$ (of the piezo), the various domain
1819 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
1820 amounts, and we need to tweak the portion that each extends to balance
1821 the tension amongst the active groups. Since the tension balancing is
1822 separable from the bulk of the simulation, we break it out into a
1823 separate module. The interface is defined in [[tension_balance.h]],
1824 the implementation in [[tension_balance.c]], and the unit testing in
1825 [[check_tension_balance.c]]
1827 <<tension-balance.h>>=
1828 #ifndef TENSION_BALANCE_H
1829 #define TENSION_BALANCE_H
1831 <<tension balance function declaration>>
1832 #endif /* TENSION_BALANCE_H */
1835 <<tension balance functions>>=
1836 <<one dimensional bracket>>
1837 <<one dimensional solve>>
1839 <<tension balance function>>
1842 <<tension balance module makefile lines>>=
1843 NW_SAWSIM_MODS += tension_balance
1844 CHECK_BINS += check_tension_balance
1845 check_tension_balance_MODS =
1848 The entire force balancing problem reduces to a solving two nested
1849 one-dimensional root-finding problems. First we define the one
1850 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
1851 non-empty group). $x(x_0)$ is also strictly monotonically increasing,
1852 so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
1853 stores the last successful value of $x$, since we expect to be taking
1854 small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
1855 indicates that the groups have changed.
1856 <<tension balance function declaration>>=
1857 double tension_balance(int num_tension_groups,
1858 one_dim_fn_t **tension_handlers,
1863 <<tension balance function>>=
1864 double tension_balance(int num_tension_groups,
1865 one_dim_fn_t **tension_handlers,
1870 static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL};
1871 double F, xo_guess, xo, lb, ub;
1872 double min_relx=1e-6, min_rely=1e-6;
1873 int max_steps=200, i;
1875 assert(num_tension_groups > 0);
1876 assert(tension_handlers != NULL);
1877 assert(params != NULL);
1878 assert(num_tension_groups > 0);
1880 if (num_tension_groups == 1) { /* only one group, no balancing required */
1883 //fprintf(stderr, "balancing for x = %g with ", x);
1884 //for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
1885 //fprintf(stderr, "\n");
1886 if (last_x == -1) { /* new group setup, reset x_xo_data */
1887 /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
1888 if (x_xo_data.xi != NULL)
1890 /* construct new x_xo_data */
1891 x_xo_data.n_groups = num_tension_groups;
1892 x_xo_data.tension_handler = tension_handlers;
1893 x_xo_data.handler_data = params;
1894 x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
1895 for (i=0; i<num_tension_groups; i++)
1896 x_xo_data.xi[i] = -1.0;
1898 if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
1899 if (x == 0) xo_guess = length_scale;
1900 else xo_guess = x/num_tension_groups;
1902 fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
1904 oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
1905 } else { /* work off the last balanced point */
1907 fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, last_x, x-last_x, x);
1910 lb = x_xo_data.xi[0];
1911 ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */
1912 } else if (x < last_x) {
1913 lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
1914 ub = x_xo_data.xi[0];
1915 } else { /* x == last_x */
1916 fprintf(stderr, "not moving\n");
1917 lb= x_xo_data.xi[0];
1918 ub= x_xo_data.xi[0];
1922 fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
1923 __FUNCTION__, x, lb, ub);
1925 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
1926 F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
1928 if (num_tension_groups > 1) {
1929 fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
1930 for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
1931 fprintf(stderr, "\n");
1935 F = (*tension_handlers[0])(xo, params[0]);
1940 <<tension balance internal definitions>>=
1941 <<x of x0 definitions>>
1944 <<x of x0 definitions>>=
1945 typedef struct x_of_xo_data_struct {
1947 one_dim_fn_t **tension_handler; /* array of fn pointers */
1948 void **handler_data; /* array of void* pointers */
1954 double x_of_xo(double xo, void *pdata)
1956 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
1957 double F, x=0, xi, xi_guess, lb, ub;
1959 double min_relx=1e-6, min_rely=1e-6;
1961 assert(data->n_groups > 0);
1963 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
1965 fprintf(stderr, "%s: finding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
1968 if (data->xi) data->xi[0] = xo;
1969 for (i=1; i < data->n_groups; i++) {
1970 if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
1971 else xi_guess = data->xi[i];
1973 fprintf(stderr, "%s: searching for proper x[%d] for F[%d](x[%d]) = %g N\n", __FUNCTION__, i, i, i, F);
1975 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
1976 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
1978 fprintf(stderr, "%s: found F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
1982 if (data->xi) data->xi[i] = xi;
1985 fprintf(stderr, "%s: found x(x0=%g) = %g\n", __FUNCTION__, xo, x);
1991 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
1992 number of steps for monotonically increasing $f(x)$. Simple
1993 bisection, so it's robust and converges fairly quickly. You can set
1994 the maximum number of steps to take, but if you have enough processor
1995 power, it's better to limit your precision with the relative limits
1996 [[min_relx]] and [[y]]. These ensure that the width of the bracket is
1997 small on the length scale of it's larger side. Note that \emph{both}
1998 relative limits must be satisfied for the search to stop.
1999 <<one dimensional function definition>>=
2000 /* equivalent to gsl_func_t */
2001 typedef double one_dim_fn_t(double x, void *params);
2003 <<one dimensional solve declaration>>=
2004 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2005 double min_relx, double min_rely, int max_steps,
2009 <<one dimensional solve>>=
2010 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2011 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2012 double min_relx, double min_rely, int max_steps,
2015 double yx, ylb, yub, x;
2018 ylb = (*f)(lb, params);
2019 yub = (*f)(ub, params);
2021 /* check some simple cases */
2023 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bounded */
2024 else return lb; /* any x on the range [lb, ub] would work */
2026 if (ylb == y) { x=lb; yx=ylb; goto end; }
2027 if (yub == y) { x=ub; yx=yub; goto end; }
2030 fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2032 assert(ylb < y); assert(yub > y);
2033 for (i=0; i<max_steps; i++) {
2035 yx = (*f)(x, params);
2037 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);
2039 if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2040 else if (yx > y) { ub=x; yub=yx; }
2041 else /* yx < y */ { lb=x; ylb=yx; }
2046 fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2052 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2053 Generate bracketing $x$ values through bisection or exponential growth.
2054 <<one dimensional bracket declaration>>=
2055 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2058 <<one dimensional bracket>>=
2059 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2061 double yx, step, x=xguess;
2063 yx = (*f)(x, params);
2065 fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2072 fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2076 if (x == 0) x = length_scale/2.0; /* HACK */
2079 yx = (*f)(x, params);
2081 fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2086 fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2089 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2093 \subsection{Balancing implementation}
2095 <<tension-balance.c>>=
2098 <<tension balance includes>>
2099 #include "tension_balance.h"
2101 double length_scale = 1e-10; /* HACK */
2103 <<tension balance internal definitions>>
2104 <<tension balance functions>>
2107 <<tension balance includes>>=
2108 #include <assert.h> /* assert() */
2109 #include <stdlib.h> /* NULL */
2110 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2111 #include <stdio.h> /* fprintf(), stdout */
2115 \subsection{Balancing unit tests}
2117 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2118 <<check-tension-balance.c>>=
2119 <<tension balance check includes>>
2120 <<tension balance test suite>>
2121 <<main check program>>
2124 <<tension balance check includes>>=
2125 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2126 #include <assert.h> /* assert() */
2129 #include "tension_balance.h"
2132 <<tension balance test suite>>=
2133 <<tension balance function tests>>
2134 <<tension balance suite definition>>
2137 <<tension balance suite definition>>=
2138 Suite *test_suite (void)
2140 Suite *s = suite_create ("tension balance");
2141 <<tension balance function test case defs>>
2143 <<tension balance function test case adds>>
2148 <<tension balance function tests>>=
2149 <<check relative error>>
2151 double hooke(double x, void *pK)
2154 return *((double*)pK) * x;
2157 START_TEST(test_single_function)
2159 double k=5, x=3, last_x=2.0, F;
2160 one_dim_fn_t *handlers[] = {&hooke};
2161 void *data[] = {&k};
2162 F = tension_balance(1, handlers, data, last_x, x);
2163 fail_unless(F = k*x, NULL);
2168 We can also test balancing two springs with different spring constants.
2169 The analytic solution for a general number of springs is given in Appendix \ref{app.math_hooke}.
2170 <<tension balance function tests>>=
2171 START_TEST(test_double_hooke)
2173 double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e;
2174 one_dim_fn_t *handlers[] = {&hooke, &hooke};
2175 void *data[] = {&k1, &k2};
2176 F = tension_balance(2, handlers, data, last_x, x);
2180 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2181 //CHECK_ERR(1e-6, x1e, xi[0]);
2182 //CHECK_ERR(1e-6, x2e, xi[1]);
2183 CHECK_ERR(1e-6, Fe, F);
2187 TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above.
2189 <<tension balance function test case defs>>=
2190 TCase *tc_tbfunc = tcase_create("tension balance function");
2193 <<tension balance function test case adds>>=
2194 tcase_add_test(tc_tbfunc, test_single_function);
2195 tcase_add_test(tc_tbfunc, test_double_hooke);
2196 suite_add_tcase(s, tc_tbfunc);
2201 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2202 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2203 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2209 <<list definitions>>
2210 <<list declarations>>
2214 <<list declarations>>=
2215 <<head and tail declarations>>
2216 <<list length declaration>>
2217 <<push declaration>>
2219 <<list destroy declaration>>
2220 <<list to array declaration>>
2221 <<move declaration>>
2222 <<sort declaration>>
2223 <<unique declaration>>
2224 <<list print declaration>>
2228 <<create and destroy node>>
2241 <<list module makefile lines>>=
2242 NW_SAWSIM_MODS += list
2243 CHECK_BINS += check_list
2247 <<list definitions>>=
2248 typedef struct list_struct {
2249 struct list_struct *next;
2250 struct list_struct *prev;
2254 <<comparison function typedef>>
2257 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2258 <<head and tail declarations>>=
2259 list_t *head(list_t *list);
2260 list_t *tail(list_t *list);
2263 list_t *head(list_t *list)
2267 while (list->prev != NULL) {
2273 list_t *tail(list_t *list)
2277 while (list->next != NULL) {
2284 <<list length declaration>>=
2285 int list_length(list_t *list);
2288 int list_length(list_t *list)
2295 while (list->next != NULL) {
2303 [[push]] inserts a node relative to the current position in the list
2304 without changing the current position.
2305 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2306 so we set it to that of the pushed domain.
2307 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2308 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2309 <<push declaration>>=
2310 void push(list_t **pList, void *data, int below);
2313 void push(list_t **pList, void *data, int below)
2315 list_t *list, *new_node;
2316 assert(pList != NULL);
2318 new_node = create_node(data);
2321 else if (below==0) { /* insert above */
2322 if (list->prev != NULL)
2323 list->prev->next = new_node;
2324 new_node->prev = list->prev;
2325 list->prev = new_node;
2326 new_node->next = list;
2327 } else { /* insert below */
2328 if (list->next != NULL)
2329 list->next->prev = new_node;
2330 new_node->next = list->next;
2331 list->next = new_node;
2332 new_node->prev = list;
2337 [[pop]] removes the current domain node, moving the current position
2338 to the node after it, unless that node is [[NULL]], in which case move
2339 the current position to the node before it.
2340 <<pop declaration>>=
2341 void *pop(list_t **pList);
2344 void *pop(list_t **pList)
2346 list_t *list, *popped;
2348 assert(pList != NULL);
2350 assert(list != NULL); /* not an empty list */
2352 /* bypass the popped node */
2353 if (list->prev != NULL)
2354 list->prev->next = list->next;
2355 if (list->next != NULL) {
2356 list->next->prev = list->prev;
2357 *pList = list->next;
2359 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
2361 destroy_node(popped);
2366 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
2367 If the cleanup function is [[NULL]], no cleanup function is called.
2368 The nodes are not popped in any particular order.
2369 <<list destroy declaration>>=
2370 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
2373 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
2377 assert(pList != NULL);
2380 assert(list != NULL); /* not an empty list */
2381 while (list != NULL) {
2383 if (destroy != NULL)
2389 [[list_to_array]] converts a list to an array of pointers, preserving order.
2390 <<list to array declaration>>=
2391 void list_to_array(list_t *list, int *length, void ***array);
2394 void list_to_array(list_t *list, int *pLength, void ***pArray)
2398 assert(list != NULL);
2399 assert(pLength != NULL);
2400 assert(pArray != NULL);
2402 length = list_length(list);
2403 array = (void **)malloc(sizeof(void *)*length);
2404 assert(array != NULL);
2407 while (list != NULL) {
2408 array[i++] = list->d;
2416 [[move]] moves the current element along the list in either direction.
2417 <<move declaration>>=
2418 void move(list_t **pList, int downstream);
2421 void move(list_t **pList, int downstream)
2423 list_t *list, *flipper;
2425 assert(pList != NULL);
2426 list = *pList; /* a>B>c>d */
2427 assert(list != NULL); /* not an empty list */
2428 if (downstream == 0)
2429 flipper = list->prev; /* flipper is A */
2431 flipper = list->next; /* flipper is C */
2432 /* ensure there is somewhere to go in the direction we're heading */
2433 assert(flipper != NULL);
2434 /* remove the current element */
2435 data = pop(&list); /* data=B, a>C>d */
2436 /* and put it back in in the correct spot */
2438 if (downstream == 0) {
2439 push(&list, data, 0); /* b>A>c>d */
2440 list = list->prev; /* B>a>c>d */
2442 push(&list, data, 1); /* a>C>b>d */
2443 list = list->next; /* a>c>B>d */
2449 [[sort]] sorts a list from smallest to largest according to a user
2450 specified comparison.
2451 <<comparison function typedef>>=
2452 typedef int (comparison_fn_t)(void *A, void *B);
2455 <<int comparison function>>=
2456 int int_comparison(void *A, void *B)
2461 if (a > b) return 1;
2462 else if (a == b) return 0;
2466 <<sort declaration>>=
2467 void sort(list_t **pList, comparison_fn_t *comp);
2469 Here we use the bubble sort algorithm.
2471 void sort(list_t **pList, comparison_fn_t *comp)
2475 assert(pList != NULL);
2477 assert(list != NULL); /* not an empty list */
2478 while (stable == 0) {
2481 while (list->next != NULL) {
2482 if ((*comp)(list->d, list->next->d) > 0) {
2484 move(&list, 1 /* downstream */);
2493 [[unique]] removes duplicates from a sorted list according to a user
2494 specified comparison. Don't do this unless you have other pointers
2495 any dynamically allocated data in your list, because [[unique]] just
2496 drops any non-unique data without freeing it.
2497 <<unique declaration>>=
2498 void unique(list_t **pList, comparison_fn_t *comp);
2501 void unique(list_t **pList, comparison_fn_t *comp)
2504 assert(pList != NULL);
2506 assert(list != NULL); /* not an empty list */
2508 while (list->next != NULL) {
2509 if ((*comp)(list->d, list->next->d) == 0) {
2518 [[list_print]] prints a list to a [[FILE *]] stream.
2519 <<list print declaration>>=
2520 void list_print(FILE *file, list_t *list, char *list_name);
2523 void list_print(FILE *file, list_t *list, char *list_name)
2525 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
2527 while (list != NULL) {
2528 fprintf(file, " %p", list->d);
2531 fprintf(file, "\n");
2535 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
2536 <<create and destroy node>>=
2537 list_t *create_node(void *data)
2539 list_t *ret = (list_t *)malloc(sizeof(list_t));
2540 assert(ret != NULL);
2547 void destroy_node(list_t *node)
2553 The user must free the data pointed to by the node on their own.
2555 \subsection{List implementation}
2565 #include <assert.h> /* assert() */
2566 #include <stdlib.h> /* malloc(), free(), rand() */
2567 #include <stdio.h> /* fprintf(), stdout, FILE */
2568 #include "global.h" /* destroy_data_func_t */
2571 \subsection{List unit tests}
2573 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2575 <<list check includes>>
2577 <<main check program>>
2580 <<list check includes>>=
2581 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2582 #include <stdio.h> /* FILE */
2588 <<list test suite>>=
2591 <<list suite definition>>
2594 <<list suite definition>>=
2595 Suite *test_suite (void)
2597 Suite *s = suite_create ("list");
2598 <<push test case defs>>
2599 <<pop test case defs>>
2601 <<push test case adds>>
2602 <<pop test case adds>>
2608 START_TEST(test_push)
2613 push(&list, (void *)i, 1);
2614 fail_unless(list == head(list), NULL);
2615 fail_unless( (int)list->d == 0 );
2616 for (i=0; i<n; i++) {
2617 /* we expect 0, n-1, n-2, ..., 1 */
2620 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
2626 <<push test case defs>>=
2627 TCase *tc_push = tcase_create("push");
2630 <<push test case adds>>=
2631 tcase_add_test(tc_push, test_push);
2632 suite_add_tcase(s, tc_push);
2637 <<pop test case defs>>=
2639 <<pop test case adds>>=
2642 \section{Function string evaluation}
2644 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).
2645 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
2646 The solution is to run a scripting language as a subprocess accessed via pipes.
2647 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
2649 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
2650 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
2651 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.
2652 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
2654 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.
2655 Then you can either statically or dynamically link to those hardcoded functions.
2656 While much less flexible, this approach would be a fairly simple-to-implement fallback.
2658 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
2659 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
2662 #ifndef STRING_EVAL_H
2663 #define STRING_EVAL_H
2665 <<string eval setup declaration>>
2666 <<string eval function declaration>>
2667 <<string eval teardown declaration>>
2668 #endif /* STRING_EVAL_H */
2671 <<string eval module makefile lines>>=
2672 NW_SAWSIM_MODS += string_eval
2673 CHECK_BINS += check_string_eval
2674 check_string_eval_MODS =
2677 For an introduction to POSIX process control, see\\
2678 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
2679 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
2680 and of course, the relavent [[man]] pages.
2682 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
2683 [[execvp]] replaces the calling process' program with a new program.
2684 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
2685 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
2686 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
2687 The new program needs command line arguments, just like it would if you were running it from a shell.
2688 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
2689 with the final array entry being a [[NULL]] pointer.
2691 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
2692 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
2693 <<string eval subprocess definitions>>=
2694 #define SUBPROCESS "python"
2695 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
2696 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
2697 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
2699 The [[i]] option lets Python know that it should run in interactive mode.
2700 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
2701 In interactive mode, python acts on every instruction as soon as it is recieved.
2702 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.
2703 %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.
2705 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
2706 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
2707 [[fork]] returns two copies of the same program, executing the original code.
2708 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
2709 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
2711 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.
2712 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
2713 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
2714 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
2715 <<string eval pipe definitions>>=
2716 #define PIPE_READ 0 /* the end you read from */
2717 #define PIPE_WRITE 1 /* the end you write to */
2719 #define STDIN 0 /* initial index of stdin pair */
2720 #define STDOUT 2 /* initial index of stdout pair */
2723 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
2725 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.
2727 <<string eval setup declaration>>=
2728 extern void string_eval_setup(FILE **pIn, FILE **pOut);
2730 <<string eval setup definition>>=
2731 void string_eval_setup(FILE **pIn, FILE **pOut)
2734 int pfd[4]; /* file descriptors for stdin and stdout */
2736 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
2737 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
2739 pid = fork(); /* split process into two copies */
2741 if (pid == -1) { /* fork error */
2742 perror("fork error");
2744 } else if (pid == 0) { /* child */
2745 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
2746 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
2747 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
2748 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
2749 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
2750 perror("exec error"); /* exec shouldn't return */
2752 } else { /* parent */
2753 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
2754 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
2755 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
2756 if ( *pIn == NULL ) {
2757 perror("fdopen (in)");
2760 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
2761 if ( *pOut == NULL ) {
2762 perror("fdopen (out)");
2769 To use the evaluating subprocess, we just pipe in our command, and read out the results.
2770 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
2771 <<string eval function declaration>>=
2772 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
2774 <<string eval function definition>>=
2775 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
2778 rc = fprintf(in, "%s", input);
2779 assert(rc == strlen(input));
2782 alarm(1); /* set a one second timeout on the read */
2783 assert( fgets(output, buflen, out) != NULL );
2784 alarm(0); /* cancel the timeout */
2785 //fprintf(stderr, "eval: %s ----> %s", input, output);
2788 The [[alarm]] calls set and clear a timeout on the returned output.
2789 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.
2790 This protects against invalid input for which a line of output is not printed to [[stdout]].
2791 Other invalid inputs (e.g. those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
2792 If you are getting strange results, check your python code seperately. TODO, better error handling.
2794 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
2795 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
2796 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
2797 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
2798 <<string eval teardown declaration>>=
2799 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
2802 <<string eval teardown definition>>=
2803 void string_eval_teardown(FILE **pIn, FILE **pOut)
2808 /* redirect Python's stderr */
2809 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
2813 assert( fclose(*pIn) == 0 );
2815 assert( fclose(*pOut) == 0 );
2818 /* wait for python to exit */
2820 pid = wait(&stat_loc);
2827 if (WIFEXITED(stat_loc)) {
2828 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
2829 } else if (WIFSIGNALED(stat_loc)) {
2830 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
2835 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
2837 \subsection{String evaluation implementation}
2841 <<string eval includes>>
2842 #include "string_eval.h"
2843 <<string eval internal definitions>>
2844 <<string eval functions>>
2847 <<string eval includes>>=
2848 #include <assert.h> /* assert() */
2849 #include <stdlib.h> /* NULL */
2850 #include <stdio.h> /* fprintf(), stdout, fdopen() */
2851 #include <string.h> /* strlen() */
2852 #include <sys/types.h> /* pid_t */
2853 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
2854 #include <wait.h> /* wait() */
2857 <<string eval internal definitions>>=
2858 <<string eval subprocess definitions>>
2859 <<string eval pipe definitions>>
2862 <<string eval functions>>=
2863 <<string eval setup definition>>
2864 <<string eval function definition>>
2865 <<string eval teardown definition>>
2868 \subsection{String evaluation unit tests}
2870 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2871 <<check-string-eval.c>>=
2872 <<string eval check includes>>
2873 <<string comparison definition and globals>>
2874 <<string eval test suite>>
2875 <<main check program>>
2878 <<string eval check includes>>=
2879 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
2880 #include <stdio.h> /* printf() */
2881 #include <assert.h> /* assert() */
2882 #include <string.h> /* strlen() */
2883 #include <signal.h> /* SIGKILL */
2885 #include "string_eval.h"
2888 <<string eval test suite>>=
2889 <<string eval tests>>
2890 <<string eval suite definition>>
2893 <<string eval suite definition>>=
2894 Suite *test_suite (void)
2896 Suite *s = suite_create ("string eval");
2897 <<string eval test case defs>>
2899 <<string eval test case add>>
2904 <<string eval tests>>=
2905 START_TEST(test_setup_teardown)
2908 string_eval_setup(&in, &out);
2909 string_eval_teardown(&in, &out);
2912 START_TEST(test_invalid_command)
2915 char input[80], output[80]={};
2916 string_eval_setup(&in, &out);
2917 sprintf(input, "print ABCDefg\n");
2918 string_eval(in, out, input, 80, output);
2919 string_eval_teardown(&in, &out);
2922 START_TEST(test_simple_eval)
2925 char input[80], output[80]={};
2926 string_eval_setup(&in, &out);
2927 sprintf(input, "print 3+4*5\n");
2928 string_eval(in, out, input, 80, output);
2929 fail_unless(STRMATCH(output,"23\n"), NULL);
2930 string_eval_teardown(&in, &out);
2933 START_TEST(test_multiple_evals)
2936 char input[80], output[80]={};
2937 string_eval_setup(&in, &out);
2938 sprintf(input, "print 3+4*5\n");
2939 string_eval(in, out, input, 80, output);
2940 fail_unless(STRMATCH(output,"23\n"), NULL);
2941 sprintf(input, "print (3**2 + 4**2)**0.5\n");
2942 string_eval(in, out, input, 80, output);
2943 fail_unless(STRMATCH(output,"5.0\n"), NULL);
2944 string_eval_teardown(&in, &out);
2947 START_TEST(test_eval_with_state)
2950 char input[80], output[80]={};
2951 string_eval_setup(&in, &out);
2952 sprintf(input, "print 3+4*5\n");
2953 fprintf(in, "A = 3\n");
2954 sprintf(input, "print A*3\n");
2955 string_eval(in, out, input, 80, output);
2956 fail_unless(STRMATCH(output,"9\n"), NULL);
2957 string_eval_teardown(&in, &out);
2961 <<string eval test case defs>>=
2962 TCase *tc_string_eval = tcase_create("string_eval");
2964 <<string eval test case add>>=
2965 tcase_add_test(tc_string_eval, test_setup_teardown);
2966 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
2967 tcase_add_test(tc_string_eval, test_simple_eval);
2968 tcase_add_test(tc_string_eval, test_multiple_evals);
2969 tcase_add_test(tc_string_eval, test_eval_with_state);
2970 suite_add_tcase(s, tc_string_eval);
2974 \section{Accelerating function evaluation}
2976 My first version-0.3 code was running very slowly.
2977 With the options suggested in the help
2978 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
2979 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
2980 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
2981 That's only 15 calls per solution, so the search algorithm seems reasonable.
2982 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
2987 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
2989 #endif /* ACCEL_K_H */
2992 <<accel k module makefile lines>>=
2993 NW_SAWSIM_MODS += accel_k
2994 #CHECK_BINS += check_accel_k
2995 check_accel_k_MODS =
2999 #include <assert.h> /* assert() */
3000 #include <stdlib.h> /* realloc(), free(), NULL */
3001 #include "global.h" /* environment_t */
3002 #include "k_model.h" /* k_func_t */
3003 #include "interp.h" /* interp_* */
3004 #include "accel_k.h"
3006 typedef struct accel_k_struct {
3007 interp_table_t *itable;
3013 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3014 static int num_accels = 0;
3015 static accel_k_t *accels=NULL;
3017 /* Wrap k in the standard f(x) acceleration form */
3018 static double k_wrap(double F, void *params)
3020 accel_k_t *p = (accel_k_t *)params;
3021 return (*p->k)(F, p->env, p->params);
3024 static int k_tol(double FA, double kA, double FB, double kB)
3027 if (FB-FA > 1e-12) {
3028 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3029 return 1; /* unnacceptable */
3031 //printf("acceptable tol\n");
3032 return 0; /* acceptable */
3036 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3039 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3040 assert(accels != NULL);
3041 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3043 accels[i].env = env;
3044 accels[i].params = params;
3051 for (i=0; i<num_accels; i++)
3052 interp_table_free(accels[i].itable);
3054 if (accels) free(accels);
3058 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3061 for (i=0; i<num_accels; i++) {
3062 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3063 /* We've already setup this function.
3064 * WARNING: we're assuming param and env are the same. */
3065 return interp_table_eval(accels[i].itable, accels+i, F);
3068 /* set up a new acceleration environment */
3069 i = add_accel_k(k, env, params);
3070 return interp_table_eval(accels[i].itable, accels+i, F);
3074 \section{Tension models}
3075 \label{sec.tension_models}
3077 TODO: tension model intro.
3078 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.
3080 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3081 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]].
3083 <<tension-model.h>>=
3084 #ifndef TENSION_MODEL_H
3085 #define TENSION_MODEL_H
3087 <<tension handler types>>
3088 <<constant tension model declarations>>
3089 <<hooke tension model declarations>>
3090 <<worm-like chain tension model declarations>>
3091 <<freely-jointed chain tension model declarations>>
3092 <<find tension definitions>>
3093 <<tension model global declarations>>
3094 #endif /* TENSION_MODEL_H */
3097 <<tension model module makefile lines>>=
3098 NW_SAWSIM_MODS += tension_model
3099 #CHECK_BINS += check_tension_model
3100 check_tension_model_MODS =
3102 <<tension model utils makefile lines>>=
3103 TENSION_MODEL_MODS = tension_model parse list tension_balance
3104 TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
3105 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \
3106 $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3107 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw $(BUILD_DIR)
3108 notangle -Rtension-model-utils.c $< > $@
3109 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) $(BIN_DIR)
3110 gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3111 clean_tension_model_utils :
3112 rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
3116 \label{sec.null_tension}
3118 For unstretchable domains.
3120 <<null tension model getopt>>=
3121 {NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3124 \subsection{Constant}
3125 \label{sec.const_tension}
3127 <<constant tension functions>>=
3128 <<constant tension function>>
3129 <<constant tension parameter create/destroy functions>>
3132 <<constant tension model declarations>>=
3133 <<constant tension function declaration>>
3134 <<constant tension parameter create/destroy function declarations>>
3135 <<constant tension model global declarations>>
3139 An infinitely stretchable domain providing a constant tension.
3141 <<constant tension function declaration>>=
3142 extern double const_tension_handler(double x, void *pdata);
3144 <<constant tension function>>=
3145 double const_tension_handler(double x, void *pdata)
3147 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3148 list_t *list = data->group;
3153 assert(list != NULL); /* empty active group?! */
3154 F = ((const_tension_param_t *)list->d)->F;
3156 fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3158 while (list != NULL) {
3159 assert(((const_tension_param_t *)list->d)->F == F);
3167 <<constant tension structure>>=
3168 typedef struct const_tension_param_struct {
3169 double F; /* tension (force) in N */
3170 } const_tension_param_t;
3174 <<constant tension parameter create/destroy function declarations>>=
3175 extern void *string_create_const_tension_param_t(char **param_strings);
3176 extern void destroy_const_tension_param_t(void *p);
3179 <<constant tension parameter create/destroy functions>>=
3180 const_tension_param_t *create_const_tension_param_t(double F)
3182 const_tension_param_t *ret
3183 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3184 assert(ret != NULL);
3189 void *string_create_const_tension_param_t(char **param_strings)
3191 return create_const_tension_param_t(atof(param_strings[0]));
3194 void destroy_const_tension_param_t(void *p)
3202 <<constant tension model global declarations>>=
3203 extern int num_const_tension_params;
3204 extern char *const_tension_param_descriptions[];
3205 extern char const_tension_param_string[];
3208 <<constant tension model globals>>=
3209 int num_const_tension_params = 1;
3210 char *const_tension_param_descriptions[] = {"tension F, N"};
3211 char const_tension_param_string[] = "0";
3214 <<constant tension model getopt>>=
3215 {&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}
3221 <<hooke functions>>=
3223 <<hooke parameter create/destroy functions>>
3226 <<hooke tension model declarations>>=
3227 <<hooke tension function declaration>>
3228 <<hooke tension parameter create/destroy function declarations>>
3229 <<hooke tension model global declarations>>
3233 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
3234 The behavior of a series of springs $k_i$ in series is given by
3236 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
3238 For a simple proof, see Appendix \ref{app.math_hooke}.
3240 <<hooke tension function declaration>>=
3241 extern double hooke_handler(double x, void *pdata);
3245 double hooke_handler(double x, void *pdata)
3247 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3248 list_t *list = data->group;
3253 assert(list != NULL); /* empty active group?! */
3254 while (list != NULL) {
3255 assert( ((hooke_param_t *)list->d)->k > 0 );
3256 k += 1.0/ ((hooke_param_t *)list->d)->k;
3258 fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
3259 __FUNCTION__, ((hooke_param_t *)list->d)->k);
3265 fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N\n",
3266 __FUNCTION__, k, x, k*x);
3272 <<hooke structure>>=
3273 typedef struct hooke_param_struct {
3274 double k; /* spring constant in N/m */
3278 <<hooke tension parameter create/destroy function declarations>>=
3279 extern void *string_create_hooke_param_t(char **param_strings);
3280 extern void destroy_hooke_param_t(void *p);
3283 <<hooke parameter create/destroy functions>>=
3284 hooke_param_t *create_hooke_param_t(double k)
3286 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
3287 assert(ret != NULL);
3292 void *string_create_hooke_param_t(char **param_strings)
3294 return create_hooke_param_t(atof(param_strings[0]));
3297 void destroy_hooke_param_t(void *p)
3304 <<hooke tension model global declarations>>=
3305 extern int num_hooke_params;
3306 extern char *hooke_param_descriptions[];
3307 extern char hooke_param_string[];
3310 <<hooke tension model globals>>=
3311 int num_hooke_params = 1;
3312 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
3313 char hooke_param_string[]="0.05";
3316 <<hooke tension model getopt>>=
3317 {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}
3320 \subsection{Worm-like chain}
3323 We can model several unfolded domains with a single worm-like chain.
3324 <<worm-like chain functions>>=
3325 <<worm-like chain function>>
3326 <<worm-like chain handler>>
3327 <<worm-like chain create/destroy functions>>
3330 <<worm-like chain tension model declarations>>=
3331 <<worm-like chain handler declaration>>
3332 <<worm-like chain create/destroy function declarations>>
3333 <<worm-like chain tension model global declarations>>
3337 The combination of all unfolded domains is modeled as a worm like chain,
3338 with the force $F_{\text{WLC}}$ approximately given by
3340 F_{\text{WLC}} \approx \frac{k_B T}{p}
3341 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
3343 where $x$ is the distance between the fixed ends,
3344 $k_B$ is Boltzmann's constant,
3345 $T$ is the absolute temperature,
3346 $p$ is the persistence length, and
3347 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
3348 <<worm-like chain function>>=
3349 double wlc(double x, double T, double p, double L)
3351 double xL=0.0; /* uses global k_B */
3353 assert(T > 0); assert(p > 0); assert(L > 0);
3354 if (x >= L) return HUGE_VAL;
3355 xL = x/L; /* unitless */
3356 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
3357 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
3358 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
3361 This model assumes that each unfolded domain has the same persistence length.
3363 <<worm-like chain handler declaration>>=
3364 extern double wlc_handler(double x, void *pdata);
3367 <<worm-like chain handler>>=
3368 double wlc_handler(double x, void *pdata)
3370 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3371 list_t *list = data->group;
3375 assert(list != NULL); /* empty active group?! */
3376 p = ((wlc_param_t *)list->d)->p;
3377 while (list != NULL) {
3378 /* enforce consistent p */
3379 assert( ((wlc_param_t *)list->d)->p == p);
3380 L += ((wlc_param_t *)list->d)->L;
3382 fprintf(stderr, "%s: WLC section %g m long in series\n",
3383 __FUNCTION__, ((wlc_param_t *)list->d)->L);
3388 fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
3389 __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
3391 return wlc(x, data->env->T, p, L);
3395 <<worm-like chain structure>>=
3396 typedef struct wlc_param_struct {
3397 double p; /* WLC persistence length (m) */
3398 double L; /* WLC contour length (m) */
3402 <<worm-like chain create/destroy function declarations>>=
3403 extern void *string_create_wlc_param_t(char **param_strings);
3404 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
3407 <<worm-like chain create/destroy functions>>=
3408 wlc_param_t *create_wlc_param_t(double p, double L)
3410 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
3411 assert(ret != NULL);
3417 void *string_create_wlc_param_t(char **param_strings)
3419 return create_wlc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3422 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
3429 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.
3430 TODO (now we copy the string before parsing, but still do this for clarity).
3432 <<valid string write code>>=
3433 char string[] = "abc";
3436 <<valid string write code>>=
3437 char *string = "abc";
3441 <<worm-like chain tension model global declarations>>=
3442 extern int num_wlc_params;
3443 extern char *wlc_param_descriptions[];
3444 extern char wlc_param_string[];
3446 <<worm-like chain tension model globals>>=
3447 int num_wlc_params = 2;
3448 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
3449 char wlc_param_string[]="0.39e-9,28.4e-9";
3452 <<worm-like chain tension model getopt>>=
3453 {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}
3455 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
3457 \subsection{Freely-jointed chain}
3460 <<freely-jointed chain functions>>=
3461 <<freely-jointed chain function>>
3462 <<freely-jointed chain handler>>
3463 <<freely-jointed chain create/destroy functions>>
3466 <<freely-jointed chain tension model declarations>>=
3467 <<freely-jointed chain handler declaration>>
3468 <<freely-jointed chain create/destroy function declarations>>
3469 <<freely-jointed chain tension model global declarations>>
3473 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
3474 The tension of a chain of $N$ such links is given by
3476 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
3478 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}.
3479 <<freely-jointed chain function>>=
3480 double langevin(double x, void *params)
3483 return 1.0/tanh(x) - 1/x;
3485 <<one dimensional bracket declaration>>
3486 <<one dimensional solve declaration>>
3487 double inv_langevin(double x)
3489 double lb=0.0, ub=1.0, ret;
3490 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
3491 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
3495 double fjc(double x, double T, double l, int N)
3497 double L = l*(double)N;
3498 if (x == 0) return 0;
3499 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
3500 return k_B*T/l * inv_langevin(x/L);
3504 <<freely-jointed chain handler declaration>>=
3505 extern double fjc_handler(double x, void *pdata);
3507 <<freely-jointed chain handler>>=
3508 double fjc_handler(double x, void *pdata)
3510 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3511 list_t *list = data->group;
3516 assert(list != NULL); /* empty active group?! */
3517 l = ((fjc_param_t *)list->d)->l;
3518 while (list != NULL) {
3519 /* enforce consistent l */
3520 assert( ((fjc_param_t *)list->d)->l == l);
3521 N += ((fjc_param_t *)list->d)->N;
3523 fprintf(stderr, "%s: FJC section %d links long in series\n",
3524 __FUNCTION__, ((fjc_param_t *)list->d)->N);
3529 fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
3530 __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
3532 return fjc(x, data->env->T, l, N);
3536 <<freely-jointed chain structure>>=
3537 typedef struct fjc_param_struct {
3538 double l; /* FJC link length (m) */
3539 int N; /* FJC number of links */
3543 <<freely-jointed chain create/destroy function declarations>>=
3544 extern void *string_create_fjc_param_t(char **param_strings);
3545 extern void destroy_fjc_param_t(void *p);
3548 <<freely-jointed chain create/destroy functions>>=
3549 fjc_param_t *create_fjc_param_t(double l, double N)
3551 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
3552 assert(ret != NULL);
3558 void *string_create_fjc_param_t(char **param_strings)
3560 return create_fjc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3563 void destroy_fjc_param_t(void *p)
3570 <<freely-jointed chain tension model global declarations>>=
3571 extern int num_fjc_params;
3572 extern char *fjc_param_descriptions[];
3573 extern char fjc_param_string[];
3576 <<freely-jointed chain tension model globals>>=
3577 int num_fjc_params = 2;
3578 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
3579 char fjc_param_string[]="0.5e-9,200";
3582 <<freely-jointed chain tension model getopt>>=
3583 {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}
3586 \subsection{Tension model implementation}
3588 <<tension-model.c>>=
3591 <<tension model includes>>
3592 #include "tension_model.h"
3593 <<tension model internal definitions>>
3594 <<tension model globals>>
3595 <<tension model functions>>
3598 <<tension model includes>>=
3599 #include <assert.h> /* assert() */
3600 #include <stdlib.h> /* NULL */
3601 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
3602 #include <stdio.h> /* fprintf(), stdout */
3605 #include "tension_balance.h" /* oneD_*() */
3608 <<tension model internal definitions>>=
3609 <<constant tension structure>>
3611 <<worm-like chain structure>>
3612 <<freely-jointed chain structure>>
3615 <<tension model globals>>=
3616 <<tension handler array global>>
3617 <<constant tension model globals>>
3618 <<hooke tension model globals>>
3619 <<worm-like chain tension model globals>>
3620 <<freely-jointed chain tension model globals>>
3623 <<tension model functions>>=
3624 <<constant tension functions>>
3626 <<worm-like chain functions>>
3627 <<freely-jointed chain functions>>
3631 \subsection{Utilities}
3633 The tension models can get complicated, and users may want to reassure
3634 themselves that this portion of the input physics are functioning properly. The
3635 stand-alone program [[tension_model_utils]] demonstrates the tension model
3636 interface without the overhead of having to understand a full simulation.
3637 [[tension_model_utils]] takes command line model arguments like the full
3638 simulation, and prints $F(x)$ data to the screen for the selected model over a
3641 <<tension-model-utils.c>>=
3643 <<tension model utility includes>>
3644 <<tension model utility definitions>>
3645 <<tension model utility getopt functions>>
3647 <<tension model utility main>>
3650 <<tension model utility main>>=
3651 int main(int argc, char **argv)
3653 <<tension model getopt array definition>>
3654 tension_model_getopt_t *model = NULL;
3658 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
3659 tension_handler_data_t tdata;
3660 int num_param_args; /* for INIT_MODEL() */
3661 char **param_args; /* for INIT_MODEL() */
3662 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model, &flags);
3664 if (flags & VFLAG) {
3665 printf("#initializing model %s with parameters %s\n", model->name, model->params);
3667 INIT_MODEL("utility", model, model->params, params);
3669 push(&tdata.group, params, 1);
3671 tdata.persist = NULL;
3672 if (model->handler == NULL) {
3673 printf("No tension function for model %s\n", model->name);
3677 double dx=1e-10, x=0, F=0;
3678 printf("#F (N)\tk (%% pop. per s)\n");
3679 while (F >= 0 && F < 1e5 && x < 1e-6) {
3680 F = (*model->handler)(x, &tdata);
3681 printf("%g\t%g\n", x, F);
3685 params = pop(&tdata.group);
3686 if (model->destructor)
3687 (*model->destructor)(params);
3692 <<tension model utility includes>>=
3695 #include <string.h> /* strlen() */
3696 #include <assert.h> /* assert() */
3700 #include "tension_model.h"
3703 <<tension model utility definitions>>=
3704 <<version definition>>
3705 #define VFLAG 1 /* verbose */
3706 <<string comparison definition and globals>>
3707 <<tension model getopt definitions>>
3708 <<initialize model definition>>
3712 <<tension model utility getopt functions>>=
3713 <<index tension model>>
3714 <<tension model utility help>>
3715 <<tension model utility get options>>
3718 <<tension model utility help>>=
3719 void help(char *prog_name,
3721 int n_tension_models, tension_model_getopt_t *tension_models,
3725 printf("usage: %s [options]\n", prog_name);
3726 printf("Version %s\n\n", VERSION);
3727 printf("A utility for understanding the available tension models\n\n");
3728 printf("Environmental options:\n");
3729 printf("-T T\tTemperature T (currently %g K)\n", env->T);
3730 printf("-C T\tYou can also set the temperature T in Celsius\n");
3731 printf("Model options:\n");
3732 printf("-m model\tFolded tension model (currently %s)\n",
3733 tension_models[tension_model].name);
3734 printf("-a args\tFolded tension model argument string (currently %s)\n",
3735 tension_models[tension_model].params);
3736 printf("F(x) is calculated for a range of x and printed\n");
3737 printf("For example:\n");
3738 printf(" #Distance (x)\tForce (N)\n");
3739 printf(" 123.456\t7.89\n");
3741 printf("-V\tChange output to verbose mode\n");
3742 printf("-h\tPrint this help and exit\n");
3744 printf("Tension models:\n");
3745 for (i=0; i<n_tension_models; i++) {
3746 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
3747 for (j=0; j < tension_models[i].num_params ; j++)
3748 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
3749 printf("\t\tdefaults: %s\n", tension_models[i].params);
3754 <<tension model utility get options>>=
3755 void get_options(int argc, char **argv, environment_t *env,
3756 int n_tension_models, tension_model_getopt_t *tension_models,
3757 tension_model_getopt_t **model,
3758 unsigned int *flags)
3760 char *prog_name = NULL;
3761 char c, options[] = "T:C:m:a:Vh";
3762 int tension_model=0, num_strings;
3763 extern char *optarg;
3764 extern int optind, optopt, opterr;
3768 /* setup defaults */
3770 prog_name = argv[0];
3771 env->T = 300.0; /* K */
3773 *model = tension_models;
3775 while ((c=getopt(argc, argv, options)) != -1) {
3777 case 'T': env->T = atof(optarg); break;
3778 case 'C': env->T = atof(optarg)+273.15; break;
3780 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
3781 *model = tension_models+tension_model;
3784 tension_models[tension_model].params = optarg;
3786 case 'V': *flags |= VFLAG; break;
3788 fprintf(stderr, "unrecognized option '%c'\n", optopt);
3789 /* fall through to default case */
3791 help(prog_name, env, n_tension_models, tension_models, tension_model);
3800 \section{Unfolding rate models}
3801 \label{sec.k_models}
3803 We have two main choices for determining $k$: Bell's model, which assumes the
3804 folded domains exist in a single `folded' state and make instantaneous,
3805 stochastic transitions over a free energy barrier; and Kramers' model, which
3806 treats the unfolding as a thermalized diffusion process.
3807 We deal with these options like we dealt with the different tension models: by bundling all model-specific
3808 parameters into structures, with model specific functions for determining $k$.
3810 <<k func definition>>=
3811 typedef double k_func_t(double F, environment_t *env, void *params);
3814 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3815 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]].
3817 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
3818 because \LaTeX\ doesn't like underscores outside math mode.
3819 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
3820 You could use spaces instead (`\verb+ +'),
3821 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
3822 which I don't like as much.
3828 <<k func definition>>
3829 <<null k declarations>>
3830 <<const k declarations>>
3831 <<bell k declarations>>
3832 <<kramers k declarations>>
3833 <<kramers integ k declarations>>
3834 #endif /* K_MODEL_H */
3837 <<k model module makefile lines>>=
3838 NW_SAWSIM_MODS += k_model
3839 CHECK_BINS += check_k_model
3840 check_k_model_MODS = parse string_eval
3842 [[check_k_model]] also depends on the parse module.
3844 <<k model utils makefile lines>>=
3845 K_MODEL_MODS = k_model parse string_eval
3846 K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
3847 $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h)
3848 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw $(BUILD)
3849 notangle -Rk-model-utils.c $< > $@
3850 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) $(BIN_DIR)
3851 gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
3852 clean_k_model_utils :
3853 rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
3859 For domains that are never folded.
3861 <<null k declarations>>=
3863 <<null k functions>>=
3868 <<null k model getopt>>=
3869 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
3872 \subsection{Constant rate model}
3875 This is a pretty boring model, but it is usefull for testing the engine.
3879 where $k_0$ is the constant unfolding rate.
3881 <<const k functions>>=
3882 <<const k function>>
3883 <<const k structure create/destroy functions>>
3886 <<const k declarations>>=
3887 <<const k function declaration>>
3888 <<const k structure create/destroy function declarations>>
3889 <<const k global declarations>>
3893 <<const k structure>>=
3894 typedef struct const_k_param_struct {
3895 double knot; /* transition rate k_0 (frac population per s) */
3899 <<const k function declaration>>=
3900 double const_k(double F, environment_t *env, void *const_k_params);
3902 <<const k function>>=
3903 double const_k(double F, environment_t *env, void *const_k_params)
3904 { /* Returns the rate constant k in frac pop per s. */
3905 const_k_param_t *p = (const_k_param_t *) const_k_params;
3907 assert(p->knot > 0);
3912 <<const k structure create/destroy function declarations>>=
3913 void *string_create_const_k_param_t(char **param_strings);
3914 void destroy_const_k_param_t(void *p);
3917 <<const k structure create/destroy functions>>=
3918 const_k_param_t *create_const_k_param_t(double knot)
3920 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
3921 assert(ret != NULL);
3926 void *string_create_const_k_param_t(char **param_strings)
3928 return create_const_k_param_t(atof(param_strings[0]));
3931 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
3938 <<const k global declarations>>=
3939 extern int num_const_k_params;
3940 extern char *const_k_param_descriptions[];
3941 extern char const_k_param_string[];
3944 <<const k globals>>=
3945 int num_const_k_params = 1;
3946 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
3947 char const_k_param_string[]="1";
3950 <<const k model getopt>>=
3951 {"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}
3954 \subsection{Bell's model}
3957 Quantitatively, Bell's model gives $k$ as
3959 k = k_0 \cdot e^{F dx / k_B T} \;,
3961 where $k_0$ is the unforced unfolding rate,
3962 $F$ is the applied unfolding force,
3963 $dx$ is the distance to the transition state, and
3964 $k_B T$ is the thermal energy\citep{schlierf06}.
3966 <<bell k functions>>=
3968 <<bell k structure create/destroy functions>>
3971 <<bell k declarations>>=
3972 <<bell k function declaration>>
3973 <<bell k structure create/destroy function declarations>>
3974 <<bell k global declarations>>
3978 <<bell k structure>>=
3979 typedef struct bell_param_struct {
3980 double knot; /* transition rate k_0 (frac population per s) */
3981 double dx; /* `distance to transition state' (nm) */
3985 <<bell k function declaration>>=
3986 double bell_k(double F, environment_t *env, void *bell_params);
3988 <<bell k function>>=
3989 double bell_k(double F, environment_t *env, void *bell_params)
3990 { /* Returns the rate constant k in frac pop per s.
3991 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
3992 * uses global k_B in J/K */
3993 bell_param_t *p = (bell_param_t *) bell_params;
3994 assert(F >= 0); assert(env->T > 0);
3996 assert(p->knot > 0); assert(p->dx >= 0);
3998 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
3999 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
4000 p->knot * exp(F*p->dx / (k_B*env->T) ));
4001 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
4003 return p->knot * exp(F*p->dx / (k_B*env->T) );
4007 <<bell k structure create/destroy function declarations>>=
4008 void *string_create_bell_param_t(char **param_strings);
4009 void destroy_bell_param_t(void *p);
4012 <<bell k structure create/destroy functions>>=
4013 bell_param_t *create_bell_param_t(double knot, double dx)
4015 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
4016 assert(ret != NULL);
4022 void *string_create_bell_param_t(char **param_strings)
4024 return create_bell_param_t(atof(param_strings[0]), atof(param_strings[1]));
4027 void destroy_bell_param_t(void *p /* bell_param_t * */)
4034 <<bell k global declarations>>=
4035 extern int num_bell_params;
4036 extern char *bell_param_descriptions[];
4037 extern char bell_param_string[];
4041 int num_bell_params = 2;
4042 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4043 char bell_param_string[]="3.3e-4,0.25e-9";
4046 <<bell k model getopt>>=
4047 {"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}
4049 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4050 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4052 \subsection{Kramer's model}
4055 Kramer's model gives $k$ as
4057 % \frac{1}{k} = \frac{1}{D}
4058 % \int_{x_\text{min}}^{x_\text{max}}
4059 % e^{\frac{-U_F(x)}{k_B T}}
4060 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4063 %where $D$ is the diffusion constant,
4064 %$U_F(x)$ is the free energy along the unfolding coordinate, and
4065 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4066 % before the minimum and after the maximum, respectively \citep{schlierf06}.
4068 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
4070 where $D$ is the diffusion constant,
4072 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
4074 are length scales where
4075 $x_c(F)$ is the location of the bound state and
4076 $x_{ts}(F)$ is the location of the transition state,
4077 $E(F,x)$ is the free energy, and
4078 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
4079 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
4080 \citet{evans97} Eqn.~3).
4082 <<kramers k functions>>=
4083 <<kramers k function>>
4084 <<kramers k structure create/destroy functions>>
4087 <<kramers k declarations>>=
4088 <<kramers k function declaration>>
4089 <<kramers k structure create/destroy function declarations>>
4090 <<kramers k global declarations>>
4094 <<kramers k structure>>=
4095 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
4096 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
4098 typedef struct kramers_param_struct {
4099 double D; /* diffusion rate (um/s) */
4106 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
4107 //kramers_x_func_t *xts; /* function returning transition x (nm) */
4108 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
4109 //kramers_E_func_t *E; /* function returning E (J) */
4110 //double *E_params; /* parametrize the energy landscape */
4111 //int n_E_params; /* length of E_params array */
4115 <<kramers k function declaration>>=
4116 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
4117 extern double kramers_k(double F, environment_t *env, void *kramers_params);
4119 <<kramers k function>>=
4120 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
4122 static char input[160]={0};
4123 static char output[80]={0};
4124 /* setup the environment */
4125 fprintf(in, "F = %g; T = %g\n", F, T);
4126 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4127 string_eval(in, out, input, 80, output);
4128 return atof(output);
4131 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
4133 static char input[160]={0};
4134 static char output[80]={0};
4135 /* setup the environment */
4136 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
4137 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4138 string_eval(in, out, input, 80, output);
4139 return atof(output);
4142 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
4144 kramers_param_t *p = (kramers_param_t *) kramers_params;
4145 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
4148 double kramers_k(double F, environment_t *env, void *kramers_params)
4149 { /* Returns the rate constant k in frac pop per s.
4150 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4151 * uses global k_B in J/K */
4152 kramers_param_t *p = (kramers_param_t *) kramers_params;
4153 double kbT, xc, xts, lc, lts, Eb;
4154 assert(F >= 0); assert(env->T > 0);
4157 assert(p->xc != NULL); assert(p->xts != NULL);
4158 assert(p->ddEdxx != NULL); assert(p->E != NULL);
4159 assert(p->in != NULL); assert(p->out != NULL);
4161 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
4162 if (xc == -1.0) return -1.0; /* error (Too much force) */
4163 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
4164 if (xc == -1.0) return -1.0; /* error (Too much force) */
4165 lc = sqrt(2.0*M_PI*kbT /
4166 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
4167 lts = sqrt(-2.0*M_PI*kbT/
4168 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
4169 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
4170 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
4171 //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));
4172 return p->D/(lc*lts) * exp(-Eb/kbT);
4176 <<kramers k structure create/destroy function declarations>>=
4177 void *string_create_kramers_param_t(char **param_strings);
4178 void destroy_kramers_param_t(void *p);
4181 <<kramers k structure create/destroy functions>>=
4182 kramers_param_t *create_kramers_param_t(double D,
4183 char *xc_fn, char *xts_fn,
4184 char *E_fn, char *ddEdxx_fn,
4185 char *global_define_string)
4186 // kramers_x_func_t *xc, kramers_x_func_t *xts,
4187 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
4188 // double *E_params, int n_E_params)
4190 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
4191 assert(ret != NULL);
4192 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4193 assert(ret->xc != NULL);
4194 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4195 assert(ret->xts != NULL);
4196 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
4197 assert(ret->E != NULL);
4198 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
4199 assert(ret->ddEdxx != NULL);
4201 strcpy(ret->xc, xc_fn);
4202 strcpy(ret->xts, xts_fn);
4203 strcpy(ret->E, E_fn);
4204 strcpy(ret->ddEdxx, ddEdxx_fn);
4205 //ret->E_params = E_params;
4206 //ret->n_E_params = n_E_params;
4207 string_eval_setup(&ret->in, &ret->out);
4208 fprintf(ret->in, "kB = %g\n", k_B);
4209 fprintf(ret->in, "%s\n", global_define_string);
4213 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
4214 void *string_create_kramers_param_t(char **param_strings)
4216 return create_kramers_param_t(atof(param_strings[0]),
4224 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
4226 kramers_param_t *pk = (kramers_param_t *)p;
4228 string_eval_teardown(&pk->in, &pk->out);
4230 // free(pk->E_params);
4236 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
4237 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.
4238 We follow \citet{shillcock98} and use
4240 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
4241 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
4244 where TODO, variable meanings.%\citep{shillcock98}.
4245 The first and second derivatives of $E(F,E_b,x,x_s)$ are
4247 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\\
4248 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
4250 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
4251 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
4252 The roots are therefor at
4254 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
4255 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
4256 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
4259 <<kramers k global declarations>>=
4260 extern int num_kramers_params;
4261 extern char *kramers_param_descriptions[];
4262 extern char kramers_param_string[];
4265 <<kramers k globals>>=
4266 int num_kramers_params = 6;
4267 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"};
4268 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)}";
4270 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
4271 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
4272 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
4273 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.
4274 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
4275 It works so long as [[val_1]] is not [[false]].
4277 <<kramers k model getopt>>=
4278 {"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}
4281 \subsection{Kramer's model (natural cubic splines)}
4282 \label{sec.kramers_integ}
4284 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
4285 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
4286 \citet{schlierf06} Eqn.~2)
4288 \frac{1}{k} = \frac{1}{D}
4289 \int_{x_\text{min}}^{x_\text{max}}
4290 e^{\frac{U_F(x)}{k_B T}}
4291 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4294 where $D$ is the diffusion constant,
4295 $U_F(x)$ is the free energy along the unfolding coordinate, and
4296 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4297 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
4299 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
4300 interpolating between them with cubic splines.
4301 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
4303 <<kramers integ k functions>>=
4304 <<spline functions>>
4305 <<kramers integ A k functions>>
4306 <<kramers integ B k functions>>
4307 <<kramers integ k function>>
4308 <<kramers integ k structure create/destroy functions>>
4311 <<kramers integ k declarations>>=
4312 <<kramers integ k function declaration>>
4313 <<kramers integ k structure create/destroy function declarations>>
4314 <<kramers integ k global declarations>>
4318 <<kramers integ k structure>>=
4319 typedef double func_t(double x, void *params);
4320 typedef struct kramers_integ_param_struct {
4321 double D; /* diffusion rate (um/s) */
4322 double x_min; /* integration bounds */
4324 func_t *E; /* function returning E (J) */
4325 void *E_params; /* parametrize the energy landscape */
4326 destroy_data_func_t *destroy_E_params;
4328 double F; /* for passing into GSL evaluations */
4330 } kramers_integ_param_t;
4333 <<spline param structure>>=
4334 typedef struct spline_param_struct {
4335 int n_knots; /* natural cubic spline knots for U(x) */
4338 gsl_spline *spline; /* GSL spline parameters */
4339 gsl_interp_accel *acc; /* GSL spline acceleration data */
4343 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
4345 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4347 <<kramers integ A k functions>>=
4348 double kramers_integ_k_integrandA(double x, void *params)
4350 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4351 double U = (*p->E)(x, p->E_params);
4352 double Fx = p->F * x;
4353 double kbT = k_B * p->env->T;
4354 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
4355 return exp(-(U-Fx)/kbT);
4357 double kramers_integ_k_integralA(double x, void *params)
4359 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4361 double result, abserr;
4362 size_t neval; /* for qng */
4363 static gsl_integration_workspace *w = NULL;
4364 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4365 f.function = &kramers_integ_k_integrandA;
4367 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4368 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4369 //fprintf(stderr, "integralA = %g\n", result);
4375 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
4376 \int_{x_\text{min}}^{x_\text{max}}
4377 e^{\frac{U_F(x)}{k_B T}}
4378 \text{Integral}_A(x)
4381 <<kramers integ B k functions>>=
4382 double kramers_integ_k_integrandB(double x, void *params)
4384 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4385 double U = (*p->E)(x, p->E_params);
4386 double Fx = p->F * x;
4387 double kbT = k_B * p->env->T;
4388 double intA = kramers_integ_k_integralA(x, params);
4389 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
4390 return exp((U-Fx)/kbT)*intA;
4392 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
4394 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4396 double result, abserr;
4397 size_t neval; /* for qng */
4398 static gsl_integration_workspace *w = NULL;
4399 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4400 f.function = &kramers_integ_k_integrandB;
4404 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4405 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4406 //fprintf(stderr, "integralB = %g\n", result);
4411 With the integrals taken care of, evaluating $k$ is simply
4413 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
4415 <<kramers integ k function declaration>>=
4416 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
4418 <<kramers integ k function>>=
4419 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
4420 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
4421 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
4422 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
4426 <<kramers integ k structure create/destroy function declarations>>=
4427 void *string_create_kramers_integ_param_t(char **param_strings);
4428 void destroy_kramers_integ_param_t(void *p);
4431 <<kramers integ k structure create/destroy functions>>=
4432 kramers_integ_param_t *create_kramers_integ_param_t(double D,
4433 double xmin, double xmax, func_t *E,
4435 destroy_data_func_t *destroy_E_params)
4437 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
4438 assert(ret != NULL);
4443 ret->E_params = E_params;
4444 ret->destroy_E_params = destroy_E_params;
4448 void *string_create_kramers_integ_param_t(char **param_strings)
4450 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
4451 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
4452 void *E_params = string_create_spline_param_t(param_strings+1);
4453 return create_kramers_integ_param_t(atof(param_strings[0]),
4454 atof(param_strings[2]),
4455 atof(param_strings[3]),
4456 &spline_eval, E_params,
4457 destroy_spline_param_t);
4460 void destroy_kramers_integ_param_t(void *params)
4462 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4465 (*p->destroy_E_params)(p->E_params);
4471 Finally we have the GSL spline wrappers:
4472 <<spline functions>>=
4474 <<create/destroy spline>>
4477 <<spline function>>=
4478 double spline_eval(double x, void *spline_params)
4480 spline_param_t *p = (spline_param_t *)(spline_params);
4481 return gsl_spline_eval(p->spline, x, p->acc);
4485 <<create/destroy spline>>=
4486 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
4488 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
4489 assert(ret != NULL);
4490 ret->n_knots = n_knots;
4493 ret->acc = gsl_interp_accel_alloc();
4494 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
4495 gsl_spline_init(ret->spline, x, y, n_knots);
4499 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
4500 void *string_create_spline_param_t(char **param_strings)
4504 double *x=NULL, *y=NULL;
4505 /* break into ordered pairs */
4506 parse_list_string(param_strings[0], SEP, '(', ')',
4507 &num_knots, &knot_coords);
4508 x = (double *)malloc(sizeof(double)*num_knots);
4510 y = (double *)malloc(sizeof(double)*num_knots);
4512 for (i=0; i<num_knots; i++) {
4513 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
4514 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
4517 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
4519 free_parsed_list(num_knots, knot_coords);
4520 return create_spline_param_t(num_knots, x, y);
4523 void destroy_spline_param_t(void *params)
4525 spline_param_t *p = (spline_param_t *)params;
4528 gsl_spline_free(p->spline);
4530 gsl_interp_accel_free(p->acc);
4540 <<kramers integ k global declarations>>=
4541 extern int num_kramers_integ_params;
4542 extern char *kramers_integ_param_descriptions[];
4543 extern char kramers_integ_param_string[];
4546 <<kramers integ k globals>>=
4547 int num_kramers_integ_params = 4;
4548 char *kramers_integ_param_descriptions[] = {
4549 "diffusion D, m^2 / s",
4550 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
4551 "starting integration bound x_min, m",
4552 "ending integration bound x_max, m"};
4553 //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";
4554 //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";
4555 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";
4556 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
4557 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
4558 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
4561 <<kramers integ k model getopt>>=
4562 {"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}
4564 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4565 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4567 \subsection{Unfolding model implementation}
4571 <<k model includes>>
4572 #include "k_model.h"
4573 <<k model internal definitions>>
4575 <<k model functions>>
4578 <<k model includes>>=
4579 #include <assert.h> /* assert() */
4580 #include <stdlib.h> /* NULL, malloc() */
4581 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4582 #include <stdio.h> /* fprintf(), stdout */
4583 #include <string.h> /* strlen(), strcpy() */
4584 #include <gsl/gsl_integration.h>
4585 #include <gsl/gsl_spline.h>
4590 <<k model internal definitions>>=
4591 <<const k structure>>
4592 <<bell k structure>>
4593 <<kramers k structure>>
4594 <<kramers integ k structure>>
4595 <<spline param structure>>
4598 <<k model globals>>=
4602 <<kramers k globals>>
4603 <<kramers integ k globals>>
4606 <<k model functions>>=
4607 <<null k functions>>
4608 <<const k functions>>
4609 <<bell k functions>>
4610 <<kramers k functions>>
4611 <<kramers integ k functions>>
4614 \subsection{Unfolding model unit tests}
4616 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
4617 <<check-k-model.c>>=
4618 <<k model check includes>>
4619 <<check relative error>>
4621 <<k model test suite>>
4622 <<main check program>>
4625 <<k model check includes>>=
4626 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
4627 #include <stdio.h> /* sprintf() */
4628 #include <assert.h> /* assert() */
4629 #include <math.h> /* exp() */
4632 #include "k_model.h"
4635 <<k model test suite>>=
4638 <<k model suite definition>>
4641 <<k model suite definition>>=
4642 Suite *test_suite (void)
4644 Suite *s = suite_create ("k model");
4645 <<const k test case defs>>
4646 <<bell k test case defs>>
4648 <<const k test case adds>>
4649 <<bell k test case adds>>
4654 \subsubsection{Constant}
4656 <<const k test case defs>>=
4657 TCase *tc_const_k = tcase_create("Constant k");
4660 <<const k test case adds>>=
4661 tcase_add_test(tc_const_k, test_const_k_create_destroy);
4662 tcase_add_test(tc_const_k, test_const_k_over_range);
4663 suite_add_tcase(s, tc_const_k);
4668 START_TEST(test_const_k_create_destroy)
4670 char *knot[] = {"1","2","3","4","5","6"};
4671 char *params[] = {knot[0]};
4674 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4675 params[0] = knot[i];
4676 p = string_create_const_k_param_t(params);
4677 destroy_const_k_param_t(p);
4682 START_TEST(test_const_k_over_range)
4684 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
4685 char *knot[] = {"1","2","3","4","5","6"};
4686 char *params[] = {knot[0]};
4689 char param_string[80];
4692 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4693 params[0] = knot[i];
4694 p = string_create_const_k_param_t(params);
4695 for ( F=Fm; F<FM; F+=dF ) {
4696 fail_unless(const_k(F, &env, p)==atof(knot[i]),
4697 "const_k(%g, %g, &{%s}) = %g != %s",
4698 F, T, knot[i], const_k(F, &env, p), knot[i]);
4700 destroy_const_k_param_t(p);
4706 \subsubsection{Bell}
4708 <<bell k test case defs>>=
4709 TCase *tc_bell = tcase_create("Bell's k");
4712 <<bell k test case adds>>=
4713 tcase_add_test(tc_bell, test_bell_k_create_destroy);
4714 tcase_add_test(tc_bell, test_bell_k_at_zero);
4715 tcase_add_test(tc_bell, test_bell_k_at_one);
4716 suite_add_tcase(s, tc_bell);
4720 START_TEST(test_bell_k_create_destroy)
4722 char *knot[] = {"1","2","3","4","5","6"};
4724 char *params[] = {knot[0], dx};
4727 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4728 params[0] = knot[i];
4729 p = string_create_bell_param_t(params);
4730 destroy_bell_param_t(p);
4735 START_TEST(test_bell_k_at_zero)
4737 double F=0.0, T=300.0;
4738 char *knot[] = {"1","2","3","4","5","6"};
4740 char *params[] = {knot[0], dx};
4743 char param_string[80];
4746 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4747 params[0] = knot[i];
4748 p = string_create_bell_param_t(params);
4749 fail_unless(bell_k(F, &env, p)==atof(knot[i]),
4750 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
4751 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
4752 destroy_bell_param_t(p);
4757 START_TEST(test_bell_k_at_one)
4760 char *knot[] = {"1","2","3","4","5","6"};
4762 char *params[] = {knot[0], dx};
4763 double F= k_B*T/atof(dx);
4766 char param_string[80];
4769 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4770 params[0] = knot[i];
4771 p = string_create_bell_param_t(params);
4772 CHECK_ERR(1e-6, atof(knot[i])*exp(1.0), bell_k(F, &env, p));
4773 //fail_unless(bell_k(F, &env, p)==atof(knot[i])*exp(1.0),
4774 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
4775 // F, T, knot[i], dx, bell_k(F, &env, p), atof(knot[i])*exp(1.0));
4776 destroy_bell_param_t(p);
4782 <<kramers k tests>>=
4785 <<kramers k test case def>>=
4788 <<kramers k test case add>>=
4791 <<k model function tests>>=
4792 <<check relative error>>
4794 START_TEST(test_something)
4796 double k=5, x=3, last_x=2.0, F;
4797 one_dim_fn_t *handlers[] = {&hooke};
4798 void *data[] = {&k};
4799 double xi[] = {0.0};
4801 int new_active_groups = 1;
4802 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
4803 fail_unless(F = k*x, NULL);
4808 \subsection{Utilities}
4810 The unfolding models can get complicated, and are sometimes hard to explain to others.
4811 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
4812 the overhead of having to understand a full simulation.
4813 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
4814 or special mode, where the operation depends on the specific model selected.
4816 <<k-model-utils.c>>=
4818 <<k model utility includes>>
4819 <<k model utility definitions>>
4820 <<k model utility getopt functions>>
4821 <<k model utility multi model E>>
4822 <<k model utility main>>
4825 <<k model utility main>>=
4826 int main(int argc, char **argv)
4828 <<k model getopt array definition>>
4829 k_model_getopt_t *model = NULL;
4834 int num_param_args; /* for INIT_MODEL() */
4835 char **param_args; /* for INIT_MODEL() */
4837 double special_xmin;
4838 double special_xmax;
4839 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
4840 &Fmax, &special_xmin, &special_xmax, &flags);
4841 if (flags & VFLAG) {
4842 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4844 INIT_MODEL("utility", model, model->params, params);
4847 if (model->k == NULL) {
4848 printf("No k function for model %s\n", model->name);
4852 printf("Fmax = %g <= 0\n", Fmax);
4858 printf("#F (N)\tk (%% pop. per s)\n");
4859 for (i=0; i<=N; i++) {
4860 F = Fmax*i/(double)N;
4861 k = (*model->k)(F, &env, params);
4863 printf("%g\t%g\n", F, k);
4868 if (model->k == NULL || model->k == &bell_k) {
4869 printf("No special mode for model %s\n", model->name);
4872 if (special_xmax <= special_xmin) {
4873 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
4877 double Xrng=(special_xmax-special_xmin), x, E;
4879 printf("#x (m)\tE (J)\n");
4880 for (i=0; i<=N; i++) {
4881 x = special_xmin + Xrng*i/(double)N;
4882 E = multi_model_E(model->k, params, &env, x);
4883 printf("%g\t%g\n", x, E);
4890 if (model->destructor)
4891 (*model->destructor)(params);
4896 <<k model utility multi model E>>=
4897 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
4899 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
4901 if (k == kramers_integ_k)
4902 E = (*p->E)(x, p->E_params);
4904 E = kramers_E(0, env, model_params, x);
4910 <<k model utility includes>>=
4913 #include <string.h> /* strlen() */
4914 #include <assert.h> /* assert() */
4917 #include "string_eval.h"
4918 #include "k_model.h"
4921 <<k model utility definitions>>=
4922 <<version definition>>
4923 #define VFLAG 1 /* verbose */
4924 enum mode_t {M_K_OF_F, M_SPECIAL};
4925 <<string comparison definition and globals>>
4926 <<k model getopt definitions>>
4927 <<initialize model definition>>
4928 <<kramers k structure>>
4929 <<kramers integ k structure>>
4933 <<k model utility getopt functions>>=
4935 <<k model utility help>>
4936 <<k model utility get options>>
4939 <<k model utility help>>=
4940 void help(char *prog_name,
4942 int n_k_models, k_model_getopt_t *k_models,
4946 printf("usage: %s [options]\n", prog_name);
4947 printf("Version %s\n\n", VERSION);
4948 printf("A utility for understanding the available unfolding models\n\n");
4949 printf("Environmental options:\n");
4950 printf("-T T\tTemperature T (currently %g K)\n", env->T);
4951 printf("-C T\tYou can also set the temperature T in Celsius\n");
4952 printf("Model options:\n");
4953 printf("-k model\tTransition rate model (currently %s)\n",
4954 k_models[k_model].name);
4955 printf("-K args\tTransition rate model argument string (currently %s)\n",
4956 k_models[k_model].params);
4957 printf("There are two output modes. In standard mode, k(F) is printed\n");
4958 printf("For example:\n");
4959 printf(" #Force (N)\tk (% pop. per s)\n");
4960 printf(" 123.456\t7.89\n");
4962 printf("In special mode, the output depends on the model.\n");
4963 printf("For models defining an energy function E(x), that function is printed\n");
4964 printf(" #Position (m)\tE (J)\n");
4965 printf(" 0.0012\t0.0034\n");
4967 printf("-m\tChange output to standard mode\n");
4968 printf("-M\tChange output to special mode\n");
4969 printf("-F\tSet the maximum F value for the standard mode k(F) output\n");
4970 printf("-x\tSet the minimum x value for the special mode E(x) output\n");
4971 printf("-X\tSet the maximum x value for the special mode E(x) output\n");
4972 printf("-V\tChange output to verbose mode\n");
4973 printf("-h\tPrint this help and exit\n");
4975 printf("Unfolding rate models:\n");
4976 for (i=0; i<n_k_models; i++) {
4977 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
4978 for (j=0; j < k_models[i].num_params ; j++)
4979 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
4980 printf("\t\tdefaults: %s\n", k_models[i].params);
4985 <<k model utility get options>>=
4986 void get_options(int argc, char **argv, environment_t *env,
4987 int n_k_models, k_model_getopt_t *k_models,
4988 enum mode_t *mode, k_model_getopt_t **model,
4989 double *Fmax, double *special_xmin, double *special_xmax,
4990 unsigned int *flags)
4992 char *prog_name = NULL;
4993 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
4995 extern char *optarg;
4996 extern int optind, optopt, opterr;
5000 /* setup defaults */
5002 prog_name = argv[0];
5003 env->T = 300.0; /* K */
5008 *special_xmax = 1e-8;
5009 *special_xmin = 0.0;
5011 while ((c=getopt(argc, argv, options)) != -1) {
5013 case 'T': env->T = atof(optarg); break;
5014 case 'C': env->T = atof(optarg)+273.15; break;
5016 k_model = index_k_model(n_k_models, k_models, optarg);
5017 *model = k_models+k_model;
5020 k_models[k_model].params = optarg;
5022 case 'm': *mode = M_K_OF_F; break;
5023 case 'M': *mode = M_SPECIAL; break;
5024 case 'F': *Fmax = atof(optarg); break;
5025 case 'x': *special_xmin = atof(optarg); break;
5026 case 'X': *special_xmax = atof(optarg); break;
5027 case 'V': *flags |= VFLAG; break;
5029 fprintf(stderr, "unrecognized option '%c'\n", optopt);
5030 /* fall through to default case */
5032 help(prog_name, env, n_k_models, k_models, k_model);
5041 \section{\LaTeX\ documentation}
5043 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
5044 The comment blocks are extracted (with nicely formatted code blocks), using
5045 <<latex makefile lines>>=
5046 $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw $(BUILD_DIR)
5047 noweave -latex -index -delay $< > $@
5048 $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib $(BUILD_DIR)
5052 <<latex makefile lines>>=
5053 $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
5055 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5056 $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
5057 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5058 $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
5059 mv $(BUILD_DIR)/sawsim.pdf $@
5061 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
5062 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
5063 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
5065 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
5066 <<latex makefile lines>>=
5068 rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \
5069 $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \
5070 $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.tex \
5071 $(BUILD_DIR)/sawsim.bib
5072 clean_latex : semi-clean_latex
5073 rm -f $(DOC_DIR)/sawsim.pdf
5078 [[make]] is a common utility on *nix systems for managing dependencies while building software.
5079 In this case, we have several \emph{targets} that we'd like to build:
5080 the main simulation program \prog;
5081 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
5082 the documentation [[sawsim.pdf]];
5083 and the [[Makefile]] itself.
5084 Besides the generated files, there is the utility target
5085 [[clean]] for removing all generated files except [[Makefile]].
5087 % [[dist]] for generating a distributable tar file.
5089 Extract the makefile with
5090 `[[notangle -Rmakefile src/sawsim.nw | sed 's/ /\t/' > Makefile]]'.
5091 The sed is needed because notangle replaces tabs with eight spaces,
5092 but [[make]] needs tabs.
5094 # Decide what directories to put things in
5099 # Note: Cannot use, for example, `./build', because make eats the `./'
5100 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5102 # Modules (X.c, X.h) defined in the noweb file
5105 # Modules defined outside the noweb file
5106 FREE_SAWSIM_MODS = interp tavl
5108 <<list module makefile lines>>
5109 <<tension balance module makefile lines>>
5110 <<tension model module makefile lines>>
5111 <<k model module makefile lines>>
5112 <<parse module makefile lines>>
5113 <<string eval module makefile lines>>
5114 <<accel k module makefile lines>>
5116 SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS)
5118 # Everything needed for sawsim under one roof. sawsim.c must come first
5119 SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \
5120 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h)
5121 # Libraries needed to compile sawsim
5122 LIBS = gsl gslcblas m
5123 CHECK_LIBS = gsl gslcblas m check
5124 # The non-check binaries generated
5125 BINS = sawsim tension_model_utils k_model_utils sawsim_profile
5128 # Define the major targets
5129 all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
5131 view : $(DOC_DIR)/sawsim.pdf
5133 profile : $(BIN_DIR)/sawsim_profile $(BIN_DIR)
5134 $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \
5135 -mnull -Mwlc -A0.39e-9,28e-9 -F8
5136 gprof $(BIN_DIR)/sawsim_profile gmon.out > $@
5137 check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim
5138 $(SHELL) -e -c 'for B in $^; do ./$$B; done'
5139 clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \
5140 clean_tension_model_utils \
5141 clean_k_model_utils clean_latex clean_check_sawsim
5142 rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \
5143 $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \
5144 $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
5145 $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
5146 $(BUILD_DIR)/global.h ./gmon.out
5147 $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
5149 # Various builds of sawsim
5150 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) $(BIN_DIR)
5151 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5152 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) $(BIN_DIR)
5153 gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5154 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) $(BIN_DIR)
5155 gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
5157 # Create the directories
5158 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
5161 # Copy over the source external to sawsim
5162 # Note: Cannot use, for example, `./build', because make eats the `./'
5163 # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...)
5165 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
5169 $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\
5173 ## Basic source generated with noweb
5174 # The central files sawsim.c and global.h...
5175 $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw $(BUILD_DIR)
5177 $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw $(BUILD_DIR)
5178 notangle -Rglobal.h $< > $@
5179 # ... and the modules
5180 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw $(BUILD_DIR)
5181 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5182 $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw $(BUILD_DIR)
5183 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5184 $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw $(BUILD_DIR)
5185 notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@
5186 # Note: I use `_' as a space in my file names, but noweb doesn't like
5187 # that so we use `-' to name the noweb roots and substitute here.
5189 ## Building the unit-test programs
5191 $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
5192 notangle -Rchecks $< > $@
5193 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
5194 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
5195 $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) $(BIN_DIR)
5196 gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
5197 clean_check_sawsim :
5198 rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
5199 # ... and the modules
5201 $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
5202 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \
5203 $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \
5204 $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
5205 $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
5206 $$(BUILD_DIR)/global.h $(BIN_DIR)
5207 gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
5208 $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
5210 # todo: clean up the required modules to
5211 $(CHECK_BINS:%=clean_%) :
5212 rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
5214 # Cleaning up the modules
5216 $(SAWSIM_MODS:%=clean_%) :
5217 rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h)
5219 <<tension model utils makefile lines>>
5220 <<k model utils makefile lines>>
5221 <<latex makefile lines>>
5223 Makefile : $(SRC_DIR)/sawsim.nw
5224 notangle -Rmakefile $< | sed 's/ /\t/' > $@
5226 This is fairly self-explanatory. We compile a dynamically linked
5227 version ([[sawsim]]) and a statically linked version
5228 ([[sawsim_static]]). The static version is about 9 times as big, but
5229 it runs without needing dynamic GSL libraries to link to, so it's a
5230 better format for distributing to a cluster for parallel evaluation.
5234 \subsection{Hookean springs in series}
5235 \label{app.math_hooke}
5238 The effective spring constant for $n$ springs $k_i$ in series is given by
5240 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5246 F &= k_i x_i &&\forall i \le n \\
5247 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
5248 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
5249 F &= k_1 x_1 = k_\text{eff} x \\
5250 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
5251 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5256 \addcontentsline{toc}{section}{References}
5257 \bibliography{sawsim}