1 %%%%%%%%%%%%%%%%%%%%%%%%% Start LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%
2 % we have our own preamble,
3 % use `noweave -delay` to not print noweb's automatic one
4 % -*- mode: Noweb; noweb-code-mode: c-mode -*-
5 \documentclass[letterpaper, 11pt]{article}
8 \noweboptions{smallcode}
10 \usepackage{url} % needed for natbib for underscores in urls?
11 \usepackage[breaklinks,colorlinks]{hyperref} % for internal links, \phantomsection
12 % breaklinks breaks long links
13 % colorlinks colors link text instead of boxing it
14 \usepackage{amsmath} % for \text{}, \xrightarrow[sub]{super}
15 \usepackage{amsthm} % for theorem and proof environments
16 \newtheorem{theorem}{Theorem}
18 \usepackage{doc} % BibTeX
19 \usepackage[super,sort&compress]{natbib} % fancy citation extensions
20 % super selects citations in superscript mode
21 % sort&compress automatically sorts and compresses compound citations (\cite{a,b,...})
23 %\bibliographystyle{ieeetr} % pick the bibliography style, short and sweet
24 %\bibliographystyle{plain} % pick the bibliography style, includes dates
25 \bibliographystyle{plainnat}
26 \defcitealias{sw:noweb}{{\tt noweb}}
27 \defcitealias{galassi05}{GSL}
28 \defcitealias{sw:check}{{\tt check}}
29 \defcitealias{sw:python}{Python}
34 \textheight 9.5in % leave a bit of extra room for page numbers (at top of page)
37 \setlength{\parindent}{0pt}
38 \setlength{\parskip}{5pt}
40 % For some reason, the \maketitle wants to go on it's own page
41 % so we'll just hardcode the following in our first page.
42 %\title{Sawsim: a sawtooth protein unfolding simulator}
43 %\author{W.~Trevor King}
46 \newcommand{\prog}{[[sawsim]]}
47 \newcommand{\lang}{[[C]]}
48 \newcommand{\p}[3]{\left#1 #2 \right#3} % e.g. \p({complicated expr})
49 \newcommand{\Th}{\ensuremath{\sup{\text{th}}}} % ordinal th
50 \newcommand{\U}[1]{\ensuremath{\text{ #1}}} % units: $1\U{m/s}$
52 % single spaced lists, from Alvin Alexander
53 % http://www.devdaily.com/blog/post/latex/control-line-spacing-in-itemize-enumerate-tags/
54 \newenvironment{packed_item}{
56 \setlength{\itemsep}{1pt}
57 \setlength{\parskip}{0pt}
58 \setlength{\parsep}{0pt}
62 \nwfilename{sawsim.nw}
67 Sawsim: a sawtooth protein unfolding simulator \\
72 %%%%%%%%%%%%%%%%%%%%%%%%%%%% End LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%%
74 \section{Introduction}
76 The unfolding of multi-globular protein strings is a stochastic
77 process. In the \prog\ program, we use Monte Carlo methods to
78 simulate this unfolding through an explicit time-stepping approach.
79 We develop a program in \lang\ to simulate probability of unfolding
80 under a constant extension velocity of a chain of globular domains.
82 \subsection{Related work}
86 \subsection{About this document}
88 This document was written in \citetalias{sw:noweb} with code blocks in
89 \lang\ and comment blocks in \LaTeX.
91 \subsection{Change Log}
93 Version 0.0 used only the unfolded domain WLC to determine the tension
94 and had a fixed timestep $dt$.
96 Version 0.1 added adaptive timesteps, adjusting so that the unfolding
97 probability for a given domain was constant.
99 Version 0.2 added Kramers' model unfolding rate calculations, and a
100 switch to select Bell's or Kramers' model. This induced a major
101 rewrite, introducing the tension group abstraction, and a split to
102 multiple \lang\ source files. Also added a unit-test suites for
103 sanity checking, and switched to SI units throughout.
105 Version 0.3 added integrated Kramers' model unfolding rate
106 calculations. Also replaced some of my hand-coded routines with GNU
107 Scientific Library \citepalias{galassi05} calls.
109 Version 0.4 added Python evaluation for the saddle-point Kramers'
110 model unfolding rate calculations, so that the model functions could
111 be controlled from the command line. Also added interpolating lookup
112 tables to accelerate slow unfolding rate model evaluation, and command
113 line control of tension groups.
115 <<version definition>>=
116 #define VERSION "0.4"
122 sawsim - program for simulating single molecule mechanical unfolding.
123 Copyright (C) 2008, William Trevor King
125 This program is free software; you can redistribute it and/or
126 modify it under the terms of the GNU General Public License as
127 published by the Free Software Foundation; either version 3 of the
128 License, or (at your option) any later version.
130 This program is distributed in the hope that it will be useful, but
131 WITHOUT ANY WARRANTY; without even the implied warranty of
132 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
133 See the GNU General Public License for more details.
135 You should have received a copy of the GNU General Public License
136 along with this program; if not, write to the Free Software
137 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
140 The author may be contacted at <wking@drexel.edu> on the Internet, or
141 write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
142 Philadelphia PA 19104, USA.
154 The unfolding system is stored as a chain of \emph{domains}. Domains
155 can be folded, globular protein domains, unfolded protein linkers, AFM
156 cantilevers, or other stretchable system components.
158 Each domain contributes to the total tension, which depends on the
159 chain extension. The particular model for the tension as a function
160 of extension is handled generally with function pointers. So far the
161 following models have been implemented:
163 \item Null (Appendix \ref{sec.null_tension}),
164 \item Constant (Appendix \ref{sec.const_tension}),
165 \item Hookean spring (Appendix \ref{sec.hooke}),
166 \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
167 \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
170 The tension model and parameters used for a particular domain depend
171 on whether the domain is folded or unfolded. The transition rate from
172 the folded state to the unfolded state is also handled generally with
173 function pointers, with implementations for
175 \item Null (Appendix \ref{sec.null_k}),
176 \item Bell's model (Appendix \ref{sec.bell}),
177 \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
178 \item Kramers' saddle point model (Appendix \ref{sec.kramers}).
181 The unfolding of the chain domains is modeled in two stages. First
182 the tension of the chain is calculated. Then the tension (assumed to
183 be constant over the length of the chain), is applied to each folded
184 domain, and used to calculate the probability of that subdomain
185 unfolding in the next timestep $dt$. Then the time is stepped
186 forward, and the process is repeated.
189 int main(int argc, char **argv)
191 <<tension model getopt array definition>>
192 <<k model getopt array definition>>
193 list_t *domain_list=NULL; /* variables initialized in get_options() */
194 environment_t env={0};
196 int num_folded=0, unfolding;
197 double x=0, dt, F; /* variables used in the simulation loop */
198 get_options(argc, argv, &P, &dt_max, &v, &env, NUM_TENSION_MODELS,
199 tension_models, NUM_K_MODELS, k_models, &domain_list, &num_folded);
201 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Position (m)\tForce (N)\n");
202 if (flags & FLAG_OUTPUT_UNFOLDING_FORCES) printf("#Unfolding Force (N)\n");
203 while (num_folded > 0) {
204 F = find_tension(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, unfolding);
205 dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, v);
206 unfolding=random_unfoldings(domain_list, F, dt, &env, &num_folded);
207 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\n", x, F);
210 destroy_domain_list(domain_list);
214 @ The meat of the simulation is bundled into the three functions
215 [[find_tension]], [[determine_dt]], and [[random_unfoldings]].
216 [[find_tension]] is discussed in Section \ref{sec.find_tension},
217 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
218 [[random_unfoldings]] in Section \ref{sec.unfolding_rate}.
220 Environmental parameters important in determining reaction rates and
221 tensions (e.g. temperature, pH) are stored in a single structure to
222 facilitate extension to more complicated models in the future.
223 <<environment definition>>=
224 typedef struct environment_struct {
230 <<environment definition>>
231 <<one dimensional function definition>>
232 <<create/destroy definitions>>
236 \section{Simulation functions}
238 <<simulation functions>>=
242 <<does domain unfold>>
243 <<randomly unfold domains>>
247 \label{sec.find_tension}
249 Because the stretched system may be made up of several parts (folded
250 domains, unfolded domains, spring-like cantilever, \ldots), we will
251 assign the domains to models and groups. The models are used to
252 determine a domain's tension (Hookean spring, WLC, \ldots). Groups
253 will represent collections of domains which the model should treat as
254 a single entity. A domain may belong to different groups or models
255 depending on its state. For example, a domain might be modeled as a
256 freely-jointed chain when it iss folded and as a worm-like chain class
257 when it is unfolded. The domains are assumed to be commutative, so
258 ordering is ignored. The interactions between the groups are assumed
259 to be linear, but the interactions between domains of the same group
260 need not be. This allows for non-linear group models such as th
261 worm-like or freely-jointed chains. Each model has a tension handler
262 function, which gives the tension $F$ needed to stretch a given group
263 of domains a distance $x$.
265 To understand the purpose of groups, consider a chain of proteins
266 where two folded proteins $A$ and $B$ are modeled as WLCs with
267 persistence length $p_f$ and two unfolded proteins $C$ and $D$ are
268 modeled as WLCs with persistence length $p_u$. The proteins are
269 attached to a cantilever $E$ qof spring constant $κ$. The string
270 would get split into three lists:
272 \begin{tabular}{llll}
273 Model & Group & List & Tension \\
274 WLC & 0 & $[A,B]$ & $F_\text{WLC}(x, p_f, L_A+L_B}$ \\
275 WLC & 1 & $[C,D]$ & $F_\text{WLC}(x, p_u, L_C+L_D}$ \\
276 Hooke & 0 & $[E]$ & $F_\text{Hooke}(x, κ}$ \\
279 Note that group numbers only matter \emph{within} model classes, since
280 grouping domains with seperate models doesn't make sense.
282 <<find tension definitions>>=
283 #define NUM_TENSION_MODELS 5
287 <<tension handler array global declaration>>=
288 extern one_dim_fn_t *tension_handlers[];
291 <<tension handler array global>>=
292 one_dim_fn_t *tension_handlers[] = {
294 &const_tension_handler,
302 <<tension model global declarations>>=
303 <<tension handler array global declaration>>
306 <<tension handler types>>=
307 typedef struct tension_handler_data_struct {
308 /* int num_domains; */
309 /* domain_t *domains; */
313 } tension_handler_data_t;
317 After sorting the chain into separate groups $G_i$, with tension
318 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
319 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
320 \forall i,j$. Note that there may be several groups within each model
321 class. [[find_tension]] is basically a wrapper to feed proper input
322 into the [[tension_balance]] function. [[unfolding]] should be set to
323 0 if there were no unfoldings since the previous call to
326 double find_tension(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list, environment_t *env, double x, int unfolding)
328 static int num_active_groups=0;
329 static one_dim_fn_t **per_group_handlers = NULL;
330 static void **per_group_params = NULL;
331 static double last_x;
332 tension_handler_data_t *tdata;
337 fprintf(stderr, "finding tension at distance %g%s\n", x, unfolding ? " after an unfolding" : "");
338 list_print(stderr, domain_list, "domain list");
341 if (unfolding != 0 || num_active_groups == 0) { /* setup statics */
342 /* free old statics */
343 if (per_group_handlers != NULL)
344 free(per_group_handlers);
345 if (per_group_params != NULL) {
346 for (i=0; i < num_active_groups; i++) {
347 assert(per_group_params[i] != NULL);
348 free(per_group_params[i]);
350 free(per_group_params);
353 /* get new statics */
354 get_active_groups(domain_list, num_tension_handlers, tension_handlers, &num_active_groups, &per_group_handlers, &per_group_params);
356 /* fill in the group handlers and params */
357 for (i=0; i<num_active_groups; i++) {
358 tdata = (tension_handler_data_t *) per_group_params[i];
360 fprintf(stderr, "active group %d of %d", i, num_active_groups);
361 list_print(stderr, tdata->group, " parameters");
364 /* tdata->persist continues unchanged */
367 } /* else roll with the current statics */
369 F = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_params, last_x, x);
373 @ For the moment, we will restrict our group boundaries to a single
374 dimension, so $\sum_i x_i = x$, our total extension (it is this
375 restriction that causes the groups to interact linearly). We'll also
376 restrict our extensions to all be positive. With these restrictions,
377 the problem of balancing the tensions reduces to solving a system of
378 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
379 the number of active groups. In general this can be fairly
380 complicated, but without much loss of practicality we can restrict
381 ourselves to strictly monotonically increasing, non-negative tension
382 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
383 simpler. For example, it guarantees the existence of a unique, real
384 solution for finite forces. See Appendix \ref{app.tension_balance}
385 for the tension balancing implementation.
387 Each group must define a parameter structure if the tension model is
389 a function to create the parameter structure,
390 a function to destroy the parameter structure, and
391 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as it's data argument.
392 The tension-specific parameter structures are contained in the domain
393 group field. See the Section \ref{app.model_selection} for a
394 discussion on the command line framework. See the worm-like chain
395 implementation for an example (Section \ref{sec.wlc}).
397 The major limitation of our approach is that all of our tension models
398 are Markovian (memory-less), which is not really a problem for the
399 chains (see \citet{evans99} for WLC thermalization rate calculations).
401 \subsection{Unfolding rate}
402 \label{sec.unfolding_rate}
404 Each folded domain is modeled as unimolecular, first order reaction
406 \text{Folded} \xrightarrow{k} % in tex: X atop Y
409 With the rate constant $k$ defined by
411 \frac{d[\text{Folded}]}{dt} = -k [\text{Folded}].
413 So the probability of a given protein unfolding in a short time $dt$
419 <<does domain unfold>>=
420 int domain_unfolds(double F, double dt, environment_t *env, domain_t *domain)
421 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to a folded domain */
423 k = accel_k(domain->k, F, env, domain->k_params);
424 //(*domain->k)(F, env, domain->k_params);
425 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
426 return happens(k*dt); /* dice roll for prob. k*dt event */
428 @ [[happens]] is a random decision making function defined in Appendix \ref{app.utils}.
430 Only one domain can unfold in each timestep, because the timescale of
431 a domain unfolding $dt_u$ is assumed to be less than the simulation
432 timestep $dt$, so a domain will completely unfold in a single
433 timestep. We adapt our timesteps to keep the probability of a single
434 domain unfolding low, so the probability of two domains unfolding in
435 the same timestep is negligible. This approach breaks down as the
436 adaptive timestep scheme approaches $dt \sim dt_u$, but $dt_u \sim
437 1\U{$\mu$s}$ for Ig27-like proteins \citep{klimov00}, so this
438 shouldn't be a problem. To reassure yourself, you can ask the
439 simulator to print the smallest timestep that was used with TODO.
440 <<randomly unfold domains>>=
441 int random_unfoldings(list_t *domain_list, double tension,
442 double dt, environment_t *env,
443 int *num_folded_domains)
444 { /* returns 1 if there was an unfolding and 0 otherwise */
445 while (domain_list != NULL) {
446 if (D(domain_list)->state == DS_FOLDED
447 && domain_unfolds(tension, dt, env, D(domain_list))) {
448 if (flags & FLAG_OUTPUT_UNFOLDING_FORCES)
449 fprintf(stdout, "%g\n", tension);
450 D(domain_list)->state = DS_UNFOLDED;
451 (*num_folded_domains)--;
452 return 1; /* our one domain has unfolded, stop looking for others */
454 domain_list = domain_list->next;
460 \subsection{Adaptive timesteps}
461 \label{sec.adaptive_dt}
463 We'd like to pick $dt$ so the probability of unfolding in the next
464 timestep is small. If we don't adapt our timestep, we risk breaking
465 our approximation that $F(x' \in [x, x+v \cdot dt]) = F(x)$ or that
466 only one domain unfolds in a given timestep. Because $F(x)$ is
467 monotonically increasing, excessively large timesteps will lead to
468 erroneously large unfolding forces. The simulation would have been
469 accurate for sufficiently small fixed timesteps, but adaptive
470 timesteps will allow us to move through low tension regions in fewer
471 steps, leading to a more efficient simulation.
473 The actual adaptive timestep implementation is not particularly
474 interesting, since we are only required to reduce $dt$ to somewhere
475 below a set threshold, so I've removed it to Appendix
476 \ref{app.adaptive_dt}.
482 The models provide the physics of the simulation, but the simulation
483 allows interchangeable models, and we are currently focusing on the
484 simulation itself, so we remove the actual model implementation to the
487 The tension models are defined in Appendix \ref{sec.tension_models}
488 and unfolding models are defined in Appendix \ref{sec.k_models}.
491 #define k_B 1.3806503e-23 /* J/K */
495 \section{Command line interface}
497 <<get option functions>>=
499 <<index tension model>>
505 \subsection{Model selection}
506 \label{app.model_selection}
508 The main difficulty with the command line interface in \prog\ is
509 developing an intuitive method for accessing the possibly large number
510 of available models. We'll treat this generally by defining an array
511 of available models, containing their important parameters.
513 <<get options definitions>>=
514 <<model getopt definitions>>
517 <<create/destroy definitions>>=
518 typedef void *create_data_func_t(char **param_strings);
519 typedef void destroy_data_func_t(void *param_struct);
522 <<model getopt definitions>>=
523 <<tension model getopt definitions>>
524 <<k model getopt definitions>>
528 \subsubsection{Tension}
530 To access the various tension models from the command line, we define
531 a structure that contains all the neccessary information about the
533 <<tension model getopt definitions>>=
534 typedef struct tension_model_getopt_struct {
535 one_dim_fn_t *handler; /* fn implementing the model on a group */
536 char *name; /* name string identifying the model */
537 char *description; /* a brief description of the model */
538 int num_params; /* number of extra params for the model */
539 char **param_descriptions; /* descriptions of extra params */
540 char *params; /* default values for the extra params */
541 create_data_func_t *creator; /* param-string -> model-data-struct fn */
542 destroy_data_func_t *destructor; /* fn to free the model data structure */
543 } tension_model_getopt_t;
546 <<tension model getopt array definition>>=
547 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
548 <<null tension model getopt>>,
549 <<constant tension model getopt>>,
550 <<hooke tension model getopt>>,
551 <<worm-like chain tension model getopt>>,
552 <<freely-jointed chain tension model getopt>>
556 \subsubsection{Unfolding rate}
558 <<k model getopt definitions>>=
559 #define NUM_K_MODELS 5
561 typedef struct k_model_getopt_struct {
566 char **param_descriptions;
568 create_data_func_t *creator;
569 destroy_data_func_t *destructor;
573 <<k model getopt array definition>>=
574 k_model_getopt_t k_models[NUM_K_MODELS] = {
575 <<null k model getopt>>,
576 <<const k model getopt>>,
577 <<bell k model getopt>>,
578 <<kramers k model getopt>>,
579 <<kramers integ k model getopt>>
586 void help(char *prog_name, double P, double dt_max, double v,
588 int n_tension_models, tension_model_getopt_t *tension_models,
589 int folded_tension_model, int unfolded_tension_model,
590 int n_k_models, k_model_getopt_t *k_models,
594 printf("usage: %s [options]\n", prog_name);
595 printf("Version %s\n\n", VERSION);
596 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
597 printf("Simulation options:\n");
598 printf("-P P\tTarget probability for dt (currently %g)\n", P);
599 printf("-t dt\tMaximum allowed timestep dt (currently %g)\n", dt_max);
600 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
601 printf("Environmental options:\n");
602 printf("-T T\tTemperature T (currently %g K)\n", env->T);
603 printf("-C T\tYou can also set the temperature T in Celsius\n");
604 printf("Model options:\n");
605 printf("The domains exist in either the folded or unfolded state\n");
606 printf("The following options change the default behavior in each state.\n");
607 printf("-m model[,group]\tFolded tension model (currently %s)\n",
608 tension_models[folded_tension_model].name);
609 printf("-a args\tFolded tension model argument string (currently %s)\n",
610 tension_models[folded_tension_model].params);
611 printf("-M model[,group]\tUnfolded tension model (currently %s)\n",
612 tension_models[unfolded_tension_model].name);
613 printf("-A args\tUnfolded tension model argument string (currently %s)\n",
614 tension_models[unfolded_tension_model].params);
615 printf("The following options change the unfolding rate.\n");
616 printf("-k model\tTransition rate model (currently %s)\n",
617 k_models[k_model].name);
618 printf("-K args\tTransition rate model argument string (currently %s)\n",
619 k_models[k_model].params);
620 printf("Domain creation options:\n");
621 printf("Once you've set up the appropriate default models, you need to add the domains\n");
622 printf("-F n\tAdd n folded domains with the current models\n");
623 printf("-U n\tAdd n unfolded domains with the current models\n");
624 printf("Output mode options:\n");
625 printf("There are two output modes. In standard mode, only the unfolding\n");
626 printf("events are printed. For example:\n");
627 printf(" #Unfolding Force (N)\n");
628 printf(" 123.456e-12\n");
630 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
631 printf(" #Position (m)\tForce (N)\n");
632 printf(" 0.001\t0.0023\n");
634 printf("-V\tChange output to verbose mode\n");
635 printf("-h\tPrint this help and exit\n");
637 printf("Tension models:\n");
638 for (i=0; i<n_tension_models; i++) {
639 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
640 for (j=0; j < tension_models[i].num_params ; j++)
641 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
642 printf("\t\tdefaults: %s\n", tension_models[i].params);
644 printf("Unfolding rate models:\n");
645 for (i=0; i<n_k_models; i++) {
646 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
647 for (j=0; j < k_models[i].num_params ; j++)
648 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
649 printf("\t\tdefaults: %s\n", k_models[i].params);
651 printf("Examples:\n");
652 printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded\n");
653 printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8\n", prog_name);
654 printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded\n");
655 printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K5e-5,.225e-9 -mwlc -a1e-8,4e-10 -Mwlc,1 -A0.4e-9,28.1e-9 -F8\n", prog_name);
659 \subsection{Get options}
662 void get_options(int argc, char **argv,
663 double *pP, double *pDt_max, double *pV,
665 int n_tension_models, tension_model_getopt_t *tension_models,
666 int n_k_models, k_model_getopt_t *k_models,
667 list_t **pDomain_list, int *pNum_folded)
669 char *prog_name = NULL;
670 char c, options[] = "P:t:v:T:C:m:a:M:A:k:K:F:U:Vh";
671 int ftension_model=0, utension_model=0, k_model=0;
672 char *ftension_params=NULL, *utension_params=NULL;
673 int i, n, ftension_group=0, utension_group=0;
677 extern int optind, optopt, opterr;
682 for (i=0; i<n_tension_models; i++) {
683 fprintf(stderr, "tension model %i %s:\t model's handler %p\n",
684 i, tension_models[i].name, tension_models[i].handler);
685 assert(tension_models[i].handler == tension_handlers[i]);
690 flags = FLAG_OUTPUT_UNFOLDING_FORCES;
692 *pP = 1e-3; /* % pop per s */
693 *pDt_max = 0.001; /* s */
694 *pV = 1e-6; /* m/s */
695 env->T = 300.0; /* K */
697 while ((c=getopt(argc, argv, options)) != -1) {
699 case 'P': *pP = atof(optarg); break;
700 case 't': *pDt_max = atof(optarg); break;
701 case 'v': *pV = atof(optarg); break;
702 case 'T': env->T = atof(optarg); break;
703 case 'C': env->T = atof(optarg)+273.15; break;
705 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
706 assert(num_strings == 1 || num_strings == 2);
707 if (num_strings == 1)
710 ftension_group = atoi(strings[1]);
711 ftension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
712 ftension_params = tension_models[ftension_model].params;
713 free_parsed_list(num_strings, strings);
716 ftension_params = optarg;
719 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
720 assert(num_strings == 1 || num_strings == 2);
721 if (num_strings == 1)
724 utension_group = atoi(strings[1]);
725 utension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
726 utension_params = tension_models[utension_model].params;
727 free_parsed_list(num_strings, strings);
730 utension_params = optarg;
733 k_model = index_k_model(n_k_models, k_models, optarg);
736 k_models[k_model].params = optarg;
741 for (i=0; i<n; i++) {
742 push(pDomain_list, generate_domain(DS_FOLDED,
743 tension_models+ftension_model,
746 tension_models+utension_model,
749 k_models+k_model), 1);
756 for (i=0; i<n; i++) {
757 push(pDomain_list, generate_domain(DS_UNFOLDED,
758 tension_models+ftension_model,
761 tension_models+utension_model,
764 k_models+k_model), 1);
768 flags = FLAG_OUTPUT_FULL_CURVE;
771 fprintf(stderr, "unrecognized option '%c'\n", optopt);
772 /* fall through to default case */
774 help(prog_name, *pP, *pDt_max, *pV, env, n_tension_models, tension_models, ftension_model, utension_model, n_k_models, k_models, k_model);
778 /* check the arguments */
779 assert(*pDomain_list != NULL); /* no domains! */
780 assert(*pP > 0.0); assert(*pP < 1.0);
781 assert(*pDt_max > 0.0);
783 assert(env->T > 0.0);
788 <<index tension model>>=
789 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
792 for (i=0; i<n_models; i++)
793 if (STRMATCH(models[i].name, name))
796 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
804 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
807 for (i=0; i<n_models; i++)
808 if (STRMATCH(models[i].name, name))
811 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
818 <<initialize model definition>>=
819 /* requires int num_param_args and char **param_args in the current scope
821 * INIT_MODEL("folded", folded_model, folded_params);
822 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
824 #define INIT_MODEL(role, model, param_string, param_pointer) \
826 parse_list_string(param_string, SEP, '{', '}', \
827 &num_param_args, ¶m_args); \
828 if (num_param_args != model->num_params) { \
830 "%s model %s expected %d params," \
831 role, model->name, model->num_params); \
833 "not the %d params in '%s'\n", \
834 num_param_args, param_string); \
835 assert(num_param_args == model->num_params); \
837 if (model->creator) \
838 param_pointer = (*model->creator)(param_args); \
840 param_pointer = NULL; \
841 free_parsed_list(num_param_args, param_args); \
846 void *generate_domain(enum domain_state_t state,
847 tension_model_getopt_t *folded_model,
849 char *folded_param_string,
850 tension_model_getopt_t *unfolded_model,
852 char *unfolded_param_string,
853 k_model_getopt_t *k_model)
855 void *folded_params, *unfolded_params, *k_params;
856 int num_param_args; /* for INIT_MODEL() */
857 char **param_args; /* for INIT_MODEL() */
860 fprintf(stderr, "generating %s ", state==DS_FOLDED ? "folded" : "unfolded");
861 fprintf(stderr, "domain (k: \"%s\", \"%s\", f: \"%s\" (%p) group %d \"%s\", u: \"%s\" (%p) group %d \"%s\")\n",
862 k_model->name, k_model->params,
863 folded_model->name, folded_model->handler, folded_group, folded_param_string,
864 unfolded_model->name, unfolded_model->handler, unfolded_group, unfolded_param_string);
866 INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
867 INIT_MODEL("unfolded", unfolded_model, unfolded_param_string, unfolded_params);
868 INIT_MODEL("k", k_model, k_model->params, k_params);
869 return create_domain(state,
870 k_model->k, k_params, k_model->destructor,
871 folded_model->handler, folded_group, folded_params, folded_model->destructor,
872 unfolded_model->handler, unfolded_group, unfolded_params, unfolded_model->destructor
879 \addcontentsline{toc}{section}{Appendicies}
881 \section{sawsim.c details}
885 The general layout of our simulation code is:
895 We include [[math.h]], so don't forget to link to the libm with `-lm'.
897 #include <assert.h> /* assert() */
898 #include <stdlib.h> /* malloc(), free(), rand() */
899 #include <stdio.h> /* fprintf(), stdout */
900 #include <string.h> /* strlen, strtok() */
901 #include <math.h> /* exp(), M_PI, sqrt() */
902 #include <sys/types.h> /* pid_t (returned by getpid()) */
903 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
906 #include "tension_balance.h"
908 #include "tension_model.h"
914 <<version definition>>
916 <<max/min definitions>>
917 <<string comparison definition and globals>>
918 <<initialize model definition>>
919 <<get options definitions>>
920 <<domain definitions>>
929 <<create/destroy domain>>
930 <<destroy domain list>>
932 <<simulation functions>>
933 <<boilerplate functions>>
936 <<boilerplate functions>>=
938 <<get option functions>>
944 srand(getpid()*time(NULL)); /* seed rand() */
948 <<flag definitions>>=
949 /* in octal b/c of prefixed '0' */
950 #define FLAG_OUTPUT_FULL_CURVE 01
951 #define FLAG_OUTPUT_UNFOLDING_FORCES 02
955 static unsigned long int flags = 0;
958 \subsection{Utilities}
961 <<max/min definitions>>=
962 #define MAX(a,b) ((a)>(b) ? (a) : (b))
963 #define MIN(a,b) ((a)<(b) ? (a) : (b))
966 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
967 <<string comparison definition and globals>>=
968 // Check if two strings match, return 1 if they do
969 static char *temp_string_A;
970 static char *temp_string_B;
971 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
972 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
973 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
974 /* +1 to also compare the '\0' */
977 We also define a macro for our [[check]] unit testing
978 <<check relative error>>=
979 #define CHECK_ERR(max_err, expected, received) \
981 fail_unless( (received-expected)/expected < max_err, \
982 "relative error %g >= %g in %s (Expected %g, received %g)", \
983 (received-expected)/expected, max_err, #received, \
984 expected, received); \
985 fail_unless(-(received-expected)/expected < max_err, \
986 "relative error %g >= %g in %s (Expected %g, received %g)", \
987 -(received-expected)/expected, max_err, #received, \
988 expected, received); \
993 int happens(double probability)
995 assert(probability >= 0.0); assert(probability <= 1.0);
996 return (double)rand()/RAND_MAX < probability; /* TODO: replace with GSL rand http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html x*/
1001 \subsection{Adaptive timesteps}
1002 \label{app.adaptive_dt}
1004 $F(x)$ increases with $x$, possibly exploding, as in the worm-like chain model,
1005 so basing the timestep on the the unfolding probability at the current tension
1006 is dangerous, and we need to search for a $dt$ for which
1007 $P(F(x+v*dt)) < P_\text{target}$.
1008 There are two cases to consider.
1009 In the most common, no domains have unfolded since the last step,
1010 and we expect the next step to be slightly shorter than the previous one.
1011 In the less common, domains did unfold in the last step,
1012 and we expect the next step to be considerably longer than the previous one.
1014 double search_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers,
1015 list_t *domain_list,
1016 environment_t *env, double x,
1017 double target_prob, double max_dt, double v)
1018 { /* Returns the timestep dt in seconds for the current folded domain.
1019 * Takes a list of tension handlers, the list of domains,
1020 * a pointer env to the environmental data, a starting separation x in m,
1021 * a target_prob between 0 and 1,
1022 * max_dt in s, stretching velocity v in m/s.
1024 double F, k, dtCur, dtU, dtUCur, dtL, dt;
1026 /* get upper bound using the current position */
1027 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x, 0); /* BUG. repeated calculation */
1028 //printf("Start with x = %g (F = %g)\n", x, F);
1029 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1030 //printf("x %g\tF %g\tk %g\n", x, F, k);
1031 dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */
1033 //printf("overshot max_dt\n");
1036 /* set a lower bound on dt too */
1039 /* The dt determined above may produce illegitimate forces or ks.
1040 * Reduce the upper bound until we have valid ks. */
1042 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1043 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1046 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1048 //printf("Try for dt = %g (F = %g)\n", dt, F);
1049 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1050 /* returned k may be -1.0 */
1051 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1052 while (k == -1.0) { /* reduce step until we hit a valid k */
1054 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1055 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1056 //printf("Try for dt = %g (F = %g)\n", dt, F);
1057 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1058 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1060 assert(dtU > 1e-14); /* timestep to valid k too small */
1061 dtUCur = target_prob / k; /* safe timestep back from x+dtU */
1063 return dt; /* dtU is safe. */
1065 /* dtU wasn't safe, lets see what would be. */
1066 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1067 dt = (dtU + dtL) / 2.0;
1068 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1069 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1070 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1071 dtCur = target_prob / k;
1072 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1073 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1075 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1076 dtU = dt; /* ... stepping out only dtCur would be safe */
1079 break; /* dtCur = dt */
1081 return MAX(dtUCur, dtL);
1085 To determine $dt$ for an array of potentially different folded domains,
1086 we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
1089 double determine_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list,
1090 environment_t *env, double x,
1091 double target_prob, double dt_max, double v)
1092 { /* Returns the timestep dt in seconds.
1093 * Takes the list of folded domains, target_prob between 0 and 1,
1094 * F in N, and T in K. */
1095 double dt=dt_max, new_dt;
1096 assert(target_prob > 0.0); assert(target_prob < 1.0);
1097 assert(dt_max > 0.0);
1099 /* .5 nm steps = v * dt */
1101 while (domain_list != NULL) {
1102 if (D(domain_list)->state == DS_FOLDED) {
1103 new_dt = search_dt(num_tension_handlers, tension_handlers, domain_list, env, x, target_prob, dt, v);
1104 dt = MIN(dt, new_dt);
1106 domain_list = domain_list->next;
1112 \subsection{Domain data}
1114 Currently domains exist in two states, folded and unfolded, and the
1115 only allowed transitions are folded $\rightarrow$ unfolded. Of
1116 course, it wouldn't be too complicated to extent this to a multi-state
1117 system, with an array containing the domains group for each possible
1118 state, and a matrix of transition-rate-calculating functions.
1119 However, at this point such generality seems unnecessary at this
1121 <<domain definitions>>=
1122 enum domain_state_t {DS_FOLDED,
1126 typedef struct domain_struct {
1127 enum domain_state_t state;
1128 one_dim_fn_t *folded_handler;
1130 one_dim_fn_t *unfolded_handler;
1132 k_func_t *k; /* function returning unfolding rate */
1133 void *folded_params; /* pointer to folded parameters */
1134 void *unfolded_params; /* pointer to unfolded parameters */
1135 void *k_params; /* pointer to k parameters */
1136 destroy_data_func_t *destroy_folded;
1137 destroy_data_func_t *destroy_unfolded;
1138 destroy_data_func_t *destroy_k;
1141 /* get the domain data for the current list node */
1142 #define D(list) ((domain_t *)(list)->d)
1143 /* get the tension params for the current list node */
1144 #define D_TP(list) (((domain_t *)(list)->d)->state == DS_FOLDED \
1145 ? ((domain_t *)(list)->d)->folded_params \
1146 : ((domain_t *)(list)->d)->unfolded_params)
1148 [[k]] is a pointer to the function determining the unfolding rate for a given tension.
1149 [[folded_params]] is a pointer to the parameters used by the function pointed to by [[k]].
1150 [[unfolded_params]] is a pointer to the parameters used by the group-appropriate handler function when determining the tension.
1151 The [[destroy_*]] pointers point to functions for freeing the memory [[*_params]].
1152 We store them with the domain data so that [[destroy_domain]] doesn't have to know which type of domain it's cleaning up after.
1154 [[create_]] and [[destroy_domain]] are simple wrappers around [[malloc]] and [[free]].
1155 <<create/destroy domain>>=
1156 domain_t *create_domain(enum domain_state_t state,
1159 destroy_data_func_t *destroy_k,
1160 one_dim_fn_t *folded_handler,
1162 void *folded_params,
1163 destroy_data_func_t *destroy_folded,
1164 one_dim_fn_t *unfolded_handler,
1166 void *unfolded_params,
1167 destroy_data_func_t *destroy_unfolded)
1169 domain_t *ret = (domain_t *)malloc(sizeof(domain_t));
1170 assert(ret != NULL);
1171 if (state == DS_FOLDED) {
1172 assert(k != NULL); /* the pointer points somewhere valid */
1173 assert(*k != NULL); /* and there is something useful there */
1175 assert(state == DS_UNFOLDED);
1177 ret->folded_handler = folded_handler;
1178 ret->folded_group = folded_group;
1179 ret->unfolded_handler = unfolded_handler;
1180 ret->unfolded_group = unfolded_group;
1182 ret->folded_params = folded_params;
1183 ret->unfolded_params = unfolded_params;
1184 ret->k_params = k_params;
1185 ret->destroy_folded = destroy_folded;
1186 ret->destroy_unfolded = destroy_unfolded;
1187 ret->destroy_k = destroy_k;
1189 fprintf(stderr, "generated domain %p. state %d\n", ret, ret->state);
1194 void destroy_domain(domain_t *domain)
1197 //printf("domain %p & %p\n", *domain, domain);
1198 if (domain->destroy_folded)
1199 (*domain->destroy_folded)(domain->folded_params);
1200 if (domain->destroy_unfolded)
1201 (*domain->destroy_unfolded)(domain->unfolded_params);
1202 if (domain->destroy_k)
1203 (*domain->destroy_k)(domain->k_params);
1209 <<destroy domain list>>=
1210 void destroy_domain_list(list_t *domain_list)
1212 domain_list = head(domain_list);
1213 while (domain_list != NULL)
1214 destroy_domain((domain_t *) pop(&domain_list));
1218 \subsection{Domain model and group handling}
1220 <<group functions>>=
1221 <<get tension handler>>
1223 <<int comparison function>>
1224 <<find possible groups>>
1227 <<get active groups>>
1230 <<get tension handler>>=
1231 one_dim_fn_t *get_tension_handler(domain_t *domain)
1233 if (domain->state == DS_FOLDED)
1234 return domain->folded_handler;
1236 if (domain->state != DS_UNFOLDED)
1237 fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
1238 assert(domain->state == DS_UNFOLDED);
1239 return domain->unfolded_handler;
1245 int get_group(domain_t *domain)
1247 if (domain->state == DS_FOLDED)
1248 return domain->folded_group;
1250 if (domain->state != DS_UNFOLDED)
1251 fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
1252 assert(domain->state == DS_UNFOLDED);
1253 return domain->unfolded_group;
1258 We already know all possible tension classes, since they are all
1259 hardcoded into \prog. However, there may be any number of possible
1260 groups. We define a function [[find_possible_groups]] to search for
1261 possible (not neccessarily active) groups. It can search for
1262 subgroups of a single [[handler]], or by setting [[handler]] to
1263 [[NULL]], subgroups of \emph{any} handler. It returns a list of group
1265 <<find possible groups>>=
1266 list_t *find_possible_groups(list_t *list, one_dim_fn_t *handler) {
1269 while (list != NULL) {
1270 if (handler == NULL || get_tension_handler(D(list)) == handler) {
1271 push(&ret, &D(list)->folded_group, 1);
1272 push(&ret, &D(list)->unfolded_group, 1);
1278 return ret; /* no groups with this handler, no processing needed */
1280 /* sort the ret list, and remove duplicates */
1281 sort(&ret, &int_comparison);
1282 unique(&ret, &int_comparison);
1284 fprintf(stderr, "handler %p possible groups:", handler);
1286 while (list != NULL) {
1287 fprintf(stderr, " %d", *((int *)list->d));
1290 fprintf(stderr, "\n");
1296 Search a [[list]] of domains, and determine whether a particular model
1297 class and group number combination is active.
1298 <<is group active>>=
1299 int is_group_active(list_t *list, one_dim_fn_t *handler, int group)
1302 while (list != NULL) {
1303 if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group)
1311 Search a [[list]] of domains, and return all domains matching a
1312 particular model class and group number.
1314 list_t *get_group_list(list_t *list, one_dim_fn_t *handler, int group)
1318 while (list != NULL) {
1319 if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group) {
1320 push(&ret, D_TP(list), 1); /* add a pointer to the appropriate tension parameters to our new list. */
1322 fprintf(stderr, "push domain %p %s tension parameters %p onto active group %p %d list %p\n", D(list), D(list)->state == DS_FOLDED ? "folded" : "unfolded", D_TP(list), handler, group, ret);
1330 Because all the node data in lists returned by [[get_group_list]] is
1331 also in the main domain list, you shouldn't destroy the node data
1332 popped off when destroying the group lists. It will all get cleaned
1333 up when the main domain list is destroyed.
1335 Put all this together to scan through a list of domains and construct
1336 an array of all the active groups.
1337 <<get active groups>>=
1338 void get_active_groups(list_t *list,
1339 int num_tension_handlers, one_dim_fn_t **tension_handlers,
1340 int *pNum_active_groups, one_dim_fn_t ***pPer_group_handlers, void ***pActive_groups)
1342 void **active_groups=NULL;
1343 one_dim_fn_t *handler, **per_group_handlers=NULL;
1344 int i, num_active_groups, *group;
1345 list_t *possible_group_numbers, *active_group_list, *active_groups_list=NULL, *per_group_handlers_list=NULL;
1347 /* get all the active groups in a list */
1348 for (i=0; i<num_tension_handlers; i++) {
1349 handler = tension_handlers[i];
1350 if (handler == NULL) continue; /* tensionless handler */
1351 possible_group_numbers = head(find_possible_groups(list, handler));
1352 while (possible_group_numbers != NULL) {
1353 group = (int *)pop(&possible_group_numbers);
1354 if (is_group_active(list, handler, *group) == 1) {
1355 active_group_list = get_group_list(list, handler, *group);
1356 push(&active_groups_list, active_group_list, 1);
1357 push(&per_group_handlers_list, handler, 1);
1359 fprintf(stderr, "active group %p for %p %d headed by tension parameters %p\n", active_group_list, handler, *group, head(active_group_list)->d);
1360 list_print(stderr, active_group_list, "active group");
1366 /* allocate the array we'll be returning */
1367 num_active_groups = list_length(active_groups_list);
1368 active_groups = (void **) malloc(sizeof(void *)*num_active_groups);
1369 assert(active_groups != NULL);
1370 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1371 assert(per_group_handlers != NULL);
1373 /* populate the array we'll be returning */
1374 active_groups_list = head(active_groups_list);
1375 for (i=0; i<num_active_groups; i++) {
1376 handler = pop(&per_group_handlers_list);
1377 per_group_handlers[i] = handler;
1379 active_groups[i] = (void *) malloc(sizeof(tension_handler_data_t));
1380 assert(active_groups[i] != NULL);
1381 active_group_list = pop(&active_groups_list);
1382 ((tension_handler_data_t *)active_groups[i])->group = active_group_list;
1383 ((tension_handler_data_t *)active_groups[i])->env = NULL;
1384 ((tension_handler_data_t *)active_groups[i])->persist = NULL;
1386 assert(active_groups_list == NULL);
1387 assert(per_group_handlers_list == NULL);
1389 *pNum_active_groups = num_active_groups;
1390 *pPer_group_handlers = per_group_handlers;
1391 *pActive_groups = active_groups;
1396 \section{String parsing}
1398 For handling command line arguments, we need parse delimited strings (e.g. [["param1,param2,param3..."]]).
1399 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1400 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1401 need to take care of parsing those parameters themselves.
1402 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1406 <<parse definitions>>
1407 <<parse declarations>>
1410 <<parse module makefile lines>>=
1412 notangle -Rparse.c $^ > $@
1414 notangle -Rparse.h $^ > $@
1415 check_parse.c : sawsim.nw
1416 notangle -Rcheck-parse.c $^ > $@
1417 check_parse : check_parse.c parse.c parse.h
1418 gcc -g -o $@ $< parse.c -lcheck
1420 rm -f parse.c parse.h check_parse.c check_parse
1423 <<parse definitions>>=
1424 #define SEP ',' /* argument separator character */
1427 <<parse declarations>>=
1428 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1429 int *num, char ***string_array);
1430 extern void free_parsed_list(int num, char **string_array);
1433 [[parse_list_string]] allocates memory, don't forget to free it
1434 afterward with [[free_parsed_list]]. It does not alter the original.
1436 The string may start off with a [[deeper]] character
1437 (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave
1438 [[x,y]], where the pointer is one character in on the copied string.
1439 However, when we go to free the memory, we need a pointer to the
1440 beginning of the string. In order to accommodate this for a string
1441 with $N$ argument, allocate a pointer array with $N+1$ elements, let
1442 the first $N$ elements point to the separated arguments, and let the
1443 last element point to the start of the copied string regardless of
1445 <<parse delimited list string functions>>=
1446 /* TODO, split out into parse.hc */
1447 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1450 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1451 if (string[i] == deeper) {depth++;}
1452 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1458 void parse_list_string(char *string, char sep, char deeper, char shallower,
1459 int *num, char ***string_array)
1461 char *str=NULL, **ret=NULL;
1463 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
1465 *string_array = NULL;
1468 /* make a copy of the string, so we don't change the original */
1469 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1470 assert(str != NULL);
1471 strcpy(str, string); /* we know str is long enough */
1472 /* count the number of regions, so we can allocate pointers to them */
1475 n++; i++; /* move on to next argument */
1476 i += next_delim_index(str+i, sep, deeper, shallower);
1477 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1478 } while (str[i] != '\0');
1479 ret = (char **)malloc(sizeof(char *)*(n+1));
1480 assert(ret != NULL);
1481 /* replace the separators with '\0' & assign pointers */
1482 ret[n] = str; /* point to the front of the copied string */
1485 for(i=1; i<n; i++) {
1486 j += next_delim_index(str+j, sep, deeper, shallower);
1488 ret[i] = str+j; /* point to the separated arguments */
1490 /* strip off bounding braces, if any */
1491 for(i=0; i<n; i++) {
1492 if (ret[i][0] == deeper) {
1496 if (ret[i][strlen(ret[i])-1] == shallower)
1497 ret[i][strlen(ret[i])-1] = '\0';
1500 *string_array = ret;
1503 void free_parsed_list(int num, char **string_array)
1506 /* frees all items, since they were allocated as a single string*/
1507 free(string_array[num]);
1508 /* and free the array of pointers */
1514 \subsection{Parsing implementation}
1520 <<parse delimited list string functions>>
1524 #include <assert.h> /* assert() */
1525 #include <stdlib.h> /* NULL */
1526 #include <stdio.h> /* fprintf(), stdout *//*!!*/
1527 #include <string.h> /* strlen() */
1531 \subsection{Parsing unit tests}
1533 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1535 <<parse check includes>>
1536 <<string comparison definition and globals>>
1537 <<check relative error>>
1538 <<parse test suite>>
1539 <<main check program>>
1542 <<parse check includes>>=
1543 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
1544 #include <stdio.h> /* printf() */
1545 #include <assert.h> /* assert() */
1546 #include <string.h> /* strlen() */
1551 <<parse test suite>>=
1552 <<parse list string tests>>
1553 <<parse suite definition>>
1556 <<parse suite definition>>=
1557 Suite *test_suite (void)
1559 Suite *s = suite_create ("k model");
1560 <<parse list string test case defs>>
1562 <<parse list string test case add>>
1567 <<parse list string tests>>=
1570 START_TEST(test_next_delim_index)
1572 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
1573 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
1574 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
1575 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
1576 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
1580 START_TEST(test_parse_list_null)
1584 parse_list_string(NULL, SEP, '{', '}',
1585 &num_param_args, ¶m_args);
1586 fail_unless(num_param_args == 0, NULL);
1587 fail_unless(param_args == NULL, NULL);
1590 START_TEST(test_parse_list_single_simple)
1594 parse_list_string("arg", SEP, '{', '}',
1595 &num_param_args, ¶m_args);
1596 fail_unless(num_param_args == 1, NULL);
1597 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
1600 START_TEST(test_parse_list_single_compound)
1604 parse_list_string("{x,y,z}", SEP, '{', '}',
1605 &num_param_args, ¶m_args);
1606 fail_unless(num_param_args == 1, NULL);
1607 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
1610 START_TEST(test_parse_list_double_simple)
1614 parse_list_string("abc,def", SEP, '{', '}',
1615 &num_param_args, ¶m_args);
1616 fail_unless(num_param_args == 2, NULL);
1617 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
1618 fail_unless(STRMATCH(param_args[1],"def"), NULL);
1622 <<parse list string test case defs>>=
1623 TCase *tc_parse_list_string = tcase_create("parse list string");
1625 <<parse list string test case add>>=
1626 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
1627 tcase_add_test(tc_parse_list_string, test_parse_list_null);
1628 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
1629 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
1630 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
1631 suite_add_tcase(s, tc_parse_list_string);
1635 \section{Unit tests}
1637 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1646 <<main check program>>
1658 <<determine dt tests>>
1660 <<does domain unfold tests>>
1661 <<randomly unfold domains tests>>
1662 <<suite definition>>
1665 <<suite definition>>=
1666 Suite *test_suite (void)
1668 Suite *s = suite_create ("sawsim");
1669 <<F test case defs>>
1670 <<determine dt test case defs>>
1671 <<happens test case defs>>
1672 <<does domain unfold test case defs>>
1673 <<randomly unfold domains test case defs>>
1676 <<determine dt test case add>>
1677 <<happens test case add>>
1678 <<does domain unfold test case add>>
1679 <<randomly unfold domains test case add>>
1682 tcase_add_checked_fixture(tc_strip_address,
1683 setup_strip_address,
1684 teardown_strip_address);
1690 <<main check program>>=
1694 Suite *s = test_suite();
1695 SRunner *sr = srunner_create(s);
1696 srunner_run_all(sr, CK_ENV);
1697 number_failed = srunner_ntests_failed(sr);
1699 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
1703 \subsection{F tests}
1705 <<worm-like chain tests>>
1707 <<F test case defs>>=
1708 <<worm-like chain test case def>>
1710 <<F test case add>>=
1711 <<worm-like chain test case add>>
1714 <<worm-like chain tests>>=
1715 START_TEST(test_wlc_at_zero)
1717 double T=1.0, L=1.0, p=0.1, x=0.0;
1718 fail_unless(wlc(x, T, p, L)==0, NULL);
1722 START_TEST(test_wlc_at_half)
1724 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
1725 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
1726 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
1728 * linear term = x/L = 0.5
1729 * nonlinear + linear = 0.75 + 0.5 = 1.25
1730 * wlc = 10e21*1.25 = 12.5e21
1732 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
1733 "wlc(%g, %g, %g, %g) = %g != %g",
1734 x, T, p, L, wlc(x, T, p, L), 12.5e21);
1738 <<worm-like chain test case def>>=
1739 TCase *tc_wlc = tcase_create("WLC");
1742 <<worm-like chain test case add>>=
1743 tcase_add_test(tc_wlc, test_wlc_at_zero);
1744 tcase_add_test(tc_wlc, test_wlc_at_half);
1745 suite_add_tcase(s, tc_wlc);
1748 \subsection{Model tests}
1750 Check the searching with [[linear_k]].
1751 Check overwhelming force treatment with the heavyside-esque step [[?]].
1752 <<determine dt tests>>=
1753 double linear_k(double F, environment_t *env, void *params)
1755 double Fnot = *(double *)params;
1759 START_TEST(test_determine_dt_linear_k)
1762 double dt_max=3.0, Fnot=3.0;
1763 double F[]={0,1,2,3,4,5,6};
1764 domain_t dom; /* use both parts at once for folded/unfolded */
1768 dom->next = dom->prev = NULL;
1769 dom->k_func_t = linear_k;
1770 dom->folded_params = &Fnot;
1771 dom->unfolded_params = !!!!!!!!!
1772 dom->destroy_folded = dom->destroy_unfolded = NULL;
1773 for( i=0; i < sizeof(F)/sizeof(double); i++) {
1774 fail_unless(determine_dt(folded, unfolded, targ_p, dt_max, v, x, T)==1, NULL);
1780 <<determine dt test case defs>>=
1781 TCase *tc_determine_dt = tcase_create("Determine dt");
1783 <<determine dt test case add>>=
1784 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
1785 suite_add_tcase(s, tc_determine_dt);
1790 <<happens test case defs>>=
1792 <<happens test case add>>=
1795 <<does domain unfold tests>>=
1797 <<does domain unfold test case defs>>=
1799 <<does domain unfold test case add>>=
1802 <<randomly unfold domains tests>>=
1804 <<randomly unfold domains test case defs>>=
1806 <<randomly unfold domains test case add>>=
1810 \section{Balancing group extensions}
1811 \label{app.tension_balance}
1813 For a given total extension $x$ (of the piezo), the various domain
1814 models (WLC, FJC, Hookean springs, \ldots) and groups extend different
1815 amounts, and we need to tweak the portion that each extends to balance
1816 the tension amongst the active groups. Since the tension balancing is
1817 separable from the bulk of the simulation, we break it out into a
1818 separate module. The interface is defined in [[tension_balance.h]],
1819 the implementation in [[tension_balance.c]], and the unit testing in
1820 [[check_tension_balance.c]]
1822 <<tension-balance.h>>=
1824 <<tension balance function declaration>>
1827 <<tension balance functions>>=
1828 <<one dimensional bracket>>
1829 <<one dimensional solve>>
1831 <<tension balance function>>
1834 <<tension balance module makefile lines>>=
1835 tension_balance.c : sawsim.nw
1836 notangle -Rtension-balance.c $^ > $@
1837 tension_balance.h : sawsim.nw
1838 notangle -Rtension-balance.h $^ > $@
1839 check_tension_balance.c : sawsim.nw
1840 notangle -Rcheck-tension-balance.c $^ > $@
1841 check_tension_balance : check_tension_balance.c global.h tension_balance.c tension_balance.h
1842 gcc -g -o $@ $< tension_balance.c -lcheck
1844 rm -f tension_balance.c tension_balance.h check_tension_balance.c check_tension_balance
1847 The entire force balancing problem reduces to a solving two nested
1848 one-dimensional root-finding problems. First we define the one
1849 dimensional tension $F(x_0)$ (where $i=0$ is the index of the first
1850 non-empty group). $x(x_0)$ is also strictly monotonically increasing,
1851 so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]]
1852 stores the last successful value of $x$, since we expect to be taking
1853 small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]]
1854 indicates that the groups have changed.
1855 <<tension balance function declaration>>=
1856 double tension_balance(int num_tension_groups,
1857 one_dim_fn_t **tension_handlers,
1862 <<tension balance function>>=
1863 double tension_balance(int num_tension_groups,
1864 one_dim_fn_t **tension_handlers,
1869 static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL};
1870 double F, x_guess, xo, lb, ub;
1871 double min_dx=1e-10, min_dy=1e-15;
1872 int max_steps=100, i;
1874 assert(num_tension_groups > 0);
1875 assert(tension_handlers != NULL);
1876 assert(params != NULL);
1877 assert(num_tension_groups > 0);
1879 if (num_tension_groups == 1) { /* only one group, no balancing required */
1882 //printf("balancing for x = %g with ", x);
1883 //for(i=0; i<num_tension_groups; i++) printf(" %d: %g\t", i, x_xo_data.xi[i]);
1885 if (last_x == -1) { /* new group setup, reset x_xo_data */
1886 /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */
1887 if (x_xo_data.xi != NULL)
1889 /* construct new x_xo_data */
1890 x_xo_data.n_groups = num_tension_groups;
1891 x_xo_data.tension_handler = tension_handlers;
1892 x_xo_data.handler_data = params;
1893 x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
1894 for (i=0; i<num_tension_groups; i++)
1895 x_xo_data.xi[i] = -1.0;
1897 if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
1898 if (x == 0) x_guess = 1.0;
1899 else x_guess = x/num_tension_groups;
1901 fprintf(stderr, "search bracket for %g with guess of %g\n", x, x_guess);
1903 oneD_bracket(x_of_xo, &x_xo_data, x, x_guess, &lb, &ub);
1904 } else { /* work off the last balanced point */
1906 fprintf(stderr, "step off the old bracketing x %g + %g = %g\n", last_x, x-last_x, x);
1909 lb = x_xo_data.xi[0];
1910 ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */
1911 } else if (x < last_x) {
1912 lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
1913 ub = x_xo_data.xi[0];
1914 } else { /* x == last_x */
1915 printf("not moving\n");
1916 lb= x_xo_data.xi[0];
1917 ub= x_xo_data.xi[0];
1920 //printf("lb %g,\tub %g\n", lb, ub);
1921 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_dx, min_dy, max_steps, NULL);
1922 if (num_tension_groups > 1 && 0) {
1923 printf("balanced for x = %g with ", x);
1924 for(i=0; i<num_tension_groups; i++) printf(" %d: %g\t", i, x_xo_data.xi[i]);
1928 F = (*tension_handlers[0])(xo, params[0]);
1933 <<tension balance internal definitions>>=
1934 <<x of x0 definitions>>
1937 <<x of x0 definitions>>=
1938 typedef struct x_of_xo_data_struct {
1940 one_dim_fn_t **tension_handler; /* array of fn pointers */
1941 void **handler_data; /* array of void* pointers */
1947 double x_of_xo(double xo, void *pdata)
1949 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
1950 double F, x=0, xi, xi_guess, lb, ub;
1952 double min_dx=1e-10, min_dy=1e-14;
1954 assert(data->n_groups > 0);
1956 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
1958 if (data->xi) data->xi[0] = xo;
1959 for (i=1; i < data->n_groups; i++) {
1960 if (data->xi[i] == -1) xi_guess = 1.0;
1961 else xi_guess = data->xi[i];
1962 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
1963 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_dx, min_dy, max_steps, NULL);
1966 if (data->xi) data->xi[i] = xi;
1972 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited number of steps for monotonically increasing $f(x)$.
1973 Simple bisection, so it's robust and converges fairly quickly.
1974 <<one dimensional function definition>>=
1975 /* equivalent to gsl_func_t */
1976 typedef double one_dim_fn_t(double x, void *params);
1978 <<one dimensional solve declaration>>=
1979 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
1980 double min_dx, double min_dy, int max_steps,
1984 <<one dimensional solve>>=
1985 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
1986 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
1987 double min_dx, double min_dy, int max_steps,
1990 double yx, ylb, yub, x;
1993 ylb = (*f)(lb, params);
1994 yub = (*f)(ub, params);
1996 /* check some simple cases */
1998 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bounded */
1999 else return lb; /* any x on the range [lb, ub] would work */
2001 if (ylb == y) { x=lb; yx=ylb; goto end; }
2002 if (yub == y) { x=ub; yx=yub; goto end; }
2004 //printf("lb %g, x %g, ub %g\tylb %g, y %g, yub %g\n", lb, x, ub, ylb, y, yub);
2005 assert(ylb < y); assert(yub > y);
2006 for (i=0; i<max_steps; i++) {
2008 yx = (*f)(x, params);
2009 if (ub-lb < min_dx || yub-ylb < min_dy || yx == y) break;
2010 else if (yx > y) { ub=x; yub=yx; }
2011 else /* yx < y */ { lb=x; ylb=yx; }
2019 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2020 Generate bracketing $x$ values through bisection or exponential growth.
2021 <<one dimensional bracket declaration>>=
2022 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2025 <<one dimensional bracket>>=
2026 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2028 double yx, step, x=xguess;
2030 yx = (*f)(x, params);
2031 //fprintf(stdout, "bracketing %g, start at f(%g) = %g\n", y, x, yx);
2038 if (x == 0) x = 0.5; /* guess a scale of 1.0 */
2041 yx = (*f)(x, params);
2042 //fprintf(stdout, "increasing to f(%g) = %g\n", x, yx);
2046 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2050 \subsection{Balancing implementation}
2052 <<tension-balance.c>>=
2054 <<tension balance includes>>
2055 #include "tension_balance.h"
2056 <<tension balance internal definitions>>
2057 <<tension balance functions>>
2060 <<tension balance includes>>=
2061 #include <assert.h> /* assert() */
2062 #include <stdlib.h> /* NULL */
2063 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2064 #include <stdio.h> /* fprintf(), stdout */
2068 \subsection{Balancing unit tests}
2070 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2071 <<check-tension-balance.c>>=
2072 <<tension balance check includes>>
2073 <<tension balance test suite>>
2074 <<main check program>>
2077 <<tension balance check includes>>=
2078 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2079 #include <assert.h> /* assert() */
2082 #include "tension_balance.h"
2085 <<tension balance test suite>>=
2086 <<tension balance function tests>>
2087 <<tension balance suite definition>>
2090 <<tension balance suite definition>>=
2091 Suite *test_suite (void)
2093 Suite *s = suite_create ("tension balance");
2094 <<tension balance function test case defs>>
2096 <<tension balance function test case adds>>
2101 <<tension balance function tests>>=
2102 <<check relative error>>
2104 double hooke(void *pK, double x)
2107 return *((double*)pK) * x;
2110 START_TEST(test_single_function)
2112 double k=5, x=3, last_x=2.0, F;
2113 one_dim_fn_t *handlers[] = {&hooke};
2114 void *data[] = {&k};
2115 double xi[] = {0.0};
2117 int new_active_groups = 1;
2118 F = tension_balance(1, handlers, data, active, new_active_groups, xi, last_x, x);
2119 fail_unless(F = k*x, NULL);
2124 We can also test balancing two springs with different spring constants.
2125 The analytic solution for a general number of springs is given in Appendix \ref{app.math_hooke}.
2126 <<tension balance function tests>>=
2127 START_TEST(test_double_hooke)
2129 double k1=5, k2=4, x=3, last_x=2.0, F, Fe, x1e, x2e;
2130 one_dim_fn_t *handlers[] = {&hooke, &hooke, NULL};
2131 void *data[] = {&k1, &k2, NULL};
2132 double xi[] = {0.0, 0.0, 0.0};
2133 int active[] = {1, 1, 0};
2134 int new_active_groups = 1;
2135 F = tension_balance(3, handlers, data, active, new_active_groups, xi, last_x, x);
2139 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2140 CHECK_ERR(1e-6, x1e, xi[0]);
2141 CHECK_ERR(1e-6, x2e, xi[1]);
2142 CHECK_ERR(1e-6, Fe, F);
2147 <<tension balance function test case defs>>=
2148 TCase *tc_tbfunc = tcase_create("tension balance function");
2151 <<tension balance function test case adds>>=
2152 tcase_add_test(tc_tbfunc, test_single_function);
2153 tcase_add_test(tc_tbfunc, test_double_hooke);
2154 suite_add_tcase(s, tc_tbfunc);
2159 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2160 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2161 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2165 <<list definitions>>
2166 <<list declarations>>
2169 <<list declarations>>=
2170 <<head and tail declarations>>
2171 <<list length declaration>>
2172 <<push declaration>>
2174 <<list destroy declaration>>
2175 <<list to array declaration>>
2176 <<move declaration>>
2177 <<sort declaration>>
2178 <<unique declaration>>
2179 <<list print declaration>>
2183 <<create and destroy node>>
2196 <<list module makefile lines>>=
2198 notangle -Rlist.c $^ > $@
2200 notangle -Rlist.h $^ > $@
2201 check_list.c : sawsim.nw
2202 notangle -Rcheck-list.c $^ > $@
2203 check_list : check_list.c global.h list.c list.h
2204 gcc -g -o $@ $< list.c -lcheck
2206 rm -f list.c list.h check_list.c check_list
2209 <<list definitions>>=
2210 typedef struct list_struct {
2211 struct list_struct *next;
2212 struct list_struct *prev;
2216 <<comparison function typedef>>
2219 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2220 <<head and tail declarations>>=
2221 list_t *head(list_t *list);
2222 list_t *tail(list_t *list);
2225 list_t *head(list_t *list)
2229 while (list->prev != NULL) {
2235 list_t *tail(list_t *list)
2239 while (list->next != NULL) {
2246 <<list length declaration>>=
2247 int list_length(list_t *list);
2250 int list_length(list_t *list)
2257 while (list->next != NULL) {
2265 [[push]] inserts a node relative to the current position in the list
2266 without changing the current position.
2267 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2268 so we set it to that of the pushed domain.
2269 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2270 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2271 <<push declaration>>=
2272 void push(list_t **pList, void *data, int below);
2275 void push(list_t **pList, void *data, int below)
2277 list_t *list, *new_node;
2278 assert(pList != NULL);
2280 new_node = create_node(data);
2283 else if (below==0) { /* insert above */
2284 if (list->prev != NULL)
2285 list->prev->next = new_node;
2286 new_node->prev = list->prev;
2287 list->prev = new_node;
2288 new_node->next = list;
2289 } else { /* insert below */
2290 if (list->next != NULL)
2291 list->next->prev = new_node;
2292 new_node->next = list->next;
2293 list->next = new_node;
2294 new_node->prev = list;
2299 [[pop]] removes the current domain node, moving the current position
2300 to the node after it, unless that node is [[NULL]], in which case move
2301 the current position to the node before it.
2302 <<pop declaration>>=
2303 void *pop(list_t **pList);
2306 void *pop(list_t **pList)
2308 list_t *list, *popped;
2310 assert(pList != NULL);
2312 assert(list != NULL); /* not an empty list */
2314 /* bypass the popped node */
2315 if (list->prev != NULL)
2316 list->prev->next = list->next;
2317 if (list->next != NULL) {
2318 list->next->prev = list->prev;
2319 *pList = list->next;
2321 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
2323 destroy_node(popped);
2328 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
2329 If the cleanup function is [[NULL]], no cleanup function is called.
2330 The nodes are not popped in any particular order.
2331 <<list destroy declaration>>=
2332 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
2335 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
2339 assert(pList != NULL);
2342 assert(list != NULL); /* not an empty list */
2343 while (list != NULL) {
2345 if (destroy != NULL)
2351 [[list_to_array]] converts a list to an array of pointers, preserving order.
2352 <<list to array declaration>>=
2353 void list_to_array(list_t *list, int *length, void ***array);
2356 void list_to_array(list_t *list, int *pLength, void ***pArray)
2360 assert(list != NULL);
2361 assert(pLength != NULL);
2362 assert(pArray != NULL);
2364 length = list_length(list);
2365 array = (void **)malloc(sizeof(void *)*length);
2366 assert(array != NULL);
2369 while (list != NULL) {
2370 array[i++] = list->d;
2378 [[move]] moves the current element along the list in either direction.
2379 <<move declaration>>=
2380 void move(list_t **pList, int downstream);
2383 void move(list_t **pList, int downstream)
2385 list_t *list, *flipper;
2387 assert(pList != NULL);
2388 list = *pList; /* a>B>c>d */
2389 assert(list != NULL); /* not an empty list */
2390 if (downstream == 0)
2391 flipper = list->prev; /* flipper is A */
2393 flipper = list->next; /* flipper is C */
2394 /* ensure there is somewhere to go in the direction we're heading */
2395 assert(flipper != NULL);
2396 /* remove the current element */
2397 data = pop(&list); /* data=B, a>C>d */
2398 /* and put it back in in the correct spot */
2400 if (downstream == 0) {
2401 push(&list, data, 0); /* b>A>c>d */
2402 list = list->prev; /* B>a>c>d */
2404 push(&list, data, 1); /* a>C>b>d */
2405 list = list->next; /* a>c>B>d */
2411 [[sort]] sorts a list from smallest to largest according to a user
2412 specified comparison.
2413 <<comparison function typedef>>=
2414 typedef int (comparison_fn_t)(void *A, void *B);
2417 <<int comparison function>>=
2418 int int_comparison(void *A, void *B)
2423 if (a > b) return 1;
2424 else if (a == b) return 0;
2428 <<sort declaration>>=
2429 void sort(list_t **pList, comparison_fn_t *comp);
2431 Here we use the bubble sort algorithm.
2433 void sort(list_t **pList, comparison_fn_t *comp)
2437 assert(pList != NULL);
2439 assert(list != NULL); /* not an empty list */
2440 while (stable == 0) {
2443 while (list->next != NULL) {
2444 if ((*comp)(list->d, list->next->d) > 0) {
2446 move(&list, 1 /* downstream */);
2455 [[unique]] removes duplicates from a sorted list according to a user
2456 specified comparison. Don't do this unless you have other pointers
2457 any dynamically allocated data in your list, because [[unique]] just
2458 drops any non-unique data without freeing it.
2459 <<unique declaration>>=
2460 void unique(list_t **pList, comparison_fn_t *comp);
2463 void unique(list_t **pList, comparison_fn_t *comp)
2466 assert(pList != NULL);
2468 assert(list != NULL); /* not an empty list */
2470 while (list->next != NULL) {
2471 if ((*comp)(list->d, list->next->d) == 0) {
2480 [[list_print]] prints a list to a [[FILE *]] stream.
2481 <<list print declaration>>=
2482 void list_print(FILE *file, list_t *list, char *list_name);
2485 void list_print(FILE *file, list_t *list, char *list_name)
2487 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
2489 while (list != NULL) {
2490 fprintf(file, " %p", list->d);
2493 fprintf(file, "\n");
2497 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
2498 <<create and destroy node>>=
2499 list_t *create_node(void *data)
2501 list_t *ret = (list_t *)malloc(sizeof(list_t));
2502 assert(ret != NULL);
2509 void destroy_node(list_t *node)
2515 The user must free the data pointed to by the node on their own.
2517 \subsection{List implementation}
2527 #include <assert.h> /* assert() */
2528 #include <stdlib.h> /* malloc(), free(), rand() */
2529 #include <stdio.h> /* fprintf(), stdout, FILE */
2530 #include "global.h" /* destroy_data_func_t */
2533 \subsection{List unit tests}
2535 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2537 <<list check includes>>
2539 <<main check program>>
2542 <<list check includes>>=
2543 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2548 <<list test suite>>=
2551 <<list suite definition>>
2554 <<list suite definition>>=
2555 Suite *test_suite (void)
2557 Suite *s = suite_create ("list");
2558 <<push test case defs>>
2559 <<pop test case defs>>
2561 <<push test case adds>>
2562 <<pop test case adds>>
2568 START_TEST(test_push)
2573 push(&list, (void *)i, 1);
2574 fail_unless(list == head(list), NULL);
2575 fail_unless( (int)list->d == 0 );
2576 for (i=0; i<n; i++) {
2577 /* we expect 0, n-1, n-2, ..., 1 */
2580 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
2586 <<push test case defs>>=
2587 TCase *tc_push = tcase_create("push");
2590 <<push test case adds>>=
2591 tcase_add_test(tc_push, test_push);
2592 suite_add_tcase(s, tc_push);
2597 <<pop test case defs>>=
2599 <<pop test case adds>>=
2602 \section{Function string evaluation}
2604 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).
2605 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
2606 The solution is to run a scripting language as a subprocess accessed via pipes.
2607 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
2609 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
2610 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
2611 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.
2612 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
2614 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.
2615 Then you can either statically or dynamically link to those hardcoded functions.
2616 While much less flexible, this approach would be a fairly simple-to-implement fallback.
2618 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
2619 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
2623 <<string eval setup declaration>>
2624 <<string eval function declaration>>
2625 <<string eval teardown declaration>>
2628 <<string eval module makefile lines>>=
2629 string_eval.c : sawsim.nw
2630 notangle -Rstring-eval.c $^ > $@
2631 string_eval.h : sawsim.nw
2632 notangle -Rstring-eval.h $^ > $@
2633 check_string_eval.c : sawsim.nw
2634 notangle -Rcheck-string-eval.c $^ > $@
2635 check_string_eval : check_string_eval.c string_eval.c string_eval.h
2636 gcc -g -o $@ $< string_eval.c -lcheck -lgsl -lgslcblas -lm
2638 rm -f string_eval.c string_eval.h check_string_eval.c check_string_eval
2641 For an introduction to POSIX process control, see\\
2642 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
2643 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
2644 and of course, the relavent [[man]] pages.
2646 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
2647 [[execvp]] replaces the calling process' program with a new program.
2648 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
2649 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
2650 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
2651 The new program needs command line arguments, just like it would if you were running it from a shell.
2652 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
2653 with the final array entry being a [[NULL]] pointer.
2655 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
2656 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
2657 <<string eval subprocess definitions>>=
2658 #define SUBPROCESS "python"
2659 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
2660 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
2661 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
2663 The [[i]] option lets Python know that it should run in interactive mode.
2664 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
2665 In interactive mode, python acts on every instruction as soon as it is recieved.
2666 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.
2667 %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.
2669 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
2670 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
2671 [[fork]] returns two copies of the same program, executing the original code.
2672 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
2673 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
2675 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.
2676 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
2677 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
2678 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
2679 <<string eval pipe definitions>>=
2680 #define PIPE_READ 0 /* the end you read from */
2681 #define PIPE_WRITE 1 /* the end you write to */
2683 #define STDIN 0 /* initial index of stdin pair */
2684 #define STDOUT 2 /* initial index of stdout pair */
2687 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
2689 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.
2691 <<string eval setup declaration>>=
2692 extern void string_eval_setup(FILE **pIn, FILE **pOut);
2694 <<string eval setup definition>>=
2695 void string_eval_setup(FILE **pIn, FILE **pOut)
2698 int pfd[4]; /* file descriptors for stdin and stdout */
2700 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
2701 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
2703 pid = fork(); /* split process into two copies */
2705 if (pid == -1) { /* fork error */
2706 perror("fork error");
2708 } else if (pid == 0) { /* child */
2709 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
2710 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
2711 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
2712 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
2713 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
2714 perror("exec error"); /* exec shouldn't return */
2716 } else { /* parent */
2717 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
2718 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
2719 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
2720 if ( *pIn == NULL ) {
2721 perror("fdopen (in)");
2724 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
2725 if ( *pOut == NULL ) {
2726 perror("fdopen (out)");
2733 To use the evaluating subprocess, we just pipe in our command, and read out the results.
2734 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
2735 <<string eval function declaration>>=
2736 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
2738 <<string eval function definition>>=
2739 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
2742 rc = fprintf(in, "%s", input);
2743 assert(rc == strlen(input));
2746 alarm(1); /* set a one second timeout on the read */
2747 assert( fgets(output, buflen, out) != NULL );
2748 alarm(0); /* cancel the timeout */
2749 //fprintf(stderr, "eval: %s ----> %s", input, output);
2752 The [[alarm]] calls set and clear a timeout on the returned output.
2753 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.
2754 This protects against invalid input for which a line of output is not printed to [[stdout]].
2755 Other invalid inputs (e.g. those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
2756 If you are getting strange results, check your python code seperately. TODO, better error handling.
2758 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
2759 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
2760 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
2761 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
2762 <<string eval teardown declaration>>=
2763 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
2766 <<string eval teardown definition>>=
2767 void string_eval_teardown(FILE **pIn, FILE **pOut)
2772 /* redirect Python's stderr */
2773 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
2777 assert( fclose(*pIn) == 0 );
2779 assert( fclose(*pOut) == 0 );
2782 /* wait for python to exit */
2784 pid = wait(&stat_loc);
2791 if (WIFEXITED(stat_loc)) {
2792 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
2793 } else if (WIFSIGNALED(stat_loc)) {
2794 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
2799 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
2801 \subsection{String evaluation implementation}
2805 <<string eval includes>>
2806 #include "string_eval.h"
2807 <<string eval internal definitions>>
2808 <<string eval functions>>
2811 <<string eval includes>>=
2812 #include <assert.h> /* assert() */
2813 #include <stdlib.h> /* NULL */
2814 #include <stdio.h> /* fprintf(), stdout, fdopen() */
2815 #include <string.h> /* strlen() */
2816 #include <sys/types.h> /* pid_t */
2817 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
2818 #include <wait.h> /* wait() */
2821 <<string eval internal definitions>>=
2822 <<string eval subprocess definitions>>
2823 <<string eval pipe definitions>>
2826 <<string eval functions>>=
2827 <<string eval setup definition>>
2828 <<string eval function definition>>
2829 <<string eval teardown definition>>
2832 \subsection{String evaluation unit tests}
2834 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2835 <<check-string-eval.c>>=
2836 <<string eval check includes>>
2837 <<string comparison definition and globals>>
2838 <<string eval test suite>>
2839 <<main check program>>
2842 <<string eval check includes>>=
2843 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
2844 #include <stdio.h> /* printf() */
2845 #include <assert.h> /* assert() */
2846 #include <string.h> /* strlen() */
2847 #include <signal.h> /* SIGKILL */
2849 #include "string_eval.h"
2852 <<string eval test suite>>=
2853 <<string eval tests>>
2854 <<string eval suite definition>>
2857 <<string eval suite definition>>=
2858 Suite *test_suite (void)
2860 Suite *s = suite_create ("string eval");
2861 <<string eval test case defs>>
2863 <<string eval test case add>>
2868 <<string eval tests>>=
2869 START_TEST(test_setup_teardown)
2872 string_eval_setup(&in, &out);
2873 string_eval_teardown(&in, &out);
2876 START_TEST(test_invalid_command)
2879 char input[80], output[80]={};
2880 string_eval_setup(&in, &out);
2881 sprintf(input, "print ABCDefg\n");
2882 string_eval(in, out, input, 80, output);
2883 string_eval_teardown(&in, &out);
2886 START_TEST(test_simple_eval)
2889 char input[80], output[80]={};
2890 string_eval_setup(&in, &out);
2891 sprintf(input, "print 3+4*5\n");
2892 string_eval(in, out, input, 80, output);
2893 fail_unless(STRMATCH(output,"23\n"), NULL);
2894 string_eval_teardown(&in, &out);
2897 START_TEST(test_multiple_evals)
2900 char input[80], output[80]={};
2901 string_eval_setup(&in, &out);
2902 sprintf(input, "print 3+4*5\n");
2903 string_eval(in, out, input, 80, output);
2904 fail_unless(STRMATCH(output,"23\n"), NULL);
2905 sprintf(input, "print (3**2 + 4**2)**0.5\n");
2906 string_eval(in, out, input, 80, output);
2907 fail_unless(STRMATCH(output,"5.0\n"), NULL);
2908 string_eval_teardown(&in, &out);
2911 START_TEST(test_eval_with_state)
2914 char input[80], output[80]={};
2915 string_eval_setup(&in, &out);
2916 sprintf(input, "print 3+4*5\n");
2917 fprintf(in, "A = 3\n");
2918 sprintf(input, "print A*3\n");
2919 string_eval(in, out, input, 80, output);
2920 fail_unless(STRMATCH(output,"9\n"), NULL);
2921 string_eval_teardown(&in, &out);
2925 <<string eval test case defs>>=
2926 TCase *tc_string_eval = tcase_create("string_eval");
2928 <<string eval test case add>>=
2929 tcase_add_test(tc_string_eval, test_setup_teardown);
2930 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
2931 tcase_add_test(tc_string_eval, test_simple_eval);
2932 tcase_add_test(tc_string_eval, test_multiple_evals);
2933 tcase_add_test(tc_string_eval, test_eval_with_state);
2934 suite_add_tcase(s, tc_string_eval);
2938 \section{Accelerating function evaluation}
2940 My first version-0.3 code was running very slowly.
2941 With the options suggested in the help
2942 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
2943 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
2944 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
2945 That's only 15 calls per solution, so the search algorithm seems reasonable.
2946 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
2949 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
2953 <<accel k module makefile lines>>=
2954 accel_k.c : sawsim.nw
2955 notangle -Raccel-k.c $^ > $@
2956 accel_k.h : sawsim.nw
2957 notangle -Raccel-k.h $^ > $@
2958 check_accel_k.c : sawsim.nw
2959 notangle -Rcheck-accel_k.c $^ > $@
2960 check_accel_k : check_accel_k.c global.h
2961 gcc -g -o $@ $< accel_k.c -lcheck -lgsl -lgslcblas -lm
2963 rm -f accel_k.c accel_k.h check_accel_k.c check_accel_k
2967 #include <assert.h> /* assert() */
2968 #include <stdlib.h> /* realloc(), free(), NULL */
2969 #include "global.h" /* environment_t */
2970 #include "k_model.h" /* k_func_t */
2971 #include "interp.h" /* interp_* */
2972 #include "accel_k.h"
2974 typedef struct accel_k_struct {
2975 interp_table_t *itable;
2981 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
2982 static int num_accels = 0;
2983 static accel_k_t *accels=NULL;
2985 /* Wrap k in the standard f(x) acceleration form */
2986 static double k_wrap(double F, void *params)
2988 accel_k_t *p = (accel_k_t *)params;
2989 return (*p->k)(F, p->env, p->params);
2992 static int k_tol(double FA, double kA, double FB, double kB)
2995 if (FB-FA > 1e-12) {
2996 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
2997 return 1; /* unnacceptable */
2999 //printf("acceptable tol\n");
3000 return 0; /* acceptable */
3004 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3007 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3008 assert(accels != NULL);
3009 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3011 accels[i].env = env;
3012 accels[i].params = params;
3019 for (i=0; i<num_accels; i++)
3020 interp_table_free(accels[i].itable);
3022 if (accels) free(accels);
3026 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3029 for (i=0; i<num_accels; i++) {
3030 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3031 /* We've already setup this function.
3032 * WARNING: we're assuming param and env are the same. */
3033 return interp_table_eval(accels[i].itable, accels+i, F);
3036 /* set up a new acceleration environment */
3037 i = add_accel_k(k, env, params);
3038 return interp_table_eval(accels[i].itable, accels+i, F);
3042 \section{Tension models}
3043 \label{sec.tension_models}
3045 TODO: tension model intro.
3046 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.
3048 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3049 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]].
3051 <<tension-model.h>>=
3053 <<tension handler types>>
3054 <<constant tension model declarations>>
3055 <<hooke tension model declarations>>
3056 <<worm-like chain tension model declarations>>
3057 <<freely-jointed chain tension model declarations>>
3058 <<find tension definitions>>
3059 <<tension model global declarations>>
3062 <<tension model module makefile lines>>=
3063 tension_model.c : sawsim.nw
3064 notangle -Rtension-model.c $^ > $@
3065 tension_model.h : sawsim.nw
3066 notangle -Rtension-model.h $^ > $@
3067 check_tension_model.c : sawsim.nw
3068 notangle -Rcheck-tension-model.c $^ > $@
3069 check_tension_model : check_tension_model.c global.h tension_model.c tension_model.h
3070 gcc -g -o $@ $< tension_model.c -lcheck -lgsl -lgslcblas -lm
3071 clean_tension_model : clean_tension_model_utils
3072 rm -f tension_model.c tension_model.h check_tension_model.c check_tension_model
3073 tension_model_utils.c : sawsim.nw
3074 notangle -Rtension-model-utils.c $^ > $@
3075 tension_model_utils : tension_model_utils.c global.h tension_model.c tension_model.h parse.c parse.h \
3076 list.c list.h tension_balance.c tension_balance.h
3077 gcc -g -o $@ $< tension_model.c parse.c list.c tension_balance.c -lgsl -lgslcblas -lm
3078 tension_model_utils_static : tension_model_utils.c global.h tension_model.c tension_model.h parse.c parse.h \
3079 list.c list.h tension_balance.c tension_balance.h
3080 gcc -g -static -o $@ $< tension_model.c parse.c list.c tension_balance.c -lgsl -lgslcblas -lm
3081 clean_tension_model_utils :
3082 rm -f tension_model_utils.c tension_model_utils
3086 \label{sec.null_tension}
3088 For unstretchable domains.
3090 <<null tension model getopt>>=
3091 {NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3094 \subsection{Constant}
3095 \label{sec.const_tension}
3097 <<constant tension functions>>=
3098 <<constant tension function>>
3099 <<constant tension parameter create/destroy functions>>
3102 <<constant tension model declarations>>=
3103 <<constant tension function declaration>>
3104 <<constant tension parameter create/destroy function declarations>>
3105 <<constant tension model global declarations>>
3109 An infinitely stretchable domain providing a constant tension.
3111 <<constant tension function declaration>>=
3112 extern double const_tension_handler(double x, void *pdata);
3114 <<constant tension function>>=
3115 double const_tension_handler(double x, void *pdata)
3117 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3118 list_t *list = data->group;
3123 assert(list != NULL); /* empty active group?! */
3124 F = ((const_tension_param_t *)list->d)->F;
3125 while (list != NULL) {
3126 assert(((const_tension_param_t *)list->d)->F == F);
3134 <<constant tension structure>>=
3135 typedef struct const_tension_param_struct {
3136 double F; /* tension (force) in N */
3137 } const_tension_param_t;
3141 <<constant tension parameter create/destroy function declarations>>=
3142 extern void *string_create_const_tension_param_t(char **param_strings);
3143 extern void destroy_const_tension_param_t(void *p);
3146 <<constant tension parameter create/destroy functions>>=
3147 const_tension_param_t *create_const_tension_param_t(double F)
3149 const_tension_param_t *ret
3150 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3151 assert(ret != NULL);
3156 void *string_create_const_tension_param_t(char **param_strings)
3158 return create_const_tension_param_t(atof(param_strings[0]));
3161 void destroy_const_tension_param_t(void *p)
3169 <<constant tension model global declarations>>=
3170 extern int num_const_tension_params;
3171 extern char *const_tension_param_descriptions[];
3172 extern char const_tension_param_string[];
3175 <<constant tension model globals>>=
3176 int num_const_tension_params = 1;
3177 char *const_tension_param_descriptions[] = {"tension F, N"};
3178 char const_tension_param_string[] = "0";
3181 <<constant tension model getopt>>=
3182 {&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}
3188 <<hooke functions>>=
3190 <<hooke parameter create/destroy functions>>
3193 <<hooke tension model declarations>>=
3194 <<hooke tension function declaration>>
3195 <<hooke tension parameter create/destroy function declarations>>
3196 <<hooke tension model global declarations>>
3200 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
3201 The behavior of a series of springs $k_i$ in series is given by
3203 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
3205 For a simple proof, see Appendix \ref{app.math_hooke}.
3207 <<hooke tension function declaration>>=
3208 extern double hooke_handler(double x, void *pdata);
3212 double hooke_handler(double x, void *pdata)
3214 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3215 list_t *list = data->group;
3220 assert(list != NULL); /* empty active group?! */
3221 while (list != NULL) {
3222 assert( ((hooke_param_t *)list->d)->k > 0 );
3223 k += 1.0/ ((hooke_param_t *)list->d)->k;
3231 <<hooke structure>>=
3232 typedef struct hooke_param_struct {
3233 double k; /* spring constant in N/m */
3237 <<hooke tension parameter create/destroy function declarations>>=
3238 extern void *string_create_hooke_param_t(char **param_strings);
3239 extern void destroy_hooke_param_t(void *p);
3242 <<hooke parameter create/destroy functions>>=
3243 hooke_param_t *create_hooke_param_t(double k)
3245 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
3246 assert(ret != NULL);
3251 void *string_create_hooke_param_t(char **param_strings)
3253 return create_hooke_param_t(atof(param_strings[0]));
3256 void destroy_hooke_param_t(void *p)
3263 <<hooke tension model global declarations>>=
3264 extern int num_hooke_params;
3265 extern char *hooke_param_descriptions[];
3266 extern char hooke_param_string[];
3269 <<hooke tension model globals>>=
3270 int num_hooke_params = 1;
3271 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
3272 char hooke_param_string[]="0.05";
3275 <<hooke tension model getopt>>=
3276 {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}
3279 \subsection{Worm-like chain}
3282 We can model several unfolded domains with a single worm-like chain.
3283 <<worm-like chain functions>>=
3284 <<worm-like chain function>>
3285 <<worm-like chain handler>>
3286 <<worm-like chain create/destroy functions>>
3289 <<worm-like chain tension model declarations>>=
3290 <<worm-like chain handler declaration>>
3291 <<worm-like chain create/destroy function declarations>>
3292 <<worm-like chain tension model global declarations>>
3296 The combination of all unfolded domains is modeled as a worm like chain,
3297 with the force $F_{\text{WLC}}$ approximately given by
3299 F_{\text{WLC}} \approx \frac{k_B T}{p}
3300 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
3302 where $x$ is the distance between the fixed ends,
3303 $k_B$ is Boltzmann's constant,
3304 $T$ is the absolute temperature,
3305 $p$ is the persistence length, and
3306 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
3307 <<worm-like chain function>>=
3308 double wlc(double x, double T, double p, double L)
3310 double xL=0.0; /* uses global k_B */
3312 assert(T > 0); assert(p > 0); assert(L > 0);
3313 if (x >= L) return HUGE_VAL;
3314 xL = x/L; /* unitless */
3315 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
3316 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
3317 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
3320 This model assumes that each unfolded domain has the same persistence length.
3322 <<worm-like chain handler declaration>>=
3323 extern double wlc_handler(double x, void *pdata);
3326 <<worm-like chain handler>>=
3327 double wlc_handler(double x, void *pdata)
3329 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3330 list_t *list = data->group;
3334 assert(list != NULL); /* empty active group?! */
3335 p = ((wlc_param_t *)list->d)->p;
3336 while (list != NULL) {
3337 /* enforce consistent p */
3338 assert( ((wlc_param_t *)list->d)->p == p);
3339 L += ((wlc_param_t *)list->d)->L;
3342 return wlc(x, data->env->T, p, L);
3346 <<worm-like chain structure>>=
3347 typedef struct wlc_param_struct {
3348 double p; /* WLC persistence length (m) */
3349 double L; /* WLC contour length (m) */
3353 <<worm-like chain create/destroy function declarations>>=
3354 extern void *string_create_wlc_param_t(char **param_strings);
3355 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
3358 <<worm-like chain create/destroy functions>>=
3359 wlc_param_t *create_wlc_param_t(double p, double L)
3361 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
3362 assert(ret != NULL);
3368 void *string_create_wlc_param_t(char **param_strings)
3370 return create_wlc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3373 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
3380 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.
3381 TODO (now we copy the string before parsing, but still do this for clarity).
3383 <<valid string write code>>=
3384 char string[] = "abc";
3387 <<valid string write code>>=
3388 char *string = "abc";
3392 <<worm-like chain tension model global declarations>>=
3393 extern int num_wlc_params;
3394 extern char *wlc_param_descriptions[];
3395 extern char wlc_param_string[];
3397 <<worm-like chain tension model globals>>=
3398 int num_wlc_params = 2;
3399 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
3400 char wlc_param_string[]="0.39e-9,28.4e-9";
3403 <<worm-like chain tension model getopt>>=
3404 {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}
3406 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
3408 \subsection{Freely-jointed chain}
3411 <<freely-jointed chain functions>>=
3412 <<freely-jointed chain function>>
3413 <<freely-jointed chain handler>>
3414 <<freely-jointed chain create/destroy functions>>
3417 <<freely-jointed chain tension model declarations>>=
3418 <<freely-jointed chain handler declaration>>
3419 <<freely-jointed chain create/destroy function declarations>>
3420 <<freely-jointed chain tension model global declarations>>
3424 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
3425 The tension of a chain of $N$ such links is given by
3427 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
3429 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}.
3430 <<freely-jointed chain function>>=
3431 double langevin(double x, void *params)
3434 return 1.0/tanh(x) - 1/x;
3436 <<one dimensional bracket declaration>>
3437 <<one dimensional solve declaration>>
3438 double inv_langevin(double x)
3440 double lb=0.0, ub=1.0, ret;
3441 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
3442 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
3446 double fjc(double x, double T, double l, int N)
3448 double L = l*(double)N;
3449 if (x == 0) return 0;
3450 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
3451 return k_B*T/l * inv_langevin(x/L);
3455 <<freely-jointed chain handler declaration>>=
3456 extern double fjc_handler(double x, void *pdata);
3458 <<freely-jointed chain handler>>=
3459 double fjc_handler(double x, void *pdata)
3461 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3462 list_t *list = data->group;
3467 assert(list != NULL); /* empty active group?! */
3468 l = ((fjc_param_t *)list->d)->l;
3469 while (list != NULL) {
3470 /* enforce consistent l */
3471 assert( ((fjc_param_t *)list->d)->l == l);
3472 N += ((fjc_param_t *)list->d)->N;
3475 return fjc(x, data->env->T, l, N);
3479 <<freely-jointed chain structure>>=
3480 typedef struct fjc_param_struct {
3481 double l; /* FJC link length (m) */
3482 int N; /* FJC number of links */
3486 <<freely-jointed chain create/destroy function declarations>>=
3487 extern void *string_create_fjc_param_t(char **param_strings);
3488 extern void destroy_fjc_param_t(void *p);
3491 <<freely-jointed chain create/destroy functions>>=
3492 fjc_param_t *create_fjc_param_t(double l, double N)
3494 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
3495 assert(ret != NULL);
3501 void *string_create_fjc_param_t(char **param_strings)
3503 return create_fjc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3506 void destroy_fjc_param_t(void *p)
3513 <<freely-jointed chain tension model global declarations>>=
3514 extern int num_fjc_params;
3515 extern char *fjc_param_descriptions[];
3516 extern char fjc_param_string[];
3519 <<freely-jointed chain tension model globals>>=
3520 int num_fjc_params = 2;
3521 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
3522 char fjc_param_string[]="0.5e-9,200";
3525 <<freely-jointed chain tension model getopt>>=
3526 {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}
3529 \subsection{Tension model implementation}
3531 <<tension-model.c>>=
3533 <<tension model includes>>
3534 #include "tension_model.h"
3535 <<tension model internal definitions>>
3536 <<tension model globals>>
3537 <<tension model functions>>
3540 <<tension model includes>>=
3541 #include <assert.h> /* assert() */
3542 #include <stdlib.h> /* NULL */
3543 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
3544 #include <stdio.h> /* fprintf(), stdout */
3547 #include "tension_balance.h" /* oneD_*() */
3550 <<tension model internal definitions>>=
3551 <<constant tension structure>>
3553 <<worm-like chain structure>>
3554 <<freely-jointed chain structure>>
3557 <<tension model globals>>=
3558 <<tension handler array global>>
3559 <<constant tension model globals>>
3560 <<hooke tension model globals>>
3561 <<worm-like chain tension model globals>>
3562 <<freely-jointed chain tension model globals>>
3565 <<tension model functions>>=
3566 <<constant tension functions>>
3568 <<worm-like chain functions>>
3569 <<freely-jointed chain functions>>
3573 \subsection{Utilities}
3575 The tension models can get complicated, and users may want to reassure
3576 themselves that this portion of the input physics are functioning properly. The
3577 stand-alone program [[tension_model_utils]] demonstrates the tension model
3578 interface without the overhead of having to understand a full simulation.
3579 [[tension_model_utils]] takes command line model arguments like the full
3580 simulation, and prints $F(x)$ data to the screen for the selected model over a
3583 <<tension-model-utils.c>>=
3585 <<tension model utility includes>>
3586 <<tension model utility definitions>>
3587 <<tension model utility getopt functions>>
3589 <<tension model utility main>>
3592 <<tension model utility main>>=
3593 int main(int argc, char **argv)
3595 <<tension model getopt array definition>>
3596 tension_model_getopt_t *model = NULL;
3600 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
3601 tension_handler_data_t tdata;
3602 int num_param_args; /* for INIT_MODEL() */
3603 char **param_args; /* for INIT_MODEL() */
3604 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model, &flags);
3605 setup(tension_handler);
3606 if (flags & VFLAG) {
3607 printf("#initializing model %s with parameters %s\n", model->name, model->params);
3609 INIT_MODEL("utility", model, params);
3611 push(&tdata.group, params, 1);
3613 tdata.persist = NULL;
3614 if (tension_handler[model->tg_group] == NULL) {
3615 printf("No tension function for model %s\n", model->name);
3619 double dx=1e-10, x=0, F=0;
3620 printf("#F (N)\tk (%% pop. per s)\n");
3621 while (F >= 0 && F < 1e5 && x < 1e-6) {
3622 F = (*tension_handler[model->tg_group])(x, &tdata);
3623 printf("%g\t%g\n", x, F);
3627 params = pop(&tdata.group);
3628 if (model->destructor)
3629 (*model->destructor)(params);
3634 <<tension model utility includes>>=
3637 #include <string.h> /* strlen() */
3638 #include <assert.h> /* assert() */
3642 #include "tension_model.h"
3645 <<tension model utility definitions>>=
3646 <<version definition>>
3647 #define VFLAG 1 /* verbose */
3648 <<string comparison definition and globals>>
3649 <<tension model getopt definitions>>
3650 <<initialize model definition>>
3654 <<tension model utility getopt functions>>=
3655 <<index tension model>>
3656 <<tension model utility help>>
3657 <<tension model utility get options>>
3660 <<tension model utility help>>=
3661 void help(char *prog_name,
3663 int n_tension_models, tension_model_getopt_t *tension_models,
3667 printf("usage: %s [options]\n", prog_name);
3668 printf("Version %s\n\n", VERSION);
3669 printf("A utility for understanding the available tension models\n\n");
3670 printf("Environmental options:\n");
3671 printf("-T T\tTemperature T (currently %g K)\n", env->T);
3672 printf("-C T\tYou can also set the temperature T in Celsius\n");
3673 printf("Model options:\n");
3674 printf("-m model\tFolded tension model (currently %s)\n",
3675 tension_models[tension_model].name);
3676 printf("-a args\tFolded tension model argument string (currently %s)\n",
3677 tension_models[tension_model].params);
3678 printf("F(x) is calculated for a range of x and printed\n");
3679 printf("For example:\n");
3680 printf(" #Distance (x)\tForce (N)\n");
3681 printf(" 123.456\t7.89\n");
3683 printf("-V\tChange output to verbose mode\n");
3684 printf("-h\tPrint this help and exit\n");
3686 printf("Tension models:\n");
3687 for (i=0; i<n_tension_models; i++) {
3688 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
3689 for (j=0; j < tension_models[i].num_params ; j++)
3690 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
3691 printf("\t\tdefaults: %s\n", tension_models[i].params);
3696 <<tension model utility get options>>=
3697 void get_options(int argc, char **argv, environment_t *env,
3698 int n_tension_models, tension_model_getopt_t *tension_models,
3699 tension_model_getopt_t **model,
3700 unsigned int *flags)
3702 char *prog_name = NULL;
3703 char c, options[] = "T:C:m:a:Vh";
3704 int tension_model=0, num_strings;
3705 extern char *optarg;
3706 extern int optind, optopt, opterr;
3710 /* setup defaults */
3712 prog_name = argv[0];
3713 env->T = 300.0; /* K */
3715 *model = tension_models;
3717 while ((c=getopt(argc, argv, options)) != -1) {
3719 case 'T': env->T = atof(optarg); break;
3720 case 'C': env->T = atof(optarg)+273.15; break;
3722 parse_list_string(optarg, ',', NULL, NULL, &num_strings, string_array);
3723 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
3724 *model = tension_models+tension_model;
3727 tension_models[tension_model].params = optarg;
3729 case 'V': *flags |= VFLAG; break;
3731 fprintf(stderr, "unrecognized option '%c'\n", optopt);
3732 /* fall through to default case */
3734 help(prog_name, env, n_tension_models, tension_models, tension_model);
3743 \section{Unfolding rate models}
3744 \label{sec.k_models}
3746 We have two main choices for determining $k$: Bell's model, which assumes the
3747 folded domains exist in a single `folded' state and make instantaneous,
3748 stochastic transitions over a free energy barrier; and Kramers' model, which
3749 treats the unfolding as a thermalized diffusion process.
3750 We deal with these options like we dealt with the different tension models: by bundling all model-specific
3751 parameters into structures, with model specific functions for determining $k$.
3753 <<k func definition>>=
3754 typedef double k_func_t(double F, environment_t *env, void *params);
3757 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3758 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]].
3760 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
3761 because \LaTeX\ doesn't like underscores outside math mode.
3762 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
3763 You could use spaces instead (`\verb+ +'),
3764 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
3765 which I don't like as much.
3769 <<k func definition>>
3770 <<null k declarations>>
3771 <<const k declarations>>
3772 <<bell k declarations>>
3773 <<kramers k declarations>>
3774 <<kramers integ k declarations>>
3777 <<k model module makefile lines>>=
3778 k_model.c : sawsim.nw
3779 notangle -Rk-model.c $^ > $@
3780 k_model.h : sawsim.nw
3781 notangle -Rk-model.h $^ > $@
3782 check_k_model.c : sawsim.nw
3783 notangle -Rcheck-k-model.c $^ > $@
3784 check_k_model : check_k_model.c global.h k_model.c k_model.h
3785 gcc -g -o $@ $< k_model.c -lcheck -lgsl -lgslcblas -lm
3786 clean_k_model : clean_k_model_utils
3787 rm -f k_model.c k_model.h check_k_model.c check_k_model
3788 k_model_utils.c : sawsim.nw
3789 notangle -Rk-model-utils.c $^ > $@
3790 k_model_utils : k_model_utils.c global.h k_model.c k_model.h parse.c parse.h string_eval.c string_eval.h
3791 gcc -g -o $@ $< k_model.c parse.c string_eval.c -lgsl -lgslcblas -lm
3792 k_model_utils_static : k_model_utils.c global.h k_model.c k_model.h parse.c parse.h string_eval.c string_eval.h
3793 gcc -g -static -o $@ $< k_model.c parse.c string_eval.c -lgsl -lgslcblas -lm
3794 clean_k_model_utils :
3795 rm -f k_model_utils.c k_model_utils
3801 For domains that are never folded.
3803 <<null k declarations>>=
3805 <<null k functions>>=
3810 <<null k model getopt>>=
3811 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
3814 \subsection{Constant rate model}
3817 This is a pretty boring model, but it is usefull for testing the engine.
3821 where $k_0$ is the constant unfolding rate.
3823 <<const k functions>>=
3824 <<const k function>>
3825 <<const k structure create/destroy functions>>
3828 <<const k declarations>>=
3829 <<const k function declaration>>
3830 <<const k structure create/destroy function declarations>>
3831 <<const k global declarations>>
3835 <<const k structure>>=
3836 typedef struct const_k_param_struct {
3837 double knot; /* transition rate k_0 (frac population per s) */
3841 <<const k function declaration>>=
3842 double const_k(double F, environment_t *env, void *const_k_params);
3844 <<const k function>>=
3845 double const_k(double F, environment_t *env, void *const_k_params)
3846 { /* Returns the rate constant k in frac pop per s. */
3847 const_k_param_t *p = (const_k_param_t *) const_k_params;
3849 assert(p->knot > 0);
3854 <<const k structure create/destroy function declarations>>=
3855 void *string_create_const_k_param_t(char **param_strings);
3856 void destroy_const_k_param_t(void *p);
3859 <<const k structure create/destroy functions>>=
3860 const_k_param_t *create_const_k_param_t(double knot)
3862 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
3863 assert(ret != NULL);
3868 void *string_create_const_k_param_t(char **param_strings)
3870 return create_const_k_param_t(atof(param_strings[0]));
3873 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
3880 <<const k global declarations>>=
3881 extern int num_const_k_params;
3882 extern char *const_k_param_descriptions[];
3883 extern char const_k_param_string[];
3886 <<const k globals>>=
3887 int num_const_k_params = 1;
3888 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
3889 char const_k_param_string[]="1";
3892 <<const k model getopt>>=
3893 {"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}
3896 \subsection{Bell's model}
3899 Quantitatively, Bell's model gives $k$ as
3901 k = k_0 \cdot e^{F dx / k_B T} \;,
3903 where $k_0$ is the unforced unfolding rate,
3904 $F$ is the applied unfolding force,
3905 $dx$ is the distance to the transition state, and
3906 $k_B T$ is the thermal energy\citep{schlierf06}.
3908 <<bell k functions>>=
3910 <<bell k structure create/destroy functions>>
3913 <<bell k declarations>>=
3914 <<bell k function declaration>>
3915 <<bell k structure create/destroy function declarations>>
3916 <<bell k global declarations>>
3920 <<bell k structure>>=
3921 typedef struct bell_param_struct {
3922 double knot; /* transition rate k_0 (frac population per s) */
3923 double dx; /* `distance to transition state' (nm) */
3927 <<bell k function declaration>>=
3928 double bell_k(double F, environment_t *env, void *bell_params);
3930 <<bell k function>>=
3931 double bell_k(double F, environment_t *env, void *bell_params)
3932 { /* Returns the rate constant k in frac pop per s.
3933 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
3934 * uses global k_B in J/K */
3935 bell_param_t *p = (bell_param_t *) bell_params;
3936 assert(F >= 0); assert(env->T > 0);
3938 assert(p->knot > 0); assert(p->dx >= 0);
3940 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
3941 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
3942 p->knot * exp(F*p->dx / (k_B*env->T) ));
3943 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
3945 return p->knot * exp(F*p->dx / (k_B*env->T) );
3949 <<bell k structure create/destroy function declarations>>=
3950 void *string_create_bell_param_t(char **param_strings);
3951 void destroy_bell_param_t(void *p);
3954 <<bell k structure create/destroy functions>>=
3955 bell_param_t *create_bell_param_t(double knot, double dx)
3957 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
3958 assert(ret != NULL);
3964 void *string_create_bell_param_t(char **param_strings)
3966 return create_bell_param_t(atof(param_strings[0]), atof(param_strings[1]));
3969 void destroy_bell_param_t(void *p /* bell_param_t * */)
3976 <<bell k global declarations>>=
3977 extern int num_bell_params;
3978 extern char *bell_param_descriptions[];
3979 extern char bell_param_string[];
3983 int num_bell_params = 2;
3984 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
3985 char bell_param_string[]="3.3e-4,0.25e-9";
3988 <<bell k model getopt>>=
3989 {"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}
3991 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
3992 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
3994 \subsection{Kramer's model}
3997 Kramer's model gives $k$ as
3999 % \frac{1}{k} = \frac{1}{D}
4000 % \int_{x_\text{min}}^{x_\text{max}}
4001 % e^{\frac{-U_F(x)}{k_B T}}
4002 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4005 %where $D$ is the diffusion constant,
4006 %$U_F(x)$ is the free energy along the unfolding coordinate, and
4007 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4008 % before the minimum and after the maximum, respectively \citep{schlierf06}.
4010 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
4012 where $D$ is the diffusion constant,
4014 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
4016 are length scales where
4017 $x_c(F)$ is the location of the bound state and
4018 $x_{ts}(F)$ is the location of the transition state,
4019 $E(F,x)$ is the free energy, and
4020 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
4021 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
4022 \citet{evans97} Eqn.~3).
4024 <<kramers k functions>>=
4025 <<kramers k function>>
4026 <<kramers k structure create/destroy functions>>
4029 <<kramers k declarations>>=
4030 <<kramers k function declaration>>
4031 <<kramers k structure create/destroy function declarations>>
4032 <<kramers k global declarations>>
4036 <<kramers k structure>>=
4037 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
4038 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
4040 typedef struct kramers_param_struct {
4041 double D; /* diffusion rate (um/s) */
4048 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
4049 //kramers_x_func_t *xts; /* function returning transition x (nm) */
4050 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
4051 //kramers_E_func_t *E; /* function returning E (J) */
4052 //double *E_params; /* parametrize the energy landscape */
4053 //int n_E_params; /* length of E_params array */
4057 <<kramers k function declaration>>=
4058 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
4059 extern double kramers_k(double F, environment_t *env, void *kramers_params);
4061 <<kramers k function>>=
4062 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
4064 static char input[160]={0};
4065 static char output[80]={0};
4066 /* setup the environment */
4067 fprintf(in, "F = %g; T = %g\n", F, T);
4068 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4069 string_eval(in, out, input, 80, output);
4070 return atof(output);
4073 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
4075 static char input[160]={0};
4076 static char output[80]={0};
4077 /* setup the environment */
4078 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
4079 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4080 string_eval(in, out, input, 80, output);
4081 return atof(output);
4084 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
4086 kramers_param_t *p = (kramers_param_t *) kramers_params;
4087 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
4090 double kramers_k(double F, environment_t *env, void *kramers_params)
4091 { /* Returns the rate constant k in frac pop per s.
4092 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4093 * uses global k_B in J/K */
4094 kramers_param_t *p = (kramers_param_t *) kramers_params;
4095 double kbT, xc, xts, lc, lts, Eb;
4096 assert(F >= 0); assert(env->T > 0);
4099 assert(p->xc != NULL); assert(p->xts != NULL);
4100 assert(p->ddEdxx != NULL); assert(p->E != NULL);
4101 assert(p->in != NULL); assert(p->out != NULL);
4103 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
4104 if (xc == -1.0) return -1.0; /* error (Too much force) */
4105 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
4106 if (xc == -1.0) return -1.0; /* error (Too much force) */
4107 lc = sqrt(2.0*M_PI*kbT /
4108 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
4109 lts = sqrt(-2.0*M_PI*kbT/
4110 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
4111 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
4112 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
4113 //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));
4114 return p->D/(lc*lts) * exp(-Eb/kbT);
4118 <<kramers k structure create/destroy function declarations>>=
4119 void *string_create_kramers_param_t(char **param_strings);
4120 void destroy_kramers_param_t(void *p);
4123 <<kramers k structure create/destroy functions>>=
4124 kramers_param_t *create_kramers_param_t(double D,
4125 char *xc_fn, char *xts_fn,
4126 char *E_fn, char *ddEdxx_fn,
4127 char *global_define_string)
4128 // kramers_x_func_t *xc, kramers_x_func_t *xts,
4129 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
4130 // double *E_params, int n_E_params)
4132 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
4133 assert(ret != NULL);
4134 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4135 assert(ret->xc != NULL);
4136 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4137 assert(ret->xts != NULL);
4138 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
4139 assert(ret->E != NULL);
4140 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
4141 assert(ret->ddEdxx != NULL);
4143 strcpy(ret->xc, xc_fn);
4144 strcpy(ret->xts, xts_fn);
4145 strcpy(ret->E, E_fn);
4146 strcpy(ret->ddEdxx, ddEdxx_fn);
4147 //ret->E_params = E_params;
4148 //ret->n_E_params = n_E_params;
4149 string_eval_setup(&ret->in, &ret->out);
4150 fprintf(ret->in, "kB = %g\n", k_B);
4151 fprintf(ret->in, "%s\n", global_define_string);
4155 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
4156 void *string_create_kramers_param_t(char **param_strings)
4158 return create_kramers_param_t(atof(param_strings[0]),
4166 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
4168 kramers_param_t *pk = (kramers_param_t *)p;
4170 string_eval_teardown(&pk->in, &pk->out);
4172 // free(pk->E_params);
4178 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
4179 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.
4180 We follow \citet{shillcock98} and use
4182 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
4183 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
4186 where TODO, variable meanings.%\citep{shillcock98}.
4187 The first and second derivatives of $E(F,E_b,x,x_s)$ are
4189 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\\
4190 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
4192 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
4193 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
4194 The roots are therefor at
4196 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
4197 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
4198 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
4201 <<kramers k global declarations>>=
4202 extern int num_kramers_params;
4203 extern char *kramers_param_descriptions[];
4204 extern char kramers_param_string[];
4207 <<kramers k globals>>=
4208 int num_kramers_params = 6;
4209 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"};
4210 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)}";
4212 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
4213 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
4214 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
4215 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.
4216 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
4217 It works so long as [[val_1]] is not [[false]].
4219 <<kramers k model getopt>>=
4220 {"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}
4223 \subsection{Kramer's model (natural cubic splines)}
4224 \label{sec.kramers_integ}
4226 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
4227 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
4228 \citet{schlierf06} Eqn.~2)
4230 \frac{1}{k} = \frac{1}{D}
4231 \int_{x_\text{min}}^{x_\text{max}}
4232 e^{\frac{U_F(x)}{k_B T}}
4233 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4236 where $D$ is the diffusion constant,
4237 $U_F(x)$ is the free energy along the unfolding coordinate, and
4238 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4239 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
4241 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
4242 interpolating between them with cubic splines.
4243 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
4245 <<kramers integ k functions>>=
4246 <<spline functions>>
4247 <<kramers integ A k functions>>
4248 <<kramers integ B k functions>>
4249 <<kramers integ k function>>
4250 <<kramers integ k structure create/destroy functions>>
4253 <<kramers integ k declarations>>=
4254 <<kramers integ k function declaration>>
4255 <<kramers integ k structure create/destroy function declarations>>
4256 <<kramers integ k global declarations>>
4260 <<kramers integ k structure>>=
4261 typedef double func_t(double x, void *params);
4262 typedef struct kramers_integ_param_struct {
4263 double D; /* diffusion rate (um/s) */
4264 double x_min; /* integration bounds */
4266 func_t *E; /* function returning E (J) */
4267 void *E_params; /* parametrize the energy landscape */
4268 destroy_data_func_t *destroy_E_params;
4270 double F; /* for passing into GSL evaluations */
4272 } kramers_integ_param_t;
4275 <<spline param structure>>=
4276 typedef struct spline_param_struct {
4277 int n_knots; /* natural cubic spline knots for U(x) */
4280 gsl_spline *spline; /* GSL spline parameters */
4281 gsl_interp_accel *acc; /* GSL spline acceleration data */
4285 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
4287 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4289 <<kramers integ A k functions>>=
4290 double kramers_integ_k_integrandA(double x, void *params)
4292 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4293 double U = (*p->E)(x, p->E_params);
4294 double Fx = p->F * x;
4295 double kbT = k_B * p->env->T;
4296 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
4297 return exp(-(U-Fx)/kbT);
4299 double kramers_integ_k_integralA(double x, void *params)
4301 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4303 double result, abserr;
4304 size_t neval; /* for qng */
4305 static gsl_integration_workspace *w = NULL;
4306 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4307 f.function = &kramers_integ_k_integrandA;
4309 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4310 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4311 //fprintf(stderr, "integralA = %g\n", result);
4317 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
4318 \int_{x_\text{min}}^{x_\text{max}}
4319 e^{\frac{U_F(x)}{k_B T}}
4320 \text{Integral}_A(x)
4323 <<kramers integ B k functions>>=
4324 double kramers_integ_k_integrandB(double x, void *params)
4326 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4327 double U = (*p->E)(x, p->E_params);
4328 double Fx = p->F * x;
4329 double kbT = k_B * p->env->T;
4330 double intA = kramers_integ_k_integralA(x, params);
4331 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
4332 return exp((U-Fx)/kbT)*intA;
4334 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
4336 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4338 double result, abserr;
4339 size_t neval; /* for qng */
4340 static gsl_integration_workspace *w = NULL;
4341 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4342 f.function = &kramers_integ_k_integrandB;
4346 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4347 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4348 //fprintf(stderr, "integralB = %g\n", result);
4353 With the integrals taken care of, evaluating $k$ is simply
4355 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
4357 <<kramers integ k function declaration>>=
4358 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
4360 <<kramers integ k function>>=
4361 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
4362 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
4363 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
4364 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
4368 <<kramers integ k structure create/destroy function declarations>>=
4369 void *string_create_kramers_integ_param_t(char **param_strings);
4370 void destroy_kramers_integ_param_t(void *p);
4373 <<kramers integ k structure create/destroy functions>>=
4374 kramers_integ_param_t *create_kramers_integ_param_t(double D,
4375 double xmin, double xmax, func_t *E,
4377 destroy_data_func_t *destroy_E_params)
4379 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
4380 assert(ret != NULL);
4385 ret->E_params = E_params;
4386 ret->destroy_E_params = destroy_E_params;
4390 void *string_create_kramers_integ_param_t(char **param_strings)
4392 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
4393 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
4394 void *E_params = string_create_spline_param_t(param_strings+1);
4395 return create_kramers_integ_param_t(atof(param_strings[0]),
4396 atof(param_strings[2]),
4397 atof(param_strings[3]),
4398 &spline_eval, E_params,
4399 destroy_spline_param_t);
4402 void destroy_kramers_integ_param_t(void *params)
4404 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4407 (*p->destroy_E_params)(p->E_params);
4413 Finally we have the GSL spline wrappers:
4414 <<spline functions>>=
4416 <<create/destroy spline>>
4419 <<spline function>>=
4420 double spline_eval(double x, void *spline_params)
4422 spline_param_t *p = (spline_param_t *)(spline_params);
4423 return gsl_spline_eval(p->spline, x, p->acc);
4427 <<create/destroy spline>>=
4428 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
4430 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
4431 assert(ret != NULL);
4432 ret->n_knots = n_knots;
4435 ret->acc = gsl_interp_accel_alloc();
4436 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
4437 gsl_spline_init(ret->spline, x, y, n_knots);
4441 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
4442 void *string_create_spline_param_t(char **param_strings)
4446 double *x=NULL, *y=NULL;
4447 /* break into ordered pairs */
4448 parse_list_string(param_strings[0], SEP, '(', ')',
4449 &num_knots, &knot_coords);
4450 x = (double *)malloc(sizeof(double)*num_knots);
4452 y = (double *)malloc(sizeof(double)*num_knots);
4454 for (i=0; i<num_knots; i++) {
4455 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
4456 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
4459 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
4461 free_parsed_list(num_knots, knot_coords);
4462 return create_spline_param_t(num_knots, x, y);
4465 void destroy_spline_param_t(void *params)
4467 spline_param_t *p = (spline_param_t *)params;
4470 gsl_spline_free(p->spline);
4472 gsl_interp_accel_free(p->acc);
4482 <<kramers integ k global declarations>>=
4483 extern int num_kramers_integ_params;
4484 extern char *kramers_integ_param_descriptions[];
4485 extern char kramers_integ_param_string[];
4488 <<kramers integ k globals>>=
4489 int num_kramers_integ_params = 4;
4490 char *kramers_integ_param_descriptions[] = {
4491 "diffusion D, m^2 / s",
4492 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
4493 "starting integration bound x_min, m",
4494 "ending integration bound x_max, m"};
4495 //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";
4496 //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";
4497 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";
4498 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
4499 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
4500 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
4503 <<kramers integ k model getopt>>=
4504 {"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}
4506 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4507 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4509 \subsection{Unfolding model implementation}
4513 <<k model includes>>
4514 #include "k_model.h"
4515 <<k model internal definitions>>
4517 <<k model functions>>
4520 <<k model includes>>=
4521 #include <assert.h> /* assert() */
4522 #include <stdlib.h> /* NULL, malloc() */
4523 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4524 #include <stdio.h> /* fprintf(), stdout */
4525 #include <string.h> /* strlen(), strcpy() */
4526 #include <gsl/gsl_integration.h>
4527 #include <gsl/gsl_spline.h>
4532 <<k model internal definitions>>=
4533 <<const k structure>>
4534 <<bell k structure>>
4535 <<kramers k structure>>
4536 <<kramers integ k structure>>
4537 <<spline param structure>>
4540 <<k model globals>>=
4544 <<kramers k globals>>
4545 <<kramers integ k globals>>
4548 <<k model functions>>=
4549 <<null k functions>>
4550 <<const k functions>>
4551 <<bell k functions>>
4552 <<kramers k functions>>
4553 <<kramers integ k functions>>
4556 \subsection{Unfolding model unit tests}
4558 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
4559 <<check-k-model.c>>=
4560 <<k model check includes>>
4561 <<check relative error>>
4563 <<k model test suite>>
4564 <<main check program>>
4567 <<k model check includes>>=
4568 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
4569 #include <stdio.h> /* sprintf() */
4570 #include <assert.h> /* assert() */
4571 #include <math.h> /* exp() */
4574 #include "k_model.h"
4577 <<k model test suite>>=
4580 <<k model suite definition>>
4583 <<k model suite definition>>=
4584 Suite *test_suite (void)
4586 Suite *s = suite_create ("k model");
4587 <<const k test case defs>>
4588 <<bell k test case defs>>
4590 <<const k test case adds>>
4591 <<bell k test case adds>>
4596 \subsubsection{Constant}
4598 <<const k test case defs>>=
4599 TCase *tc_const_k = tcase_create("Constant k");
4602 <<const k test case adds>>=
4603 tcase_add_test(tc_const_k, test_const_k_create_destroy);
4604 tcase_add_test(tc_const_k, test_const_k_over_range);
4605 suite_add_tcase(s, tc_const_k);
4610 START_TEST(test_const_k_create_destroy)
4612 char *knot[] = {"1","2","3","4","5","6"};
4613 char *params[] = {knot[0]};
4616 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4617 params[0] = knot[i];
4618 p = string_create_const_param_t(params);
4619 destroy_const_param_t(p);
4624 START_TEST(test_const_k_over_range)
4626 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
4627 char *knot[] = {"1","2","3","4","5","6"};
4628 char *params[] = {knot[0]};
4631 char param_string[80];
4634 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4635 params[0] = knot[i];
4636 p = string_create_const_param_t(params);
4637 for ( F=Fm, F<FM, F+=dF ) {
4638 fail_unless(const_k(F, &env, p)==atof(knot[i]),
4639 "const_k(%g, %g, &{%s,%s}) = %g != %s",
4640 F, T, knot[i], dx, const_k(F, &env, p), knot[i]);
4642 destroy_const_param_t(p);
4648 \subsubsection{Bell}
4650 <<bell k test case defs>>=
4651 TCase *tc_bell = tcase_create("Bell's k");
4654 <<bell k test case adds>>=
4655 tcase_add_test(tc_bell, test_bell_k_create_destroy);
4656 tcase_add_test(tc_bell, test_bell_k_at_zero);
4657 tcase_add_test(tc_bell, test_bell_k_at_one);
4658 suite_add_tcase(s, tc_bell);
4662 START_TEST(test_bell_k_create_destroy)
4664 char *knot[] = {"1","2","3","4","5","6"};
4666 char *params[] = {knot[0], dx};
4669 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4670 params[0] = knot[i];
4671 p = string_create_bell_param_t(params);
4672 destroy_bell_param_t(p);
4677 START_TEST(test_bell_k_at_zero)
4679 double F=0.0, T=300.0;
4680 char *knot[] = {"1","2","3","4","5","6"};
4682 char *params[] = {knot[0], dx};
4685 char param_string[80];
4688 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4689 params[0] = knot[i];
4690 p = string_create_bell_param_t(params);
4691 fail_unless(bell_k(F, &env, p)==atof(knot[i]),
4692 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
4693 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
4694 destroy_bell_param_t(p);
4699 START_TEST(test_bell_k_at_one)
4702 char *knot[] = {"1","2","3","4","5","6"};
4704 char *params[] = {knot[0], dx};
4705 double F= k_B*T/atof(dx);
4708 char param_string[80];
4711 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4712 params[0] = knot[i];
4713 p = string_create_bell_param_t(params);
4714 CHECK_ERR(1e-6, atof(knot[i])*exp(1.0), bell_k(F, &env, p));
4715 //fail_unless(bell_k(F, &env, p)==atof(knot[i])*exp(1.0),
4716 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
4717 // F, T, knot[i], dx, bell_k(F, &env, p), atof(knot[i])*exp(1.0));
4718 destroy_bell_param_t(p);
4724 <<kramers k tests>>=
4727 <<kramers k test case def>>=
4730 <<kramers k test case add>>=
4733 <<k model function tests>>=
4734 <<check relative error>>
4736 START_TEST(test_something)
4738 double k=5, x=3, last_x=2.0, F;
4739 one_dim_fn_t *handlers[] = {&hooke};
4740 void *data[] = {&k};
4741 double xi[] = {0.0};
4743 int new_active_groups = 1;
4744 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
4745 fail_unless(F = k*x, NULL);
4750 \subsection{Utilities}
4752 The unfolding models can get complicated, and are sometimes hard to explain to others.
4753 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
4754 the overhead of having to understand a full simulation.
4755 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
4756 or special mode, where the operation depends on the specific model selected.
4758 <<k-model-utils.c>>=
4760 <<k model utility includes>>
4761 <<k model utility definitions>>
4762 <<k model utility getopt functions>>
4763 <<k model utility multi model E>>
4764 <<k model utility main>>
4767 <<k model utility main>>=
4768 int main(int argc, char **argv)
4770 <<k model getopt array definition>>
4771 k_model_getopt_t *model = NULL;
4776 int num_param_args; /* for INIT_MODEL() */
4777 char **param_args; /* for INIT_MODEL() */
4779 double special_xmin;
4780 double special_xmax;
4781 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
4782 &Fmax, &special_xmin, &special_xmax, &flags);
4783 if (flags & VFLAG) {
4784 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4786 INIT_MODEL("utility", model, params);
4789 if (model->k == NULL) {
4790 printf("No k function for model %s\n", model->name);
4794 printf("Fmax = %g <= 0\n", Fmax);
4800 printf("#F (N)\tk (%% pop. per s)\n");
4801 for (i=0; i<=N; i++) {
4802 F = Fmax*i/(double)N;
4803 k = (*model->k)(F, &env, params);
4805 printf("%g\t%g\n", F, k);
4810 if (model->k == NULL || model->k == &bell_k) {
4811 printf("No special mode for model %s\n", model->name);
4814 if (special_xmax <= special_xmin) {
4815 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
4819 double Xrng=(special_xmax-special_xmin), x, E;
4821 printf("#x (m)\tE (J)\n");
4822 for (i=0; i<=N; i++) {
4823 x = special_xmin + Xrng*i/(double)N;
4824 E = multi_model_E(model->k, params, &env, x);
4825 printf("%g\t%g\n", x, E);
4832 if (model->destructor)
4833 (*model->destructor)(params);
4838 <<k model utility multi model E>>=
4839 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
4841 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
4843 if (k == kramers_integ_k)
4844 E = (*p->E)(x, p->E_params);
4846 E = kramers_E(0, env, model_params, x);
4852 <<k model utility includes>>=
4855 #include <string.h> /* strlen() */
4856 #include <assert.h> /* assert() */
4859 #include "string_eval.h"
4860 #include "k_model.h"
4863 <<k model utility definitions>>=
4864 <<version definition>>
4865 #define VFLAG 1 /* verbose */
4866 enum mode_t {M_K_OF_F, M_SPECIAL};
4867 <<string comparison definition and globals>>
4868 <<k model getopt definitions>>
4869 <<initialize model definition>>
4870 <<kramers k structure>>
4871 <<kramers integ k structure>>
4875 <<k model utility getopt functions>>=
4877 <<k model utility help>>
4878 <<k model utility get options>>
4881 <<k model utility help>>=
4882 void help(char *prog_name,
4884 int n_k_models, k_model_getopt_t *k_models,
4888 printf("usage: %s [options]\n", prog_name);
4889 printf("Version %s\n\n", VERSION);
4890 printf("A utility for understanding the available unfolding models\n\n");
4891 printf("Environmental options:\n");
4892 printf("-T T\tTemperature T (currently %g K)\n", env->T);
4893 printf("-C T\tYou can also set the temperature T in Celsius\n");
4894 printf("Model options:\n");
4895 printf("-k model\tTransition rate model (currently %s)\n",
4896 k_models[k_model].name);
4897 printf("-K args\tTransition rate model argument string (currently %s)\n",
4898 k_models[k_model].params);
4899 printf("There are two output modes. In standard mode, k(F) is printed\n");
4900 printf("For example:\n");
4901 printf(" #Force (N)\tk (% pop. per s)\n");
4902 printf(" 123.456\t7.89\n");
4904 printf("In special mode, the output depends on the model.\n");
4905 printf("For models defining an energy function E(x), that function is printed\n");
4906 printf(" #Position (m)\tE (J)\n");
4907 printf(" 0.0012\t0.0034\n");
4909 printf("-m\tChange output to standard mode\n");
4910 printf("-M\tChange output to special mode\n");
4911 printf("-F\tSet the maximum F value for the standard mode k(F) output\n");
4912 printf("-x\tSet the minimum x value for the special mode E(x) output\n");
4913 printf("-X\tSet the maximum x value for the special mode E(x) output\n");
4914 printf("-V\tChange output to verbose mode\n");
4915 printf("-h\tPrint this help and exit\n");
4917 printf("Unfolding rate models:\n");
4918 for (i=0; i<n_k_models; i++) {
4919 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
4920 for (j=0; j < k_models[i].num_params ; j++)
4921 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
4922 printf("\t\tdefaults: %s\n", k_models[i].params);
4927 <<k model utility get options>>=
4928 void get_options(int argc, char **argv, environment_t *env,
4929 int n_k_models, k_model_getopt_t *k_models,
4930 enum mode_t *mode, k_model_getopt_t **model,
4931 double *Fmax, double *special_xmin, double *special_xmax,
4932 unsigned int *flags)
4934 char *prog_name = NULL;
4935 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
4937 extern char *optarg;
4938 extern int optind, optopt, opterr;
4942 /* setup defaults */
4944 prog_name = argv[0];
4945 env->T = 300.0; /* K */
4950 *special_xmax = 1e-8;
4951 *special_xmin = 0.0;
4953 while ((c=getopt(argc, argv, options)) != -1) {
4955 case 'T': env->T = atof(optarg); break;
4956 case 'C': env->T = atof(optarg)+273.15; break;
4958 k_model = index_k_model(n_k_models, k_models, optarg);
4959 *model = k_models+k_model;
4962 k_models[k_model].params = optarg;
4964 case 'm': *mode = M_K_OF_F; break;
4965 case 'M': *mode = M_SPECIAL; break;
4966 case 'F': *Fmax = atof(optarg); break;
4967 case 'x': *special_xmin = atof(optarg); break;
4968 case 'X': *special_xmax = atof(optarg); break;
4969 case 'V': *flags |= VFLAG; break;
4971 fprintf(stderr, "unrecognized option '%c'\n", optopt);
4972 /* fall through to default case */
4974 help(prog_name, env, n_k_models, k_models, k_model);
4983 \section{\LaTeX\ documentation}
4985 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
4986 The comment blocks are extracted (with nicely formatted code blocks), using
4987 <<latex makefile lines>>=
4988 sawsim.tex : sawsim.nw
4989 noweave -latex -index -delay $^ > $@
4992 <<latex makefile lines>>=
4993 sawsim.pdf : sawsim.tex sawsim.bib
4999 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
5000 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
5001 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
5003 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
5004 <<latex makefile lines>>=
5006 rm -f sawsim.aux sawsim.bbl sawsim.blg sawsim.log sawsim.out sawsim.tex
5007 clean_latex : semi-clean_latex
5013 [[make]] is a common utility on *nix systems for managing dependencies while building software.
5014 In this case, we have several \emph{targets} that we'd like to build:
5015 the main simulation program \prog;
5016 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
5017 the documentation [[sawsim.pdf]];
5018 and the [[Makefile]] itself.
5019 Besides the generated files, there is the utility target
5020 [[clean]] for removing all generated files except [[Makefile]].
5022 % [[dist]] for generating a distributable tar file.
5024 Extract the makefile with
5025 `[[notangle -Rmakefile sawsim.nw | sed 's/ /\t/' > Makefile]]'.
5026 The sed is needed because notangle replaces tabs with eight spaces,
5027 but [[make]] needs tabs.
5029 all : sawsim sawsim.pdf Makefile tension_model_utils k_model_utils
5033 profile : sawsim_profile
5034 sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ -mnull -Mwlc -A0.39e-9,28e-9 -F8
5035 gprof sawsim_profile gmon.out > $@
5037 clean : clean_check clean_list clean_tension clean_k_model clean_parse clean_string_eval clean_latex
5038 rm -f sawsim.c sawsim sawsim_static sawsim_profile gmon.out global.h
5040 sawsim.c : sawsim.nw
5042 global.h : sawsim.nw
5043 notangle -Rglobal.h $^ > $@
5044 sawsim : sawsim.c global.h list.c list.h tension_balance.c tension_balance.h \
5045 tension_model.c tension_model.h k_model.c k_model.h parse.c parse.h string_eval.c string_eval.h \
5047 gcc -g -o $@ $< list.c tension_balance.c tension_model.c k_model.c parse.c string_eval.c accel_k.c interp.c tavl.c -lgsl -lgslcblas -lm
5048 sawsim_static : sawsim
5049 gcc -g -static -o $@ sawsim.c list.c tension_balance.c tension_model.c k_model.c parse.c string_eval.c accel_k.c interp.c tavl.c -lgsl -lgslcblas -lm
5050 sawsim_profile : sawsim
5051 gcc -g -pg -o $@ sawsim.c list.c tension_balance.c tension_model.c k_model.c parse.c string_eval.c accel_k.c interp.c tavl.c -lgsl -lgslcblas -lm
5054 check_sawsim.c : sawsim.nw
5055 notangle -Rchecks $^ > $@
5056 check_sawsim : check_sawsim.c global.h list.c list.h tension_balance.c tension_balance.h \
5057 k_model.c k_model.h parse.c parse.h string_eval.c string_eval.h accel_k.c accel_k.h
5058 gcc -g -o $@ $< list.c tension_balance.c k_model.c parse.c string_eval.c accel_k.c interp.c tavl.c -lgsl -lgslcblas -lm -lcheck
5060 rm -f check_sawsim check_sawsim.c
5061 check : check_list check_tension_balance check_k_model check_parse check_string_eval check_accel_k check_sawsim
5063 check_tension_balance
5070 <<list module makefile lines>>
5071 <<tension balance module makefile lines>>
5072 <<tension model module makefile lines>>
5073 <<k model module makefile lines>>
5074 <<parse module makefile lines>>
5075 <<string eval module makefile lines>>
5076 <<accel k module makefile lines>>
5077 <<latex makefile lines>>
5079 Makefile : sawsim.nw
5080 notangle -Rmakefile $^ | sed 's/ /\t/' > $@
5082 This is fairly self-explanatory.
5083 We compile a dynamically linked version ([[sawsim]]) and a statically linked version ([[sawsim_static]]).
5084 The static version is about 9 times as big, but it runs without needing dynamic GSL libraries to link to, so it's a better format for distributing to a cluster for parallel evaluation.
5089 \subsection{Hookean springs in series}
5090 \label{app.math_hooke}
5093 The effective spring constant for $n$ springs $k_i$ in series is given by
5095 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5101 F &= k_i x_i &&\forall i \le n \\
5102 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
5103 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
5104 F &= k_1 x_1 = k_\text{eff} x \\
5105 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
5106 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5111 \addcontentsline{toc}{section}{References}
5112 \bibliography{sawsim}