1 %%%%%%%%%%%%%%%%%%%%%%%%% Start LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%
2 % we have our own preamble,
3 % use `noweave -delay` to not print noweb's automatic one
4 % -*- mode: Noweb; noweb-code-mode: c-mode -*-
5 \documentclass[letterpaper, 11pt]{article}
8 \noweboptions{smallcode}
10 \usepackage{url} % needed for natbib for underscores in urls?
11 \usepackage[breaklinks,colorlinks]{hyperref} % for internal links, \phantomsection
12 % breaklinks breaks long links
13 % colorlinks colors link text instead of boxing it
14 \usepackage{amsmath} % for \text{}, \xrightarrow[sub]{super}
15 \usepackage{amsthm} % for theorem and proof environments
16 \newtheorem{theorem}{Theorem}
18 \usepackage{doc} % BibTeX
19 \usepackage[super,sort&compress]{natbib} % fancy citation extensions
20 % super selects citations in superscript mode
21 % sort&compress automatically sorts and compresses compound citations (\cite{a,b,...})
23 %\bibliographystyle{ieeetr} % pick the bibliography style, short and sweet
24 %\bibliographystyle{plain} % pick the bibliography style, includes dates
25 \bibliographystyle{plainnat}
26 \defcitealias{sw:noweb}{{\tt noweb}}
27 \defcitealias{galassi05}{GSL}
28 \defcitealias{sw:check}{{\tt check}}
29 \defcitealias{sw:python}{Python}
34 \textheight 9.5in % leave a bit of extra room for page numbers (at top of page)
37 \setlength{\parindent}{0pt}
38 \setlength{\parskip}{5pt}
40 % For some reason, the \maketitle wants to go on it's own page
41 % so we'll just hardcode the following in our first page.
42 %\title{Sawsim: a sawtooth protein unfolding simulator}
43 %\author{W.~Trevor King}
46 \newcommand{\prog}{[[sawsim]]}
47 \newcommand{\lang}{[[C]]}
48 \newcommand{\p}[3]{\left#1 #2 \right#3} % e.g. \p({complicated expr})
49 \newcommand{\Th}{\ensuremath{\sup{\text{th}}}} % ordinal th
50 \newcommand{\U}[1]{\ensuremath{\text{ #1}}} % units: $1\U{m/s}$
52 % single spaced lists, from Alvin Alexander
53 % http://www.devdaily.com/blog/post/latex/control-line-spacing-in-itemize-enumerate-tags/
54 \newenvironment{packed_item}{
56 \setlength{\itemsep}{1pt}
57 \setlength{\parskip}{0pt}
58 \setlength{\parsep}{0pt}
62 \nwfilename{sawsim.nw}
67 Sawsim: a sawtooth protein unfolding simulator \\
72 %%%%%%%%%%%%%%%%%%%%%%%%%%%% End LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%%
74 \section{Introduction}
76 The unfolding of multi-globular protein strings is a stochastic
77 process. In the \prog\ program, we use Monte Carlo methods to
78 simulate this unfolding through an explicit time-stepping approach.
79 We develop a program in \lang\ to simulate probability of unfolding
80 under a constant extension velocity of a chain of globular domains.
82 \subsection{Related work}
86 \subsection{About this document}
88 This document was written in \citetalias{sw:noweb} with code blocks in
89 \lang\ and comment blocks in \LaTeX.
91 \subsection{Change Log}
93 Version 0.0 used only the unfolded domain WLC to determine the tension
94 and had a fixed timestep $dt$.
96 Version 0.1 added adaptive timesteps, adjusting so that the unfolding
97 probability for a given domain was constant.
99 Version 0.2 added Kramers' model unfolding rate calculations, and a
100 switch to select Bell's or Kramers' model. This induced a major
101 rewrite, introducing the tension group abstraction, and a split to
102 multiple \lang\ source files. Also added a unit-test suites for
103 sanity checking, and switched to SI units throughout.
105 Version 0.3 added integrated Kramers' model unfolding rate
106 calculations. Also replaced some of my hand-coded routines with GNU
107 Scientific Library \citepalias{galassi05} calls.
109 Version 0.4 added Python evaluation for the saddle-point Kramers'
110 model unfolding rate calculations, so that the model functions could
111 be controlled from the command line. Also added interpolating lookup
112 tables to accelerate slow unfolding rate model evaluation, and command
113 line control of tension groups.
115 <<version definition>>=
116 #define VERSION "0.4"
122 sawsim - program for simulating single molecule mechanical unfolding.
123 Copyright (C) 2008, William Trevor King
125 This program is free software; you can redistribute it and/or
126 modify it under the terms of the GNU General Public License as
127 published by the Free Software Foundation; either version 3 of the
128 License, or (at your option) any later version.
130 This program is distributed in the hope that it will be useful, but
131 WITHOUT ANY WARRANTY; without even the implied warranty of
132 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
133 See the GNU General Public License for more details.
135 You should have received a copy of the GNU General Public License
136 along with this program; if not, write to the Free Software
137 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
140 The author may be contacted at <wking@drexel.edu> on the Internet, or
141 write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
142 Philadelphia PA 19104, USA.
154 The unfolding system is stored as a chain of \emph{domains}. Domains
155 can be folded, globular protein domains, unfolded protein linkers, AFM
156 cantilevers, or other stretchable system components.
158 Each domain contributes to the total tension, which depends on the
159 chain extension. The particular model for the tension as a function
160 of extension is handled generally with function pointers. So far the
161 following models have been implemented:
163 \item Null (Appendix \ref{sec.null_tension}),
164 \item Constant (Appendix \ref{sec.const_tension}),
165 \item Hookean spring (Appendix \ref{sec.hooke}),
166 \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
167 \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
170 The tension model and parameters used for a particular domain depend
171 on whether the domain is folded or unfolded. The transition rate from
172 the folded state to the unfolded state is also handled generally with
173 function pointers, with implementations for
175 \item Null (Appendix \ref{sec.null_k}),
176 \item Bell's model (Appendix \ref{sec.bell}),
177 \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
178 \item Kramers' saddle point model (Appendix \ref{sec.kramers}).
181 The unfolding of the chain domains is modeled in two stages. First
182 the tension of the chain is calculated. Then the tension (assumed to
183 be constant over the length of the chain), is applied to each folded
184 domain, and used to calculate the probability of that subdomain
185 unfolding in the next timestep $dt$. Then the time is stepped
186 forward, and the process is repeated.
189 int main(int argc, char **argv)
191 <<tension model getopt array definition>>
192 <<k model getopt array definition>>
193 list_t *domain_list=NULL; /* variables initialized in get_options() */
194 environment_t env={0};
196 int num_folded=0, unfolding;
197 double x=0, dt, F; /* variables used in the simulation loop */
198 get_options(argc, argv, &P, &dt_max, &v, &env, NUM_TENSION_MODELS,
199 tension_models, NUM_K_MODELS, k_models, &domain_list, &num_folded);
201 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Position (m)\tForce (N)\n");
202 if (flags & FLAG_OUTPUT_UNFOLDING_FORCES) printf("#Unfolding Force (N)\n");
203 while (num_folded > 0) {
204 F = find_tension(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, unfolding);
205 dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, v);
206 unfolding=random_unfoldings(domain_list, F, dt, &env, &num_folded);
207 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\n", x, F);
210 destroy_domain_list(domain_list);
214 @ The meat of the simulation is bundled into the three functions
215 [[find_tension]], [[determine_dt]], and [[random_unfoldings]].
216 [[find_tension]] is discussed in Section \ref{sec.find_tension},
217 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
218 [[random_unfoldings]] in Section \ref{sec.unfolding_rate}.
220 Environmental parameters important in determining reaction rates and
221 tensions (e.g. temperature, pH) are stored in a single structure to
222 facilitate extension to more complicated models in the future.
223 <<environment definition>>=
224 typedef struct environment_struct {
230 <<environment definition>>
231 <<one dimensional function definition>>
232 <<create/destroy definitions>>
236 \section{Simulation functions}
238 <<simulation functions>>=
242 <<does domain unfold>>
243 <<randomly unfold domains>>
247 \label{sec.find_tension}
249 Because the stretched system may be made up of several parts (folded
250 domains, unfolded domains, spring-like cantilever, \ldots), we will
251 assign the domains to models and groups. The models are used to
252 determine a domain's tension (Hookean spring, WLC, \ldots). Groups
253 will represent collections of domains which the model should treat as
254 a single entity. A domain may belong to different groups or models
255 depending on its state. For example, a domain might be modeled as a
256 freely-jointed chain when it iss folded and as a worm-like chain class
257 when it is unfolded. The domains are assumed to be commutative, so
258 ordering is ignored. The interactions between the groups are assumed
259 to be linear, but the interactions between domains of the same group
260 need not be. This allows for non-linear group models such as th
261 worm-like or freely-jointed chains. Each model has a tension handler
262 function, which gives the tension $F$ needed to stretch a given group
263 of domains a distance $x$.
265 To understand the purpose of groups, consider a chain of proteins
266 where two folded proteins $A$ and $B$ are modeled as WLCs with
267 persistence length $p_f$ and two unfolded proteins $C$ and $D$ are
268 modeled as WLCs with persistence length $p_u$. The proteins are
269 attached to a cantilever $E$ qof spring constant $κ$. The string
270 would get split into three lists:
272 \begin{tabular}{llll}
273 Model & Group & List & Tension \\
274 WLC & 0 & $[A,B]$ & $F_\text{WLC}(x, p_f, L_A+L_B)$ \\
275 WLC & 1 & $[C,D]$ & $F_\text{WLC}(x, p_u, L_C+L_D)$ \\
276 Hooke & 0 & $[E]$ & $F_\text{Hooke}(x, κ)$ \\
279 Note that group numbers only matter \emph{within} model classes, since
280 grouping domains with seperate models doesn't make sense.
282 <<find tension definitions>>=
283 #define NUM_TENSION_MODELS 5
287 <<tension handler array global declaration>>=
288 extern one_dim_fn_t *tension_handlers[];
291 <<tension handler array global>>=
292 one_dim_fn_t *tension_handlers[] = {
294 &const_tension_handler,
302 <<tension model global declarations>>=
303 <<tension handler array global declaration>>
306 <<tension handler types>>=
307 typedef struct tension_handler_data_struct {
308 /* int num_domains; */
309 /* domain_t *domains; */
313 } tension_handler_data_t;
317 After sorting the chain into separate groups $G_i$, with tension
318 handlers $F_i(G_i; x_i)$, we need to balance the extension $x_i$ of
319 each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\;
320 \forall i,j$. Note that there may be several groups within each model
321 class. [[find_tension]] is basically a wrapper to feed proper input
322 into the [[tension_balance]] function. [[unfolding]] should be set to
323 0 if there were no unfoldings since the previous call to
326 double find_tension(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list, environment_t *env, double x, int unfolding)
328 static int num_active_groups=0;
329 static one_dim_fn_t **per_group_handlers = NULL;
330 static void **per_group_params = NULL;
331 static double last_x;
332 tension_handler_data_t *tdata;
337 fprintf(stderr, "finding tension at distance %g%s\n", x, unfolding ? " after an unfolding" : "");
338 list_print(stderr, domain_list, "domain list");
341 if (unfolding != 0 || num_active_groups == 0) { /* setup statics */
342 /* free old statics */
343 if (per_group_handlers != NULL)
344 free(per_group_handlers);
345 if (per_group_params != NULL) {
346 for (i=0; i < num_active_groups; i++) {
347 assert(per_group_params[i] != NULL);
348 free(per_group_params[i]);
350 free(per_group_params);
353 /* get new statics */
354 get_active_groups(domain_list, num_tension_handlers, tension_handlers, &num_active_groups, &per_group_handlers, &per_group_params);
356 /* fill in the group handlers and params */
357 for (i=0; i<num_active_groups; i++) {
358 tdata = (tension_handler_data_t *) per_group_params[i];
360 fprintf(stderr, "active group %d of %d", i, num_active_groups);
361 list_print(stderr, tdata->group, " parameters");
364 /* tdata->persist continues unchanged */
367 } /* else roll with the current statics */
369 F = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_params, last_x, x);
373 @ For the moment, we will restrict our group boundaries to a single
374 dimension, so $\sum_i x_i = x$, our total extension (it is this
375 restriction that causes the groups to interact linearly). We'll also
376 restrict our extensions to all be positive. With these restrictions,
377 the problem of balancing the tensions reduces to solving a system of
378 $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is
379 the number of active groups. In general this can be fairly
380 complicated, but without much loss of practicality we can restrict
381 ourselves to strictly monotonically increasing, non-negative tension
382 functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much
383 simpler. For example, it guarantees the existence of a unique, real
384 solution for finite forces. See Appendix \ref{app.tension_balance}
385 for the tension balancing implementation.
387 Each group must define a parameter structure if the tension model is
389 a function to create the parameter structure,
390 a function to destroy the parameter structure, and
391 a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as it's data argument.
392 The tension-specific parameter structures are contained in the domain
393 group field. See the Section \ref{app.model_selection} for a
394 discussion on the command line framework. See the worm-like chain
395 implementation for an example (Section \ref{sec.wlc}).
397 The major limitation of our approach is that all of our tension models
398 are Markovian (memory-less), which is not really a problem for the
399 chains (see \citet{evans99} for WLC thermalization rate calculations).
401 \subsection{Unfolding rate}
402 \label{sec.unfolding_rate}
404 Each folded domain is modeled as unimolecular, first order reaction
406 \text{Folded} \xrightarrow{k} % in tex: X atop Y
409 With the rate constant $k$ defined by
411 \frac{d[\text{Folded}]}{dt} = -k [\text{Folded}].
413 So the probability of a given protein unfolding in a short time $dt$
419 <<does domain unfold>>=
420 int domain_unfolds(double F, double dt, environment_t *env, domain_t *domain)
421 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to a folded domain */
423 k = accel_k(domain->k, F, env, domain->k_params);
424 //(*domain->k)(F, env, domain->k_params);
425 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
426 return happens(k*dt); /* dice roll for prob. k*dt event */
428 @ [[happens]] is a random decision making function defined in Appendix \ref{app.utils}.
430 Only one domain can unfold in each timestep, because the timescale of
431 a domain unfolding $dt_u$ is assumed to be less than the simulation
432 timestep $dt$, so a domain will completely unfold in a single
433 timestep. We adapt our timesteps to keep the probability of a single
434 domain unfolding low, so the probability of two domains unfolding in
435 the same timestep is negligible. This approach breaks down as the
436 adaptive timestep scheme approaches $dt \sim dt_u$, but $dt_u \sim
437 1\U{$\mu$s}$ for Ig27-like proteins \citep{klimov00}, so this
438 shouldn't be a problem. To reassure yourself, you can ask the
439 simulator to print the smallest timestep that was used with TODO.
440 <<randomly unfold domains>>=
441 int random_unfoldings(list_t *domain_list, double tension,
442 double dt, environment_t *env,
443 int *num_folded_domains)
444 { /* returns 1 if there was an unfolding and 0 otherwise */
445 while (domain_list != NULL) {
446 if (D(domain_list)->state == DS_FOLDED
447 && domain_unfolds(tension, dt, env, D(domain_list))) {
448 if (flags & FLAG_OUTPUT_UNFOLDING_FORCES)
449 fprintf(stdout, "%g\n", tension);
450 D(domain_list)->state = DS_UNFOLDED;
451 (*num_folded_domains)--;
452 return 1; /* our one domain has unfolded, stop looking for others */
454 domain_list = domain_list->next;
460 \subsection{Adaptive timesteps}
461 \label{sec.adaptive_dt}
463 We'd like to pick $dt$ so the probability of unfolding in the next
464 timestep is small. If we don't adapt our timestep, we risk breaking
465 our approximation that $F(x' \in [x, x+v \cdot dt]) = F(x)$ or that
466 only one domain unfolds in a given timestep. Because $F(x)$ is
467 monotonically increasing, excessively large timesteps will lead to
468 erroneously large unfolding forces. The simulation would have been
469 accurate for sufficiently small fixed timesteps, but adaptive
470 timesteps will allow us to move through low tension regions in fewer
471 steps, leading to a more efficient simulation.
473 The actual adaptive timestep implementation is not particularly
474 interesting, since we are only required to reduce $dt$ to somewhere
475 below a set threshold, so I've removed it to Appendix
476 \ref{app.adaptive_dt}.
482 The models provide the physics of the simulation, but the simulation
483 allows interchangeable models, and we are currently focusing on the
484 simulation itself, so we remove the actual model implementation to the
487 The tension models are defined in Appendix \ref{sec.tension_models}
488 and unfolding models are defined in Appendix \ref{sec.k_models}.
491 #define k_B 1.3806503e-23 /* J/K */
495 \section{Command line interface}
497 <<get option functions>>=
499 <<index tension model>>
505 \subsection{Model selection}
506 \label{app.model_selection}
508 The main difficulty with the command line interface in \prog\ is
509 developing an intuitive method for accessing the possibly large number
510 of available models. We'll treat this generally by defining an array
511 of available models, containing their important parameters.
513 <<get options definitions>>=
514 <<model getopt definitions>>
517 <<create/destroy definitions>>=
518 typedef void *create_data_func_t(char **param_strings);
519 typedef void destroy_data_func_t(void *param_struct);
522 <<model getopt definitions>>=
523 <<tension model getopt definitions>>
524 <<k model getopt definitions>>
528 \subsubsection{Tension}
530 To access the various tension models from the command line, we define
531 a structure that contains all the neccessary information about the
533 <<tension model getopt definitions>>=
534 typedef struct tension_model_getopt_struct {
535 one_dim_fn_t *handler; /* fn implementing the model on a group */
536 char *name; /* name string identifying the model */
537 char *description; /* a brief description of the model */
538 int num_params; /* number of extra params for the model */
539 char **param_descriptions; /* descriptions of extra params */
540 char *params; /* default values for the extra params */
541 create_data_func_t *creator; /* param-string -> model-data-struct fn */
542 destroy_data_func_t *destructor; /* fn to free the model data structure */
543 } tension_model_getopt_t;
546 <<tension model getopt array definition>>=
547 tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
548 <<null tension model getopt>>,
549 <<constant tension model getopt>>,
550 <<hooke tension model getopt>>,
551 <<worm-like chain tension model getopt>>,
552 <<freely-jointed chain tension model getopt>>
556 \subsubsection{Unfolding rate}
558 <<k model getopt definitions>>=
559 #define NUM_K_MODELS 5
561 typedef struct k_model_getopt_struct {
566 char **param_descriptions;
568 create_data_func_t *creator;
569 destroy_data_func_t *destructor;
573 <<k model getopt array definition>>=
574 k_model_getopt_t k_models[NUM_K_MODELS] = {
575 <<null k model getopt>>,
576 <<const k model getopt>>,
577 <<bell k model getopt>>,
578 <<kramers k model getopt>>,
579 <<kramers integ k model getopt>>
586 void help(char *prog_name, double P, double dt_max, double v,
588 int n_tension_models, tension_model_getopt_t *tension_models,
589 int folded_tension_model, int unfolded_tension_model,
590 int n_k_models, k_model_getopt_t *k_models,
594 printf("usage: %s [options]\n", prog_name);
595 printf("Version %s\n\n", VERSION);
596 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
597 printf("Simulation options:\n");
598 printf("-P P\tTarget probability for dt (currently %g)\n", P);
599 printf("-t dt\tMaximum allowed timestep dt (currently %g)\n", dt_max);
600 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
601 printf("Environmental options:\n");
602 printf("-T T\tTemperature T (currently %g K)\n", env->T);
603 printf("-C T\tYou can also set the temperature T in Celsius\n");
604 printf("Model options:\n");
605 printf("The domains exist in either the folded or unfolded state\n");
606 printf("The following options change the default behavior in each state.\n");
607 printf("-m model[,group]\tFolded tension model (currently %s)\n",
608 tension_models[folded_tension_model].name);
609 printf("-a args\tFolded tension model argument string (currently %s)\n",
610 tension_models[folded_tension_model].params);
611 printf("-M model[,group]\tUnfolded tension model (currently %s)\n",
612 tension_models[unfolded_tension_model].name);
613 printf("-A args\tUnfolded tension model argument string (currently %s)\n",
614 tension_models[unfolded_tension_model].params);
615 printf("The following options change the unfolding rate.\n");
616 printf("-k model\tTransition rate model (currently %s)\n",
617 k_models[k_model].name);
618 printf("-K args\tTransition rate model argument string (currently %s)\n",
619 k_models[k_model].params);
620 printf("Domain creation options:\n");
621 printf("Once you've set up the appropriate default models, you need to add the domains\n");
622 printf("-F n\tAdd n folded domains with the current models\n");
623 printf("-U n\tAdd n unfolded domains with the current models\n");
624 printf("Output mode options:\n");
625 printf("There are two output modes. In standard mode, only the unfolding\n");
626 printf("events are printed. For example:\n");
627 printf(" #Unfolding Force (N)\n");
628 printf(" 123.456e-12\n");
630 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
631 printf(" #Position (m)\tForce (N)\n");
632 printf(" 0.001\t0.0023\n");
634 printf("-V\tChange output to verbose mode\n");
635 printf("-h\tPrint this help and exit\n");
637 printf("Tension models:\n");
638 for (i=0; i<n_tension_models; i++) {
639 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
640 for (j=0; j < tension_models[i].num_params ; j++)
641 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
642 printf("\t\tdefaults: %s\n", tension_models[i].params);
644 printf("Unfolding rate models:\n");
645 for (i=0; i<n_k_models; i++) {
646 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
647 for (j=0; j < k_models[i].num_params ; j++)
648 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
649 printf("\t\tdefaults: %s\n", k_models[i].params);
651 printf("Examples:\n");
652 printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded\n");
653 printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8\n", prog_name);
654 printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded\n");
655 printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K5e-5,.225e-9 -mwlc -a1e-8,4e-10 -Mwlc,1 -A0.4e-9,28.1e-9 -F8\n", prog_name);
659 \subsection{Get options}
662 void get_options(int argc, char **argv,
663 double *pP, double *pDt_max, double *pV,
665 int n_tension_models, tension_model_getopt_t *tension_models,
666 int n_k_models, k_model_getopt_t *k_models,
667 list_t **pDomain_list, int *pNum_folded)
669 char *prog_name = NULL;
670 char c, options[] = "P:t:v:T:C:m:a:M:A:k:K:F:U:Vh";
671 int ftension_model=0, utension_model=0, k_model=0;
672 char *ftension_params=NULL, *utension_params=NULL;
673 int i, n, ftension_group=0, utension_group=0;
677 extern int optind, optopt, opterr;
682 for (i=0; i<n_tension_models; i++) {
683 fprintf(stderr, "tension model %i %s:\t model's handler %p\n",
684 i, tension_models[i].name, tension_models[i].handler);
685 assert(tension_models[i].handler == tension_handlers[i]);
690 flags = FLAG_OUTPUT_UNFOLDING_FORCES;
692 *pP = 1e-3; /* % pop per s */
693 *pDt_max = 0.001; /* s */
694 *pV = 1e-6; /* m/s */
695 env->T = 300.0; /* K */
697 while ((c=getopt(argc, argv, options)) != -1) {
699 case 'P': *pP = atof(optarg); break;
700 case 't': *pDt_max = atof(optarg); break;
701 case 'v': *pV = atof(optarg); break;
702 case 'T': env->T = atof(optarg); break;
703 case 'C': env->T = atof(optarg)+273.15; break;
705 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
706 assert(num_strings == 1 || num_strings == 2);
707 if (num_strings == 1)
710 ftension_group = atoi(strings[1]);
711 ftension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
712 ftension_params = tension_models[ftension_model].params;
713 free_parsed_list(num_strings, strings);
716 ftension_params = optarg;
719 parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
720 assert(num_strings == 1 || num_strings == 2);
721 if (num_strings == 1)
724 utension_group = atoi(strings[1]);
725 utension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
726 utension_params = tension_models[utension_model].params;
727 free_parsed_list(num_strings, strings);
730 utension_params = optarg;
733 k_model = index_k_model(n_k_models, k_models, optarg);
736 k_models[k_model].params = optarg;
741 for (i=0; i<n; i++) {
742 push(pDomain_list, generate_domain(DS_FOLDED,
743 tension_models+ftension_model,
746 tension_models+utension_model,
749 k_models+k_model), 1);
756 for (i=0; i<n; i++) {
757 push(pDomain_list, generate_domain(DS_UNFOLDED,
758 tension_models+ftension_model,
761 tension_models+utension_model,
764 k_models+k_model), 1);
768 flags = FLAG_OUTPUT_FULL_CURVE;
771 fprintf(stderr, "unrecognized option '%c'\n", optopt);
772 /* fall through to default case */
774 help(prog_name, *pP, *pDt_max, *pV, env, n_tension_models, tension_models, ftension_model, utension_model, n_k_models, k_models, k_model);
778 /* check the arguments */
779 assert(*pDomain_list != NULL); /* no domains! */
780 assert(*pP > 0.0); assert(*pP < 1.0);
781 assert(*pDt_max > 0.0);
783 assert(env->T > 0.0);
788 <<index tension model>>=
789 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
792 for (i=0; i<n_models; i++)
793 if (STRMATCH(models[i].name, name))
796 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
804 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
807 for (i=0; i<n_models; i++)
808 if (STRMATCH(models[i].name, name))
811 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
818 <<initialize model definition>>=
819 /* requires int num_param_args and char **param_args in the current scope
821 * INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
822 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
824 #define INIT_MODEL(role, model, param_string, param_pointer) \
826 parse_list_string(param_string, SEP, '{', '}', \
827 &num_param_args, ¶m_args); \
828 if (num_param_args != model->num_params) { \
830 "%s model %s expected %d params,\n", \
831 role, model->name, model->num_params); \
833 "not the %d params in '%s'\n", \
834 num_param_args, param_string); \
835 assert(num_param_args == model->num_params); \
837 if (model->creator) \
838 param_pointer = (*model->creator)(param_args); \
840 param_pointer = NULL; \
841 free_parsed_list(num_param_args, param_args); \
846 void *generate_domain(enum domain_state_t state,
847 tension_model_getopt_t *folded_model,
849 char *folded_param_string,
850 tension_model_getopt_t *unfolded_model,
852 char *unfolded_param_string,
853 k_model_getopt_t *k_model)
855 void *folded_params, *unfolded_params, *k_params;
856 int num_param_args; /* for INIT_MODEL() */
857 char **param_args; /* for INIT_MODEL() */
860 fprintf(stderr, "generating %s ", state==DS_FOLDED ? "folded" : "unfolded");
861 fprintf(stderr, "domain (k: \"%s\", \"%s\", f: \"%s\" (%p) group %d \"%s\", u: \"%s\" (%p) group %d \"%s\")\n",
862 k_model->name, k_model->params,
863 folded_model->name, folded_model->handler, folded_group, folded_param_string,
864 unfolded_model->name, unfolded_model->handler, unfolded_group, unfolded_param_string);
866 INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
867 INIT_MODEL("unfolded", unfolded_model, unfolded_param_string, unfolded_params);
868 INIT_MODEL("k", k_model, k_model->params, k_params);
869 return create_domain(state,
870 k_model->k, k_params, k_model->destructor,
871 folded_model->handler, folded_group, folded_params, folded_model->destructor,
872 unfolded_model->handler, unfolded_group, unfolded_params, unfolded_model->destructor
879 \addcontentsline{toc}{section}{Appendicies}
881 \section{sawsim.c details}
885 The general layout of our simulation code is:
895 We include [[math.h]], so don't forget to link to the libm with `-lm'.
897 #include <assert.h> /* assert() */
898 #include <stdlib.h> /* malloc(), free(), rand() */
899 #include <stdio.h> /* fprintf(), stdout */
900 #include <string.h> /* strlen, strtok() */
901 #include <math.h> /* exp(), M_PI, sqrt() */
902 #include <sys/types.h> /* pid_t (returned by getpid()) */
903 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
906 #include "tension_balance.h"
908 #include "tension_model.h"
914 <<version definition>>
916 <<max/min definitions>>
917 <<string comparison definition and globals>>
918 <<initialize model definition>>
919 <<get options definitions>>
920 <<domain definitions>>
929 <<create/destroy domain>>
930 <<destroy domain list>>
932 <<simulation functions>>
933 <<boilerplate functions>>
936 <<boilerplate functions>>=
938 <<get option functions>>
944 srand(getpid()*time(NULL)); /* seed rand() */
948 <<flag definitions>>=
949 /* in octal b/c of prefixed '0' */
950 #define FLAG_OUTPUT_FULL_CURVE 01
951 #define FLAG_OUTPUT_UNFOLDING_FORCES 02
955 static unsigned long int flags = 0;
958 \subsection{Utilities}
961 <<max/min definitions>>=
962 #define MAX(a,b) ((a)>(b) ? (a) : (b))
963 #define MIN(a,b) ((a)<(b) ? (a) : (b))
966 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
967 <<string comparison definition and globals>>=
968 // Check if two strings match, return 1 if they do
969 static char *temp_string_A;
970 static char *temp_string_B;
971 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
972 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
973 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
974 /* +1 to also compare the '\0' */
977 We also define a macro for our [[check]] unit testing
978 <<check relative error>>=
979 #define CHECK_ERR(max_err, expected, received) \
981 fail_unless( (received-expected)/expected < max_err, \
982 "relative error %g >= %g in %s (Expected %g, received %g)", \
983 (received-expected)/expected, max_err, #received, \
984 expected, received); \
985 fail_unless(-(received-expected)/expected < max_err, \
986 "relative error %g >= %g in %s (Expected %g, received %g)", \
987 -(received-expected)/expected, max_err, #received, \
988 expected, received); \
993 int happens(double probability)
995 assert(probability >= 0.0); assert(probability <= 1.0);
996 return (double)rand()/RAND_MAX < probability; /* TODO: replace with GSL rand http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html x*/
1001 \subsection{Adaptive timesteps}
1002 \label{app.adaptive_dt}
1004 $F(x)$ increases with $x$, possibly exploding, as in the worm-like chain model,
1005 so basing the timestep on the the unfolding probability at the current tension
1006 is dangerous, and we need to search for a $dt$ for which
1007 $P(F(x+v*dt)) < P_\text{target}$.
1008 There are two cases to consider.
1009 In the most common, no domains have unfolded since the last step,
1010 and we expect the next step to be slightly shorter than the previous one.
1011 In the less common, domains did unfold in the last step,
1012 and we expect the next step to be considerably longer than the previous one.
1014 double search_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers,
1015 list_t *domain_list,
1016 environment_t *env, double x,
1017 double target_prob, double max_dt, double v)
1018 { /* Returns the timestep dt in seconds for the current folded domain.
1019 * Takes a list of tension handlers, the list of domains,
1020 * a pointer env to the environmental data, a starting separation x in m,
1021 * a target_prob between 0 and 1,
1022 * max_dt in s, stretching velocity v in m/s.
1024 double F, k, dtCur, dtU, dtUCur, dtL, dt;
1026 /* get upper bound using the current position */
1027 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x, 0); /* BUG. repeated calculation */
1028 //printf("Start with x = %g (F = %g)\n", x, F);
1029 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1030 //printf("x %g\tF %g\tk %g\n", x, F, k);
1031 dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */
1033 //printf("overshot max_dt\n");
1036 /* set a lower bound on dt too */
1039 /* The dt determined above may produce illegitimate forces or ks.
1040 * Reduce the upper bound until we have valid ks. */
1042 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1043 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
1046 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1048 //printf("Try for dt = %g (F = %g)\n", dt, F);
1049 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1050 /* returned k may be -1.0 */
1051 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1052 while (k == -1.0) { /* reduce step until we hit a valid k */
1054 dt = dtU; /* hopefully, we can use the max dt, see if it works */
1055 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1056 //printf("Try for dt = %g (F = %g)\n", dt, F);
1057 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1058 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
1060 assert(dtU > 1e-14); /* timestep to valid k too small */
1061 dtUCur = target_prob / k; /* safe timestep back from x+dtU */
1063 return dt; /* dtU is safe. */
1065 /* dtU wasn't safe, lets see what would be. */
1066 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
1067 dt = (dtU + dtL) / 2.0;
1068 F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
1069 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
1070 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
1071 dtCur = target_prob / k;
1072 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
1073 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
1075 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
1076 dtU = dt; /* ... stepping out only dtCur would be safe */
1079 break; /* dtCur = dt */
1081 return MAX(dtUCur, dtL);
1085 To determine $dt$ for an array of potentially different folded domains,
1086 we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
1089 double determine_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list,
1090 environment_t *env, double x,
1091 double target_prob, double dt_max, double v)
1092 { /* Returns the timestep dt in seconds.
1093 * Takes the list of folded domains, target_prob between 0 and 1,
1094 * F in N, and T in K. */
1095 double dt=dt_max, new_dt;
1096 assert(target_prob > 0.0); assert(target_prob < 1.0);
1097 assert(dt_max > 0.0);
1099 /* .5 nm steps = v * dt */
1101 while (domain_list != NULL) {
1102 if (D(domain_list)->state == DS_FOLDED) {
1103 new_dt = search_dt(num_tension_handlers, tension_handlers, domain_list, env, x, target_prob, dt, v);
1104 dt = MIN(dt, new_dt);
1106 domain_list = domain_list->next;
1112 \subsection{Domain data}
1114 Currently domains exist in two states, folded and unfolded, and the
1115 only allowed transitions are folded $\rightarrow$ unfolded. Of
1116 course, it wouldn't be too complicated to extent this to a multi-state
1117 system, with an array containing the domains group for each possible
1118 state, and a matrix of transition-rate-calculating functions.
1119 However, at this point such generality seems unnecessary at this
1121 <<domain definitions>>=
1122 enum domain_state_t {DS_FOLDED,
1126 typedef struct domain_struct {
1127 enum domain_state_t state;
1128 one_dim_fn_t *folded_handler;
1130 one_dim_fn_t *unfolded_handler;
1132 k_func_t *k; /* function returning unfolding rate */
1133 void *folded_params; /* pointer to folded parameters */
1134 void *unfolded_params; /* pointer to unfolded parameters */
1135 void *k_params; /* pointer to k parameters */
1136 destroy_data_func_t *destroy_folded;
1137 destroy_data_func_t *destroy_unfolded;
1138 destroy_data_func_t *destroy_k;
1141 /* get the domain data for the current list node */
1142 #define D(list) ((domain_t *)(list)->d)
1143 /* get the tension params for the current list node */
1144 #define D_TP(list) (((domain_t *)(list)->d)->state == DS_FOLDED \
1145 ? ((domain_t *)(list)->d)->folded_params \
1146 : ((domain_t *)(list)->d)->unfolded_params)
1148 [[k]] is a pointer to the function determining the unfolding rate for a given tension.
1149 [[folded_params]] is a pointer to the parameters used by the function pointed to by [[k]].
1150 [[unfolded_params]] is a pointer to the parameters used by the group-appropriate handler function when determining the tension.
1151 The [[destroy_*]] pointers point to functions for freeing the memory [[*_params]].
1152 We store them with the domain data so that [[destroy_domain]] doesn't have to know which type of domain it's cleaning up after.
1154 [[create_]] and [[destroy_domain]] are simple wrappers around [[malloc]] and [[free]].
1155 <<create/destroy domain>>=
1156 domain_t *create_domain(enum domain_state_t state,
1159 destroy_data_func_t *destroy_k,
1160 one_dim_fn_t *folded_handler,
1162 void *folded_params,
1163 destroy_data_func_t *destroy_folded,
1164 one_dim_fn_t *unfolded_handler,
1166 void *unfolded_params,
1167 destroy_data_func_t *destroy_unfolded)
1169 domain_t *ret = (domain_t *)malloc(sizeof(domain_t));
1170 assert(ret != NULL);
1171 if (state == DS_FOLDED) {
1172 assert(k != NULL); /* the pointer points somewhere valid */
1173 assert(*k != NULL); /* and there is something useful there */
1175 assert(state == DS_UNFOLDED);
1177 ret->folded_handler = folded_handler;
1178 ret->folded_group = folded_group;
1179 ret->unfolded_handler = unfolded_handler;
1180 ret->unfolded_group = unfolded_group;
1182 ret->folded_params = folded_params;
1183 ret->unfolded_params = unfolded_params;
1184 ret->k_params = k_params;
1185 ret->destroy_folded = destroy_folded;
1186 ret->destroy_unfolded = destroy_unfolded;
1187 ret->destroy_k = destroy_k;
1189 fprintf(stderr, "generated domain %p. state %d\n", ret, ret->state);
1194 void destroy_domain(domain_t *domain)
1197 //printf("domain %p & %p\n", *domain, domain);
1198 if (domain->destroy_folded)
1199 (*domain->destroy_folded)(domain->folded_params);
1200 if (domain->destroy_unfolded)
1201 (*domain->destroy_unfolded)(domain->unfolded_params);
1202 if (domain->destroy_k)
1203 (*domain->destroy_k)(domain->k_params);
1209 <<destroy domain list>>=
1210 void destroy_domain_list(list_t *domain_list)
1212 domain_list = head(domain_list);
1213 while (domain_list != NULL)
1214 destroy_domain((domain_t *) pop(&domain_list));
1218 \subsection{Domain model and group handling}
1220 <<group functions>>=
1221 <<get tension handler>>
1223 <<int comparison function>>
1224 <<find possible groups>>
1227 <<get active groups>>
1230 <<get tension handler>>=
1231 one_dim_fn_t *get_tension_handler(domain_t *domain)
1233 if (domain->state == DS_FOLDED)
1234 return domain->folded_handler;
1236 if (domain->state != DS_UNFOLDED)
1237 fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
1238 assert(domain->state == DS_UNFOLDED);
1239 return domain->unfolded_handler;
1245 int get_group(domain_t *domain)
1247 if (domain->state == DS_FOLDED)
1248 return domain->folded_group;
1250 if (domain->state != DS_UNFOLDED)
1251 fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
1252 assert(domain->state == DS_UNFOLDED);
1253 return domain->unfolded_group;
1258 We already know all possible tension classes, since they are all
1259 hardcoded into \prog. However, there may be any number of possible
1260 groups. We define a function [[find_possible_groups]] to search for
1261 possible (not neccessarily active) groups. It can search for
1262 subgroups of a single [[handler]], or by setting [[handler]] to
1263 [[NULL]], subgroups of \emph{any} handler. It returns a list of group
1265 <<find possible groups>>=
1266 list_t *find_possible_groups(list_t *list, one_dim_fn_t *handler) {
1269 while (list != NULL) {
1270 if (handler == NULL || get_tension_handler(D(list)) == handler) {
1271 push(&ret, &D(list)->folded_group, 1);
1272 push(&ret, &D(list)->unfolded_group, 1);
1278 return ret; /* no groups with this handler, no processing needed */
1280 /* sort the ret list, and remove duplicates */
1281 sort(&ret, &int_comparison);
1282 unique(&ret, &int_comparison);
1284 fprintf(stderr, "handler %p possible groups:", handler);
1286 while (list != NULL) {
1287 fprintf(stderr, " %d", *((int *)list->d));
1290 fprintf(stderr, "\n");
1296 Search a [[list]] of domains, and determine whether a particular model
1297 class and group number combination is active.
1298 <<is group active>>=
1299 int is_group_active(list_t *list, one_dim_fn_t *handler, int group)
1302 while (list != NULL) {
1303 if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group)
1311 Search a [[list]] of domains, and return all domains matching a
1312 particular model class and group number.
1314 list_t *get_group_list(list_t *list, one_dim_fn_t *handler, int group)
1318 while (list != NULL) {
1319 if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group) {
1320 push(&ret, D_TP(list), 1); /* add a pointer to the appropriate tension parameters to our new list. */
1322 fprintf(stderr, "push domain %p %s tension parameters %p onto active group %p %d list %p\n", D(list), D(list)->state == DS_FOLDED ? "folded" : "unfolded", D_TP(list), handler, group, ret);
1330 Because all the node data in lists returned by [[get_group_list]] is
1331 also in the main domain list, you shouldn't destroy the node data
1332 popped off when destroying the group lists. It will all get cleaned
1333 up when the main domain list is destroyed.
1335 Put all this together to scan through a list of domains and construct
1336 an array of all the active groups.
1337 <<get active groups>>=
1338 void get_active_groups(list_t *list,
1339 int num_tension_handlers, one_dim_fn_t **tension_handlers,
1340 int *pNum_active_groups, one_dim_fn_t ***pPer_group_handlers, void ***pActive_groups)
1342 void **active_groups=NULL;
1343 one_dim_fn_t *handler, **per_group_handlers=NULL;
1344 int i, num_active_groups, *group;
1345 list_t *possible_group_numbers, *active_group_list, *active_groups_list=NULL, *per_group_handlers_list=NULL;
1347 /* get all the active groups in a list */
1348 for (i=0; i<num_tension_handlers; i++) {
1349 handler = tension_handlers[i];
1350 if (handler == NULL) continue; /* tensionless handler */
1351 possible_group_numbers = head(find_possible_groups(list, handler));
1352 while (possible_group_numbers != NULL) {
1353 group = (int *)pop(&possible_group_numbers);
1354 if (is_group_active(list, handler, *group) == 1) {
1355 active_group_list = get_group_list(list, handler, *group);
1356 push(&active_groups_list, active_group_list, 1);
1357 push(&per_group_handlers_list, handler, 1);
1359 fprintf(stderr, "active group %p for %p %d headed by tension parameters %p\n", active_group_list, handler, *group, head(active_group_list)->d);
1360 list_print(stderr, active_group_list, "active group");
1366 /* allocate the array we'll be returning */
1367 num_active_groups = list_length(active_groups_list);
1368 active_groups = (void **) malloc(sizeof(void *)*num_active_groups);
1369 assert(active_groups != NULL);
1370 per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups);
1371 assert(per_group_handlers != NULL);
1373 /* populate the array we'll be returning */
1374 active_groups_list = head(active_groups_list);
1375 for (i=0; i<num_active_groups; i++) {
1376 handler = pop(&per_group_handlers_list);
1377 per_group_handlers[i] = handler;
1379 active_groups[i] = (void *) malloc(sizeof(tension_handler_data_t));
1380 assert(active_groups[i] != NULL);
1381 active_group_list = pop(&active_groups_list);
1382 ((tension_handler_data_t *)active_groups[i])->group = active_group_list;
1383 ((tension_handler_data_t *)active_groups[i])->env = NULL;
1384 ((tension_handler_data_t *)active_groups[i])->persist = NULL;
1386 assert(active_groups_list == NULL);
1387 assert(per_group_handlers_list == NULL);
1389 *pNum_active_groups = num_active_groups;
1390 *pPer_group_handlers = per_group_handlers;
1391 *pActive_groups = active_groups;
1396 \section{String parsing}
1398 For handling command line arguments, we need parse delimited strings (e.g. [["param1,param2,param3..."]]).
1399 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1400 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1401 need to take care of parsing those parameters themselves.
1402 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1406 <<parse definitions>>
1407 <<parse declarations>>
1410 <<parse module makefile lines>>=
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, xo_guess, xo, lb, ub;
1871 double min_relx=1e-6, min_rely=1e-6;
1872 int max_steps=200, 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 //fprintf(stderr, "balancing for x = %g with ", x);
1883 //for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
1884 //fprintf(stderr, "\n");
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) xo_guess = length_scale;
1899 else xo_guess = x/num_tension_groups;
1901 fprintf(stderr, "%s:\tfind a bracket for x = %g with guess of x0 = %g\n", __FUNCTION__, x, xo_guess);
1903 oneD_bracket(x_of_xo, &x_xo_data, x, xo_guess, &lb, &ub);
1904 } else { /* work off the last balanced point */
1906 fprintf(stderr, "%s:\tstep off the old bracketing x %g + %g = %g\n", __FUNCTION__, 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 fprintf(stderr, "not moving\n");
1916 lb= x_xo_data.xi[0];
1917 ub= x_xo_data.xi[0];
1921 fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n",
1922 __FUNCTION__, x, lb, ub);
1924 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL);
1925 F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */
1927 if (num_tension_groups > 1) {
1928 fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x);
1929 for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
1930 fprintf(stderr, "\n");
1934 F = (*tension_handlers[0])(xo, params[0]);
1939 <<tension balance internal definitions>>=
1940 <<x of x0 definitions>>
1943 <<x of x0 definitions>>=
1944 typedef struct x_of_xo_data_struct {
1946 one_dim_fn_t **tension_handler; /* array of fn pointers */
1947 void **handler_data; /* array of void* pointers */
1953 double x_of_xo(double xo, void *pdata)
1955 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
1956 double F, x=0, xi, xi_guess, lb, ub;
1958 double min_relx=1e-6, min_rely=1e-6;
1960 assert(data->n_groups > 0);
1962 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
1964 fprintf(stderr, "%s: finding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
1967 if (data->xi) data->xi[0] = xo;
1968 for (i=1; i < data->n_groups; i++) {
1969 if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */
1970 else xi_guess = data->xi[i];
1972 fprintf(stderr, "%s: searching for proper x[%d] for F[%d](x[%d]) = %g N\n", __FUNCTION__, i, i, i, F);
1974 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub);
1975 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL);
1977 fprintf(stderr, "%s: found F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
1981 if (data->xi) data->xi[i] = xi;
1984 fprintf(stderr, "%s: found x(x0=%g) = %g\n", __FUNCTION__, xo, x);
1990 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited
1991 number of steps for monotonically increasing $f(x)$. Simple
1992 bisection, so it's robust and converges fairly quickly. You can set
1993 the maximum number of steps to take, but if you have enough processor
1994 power, it's better to limit your precision with the relative limits
1995 [[min_relx]] and [[y]]. These ensure that the width of the bracket is
1996 small on the length scale of it's larger side. Note that \emph{both}
1997 relative limits must be satisfied for the search to stop.
1998 <<one dimensional function definition>>=
1999 /* equivalent to gsl_func_t */
2000 typedef double one_dim_fn_t(double x, void *params);
2002 <<one dimensional solve declaration>>=
2003 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2004 double min_relx, double min_rely, int max_steps,
2008 <<one dimensional solve>>=
2009 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
2010 double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub,
2011 double min_relx, double min_rely, int max_steps,
2014 double yx, ylb, yub, x;
2017 ylb = (*f)(lb, params);
2018 yub = (*f)(ub, params);
2020 /* check some simple cases */
2022 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bounded */
2023 else return lb; /* any x on the range [lb, ub] would work */
2025 if (ylb == y) { x=lb; yx=ylb; goto end; }
2026 if (yub == y) { x=ub; yx=yub; goto end; }
2029 fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub);
2031 assert(ylb < y); assert(yub > y);
2032 for (i=0; i<max_steps; i++) {
2034 yx = (*f)(x, params);
2036 fprintf(stderr, "%s:\tbracket: lb %g, x %g, ub %g\tylb %g, yx %g, yub %g\t(y %g)\n", __FUNCTION__, lb, x, ub, ylb, yx, yub, y);
2038 if (yx == y || ((ub-lb)/ub < min_relx && (yub-ylb)/yub < min_rely)) break;
2039 else if (yx > y) { ub=x; yub=yx; }
2040 else /* yx < y */ { lb=x; ylb=yx; }
2045 fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y);
2051 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
2052 Generate bracketing $x$ values through bisection or exponential growth.
2053 <<one dimensional bracket declaration>>=
2054 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
2057 <<one dimensional bracket>>=
2058 void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub)
2060 double yx, step, x=xguess;
2062 yx = (*f)(x, params);
2064 fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params);
2071 fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x);
2075 if (x == 0) x = length_scale/2.0; /* HACK */
2078 yx = (*f)(x, params);
2080 fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx);
2085 fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x);
2088 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
2092 \subsection{Balancing implementation}
2094 <<tension-balance.c>>=
2097 <<tension balance includes>>
2098 #include "tension_balance.h"
2100 double length_scale = 1e-10; /* HACK */
2102 <<tension balance internal definitions>>
2103 <<tension balance functions>>
2106 <<tension balance includes>>=
2107 #include <assert.h> /* assert() */
2108 #include <stdlib.h> /* NULL */
2109 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
2110 #include <stdio.h> /* fprintf(), stdout */
2114 \subsection{Balancing unit tests}
2116 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2117 <<check-tension-balance.c>>=
2118 <<tension balance check includes>>
2119 <<tension balance test suite>>
2120 <<main check program>>
2123 <<tension balance check includes>>=
2124 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2125 #include <assert.h> /* assert() */
2128 #include "tension_balance.h"
2131 <<tension balance test suite>>=
2132 <<tension balance function tests>>
2133 <<tension balance suite definition>>
2136 <<tension balance suite definition>>=
2137 Suite *test_suite (void)
2139 Suite *s = suite_create ("tension balance");
2140 <<tension balance function test case defs>>
2142 <<tension balance function test case adds>>
2147 <<tension balance function tests>>=
2148 <<check relative error>>
2150 double hooke(double x, void *pK)
2153 return *((double*)pK) * x;
2156 START_TEST(test_single_function)
2158 double k=5, x=3, last_x=2.0, F;
2159 one_dim_fn_t *handlers[] = {&hooke};
2160 void *data[] = {&k};
2161 F = tension_balance(1, handlers, data, last_x, x);
2162 fail_unless(F = k*x, NULL);
2167 We can also test balancing two springs with different spring constants.
2168 The analytic solution for a general number of springs is given in Appendix \ref{app.math_hooke}.
2169 <<tension balance function tests>>=
2170 START_TEST(test_double_hooke)
2172 double k1=5, k2=4, x=3, last_x=2.0, F, Fe, x1e, x2e;
2173 one_dim_fn_t *handlers[] = {&hooke, &hooke};
2174 void *data[] = {&k1, &k2};
2175 F = tension_balance(2, handlers, data, last_x, x);
2179 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
2180 //CHECK_ERR(1e-6, x1e, xi[0]);
2181 //CHECK_ERR(1e-6, x2e, xi[1]);
2182 CHECK_ERR(1e-6, Fe, F);
2187 <<tension balance function test case defs>>=
2188 TCase *tc_tbfunc = tcase_create("tension balance function");
2191 <<tension balance function test case adds>>=
2192 tcase_add_test(tc_tbfunc, test_single_function);
2193 tcase_add_test(tc_tbfunc, test_double_hooke);
2194 suite_add_tcase(s, tc_tbfunc);
2199 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
2200 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
2201 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
2205 <<list definitions>>
2206 <<list declarations>>
2209 <<list declarations>>=
2210 <<head and tail declarations>>
2211 <<list length declaration>>
2212 <<push declaration>>
2214 <<list destroy declaration>>
2215 <<list to array declaration>>
2216 <<move declaration>>
2217 <<sort declaration>>
2218 <<unique declaration>>
2219 <<list print declaration>>
2223 <<create and destroy node>>
2236 <<list module makefile lines>>=
2238 notangle -Rlist.c $^ > $@
2240 notangle -Rlist.h $^ > $@
2241 check_list.c : sawsim.nw
2242 notangle -Rcheck-list.c $^ > $@
2243 check_list : check_list.c global.h list.c list.h
2244 gcc -g -o $@ $< list.c -lcheck
2246 rm -f list.c list.h check_list.c check_list
2249 <<list definitions>>=
2250 typedef struct list_struct {
2251 struct list_struct *next;
2252 struct list_struct *prev;
2256 <<comparison function typedef>>
2259 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
2260 <<head and tail declarations>>=
2261 list_t *head(list_t *list);
2262 list_t *tail(list_t *list);
2265 list_t *head(list_t *list)
2269 while (list->prev != NULL) {
2275 list_t *tail(list_t *list)
2279 while (list->next != NULL) {
2286 <<list length declaration>>=
2287 int list_length(list_t *list);
2290 int list_length(list_t *list)
2297 while (list->next != NULL) {
2305 [[push]] inserts a node relative to the current position in the list
2306 without changing the current position.
2307 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
2308 so we set it to that of the pushed domain.
2309 If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c),
2310 otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
2311 <<push declaration>>=
2312 void push(list_t **pList, void *data, int below);
2315 void push(list_t **pList, void *data, int below)
2317 list_t *list, *new_node;
2318 assert(pList != NULL);
2320 new_node = create_node(data);
2323 else if (below==0) { /* insert above */
2324 if (list->prev != NULL)
2325 list->prev->next = new_node;
2326 new_node->prev = list->prev;
2327 list->prev = new_node;
2328 new_node->next = list;
2329 } else { /* insert below */
2330 if (list->next != NULL)
2331 list->next->prev = new_node;
2332 new_node->next = list->next;
2333 list->next = new_node;
2334 new_node->prev = list;
2339 [[pop]] removes the current domain node, moving the current position
2340 to the node after it, unless that node is [[NULL]], in which case move
2341 the current position to the node before it.
2342 <<pop declaration>>=
2343 void *pop(list_t **pList);
2346 void *pop(list_t **pList)
2348 list_t *list, *popped;
2350 assert(pList != NULL);
2352 assert(list != NULL); /* not an empty list */
2354 /* bypass the popped node */
2355 if (list->prev != NULL)
2356 list->prev->next = list->next;
2357 if (list->next != NULL) {
2358 list->next->prev = list->prev;
2359 *pList = list->next;
2361 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
2363 destroy_node(popped);
2368 [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data.
2369 If the cleanup function is [[NULL]], no cleanup function is called.
2370 The nodes are not popped in any particular order.
2371 <<list destroy declaration>>=
2372 void list_destroy(list_t **pList, destroy_data_func_t *destroy);
2375 void list_destroy(list_t **pList, destroy_data_func_t *destroy)
2379 assert(pList != NULL);
2382 assert(list != NULL); /* not an empty list */
2383 while (list != NULL) {
2385 if (destroy != NULL)
2391 [[list_to_array]] converts a list to an array of pointers, preserving order.
2392 <<list to array declaration>>=
2393 void list_to_array(list_t *list, int *length, void ***array);
2396 void list_to_array(list_t *list, int *pLength, void ***pArray)
2400 assert(list != NULL);
2401 assert(pLength != NULL);
2402 assert(pArray != NULL);
2404 length = list_length(list);
2405 array = (void **)malloc(sizeof(void *)*length);
2406 assert(array != NULL);
2409 while (list != NULL) {
2410 array[i++] = list->d;
2418 [[move]] moves the current element along the list in either direction.
2419 <<move declaration>>=
2420 void move(list_t **pList, int downstream);
2423 void move(list_t **pList, int downstream)
2425 list_t *list, *flipper;
2427 assert(pList != NULL);
2428 list = *pList; /* a>B>c>d */
2429 assert(list != NULL); /* not an empty list */
2430 if (downstream == 0)
2431 flipper = list->prev; /* flipper is A */
2433 flipper = list->next; /* flipper is C */
2434 /* ensure there is somewhere to go in the direction we're heading */
2435 assert(flipper != NULL);
2436 /* remove the current element */
2437 data = pop(&list); /* data=B, a>C>d */
2438 /* and put it back in in the correct spot */
2440 if (downstream == 0) {
2441 push(&list, data, 0); /* b>A>c>d */
2442 list = list->prev; /* B>a>c>d */
2444 push(&list, data, 1); /* a>C>b>d */
2445 list = list->next; /* a>c>B>d */
2451 [[sort]] sorts a list from smallest to largest according to a user
2452 specified comparison.
2453 <<comparison function typedef>>=
2454 typedef int (comparison_fn_t)(void *A, void *B);
2457 <<int comparison function>>=
2458 int int_comparison(void *A, void *B)
2463 if (a > b) return 1;
2464 else if (a == b) return 0;
2468 <<sort declaration>>=
2469 void sort(list_t **pList, comparison_fn_t *comp);
2471 Here we use the bubble sort algorithm.
2473 void sort(list_t **pList, comparison_fn_t *comp)
2477 assert(pList != NULL);
2479 assert(list != NULL); /* not an empty list */
2480 while (stable == 0) {
2483 while (list->next != NULL) {
2484 if ((*comp)(list->d, list->next->d) > 0) {
2486 move(&list, 1 /* downstream */);
2495 [[unique]] removes duplicates from a sorted list according to a user
2496 specified comparison. Don't do this unless you have other pointers
2497 any dynamically allocated data in your list, because [[unique]] just
2498 drops any non-unique data without freeing it.
2499 <<unique declaration>>=
2500 void unique(list_t **pList, comparison_fn_t *comp);
2503 void unique(list_t **pList, comparison_fn_t *comp)
2506 assert(pList != NULL);
2508 assert(list != NULL); /* not an empty list */
2510 while (list->next != NULL) {
2511 if ((*comp)(list->d, list->next->d) == 0) {
2520 [[list_print]] prints a list to a [[FILE *]] stream.
2521 <<list print declaration>>=
2522 void list_print(FILE *file, list_t *list, char *list_name);
2525 void list_print(FILE *file, list_t *list, char *list_name)
2527 fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name);
2529 while (list != NULL) {
2530 fprintf(file, " %p", list->d);
2533 fprintf(file, "\n");
2537 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
2538 <<create and destroy node>>=
2539 list_t *create_node(void *data)
2541 list_t *ret = (list_t *)malloc(sizeof(list_t));
2542 assert(ret != NULL);
2549 void destroy_node(list_t *node)
2555 The user must free the data pointed to by the node on their own.
2557 \subsection{List implementation}
2567 #include <assert.h> /* assert() */
2568 #include <stdlib.h> /* malloc(), free(), rand() */
2569 #include <stdio.h> /* fprintf(), stdout, FILE */
2570 #include "global.h" /* destroy_data_func_t */
2573 \subsection{List unit tests}
2575 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2577 <<list check includes>>
2579 <<main check program>>
2582 <<list check includes>>=
2583 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2588 <<list test suite>>=
2591 <<list suite definition>>
2594 <<list suite definition>>=
2595 Suite *test_suite (void)
2597 Suite *s = suite_create ("list");
2598 <<push test case defs>>
2599 <<pop test case defs>>
2601 <<push test case adds>>
2602 <<pop test case adds>>
2608 START_TEST(test_push)
2613 push(&list, (void *)i, 1);
2614 fail_unless(list == head(list), NULL);
2615 fail_unless( (int)list->d == 0 );
2616 for (i=0; i<n; i++) {
2617 /* we expect 0, n-1, n-2, ..., 1 */
2620 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
2626 <<push test case defs>>=
2627 TCase *tc_push = tcase_create("push");
2630 <<push test case adds>>=
2631 tcase_add_test(tc_push, test_push);
2632 suite_add_tcase(s, tc_push);
2637 <<pop test case defs>>=
2639 <<pop test case adds>>=
2642 \section{Function string evaluation}
2644 For the saddle-point approximated Kramers' model (Section \ref{sec.kramers}) we need the ability to evaluate user-supplied functions ($E(x)$, $x_{ts}(F)$, \ldots).
2645 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
2646 The solution is to run a scripting language as a subprocess accessed via pipes.
2647 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
2649 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
2650 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
2651 Most of this is also POSIX-specific (as far as I know), so you'll have to do some legwork to port it to non-POSIX systems.
2652 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
2654 If you feel the benefits do \emph{not} outweigh the costs, or you are on a non-POSIX system (e.g. MS Windows without Cygwin), you should probably hardcode your functions in \lang.
2655 Then you can either statically or dynamically link to those hardcoded functions.
2656 While much less flexible, this approach would be a fairly simple-to-implement fallback.
2658 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
2659 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
2663 <<string eval setup declaration>>
2664 <<string eval function declaration>>
2665 <<string eval teardown declaration>>
2668 <<string eval module makefile lines>>=
2669 string_eval.c : sawsim.nw
2670 notangle -Rstring-eval.c $^ > $@
2671 string_eval.h : sawsim.nw
2672 notangle -Rstring-eval.h $^ > $@
2673 check_string_eval.c : sawsim.nw
2674 notangle -Rcheck-string-eval.c $^ > $@
2675 check_string_eval : check_string_eval.c string_eval.c string_eval.h
2676 gcc -g -o $@ $< string_eval.c -lcheck -lgsl -lgslcblas -lm
2678 rm -f string_eval.c string_eval.h check_string_eval.c check_string_eval
2681 For an introduction to POSIX process control, see\\
2682 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
2683 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
2684 and of course, the relavent [[man]] pages.
2686 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
2687 [[execvp]] replaces the calling process' program with a new program.
2688 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
2689 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
2690 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
2691 The new program needs command line arguments, just like it would if you were running it from a shell.
2692 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
2693 with the final array entry being a [[NULL]] pointer.
2695 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
2696 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
2697 <<string eval subprocess definitions>>=
2698 #define SUBPROCESS "python"
2699 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
2700 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
2701 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
2703 The [[i]] option lets Python know that it should run in interactive mode.
2704 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
2705 In interactive mode, python acts on every instruction as soon as it is recieved.
2706 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.
2707 %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.
2709 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
2710 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
2711 [[fork]] returns two copies of the same program, executing the original code.
2712 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
2713 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
2715 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.
2716 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
2717 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
2718 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
2719 <<string eval pipe definitions>>=
2720 #define PIPE_READ 0 /* the end you read from */
2721 #define PIPE_WRITE 1 /* the end you write to */
2723 #define STDIN 0 /* initial index of stdin pair */
2724 #define STDOUT 2 /* initial index of stdout pair */
2727 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
2729 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.
2731 <<string eval setup declaration>>=
2732 extern void string_eval_setup(FILE **pIn, FILE **pOut);
2734 <<string eval setup definition>>=
2735 void string_eval_setup(FILE **pIn, FILE **pOut)
2738 int pfd[4]; /* file descriptors for stdin and stdout */
2740 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
2741 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
2743 pid = fork(); /* split process into two copies */
2745 if (pid == -1) { /* fork error */
2746 perror("fork error");
2748 } else if (pid == 0) { /* child */
2749 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
2750 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
2751 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
2752 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
2753 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
2754 perror("exec error"); /* exec shouldn't return */
2756 } else { /* parent */
2757 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
2758 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
2759 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
2760 if ( *pIn == NULL ) {
2761 perror("fdopen (in)");
2764 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
2765 if ( *pOut == NULL ) {
2766 perror("fdopen (out)");
2773 To use the evaluating subprocess, we just pipe in our command, and read out the results.
2774 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
2775 <<string eval function declaration>>=
2776 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
2778 <<string eval function definition>>=
2779 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
2782 rc = fprintf(in, "%s", input);
2783 assert(rc == strlen(input));
2786 alarm(1); /* set a one second timeout on the read */
2787 assert( fgets(output, buflen, out) != NULL );
2788 alarm(0); /* cancel the timeout */
2789 //fprintf(stderr, "eval: %s ----> %s", input, output);
2792 The [[alarm]] calls set and clear a timeout on the returned output.
2793 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.
2794 This protects against invalid input for which a line of output is not printed to [[stdout]].
2795 Other invalid inputs (e.g. those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
2796 If you are getting strange results, check your python code seperately. TODO, better error handling.
2798 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
2799 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
2800 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
2801 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
2802 <<string eval teardown declaration>>=
2803 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
2806 <<string eval teardown definition>>=
2807 void string_eval_teardown(FILE **pIn, FILE **pOut)
2812 /* redirect Python's stderr */
2813 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
2817 assert( fclose(*pIn) == 0 );
2819 assert( fclose(*pOut) == 0 );
2822 /* wait for python to exit */
2824 pid = wait(&stat_loc);
2831 if (WIFEXITED(stat_loc)) {
2832 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
2833 } else if (WIFSIGNALED(stat_loc)) {
2834 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
2839 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
2841 \subsection{String evaluation implementation}
2845 <<string eval includes>>
2846 #include "string_eval.h"
2847 <<string eval internal definitions>>
2848 <<string eval functions>>
2851 <<string eval includes>>=
2852 #include <assert.h> /* assert() */
2853 #include <stdlib.h> /* NULL */
2854 #include <stdio.h> /* fprintf(), stdout, fdopen() */
2855 #include <string.h> /* strlen() */
2856 #include <sys/types.h> /* pid_t */
2857 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
2858 #include <wait.h> /* wait() */
2861 <<string eval internal definitions>>=
2862 <<string eval subprocess definitions>>
2863 <<string eval pipe definitions>>
2866 <<string eval functions>>=
2867 <<string eval setup definition>>
2868 <<string eval function definition>>
2869 <<string eval teardown definition>>
2872 \subsection{String evaluation unit tests}
2874 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2875 <<check-string-eval.c>>=
2876 <<string eval check includes>>
2877 <<string comparison definition and globals>>
2878 <<string eval test suite>>
2879 <<main check program>>
2882 <<string eval check includes>>=
2883 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
2884 #include <stdio.h> /* printf() */
2885 #include <assert.h> /* assert() */
2886 #include <string.h> /* strlen() */
2887 #include <signal.h> /* SIGKILL */
2889 #include "string_eval.h"
2892 <<string eval test suite>>=
2893 <<string eval tests>>
2894 <<string eval suite definition>>
2897 <<string eval suite definition>>=
2898 Suite *test_suite (void)
2900 Suite *s = suite_create ("string eval");
2901 <<string eval test case defs>>
2903 <<string eval test case add>>
2908 <<string eval tests>>=
2909 START_TEST(test_setup_teardown)
2912 string_eval_setup(&in, &out);
2913 string_eval_teardown(&in, &out);
2916 START_TEST(test_invalid_command)
2919 char input[80], output[80]={};
2920 string_eval_setup(&in, &out);
2921 sprintf(input, "print ABCDefg\n");
2922 string_eval(in, out, input, 80, output);
2923 string_eval_teardown(&in, &out);
2926 START_TEST(test_simple_eval)
2929 char input[80], output[80]={};
2930 string_eval_setup(&in, &out);
2931 sprintf(input, "print 3+4*5\n");
2932 string_eval(in, out, input, 80, output);
2933 fail_unless(STRMATCH(output,"23\n"), NULL);
2934 string_eval_teardown(&in, &out);
2937 START_TEST(test_multiple_evals)
2940 char input[80], output[80]={};
2941 string_eval_setup(&in, &out);
2942 sprintf(input, "print 3+4*5\n");
2943 string_eval(in, out, input, 80, output);
2944 fail_unless(STRMATCH(output,"23\n"), NULL);
2945 sprintf(input, "print (3**2 + 4**2)**0.5\n");
2946 string_eval(in, out, input, 80, output);
2947 fail_unless(STRMATCH(output,"5.0\n"), NULL);
2948 string_eval_teardown(&in, &out);
2951 START_TEST(test_eval_with_state)
2954 char input[80], output[80]={};
2955 string_eval_setup(&in, &out);
2956 sprintf(input, "print 3+4*5\n");
2957 fprintf(in, "A = 3\n");
2958 sprintf(input, "print A*3\n");
2959 string_eval(in, out, input, 80, output);
2960 fail_unless(STRMATCH(output,"9\n"), NULL);
2961 string_eval_teardown(&in, &out);
2965 <<string eval test case defs>>=
2966 TCase *tc_string_eval = tcase_create("string_eval");
2968 <<string eval test case add>>=
2969 tcase_add_test(tc_string_eval, test_setup_teardown);
2970 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
2971 tcase_add_test(tc_string_eval, test_simple_eval);
2972 tcase_add_test(tc_string_eval, test_multiple_evals);
2973 tcase_add_test(tc_string_eval, test_eval_with_state);
2974 suite_add_tcase(s, tc_string_eval);
2978 \section{Accelerating function evaluation}
2980 My first version-0.3 code was running very slowly.
2981 With the options suggested in the help
2982 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
2983 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
2984 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
2985 That's only 15 calls per solution, so the search algorithm seems reasonable.
2986 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
2989 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
2993 <<accel k module makefile lines>>=
2994 accel_k.c : sawsim.nw
2995 notangle -Raccel-k.c $^ > $@
2996 accel_k.h : sawsim.nw
2997 notangle -Raccel-k.h $^ > $@
2998 check_accel_k.c : sawsim.nw
2999 notangle -Rcheck-accel_k.c $^ > $@
3000 check_accel_k : check_accel_k.c global.h
3001 gcc -g -o $@ $< accel_k.c -lcheck -lgsl -lgslcblas -lm
3003 rm -f accel_k.c accel_k.h check_accel_k.c check_accel_k
3007 #include <assert.h> /* assert() */
3008 #include <stdlib.h> /* realloc(), free(), NULL */
3009 #include "global.h" /* environment_t */
3010 #include "k_model.h" /* k_func_t */
3011 #include "interp.h" /* interp_* */
3012 #include "accel_k.h"
3014 typedef struct accel_k_struct {
3015 interp_table_t *itable;
3021 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
3022 static int num_accels = 0;
3023 static accel_k_t *accels=NULL;
3025 /* Wrap k in the standard f(x) acceleration form */
3026 static double k_wrap(double F, void *params)
3028 accel_k_t *p = (accel_k_t *)params;
3029 return (*p->k)(F, p->env, p->params);
3032 static int k_tol(double FA, double kA, double FB, double kB)
3035 if (FB-FA > 1e-12) {
3036 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
3037 return 1; /* unnacceptable */
3039 //printf("acceptable tol\n");
3040 return 0; /* acceptable */
3044 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
3047 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
3048 assert(accels != NULL);
3049 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
3051 accels[i].env = env;
3052 accels[i].params = params;
3059 for (i=0; i<num_accels; i++)
3060 interp_table_free(accels[i].itable);
3062 if (accels) free(accels);
3066 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
3069 for (i=0; i<num_accels; i++) {
3070 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
3071 /* We've already setup this function.
3072 * WARNING: we're assuming param and env are the same. */
3073 return interp_table_eval(accels[i].itable, accels+i, F);
3076 /* set up a new acceleration environment */
3077 i = add_accel_k(k, env, params);
3078 return interp_table_eval(accels[i].itable, accels+i, F);
3082 \section{Tension models}
3083 \label{sec.tension_models}
3085 TODO: tension model intro.
3086 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.
3088 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3089 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]].
3091 <<tension-model.h>>=
3093 <<tension handler types>>
3094 <<constant tension model declarations>>
3095 <<hooke tension model declarations>>
3096 <<worm-like chain tension model declarations>>
3097 <<freely-jointed chain tension model declarations>>
3098 <<find tension definitions>>
3099 <<tension model global declarations>>
3102 <<tension model module makefile lines>>=
3103 tension_model.c : sawsim.nw
3104 notangle -Rtension-model.c $^ > $@
3105 tension_model.h : sawsim.nw
3106 notangle -Rtension-model.h $^ > $@
3107 check_tension_model.c : sawsim.nw
3108 notangle -Rcheck-tension-model.c $^ > $@
3109 check_tension_model : check_tension_model.c global.h tension_model.c tension_model.h
3110 gcc -g -o $@ $< tension_model.c -lcheck -lgsl -lgslcblas -lm
3111 clean_tension_model : clean_tension_model_utils
3112 rm -f tension_model.c tension_model.h check_tension_model.c check_tension_model
3113 tension_model_utils.c : sawsim.nw
3114 notangle -Rtension-model-utils.c $^ > $@
3115 tension_model_utils : tension_model_utils.c global.h tension_model.c tension_model.h parse.c parse.h \
3116 list.c list.h tension_balance.c tension_balance.h
3117 gcc -g -o $@ $< tension_model.c parse.c list.c tension_balance.c -lgsl -lgslcblas -lm
3118 tension_model_utils_static : tension_model_utils.c global.h tension_model.c tension_model.h parse.c parse.h \
3119 list.c list.h tension_balance.c tension_balance.h
3120 gcc -g -static -o $@ $< tension_model.c parse.c list.c tension_balance.c -lgsl -lgslcblas -lm
3121 clean_tension_model_utils :
3122 rm -f tension_model_utils.c tension_model_utils
3126 \label{sec.null_tension}
3128 For unstretchable domains.
3130 <<null tension model getopt>>=
3131 {NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
3134 \subsection{Constant}
3135 \label{sec.const_tension}
3137 <<constant tension functions>>=
3138 <<constant tension function>>
3139 <<constant tension parameter create/destroy functions>>
3142 <<constant tension model declarations>>=
3143 <<constant tension function declaration>>
3144 <<constant tension parameter create/destroy function declarations>>
3145 <<constant tension model global declarations>>
3149 An infinitely stretchable domain providing a constant tension.
3151 <<constant tension function declaration>>=
3152 extern double const_tension_handler(double x, void *pdata);
3154 <<constant tension function>>=
3155 double const_tension_handler(double x, void *pdata)
3157 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3158 list_t *list = data->group;
3163 assert(list != NULL); /* empty active group?! */
3164 F = ((const_tension_param_t *)list->d)->F;
3166 fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F);
3168 while (list != NULL) {
3169 assert(((const_tension_param_t *)list->d)->F == F);
3177 <<constant tension structure>>=
3178 typedef struct const_tension_param_struct {
3179 double F; /* tension (force) in N */
3180 } const_tension_param_t;
3184 <<constant tension parameter create/destroy function declarations>>=
3185 extern void *string_create_const_tension_param_t(char **param_strings);
3186 extern void destroy_const_tension_param_t(void *p);
3189 <<constant tension parameter create/destroy functions>>=
3190 const_tension_param_t *create_const_tension_param_t(double F)
3192 const_tension_param_t *ret
3193 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
3194 assert(ret != NULL);
3199 void *string_create_const_tension_param_t(char **param_strings)
3201 return create_const_tension_param_t(atof(param_strings[0]));
3204 void destroy_const_tension_param_t(void *p)
3212 <<constant tension model global declarations>>=
3213 extern int num_const_tension_params;
3214 extern char *const_tension_param_descriptions[];
3215 extern char const_tension_param_string[];
3218 <<constant tension model globals>>=
3219 int num_const_tension_params = 1;
3220 char *const_tension_param_descriptions[] = {"tension F, N"};
3221 char const_tension_param_string[] = "0";
3224 <<constant tension model getopt>>=
3225 {&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}
3231 <<hooke functions>>=
3233 <<hooke parameter create/destroy functions>>
3236 <<hooke tension model declarations>>=
3237 <<hooke tension function declaration>>
3238 <<hooke tension parameter create/destroy function declarations>>
3239 <<hooke tension model global declarations>>
3243 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
3244 The behavior of a series of springs $k_i$ in series is given by
3246 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
3248 For a simple proof, see Appendix \ref{app.math_hooke}.
3250 <<hooke tension function declaration>>=
3251 extern double hooke_handler(double x, void *pdata);
3255 double hooke_handler(double x, void *pdata)
3257 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3258 list_t *list = data->group;
3263 assert(list != NULL); /* empty active group?! */
3264 while (list != NULL) {
3265 assert( ((hooke_param_t *)list->d)->k > 0 );
3266 k += 1.0/ ((hooke_param_t *)list->d)->k;
3268 fprintf(stderr, "%s: Hooke spring %g N/m in series\n",
3269 __FUNCTION__, ((hooke_param_t *)list->d)->k);
3275 fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N\n",
3276 __FUNCTION__, k, x, k*x);
3282 <<hooke structure>>=
3283 typedef struct hooke_param_struct {
3284 double k; /* spring constant in N/m */
3288 <<hooke tension parameter create/destroy function declarations>>=
3289 extern void *string_create_hooke_param_t(char **param_strings);
3290 extern void destroy_hooke_param_t(void *p);
3293 <<hooke parameter create/destroy functions>>=
3294 hooke_param_t *create_hooke_param_t(double k)
3296 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
3297 assert(ret != NULL);
3302 void *string_create_hooke_param_t(char **param_strings)
3304 return create_hooke_param_t(atof(param_strings[0]));
3307 void destroy_hooke_param_t(void *p)
3314 <<hooke tension model global declarations>>=
3315 extern int num_hooke_params;
3316 extern char *hooke_param_descriptions[];
3317 extern char hooke_param_string[];
3320 <<hooke tension model globals>>=
3321 int num_hooke_params = 1;
3322 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
3323 char hooke_param_string[]="0.05";
3326 <<hooke tension model getopt>>=
3327 {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}
3330 \subsection{Worm-like chain}
3333 We can model several unfolded domains with a single worm-like chain.
3334 <<worm-like chain functions>>=
3335 <<worm-like chain function>>
3336 <<worm-like chain handler>>
3337 <<worm-like chain create/destroy functions>>
3340 <<worm-like chain tension model declarations>>=
3341 <<worm-like chain handler declaration>>
3342 <<worm-like chain create/destroy function declarations>>
3343 <<worm-like chain tension model global declarations>>
3347 The combination of all unfolded domains is modeled as a worm like chain,
3348 with the force $F_{\text{WLC}}$ approximately given by
3350 F_{\text{WLC}} \approx \frac{k_B T}{p}
3351 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
3353 where $x$ is the distance between the fixed ends,
3354 $k_B$ is Boltzmann's constant,
3355 $T$ is the absolute temperature,
3356 $p$ is the persistence length, and
3357 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
3358 <<worm-like chain function>>=
3359 double wlc(double x, double T, double p, double L)
3361 double xL=0.0; /* uses global k_B */
3363 assert(T > 0); assert(p > 0); assert(L > 0);
3364 if (x >= L) return HUGE_VAL;
3365 xL = x/L; /* unitless */
3366 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
3367 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
3368 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
3371 This model assumes that each unfolded domain has the same persistence length.
3373 <<worm-like chain handler declaration>>=
3374 extern double wlc_handler(double x, void *pdata);
3377 <<worm-like chain handler>>=
3378 double wlc_handler(double x, void *pdata)
3380 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3381 list_t *list = data->group;
3385 assert(list != NULL); /* empty active group?! */
3386 p = ((wlc_param_t *)list->d)->p;
3387 while (list != NULL) {
3388 /* enforce consistent p */
3389 assert( ((wlc_param_t *)list->d)->p == p);
3390 L += ((wlc_param_t *)list->d)->L;
3392 fprintf(stderr, "%s: WLC section %g m long in series\n",
3393 __FUNCTION__, ((wlc_param_t *)list->d)->L);
3398 fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n",
3399 __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L));
3401 return wlc(x, data->env->T, p, L);
3405 <<worm-like chain structure>>=
3406 typedef struct wlc_param_struct {
3407 double p; /* WLC persistence length (m) */
3408 double L; /* WLC contour length (m) */
3412 <<worm-like chain create/destroy function declarations>>=
3413 extern void *string_create_wlc_param_t(char **param_strings);
3414 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
3417 <<worm-like chain create/destroy functions>>=
3418 wlc_param_t *create_wlc_param_t(double p, double L)
3420 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
3421 assert(ret != NULL);
3427 void *string_create_wlc_param_t(char **param_strings)
3429 return create_wlc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3432 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
3439 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.
3440 TODO (now we copy the string before parsing, but still do this for clarity).
3442 <<valid string write code>>=
3443 char string[] = "abc";
3446 <<valid string write code>>=
3447 char *string = "abc";
3451 <<worm-like chain tension model global declarations>>=
3452 extern int num_wlc_params;
3453 extern char *wlc_param_descriptions[];
3454 extern char wlc_param_string[];
3456 <<worm-like chain tension model globals>>=
3457 int num_wlc_params = 2;
3458 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
3459 char wlc_param_string[]="0.39e-9,28.4e-9";
3462 <<worm-like chain tension model getopt>>=
3463 {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}
3465 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
3467 \subsection{Freely-jointed chain}
3470 <<freely-jointed chain functions>>=
3471 <<freely-jointed chain function>>
3472 <<freely-jointed chain handler>>
3473 <<freely-jointed chain create/destroy functions>>
3476 <<freely-jointed chain tension model declarations>>=
3477 <<freely-jointed chain handler declaration>>
3478 <<freely-jointed chain create/destroy function declarations>>
3479 <<freely-jointed chain tension model global declarations>>
3483 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
3484 The tension of a chain of $N$ such links is given by
3486 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
3488 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}.
3489 <<freely-jointed chain function>>=
3490 double langevin(double x, void *params)
3493 return 1.0/tanh(x) - 1/x;
3495 <<one dimensional bracket declaration>>
3496 <<one dimensional solve declaration>>
3497 double inv_langevin(double x)
3499 double lb=0.0, ub=1.0, ret;
3500 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
3501 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
3505 double fjc(double x, double T, double l, int N)
3507 double L = l*(double)N;
3508 if (x == 0) return 0;
3509 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
3510 return k_B*T/l * inv_langevin(x/L);
3514 <<freely-jointed chain handler declaration>>=
3515 extern double fjc_handler(double x, void *pdata);
3517 <<freely-jointed chain handler>>=
3518 double fjc_handler(double x, void *pdata)
3520 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
3521 list_t *list = data->group;
3526 assert(list != NULL); /* empty active group?! */
3527 l = ((fjc_param_t *)list->d)->l;
3528 while (list != NULL) {
3529 /* enforce consistent l */
3530 assert( ((fjc_param_t *)list->d)->l == l);
3531 N += ((fjc_param_t *)list->d)->N;
3533 fprintf(stderr, "%s: FJC section %d links long in series\n",
3534 __FUNCTION__, ((fjc_param_t *)list->d)->N);
3539 fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n",
3540 __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N));
3542 return fjc(x, data->env->T, l, N);
3546 <<freely-jointed chain structure>>=
3547 typedef struct fjc_param_struct {
3548 double l; /* FJC link length (m) */
3549 int N; /* FJC number of links */
3553 <<freely-jointed chain create/destroy function declarations>>=
3554 extern void *string_create_fjc_param_t(char **param_strings);
3555 extern void destroy_fjc_param_t(void *p);
3558 <<freely-jointed chain create/destroy functions>>=
3559 fjc_param_t *create_fjc_param_t(double l, double N)
3561 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
3562 assert(ret != NULL);
3568 void *string_create_fjc_param_t(char **param_strings)
3570 return create_fjc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3573 void destroy_fjc_param_t(void *p)
3580 <<freely-jointed chain tension model global declarations>>=
3581 extern int num_fjc_params;
3582 extern char *fjc_param_descriptions[];
3583 extern char fjc_param_string[];
3586 <<freely-jointed chain tension model globals>>=
3587 int num_fjc_params = 2;
3588 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
3589 char fjc_param_string[]="0.5e-9,200";
3592 <<freely-jointed chain tension model getopt>>=
3593 {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}
3596 \subsection{Tension model implementation}
3598 <<tension-model.c>>=
3601 <<tension model includes>>
3602 #include "tension_model.h"
3603 <<tension model internal definitions>>
3604 <<tension model globals>>
3605 <<tension model functions>>
3608 <<tension model includes>>=
3609 #include <assert.h> /* assert() */
3610 #include <stdlib.h> /* NULL */
3611 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
3612 #include <stdio.h> /* fprintf(), stdout */
3615 #include "tension_balance.h" /* oneD_*() */
3618 <<tension model internal definitions>>=
3619 <<constant tension structure>>
3621 <<worm-like chain structure>>
3622 <<freely-jointed chain structure>>
3625 <<tension model globals>>=
3626 <<tension handler array global>>
3627 <<constant tension model globals>>
3628 <<hooke tension model globals>>
3629 <<worm-like chain tension model globals>>
3630 <<freely-jointed chain tension model globals>>
3633 <<tension model functions>>=
3634 <<constant tension functions>>
3636 <<worm-like chain functions>>
3637 <<freely-jointed chain functions>>
3641 \subsection{Utilities}
3643 The tension models can get complicated, and users may want to reassure
3644 themselves that this portion of the input physics are functioning properly. The
3645 stand-alone program [[tension_model_utils]] demonstrates the tension model
3646 interface without the overhead of having to understand a full simulation.
3647 [[tension_model_utils]] takes command line model arguments like the full
3648 simulation, and prints $F(x)$ data to the screen for the selected model over a
3651 <<tension-model-utils.c>>=
3653 <<tension model utility includes>>
3654 <<tension model utility definitions>>
3655 <<tension model utility getopt functions>>
3657 <<tension model utility main>>
3660 <<tension model utility main>>=
3661 int main(int argc, char **argv)
3663 <<tension model getopt array definition>>
3664 tension_model_getopt_t *model = NULL;
3668 one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0};
3669 tension_handler_data_t tdata;
3670 int num_param_args; /* for INIT_MODEL() */
3671 char **param_args; /* for INIT_MODEL() */
3672 get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model, &flags);
3674 if (flags & VFLAG) {
3675 printf("#initializing model %s with parameters %s\n", model->name, model->params);
3677 INIT_MODEL("utility", model, model->params, params);
3679 push(&tdata.group, params, 1);
3681 tdata.persist = NULL;
3682 if (model->handler == NULL) {
3683 printf("No tension function for model %s\n", model->name);
3687 double dx=1e-10, x=0, F=0;
3688 printf("#F (N)\tk (%% pop. per s)\n");
3689 while (F >= 0 && F < 1e5 && x < 1e-6) {
3690 F = (*model->handler)(x, &tdata);
3691 printf("%g\t%g\n", x, F);
3695 params = pop(&tdata.group);
3696 if (model->destructor)
3697 (*model->destructor)(params);
3702 <<tension model utility includes>>=
3705 #include <string.h> /* strlen() */
3706 #include <assert.h> /* assert() */
3710 #include "tension_model.h"
3713 <<tension model utility definitions>>=
3714 <<version definition>>
3715 #define VFLAG 1 /* verbose */
3716 <<string comparison definition and globals>>
3717 <<tension model getopt definitions>>
3718 <<initialize model definition>>
3722 <<tension model utility getopt functions>>=
3723 <<index tension model>>
3724 <<tension model utility help>>
3725 <<tension model utility get options>>
3728 <<tension model utility help>>=
3729 void help(char *prog_name,
3731 int n_tension_models, tension_model_getopt_t *tension_models,
3735 printf("usage: %s [options]\n", prog_name);
3736 printf("Version %s\n\n", VERSION);
3737 printf("A utility for understanding the available tension models\n\n");
3738 printf("Environmental options:\n");
3739 printf("-T T\tTemperature T (currently %g K)\n", env->T);
3740 printf("-C T\tYou can also set the temperature T in Celsius\n");
3741 printf("Model options:\n");
3742 printf("-m model\tFolded tension model (currently %s)\n",
3743 tension_models[tension_model].name);
3744 printf("-a args\tFolded tension model argument string (currently %s)\n",
3745 tension_models[tension_model].params);
3746 printf("F(x) is calculated for a range of x and printed\n");
3747 printf("For example:\n");
3748 printf(" #Distance (x)\tForce (N)\n");
3749 printf(" 123.456\t7.89\n");
3751 printf("-V\tChange output to verbose mode\n");
3752 printf("-h\tPrint this help and exit\n");
3754 printf("Tension models:\n");
3755 for (i=0; i<n_tension_models; i++) {
3756 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
3757 for (j=0; j < tension_models[i].num_params ; j++)
3758 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
3759 printf("\t\tdefaults: %s\n", tension_models[i].params);
3764 <<tension model utility get options>>=
3765 void get_options(int argc, char **argv, environment_t *env,
3766 int n_tension_models, tension_model_getopt_t *tension_models,
3767 tension_model_getopt_t **model,
3768 unsigned int *flags)
3770 char *prog_name = NULL;
3771 char c, options[] = "T:C:m:a:Vh";
3772 int tension_model=0, num_strings;
3773 extern char *optarg;
3774 extern int optind, optopt, opterr;
3778 /* setup defaults */
3780 prog_name = argv[0];
3781 env->T = 300.0; /* K */
3783 *model = tension_models;
3785 while ((c=getopt(argc, argv, options)) != -1) {
3787 case 'T': env->T = atof(optarg); break;
3788 case 'C': env->T = atof(optarg)+273.15; break;
3790 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
3791 *model = tension_models+tension_model;
3794 tension_models[tension_model].params = optarg;
3796 case 'V': *flags |= VFLAG; break;
3798 fprintf(stderr, "unrecognized option '%c'\n", optopt);
3799 /* fall through to default case */
3801 help(prog_name, env, n_tension_models, tension_models, tension_model);
3810 \section{Unfolding rate models}
3811 \label{sec.k_models}
3813 We have two main choices for determining $k$: Bell's model, which assumes the
3814 folded domains exist in a single `folded' state and make instantaneous,
3815 stochastic transitions over a free energy barrier; and Kramers' model, which
3816 treats the unfolding as a thermalized diffusion process.
3817 We deal with these options like we dealt with the different tension models: by bundling all model-specific
3818 parameters into structures, with model specific functions for determining $k$.
3820 <<k func definition>>=
3821 typedef double k_func_t(double F, environment_t *env, void *params);
3824 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3825 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]].
3827 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
3828 because \LaTeX\ doesn't like underscores outside math mode.
3829 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
3830 You could use spaces instead (`\verb+ +'),
3831 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
3832 which I don't like as much.
3836 <<k func definition>>
3837 <<null k declarations>>
3838 <<const k declarations>>
3839 <<bell k declarations>>
3840 <<kramers k declarations>>
3841 <<kramers integ k declarations>>
3844 <<k model module makefile lines>>=
3845 k_model.c : sawsim.nw
3846 notangle -Rk-model.c $^ > $@
3847 k_model.h : sawsim.nw
3848 notangle -Rk-model.h $^ > $@
3849 check_k_model.c : sawsim.nw
3850 notangle -Rcheck-k-model.c $^ > $@
3851 check_k_model : check_k_model.c global.h k_model.c k_model.h
3852 gcc -g -o $@ $< k_model.c -lcheck -lgsl -lgslcblas -lm
3853 clean_k_model : clean_k_model_utils
3854 rm -f k_model.c k_model.h check_k_model.c check_k_model
3855 k_model_utils.c : sawsim.nw
3856 notangle -Rk-model-utils.c $^ > $@
3857 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
3858 gcc -g -o $@ $< k_model.c parse.c string_eval.c -lgsl -lgslcblas -lm
3859 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
3860 gcc -g -static -o $@ $< k_model.c parse.c string_eval.c -lgsl -lgslcblas -lm
3861 clean_k_model_utils :
3862 rm -f k_model_utils.c k_model_utils
3868 For domains that are never folded.
3870 <<null k declarations>>=
3872 <<null k functions>>=
3877 <<null k model getopt>>=
3878 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
3881 \subsection{Constant rate model}
3884 This is a pretty boring model, but it is usefull for testing the engine.
3888 where $k_0$ is the constant unfolding rate.
3890 <<const k functions>>=
3891 <<const k function>>
3892 <<const k structure create/destroy functions>>
3895 <<const k declarations>>=
3896 <<const k function declaration>>
3897 <<const k structure create/destroy function declarations>>
3898 <<const k global declarations>>
3902 <<const k structure>>=
3903 typedef struct const_k_param_struct {
3904 double knot; /* transition rate k_0 (frac population per s) */
3908 <<const k function declaration>>=
3909 double const_k(double F, environment_t *env, void *const_k_params);
3911 <<const k function>>=
3912 double const_k(double F, environment_t *env, void *const_k_params)
3913 { /* Returns the rate constant k in frac pop per s. */
3914 const_k_param_t *p = (const_k_param_t *) const_k_params;
3916 assert(p->knot > 0);
3921 <<const k structure create/destroy function declarations>>=
3922 void *string_create_const_k_param_t(char **param_strings);
3923 void destroy_const_k_param_t(void *p);
3926 <<const k structure create/destroy functions>>=
3927 const_k_param_t *create_const_k_param_t(double knot)
3929 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
3930 assert(ret != NULL);
3935 void *string_create_const_k_param_t(char **param_strings)
3937 return create_const_k_param_t(atof(param_strings[0]));
3940 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
3947 <<const k global declarations>>=
3948 extern int num_const_k_params;
3949 extern char *const_k_param_descriptions[];
3950 extern char const_k_param_string[];
3953 <<const k globals>>=
3954 int num_const_k_params = 1;
3955 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
3956 char const_k_param_string[]="1";
3959 <<const k model getopt>>=
3960 {"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}
3963 \subsection{Bell's model}
3966 Quantitatively, Bell's model gives $k$ as
3968 k = k_0 \cdot e^{F dx / k_B T} \;,
3970 where $k_0$ is the unforced unfolding rate,
3971 $F$ is the applied unfolding force,
3972 $dx$ is the distance to the transition state, and
3973 $k_B T$ is the thermal energy\citep{schlierf06}.
3975 <<bell k functions>>=
3977 <<bell k structure create/destroy functions>>
3980 <<bell k declarations>>=
3981 <<bell k function declaration>>
3982 <<bell k structure create/destroy function declarations>>
3983 <<bell k global declarations>>
3987 <<bell k structure>>=
3988 typedef struct bell_param_struct {
3989 double knot; /* transition rate k_0 (frac population per s) */
3990 double dx; /* `distance to transition state' (nm) */
3994 <<bell k function declaration>>=
3995 double bell_k(double F, environment_t *env, void *bell_params);
3997 <<bell k function>>=
3998 double bell_k(double F, environment_t *env, void *bell_params)
3999 { /* Returns the rate constant k in frac pop per s.
4000 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4001 * uses global k_B in J/K */
4002 bell_param_t *p = (bell_param_t *) bell_params;
4003 assert(F >= 0); assert(env->T > 0);
4005 assert(p->knot > 0); assert(p->dx >= 0);
4007 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
4008 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
4009 p->knot * exp(F*p->dx / (k_B*env->T) ));
4010 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
4012 return p->knot * exp(F*p->dx / (k_B*env->T) );
4016 <<bell k structure create/destroy function declarations>>=
4017 void *string_create_bell_param_t(char **param_strings);
4018 void destroy_bell_param_t(void *p);
4021 <<bell k structure create/destroy functions>>=
4022 bell_param_t *create_bell_param_t(double knot, double dx)
4024 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
4025 assert(ret != NULL);
4031 void *string_create_bell_param_t(char **param_strings)
4033 return create_bell_param_t(atof(param_strings[0]), atof(param_strings[1]));
4036 void destroy_bell_param_t(void *p /* bell_param_t * */)
4043 <<bell k global declarations>>=
4044 extern int num_bell_params;
4045 extern char *bell_param_descriptions[];
4046 extern char bell_param_string[];
4050 int num_bell_params = 2;
4051 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
4052 char bell_param_string[]="3.3e-4,0.25e-9";
4055 <<bell k model getopt>>=
4056 {"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}
4058 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4059 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4061 \subsection{Kramer's model}
4064 Kramer's model gives $k$ as
4066 % \frac{1}{k} = \frac{1}{D}
4067 % \int_{x_\text{min}}^{x_\text{max}}
4068 % e^{\frac{-U_F(x)}{k_B T}}
4069 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4072 %where $D$ is the diffusion constant,
4073 %$U_F(x)$ is the free energy along the unfolding coordinate, and
4074 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4075 % before the minimum and after the maximum, respectively \citep{schlierf06}.
4077 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
4079 where $D$ is the diffusion constant,
4081 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
4083 are length scales where
4084 $x_c(F)$ is the location of the bound state and
4085 $x_{ts}(F)$ is the location of the transition state,
4086 $E(F,x)$ is the free energy, and
4087 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
4088 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c,
4089 \citet{evans97} Eqn.~3).
4091 <<kramers k functions>>=
4092 <<kramers k function>>
4093 <<kramers k structure create/destroy functions>>
4096 <<kramers k declarations>>=
4097 <<kramers k function declaration>>
4098 <<kramers k structure create/destroy function declarations>>
4099 <<kramers k global declarations>>
4103 <<kramers k structure>>=
4104 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
4105 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
4107 typedef struct kramers_param_struct {
4108 double D; /* diffusion rate (um/s) */
4115 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
4116 //kramers_x_func_t *xts; /* function returning transition x (nm) */
4117 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
4118 //kramers_E_func_t *E; /* function returning E (J) */
4119 //double *E_params; /* parametrize the energy landscape */
4120 //int n_E_params; /* length of E_params array */
4124 <<kramers k function declaration>>=
4125 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
4126 extern double kramers_k(double F, environment_t *env, void *kramers_params);
4128 <<kramers k function>>=
4129 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
4131 static char input[160]={0};
4132 static char output[80]={0};
4133 /* setup the environment */
4134 fprintf(in, "F = %g; T = %g\n", F, T);
4135 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4136 string_eval(in, out, input, 80, output);
4137 return atof(output);
4140 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
4142 static char input[160]={0};
4143 static char output[80]={0};
4144 /* setup the environment */
4145 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
4146 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
4147 string_eval(in, out, input, 80, output);
4148 return atof(output);
4151 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
4153 kramers_param_t *p = (kramers_param_t *) kramers_params;
4154 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
4157 double kramers_k(double F, environment_t *env, void *kramers_params)
4158 { /* Returns the rate constant k in frac pop per s.
4159 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
4160 * uses global k_B in J/K */
4161 kramers_param_t *p = (kramers_param_t *) kramers_params;
4162 double kbT, xc, xts, lc, lts, Eb;
4163 assert(F >= 0); assert(env->T > 0);
4166 assert(p->xc != NULL); assert(p->xts != NULL);
4167 assert(p->ddEdxx != NULL); assert(p->E != NULL);
4168 assert(p->in != NULL); assert(p->out != NULL);
4170 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
4171 if (xc == -1.0) return -1.0; /* error (Too much force) */
4172 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
4173 if (xc == -1.0) return -1.0; /* error (Too much force) */
4174 lc = sqrt(2.0*M_PI*kbT /
4175 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
4176 lts = sqrt(-2.0*M_PI*kbT/
4177 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
4178 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
4179 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
4180 //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));
4181 return p->D/(lc*lts) * exp(-Eb/kbT);
4185 <<kramers k structure create/destroy function declarations>>=
4186 void *string_create_kramers_param_t(char **param_strings);
4187 void destroy_kramers_param_t(void *p);
4190 <<kramers k structure create/destroy functions>>=
4191 kramers_param_t *create_kramers_param_t(double D,
4192 char *xc_fn, char *xts_fn,
4193 char *E_fn, char *ddEdxx_fn,
4194 char *global_define_string)
4195 // kramers_x_func_t *xc, kramers_x_func_t *xts,
4196 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
4197 // double *E_params, int n_E_params)
4199 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
4200 assert(ret != NULL);
4201 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4202 assert(ret->xc != NULL);
4203 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
4204 assert(ret->xts != NULL);
4205 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
4206 assert(ret->E != NULL);
4207 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
4208 assert(ret->ddEdxx != NULL);
4210 strcpy(ret->xc, xc_fn);
4211 strcpy(ret->xts, xts_fn);
4212 strcpy(ret->E, E_fn);
4213 strcpy(ret->ddEdxx, ddEdxx_fn);
4214 //ret->E_params = E_params;
4215 //ret->n_E_params = n_E_params;
4216 string_eval_setup(&ret->in, &ret->out);
4217 fprintf(ret->in, "kB = %g\n", k_B);
4218 fprintf(ret->in, "%s\n", global_define_string);
4222 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
4223 void *string_create_kramers_param_t(char **param_strings)
4225 return create_kramers_param_t(atof(param_strings[0]),
4233 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
4235 kramers_param_t *pk = (kramers_param_t *)p;
4237 string_eval_teardown(&pk->in, &pk->out);
4239 // free(pk->E_params);
4245 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
4246 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.
4247 We follow \citet{shillcock98} and use
4249 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
4250 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
4253 where TODO, variable meanings.%\citep{shillcock98}.
4254 The first and second derivatives of $E(F,E_b,x,x_s)$ are
4256 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\\
4257 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
4259 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
4260 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
4261 The roots are therefor at
4263 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
4264 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
4265 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
4268 <<kramers k global declarations>>=
4269 extern int num_kramers_params;
4270 extern char *kramers_param_descriptions[];
4271 extern char kramers_param_string[];
4274 <<kramers k globals>>=
4275 int num_kramers_params = 6;
4276 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"};
4277 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)}";
4279 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
4280 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
4281 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
4282 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.
4283 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
4284 It works so long as [[val_1]] is not [[false]].
4286 <<kramers k model getopt>>=
4287 {"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}
4290 \subsection{Kramer's model (natural cubic splines)}
4291 \label{sec.kramers_integ}
4293 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
4294 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
4295 \citet{schlierf06} Eqn.~2)
4297 \frac{1}{k} = \frac{1}{D}
4298 \int_{x_\text{min}}^{x_\text{max}}
4299 e^{\frac{U_F(x)}{k_B T}}
4300 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4303 where $D$ is the diffusion constant,
4304 $U_F(x)$ is the free energy along the unfolding coordinate, and
4305 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
4306 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
4308 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
4309 interpolating between them with cubic splines.
4310 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
4312 <<kramers integ k functions>>=
4313 <<spline functions>>
4314 <<kramers integ A k functions>>
4315 <<kramers integ B k functions>>
4316 <<kramers integ k function>>
4317 <<kramers integ k structure create/destroy functions>>
4320 <<kramers integ k declarations>>=
4321 <<kramers integ k function declaration>>
4322 <<kramers integ k structure create/destroy function declarations>>
4323 <<kramers integ k global declarations>>
4327 <<kramers integ k structure>>=
4328 typedef double func_t(double x, void *params);
4329 typedef struct kramers_integ_param_struct {
4330 double D; /* diffusion rate (um/s) */
4331 double x_min; /* integration bounds */
4333 func_t *E; /* function returning E (J) */
4334 void *E_params; /* parametrize the energy landscape */
4335 destroy_data_func_t *destroy_E_params;
4337 double F; /* for passing into GSL evaluations */
4339 } kramers_integ_param_t;
4342 <<spline param structure>>=
4343 typedef struct spline_param_struct {
4344 int n_knots; /* natural cubic spline knots for U(x) */
4347 gsl_spline *spline; /* GSL spline parameters */
4348 gsl_interp_accel *acc; /* GSL spline acceleration data */
4352 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
4354 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
4356 <<kramers integ A k functions>>=
4357 double kramers_integ_k_integrandA(double x, void *params)
4359 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4360 double U = (*p->E)(x, p->E_params);
4361 double Fx = p->F * x;
4362 double kbT = k_B * p->env->T;
4363 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
4364 return exp(-(U-Fx)/kbT);
4366 double kramers_integ_k_integralA(double x, void *params)
4368 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4370 double result, abserr;
4371 size_t neval; /* for qng */
4372 static gsl_integration_workspace *w = NULL;
4373 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4374 f.function = &kramers_integ_k_integrandA;
4376 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4377 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4378 //fprintf(stderr, "integralA = %g\n", result);
4384 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
4385 \int_{x_\text{min}}^{x_\text{max}}
4386 e^{\frac{U_F(x)}{k_B T}}
4387 \text{Integral}_A(x)
4390 <<kramers integ B k functions>>=
4391 double kramers_integ_k_integrandB(double x, void *params)
4393 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4394 double U = (*p->E)(x, p->E_params);
4395 double Fx = p->F * x;
4396 double kbT = k_B * p->env->T;
4397 double intA = kramers_integ_k_integralA(x, params);
4398 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
4399 return exp((U-Fx)/kbT)*intA;
4401 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
4403 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4405 double result, abserr;
4406 size_t neval; /* for qng */
4407 static gsl_integration_workspace *w = NULL;
4408 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
4409 f.function = &kramers_integ_k_integrandB;
4413 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
4414 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
4415 //fprintf(stderr, "integralB = %g\n", result);
4420 With the integrals taken care of, evaluating $k$ is simply
4422 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
4424 <<kramers integ k function declaration>>=
4425 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
4427 <<kramers integ k function>>=
4428 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
4429 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
4430 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
4431 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
4435 <<kramers integ k structure create/destroy function declarations>>=
4436 void *string_create_kramers_integ_param_t(char **param_strings);
4437 void destroy_kramers_integ_param_t(void *p);
4440 <<kramers integ k structure create/destroy functions>>=
4441 kramers_integ_param_t *create_kramers_integ_param_t(double D,
4442 double xmin, double xmax, func_t *E,
4444 destroy_data_func_t *destroy_E_params)
4446 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
4447 assert(ret != NULL);
4452 ret->E_params = E_params;
4453 ret->destroy_E_params = destroy_E_params;
4457 void *string_create_kramers_integ_param_t(char **param_strings)
4459 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
4460 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
4461 void *E_params = string_create_spline_param_t(param_strings+1);
4462 return create_kramers_integ_param_t(atof(param_strings[0]),
4463 atof(param_strings[2]),
4464 atof(param_strings[3]),
4465 &spline_eval, E_params,
4466 destroy_spline_param_t);
4469 void destroy_kramers_integ_param_t(void *params)
4471 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
4474 (*p->destroy_E_params)(p->E_params);
4480 Finally we have the GSL spline wrappers:
4481 <<spline functions>>=
4483 <<create/destroy spline>>
4486 <<spline function>>=
4487 double spline_eval(double x, void *spline_params)
4489 spline_param_t *p = (spline_param_t *)(spline_params);
4490 return gsl_spline_eval(p->spline, x, p->acc);
4494 <<create/destroy spline>>=
4495 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
4497 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
4498 assert(ret != NULL);
4499 ret->n_knots = n_knots;
4502 ret->acc = gsl_interp_accel_alloc();
4503 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
4504 gsl_spline_init(ret->spline, x, y, n_knots);
4508 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
4509 void *string_create_spline_param_t(char **param_strings)
4513 double *x=NULL, *y=NULL;
4514 /* break into ordered pairs */
4515 parse_list_string(param_strings[0], SEP, '(', ')',
4516 &num_knots, &knot_coords);
4517 x = (double *)malloc(sizeof(double)*num_knots);
4519 y = (double *)malloc(sizeof(double)*num_knots);
4521 for (i=0; i<num_knots; i++) {
4522 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
4523 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
4526 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
4528 free_parsed_list(num_knots, knot_coords);
4529 return create_spline_param_t(num_knots, x, y);
4532 void destroy_spline_param_t(void *params)
4534 spline_param_t *p = (spline_param_t *)params;
4537 gsl_spline_free(p->spline);
4539 gsl_interp_accel_free(p->acc);
4549 <<kramers integ k global declarations>>=
4550 extern int num_kramers_integ_params;
4551 extern char *kramers_integ_param_descriptions[];
4552 extern char kramers_integ_param_string[];
4555 <<kramers integ k globals>>=
4556 int num_kramers_integ_params = 4;
4557 char *kramers_integ_param_descriptions[] = {
4558 "diffusion D, m^2 / s",
4559 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
4560 "starting integration bound x_min, m",
4561 "ending integration bound x_max, m"};
4562 //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";
4563 //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";
4564 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";
4565 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
4566 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
4567 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
4570 <<kramers integ k model getopt>>=
4571 {"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}
4573 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4574 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4576 \subsection{Unfolding model implementation}
4580 <<k model includes>>
4581 #include "k_model.h"
4582 <<k model internal definitions>>
4584 <<k model functions>>
4587 <<k model includes>>=
4588 #include <assert.h> /* assert() */
4589 #include <stdlib.h> /* NULL, malloc() */
4590 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4591 #include <stdio.h> /* fprintf(), stdout */
4592 #include <string.h> /* strlen(), strcpy() */
4593 #include <gsl/gsl_integration.h>
4594 #include <gsl/gsl_spline.h>
4599 <<k model internal definitions>>=
4600 <<const k structure>>
4601 <<bell k structure>>
4602 <<kramers k structure>>
4603 <<kramers integ k structure>>
4604 <<spline param structure>>
4607 <<k model globals>>=
4611 <<kramers k globals>>
4612 <<kramers integ k globals>>
4615 <<k model functions>>=
4616 <<null k functions>>
4617 <<const k functions>>
4618 <<bell k functions>>
4619 <<kramers k functions>>
4620 <<kramers integ k functions>>
4623 \subsection{Unfolding model unit tests}
4625 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
4626 <<check-k-model.c>>=
4627 <<k model check includes>>
4628 <<check relative error>>
4630 <<k model test suite>>
4631 <<main check program>>
4634 <<k model check includes>>=
4635 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
4636 #include <stdio.h> /* sprintf() */
4637 #include <assert.h> /* assert() */
4638 #include <math.h> /* exp() */
4641 #include "k_model.h"
4644 <<k model test suite>>=
4647 <<k model suite definition>>
4650 <<k model suite definition>>=
4651 Suite *test_suite (void)
4653 Suite *s = suite_create ("k model");
4654 <<const k test case defs>>
4655 <<bell k test case defs>>
4657 <<const k test case adds>>
4658 <<bell k test case adds>>
4663 \subsubsection{Constant}
4665 <<const k test case defs>>=
4666 TCase *tc_const_k = tcase_create("Constant k");
4669 <<const k test case adds>>=
4670 tcase_add_test(tc_const_k, test_const_k_create_destroy);
4671 tcase_add_test(tc_const_k, test_const_k_over_range);
4672 suite_add_tcase(s, tc_const_k);
4677 START_TEST(test_const_k_create_destroy)
4679 char *knot[] = {"1","2","3","4","5","6"};
4680 char *params[] = {knot[0]};
4683 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4684 params[0] = knot[i];
4685 p = string_create_const_param_t(params);
4686 destroy_const_param_t(p);
4691 START_TEST(test_const_k_over_range)
4693 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
4694 char *knot[] = {"1","2","3","4","5","6"};
4695 char *params[] = {knot[0]};
4698 char param_string[80];
4701 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4702 params[0] = knot[i];
4703 p = string_create_const_param_t(params);
4704 for ( F=Fm, F<FM, F+=dF ) {
4705 fail_unless(const_k(F, &env, p)==atof(knot[i]),
4706 "const_k(%g, %g, &{%s,%s}) = %g != %s",
4707 F, T, knot[i], dx, const_k(F, &env, p), knot[i]);
4709 destroy_const_param_t(p);
4715 \subsubsection{Bell}
4717 <<bell k test case defs>>=
4718 TCase *tc_bell = tcase_create("Bell's k");
4721 <<bell k test case adds>>=
4722 tcase_add_test(tc_bell, test_bell_k_create_destroy);
4723 tcase_add_test(tc_bell, test_bell_k_at_zero);
4724 tcase_add_test(tc_bell, test_bell_k_at_one);
4725 suite_add_tcase(s, tc_bell);
4729 START_TEST(test_bell_k_create_destroy)
4731 char *knot[] = {"1","2","3","4","5","6"};
4733 char *params[] = {knot[0], dx};
4736 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4737 params[0] = knot[i];
4738 p = string_create_bell_param_t(params);
4739 destroy_bell_param_t(p);
4744 START_TEST(test_bell_k_at_zero)
4746 double F=0.0, T=300.0;
4747 char *knot[] = {"1","2","3","4","5","6"};
4749 char *params[] = {knot[0], dx};
4752 char param_string[80];
4755 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4756 params[0] = knot[i];
4757 p = string_create_bell_param_t(params);
4758 fail_unless(bell_k(F, &env, p)==atof(knot[i]),
4759 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
4760 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
4761 destroy_bell_param_t(p);
4766 START_TEST(test_bell_k_at_one)
4769 char *knot[] = {"1","2","3","4","5","6"};
4771 char *params[] = {knot[0], dx};
4772 double F= k_B*T/atof(dx);
4775 char param_string[80];
4778 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4779 params[0] = knot[i];
4780 p = string_create_bell_param_t(params);
4781 CHECK_ERR(1e-6, atof(knot[i])*exp(1.0), bell_k(F, &env, p));
4782 //fail_unless(bell_k(F, &env, p)==atof(knot[i])*exp(1.0),
4783 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
4784 // F, T, knot[i], dx, bell_k(F, &env, p), atof(knot[i])*exp(1.0));
4785 destroy_bell_param_t(p);
4791 <<kramers k tests>>=
4794 <<kramers k test case def>>=
4797 <<kramers k test case add>>=
4800 <<k model function tests>>=
4801 <<check relative error>>
4803 START_TEST(test_something)
4805 double k=5, x=3, last_x=2.0, F;
4806 one_dim_fn_t *handlers[] = {&hooke};
4807 void *data[] = {&k};
4808 double xi[] = {0.0};
4810 int new_active_groups = 1;
4811 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
4812 fail_unless(F = k*x, NULL);
4817 \subsection{Utilities}
4819 The unfolding models can get complicated, and are sometimes hard to explain to others.
4820 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
4821 the overhead of having to understand a full simulation.
4822 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
4823 or special mode, where the operation depends on the specific model selected.
4825 <<k-model-utils.c>>=
4827 <<k model utility includes>>
4828 <<k model utility definitions>>
4829 <<k model utility getopt functions>>
4830 <<k model utility multi model E>>
4831 <<k model utility main>>
4834 <<k model utility main>>=
4835 int main(int argc, char **argv)
4837 <<k model getopt array definition>>
4838 k_model_getopt_t *model = NULL;
4843 int num_param_args; /* for INIT_MODEL() */
4844 char **param_args; /* for INIT_MODEL() */
4846 double special_xmin;
4847 double special_xmax;
4848 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
4849 &Fmax, &special_xmin, &special_xmax, &flags);
4850 if (flags & VFLAG) {
4851 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4853 INIT_MODEL("utility", model, model->params, params);
4856 if (model->k == NULL) {
4857 printf("No k function for model %s\n", model->name);
4861 printf("Fmax = %g <= 0\n", Fmax);
4867 printf("#F (N)\tk (%% pop. per s)\n");
4868 for (i=0; i<=N; i++) {
4869 F = Fmax*i/(double)N;
4870 k = (*model->k)(F, &env, params);
4872 printf("%g\t%g\n", F, k);
4877 if (model->k == NULL || model->k == &bell_k) {
4878 printf("No special mode for model %s\n", model->name);
4881 if (special_xmax <= special_xmin) {
4882 printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
4886 double Xrng=(special_xmax-special_xmin), x, E;
4888 printf("#x (m)\tE (J)\n");
4889 for (i=0; i<=N; i++) {
4890 x = special_xmin + Xrng*i/(double)N;
4891 E = multi_model_E(model->k, params, &env, x);
4892 printf("%g\t%g\n", x, E);
4899 if (model->destructor)
4900 (*model->destructor)(params);
4905 <<k model utility multi model E>>=
4906 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
4908 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
4910 if (k == kramers_integ_k)
4911 E = (*p->E)(x, p->E_params);
4913 E = kramers_E(0, env, model_params, x);
4919 <<k model utility includes>>=
4922 #include <string.h> /* strlen() */
4923 #include <assert.h> /* assert() */
4926 #include "string_eval.h"
4927 #include "k_model.h"
4930 <<k model utility definitions>>=
4931 <<version definition>>
4932 #define VFLAG 1 /* verbose */
4933 enum mode_t {M_K_OF_F, M_SPECIAL};
4934 <<string comparison definition and globals>>
4935 <<k model getopt definitions>>
4936 <<initialize model definition>>
4937 <<kramers k structure>>
4938 <<kramers integ k structure>>
4942 <<k model utility getopt functions>>=
4944 <<k model utility help>>
4945 <<k model utility get options>>
4948 <<k model utility help>>=
4949 void help(char *prog_name,
4951 int n_k_models, k_model_getopt_t *k_models,
4955 printf("usage: %s [options]\n", prog_name);
4956 printf("Version %s\n\n", VERSION);
4957 printf("A utility for understanding the available unfolding models\n\n");
4958 printf("Environmental options:\n");
4959 printf("-T T\tTemperature T (currently %g K)\n", env->T);
4960 printf("-C T\tYou can also set the temperature T in Celsius\n");
4961 printf("Model options:\n");
4962 printf("-k model\tTransition rate model (currently %s)\n",
4963 k_models[k_model].name);
4964 printf("-K args\tTransition rate model argument string (currently %s)\n",
4965 k_models[k_model].params);
4966 printf("There are two output modes. In standard mode, k(F) is printed\n");
4967 printf("For example:\n");
4968 printf(" #Force (N)\tk (% pop. per s)\n");
4969 printf(" 123.456\t7.89\n");
4971 printf("In special mode, the output depends on the model.\n");
4972 printf("For models defining an energy function E(x), that function is printed\n");
4973 printf(" #Position (m)\tE (J)\n");
4974 printf(" 0.0012\t0.0034\n");
4976 printf("-m\tChange output to standard mode\n");
4977 printf("-M\tChange output to special mode\n");
4978 printf("-F\tSet the maximum F value for the standard mode k(F) output\n");
4979 printf("-x\tSet the minimum x value for the special mode E(x) output\n");
4980 printf("-X\tSet the maximum x value for the special mode E(x) output\n");
4981 printf("-V\tChange output to verbose mode\n");
4982 printf("-h\tPrint this help and exit\n");
4984 printf("Unfolding rate models:\n");
4985 for (i=0; i<n_k_models; i++) {
4986 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
4987 for (j=0; j < k_models[i].num_params ; j++)
4988 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
4989 printf("\t\tdefaults: %s\n", k_models[i].params);
4994 <<k model utility get options>>=
4995 void get_options(int argc, char **argv, environment_t *env,
4996 int n_k_models, k_model_getopt_t *k_models,
4997 enum mode_t *mode, k_model_getopt_t **model,
4998 double *Fmax, double *special_xmin, double *special_xmax,
4999 unsigned int *flags)
5001 char *prog_name = NULL;
5002 char c, options[] = "T:C:k:K:mMF:x:X:Vh";
5004 extern char *optarg;
5005 extern int optind, optopt, opterr;
5009 /* setup defaults */
5011 prog_name = argv[0];
5012 env->T = 300.0; /* K */
5017 *special_xmax = 1e-8;
5018 *special_xmin = 0.0;
5020 while ((c=getopt(argc, argv, options)) != -1) {
5022 case 'T': env->T = atof(optarg); break;
5023 case 'C': env->T = atof(optarg)+273.15; break;
5025 k_model = index_k_model(n_k_models, k_models, optarg);
5026 *model = k_models+k_model;
5029 k_models[k_model].params = optarg;
5031 case 'm': *mode = M_K_OF_F; break;
5032 case 'M': *mode = M_SPECIAL; break;
5033 case 'F': *Fmax = atof(optarg); break;
5034 case 'x': *special_xmin = atof(optarg); break;
5035 case 'X': *special_xmax = atof(optarg); break;
5036 case 'V': *flags |= VFLAG; break;
5038 fprintf(stderr, "unrecognized option '%c'\n", optopt);
5039 /* fall through to default case */
5041 help(prog_name, env, n_k_models, k_models, k_model);
5050 \section{\LaTeX\ documentation}
5052 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
5053 The comment blocks are extracted (with nicely formatted code blocks), using
5054 <<latex makefile lines>>=
5055 sawsim.tex : sawsim.nw
5056 noweave -latex -index -delay $^ > $@
5059 <<latex makefile lines>>=
5060 sawsim.pdf : sawsim.tex sawsim.bib
5066 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
5067 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
5068 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
5070 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
5071 <<latex makefile lines>>=
5073 rm -f sawsim.aux sawsim.bbl sawsim.blg sawsim.log sawsim.out sawsim.tex
5074 clean_latex : semi-clean_latex
5080 [[make]] is a common utility on *nix systems for managing dependencies while building software.
5081 In this case, we have several \emph{targets} that we'd like to build:
5082 the main simulation program \prog;
5083 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
5084 the documentation [[sawsim.pdf]];
5085 and the [[Makefile]] itself.
5086 Besides the generated files, there is the utility target
5087 [[clean]] for removing all generated files except [[Makefile]].
5089 % [[dist]] for generating a distributable tar file.
5091 Extract the makefile with
5092 `[[notangle -Rmakefile sawsim.nw | sed 's/ /\t/' > Makefile]]'.
5093 The sed is needed because notangle replaces tabs with eight spaces,
5094 but [[make]] needs tabs.
5096 all : sawsim sawsim.pdf Makefile tension_model_utils k_model_utils
5100 profile : sawsim_profile
5101 sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ -mnull -Mwlc -A0.39e-9,28e-9 -F8
5102 gprof sawsim_profile gmon.out > $@
5104 clean : clean_check clean_list clean_tension clean_tension_model clean_tension_model_utils clean_k_model clean_parse clean_string_eval clean_latex
5105 rm -f sawsim.c sawsim sawsim_static sawsim_profile gmon.out global.h
5107 sawsim.c : sawsim.nw
5109 global.h : sawsim.nw
5110 notangle -Rglobal.h $^ > $@
5111 sawsim : sawsim.c global.h list.c list.h tension_balance.c tension_balance.h \
5112 tension_model.c tension_model.h k_model.c k_model.h parse.c parse.h string_eval.c string_eval.h \
5114 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
5115 sawsim_static : sawsim
5116 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
5117 sawsim_profile : sawsim
5118 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
5121 check_sawsim.c : sawsim.nw
5122 notangle -Rchecks $^ > $@
5123 check_sawsim : check_sawsim.c global.h list.c list.h tension_balance.c tension_balance.h \
5124 k_model.c k_model.h parse.c parse.h string_eval.c string_eval.h accel_k.c accel_k.h
5125 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
5127 rm -f check_sawsim check_sawsim.c
5128 check : check_list check_tension_balance check_k_model check_parse check_string_eval check_accel_k check_sawsim
5130 check_tension_balance
5137 <<list module makefile lines>>
5138 <<tension balance module makefile lines>>
5139 <<tension model module makefile lines>>
5140 <<k model module makefile lines>>
5141 <<parse module makefile lines>>
5142 <<string eval module makefile lines>>
5143 <<accel k module makefile lines>>
5144 <<latex makefile lines>>
5146 Makefile : sawsim.nw
5147 notangle -Rmakefile $^ | sed 's/ /\t/' > $@
5149 This is fairly self-explanatory.
5150 We compile a dynamically linked version ([[sawsim]]) and a statically linked version ([[sawsim_static]]).
5151 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.
5156 \subsection{Hookean springs in series}
5157 \label{app.math_hooke}
5160 The effective spring constant for $n$ springs $k_i$ in series is given by
5162 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5168 F &= k_i x_i &&\forall i \le n \\
5169 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
5170 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
5171 F &= k_1 x_1 = k_\text{eff} x \\
5172 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
5173 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
5178 \addcontentsline{toc}{section}{References}
5179 \bibliography{sawsim}