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 process.
77 In the \prog\ program, we use Monte Carlo methods to simulate this unfolding
78 through an explicit time-stepping approach.
79 We develop a program in \lang\ to simulate probability of unfolding under
80 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 \lang\
89 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,
100 and a switch to select Bell's or Kramers' model.
101 This induced a major rewrite, introducing the tension groups abstraction, and a split to multiple \lang\ source files.
102 Also added a unit-test suites for sanity checking, and switched to SI units throughout.
104 Version 0.3 added integrated Kramers' model unfolding rate calculations.
105 Also replaced some of my hand-coded routines with GNU Scientific Library \citepalias{galassi05} calls.
107 Version 0.4 added Python evaluation for the saddle-point Kramers' model unfolding rate calculations,
108 so that the model functions could be controlled from the command line.
109 Also added interpolating lookup tables to accelerate slow unfolding rate model evaluation.
111 <<version definition>>=
112 #define VERSION "0.4"
118 sawsim - program for simulating single molecule mechanical unfolding.
119 Copyright (C) 2008, William Trevor King
121 This program is free software; you can redistribute it and/or
122 modify it under the terms of the GNU General Public License as
123 published by the Free Software Foundation; either version 3 of the
124 License, or (at your option) any later version.
126 This program is distributed in the hope that it will be useful, but
127 WITHOUT ANY WARRANTY; without even the implied warranty of
128 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
129 See the GNU General Public License for more details.
131 You should have received a copy of the GNU General Public License
132 along with this program; if not, write to the Free Software
133 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
136 The author may be contacted at <wking@drexel.edu> on the Internet, or
137 write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
138 Philadelphia PA 19104, USA.
150 The unfolding system is stored as a chain of \emph{domains}.
151 Domains can be folded, globular protein domains, unfolded protein linkers,
152 AFM cantilevers, or other stretchable system components.
154 Each domain contributes to the total tension, which depends on the chain extension.
155 The particular model for the tension as a function of extension is handled generally with function pointers.
156 So far the following models have been implemented:
158 \item Null (Appendix \ref{sec.null_tension}),
159 \item Constant (Appendix \ref{sec.const_tension}),
160 \item Hookean spring (Appendix \ref{sec.hooke}),
161 \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
162 \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
165 The tension model and parameters used for a particular domain depend on whether the domain is folded or unfolded.
166 The transition rate from the folded state to the unfolded state is also handled generally with function pointers, with implementations for
168 \item Null (Appendix \ref{sec.null_k}),
169 \item Bell's model (Appendix \ref{sec.bell}),
170 \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and
171 \item Kramers' saddle point model (Appendix \ref{sec.kramers}).
174 The unfolding of the chain domains is modeled in two stages.
175 First the tension of the chain is calculated.
176 Then the tension (assumed to be constant over the length of the chain),
177 is applied to each folded domain, and used to calculate the probability of
178 that subdomain unfolding in the next timestep $dt$.
179 Then the time is stepped forward, and the process is repeated.
182 int main(int argc, char **argv)
184 <<tension model getopt array definition>>
185 <<k model getopt array definition>>
186 list_t *domain_list=NULL; /* variables initialized in get_options() */
187 environment_t env={0};
190 double x=0, dt, F; /* variables used in the simulation loop */
191 one_dim_func_t *tension_handler[NUM_TENSION_GROUPS] = {0};
192 get_options(argc, argv, &P, &dt_max, &v, &env, NUM_TENSION_GROUPS,
193 tension_models, NUM_K_MODELS, k_models, &domain_list, &num_folded);
194 setup(tension_handler);
195 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Position (m)\tForce (N)\n");
196 if (flags & FLAG_OUTPUT_UNFOLDING_FORCES) printf("#Unfolding Force (N)\n");
197 while (num_folded > 0) {
198 F = find_tension(tension_handler, domain_list, &env, x);
199 dt = determine_dt(tension_handler, domain_list, &env, x, P, dt_max, v);
200 random_unfoldings(domain_list, F, dt, &env, &num_folded);
201 if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\n", x, F);
204 destroy_domain_list(domain_list);
208 @ The meat of the simulation is bundled into the three functions
209 [[find_tension]], [[determine_dt]], and [[random_unfoldings]].
210 [[find_tension]] is discussed in Section \ref{sec.find_tension},
211 [[determine_dt]] in Section \ref{sec.adaptive_dt}, and
212 [[random_unfoldings]] in Section \ref{sec.unfolding_rate}.
214 Environmental parameters important in determining reaction rates and tensions
215 (e.g. temperature, pH) are stored in a single structure
216 to facilitate extension to more complicated models in the future.
217 <<environment definition>>=
218 typedef struct environment_struct {
224 <<environment definition>>
225 <<one dimensional function definition>>
226 <<create/destroy definitions>>
230 \section{Simulation functions}
232 <<simulation functions>>=
236 <<does domain unfold>>
237 <<randomly unfold domains>>
241 \label{sec.find_tension}
243 Because the stretched system may be made up of several parts (folded domains, unfolded domains, spring-like cantilever, \ldots),
244 we will assign the domains to groups. For example, a domain might be in the [[NULL]] group when it's folded and in the worm-like chain group when it is unfolded.
245 The domains are assumed to be commutative, so ordering is ignored.
246 The interactions between the groups are assumed to be linear, but the interactions between domains of the same group need not be.
247 This allows for non-linear group models such as th worm-like or freely-jointed chains.
248 Each handler function receives a list of domains matching it's group number.
249 <<find tension definitions>>=
250 enum tension_group_t {TG_NULL=0,
256 #define NUM_TENSION_GROUPS 5
260 <<tension handler types>>=
261 typedef struct tension_handler_data_struct {
265 } tension_handler_data_t;
269 After sorting the chain into separate groups $G_i$, with tension handlers $F_i(G_i; x_i)$,
270 we need to balance the extension $x_i$ of each group to get a consistent tension $F_i(x_i) = F_j(x_j) \;\; \forall i,j$.
271 For the moment, we will restrict our group boundaries to a single dimension, so $\sum_i x_i = x$, our total extension (it is this restriction that causes the groups to interact linearly).
272 We'll also restrict our extensions to all be positive.
273 With these restrictions, the problem of balancing the tensions reduces to solving a system of $N+1$ possibly non-linear equations with $N$ unknowns,
274 where $N$ is the number of active groups.
275 In general this can be fairly complicated, but without much loss of practicality we can restrict ourselves to strictly monotonically increasing, non-negative tension functions $F_i(x)$, with $F_i(0)=0$, which makes the situation much simpler.
276 For example, it guarantees the existence of a unique, real solution for finite forces.
278 double find_tension(one_dim_func_t **tension_handler, list_t *domain_list, environment_t *env, double x)
280 static tension_handler_data_t data[NUM_TENSION_GROUPS] = {0};
281 static void *pdata[NUM_TENSION_GROUPS] = {0};
282 static double xi[NUM_TENSION_GROUPS] = {0}, last_x = 0;
283 static int active_groups[NUM_TENSION_GROUPS] = {0};
284 int i, new_active_groups=0;
287 for (i=0; i<NUM_TENSION_GROUPS; i++) {
289 data[i].group = get_group_list(domain_list, i);
291 /* data[i].persist continues unchanged */
292 if (data[i].group && i != TG_NULL) {
293 /* an active group */
294 if (active_groups[i] == 0) /* newly active */
295 new_active_groups = 1;
296 active_groups[i] = 1;
298 if (active_groups[i] == 1) /* newly inactive */
299 new_active_groups = 1;
300 active_groups[i] = 0;
303 F = tension_balance(NUM_TENSION_GROUPS, tension_handler, pdata, active_groups, new_active_groups, xi, last_x, x);
308 See Appendix \ref{app.tension_balance} for the tension balancing implementation.
310 Each group must define a parameter structure if the tension model is parametrized,
311 a function to create the parameter structure,
312 a function to destroy the parameter structure, and
313 a tension function of type [[one_dim_func_t]] that receives a [[tension_handler_data_t *]] as it's data argument.
314 The tension-specific parameter structures are contained in the domain group field.
315 See the Section \ref{app.model_selection} for a discussion on the command line framework.
316 See the worm-like chain implementation for an example (Section \ref{sec.wlc}).
318 The major limitation of our approach is that all of our tension models are
319 Markovian (memory-less), which is not really a problem for the chains
320 (see \citet{evans99} for WLC thermalization rate calculations).
322 \subsection{Unfolding rate}
323 \label{sec.unfolding_rate}
325 Each folded domain is modeled as unimolecular, first order reaction
327 \text{Folded} \xrightarrow{k} % in tex: X atop Y
330 With the rate constant $k$ defined by
332 \frac{d[\text{Folded}]}{dt} = -k [\text{Folded}].
334 So the probability of a given protein unfolding in a short time $dt$ is given
340 <<does domain unfold>>=
341 int domain_unfolds(double F, double dt, environment_t *env, domain_t *domain)
342 { /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to a folded domain */
344 k = accel_k(domain->k, F, env, domain->k_params);
345 //(*domain->k)(F, env, domain->k_params);
346 //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
347 return happens(k*dt); /* dice roll for prob. k*dt event */
349 @ [[happens]] is a random decision making function defined in Appendix \ref{app.utils}.
351 Only one domain can unfold in each timestep,
352 because the timescale of a domain unfolding $dt_u$
353 is assumed to be less than the simulation timestep $dt$,
354 so a domain will completely unfold in a single timestep.
355 We adapt our timesteps to keep the probability of a single domain unfolding
356 low, so the probability of two domains unfolding in the same timestep is
358 This approach breaks down as the adaptive timestep scheme approaches
359 $dt \sim dt_u$, but $dt_u \sim 1\U{$\mu$s}$ for Ig27-like proteins
360 \citep{klimov00}, so this shouldn't be a problem.
361 To reassure yourself, you can ask the simulator to print the smallest timestep
364 <<randomly unfold domains>>=
365 void random_unfoldings(list_t *domain_list, double tension,
366 double dt, environment_t *env,
367 int *num_folded_domains)
369 while (domain_list != NULL) {
370 if (D(domain_list)->state == DS_FOLDED
371 && domain_unfolds(tension, dt, env, D(domain_list))) {
372 if (flags & FLAG_OUTPUT_UNFOLDING_FORCES)
373 fprintf(stdout, "%g\n", tension);
374 D(domain_list)->state = DS_UNFOLDED;
375 (*num_folded_domains)--;
376 break; /* our one domain has unfolded, stop looking for others */
378 domain_list = domain_list->next;
383 \subsection{Adaptive timesteps}
384 \label{sec.adaptive_dt}
386 We'd like to pick $dt$ so the probability of unfolding in the next timestep is
387 small. If we don't adapt our timestep, we risk breaking our approximation
388 that $F(x' \in [x, x+v \cdot dt]) = F(x)$
389 or that only one domain unfolds in a given timestep.
390 Because $F(x)$ is monotonically increasing, excessively large timesteps will
391 lead to erroneously large unfolding forces.
392 The simulation would have been accurate for sufficiently small fixed timesteps,
393 but adaptive timesteps will allow us to move through low tension regions in
394 fewer steps, leading to a more efficient simulation.
396 The actual adaptive timestep implementation is not particularly interesting,
397 since we are only required to reduce $dt$ to somewhere below a set threshold,
398 so I've removed it to Appendix \ref{app.adaptive_dt}.
403 The models provide the physics of the simulation,
404 but the simulation allows interchangeable models,
405 and we are currently focusing on the simulation itself,
406 so we remove the actual model implementation to the appendices.
408 The tension models are defined in Appendix \ref{sec.tension_models}
409 and unfolding models are defined in Appendix \ref{sec.k_models}.
412 #define k_B 1.3806503e-23 /* J/K */
416 \section{Command line interface}
418 <<get option functions>>=
420 <<index tension model>>
426 \subsection{Model selection}
427 \label{app.model_selection}
429 The main difficulty with the command line interface in \prog\ is developing an intuitive method for accessing the possibly large number of available models.
430 We'll treat this generally by defining an array of available models, containing their important parameters.
432 <<get options definitions>>=
433 <<model getopt definitions>>
436 <<create/destroy definitions>>=
437 typedef void *create_data_func_t(char **param_strings);
438 typedef void destroy_data_func_t(void *param_struct);
441 <<model getopt definitions>>=
442 <<tension model getopt definitions>>
443 <<k model getopt definitions>>
447 \subsubsection{Tension}
449 <<tension model getopt definitions>>=
450 typedef struct tension_model_getopt_struct {
451 enum tension_group_t tg_group;
455 char **param_descriptions;
457 create_data_func_t *creator;
458 destroy_data_func_t *destructor;
459 } tension_model_getopt_t;
462 <<tension model getopt array definition>>=
463 tension_model_getopt_t tension_models[NUM_TENSION_GROUPS] = {
464 <<null tension model getopt>>,
465 <<constant tension model getopt>>,
466 <<hooke tension model getopt>>,
467 <<worm-like chain tension model getopt>>,
468 <<freely-jointed chain tension model getopt>>
472 \subsubsection{Unfolding rate}
474 <<k model getopt definitions>>=
475 #define NUM_K_MODELS 5
477 typedef struct k_model_getopt_struct {
482 char **param_descriptions;
484 create_data_func_t *creator;
485 destroy_data_func_t *destructor;
489 <<k model getopt array definition>>=
490 k_model_getopt_t k_models[NUM_K_MODELS] = {
491 <<null k model getopt>>,
492 <<const k model getopt>>,
493 <<bell k model getopt>>,
494 <<kramers k model getopt>>,
495 <<kramers integ k model getopt>>
502 void help(char *prog_name, double P, double dt_max, double v,
504 int n_tension_models, tension_model_getopt_t *tension_models,
505 int folded_tension_model, int unfolded_tension_model,
506 int n_k_models, k_model_getopt_t *k_models,
510 printf("usage: %s [options]\n", prog_name);
511 printf("Version %s\n\n", VERSION);
512 printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n");
513 printf("Simulation options:\n");
514 printf("-P P\tTarget probability for dt (currently %g)\n", P);
515 printf("-t dt\tMaximum allowed timestep dt (currently %g)\n", dt_max);
516 printf("-v v\tPulling velocity v (currently %g nm/s)\n", v);
517 printf("Environmental options:\n");
518 printf("-T T\tTemperature T (currently %g K)\n", env->T);
519 printf("-C T\tYou can also set the temperature T in Celsius\n");
520 printf("Model options:\n");
521 printf("The domains exist in either the folded or unfolded state\n");
522 printf("The following options change the default behavior in each state.\n");
523 printf("-m model\tFolded tension model (currently %s)\n",
524 tension_models[folded_tension_model].name);
525 printf("-a args\tFolded tension model argument string (currently %s)\n",
526 tension_models[folded_tension_model].params);
527 printf("-M model\tUnfolded tension model (currently %s)\n",
528 tension_models[unfolded_tension_model].name);
529 printf("-A args\tUnfolded tension model argument string (currently %s)\n",
530 tension_models[unfolded_tension_model].params);
531 printf("The following options change the unfolding rate.\n");
532 printf("-k model\tTransition rate model (currently %s)\n",
533 k_models[k_model].name);
534 printf("-K args\tTransition rate model argument string (currently %s)\n",
535 k_models[k_model].params);
536 printf("Domain creation options:\n");
537 printf("Once you've set up the appropriate default models, you need to add the domains\n");
538 printf("-F n\tAdd n folded domains with the current models\n");
539 printf("-U n\tAdd n unfolded domains with the current models\n");
540 printf("Output mode options:\n");
541 printf("There are two output modes. In standard mode, only the unfolding\n");
542 printf("events are printed. For example:\n");
543 printf(" #Unfolding Force (N)\n");
544 printf(" 123.456e-12\n");
546 printf("In verbose mode, the entire Force-vs-distance curve is output:\n");
547 printf(" #Position (m)\tForce (N)\n");
548 printf(" 0.001\t0.0023\n");
550 printf("-V\tChange output to verbose mode\n");
551 printf("-h\tPrint this help and exit\n");
553 printf("Tension models:\n");
554 for (i=0; i<n_tension_models; i++) {
555 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
556 for (j=0; j < tension_models[i].num_params ; j++)
557 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
558 printf("\t\tdefaults: %s\n", tension_models[i].params);
560 printf("Unfolding rate models:\n");
561 for (i=0; i<n_k_models; i++) {
562 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
563 for (j=0; j < k_models[i].num_params ; j++)
564 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
565 printf("\t\tdefaults: %s\n", k_models[i].params);
567 printf("Example usage: (1 spring and 8 folded domains)\n");
568 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);
572 \subsection{Get options}
575 void get_options(int argc, char **argv,
576 double *pP, double *pDt_max, double *pV,
578 int n_tension_models, tension_model_getopt_t *tension_models,
579 int n_k_models, k_model_getopt_t *k_models,
580 list_t **pDomain_list, int *pNum_folded)
582 char *prog_name = NULL;
583 char c, options[] = "P:t:v:T:C:m:a:M:A:k:K:F:U:Vh";
584 int ftension_model=0, utension_model=0, k_model=0;
587 extern int optind, optopt, opterr;
592 flags = FLAG_OUTPUT_UNFOLDING_FORCES;
594 *pP = 1e-3; /* % pop per s */
595 *pDt_max = 0.001; /* s */
596 *pV = 1e-6; /* m/s */
597 env->T = 300.0; /* K */
599 while ((c=getopt(argc, argv, options)) != -1) {
601 case 'P': *pP = atof(optarg); break;
602 case 't': *pDt_max = atof(optarg); break;
603 case 'v': *pV = atof(optarg); break;
604 case 'T': env->T = atof(optarg); break;
605 case 'C': env->T = atof(optarg)+273.15; break;
607 ftension_model = index_tension_model(n_tension_models, tension_models, optarg);
610 tension_models[ftension_model].params = optarg;
613 utension_model = index_tension_model(n_tension_models, tension_models, optarg);
616 tension_models[utension_model].params = optarg;
619 k_model = index_k_model(n_k_models, k_models, optarg);
622 k_models[k_model].params = optarg;
627 for (i=0; i<n; i++) {
628 push(pDomain_list, generate_domain(DS_FOLDED, tension_models+ftension_model, tension_models+utension_model, k_models+k_model));
635 for (i=0; i<n; i++) {
636 push(pDomain_list, generate_domain(DS_UNFOLDED, tension_models+ftension_model, tension_models+utension_model, k_models+k_model));
640 flags = FLAG_OUTPUT_FULL_CURVE;
643 fprintf(stderr, "unrecognized option '%c'\n", optopt);
644 /* fall through to default case */
646 help(prog_name, *pP, *pDt_max, *pV, env, n_tension_models, tension_models, ftension_model, utension_model, n_k_models, k_models, k_model);
650 /* check the arguments */
651 assert(*pDomain_list != NULL); /* no domains! */
652 assert(*pP > 0.0); assert(*pP < 1.0);
653 assert(*pDt_max > 0.0);
655 assert(env->T > 0.0);
660 <<index tension model>>=
661 int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
664 for (i=0; i<n_models; i++)
665 if (STRMATCH(models[i].name, name))
668 fprintf(stderr, "Unrecognized tension model '%s'\n", name);
676 int index_k_model(int n_models, k_model_getopt_t *models, char *name)
679 for (i=0; i<n_models; i++)
680 if (STRMATCH(models[i].name, name))
683 fprintf(stderr, "Unrecognized unfolding rate model '%s'\n", name);
690 <<initialize model definition>>=
691 /* requires int num_param_args and char **param_args in the current scope
693 * INIT_MODEL("folded", folded_model, folded_params);
694 * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types.
696 #define INIT_MODEL(role, model, param_pointer) \
698 parse_list_string(model->params, SEP, '{', '}', \
699 &num_param_args, ¶m_args); \
700 if (num_param_args != model->num_params) { \
702 "%s model %s expected %d params," \
703 role, model->name, model->num_params); \
705 "not the %d params in '%s'\n", \
706 num_param_args, model->params); \
707 assert(num_param_args == model->num_params); \
709 if (model->creator) \
710 param_pointer = (*model->creator)(param_args); \
712 param_pointer = NULL; \
713 free_parsed_list(num_param_args, param_args); \
718 void *generate_domain(enum domain_state_t state,
719 tension_model_getopt_t *folded_model,
720 tension_model_getopt_t *unfolded_model,
721 k_model_getopt_t *k_model)
723 void *ftension_params, *utension_params, *k_params;
724 int num_param_args; /* for INIT_MODEL() */
725 char **param_args; /* for INIT_MODEL() */
728 fprintf(stderr, "generating %s ", state==DS_FOLDED ? "folded" : "unfolded");
729 fprintf(stderr, "domain (k: \"%s\", \"%s\", f: \"%s\", \"%s\", u: \"%s\", \"%s\")\n",
730 k_model->name, k_model->params,
731 folded_model->name, folded_model->params,
732 unfolded_model->name, unfolded_model->params);
734 INIT_MODEL("folded", folded_model, ftension_params);
735 INIT_MODEL("unfolded", unfolded_model, utension_params);
736 INIT_MODEL("k", k_model, k_params);
737 return create_domain(state,
738 k_model->k, k_params,
740 folded_model->tg_group, ftension_params,
741 folded_model->destructor,
742 unfolded_model->tg_group, utension_params,
743 unfolded_model->destructor);
749 \addcontentsline{toc}{section}{Appendicies}
751 \section{sawsim.c details}
755 The general layout of our simulation code is:
765 We include [[math.h]], so don't forget to link to the libm with `-lm'.
767 #include <assert.h> /* assert() */
768 #include <stdlib.h> /* malloc(), free(), rand() */
769 #include <stdio.h> /* fprintf(), stdout */
770 #include <string.h> /* strlen, strtok() */
771 #include <math.h> /* exp(), M_PI, sqrt() */
772 #include <sys/types.h> /* pid_t (returned by getpid()) */
773 #include <unistd.h> /* getpid() (for seeding rand()), getopt() */
776 #include "tension_balance.h"
778 #include "tension_model.h"
784 <<version definition>>
786 <<max/min definitions>>
787 <<string comparison definition and globals>>
788 <<initialize model definition>>
789 <<get options definitions>>
790 <<domain definitions>>
799 <<create/destroy domain>>
800 <<destroy domain list>>
802 <<simulation functions>>
803 <<boilerplate functions>>
806 <<boilerplate functions>>=
808 <<get option functions>>
812 void setup(one_dim_func_t **tension_handler)
814 srand(getpid()*time(NULL)); /* seed rand() */
815 tension_handler[TG_NULL] = NULL;
816 tension_handler[TG_CONST] = &const_tension_handler;
817 tension_handler[TG_HOOKE] = &hooke_handler;
818 tension_handler[TG_WLC] = &wlc_handler;
819 tension_handler[TG_FJC] = &fjc_handler;
823 <<flag definitions>>=
824 /* in octal b/c of prefixed '0' */
825 #define FLAG_OUTPUT_FULL_CURVE 01
826 #define FLAG_OUTPUT_UNFOLDING_FORCES 02
830 static unsigned long int flags = 0;
833 \subsection{Utilities}
836 <<max/min definitions>>=
837 #define MAX(a,b) ((a)>(b) ? (a) : (b))
838 #define MIN(a,b) ((a)<(b) ? (a) : (b))
841 Note that [[STRMATCH]] chokes if one of the strings is [[NULL]].
842 <<string comparison definition and globals>>=
843 // Check if two strings match, return 1 if they do
844 static char *temp_string_A;
845 static char *temp_string_B;
846 #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \
847 strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \
848 !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) )
849 /* +1 to also compare the '\0' */
852 We also define a macro for our [[check]] unit testing
853 <<check relative error>>=
854 #define CHECK_ERR(max_err, expected, received) \
856 fail_unless( (received-expected)/expected < max_err, \
857 "relative error %g >= %g in %s (Expected %g, received %g)", \
858 (received-expected)/expected, max_err, #received, \
859 expected, received); \
860 fail_unless(-(received-expected)/expected < max_err, \
861 "relative error %g >= %g in %s (Expected %g, received %g)", \
862 -(received-expected)/expected, max_err, #received, \
863 expected, received); \
868 int happens(double probability)
870 assert(probability >= 0.0); assert(probability <= 1.0);
871 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*/
876 \subsection{Adaptive timesteps}
877 \label{app.adaptive_dt}
879 $F(x)$ increases with $x$, possibly exploding, as in the worm-like chain model,
880 so basing the timestep on the the unfolding probability at the current tension
881 is dangerous, and we need to search for a $dt$ for which
882 $P(F(x+v*dt)) < P_\text{target}$.
883 There are two cases to consider.
884 In the most common, no domains have unfolded since the last step,
885 and we expect the next step to be slightly shorter than the previous one.
886 In the less common, domains did unfold in the last step,
887 and we expect the next step to be considerably longer than the previous one.
889 double search_dt(one_dim_func_t **tension_handler,
891 environment_t *env, double x,
892 double target_prob, double max_dt, double v)
893 { /* Returns the timestep dt in seconds for the current folded domain.
894 * Takes a list of tension handlers, the list of domains,
895 * a pointer env to the environmental data, a starting separation x in m,
896 * a target_prob between 0 and 1,
897 * max_dt in s, stretching velocity v in m/s.
899 double F, k, dtCur, dtU, dtUCur, dtL, dt;
901 /* get upper bound using the current position */
902 F = find_tension(tension_handler, domain_list, env, x); /* BUG. repeated calculation */
903 //printf("Start with x = %g (F = %g)\n", x, F);
904 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
905 //printf("x %g\tF %g\tk %g\n", x, F, k);
906 dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */
908 //printf("overshot max_dt\n");
911 /* set a lower bound on dt too */
914 /* The dt determined above may produce illegitimate forces or ks.
915 * Reduce the upper bound until we have valid ks. */
917 F = find_tension(tension_handler, domain_list, env, x+v*dt);
918 while (F == HUGE_VAL) { /* reduce step until we hit a valid force */
921 F = find_tension(tension_handler, domain_list, env, x+v*dt);
923 //printf("Try for dt = %g (F = %g)\n", dt, F);
924 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
925 /* returned k may be -1.0 */
926 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
927 while (k == -1.0) { /* reduce step until we hit a valid k */
929 dt = dtU; /* hopefully, we can use the max dt, see if it works */
930 F = find_tension(tension_handler, domain_list, env, x+v*dt);
931 //printf("Try for dt = %g (F = %g)\n", dt, F);
932 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
933 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
935 assert(dtU > 1e-14); /* timestep to valid k too small */
936 dtUCur = target_prob / k; /* safe timestep back from x+dtU */
938 return dt; /* dtU is safe. */
940 /* dtU wasn't safe, lets see what would be. */
941 while (dtU > 1.1*dtL) { /* until the range is reasonably small */
942 dt = (dtU + dtL) / 2.0;
943 F = find_tension(tension_handler, domain_list, env, x+v*dt);
944 //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
945 k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
946 dtCur = target_prob / k;
947 //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur);
948 if (dtCur > dt) /* safe timestep back from x+dt covers dt */
950 else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */
951 dtU = dt; /* ... stepping out only dtCur would be safe */
954 break; /* dtCur = dt */
956 return MAX(dtUCur, dtL);
960 To determine $dt$ for an array of potentially different folded domains,
961 we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
964 double determine_dt(one_dim_func_t **tension_handler, list_t *domain_list,
965 environment_t *env, double x,
966 double target_prob, double dt_max, double v)
967 { /* Returns the timestep dt in seconds.
968 * Takes the list of folded domains, target_prob between 0 and 1,
969 * F in N, and T in K. */
970 double dt=dt_max, new_dt;
971 assert(target_prob > 0.0); assert(target_prob < 1.0);
972 assert(dt_max > 0.0);
974 /* .5 nm steps = v * dt */
976 while (domain_list != NULL) {
977 if (D(domain_list)->state == DS_FOLDED) {
978 new_dt = search_dt(tension_handler, domain_list, env, x, target_prob, dt, v);
979 dt = MIN(dt, new_dt);
981 domain_list = domain_list->next;
987 \subsection{Domain data}
989 Currently domains exist in two states, folded and unfolded,
990 and the only allowed transitions are folded $\rightarrow$ unfolded.
991 Of course, it wouldn't be too complicated to extent this to a multi-state system,
992 with an array containing the domains group for each possible state,
993 and a matrix of transition-rate-calculating functions.
994 However, at this point such generality seems unnecessary at this point.
995 <<domain definitions>>=
996 enum domain_state_t {DS_FOLDED,
1000 typedef struct domain_struct {
1001 enum domain_state_t state;
1002 enum tension_group_t folded_group;
1003 enum tension_group_t unfolded_group;
1004 k_func_t *k; /* function returning unfolding rate */
1005 void *folded_params; /* pointer to folded parameters */
1006 void *unfolded_params; /* pointer to unfolded parameters */
1007 void *k_params; /* pointer to k parameters */
1008 destroy_data_func_t *destroy_folded;
1009 destroy_data_func_t *destroy_unfolded;
1010 destroy_data_func_t *destroy_k;
1013 /* get the domain data for the current list node */
1014 #define D(list) ((domain_t *)(list)->d)
1015 /* get the tension params for the current list node */
1016 #define D_TP(list) (((domain_t *)(list)->d)->state == DS_FOLDED \
1017 ? ((domain_t *)(list)->d)->folded_params \
1018 : ((domain_t *)(list)->d)->unfolded_params)
1020 [[k]] is a pointer to the function determining the unfolding rate for a given tension.
1021 [[folded_params]] is a pointer to the parameters used by the function pointed to by [[k]].
1022 [[unfolded_params]] is a pointer to the parameters used by the group-appropriate handler function when determining the tension.
1023 The [[destroy_*]] pointers point to functions for freeing the memory [[*_params]].
1024 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.
1026 [[create_]] and [[destroy_domain]] are simple wrappers around [[malloc]] and [[free]].
1027 <<create/destroy domain>>=
1028 domain_t *create_domain(enum domain_state_t state,
1031 destroy_data_func_t *destroy_k,
1032 enum tension_group_t folded_group,
1033 void *folded_params,
1034 destroy_data_func_t *destroy_folded,
1035 enum tension_group_t unfolded_group,
1036 void *unfolded_params,
1037 destroy_data_func_t *destroy_unfolded)
1039 domain_t *ret = (domain_t *)malloc(sizeof(domain_t));
1040 assert(ret != NULL);
1041 if (state == DS_FOLDED) {
1042 assert(k != NULL); /* the pointer points somewhere valid */
1043 assert(*k != NULL); /* and there is something useful there */
1046 ret->folded_group = folded_group;
1047 ret->unfolded_group = unfolded_group;
1049 ret->k_params = k_params;
1050 ret->destroy_k = destroy_k;
1051 ret->folded_params = folded_params;
1052 ret->unfolded_params = unfolded_params;
1053 ret->destroy_folded = destroy_folded;
1054 ret->destroy_unfolded = destroy_unfolded;
1058 void destroy_domain(domain_t *domain)
1061 //printf("domain %p & %p\n", *domain, domain);
1062 if (domain->destroy_folded)
1063 (*domain->destroy_folded)(domain->folded_params);
1064 if (domain->destroy_unfolded)
1065 (*domain->destroy_unfolded)(domain->unfolded_params);
1066 if (domain->destroy_k)
1067 (*domain->destroy_k)(domain->k_params);
1073 <<destroy domain list>>=
1074 void destroy_domain_list(list_t *domain_list)
1076 domain_list = head(domain_list);
1077 while (domain_list != NULL)
1078 destroy_domain((domain_t *) pop(&domain_list));
1082 \subsection{Group handling}
1084 <<group functions>>=
1090 enum tension_group_t get_group(domain_t *domain)
1092 if (domain->state == DS_FOLDED)
1093 return domain->folded_group;
1095 assert(domain->state == DS_UNFOLDED);
1096 return domain->unfolded_group;
1102 list_t *get_group_list(list_t *list, enum tension_group_t group)
1106 while (list != NULL) {
1107 if (get_group(D(list)) == group)
1108 push(&ret, D_TP(list)); /* add a pointer to the appropriate tension parameters to our new list. */
1114 Because all the node data in lists returned by [[get_group_list]] is also in the main domain list,
1115 you shouldn't destroy and node data popped off when destroying the group lists.
1116 It will all get cleaned up when the main domain list is destroyed.
1119 \section{String parsing}
1121 For handling command line arguments, we need parse delimited strings (e.g. [["param1,param2,param3..."]]).
1122 The model handling in getopt is set up to handle a fixed number of arguments for each model,
1123 so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots)
1124 need to take care of parsing those parameters themselves.
1125 We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]].
1129 <<parse definitions>>
1130 <<parse declarations>>
1133 <<parse module makefile lines>>=
1135 notangle -Rparse.c $^ > $@
1137 notangle -Rparse.h $^ > $@
1138 check_parse.c : sawsim.nw
1139 notangle -Rcheck-parse.c $^ > $@
1140 check_parse : check_parse.c parse.c parse.h
1141 gcc -g -o $@ $< parse.c -lcheck
1143 rm -f parse.c parse.h check_parse.c check_parse
1146 <<parse definitions>>=
1147 #define SEP ',' /* argument separator character */
1150 <<parse declarations>>=
1151 extern void parse_list_string(char *string, char sep, char deeper, char shallower,
1152 int *num, char ***string_array);
1153 extern void free_parsed_list(int num, char **string_array);
1156 [[parse_list_string]] allocates memory, don't forget to free it afterward with [[free_parsed_list]].
1157 It does not alter the original.
1159 The string may start off with a [[deeper]] character (i.e. [["{x,y}"]]),
1160 and when it does, brace stripping will set leave [[x,y]], where the pointer is one character in on the copied string.
1161 However, when we go to free the memory, we need a pointer to the beginning of the string.
1162 In order to accommodate this for a string with $N$ argument, allocate a pointer array with $N+1$ elements,
1163 let the first $N$ elements point to the separated arguments,
1164 and let the last element point to the start of the copied string regardless of braces.
1165 <<parse delimited list string functions>>=
1166 /* TODO, split out into parse.hc */
1167 static int next_delim_index(char *string, char sep, char deeper, char shallower)
1170 while (string[i] != '\0' && !(string[i] == sep && depth == 0)) {
1171 if (string[i] == deeper) {depth++;}
1172 else if (string[i] == shallower) {depth--; assert(depth >= 0);}
1178 void parse_list_string(char *string, char sep, char deeper, char shallower,
1179 int *num, char ***string_array)
1181 char *str=NULL, **ret=NULL;
1183 if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */
1185 *string_array = NULL;
1188 /* make a copy of the string, so we don't change the original */
1189 str = (char *)malloc(sizeof(char)*(strlen(string)+1));
1190 assert(str != NULL);
1191 strcpy(str, string); /* we know str is long enough */
1192 /* count the number of regions, so we can allocate pointers to them */
1195 n++; i++; /* move on to next argument */
1196 i += next_delim_index(str+i, sep, deeper, shallower);
1197 //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr);
1198 } while (str[i] != '\0');
1199 ret = (char **)malloc(sizeof(char *)*(n+1));
1200 assert(ret != NULL);
1201 /* replace the separators with '\0' & assign pointers */
1202 ret[n] = str; /* point to the front of the copied string */
1205 for(i=1; i<n; i++) {
1206 j += next_delim_index(str+j, sep, deeper, shallower);
1208 ret[i] = str+j; /* point to the separated arguments */
1210 /* strip off bounding braces, if any */
1211 for(i=0; i<n; i++) {
1212 if (ret[i][0] == deeper) {
1216 if (ret[i][strlen(ret[i])-1] == shallower)
1217 ret[i][strlen(ret[i])-1] = '\0';
1220 *string_array = ret;
1223 void free_parsed_list(int num, char **string_array)
1226 /* frees all items, since they were allocated as a single string*/
1227 free(string_array[num]);
1228 /* and free the array of pointers */
1234 \subsection{Parsing implementation}
1240 <<parse delimited list string functions>>
1244 #include <assert.h> /* assert() */
1245 #include <stdlib.h> /* NULL */
1246 #include <stdio.h> /* fprintf(), stdout *//*!!*/
1247 #include <string.h> /* strlen() */
1251 \subsection{Parsing unit tests}
1253 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1255 <<parse check includes>>
1256 <<string comparison definition and globals>>
1257 <<check relative error>>
1258 <<parse test suite>>
1259 <<main check program>>
1262 <<parse check includes>>=
1263 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
1264 #include <stdio.h> /* printf() */
1265 #include <assert.h> /* assert() */
1266 #include <string.h> /* strlen() */
1271 <<parse test suite>>=
1272 <<parse list string tests>>
1273 <<parse suite definition>>
1276 <<parse suite definition>>=
1277 Suite *test_suite (void)
1279 Suite *s = suite_create ("k model");
1280 <<parse list string test case defs>>
1282 <<parse list string test case add>>
1287 <<parse list string tests>>=
1290 START_TEST(test_next_delim_index)
1292 fail_unless(next_delim_index("", ',', '{', '}')==0, NULL);
1293 fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL);
1294 fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL);
1295 fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL);
1296 fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL);
1300 START_TEST(test_parse_list_null)
1304 parse_list_string(NULL, SEP, '{', '}',
1305 &num_param_args, ¶m_args);
1306 fail_unless(num_param_args == 0, NULL);
1307 fail_unless(param_args == NULL, NULL);
1310 START_TEST(test_parse_list_single_simple)
1314 parse_list_string("arg", SEP, '{', '}',
1315 &num_param_args, ¶m_args);
1316 fail_unless(num_param_args == 1, NULL);
1317 fail_unless(STRMATCH(param_args[0],"arg"), NULL);
1320 START_TEST(test_parse_list_single_compound)
1324 parse_list_string("{x,y,z}", SEP, '{', '}',
1325 &num_param_args, ¶m_args);
1326 fail_unless(num_param_args == 1, NULL);
1327 fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z");
1330 START_TEST(test_parse_list_double_simple)
1334 parse_list_string("abc,def", SEP, '{', '}',
1335 &num_param_args, ¶m_args);
1336 fail_unless(num_param_args == 2, NULL);
1337 fail_unless(STRMATCH(param_args[0],"abc"), NULL);
1338 fail_unless(STRMATCH(param_args[1],"def"), NULL);
1342 <<parse list string test case defs>>=
1343 TCase *tc_parse_list_string = tcase_create("parse list string");
1345 <<parse list string test case add>>=
1346 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
1347 tcase_add_test(tc_parse_list_string, test_parse_list_null);
1348 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
1349 tcase_add_test(tc_parse_list_string, test_parse_list_single_compound);
1350 tcase_add_test(tc_parse_list_string, test_parse_list_double_simple);
1351 suite_add_tcase(s, tc_parse_list_string);
1355 \section{Unit tests}
1357 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1366 <<main check program>>
1378 <<determine dt tests>>
1380 <<does domain unfold tests>>
1381 <<randomly unfold domains tests>>
1382 <<suite definition>>
1385 <<suite definition>>=
1386 Suite *test_suite (void)
1388 Suite *s = suite_create ("sawsim");
1389 <<F test case defs>>
1390 <<determine dt test case defs>>
1391 <<happens test case defs>>
1392 <<does domain unfold test case defs>>
1393 <<randomly unfold domains test case defs>>
1396 <<determine dt test case add>>
1397 <<happens test case add>>
1398 <<does domain unfold test case add>>
1399 <<randomly unfold domains test case add>>
1402 tcase_add_checked_fixture(tc_strip_address,
1403 setup_strip_address,
1404 teardown_strip_address);
1410 <<main check program>>=
1414 Suite *s = test_suite();
1415 SRunner *sr = srunner_create(s);
1416 srunner_run_all(sr, CK_ENV);
1417 number_failed = srunner_ntests_failed(sr);
1419 return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
1423 \subsection{F tests}
1425 <<worm-like chain tests>>
1427 <<F test case defs>>=
1428 <<worm-like chain test case def>>
1430 <<F test case add>>=
1431 <<worm-like chain test case add>>
1434 <<worm-like chain tests>>=
1435 START_TEST(test_wlc_at_zero)
1437 double T=1.0, L=1.0, p=0.1, x=0.0;
1438 fail_unless(wlc(x, T, p, L)==0, NULL);
1442 START_TEST(test_wlc_at_half)
1444 double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
1445 /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
1446 * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
1448 * linear term = x/L = 0.5
1449 * nonlinear + linear = 0.75 + 0.5 = 1.25
1450 * wlc = 10e21*1.25 = 12.5e21
1452 fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
1453 "wlc(%g, %g, %g, %g) = %g != %g",
1454 x, T, p, L, wlc(x, T, p, L), 12.5e21);
1458 <<worm-like chain test case def>>=
1459 TCase *tc_wlc = tcase_create("WLC");
1462 <<worm-like chain test case add>>=
1463 tcase_add_test(tc_wlc, test_wlc_at_zero);
1464 tcase_add_test(tc_wlc, test_wlc_at_half);
1465 suite_add_tcase(s, tc_wlc);
1468 \subsection{Model tests}
1470 Check the searching with [[linear_k]].
1471 Check overwhelming force treatment with the heavyside-esque step [[?]].
1472 <<determine dt tests>>=
1473 double linear_k(double F, environment_t *env, void *params)
1475 double Fnot = *(double *)params;
1479 START_TEST(test_determine_dt_linear_k)
1482 double dt_max=3.0, Fnot=3.0;
1483 double F[]={0,1,2,3,4,5,6};
1484 domain_t dom; /* use both parts at once for folded/unfolded */
1488 dom->next = dom->prev = NULL;
1489 dom->k_func_t = linear_k;
1490 dom->folded_params = &Fnot;
1491 dom->unfolded_params = !!!!!!!!!
1492 dom->destroy_folded = dom->destroy_unfolded = NULL;
1493 for( i=0; i < sizeof(F)/sizeof(double); i++) {
1494 fail_unless(determine_dt(folded, unfolded, targ_p, dt_max, v, x, T)==1, NULL);
1500 <<determine dt test case defs>>=
1501 TCase *tc_determine_dt = tcase_create("Determine dt");
1503 <<determine dt test case add>>=
1504 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
1505 suite_add_tcase(s, tc_determine_dt);
1510 <<happens test case defs>>=
1512 <<happens test case add>>=
1515 <<does domain unfold tests>>=
1517 <<does domain unfold test case defs>>=
1519 <<does domain unfold test case add>>=
1522 <<randomly unfold domains tests>>=
1524 <<randomly unfold domains test case defs>>=
1526 <<randomly unfold domains test case add>>=
1530 \section{Balancing group extensions}
1531 \label{app.tension_balance}
1533 For a given total extension $x$ (of the piezo), the various domain groups (WLC, FJC, Hookean springs, \ldots) extend different amounts, and we need to tweak the portion that each extends to balance the tension amongst the active groups.
1534 Since the tension balancing is separable from the bulk of the simulation, we break it out into a separate module.
1535 The interface is defined in [[tension_balance.h]], the implementation in [[tension_balance.c]], and the unit testing in [[check_tension_balance.c]]
1537 <<tension-balance.h>>=
1539 <<tension balance function declaration>>
1542 <<tension balance functions>>=
1543 <<one dimensional bracket>>
1544 <<one dimensional solve>>
1546 <<tension balance function>>
1549 <<tension balance module makefile lines>>=
1550 tension_balance.c : sawsim.nw
1551 notangle -Rtension-balance.c $^ > $@
1552 tension_balance.h : sawsim.nw
1553 notangle -Rtension-balance.h $^ > $@
1554 check_tension_balance.c : sawsim.nw
1555 notangle -Rcheck-tension-balance.c $^ > $@
1556 check_tension_balance : check_tension_balance.c global.h tension_balance.c tension_balance.h
1557 gcc -g -o $@ $< tension_balance.c -lcheck
1559 rm -f tension_balance.c tension_balance.h check_tension_balance.c check_tension_balance
1562 The entire force balancing problem reduces to a solving two nested one-dimensional root-finding problems.
1563 First we define the one dimensional tension $F(x_0)$ (where $i=0$ is the index of the first non-empty group).
1564 $x(x_0)$ is also strictly monotonically increasing, so we can solve for $x_0$ such that $\sum_i x_i = x$.
1565 <<tension balance function declaration>>=
1566 double tension_balance(int num_tension_groups,
1567 one_dim_func_t **tension_handler,
1568 void **params, int *active,
1569 int new_active_groups, double *xi,
1570 double last_x, double x);
1572 <<tension balance function>>=
1573 double tension_balance(int num_tension_groups,
1574 one_dim_func_t **tension_handler,
1575 void **params, int *active,
1576 int new_active_groups, double *xi,
1577 double last_x, double x)
1578 { /* xi initialized to x values for groups at last x,
1579 * returned as x values for groups at current x */
1581 one_dim_func_t **active_handlers=NULL;
1582 void **active_params=NULL;
1583 double *active_xi=NULL;
1584 int i, active_groups=0;
1585 x_of_xo_data_t x_xo_data;
1586 double min_dx=1e-10, min_dy=1e-15;
1590 assert(num_tension_groups > 0);
1591 active_handlers = (one_dim_func_t **)malloc(sizeof(one_dim_func_t *)*num_tension_groups);
1592 assert(active_handlers != NULL);
1593 active_params = (void **)malloc(sizeof(void *)*num_tension_groups);
1594 assert(active_params != NULL);
1595 active_xi = (double *)malloc(sizeof(double)*num_tension_groups);
1596 assert(active_xi != NULL);
1598 for (i=0; i<num_tension_groups; i++) { /* remove null groups */
1600 active_handlers[active_groups] = tension_handler[i];
1601 active_params[active_groups] = params[i];
1602 active_xi[active_groups] = xi[i]; /* last balanced x_i if no new active groups */
1606 x_xo_data.n_groups = active_groups;
1607 x_xo_data.tension_handler = active_handlers;
1608 x_xo_data.handler_data = active_params;
1609 x_xo_data.xi = active_xi;
1610 if (active_groups == 1) { /* only one group, no balancing required */
1614 //printf("balancing for x = %g with ", x);
1615 //for(i=0; i<num_tension_groups; i++) printf(" %d: %g\t", i, xi[i]);
1617 if (new_active_groups || 1) { /* take a rough guess and attempt boundary conditions. */
1618 //printf("search bracket\n");
1619 oneD_bracket(x_of_xo, &x_xo_data, x, x/active_groups, &lb, &ub);
1620 } else { /* work off the last balanced point */
1621 //printf("step off the old bracketing x %g + %g = %g\n", last_x, x-last_x, x);
1624 ub = active_xi[0]+ x-last_x; /* apply all change to x0 */
1625 } else if (x < last_x) {
1626 lb = active_xi[0]- (x-last_x); /* apply all change to x0 */
1628 } else { /* x == last_x */
1629 printf("not moving\n");
1634 //printf("lb %g,\tub %g\n", lb, ub);
1635 xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_dx, min_dy, max_steps, NULL);
1637 F = (*active_handlers[0])(xo, active_params[0]);
1638 /* go back through and place the active xi data in the complete xi array */
1640 for (i=0; i<num_tension_groups; i++) {
1642 xi[i] = active_xi[active_groups];
1647 if (active_groups > 1 && 0) {
1648 printf("balanced for x = %g with ", x);
1649 for(i=0; i<num_tension_groups; i++) printf(" %d: %g\t", i, xi[i]);
1652 free(active_handlers);
1653 free(active_params);
1660 <<tension balance internal definitions>>=
1661 <<x of x0 definitions>>
1664 <<x of x0 definitions>>=
1665 typedef struct x_of_xo_data_struct {
1667 one_dim_func_t **tension_handler; /* array of fn pointers */
1668 void **handler_data; /* array of void* pointers */
1674 double x_of_xo(double xo, void *pdata)
1676 x_of_xo_data_t *data = (x_of_xo_data_t *)pdata;
1677 double F, x=0, xi, lb, ub;
1679 double min_dx=1e-10, min_dy=1e-14;
1681 assert(data->n_groups > 0);
1683 F = (*data->tension_handler[0])(xo, data->handler_data[0]);
1685 if (data->xi) data->xi[0] = xo;
1686 for (i=1; i < data->n_groups; i++) {
1687 oneD_bracket(data->tension_handler[i], data->handler_data[i], F, data->xi[i], &lb, &ub);
1688 xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_dx, min_dy, max_steps, NULL);
1691 if (data->xi) data->xi[i] = xi;
1697 Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited number of steps for monotonically increasing $f(x)$.
1698 Simple bisection, so it's robust and converges fairly quickly.
1699 <<one dimensional function definition>>=
1700 /* equivalent to gsl_func_t */
1701 typedef double one_dim_func_t(double x, void *params);
1703 <<one dimensional solve declaration>>=
1704 double oneD_solve(one_dim_func_t *f, void *params, double y, double lb, double ub,
1705 double min_dx, double min_dy, int max_steps,
1709 <<one dimensional solve>>=
1710 /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */
1711 double oneD_solve(one_dim_func_t *f, void *params, double y, double lb, double ub,
1712 double min_dx, double min_dy, int max_steps,
1715 double yx, ylb, yub, x;
1718 ylb = (*f)(lb, params);
1719 yub = (*f)(ub, params);
1721 /* check some simple cases */
1723 if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bounded */
1724 else return lb; /* any x on the range [lb, ub] would work */
1726 if (ylb == y) { x=lb; yx=ylb; goto end; }
1727 if (yub == y) { x=ub; yx=yub; goto end; }
1729 //printf("lb %g, x %g, ub %g\tylb %g, y %g, yub %g\n", lb, x, ub, ylb, y, yub);
1730 assert(ylb < y); assert(yub > y);
1731 for (i=0; i<max_steps; i++) {
1733 yx = (*f)(x, params);
1734 if (ub-lb < min_dx || yub-ylb < min_dy || yx == y) break;
1735 else if (yx > y) { ub=x; yub=yx; }
1736 else /* yx < y */ { lb=x; ylb=yx; }
1744 The one dimensional solver needs a bracketed solution, which sometimes we don't have.
1745 Generate bracketing $x$ values through bisection or exponential growth.
1746 <<one dimensional bracket declaration>>=
1747 void oneD_bracket(one_dim_func_t *f, void *params, double y, double xguess, double *lb, double *ub);
1750 <<one dimensional bracket>>=
1751 void oneD_bracket(one_dim_func_t *f, void *params, double y, double xguess, double *lb, double *ub)
1753 double yx, step, x=xguess;
1755 yx = (*f)(x, params);
1756 //fprintf(stdout, "bracketing %g, start at f(%g) = %g\n", y, x, yx);
1763 if (x == 0) x = 0.5; /* guess a scale of 1.0 */
1766 yx = (*f)(x, params);
1767 //fprintf(stdout, "increasing to f(%g) = %g\n", x, yx);
1771 //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
1775 \subsection{Balancing implementation}
1777 <<tension-balance.c>>=
1779 <<tension balance includes>>
1780 #include "tension_balance.h"
1781 <<tension balance internal definitions>>
1782 <<tension balance functions>>
1785 <<tension balance includes>>=
1786 #include <assert.h> /* assert() */
1787 #include <stdlib.h> /* NULL */
1788 #include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
1789 #include <stdio.h> /* fprintf(), stdout */
1793 \subsection{Balancing unit tests}
1795 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
1796 <<check-tension-balance.c>>=
1797 <<tension balance check includes>>
1798 <<tension balance test suite>>
1799 <<main check program>>
1802 <<tension balance check includes>>=
1803 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
1804 #include <assert.h> /* assert() */
1807 #include "tension_balance.h"
1810 <<tension balance test suite>>=
1811 <<tension balance function tests>>
1812 <<tension balance suite definition>>
1815 <<tension balance suite definition>>=
1816 Suite *test_suite (void)
1818 Suite *s = suite_create ("tension balance");
1819 <<tension balance function test case defs>>
1821 <<tension balance function test case adds>>
1826 <<tension balance function tests>>=
1827 <<check relative error>>
1829 double hooke(void *pK, double x)
1832 return *((double*)pK) * x;
1835 START_TEST(test_single_function)
1837 double k=5, x=3, last_x=2.0, F;
1838 one_dim_func_t *handlers[] = {&hooke};
1839 void *data[] = {&k};
1840 double xi[] = {0.0};
1842 int new_active_groups = 1;
1843 F = tension_balance(1, handlers, data, active, new_active_groups, xi, last_x, x);
1844 fail_unless(F = k*x, NULL);
1849 We can also test balancing two springs with different spring constants.
1850 The analytic solution for a general number of springs is given in Appendix \ref{app.math_hooke}.
1851 <<tension balance function tests>>=
1852 START_TEST(test_double_hooke)
1854 double k1=5, k2=4, x=3, last_x=2.0, F, Fe, x1e, x2e;
1855 one_dim_func_t *handlers[] = {&hooke, &hooke, NULL};
1856 void *data[] = {&k1, &k2, NULL};
1857 double xi[] = {0.0, 0.0, 0.0};
1858 int active[] = {1, 1, 0};
1859 int new_active_groups = 1;
1860 F = tension_balance(3, handlers, data, active, new_active_groups, xi, last_x, x);
1864 //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F);
1865 CHECK_ERR(1e-6, x1e, xi[0]);
1866 CHECK_ERR(1e-6, x2e, xi[1]);
1867 CHECK_ERR(1e-6, Fe, F);
1872 <<tension balance function test case defs>>=
1873 TCase *tc_tbfunc = tcase_create("tension balance function");
1876 <<tension balance function test case adds>>=
1877 tcase_add_test(tc_tbfunc, test_single_function);
1878 tcase_add_test(tc_tbfunc, test_double_hooke);
1879 suite_add_tcase(s, tc_tbfunc);
1884 The globular domains (and cantilever params, etc.) are saved in bi-directional lists.
1885 Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module.
1886 The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]]
1890 <<list definitions>>
1891 <<list declarations>>
1894 <<list declarations>>=
1895 <<head and tail declarations>>
1896 <<list length declaration>>
1897 <<push declaration>>
1902 <<create and destroy node>>
1909 <<list module makefile lines>>=
1911 notangle -Rlist.c $^ > $@
1913 notangle -Rlist.h $^ > $@
1914 check_list.c : sawsim.nw
1915 notangle -Rcheck-list.c $^ > $@
1916 check_list : check_list.c global.h list.c list.h
1917 gcc -g -o $@ $< list.c -lcheck
1919 rm -f list.c list.h check_list.c check_list
1922 <<list definitions>>=
1923 typedef struct list_struct {
1924 struct list_struct *next;
1925 struct list_struct *prev;
1931 [[head]] and [[tail]] return pointers to the head and tail nodes of the list:
1932 <<head and tail declarations>>=
1933 list_t *head(list_t *list);
1934 list_t *tail(list_t *list);
1937 list_t *head(list_t *list)
1941 while (list->prev != NULL) {
1947 list_t *tail(list_t *list)
1951 while (list->next != NULL) {
1958 <<list length declaration>>=
1959 int list_length(list_t *list);
1962 int list_length(list_t *list)
1969 while (list->next != NULL) {
1977 [[push]] inserts a node after the current position in the list
1978 without changing the current position.
1979 However, if the list we're pushing onto is [[NULL]], the current position isn't defined
1980 so we set it to that of the pushed domain.
1981 <<push declaration>>=
1982 void push(list_t **pList, void *data);
1985 void push(list_t **pList, void *data)
1987 list_t *list, *new_node;
1988 assert(pList != NULL);
1990 new_node = create_node(data);
1994 if (list->next != NULL)
1995 list->next->prev = new_node;
1996 new_node->next = list->next;
1997 list->next = new_node;
1998 new_node->prev = list;
2003 [[pop]] removes the current domain node,
2004 moving the current position to the node after it,
2005 unless that node is [[NULL]],
2006 in which case move the current position to the node before it.
2007 <<pop declaration>>=
2008 void *pop(list_t **pList);
2011 void *pop(list_t **pList)
2013 list_t *list, *popped;
2015 assert(pList != NULL);
2017 assert(list != NULL); /* not an empty list */
2019 /* bypass the popped node */
2020 if (list->prev != NULL)
2021 list->prev->next = list->next;
2022 if (list->next != NULL) {
2023 list->next->prev = list->prev;
2024 *pList = list->next;
2026 *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */
2028 destroy_node(popped);
2033 [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
2034 <<create and destroy node>>=
2035 list_t *create_node(void *data)
2037 list_t *ret = (list_t *)malloc(sizeof(list_t));
2038 assert(ret != NULL);
2045 void destroy_node(list_t *node)
2051 The user must free the data pointed to by the node on their own.
2053 \subsection{List implementation}
2063 #include <assert.h> /* assert() */
2064 #include <stdlib.h> /* malloc(), free(), rand() */
2065 //#include <stdio.h> /* fprintf(), stdout */
2068 \subsection{List unit tests}
2070 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2072 <<list check includes>>
2074 <<main check program>>
2077 <<list check includes>>=
2078 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
2083 <<list test suite>>=
2086 <<list suite definition>>
2089 <<list suite definition>>=
2090 Suite *test_suite (void)
2092 Suite *s = suite_create ("list");
2093 <<push test case defs>>
2094 <<pop test case defs>>
2096 <<push test case adds>>
2097 <<pop test case adds>>
2103 START_TEST(test_push)
2108 push(&list, (void *)i);
2109 fail_unless(list == head(list), NULL);
2110 fail_unless( (int)list->d == 0 );
2111 for (i=0; i<n; i++) {
2112 /* we expect 0, n-1, n-2, ..., 1 */
2115 fail_unless((p=(int)pop(&list)) == e, "popped %d != %d expected", p, e);
2121 <<push test case defs>>=
2122 TCase *tc_push = tcase_create("push");
2125 <<push test case adds>>=
2126 tcase_add_test(tc_push, test_push);
2127 suite_add_tcase(s, tc_push);
2132 <<pop test case defs>>=
2134 <<pop test case adds>>=
2137 \section{Function string evaluation}
2139 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).
2140 We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator.
2141 The solution is to run a scripting language as a subprocess accessed via pipes.
2142 We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired.
2144 We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom.
2145 This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been.
2146 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.
2147 We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy.
2149 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.
2150 Then you can either statically or dynamically link to those hardcoded functions.
2151 While much less flexible, this approach would be a fairly simple-to-implement fallback.
2153 Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files.
2154 The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]].
2158 <<string eval setup declaration>>
2159 <<string eval function declaration>>
2160 <<string eval teardown declaration>>
2163 <<string eval module makefile lines>>=
2164 string_eval.c : sawsim.nw
2165 notangle -Rstring-eval.c $^ > $@
2166 string_eval.h : sawsim.nw
2167 notangle -Rstring-eval.h $^ > $@
2168 check_string_eval.c : sawsim.nw
2169 notangle -Rcheck-string-eval.c $^ > $@
2170 check_string_eval : check_string_eval.c string_eval.c string_eval.h
2171 gcc -g -o $@ $< string_eval.c -lcheck -lgsl -lgslcblas -lm
2173 rm -f string_eval.c string_eval.h check_string_eval.c check_string_eval
2176 For an introduction to POSIX process control, see\\
2177 \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\
2178 \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail),
2179 and of course, the relavent [[man]] pages.
2181 We start our subprocess with [[execvp]], one of the [[exec]] family of functions.
2182 [[execvp]] replaces the calling process' program with a new program.
2183 The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable
2184 (this may be a security hole if someone messes about with [[PATH]] before you run \prog,
2185 but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries).
2186 The new program needs command line arguments, just like it would if you were running it from a shell.
2187 The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings,
2188 with the final array entry being a [[NULL]] pointer.
2190 Now that we know how [[execvp]] works, we store it's arguments in some definitions,
2191 to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language.
2192 <<string eval subprocess definitions>>=
2193 #define SUBPROCESS "python"
2194 //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}
2195 static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL};
2196 //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL};
2198 The [[i]] option lets Python know that it should run in interactive mode.
2199 In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass.
2200 In interactive mode, python acts on every instruction as soon as it is recieved.
2201 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.
2202 %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.
2204 Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]].
2205 The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python.
2206 [[fork]] returns two copies of the same program, executing the original code.
2207 In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child.
2208 We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]].
2210 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.
2211 The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe.
2212 We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]].
2213 We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access.
2214 <<string eval pipe definitions>>=
2215 #define PIPE_READ 0 /* the end you read from */
2216 #define PIPE_WRITE 1 /* the end you write to */
2218 #define STDIN 0 /* initial index of stdin pair */
2219 #define STDOUT 2 /* initial index of stdout pair */
2222 So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations.
2224 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.
2226 <<string eval setup declaration>>=
2227 extern void string_eval_setup(FILE **pIn, FILE **pOut);
2229 <<string eval setup definition>>=
2230 void string_eval_setup(FILE **pIn, FILE **pOut)
2233 int pfd[4]; /* file descriptors for stdin and stdout */
2235 assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */
2236 assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */
2238 pid = fork(); /* split process into two copies */
2240 if (pid == -1) { /* fork error */
2241 perror("fork error");
2243 } else if (pid == 0) { /* child */
2244 close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */
2245 close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */
2246 dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */
2247 dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */
2248 execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */
2249 perror("exec error"); /* exec shouldn't return */
2251 } else { /* parent */
2252 close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */
2253 close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */
2254 *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */
2255 if ( *pIn == NULL ) {
2256 perror("fdopen (in)");
2259 *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r");
2260 if ( *pOut == NULL ) {
2261 perror("fdopen (out)");
2268 To use the evaluating subprocess, we just pipe in our command, and read out the results.
2269 For the simple cases we expect here, we restrict ourselves to a single line of returned text.
2270 <<string eval function declaration>>=
2271 extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output);
2273 <<string eval function definition>>=
2274 void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output)
2277 rc = fprintf(in, "%s", input);
2278 assert(rc == strlen(input));
2281 alarm(1); /* set a one second timeout on the read */
2282 assert( fgets(output, buflen, out) != NULL );
2283 alarm(0); /* cancel the timeout */
2284 //fprintf(stderr, "eval: %s ----> %s", input, output);
2287 The [[alarm]] calls set and clear a timeout on the returned output.
2288 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.
2289 This protects against invalid input for which a line of output is not printed to [[stdout]].
2290 Other invalid inputs (e.g. those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored.
2291 If you are getting strange results, check your python code seperately. TODO, better error handling.
2293 Cleanup is fairly straightforward, we just close the connecting streams from the parent side.
2294 With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes.
2295 The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies.
2296 As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
2297 <<string eval teardown declaration>>=
2298 extern void string_eval_teardown(FILE **pIn, FILE **pOut);
2301 <<string eval teardown definition>>=
2302 void string_eval_teardown(FILE **pIn, FILE **pOut)
2307 /* redirect Python's stderr */
2308 fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n");
2312 assert( fclose(*pIn) == 0 );
2314 assert( fclose(*pOut) == 0 );
2317 /* wait for python to exit */
2319 pid = wait(&stat_loc);
2326 if (WIFEXITED(stat_loc)) {
2327 printf("child exited with status %d\n", WEXITSTATUS(stat_loc));
2328 } else if (WIFSIGNALED(stat_loc)) {
2329 printf("child terminated with signal %d\n", WTERMSIG(stat_loc));
2334 The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
2336 \subsection{String evaluation implementation}
2340 <<string eval includes>>
2341 #include "string_eval.h"
2342 <<string eval internal definitions>>
2343 <<string eval functions>>
2346 <<string eval includes>>=
2347 #include <assert.h> /* assert() */
2348 #include <stdlib.h> /* NULL */
2349 #include <stdio.h> /* fprintf(), stdout, fdopen() */
2350 #include <string.h> /* strlen() */
2351 #include <sys/types.h> /* pid_t */
2352 #include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
2353 #include <wait.h> /* wait() */
2356 <<string eval internal definitions>>=
2357 <<string eval subprocess definitions>>
2358 <<string eval pipe definitions>>
2361 <<string eval functions>>=
2362 <<string eval setup definition>>
2363 <<string eval function definition>>
2364 <<string eval teardown definition>>
2367 \subsection{String evaluation unit tests}
2369 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
2370 <<check-string-eval.c>>=
2371 <<string eval check includes>>
2372 <<string comparison definition and globals>>
2373 <<string eval test suite>>
2374 <<main check program>>
2377 <<string eval check includes>>=
2378 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
2379 #include <stdio.h> /* printf() */
2380 #include <assert.h> /* assert() */
2381 #include <string.h> /* strlen() */
2382 #include <signal.h> /* SIGKILL */
2384 #include "string_eval.h"
2387 <<string eval test suite>>=
2388 <<string eval tests>>
2389 <<string eval suite definition>>
2392 <<string eval suite definition>>=
2393 Suite *test_suite (void)
2395 Suite *s = suite_create ("string eval");
2396 <<string eval test case defs>>
2398 <<string eval test case add>>
2403 <<string eval tests>>=
2404 START_TEST(test_setup_teardown)
2407 string_eval_setup(&in, &out);
2408 string_eval_teardown(&in, &out);
2411 START_TEST(test_invalid_command)
2414 char input[80], output[80]={};
2415 string_eval_setup(&in, &out);
2416 sprintf(input, "print ABCDefg\n");
2417 string_eval(in, out, input, 80, output);
2418 string_eval_teardown(&in, &out);
2421 START_TEST(test_simple_eval)
2424 char input[80], output[80]={};
2425 string_eval_setup(&in, &out);
2426 sprintf(input, "print 3+4*5\n");
2427 string_eval(in, out, input, 80, output);
2428 fail_unless(STRMATCH(output,"23\n"), NULL);
2429 string_eval_teardown(&in, &out);
2432 START_TEST(test_multiple_evals)
2435 char input[80], output[80]={};
2436 string_eval_setup(&in, &out);
2437 sprintf(input, "print 3+4*5\n");
2438 string_eval(in, out, input, 80, output);
2439 fail_unless(STRMATCH(output,"23\n"), NULL);
2440 sprintf(input, "print (3**2 + 4**2)**0.5\n");
2441 string_eval(in, out, input, 80, output);
2442 fail_unless(STRMATCH(output,"5.0\n"), NULL);
2443 string_eval_teardown(&in, &out);
2446 START_TEST(test_eval_with_state)
2449 char input[80], output[80]={};
2450 string_eval_setup(&in, &out);
2451 sprintf(input, "print 3+4*5\n");
2452 fprintf(in, "A = 3\n");
2453 sprintf(input, "print A*3\n");
2454 string_eval(in, out, input, 80, output);
2455 fail_unless(STRMATCH(output,"9\n"), NULL);
2456 string_eval_teardown(&in, &out);
2460 <<string eval test case defs>>=
2461 TCase *tc_string_eval = tcase_create("string_eval");
2463 <<string eval test case add>>=
2464 tcase_add_test(tc_string_eval, test_setup_teardown);
2465 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
2466 tcase_add_test(tc_string_eval, test_simple_eval);
2467 tcase_add_test(tc_string_eval, test_multiple_evals);
2468 tcase_add_test(tc_string_eval, test_eval_with_state);
2469 suite_add_tcase(s, tc_string_eval);
2473 \section{Accelerating function evaluation}
2475 My first version-0.3 code was running very slowly.
2476 With the options suggested in the help
2477 ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]),
2478 the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]],
2479 making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions.
2480 That's only 15 calls per solution, so the search algorithm seems reasonable.
2481 The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table.
2484 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
2488 <<accel k module makefile lines>>=
2489 accel_k.c : sawsim.nw
2490 notangle -Raccel-k.c $^ > $@
2491 accel_k.h : sawsim.nw
2492 notangle -Raccel-k.h $^ > $@
2493 check_accel_k.c : sawsim.nw
2494 notangle -Rcheck-accel_k.c $^ > $@
2495 check_accel_k : check_accel_k.c global.h
2496 gcc -g -o $@ $< accel_k.c -lcheck -lgsl -lgslcblas -lm
2498 rm -f accel_k.c accel_k.h check_accel_k.c check_accel_k
2502 #include <assert.h> /* assert() */
2503 #include <stdlib.h> /* realloc(), free(), NULL */
2504 #include "global.h" /* environment_t */
2505 #include "k_model.h" /* k_func_t */
2506 #include "interp.h" /* interp_* */
2507 #include "accel_k.h"
2509 typedef struct accel_k_struct {
2510 interp_table_t *itable;
2516 /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */
2517 static int num_accels = 0;
2518 static accel_k_t *accels=NULL;
2520 /* Wrap k in the standard f(x) acceleration form */
2521 static double k_wrap(double F, void *params)
2523 accel_k_t *p = (accel_k_t *)params;
2524 return (*p->k)(F, p->env, p->params);
2527 static int k_tol(double FA, double kA, double FB, double kB)
2530 if (FB-FA > 1e-12) {
2531 //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB));
2532 return 1; /* unnacceptable */
2534 //printf("acceptable tol\n");
2535 return 0; /* acceptable */
2539 static int add_accel_k(k_func_t *k, environment_t *env, void *params)
2542 accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
2543 assert(accels != NULL);
2544 accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
2546 accels[i].env = env;
2547 accels[i].params = params;
2554 for (i=0; i<num_accels; i++)
2555 interp_table_free(accels[i].itable);
2557 if (accels) free(accels);
2561 double accel_k(k_func_t *k, double F, environment_t *env, void *params)
2564 for (i=0; i<num_accels; i++) {
2565 if (accels[i].k == k && accels[i].env == env && accels[i].params == params) {
2566 /* We've already setup this function.
2567 * WARNING: we're assuming param and env are the same. */
2568 return interp_table_eval(accels[i].itable, accels+i, F);
2571 /* set up a new acceleration environment */
2572 i = add_accel_k(k, env, params);
2573 return interp_table_eval(accels[i].itable, accels+i, F);
2577 \section{Tension models}
2578 \label{sec.tension_models}
2580 TODO: tension model intro.
2581 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.
2583 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
2584 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]].
2586 <<tension-model.h>>=
2588 <<tension handler types>>
2589 <<find tension definitions>>
2590 <<constant tension model declarations>>
2591 <<hooke tension model declarations>>
2592 <<worm-like chain tension model declarations>>
2593 <<freely-jointed chain tension model declarations>>
2596 <<tension model module makefile lines>>=
2597 tension_model.c : sawsim.nw
2598 notangle -Rtension-model.c $^ > $@
2599 tension_model.h : sawsim.nw
2600 notangle -Rtension-model.h $^ > $@
2601 check_tension_model.c : sawsim.nw
2602 notangle -Rcheck-tension-model.c $^ > $@
2603 check_tension_model : check_tension_model.c global.h tension_model.c tension_model.h
2604 gcc -g -o $@ $< tension_model.c -lcheck -lgsl -lgslcblas -lm
2605 clean_tension_model : clean_tension_model_utils
2606 rm -f tension_model.c tension_model.h check_tension_model.c check_tension_model
2607 tension_model_utils.c : sawsim.nw
2608 notangle -Rtension-model-utils.c $^ > $@
2609 tension_model_utils : tension_model_utils.c global.h tension_model.c tension_model.h parse.c parse.h \
2610 list.c list.h tension_balance.c tension_balance.h
2611 gcc -g -o $@ $< tension_model.c parse.c list.c tension_balance.c -lgsl -lgslcblas -lm
2612 tension_model_utils_static : tension_model_utils.c global.h tension_model.c tension_model.h parse.c parse.h \
2613 list.c list.h tension_balance.c tension_balance.h
2614 gcc -g -static -o $@ $< tension_model.c parse.c list.c tension_balance.c -lgsl -lgslcblas -lm
2615 clean_tension_model_utils :
2616 rm -f tension_model_utils.c tension_model_utils
2620 \label{sec.null_tension}
2622 For unstretchable domains.
2624 <<null tension model getopt>>=
2625 {TG_NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
2628 \subsection{Constant}
2629 \label{sec.const_tension}
2631 <<constant tension functions>>=
2632 <<constant tension function>>
2633 <<constant tension parameter create/destroy functions>>
2636 <<constant tension model declarations>>=
2637 <<constant tension function declaration>>
2638 <<constant tension parameter create/destroy function declarations>>
2639 <<constant tension model global declarations>>
2643 An infinitely stretchable domain providing a constant tension.
2645 <<constant tension function declaration>>=
2646 extern double const_tension_handler(double x, void *pdata);
2648 <<constant tension function>>=
2649 double const_tension_handler(double x, void *pdata)
2651 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
2652 list_t *list = data->group;
2657 assert(list != NULL); /* empty active group?! */
2658 F = ((const_tension_param_t *)list->d)->F;
2659 while (list != NULL) {
2660 assert(((const_tension_param_t *)list->d)->F == F);
2668 <<constant tension structure>>=
2669 typedef struct const_tension_param_struct {
2670 double F; /* tension (force) in N */
2671 } const_tension_param_t;
2675 <<constant tension parameter create/destroy function declarations>>=
2676 extern void *string_create_const_tension_param_t(char **param_strings);
2677 extern void destroy_const_tension_param_t(void *p);
2680 <<constant tension parameter create/destroy functions>>=
2681 const_tension_param_t *create_const_tension_param_t(double F)
2683 const_tension_param_t *ret
2684 = (const_tension_param_t *) malloc(sizeof(const_tension_param_t));
2685 assert(ret != NULL);
2690 void *string_create_const_tension_param_t(char **param_strings)
2692 return create_const_tension_param_t(atof(param_strings[0]));
2695 void destroy_const_tension_param_t(void *p)
2703 <<constant tension model global declarations>>=
2704 extern int num_const_tension_params;
2705 extern char *const_tension_param_descriptions[];
2706 extern char const_tension_param_string[];
2709 <<constant tension model globals>>=
2710 int num_const_tension_params = 1;
2711 char *const_tension_param_descriptions[] = {"tension F, N"};
2712 char const_tension_param_string[] = "0";
2715 <<constant tension model getopt>>=
2716 {TG_CONST, "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}
2722 <<hooke functions>>=
2724 <<hooke parameter create/destroy functions>>
2727 <<hooke tension model declarations>>=
2728 <<hooke tension function declaration>>
2729 <<hooke tension parameter create/destroy function declarations>>
2730 <<hooke tension model global declarations>>
2734 The tension of a single spring is given by $F=kx$ for some spring constant $k$.
2735 The behavior of a series of springs $k_i$ in series is given by
2737 F = \p({\sum_i \frac{1}{k_i}})^{-1} x
2739 For a simple proof, see Appendix \ref{app.math_hooke}.
2741 <<hooke tension function declaration>>=
2742 extern double hooke_handler(double x, void *pdata);
2746 double hooke_handler(double x, void *pdata)
2748 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
2749 list_t *list = data->group;
2754 assert(list != NULL); /* empty active group?! */
2755 while (list != NULL) {
2756 assert( ((hooke_param_t *)list->d)->k > 0 );
2757 k += 1.0/ ((hooke_param_t *)list->d)->k;
2765 <<hooke structure>>=
2766 typedef struct hooke_param_struct {
2767 double k; /* spring constant in N/m */
2771 <<hooke tension parameter create/destroy function declarations>>=
2772 extern void *string_create_hooke_param_t(char **param_strings);
2773 extern void destroy_hooke_param_t(void *p);
2776 <<hooke parameter create/destroy functions>>=
2777 hooke_param_t *create_hooke_param_t(double k)
2779 hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t));
2780 assert(ret != NULL);
2785 void *string_create_hooke_param_t(char **param_strings)
2787 return create_hooke_param_t(atof(param_strings[0]));
2790 void destroy_hooke_param_t(void *p)
2797 <<hooke tension model global declarations>>=
2798 extern int num_hooke_params;
2799 extern char *hooke_param_descriptions[];
2800 extern char hooke_param_string[];
2803 <<hooke tension model globals>>=
2804 int num_hooke_params = 1;
2805 char *hooke_param_descriptions[] = {"spring constant k, N/m"};
2806 char hooke_param_string[]="0.05";
2809 <<hooke tension model getopt>>=
2810 {TG_HOOKE, "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}
2813 \subsection{Worm-like chain}
2816 We can model several unfolded domains with a single worm-like chain.
2817 <<worm-like chain functions>>=
2818 <<worm-like chain function>>
2819 <<worm-like chain handler>>
2820 <<worm-like chain create/destroy functions>>
2823 <<worm-like chain tension model declarations>>=
2824 <<worm-like chain handler declaration>>
2825 <<worm-like chain create/destroy function declarations>>
2826 <<worm-like chain tension model global declarations>>
2830 The combination of all unfolded domains is modeled as a worm like chain,
2831 with the force $F_{\text{WLC}}$ approximately given by
2833 F_{\text{WLC}} \approx \frac{k_B T}{p}
2834 \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}],
2836 where $x$ is the distance between the fixed ends,
2837 $k_B$ is Boltzmann's constant,
2838 $T$ is the absolute temperature,
2839 $p$ is the persistence length, and
2840 $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
2841 <<worm-like chain function>>=
2842 double wlc(double x, double T, double p, double L)
2844 double xL=0.0; /* uses global k_B */
2846 assert(T > 0); assert(p > 0); assert(L > 0);
2847 if (x >= L) return HUGE_VAL;
2848 xL = x/L; /* unitless */
2849 //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L,
2850 // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ));
2851 return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL );
2854 This model assumes that each unfolded domain has the same persistence length.
2856 <<worm-like chain handler declaration>>=
2857 extern double wlc_handler(double x, void *pdata);
2860 <<worm-like chain handler>>=
2861 double wlc_handler(double x, void *pdata)
2863 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
2864 list_t *list = data->group;
2868 assert(list != NULL); /* empty active group?! */
2869 p = ((wlc_param_t *)list->d)->p;
2870 while (list != NULL) {
2871 /* enforce consistent p */
2872 assert( ((wlc_param_t *)list->d)->p == p);
2873 L += ((wlc_param_t *)list->d)->L;
2876 return wlc(x, data->env->T, p, L);
2880 <<worm-like chain structure>>=
2881 typedef struct wlc_param_struct {
2882 double p; /* WLC persistence length (m) */
2883 double L; /* WLC contour length (m) */
2887 <<worm-like chain create/destroy function declarations>>=
2888 extern void *string_create_wlc_param_t(char **param_strings);
2889 extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
2892 <<worm-like chain create/destroy functions>>=
2893 wlc_param_t *create_wlc_param_t(double p, double L)
2895 wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t));
2896 assert(ret != NULL);
2902 void *string_create_wlc_param_t(char **param_strings)
2904 return create_wlc_param_t(atof(param_strings[0]), atof(param_strings[1]));
2907 void destroy_wlc_param_t(void *p /* wlc_param_t * */)
2914 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.
2915 TODO (now we copy the string before parsing, but still do this for clarity).
2917 <<valid string write code>>=
2918 char string[] = "abc";
2921 <<valid string write code>>=
2922 char *string = "abc";
2926 <<worm-like chain tension model global declarations>>=
2927 extern int num_wlc_params;
2928 extern char *wlc_param_descriptions[];
2929 extern char wlc_param_string[];
2931 <<worm-like chain tension model globals>>=
2932 int num_wlc_params = 2;
2933 char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"};
2934 char wlc_param_string[]="0.39e-9,28.4e-9";
2937 <<worm-like chain tension model getopt>>=
2938 {TG_WLC, "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}
2940 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
2942 \subsection{Freely-jointed chain}
2945 <<freely-jointed chain functions>>=
2946 <<freely-jointed chain function>>
2947 <<freely-jointed chain handler>>
2948 <<freely-jointed chain create/destroy functions>>
2951 <<freely-jointed chain tension model declarations>>=
2952 <<freely-jointed chain handler declaration>>
2953 <<freely-jointed chain create/destroy function declarations>>
2954 <<freely-jointed chain tension model global declarations>>
2958 The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints.
2959 The tension of a chain of $N$ such links is given by
2961 F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;,
2963 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}.
2964 <<freely-jointed chain function>>=
2965 double langevin(double x, void *params)
2968 return 1.0/tanh(x) - 1/x;
2970 <<one dimensional bracket declaration>>
2971 <<one dimensional solve declaration>>
2972 double inv_langevin(double x)
2974 double lb=0.0, ub=1.0, ret;
2975 oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
2976 ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */
2980 double fjc(double x, double T, double l, int N)
2982 double L = l*(double)N;
2983 if (x == 0) return 0;
2984 //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
2985 return k_B*T/l * inv_langevin(x/L);
2989 <<freely-jointed chain handler declaration>>=
2990 extern double fjc_handler(double x, void *pdata);
2992 <<freely-jointed chain handler>>=
2993 double fjc_handler(double x, void *pdata)
2995 tension_handler_data_t *data = (tension_handler_data_t *)pdata;
2996 list_t *list = data->group;
3001 assert(list != NULL); /* empty active group?! */
3002 l = ((fjc_param_t *)list->d)->l;
3003 while (list != NULL) {
3004 /* enforce consistent l */
3005 assert( ((fjc_param_t *)list->d)->l == l);
3006 N += ((fjc_param_t *)list->d)->N;
3009 return fjc(x, data->env->T, l, N);
3013 <<freely-jointed chain structure>>=
3014 typedef struct fjc_param_struct {
3015 double l; /* FJC link length (m) */
3016 int N; /* FJC number of links */
3020 <<freely-jointed chain create/destroy function declarations>>=
3021 extern void *string_create_fjc_param_t(char **param_strings);
3022 extern void destroy_fjc_param_t(void *p);
3025 <<freely-jointed chain create/destroy functions>>=
3026 fjc_param_t *create_fjc_param_t(double l, double N)
3028 fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t));
3029 assert(ret != NULL);
3035 void *string_create_fjc_param_t(char **param_strings)
3037 return create_fjc_param_t(atof(param_strings[0]), atof(param_strings[1]));
3040 void destroy_fjc_param_t(void *p)
3047 <<freely-jointed chain tension model global declarations>>=
3048 extern int num_fjc_params;
3049 extern char *fjc_param_descriptions[];
3050 extern char fjc_param_string[];
3053 <<freely-jointed chain tension model globals>>=
3054 int num_fjc_params = 2;
3055 char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"};
3056 char fjc_param_string[]="0.5e-9,200";
3059 <<freely-jointed chain tension model getopt>>=
3060 {TG_FJC, "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}
3063 \subsection{Tension model implementation}
3065 <<tension-model.c>>=
3067 <<tension model includes>>
3068 #include "tension_model.h"
3069 <<tension model internal definitions>>
3070 <<tension model globals>>
3071 <<tension model functions>>
3074 <<tension model includes>>=
3075 #include <assert.h> /* assert() */
3076 #include <stdlib.h> /* NULL */
3077 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
3078 #include <stdio.h> /* fprintf(), stdout */
3081 #include "tension_balance.h" /* oneD_*() */
3084 <<tension model internal definitions>>=
3085 <<constant tension structure>>
3087 <<worm-like chain structure>>
3088 <<freely-jointed chain structure>>
3091 <<tension model globals>>=
3092 <<constant tension model globals>>
3093 <<hooke tension model globals>>
3094 <<worm-like chain tension model globals>>
3095 <<freely-jointed chain tension model globals>>
3098 <<tension model functions>>=
3099 <<constant tension functions>>
3101 <<worm-like chain functions>>
3102 <<freely-jointed chain functions>>
3106 \subsection{Utilities}
3108 The tension models can get complicated, and users may want to reassure
3109 themselves that this portion of the input physics are functioning properly. The
3110 stand-alone program [[tension_model_utils]] demonstrates the tension model
3111 interface without the overhead of having to understand a full simulation.
3112 [[tension_model_utils]] takes command line model arguments like the full
3113 simulation, and prints $F(x)$ data to the screen for the selected model over a
3116 <<tension-model-utils.c>>=
3118 <<tension model utility includes>>
3119 <<tension model utility definitions>>
3120 <<tension model utility getopt functions>>
3122 <<tension model utility main>>
3125 <<tension model utility main>>=
3126 int main(int argc, char **argv)
3128 <<tension model getopt array definition>>
3129 tension_model_getopt_t *model = NULL;
3133 one_dim_func_t *tension_handler[NUM_TENSION_GROUPS] = {0};
3134 tension_handler_data_t tdata;
3135 int num_param_args; /* for INIT_MODEL() */
3136 char **param_args; /* for INIT_MODEL() */
3137 get_options(argc, argv, &env, NUM_TENSION_GROUPS, tension_models, &model, &flags);
3138 setup(tension_handler);
3139 if (flags & VFLAG) {
3140 printf("#initializing model %s with parameters %s\n", model->name, model->params);
3142 INIT_MODEL("utility", model, params);
3144 push(&tdata.group, params);
3146 tdata.persist = NULL;
3147 if (tension_handler[model->tg_group] == NULL) {
3148 printf("No tension function for model %s\n", model->name);
3152 double dx=1e-10, x=0, F=0;
3153 printf("#F (N)\tk (%% pop. per s)\n");
3154 while (F >= 0 && F < 1e5 && x < 1e-6) {
3155 F = (*tension_handler[model->tg_group])(x, &tdata);
3156 printf("%g\t%g\n", x, F);
3160 params = pop(&tdata.group);
3161 if (model->destructor)
3162 (*model->destructor)(params);
3167 <<tension model utility includes>>=
3170 #include <string.h> /* strlen() */
3171 #include <assert.h> /* assert() */
3175 #include "tension_model.h"
3178 <<tension model utility definitions>>=
3179 <<version definition>>
3180 #define VFLAG 1 /* verbose */
3181 <<string comparison definition and globals>>
3182 <<tension model getopt definitions>>
3183 <<initialize model definition>>
3187 <<tension model utility getopt functions>>=
3188 <<index tension model>>
3189 <<tension model utility help>>
3190 <<tension model utility get options>>
3193 <<tension model utility help>>=
3194 void help(char *prog_name,
3196 int n_tension_models, tension_model_getopt_t *tension_models,
3200 printf("usage: %s [options]\n", prog_name);
3201 printf("Version %s\n\n", VERSION);
3202 printf("A utility for understanding the available tension models\n\n");
3203 printf("Environmental options:\n");
3204 printf("-T T\tTemperature T (currently %g K)\n", env->T);
3205 printf("-C T\tYou can also set the temperature T in Celsius\n");
3206 printf("Model options:\n");
3207 printf("-m model\tFolded tension model (currently %s)\n",
3208 tension_models[tension_model].name);
3209 printf("-a args\tFolded tension model argument string (currently %s)\n",
3210 tension_models[tension_model].params);
3211 printf("F(x) is calculated for a range of x and printed\n");
3212 printf("For example:\n");
3213 printf(" #Distance (x)\tForce (N)\n");
3214 printf(" 123.456\t7.89\n");
3216 printf("-V\tChange output to verbose mode\n");
3217 printf("-h\tPrint this help and exit\n");
3219 printf("Tension models:\n");
3220 for (i=0; i<n_tension_models; i++) {
3221 printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
3222 for (j=0; j < tension_models[i].num_params ; j++)
3223 printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
3224 printf("\t\tdefaults: %s\n", tension_models[i].params);
3229 <<tension model utility get options>>=
3230 void get_options(int argc, char **argv, environment_t *env,
3231 int n_tension_models, tension_model_getopt_t *tension_models,
3232 tension_model_getopt_t **model,
3233 unsigned int *flags)
3235 char *prog_name = NULL;
3236 char c, options[] = "T:C:m:a:Vh";
3237 int tension_model=0;
3238 extern char *optarg;
3239 extern int optind, optopt, opterr;
3243 /* setup defaults */
3245 prog_name = argv[0];
3246 env->T = 300.0; /* K */
3248 *model = tension_models;
3250 while ((c=getopt(argc, argv, options)) != -1) {
3252 case 'T': env->T = atof(optarg); break;
3253 case 'C': env->T = atof(optarg)+273.15; break;
3255 tension_model = index_tension_model(n_tension_models, tension_models, optarg);
3256 *model = tension_models+tension_model;
3259 tension_models[tension_model].params = optarg;
3261 case 'V': *flags |= VFLAG; break;
3263 fprintf(stderr, "unrecognized option '%c'\n", optopt);
3264 /* fall through to default case */
3266 help(prog_name, env, n_tension_models, tension_models, tension_model);
3275 \section{Unfolding rate models}
3276 \label{sec.k_models}
3278 We have two main choices for determining $k$: Bell's model, which assumes the
3279 folded domains exist in a single `folded' state and make instantaneous,
3280 stochastic transitions over a free energy barrier; and Kramers' model, which
3281 treats the unfolding as a thermalized diffusion process.
3282 We deal with these options like we dealt with the different tension models: by bundling all model-specific
3283 parameters into structures, with model specific functions for determining $k$.
3285 <<k func definition>>=
3286 typedef double k_func_t(double F, environment_t *env, void *params);
3289 Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files.
3290 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]].
3292 We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$,
3293 because \LaTeX\ doesn't like underscores outside math mode.
3294 Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+').
3295 You could use spaces instead (`\verb+ +'),
3296 but then your [[notangle]]s would look like [[notangle -R'k model.h']],
3297 which I don't like as much.
3301 <<k func definition>>
3302 <<null k declarations>>
3303 <<const k declarations>>
3304 <<bell k declarations>>
3305 <<kramers k declarations>>
3306 <<kramers integ k declarations>>
3309 <<k model module makefile lines>>=
3310 k_model.c : sawsim.nw
3311 notangle -Rk-model.c $^ > $@
3312 k_model.h : sawsim.nw
3313 notangle -Rk-model.h $^ > $@
3314 check_k_model.c : sawsim.nw
3315 notangle -Rcheck-k-model.c $^ > $@
3316 check_k_model : check_k_model.c global.h k_model.c k_model.h
3317 gcc -g -o $@ $< k_model.c -lcheck -lgsl -lgslcblas -lm
3318 clean_k_model : clean_k_model_utils
3319 rm -f k_model.c k_model.h check_k_model.c check_k_model
3320 k_model_utils.c : sawsim.nw
3321 notangle -Rk-model-utils.c $^ > $@
3322 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
3323 gcc -g -o $@ $< k_model.c parse.c string_eval.c -lgsl -lgslcblas -lm
3324 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
3325 gcc -g -static -o $@ $< k_model.c parse.c string_eval.c -lgsl -lgslcblas -lm
3326 clean_k_model_utils :
3327 rm -f k_model_utils.c k_model_utils
3333 For domains that are never folded.
3335 <<null k declarations>>=
3337 <<null k functions>>=
3342 <<null k model getopt>>=
3343 {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
3346 \subsection{Constant rate model}
3349 This is a pretty boring model, but it is usefull for testing the engine.
3353 where $k_0$ is the constant unfolding rate.
3355 <<const k functions>>=
3356 <<const k function>>
3357 <<const k structure create/destroy functions>>
3360 <<const k declarations>>=
3361 <<const k function declaration>>
3362 <<const k structure create/destroy function declarations>>
3363 <<const k global declarations>>
3367 <<const k structure>>=
3368 typedef struct const_k_param_struct {
3369 double knot; /* transition rate k_0 (frac population per s) */
3373 <<const k function declaration>>=
3374 double const_k(double F, environment_t *env, void *const_k_params);
3376 <<const k function>>=
3377 double const_k(double F, environment_t *env, void *const_k_params)
3378 { /* Returns the rate constant k in frac pop per s. */
3379 const_k_param_t *p = (const_k_param_t *) const_k_params;
3381 assert(p->knot > 0);
3386 <<const k structure create/destroy function declarations>>=
3387 void *string_create_const_k_param_t(char **param_strings);
3388 void destroy_const_k_param_t(void *p);
3391 <<const k structure create/destroy functions>>=
3392 const_k_param_t *create_const_k_param_t(double knot)
3394 const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t));
3395 assert(ret != NULL);
3400 void *string_create_const_k_param_t(char **param_strings)
3402 return create_const_k_param_t(atof(param_strings[0]));
3405 void destroy_const_k_param_t(void *p /* const_k_param_t * */)
3412 <<const k global declarations>>=
3413 extern int num_const_k_params;
3414 extern char *const_k_param_descriptions[];
3415 extern char const_k_param_string[];
3418 <<const k globals>>=
3419 int num_const_k_params = 1;
3420 char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"};
3421 char const_k_param_string[]="1";
3424 <<const k model getopt>>=
3425 {"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}
3428 \subsection{Bell's model}
3431 Quantitatively, Bell's model gives $k$ as
3433 k = k_0 \cdot e^{F dx / k_B T} \;,
3435 where $k_0$ is the unforced unfolding rate,
3436 $F$ is the applied unfolding force,
3437 $dx$ is the distance to the transition state, and
3438 $k_B T$ is the thermal energy\citep{schlierf06}.
3440 <<bell k functions>>=
3442 <<bell k structure create/destroy functions>>
3445 <<bell k declarations>>=
3446 <<bell k function declaration>>
3447 <<bell k structure create/destroy function declarations>>
3448 <<bell k global declarations>>
3452 <<bell k structure>>=
3453 typedef struct bell_param_struct {
3454 double knot; /* transition rate k_0 (frac population per s) */
3455 double dx; /* `distance to transition state' (nm) */
3459 <<bell k function declaration>>=
3460 double bell_k(double F, environment_t *env, void *bell_params);
3462 <<bell k function>>=
3463 double bell_k(double F, environment_t *env, void *bell_params)
3464 { /* Returns the rate constant k in frac pop per s.
3465 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
3466 * uses global k_B in J/K */
3467 bell_param_t *p = (bell_param_t *) bell_params;
3468 assert(F >= 0); assert(env->T > 0);
3470 assert(p->knot > 0); assert(p->dx >= 0);
3472 printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n",
3473 F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T),
3474 p->knot * exp(F*p->dx / (k_B*env->T) ));
3475 printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T);
3477 return p->knot * exp(F*p->dx / (k_B*env->T) );
3481 <<bell k structure create/destroy function declarations>>=
3482 void *string_create_bell_param_t(char **param_strings);
3483 void destroy_bell_param_t(void *p);
3486 <<bell k structure create/destroy functions>>=
3487 bell_param_t *create_bell_param_t(double knot, double dx)
3489 bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t));
3490 assert(ret != NULL);
3496 void *string_create_bell_param_t(char **param_strings)
3498 return create_bell_param_t(atof(param_strings[0]), atof(param_strings[1]));
3501 void destroy_bell_param_t(void *p /* bell_param_t * */)
3508 <<bell k global declarations>>=
3509 extern int num_bell_params;
3510 extern char *bell_param_descriptions[];
3511 extern char bell_param_string[];
3515 int num_bell_params = 2;
3516 char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"};
3517 char bell_param_string[]="3.3e-4,0.25e-9";
3520 <<bell k model getopt>>=
3521 {"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}
3523 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
3524 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
3526 \subsection{Kramer's model}
3529 Kramer's model gives $k$ as
3531 % \frac{1}{k} = \frac{1}{D}
3532 % \int_{x_\text{min}}^{x_\text{max}}
3533 % e^{\frac{-U_F(x)}{k_B T}}
3534 % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
3537 %where $D$ is the diffusion constant,
3538 %$U_F(x)$ is the free energy along the unfolding coordinate, and
3539 %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
3540 % before the minimum and after the maximum, respectively \citep{schlierf06}.
3542 k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T},
3544 where $D$ is the diffusion constant,
3546 l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}},
3548 are length scales where
3549 $x_c(F)$ is the location of the bound state and
3550 $x_{ts}(F)$ is the location of the transition state,
3551 $E(F,x)$ is the free energy, and
3552 $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy
3553 (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanngi90} Eqn.~4.56c,
3554 \citet{evans97} Eqn.~3).
3556 <<kramers k functions>>=
3557 <<kramers k function>>
3558 <<kramers k structure create/destroy functions>>
3561 <<kramers k declarations>>=
3562 <<kramers k function declaration>>
3563 <<kramers k structure create/destroy function declarations>>
3564 <<kramers k global declarations>>
3568 <<kramers k structure>>=
3569 //typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
3570 //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x);
3572 typedef struct kramers_param_struct {
3573 double D; /* diffusion rate (um/s) */
3580 //kramers_x_func_t *xc; /* function returning metastable x (nm) */
3581 //kramers_x_func_t *xts; /* function returning transition x (nm) */
3582 //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */
3583 //kramers_E_func_t *E; /* function returning E (J) */
3584 //double *E_params; /* parametrize the energy landscape */
3585 //int n_E_params; /* length of E_params array */
3589 <<kramers k function declaration>>=
3590 extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */
3591 extern double kramers_k(double F, environment_t *env, void *kramers_params);
3593 <<kramers k function>>=
3594 double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
3596 static char input[160]={0};
3597 static char output[80]={0};
3598 /* setup the environment */
3599 fprintf(in, "F = %g; T = %g\n", F, T);
3600 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
3601 string_eval(in, out, input, 80, output);
3602 return atof(output);
3605 double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
3607 static char input[160]={0};
3608 static char output[80]={0};
3609 /* setup the environment */
3610 fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x);
3611 sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */
3612 string_eval(in, out, input, 80, output);
3613 return atof(output);
3616 double kramers_E(double F, environment_t *env, void *kramers_params, double x)
3618 kramers_param_t *p = (kramers_param_t *) kramers_params;
3619 return kramers_E_fn(p->in, p->out, p->E, F, env->T, x);
3622 double kramers_k(double F, environment_t *env, void *kramers_params)
3623 { /* Returns the rate constant k in frac pop per s.
3624 * Takes F in N, T in K, knot in frac pop per s, and dx in m.
3625 * uses global k_B in J/K */
3626 kramers_param_t *p = (kramers_param_t *) kramers_params;
3627 double kbT, xc, xts, lc, lts, Eb;
3628 assert(F >= 0); assert(env->T > 0);
3631 assert(p->xc != NULL); assert(p->xts != NULL);
3632 assert(p->ddEdxx != NULL); assert(p->E != NULL);
3633 assert(p->in != NULL); assert(p->out != NULL);
3635 xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T);
3636 if (xc == -1.0) return -1.0; /* error (Too much force) */
3637 xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T);
3638 if (xc == -1.0) return -1.0; /* error (Too much force) */
3639 lc = sqrt(2.0*M_PI*kbT /
3640 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc));
3641 lts = sqrt(-2.0*M_PI*kbT/
3642 kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts));
3643 Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts)
3644 - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc);
3645 //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));
3646 return p->D/(lc*lts) * exp(-Eb/kbT);
3650 <<kramers k structure create/destroy function declarations>>=
3651 void *string_create_kramers_param_t(char **param_strings);
3652 void destroy_kramers_param_t(void *p);
3655 <<kramers k structure create/destroy functions>>=
3656 kramers_param_t *create_kramers_param_t(double D,
3657 char *xc_fn, char *xts_fn,
3658 char *E_fn, char *ddEdxx_fn,
3659 char *global_define_string)
3660 // kramers_x_func_t *xc, kramers_x_func_t *xts,
3661 // kramers_E_func_t *ddEdxx, kramers_E_func_t *E,
3662 // double *E_params, int n_E_params)
3664 kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t));
3665 assert(ret != NULL);
3666 ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
3667 assert(ret->xc != NULL);
3668 ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1));
3669 assert(ret->xts != NULL);
3670 ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1));
3671 assert(ret->E != NULL);
3672 ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1));
3673 assert(ret->ddEdxx != NULL);
3675 strcpy(ret->xc, xc_fn);
3676 strcpy(ret->xts, xts_fn);
3677 strcpy(ret->E, E_fn);
3678 strcpy(ret->ddEdxx, ddEdxx_fn);
3679 //ret->E_params = E_params;
3680 //ret->n_E_params = n_E_params;
3681 string_eval_setup(&ret->in, &ret->out);
3682 fprintf(ret->in, "kB = %g\n", k_B);
3683 fprintf(ret->in, "%s\n", global_define_string);
3687 /* takes an array of 4 functions: {"(1,2),(2,0),..."} */
3688 void *string_create_kramers_param_t(char **param_strings)
3690 return create_kramers_param_t(atof(param_strings[0]),
3698 void destroy_kramers_param_t(void *p /* kramers_param_t * */)
3700 kramers_param_t *pk = (kramers_param_t *)p;
3702 string_eval_teardown(&pk->in, &pk->out);
3704 // free(pk->E_params);
3710 For our defaults, we'll follow \citet{schlierf06} and study ddFLN4.
3711 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.
3712 We follow \citet{shillcock98} and use
3714 E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b
3715 \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}})
3718 where TODO, variable meanings.%\citep{shillcock98}.
3719 The first and second derivatives of $E(F,E_b,x,x_s)$ are
3721 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\\
3722 E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}).
3724 $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$,
3725 $b=\frac{27 E_b}{x_s^2}$, and $c=-F$.
3726 The roots are therefor at
3728 x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\
3729 &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\
3730 &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}})
3733 <<kramers k global declarations>>=
3734 extern int num_kramers_params;
3735 extern char *kramers_param_descriptions[];
3736 extern char kramers_param_string[];
3739 <<kramers k globals>>=
3740 int num_kramers_params = 6;
3741 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"};
3742 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)}";
3744 Which we initialized with parameters for ddFLN4\citep{schlierf06}.
3745 There is a bit of funny buisness in the $x_c$ definition to handle high force cases.
3746 When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down.
3747 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.
3748 The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]].
3749 It works so long as [[val_1]] is not [[false]].
3751 <<kramers k model getopt>>=
3752 {"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}
3755 \subsection{Kramer's model (natural cubic splines)}
3756 \label{sec.kramers_integ}
3758 Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model,
3759 (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1,
3760 \citet{schlierf06) Eqn.~2)
3762 \frac{1}{k} = \frac{1}{D}
3763 \int_{x_\text{min}}^{x_\text{max}}
3764 e^{\frac{U_F(x)}{k_B T}}
3765 \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
3768 where $D$ is the diffusion constant,
3769 $U_F(x)$ is the free energy along the unfolding coordinate, and
3770 $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$
3771 before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}.
3773 We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points),
3774 interpolating between them with cubic splines.
3775 We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here.
3777 <<kramers integ k functions>>=
3778 <<spline functions>>
3779 <<kramers integ A k functions>>
3780 <<kramers integ B k functions>>
3781 <<kramers integ k function>>
3782 <<kramers integ k structure create/destroy functions>>
3785 <<kramers integ k declarations>>=
3786 <<kramers integ k function declaration>>
3787 <<kramers integ k structure create/destroy function declarations>>
3788 <<kramers integ k global declarations>>
3792 <<kramers integ k structure>>=
3793 typedef double func_t(double x, void *params);
3794 typedef struct kramers_integ_param_struct {
3795 double D; /* diffusion rate (um/s) */
3796 double x_min; /* integration bounds */
3798 func_t *E; /* function returning E (J) */
3799 void *E_params; /* parametrize the energy landscape */
3800 destroy_data_func_t *destroy_E_params;
3802 double F; /* for passing into GSL evaluations */
3804 } kramers_integ_param_t;
3807 <<spline param structure>>=
3808 typedef struct spline_param_struct {
3809 int n_knots; /* natural cubic spline knots for U(x) */
3812 gsl_spline *spline; /* GSL spline parameters */
3813 gsl_interp_accel *acc; /* GSL spline acceleration data */
3817 We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting.
3819 \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx'
3821 <<kramers integ A k functions>>=
3822 double kramers_integ_k_integrandA(double x, void *params)
3824 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
3825 double U = (*p->E)(x, p->E_params);
3826 double Fx = p->F * x;
3827 double kbT = k_B * p->env->T;
3828 //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT));
3829 return exp(-(U-Fx)/kbT);
3831 double kramers_integ_k_integralA(double x, void *params)
3833 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
3835 double result, abserr;
3836 size_t neval; /* for qng */
3837 static gsl_integration_workspace *w = NULL;
3838 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
3839 f.function = &kramers_integ_k_integrandA;
3841 //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
3842 assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
3843 //fprintf(stderr, "integralA = %g\n", result);
3849 \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
3850 \int_{x_\text{min}}^{x_\text{max}}
3851 e^{\frac{U_F(x)}{k_B T}}
3852 \text{Integral}_A(x)
3855 <<kramers integ B k functions>>=
3856 double kramers_integ_k_integrandB(double x, void *params)
3858 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
3859 double U = (*p->E)(x, p->E_params);
3860 double Fx = p->F * x;
3861 double kbT = k_B * p->env->T;
3862 double intA = kramers_integ_k_integralA(x, params);
3863 //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA;
3864 return exp((U-Fx)/kbT)*intA;
3866 double kramers_integ_k_integralB(double F, environment_t *env, void *params)
3868 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
3870 double result, abserr;
3871 size_t neval; /* for qng */
3872 static gsl_integration_workspace *w = NULL;
3873 if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */
3874 f.function = &kramers_integ_k_integrandB;
3878 //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0);
3879 assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0);
3880 //fprintf(stderr, "integralB = %g\n", result);
3885 With the integrals taken care of, evaluating $k$ is simply
3887 k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})}
3889 <<kramers integ k function declaration>>=
3890 double kramers_integ_k(double F, environment_t *env, void *kramers_params);
3892 <<kramers integ k function>>=
3893 double kramers_integ_k(double F, environment_t *env, void *kramers_params)
3894 { /* Returns the rate constant k in frac pop per s. Takes F in N. */
3895 kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params;
3896 return p->D/kramers_integ_k_integralB(F, env, kramers_params);
3900 <<kramers integ k structure create/destroy function declarations>>=
3901 void *string_create_kramers_integ_param_t(char **param_strings);
3902 void destroy_kramers_integ_param_t(void *p);
3905 <<kramers integ k structure create/destroy functions>>=
3906 kramers_integ_param_t *create_kramers_integ_param_t(double D,
3907 double xmin, double xmax, func_t *E,
3909 destroy_data_func_t *destroy_E_params)
3911 kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t));
3912 assert(ret != NULL);
3917 ret->E_params = E_params;
3918 ret->destroy_E_params = destroy_E_params;
3922 void *string_create_kramers_integ_param_t(char **param_strings)
3924 //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n",
3925 // param_strings[0],param_strings[1],param_strings[2],param_strings[3]);
3926 void *E_params = string_create_spline_param_t(param_strings+1);
3927 return create_kramers_integ_param_t(atof(param_strings[0]),
3928 atof(param_strings[2]),
3929 atof(param_strings[3]),
3930 &spline_eval, E_params,
3931 destroy_spline_param_t);
3934 void destroy_kramers_integ_param_t(void *params)
3936 kramers_integ_param_t *p = (kramers_integ_param_t *)params;
3939 (*p->destroy_E_params)(p->E_params);
3945 Finally we have the GSL spline wrappers:
3946 <<spline functions>>=
3948 <<create/destroy spline>>
3951 <<spline function>>=
3952 double spline_eval(double x, void *spline_params)
3954 spline_param_t *p = (spline_param_t *)(spline_params);
3955 return gsl_spline_eval(p->spline, x, p->acc);
3959 <<create/destroy spline>>=
3960 spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
3962 spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t));
3963 assert(ret != NULL);
3964 ret->n_knots = n_knots;
3967 ret->acc = gsl_interp_accel_alloc();
3968 ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots);
3969 gsl_spline_init(ret->spline, x, y, n_knots);
3973 /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */
3974 void *string_create_spline_param_t(char **param_strings)
3978 double *x=NULL, *y=NULL;
3979 /* break into ordered pairs */
3980 parse_list_string(param_strings[0], SEP, '(', ')',
3981 &num_knots, &knot_coords);
3982 x = (double *)malloc(sizeof(double)*num_knots);
3984 y = (double *)malloc(sizeof(double)*num_knots);
3986 for (i=0; i<num_knots; i++) {
3987 if (sscanf(knot_coords[i], "%lg,%lg", x+i, y+i) != 2) {
3988 fprintf(stderr, "Could not parse pair %d, \"%s\"\n", i, knot_coords[i]);
3991 //fprintf(stderr, "%d : scanned '%s' -> (%g, %g)\n", i, knot_coords[i], x[i], y[i]);
3993 free_parsed_list(num_knots, knot_coords);
3994 return create_spline_param_t(num_knots, x, y);
3997 void destroy_spline_param_t(void *params)
3999 spline_param_t *p = (spline_param_t *)params;
4002 gsl_spline_free(p->spline);
4004 gsl_interp_accel_free(p->acc);
4014 <<kramers integ k global declarations>>=
4015 extern int num_kramers_integ_params;
4016 extern char *kramers_integ_param_descriptions[];
4017 extern char kramers_integ_param_string[];
4020 <<kramers integ k globals>>=
4021 int num_kramers_integ_params = 4;
4022 char *kramers_integ_param_descriptions[] = {
4023 "diffusion D, m^2 / s",
4024 "spline knots, {(x1,y1),(x2,y2),...}, (m,J)",
4025 "starting integration bound x_min, m",
4026 "ending integration bound x_max, m"};
4027 //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";
4028 //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";
4029 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";
4030 /* D from Schlierf 2006, Biophys 90, 4, L33, sup.
4031 * Anchor points from private communication with Schlierf, Apr 9th, 2008.
4032 * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */
4035 <<kramers integ k model getopt>>=
4036 {"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}
4038 Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
4039 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
4041 \subsection{Unfolding model implementation}
4045 <<k model includes>>
4046 #include "k_model.h"
4047 <<k model internal definitions>>
4049 <<k model functions>>
4052 <<k model includes>>=
4053 #include <assert.h> /* assert() */
4054 #include <stdlib.h> /* NULL, malloc() */
4055 #include <math.h> /* HUGE_VAL, sqrt(), exp() */
4056 #include <stdio.h> /* fprintf(), stdout */
4057 #include <string.h> /* strlen(), strcpy() */
4058 #include <gsl/gsl_integration.h>
4059 #include <gsl/gsl_spline.h>
4064 <<k model internal definitions>>=
4065 <<const k structure>>
4066 <<bell k structure>>
4067 <<kramers k structure>>
4068 <<kramers integ k structure>>
4069 <<spline param structure>>
4072 <<k model globals>>=
4076 <<kramers k globals>>
4077 <<kramers integ k globals>>
4080 <<k model functions>>=
4081 <<null k functions>>
4082 <<const k functions>>
4083 <<bell k functions>>
4084 <<kramers k functions>>
4085 <<kramers integ k functions>>
4088 \subsection{Unfolding model unit tests}
4090 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
4091 <<check-k-model.c>>=
4092 <<k model check includes>>
4093 <<check relative error>>
4095 <<k model test suite>>
4096 <<main check program>>
4099 <<k model check includes>>=
4100 #include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
4101 #include <stdio.h> /* sprintf() */
4102 #include <assert.h> /* assert() */
4103 #include <math.h> /* exp() */
4106 #include "k_model.h"
4109 <<k model test suite>>=
4112 <<k model suite definition>>
4115 <<k model suite definition>>=
4116 Suite *test_suite (void)
4118 Suite *s = suite_create ("k model");
4119 <<const k test case defs>>
4120 <<bell k test case defs>>
4122 <<const k test case adds>>
4123 <<bell k test case adds>>
4128 \subsubsection{Constant}
4130 <<const k test case defs>>=
4131 TCase *tc_const_k = tcase_create("Constant k");
4134 <<const k test case adds>>=
4135 tcase_add_test(tc_const_k, test_const_k_create_destroy);
4136 tcase_add_test(tc_const_k, test_const_k_over_range);
4137 suite_add_tcase(s, tc_const_k);
4142 START_TEST(test_const_k_create_destroy)
4144 char *knot[] = {"1","2","3","4","5","6"};
4145 char *params[] = {knot[0]};
4148 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4149 params[0] = knot[i];
4150 p = string_create_const_param_t(params);
4151 destroy_const_param_t(p);
4156 START_TEST(test_const_k_over_range)
4158 double F, Fm=0.0, FM=10, dF=.1, T=300.0;
4159 char *knot[] = {"1","2","3","4","5","6"};
4160 char *params[] = {knot[0]};
4163 char param_string[80];
4166 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4167 params[0] = knot[i];
4168 p = string_create_const_param_t(params);
4169 for ( F=Fm, F<FM, F+=dF ) {
4170 fail_unless(const_k(F, &env, p)==atof(knot[i]),
4171 "const_k(%g, %g, &{%s,%s}) = %g != %s",
4172 F, T, knot[i], dx, const_k(F, &env, p), knot[i]);
4174 destroy_const_param_t(p);
4180 \subsubsection{Bell}
4182 <<bell k test case defs>>=
4183 TCase *tc_bell = tcase_create("Bell's k");
4186 <<bell k test case adds>>=
4187 tcase_add_test(tc_bell, test_bell_k_create_destroy);
4188 tcase_add_test(tc_bell, test_bell_k_at_zero);
4189 tcase_add_test(tc_bell, test_bell_k_at_one);
4190 suite_add_tcase(s, tc_bell);
4194 START_TEST(test_bell_k_create_destroy)
4196 char *knot[] = {"1","2","3","4","5","6"};
4198 char *params[] = {knot[0], dx};
4201 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4202 params[0] = knot[i];
4203 p = string_create_bell_param_t(params);
4204 destroy_bell_param_t(p);
4209 START_TEST(test_bell_k_at_zero)
4211 double F=0.0, T=300.0;
4212 char *knot[] = {"1","2","3","4","5","6"};
4214 char *params[] = {knot[0], dx};
4217 char param_string[80];
4220 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4221 params[0] = knot[i];
4222 p = string_create_bell_param_t(params);
4223 fail_unless(bell_k(F, &env, p)==atof(knot[i]),
4224 "bell_k(%g, %g, &{%s,%s}) = %g != %s",
4225 F, T, knot[i], dx, bell_k(F, &env, p), knot[i]);
4226 destroy_bell_param_t(p);
4231 START_TEST(test_bell_k_at_one)
4234 char *knot[] = {"1","2","3","4","5","6"};
4236 char *params[] = {knot[0], dx};
4237 double F= k_B*T/atof(dx);
4240 char param_string[80];
4243 for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
4244 params[0] = knot[i];
4245 p = string_create_bell_param_t(params);
4246 CHECK_ERR(1e-6, atof(knot[i])*exp(1.0), bell_k(F, &env, p));
4247 //fail_unless(bell_k(F, &env, p)==atof(knot[i])*exp(1.0),
4248 // "bell_k(%g, %g, &{%s,%s}) = %g != %g",
4249 // F, T, knot[i], dx, bell_k(F, &env, p), atof(knot[i])*exp(1.0));
4250 destroy_bell_param_t(p);
4256 <<kramers k tests>>=
4259 <<kramers k test case def>>=
4262 <<kramers k test case add>>=
4265 <<k model function tests>>=
4266 <<check relative error>>
4268 START_TEST(test_something)
4270 double k=5, x=3, last_x=2.0, F;
4271 one_dim_func_t *handlers[] = {&hooke};
4272 void *data[] = {&k};
4273 double xi[] = {0.0};
4275 int new_active_groups = 1;
4276 F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x);
4277 fail_unless(F = k*x, NULL);
4282 \subsection{Utilities}
4284 The unfolding models can get complicated, and are sometimes hard to explain to others.
4285 The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without
4286 the overhead of having to understand a full simulation.
4287 [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen,
4288 or special mode, where the operation depends on the specific model selected.
4290 <<k-model-utils.c>>=
4292 <<k model utility includes>>
4293 <<k model utility definitions>>
4294 <<k model utility getopt functions>>
4295 <<k model utility multi model E>>
4296 <<k model utility main>>
4299 <<k model utility main>>=
4300 int main(int argc, char **argv)
4302 <<k model getopt array definition>>
4303 k_model_getopt_t *model = NULL;
4308 int num_param_args; /* for INIT_MODEL() */
4309 char **param_args; /* for INIT_MODEL() */
4310 double special_xmin;
4311 double special_xmax;
4312 get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
4313 &special_xmin, &special_xmax, &flags);
4314 if (flags & VFLAG) {
4315 printf("#initializing model %s with parameters %s\n", model->name, model->params);
4317 INIT_MODEL("utility", model, params);
4320 if (model->k == NULL) {
4321 printf("No k function for model %s\n", model->name);
4325 double dF=1e-11, F=0, k=1.0;
4326 printf("#F (N)\tk (%% pop. per s)\n");
4327 while (k > 0 && k < 1e5) {
4328 k = (*model->k)(F, &env, params);
4329 printf("%g\t%g\n", F, k);
4335 if (model->k == NULL || model->k == &bell_k) {
4336 printf("No special mode for model %s\n", model->name);
4340 double dx=(special_xmax-special_xmin)/500, x=special_xmin, E;
4341 printf("#x (m)\tE (J)\n");
4342 while (x <= special_xmax) {
4343 E = multi_model_E(model->k, params, &env, x);
4344 printf("%g\t%g\n", x, E);
4352 if (model->destructor)
4353 (*model->destructor)(params);
4358 <<k model utility multi model E>>=
4359 double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
4361 kramers_integ_param_t *p = (kramers_integ_param_t *)model_params;
4363 if (k == kramers_integ_k)
4364 E = (*p->E)(x, p->E_params);
4366 E = kramers_E(0, env, model_params, x);
4372 <<k model utility includes>>=
4375 #include <string.h> /* strlen() */
4376 #include <assert.h> /* assert() */
4379 #include "string_eval.h"
4380 #include "k_model.h"
4383 <<k model utility definitions>>=
4384 <<version definition>>
4385 #define VFLAG 1 /* verbose */
4386 enum mode_t {M_K_OF_F, M_SPECIAL};
4387 <<string comparison definition and globals>>
4388 <<k model getopt definitions>>
4389 <<initialize model definition>>
4390 <<kramers k structure>>
4391 <<kramers integ k structure>>
4395 <<k model utility getopt functions>>=
4397 <<k model utility help>>
4398 <<k model utility get options>>
4401 <<k model utility help>>=
4402 void help(char *prog_name,
4404 int n_k_models, k_model_getopt_t *k_models,
4408 printf("usage: %s [options]\n", prog_name);
4409 printf("Version %s\n\n", VERSION);
4410 printf("A utility for understanding the available unfolding models\n\n");
4411 printf("Environmental options:\n");
4412 printf("-T T\tTemperature T (currently %g K)\n", env->T);
4413 printf("-C T\tYou can also set the temperature T in Celsius\n");
4414 printf("Model options:\n");
4415 printf("-k model\tTransition rate model (currently %s)\n",
4416 k_models[k_model].name);
4417 printf("-K args\tTransition rate model argument string (currently %s)\n",
4418 k_models[k_model].params);
4419 printf("There are two output modes. In standard mode, k(F) is printed\n");
4420 printf("For example:\n");
4421 printf(" #Force (N)\tk (% pop. per s)\n");
4422 printf(" 123.456\t7.89\n");
4424 printf("In special mode, the output depends on the model.\n");
4425 printf("For models defining an energy function E(x), that function is printed\n");
4426 printf(" #Position (m)\tE (J)\n");
4427 printf(" 0.0012\t0.0034\n");
4429 printf("-m\tChange output to standard mode\n");
4430 printf("-M\tChange output to special mode\n");
4431 printf("-x\tSet the minimum x value for the special mode E(x) output\n");
4432 printf("-X\tSet the maximum x value for the special mode E(x) output\n");
4433 printf("-V\tChange output to verbose mode\n");
4434 printf("-h\tPrint this help and exit\n");
4436 printf("Unfolding rate models:\n");
4437 for (i=0; i<n_k_models; i++) {
4438 printf("\t%s\t%s\n", k_models[i].name, k_models[i].description);
4439 for (j=0; j < k_models[i].num_params ; j++)
4440 printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
4441 printf("\t\tdefaults: %s\n", k_models[i].params);
4446 <<k model utility get options>>=
4447 void get_options(int argc, char **argv, environment_t *env,
4448 int n_k_models, k_model_getopt_t *k_models,
4449 enum mode_t *mode, k_model_getopt_t **model,
4450 double *special_xmin, double *special_xmax, unsigned int *flags)
4452 char *prog_name = NULL;
4453 char c, options[] = "T:C:k:K:mMx:X:Vh";
4455 extern char *optarg;
4456 extern int optind, optopt, opterr;
4460 /* setup defaults */
4462 prog_name = argv[0];
4463 env->T = 300.0; /* K */
4467 *special_xmax = 1e-8;
4469 while ((c=getopt(argc, argv, options)) != -1) {
4471 case 'T': env->T = atof(optarg); break;
4472 case 'C': env->T = atof(optarg)+273.15; break;
4474 k_model = index_k_model(n_k_models, k_models, optarg);
4475 *model = k_models+k_model;
4478 k_models[k_model].params = optarg;
4480 case 'm': *mode = M_K_OF_F; break;
4481 case 'M': *mode = M_SPECIAL; break;
4482 case 'x': *special_xmin = atof(optarg); break;
4483 case 'X': *special_xmax = atof(optarg); break;
4484 case 'V': *flags |= VFLAG; break;
4486 fprintf(stderr, "unrecognized option '%c'\n", optopt);
4487 /* fall through to default case */
4489 help(prog_name, env, n_k_models, k_models, k_model);
4498 \section{\LaTeX\ documentation}
4500 The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations.
4501 The comment blocks are extracted (with nicely formatted code blocks), using
4502 <<latex makefile lines>>=
4503 sawsim.tex : sawsim.nw
4504 noweave -latex -index -delay $^ > $@
4507 <<latex makefile lines>>=
4508 sawsim.pdf : sawsim.tex
4514 The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information,
4515 the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and
4516 the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings.
4518 [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
4519 <<latex makefile lines>>=
4521 rm -f sawsim.aux sawsim.bbl sawsim.blg sawsim.log sawsim.out sawsim.tex
4522 clean_latex : semi-clean_latex
4528 [[make]] is a common utility on *nix systems for managing dependencies while building software.
4529 In this case, we have several \emph{targets} that we'd like to build:
4530 the main simulation program \prog;
4531 the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]];
4532 the documentation [[sawsim.pdf]];
4533 and the [[Makefile]] itself.
4534 Besides the generated files, there is the utility target
4535 [[clean]] for removing all generated files except [[Makefile]].
4537 % [[dist]] for generating a distributable tar file.
4539 Extract the makefile with
4540 `[[notangle -Rmakefile sawsim.nw | sed 's/ /\t/' > Makefile]]'.
4541 The sed is needed because notangle replaces tabs with eight spaces,
4542 but [[make]] needs tabs.
4544 all : sawsim sawsim.pdf Makefile tension_model_utils k_model_utils
4548 profile : sawsim_profile
4549 sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ -mnull -Mwlc -A0.39e-9,28e-9 -F8
4550 gprof sawsim_profile gmon.out > $@
4552 clean : clean_check clean_list clean_tension clean_k_model clean_parse clean_string_eval clean_latex
4553 rm -f sawsim.c sawsim sawsim_static sawsim_profile gmon.out global.h
4555 sawsim.c : sawsim.nw
4557 global.h : sawsim.nw
4558 notangle -Rglobal.h $^ > $@
4559 sawsim : sawsim.c global.h list.c list.h tension_balance.c tension_balance.h \
4560 tension_model.c tension_model.h k_model.c k_model.h parse.c parse.h string_eval.c string_eval.h \
4562 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
4563 sawsim_static : sawsim
4564 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
4565 sawsim_profile : sawsim
4566 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
4569 check_sawsim.c : sawsim.nw
4570 notangle -Rchecks $^ > $@
4571 check_sawsim : check_sawsim.c global.h list.c list.h tension_balance.c tension_balance.h \
4572 k_model.c k_model.h parse.c parse.h string_eval.c string_eval.h accel_k.c accel_k.h
4573 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
4575 rm -f check_sawsim check_sawsim.c
4576 check : check_list check_tension_balance check_k_model check_parse check_string_eval check_accel_k check_sawsim
4578 check_tension_balance
4585 <<list module makefile lines>>
4586 <<tension balance module makefile lines>>
4587 <<tension model module makefile lines>>
4588 <<k model module makefile lines>>
4589 <<parse module makefile lines>>
4590 <<string eval module makefile lines>>
4591 <<accel k module makefile lines>>
4592 <<latex makefile lines>>
4594 Makefile : sawsim.nw
4595 notangle -Rmakefile $^ | sed 's/ /\t/' > $@
4597 This is fairly self-explanatory.
4598 We compile a dynamically linked version ([[sawsim]]) and a statically linked version ([[sawsim_static]]).
4599 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.
4604 \subsection{Hookean springs in series}
4605 \label{app.math_hooke}
4608 The effective spring constant for $n$ springs $k_i$ in series is given by
4610 k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
4616 F &= k_i x_i &&\forall i \le n \\
4617 x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\
4618 x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\
4619 F &= k_1 x_1 = k_\text{eff} x \\
4620 k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}}
4621 = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}.
4626 \addcontentsline{toc}{section}{References}
4627 \bibliography{sawsim}