% sawsim - simulating single molecule mechanical unfolding. % Copyright (C) 2008-2010 William Trevor King % % This program is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation, either version 3 of the License, or % (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program. If not, see . % % The author may be contacted at on the Internet, or % write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St., % Philadelphia PA 19104, USA. %%%%%%%%%%%%%%%%%%%%%%%%% Start LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%% % we have our own preamble, % use `noweave -delay` to not print noweb's automatic one % -*- mode: Noweb; noweb-code-mode: c-mode -*- \documentclass[letterpaper, 11pt]{article} \usepackage{noweb} \pagestyle{noweb} \noweboptions{smallcode} \usepackage{url} % needed for natbib for underscores in urls? \usepackage[breaklinks,colorlinks]{hyperref} % for internal links, \phantomsection % breaklinks breaks long links % colorlinks colors link text instead of boxing it \usepackage{amsmath} % for \text{}, \xrightarrow[sub]{super} \usepackage{amsthm} % for theorem and proof environments \newtheorem{theorem}{Theorem} \usepackage{subfig} % Support side-by-side figures \usepackage{pgf} % Fancy graphics \usepackage{tikz} % A nice, inline PGF frontend \usetikzlibrary{automata} % Graph-theory library \usepackage{doc} % BibTeX \usepackage[super,sort&compress]{natbib} % fancy citation extensions % super selects citations in superscript mode % sort&compress automatically sorts and compresses compound citations (\cite{a,b,...}) \usepackage[mathletters]{ucs} % use math fonts to allow Greek UTF-8 in the text \usepackage[utf8x]{inputenc} % allow UTF-8 characters %\bibliographystyle{ieeetr} % pick the bibliography style, short and sweet %\bibliographystyle{plain} % pick the bibliography style, includes dates \bibliographystyle{plainnat} \defcitealias{sw:noweb}{{\tt noweb}} \defcitealias{galassi05}{GSL} \defcitealias{sw:check}{{\tt check}} \defcitealias{sw:python}{Python} \topmargin -1.0in \headheight 0.5in \headsep 0.5in \textheight 9.5in % leave a bit of extra room for page numbers (at top of page) \oddsidemargin -0.5in \textwidth 7.5in \setlength{\parindent}{0pt} \setlength{\parskip}{5pt} % For some reason, the \maketitle wants to go on it's own page % so we'll just hardcode the following in our first page. %\title{Sawsim: a sawtooth protein unfolding simulator} %\author{W.~Trevor King} %\date{\today} \newcommand{\prog}{[[sawsim]]} \newcommand{\lang}{[[C]]} \newcommand{\p}[3]{\left#1 #2 \right#3} % e.g. \p({complicated expr}) \newcommand{\Th}{\ensuremath{\sup{\text{th}}}} % ordinal th \newcommand{\U}[1]{\ensuremath{\text{ #1}}} % units: $1\U{m/s}$ % single spaced lists, from Alvin Alexander % http://www.devdaily.com/blog/post/latex/control-line-spacing-in-itemize-enumerate-tags/ \newenvironment{packed_item}{ \begin{itemize} \setlength{\itemsep}{1pt} \setlength{\parskip}{0pt} \setlength{\parsep}{0pt} }{\end{itemize}} \begin{document} \nwfilename{sawsim.nw} \nwbegindocs{0} %\maketitle \begin{centering} Sawsim: a sawtooth protein unfolding simulator \\ W.~Trevor King \\ \today \\ \end{centering} \bigskip \tableofcontents %%%%%%%%%%%%%%%%%%%%%%%%%%%% End LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Introduction} The unfolding of multi-globular protein strings is a stochastic process. In the \prog\ program, we use Monte Carlo methods to simulate this unfolding through an explicit time-stepping approach. We develop a program in \lang\ to simulate probability of unfolding under a constant extension velocity of a chain of globular domains. In order to extract physical meaning from many mechanical unfolding simulations, it is necessary to compare the experimental results against the results of simulations using different models and parameters. The model and parameters that best match the experimental data provide the ‘best’ model for that experiment. Previous mechanical unfolding studies have rolled their own simulations for modelling their experiment, and there has been little exchange and standardization of these simulations between groups. The \prog\ program attempts to fill this gap, providing a flexible, extensible, \emph{community} program, to make it easier for groups to share the burden and increase the transparency of data analysis. ‘Sawsim’ is blend of ‘sawtooth simulation’. \subsection{Related work} Sawim has been published\citep{king10}! It is, as far as I know, the only open-source Monte Carlo force spectroscopy simulator, but similar closed source simulators have been around since the seminal paper by \citet{rief97a}. In chronological order, major contributions have come from \begin{packed_item} \item \citet{rief97a} \citeyear{rief97a}: \item \citet{rief97b} \citeyear{rief97b}: \item \citet{kellermayer97} \citeyear{kellermayer97}: \item \citet{rief98} \citeyear{rief98}: first paper to focus mainly on the simulation \item \citet{oberhauser98} \citeyear{oberhauser98}: \item \citet{carrion-vazquez99a} \citeyear{carrion-vazquez99a}: \item \citet{best02} \citeyear{best02}: pointed out large fit valley \item \citet{zinober02} \citeyear{zinober02}: introduced the scaffold effect \item \citet{jollymore09} \citeyear{jollymore09}: \item \citet{king10} \citeyear{king10}: introduced sawsim \end{packed_item} \subsection{About this document} This document was written in \citetalias{sw:noweb} with code blocks in \lang\ and comment blocks in \LaTeX. \subsection{System Requirements and Dependencies} \subsection{Change Log} Version 0.0 used only the unfolded domain WLC to determine the tension and had a fixed timestep $dt$. Version 0.1 added adaptive timesteps, adjusting so that the unfolding probability for a given domain was constant. Version 0.2 added Kramers' model unfolding rate calculations, and a switch to select Bell's or Kramers' model. This induced a major rewrite, introducing the tension group abstraction, and a split to multiple \lang\ source files. Also added a unit-test suites for sanity checking, and switched to SI units throughout. Version 0.3 added integrated Kramers' model unfolding rate calculations. Also replaced some of my hand-coded routines with GNU Scientific Library \citepalias{galassi05} calls. Version 0.4 added Python evaluation for the saddle-point Kramers' model unfolding rate calculations, so that the model functions could be controlled from the command line. Also added interpolating lookup tables to accelerate slow unfolding rate model evaluation, and command line control of tension groups. Version 0.5 added the stiffness environmental parameter and it's associated unfolding models. Version 0.6 generalized from two state unfolding to arbitrary $N$-state rate reactions. This allows simulations of unfolding-refolding, metastable intermediates, etc., but required reorganizing the command line interface. Version 0.7 added tension model inverses, which seems to reduce computation time by about a factor of 2. I was expecting more, but I'll take what I can get. Version 0.8 fixed a minor bug in unfolding probability for multi-domain groups. The probability of at least one domain unfolding had been calculated as $P_N=NP_1$, rather than $P_N=1-(1-P_1)^N$. However, for small $P$ the two are equivalent. Version 0.9 added the [[-P]] option to sawsim, setting the target probability for the per-state transition $P_N$, not the per-domain transisiton $P_1$. Version 0.10 added the [[-F]] option to sawsim, setting a limit on the force change $dF$ in a single timestep for continuous pulling simulations. It also added the piston tension model (Section \ref{sec.piston_tension}), and adjusted the stiffness calculations to handle domains with undefined stiffness. <>= #define VERSION "0.10" @ \subsection{License} <>= sawsim - simulating single molecule mechanical unfolding. Copyright (C) 2008-2010, William Trevor King This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . The author may be contacted at on the Internet, or write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St., Philadelphia PA 19104, USA. @ <>= /* <> */ @ \section{Main} \label{sec.main} The unfolding system is stored as a chain of \emph{domains} (Figure \ref{fig.domain_chain}). Domains can be folded, globular protein domains, unfolded protein linkers, AFM cantilevers, or other stretchable system components. The system is modeled as graph of possible domain states with transitions between them (Figure \ref{fig.domain_states}). \begin{figure} \begin{center} \subfloat[][]{\label{fig.domain_chain} \begin{tikzpicture}[->,node distance=2cm,font=\footnotesize] \tikzstyle{every state}=[fill,draw=red!50,very thick,fill=red!20] \node[state] (A) {domain 1}; \node[state] (B) [below of=A] {domain 2}; \node[state] (C) [below of=B] {.~.~.}; \node[state] (D) [below of=C] {domain $N$}; \node (S) [below of=D] {Surface}; \node (E) [above of=A] {}; \path[-] (A) edge (B) (B) edge node [right] {Tension} (C) (C) edge (D) (D) edge (S); \path[->,green] (A) edge node [right,black] {Extension} (E); \end{tikzpicture} } \qquad \subfloat[][]{\label{fig.domain_states} \begin{tikzpicture}[->,node distance=2.5cm,shorten <=1pt,shorten >=1pt,font=\footnotesize] \tikzstyle{every state}=[fill,draw=blue!50,very thick,fill=blue!20] \node[state] (A) {cantilever}; \node[state] (C) [below of=A] {transition}; \node[state] (B) [left of=C] {folded}; \node[state] (D) [right of=C] {unfolded}; \path (B) edge [bend left] node [above] {$k_1$} (C) (C) edge [bend left] node [below] {$k_1'$} (B) edge [bend left] node [above] {$k_2$} (D) (D) edge [bend left] node [below] {$k_2'$} (C); \end{tikzpicture} } \caption{\subref{fig.domain_chain} Extending a chain of domains. One end of the chain is fixed, while the other is extended at a constant speed. The domains are coupled with rigid linkers, so the domains themselves must stretch to accomodate the extension. \subref{fig.domain_states} Each domain exists in a discrete state. At each timestep, it may transition into another state following a user-defined state matrix such as this one, showing a metastable transition state and an explicit ``cantilever'' domain.} \end{center} \end{figure} Each domain contributes to the total tension, which depends on the chain extension. The particular model for the tension as a function of extension is handled generally with function pointers. So far the following models have been implemented: \begin{packed_item} \item Null (Appendix \ref{sec.null_tension}), \item Constant (Appendix \ref{sec.const_tension}), \item Hookean spring (Appendix \ref{sec.hooke}), \item Worm-like chain (WLC, Appendix \ref{sec.wlc}), and \item Freely-jointed chain (FJC, Appendix \ref{sec.fjc}). \item Piston (Appendix \ref{sec.piston_tension}), \end{packed_item} The tension model and parameters used for a particular domain depend on the domain's current state. The transition rates between states are also handled generally with function pointers, with implementations for \begin{packed_item} \item Null (Appendix \ref{sec.null_k}), \item Bell's model (Appendix \ref{sec.bell}), \item Spring constant corrected Bell model (Appendix \ref{sec.kbell}), \item Kramers' integral model (Appendix \ref{sec.kramers_integ}), and \item Kramers' saddle point model (Appendix \ref{sec.kramers}). \end{packed_item} The unfolding of the chain domains is modeled in two stages. First the tension of the chain is calculated. Then the tension (assumed to be constant over the length of the chain, see Section \ref{sec.timescales}), is applied to each domain, and used to calculate the probability of that domain changing state in the next timestep $dt$. Then the time is stepped forward, and the process is repeated. The simulation is complete when a pre-set time limit [[t_max]] is reached or a pre-selected state [[stop_state]] is empty. <
>= int main(int argc, char **argv) { <> <> /* variables initialized in get_options() */ list_t *state_list=NULL, *transition_list=NULL; environment_t env={0}; double P, dF, t_max, dt_max, v, x_step; state_t *stop_state=NULL; /* variables used in the simulation loop */ double t=0, x=0, dt, F; /* time, distance, timestep, force */ int transition=1; /* boolean marking a transition for tension and timestep optimization */ get_options(argc, argv, &P, &dF, &t_max, &dt_max, &v, &x_step, &stop_state, &env, NUM_TENSION_MODELS, tension_models, NUM_K_MODELS, k_models, &state_list, &transition_list); setup(); if (flags & FLAG_OUTPUT_FULL_CURVE) printf("#Distance (m)\tForce (N)\tStiffness (N/m)\n"); if (flags & FLAG_OUTPUT_TRANSITION_FORCES) printf("#Force (N)\tInitial state\tFinal state\n"); while ((t_max > 0 && t <= t_max) || (stop_state != NULL && stop_state->num_domains > 0)) { F = find_tension(state_list, &env, x, transition, 0); if (x_step == 0) dt = determine_dt(state_list, transition_list, &env, F, x, dF, P, dt_max, v, transition); else dt = determine_dt(state_list, transition_list, &env, F, x, dF, P, dt_max, 0, transition); transition=random_transitions(transition_list, F, dt, &env); if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\t%g\n", x, F, env.stiffness); t += dt; if (x_step == 0) { x += v*dt; } else { if (dt == dt_max) { /* step completed */ x += x_step; dt_max = x_step / v; } else { /* still working on this step */ dt_max -= dt; } } } destroy_state_list(state_list); destroy_transition_list(transition_list); free_accels(); return 0; } @ The meat of the simulation is bundled into the three functions [[find_tension]], [[determine_dt]], and [[random_transitions]]. [[find_tension]] is discussed in Section \ref{sec.find_tension}, [[determine_dt]] in Section \ref{sec.adaptive_dt}, and [[random_transitions]] in Section \ref{sec.transition_rate}. The stretched distance is extended in one of two modes depending on [[xstep]]. If [[xstep == 0]], then the pulling is continuous, at least within the limits of the inherent discretization of a time stepped approach. At any rate, the timestep is limited by the maximum allowed timestep [[dt_max]] and the maximum allowed unfolding probability in a given step [[P]]. In the second mode [[xstep > 0]], and the end to end distance increases in discrete steps of that length. The time between steps is chosen to maintain an average unfolding velocity [[v]]. This approach is less work to simulate than smooth pulling and also closer to the reality of many experiments, but it is practically impossible to treat analytically. With both pulling styles implemented, the effects of the discretization can be easily calculated, bridging the gap between theory and experiment. Environmental parameters important in determining reaction rates and tensions (e.g.\ temperature, pH) are stored in a single structure to facilitate extension to more complicated models in the future. <>= typedef struct environment_struct { double T; double stiffness; } environment_t; @ <>= #ifndef GLOBAL_H #define GLOBAL_H #define DOUBLE_PRECISION 1e-12 <> <> <> <> #endif /* GLOBAL_H */ @ Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being included multiple times. \section{Simulation functions} <>= <> <> <> <> <> @ \subsection{Tension} \label{sec.find_tension} Because the stretched system may be made up of several parts (folded domains, unfolded domains, spring-like cantilever, \ldots), we will assign the domains to states with tension models and groups. The models are used to determine the tension of all the domain in a given state following some model (Hookean spring, WLC, \ldots). The domains are assumed to be commutative, so ordering is ignored. The interactions between the groups are assumed to be linear, but the interactions between domains of the same group need not be. Each model has a tension handler function, which gives the tension $F$ needed to stretch a given group of domains a distance $x$. Groups represent collections of states which the model should treat as a single entity. To understand the purpose of groups, consider a system with two unfolded domain states (e.g.\ for two different proteins) that were both modeled as WLCs. If both unfolded states had the same persistence length, it would be wise to assign them to the same group. This would reduce both the computational cost of balancing the tension and the error associated with the inter-group interaction linearity. Note that group numbers only matter \emph{within} model classes, since grouping domains with seperate models doesn't make sense. Within designated groups, the tension parameters for each domain are still checked for compatibility, so if you accidentally assign incompatible domains to the same group you should get an appropriate error message. <>= #define NUM_TENSION_MODELS 6 @ <>= extern one_dim_fn_t *tension_handlers[]; @ <>= one_dim_fn_t *tension_handlers[] = { NULL, &const_tension_handler, &hooke_handler, &wlc_handler, &fjc_handler, }; @ <>= <> @ <>= typedef struct tension_handler_data_struct { list_t *group_tension_params; /* one item for each domain in the group */ environment_t *env; void *persist; } tension_handler_data_t; @ After sorting the chain into separate groups $G_i$, with tension handlers $F_i(G_i; x_i)$, 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$. Note that there may be several states within each group. [[find_tension]] is basically a wrapper to feed proper input into the [[tension_balance]] function. [[transition]] should be set to 0 if there were no transitions since the previous call to [[find_tension]] to support some optimizations that assume a small increase in tension between steps (only true if no transition has occured). While we're messing with the tension handlers and if [[const_env==0]], we also grab a measure of the current linker stiffness for the environmental variable, which is needed by some of the unfolding rate models (e.g.\ the linker-stiffness corrected Bell model, Appendix \ref{sec.kbell}). <>= double find_tension(list_t *state_list, environment_t *env, double x, int transition, int const_env) { // TODO: !cont_env -> step_call, only alter env or statics if step_call==1 static int num_active_groups=0; static one_dim_fn_t **per_group_handlers = NULL; static one_dim_fn_t **per_group_inverse_handlers = NULL; static void **per_group_data = NULL; /* array of pointers to tension_handler_data_t */ static double last_x; tension_handler_data_t *tdata; double F, *pStiffness; int i; #ifdef DEBUG fprintf(stderr, "%s:\tfinding tension at distance %g%s\n", __FUNCTION__, x, transition ? " after a transition" : ""); #endif if (transition != 0 || num_active_groups == 0) { /* setup statics */ /* get new statics, freeing the old ones if they aren't NULL */ get_active_groups(state_list, &num_active_groups, &per_group_handlers, &per_group_inverse_handlers, &per_group_data); /* fill in the group handlers and params */ #ifdef DEBUG fprintf(stderr, "%s:\t%d groups, tension handler array %p, inverse tension handler array %p, handler data array %p, (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_inverse_handlers, per_group_data, per_group_data[0]); #endif for (i=0; igroup_tension_params); list_print(stderr, tdata->group_tension_params, " parameters"); #endif tdata->env = env; /* tdata->persist continues unchanged */ } last_x = -1.0; } /* else continue with the current statics */ if (const_env == 1) pStiffness = NULL; else pStiffness = &env->stiffness; F = tension_balance(num_active_groups, per_group_handlers, per_group_inverse_handlers, (void **)per_group_data, last_x, x, pStiffness); last_x = x; return F; } @ 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). We'll also restrict our extensions to all be positive. With these restrictions, the problem of balancing the tensions reduces to solving a system of $N+1$ possibly non-linear equations with $N$ unknowns, where $N$ is the number of active groups. 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. For example, it guarantees the existence of a unique, real solution for finite forces. See Appendix \ref{sec.tension_balance} for the tension balancing implementation. Each group must define a parameter structure if the tension model is parametrized, a function to create the parameter structure, a function to destroy the parameter structure, and a tension function of type [[one_dim_fn_t]] that receives a [[tension_handler_data_t *]] as its data argument. The tension-specific parameter structures are contained in the domain group field. For optimization, a group may also define an inverse tension function as an optimization. See the Section \ref{sec.model_selection} for a discussion on the command line framework and inverse function discussion. See the worm-like chain implementation for an example (Section \ref{sec.wlc}). The major limitation of our approach is that all of our tension models are Markovian (memory-less), which is not really a problem for the chains (see \citet{evans99} for WLC thermalization rate calculations). \subsection{Transition rate} \label{sec.transition_rate} Each state transition is modeled as unimolecular, first order reaction $$ \text{State 1} \xrightarrow{k} \text{State 2} $$ With the rate constant $k$ defined by $$ \frac{d[\text{State 2}]}{dt} = -\frac{d[\text{State 1}]}{dt} = k [\text{State 1}]. $$ So the probability of a given domain transitioning in a short time $dt$ is given by $$ P_1 = k \cdot dt. $$ For a group of $N$ identical domains, and therefore identical unfolding probabilities $P_1$, the probability of $n$ domains unfolding in a given timestep is given by the binomial distribution $$ P(n) = \frac{N!}{n!(N-n)!}P_1^n(1-P_1)^{(N-n)}. $$ The probability that \emph{none} of these domains unfold is then $$ P(0) = (1-P_1)^N, $$ so the probability that \emph{at least one} domain unfolds is $$ P_N = 1-(1-P_1)^N. $$ Note that for small $P$, $$ p(n>0) = 1-(1-NP_1+\mathcal{O}(P_1^2)) = NP_1 - \mathcal{O}(P^2) \approx NP_1. $$ This conversion from $P_1$ to $P_N$ is handled by the [[pN]] macro <>= /* find multi-domain transition rate from the single domain rate */ #define PN(P,N) (1.0-pow(1.0-(P), (N))) @ We can also define a macro to find the approprate $dt$ to achieve a target $P_N$ with a given $k$ and $N$. \begin{align} P_N &= 1-(1-P_1)^N \\ 1-P_1 &= (1-P_N)^{1/N} \\ P_1 &= 1-(1-P_N)^{1/N} \end{align} <>= #define P1(PN,N) (1.0-pow(1.0-(PN), 1.0/(N))) @ We take some time to discuss the meaning of $p(n>1)$ (i.e. multi-unfolding timesteps) in Section \ref{sec.timescales}. <>= int domain_transitions(double F, double dt, environment_t *env, transition_t *transition) { /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to transition structure */ double k; k = accel_k(transition->k, F, env, transition->k_params); //(*transition->k)(F, env, domain->k_params); //printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt); return happens(PN(k*dt, N_INIT(transition))); } @ [[happens]] is a random decision making function defined in Appendix \ref{sec.utils}. <>= int random_transitions(list_t *transition_list, double tension, double dt, environment_t *env) { /* returns 1 if there was a transition and 0 otherwise */ while (transition_list != NULL) { if (T(transition_list)->initial_state->num_domains > 0 && domain_transitions(tension, dt, env, T(transition_list))) { if (flags & FLAG_OUTPUT_TRANSITION_FORCES) fprintf(stdout, "%g\t%s\t%s\n", tension, T(transition_list)->initial_state->name, T(transition_list)->final_state->name); T(transition_list)->initial_state->num_domains--; T(transition_list)->final_state->num_domains++; return 1; /* our one domain has transitioned, stop looking for others */ } transition_list = transition_list->next; } return 0; } @ \subsection{Timescales} \label{sec.timescales} The simulation assumes that chain equilibration is instantaneous relative to the simulation timestep $dt$, so that tension is uniform along the chain. The quality of this assumption depends on your particular chain. For example, a damped spring thermalizes on a timescale of order $\tau = 1/\gamma$ where the damping ratio $\gamma \equiv \omega_0(-\zeta \pm \sqrt{\zeta^2-1})$, the natural angular frequency $\omega_0 \equiv \sqrt{k/m}$, $\zeta \equiv b/2\sqrt{km}$, and $b$ sets the drag $F_b = -bv$. For our cantilevers $\tau$ is on the order of milliseconds, which is longer than a timestep. However, the approximation is still reasonable, because there is rarely unfolding during the cantilever return. You could, of course, take cantilever, WLC, etc.\ response times into effect, but that would significantly complicate a simulation that seems to work well enough without bothering :p. For WLC and FJC relaxation times, see the Appendix of \citet{evans99} and \citet{kroy07}. Besides assuming our timestep is much greater than equilibration times, we also force it to be short enough so that only one domain may unfold in a given timestep. For an unfolding timescale of $dt_u = 1/k$, we adapt our timesteps to keep $dt \ll dt_u$, so the probability of two domains unfolding in the same timestep is small (see [[determine_dt]] in Section \ref{sec.adaptive_dt} and [[P]] in [[main()]] in Section \ref{sec.main} set by [[-P]] in [[get_options()]] in Section \ref{sec.get_options}). This approach breaks down as the adaptive timestep scheme approaches the transition timestep $dt_t$, but $dt_t \sim 1\U{$\mu$s}$ for Ig27-like proteins \citep{klimov00}, so this shouldn't be a problem. To reassure yourself, you can ask the simulator to print the smallest timestep that was used with TODO. Even if the likelyhood of two domains unfolding in the same timestep is small, if you run long enough simulations it will eventually occur. In this case, we assume that $dt_t \ll dt$, so even if two domains would unfold if we stuck to our unfolding rate $k$ for an entire timestep $dt$, in ``reality'' only one of those domains unfolds. Because we do not know when in the timestep the unfolding took place, we move on to the next timestep without checking for further unfoldings. This ``unchecked timestep portion'' should not contribute significant errors to the calculation, because if $dt$ was small enough that unfolding was unlikely at the pre-unfolding force, it should be even more unlikely at the post-unfolding force, especially over only a fraction of the timestep. The only dangerous cases for the current implementation are those where unfolding intermediates are much less stable than their precursor states, in which case an unfolding event during the remaining timestep may actually be likely. A simple workaround for such cases is to pick the value for [[P]] to be small enough that the $dt \ll dt_u$ assumption survives even under a transition populating the unstable intermediate. \subsection{Adaptive timesteps} \label{sec.adaptive_dt} We'd like to pick $dt$ so the probability of changing state (e.g.\ unfolding) in the next timestep is small. If we don't adapt our timestep, we also risk breaking our approximation that $F(x' \in [x, x+v \cdot dt]) = F(x)$ or that only one domain changes state in a given timestep. The simulation would have been accurate for sufficiently small fixed timesteps, but adaptive timesteps will allow us to move through low tension regions in fewer steps, leading to a more efficient simulation. Consider the two state folded $\rightarrow$ unfolded transition. Because $F(x)$ is monotonically increasing between unfoldings, excessively large timesteps will lead to erroneously small unfolding rates (an thus, higher average unfolding force). The actual adaptive timestep implementation is not particularly interesting, since we are only required to reduce $dt$ to somewhere below a set threshold, so I've removed it to Appendix \ref{sec.adaptive_dt_implementation}. \section{Models} TODO: model intro. The models provide the physics of the simulation, but the simulation allows interchangeable models, and we are currently focusing on the simulation itself, so we remove the actual model implementation to the appendices. The tension models are defined in Appendix \ref{sec.tension_models}, and unfolding models are defined in Appendix \ref{sec.k_models}. <>= #define k_B 1.3806503e-23 /* J/K */ @ \section{Command line interface} <>= <> <> <> <> @ \subsection{Model selection} \label{sec.model_selection} The main difficulty with the command line interface in \prog\ is developing an intuitive method for accessing the possibly large number of available models. We'll treat this generally by defining an array of available models, containing their important parameters. <>= <> @ <>= typedef void *create_data_func_t(char **param_strings); typedef void destroy_data_func_t(void *param_struct); @ <>= <> <> @ In order to choose models, we need a method of assembling a reaction graph such as that in Figure \ref{fig.domain_states} and population the graph with a set of domains. First we declare the domain states following \begin{center} [[-s ,[,] -S [-N ,, -K ]] \\ \end{center} e.g. \begin{center} [[-k folded,unfolded,bell -K 3.3e-4,0.25e-9]] \end{center} This creates a reaction going from the [[folded]] state to the [[unfolded]] state, with the rate constant given by the Bell model (Appendix \ref{sec.bell}). The unfolding parameters are then set to [[3.3e-4,0.25e-9]], which in the case of the Bell model are the unforced rate and transition state distance respectively. \subsubsection{Tension} To access the various tension models from the command line, we define a structure that contains all the neccessary information about the model. <>= typedef struct tension_model_getopt_struct { one_dim_fn_t *handler; /* fn implementing the model on a group */ one_dim_fn_t *inverse_handler; /* fn implementing inverse on group */ char *name; /* name string identifying the model */ char *description; /* a brief description of the model */ int num_params; /* number of extra params for the model */ char **param_descriptions; /* descriptions of extra params */ char *params; /* default values for the extra params */ create_data_func_t *creator; /* param-string -> model-data-struct fn */ destroy_data_func_t *destructor; /* fn to free the model data structure */ } tension_model_getopt_t; @ Set [[inverse_handler]] to [[NULL]] if you wish to invert using bisection. Obviously, this will be slower than computing an analytical inverse and more prone to rounding errors, but it may be easier if a closed-form inverse does not exist. <>= tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = { <>, <>, <>, <>, <>, <> }; @ \subsubsection{Unfolding rate} <>= #define NUM_K_MODELS 6 typedef struct k_model_getopt_struct { char *name; char *description; k_func_t *k; int num_params; char **param_descriptions; char *params; create_data_func_t *creator; destroy_data_func_t *destructor; } k_model_getopt_t; @ <>= k_model_getopt_t k_models[NUM_K_MODELS] = { <>, <>, <>, <>, <>, <> }; @ \subsection{help} <>= void help(char *prog_name, double P, double dF, double t_max, double dt_max, double v, double x_step, state_t *stop_state, environment_t *env, int n_tension_models, tension_model_getopt_t *tension_models, int n_k_models, k_model_getopt_t *k_models, int k_model, list_t *state_list) { state_t *state=NULL; int i, j; if (state_list != NULL) state = S(tail(state_list)); printf("usage: %s [options]\n", prog_name); printf("Version %s\n\n", VERSION); printf("Monte Carlo simulation of a multi-globular domain protein unfolding\n\n"); printf("Simulation options:\n"); printf(""); printf("-q state\tStop simulation when this state empties (- for no limit, currently %s)\n", stop_state ? stop_state->name : "-" ); printf("-t tmax\tMaximum simulated time (-1 for no limit, currently %g s)\n", t_max); printf("-d dt\tMaximum allowed timestep dt (currently %g s)\n", dt_max); printf("-P P\tMaximum transition probability for dt (currently %g)\n", P); printf("-F dF\tMaximum change in force over dt. Only relevant if dx > 0. (currently %g N)\n", dF); printf("-v v\tPulling velocity v (currently %g nm/s)\n", v); printf("-x dx\tPulling step size (currently %g m)\n", x_step); printf("\tWhen dx = 0, the pulling is continuous.\n"); printf("\tWhen dx > 0, the pulling occurs in discrete steps.\n"); printf("Environmental options:\n"); printf("-T T\tTemperature T (currently %g K)\n", env->T); printf("-C T\tYou can also set the temperature T in Celsius\n"); printf("State creation options:\n"); printf("-s name,model[[,group],params]\tCreate a new state.\n"); printf("Once you've set up the state, you may populate it with domains\n"); printf("-N n\tAdd n domains in the current state (%s)\n", state ? state->name : "-"); printf("Transition creation options:\n"); printf("-k initial_state,final_state,model[,params]\tCreate a new transition\n"); printf("To double check your reaction graph, you can produce graphviz dot output\n"); printf("describing the current states and transitions.\n"); printf("-D\tPrint dot description of states and transitions and exit.\n"); printf("Simulation output mode options:\n"); printf("There are two output modes. In standard mode, only the transition\n"); printf("events are printed. For example:\n"); printf(" #Force (N)\tInitial state\tFinal state\n"); printf(" 123.456e-12\tfolded\tunfolded\n"); printf(" ...\n"); printf("In verbose mode, the entire Force-vs-distance curve is output:\n"); printf(" #Distance (m)\tForce (N)\tStiffness (N/m)\n"); printf(" 0.001\t0.0023\t0.002\n"); printf(" ...\n"); printf("-V\tChange output to verbose mode.\n"); printf("-h\tPrint this help and exit.\n"); printf("\n"); printf("Tension models:\n"); for (i=0; i>= <> void get_options(int argc, char **argv, double *pP, double *pDF, double *pT_max, double *pDt_max, double *pV, double *pX_step, state_t **pStop_state, environment_t *env, int n_tension_models, tension_model_getopt_t *tension_models, int n_k_models, k_model_getopt_t *k_models, list_t **pState_list, list_t **pTransition_list) { char *prog_name = NULL; char c, options[] = "q:t:d:P:F:v:x:T:C:s:N:k:DVh"; int tension_model=0, k_model=0, initial_state=-1, final_state=-1; char *state_name=NULL; int i, n, tension_group=0; list_t *temp_list=NULL; int num_strings; char **strings; void *model_params; /* for INIT_MODEL() */ int num_param_args; /* for INIT_MODEL() */ char **param_args; /* for INIT_MODEL() */ extern char *optarg; extern int optind, optopt, opterr; assert(argc > 0); #ifdef DEBUG for (i=0; iT = 300.0; /* K */ *pState_list = NULL; *pTransition_list = NULL; while ((c=getopt(argc, argv, options)) != -1) { switch(c) { case 'q': if (STRMATCH(optarg, "-")) { *pStop_state = NULL; } else { temp_list = head(*pState_list); while (temp_list != NULL) { if (STRMATCH(S(temp_list)->name, optarg)) { *pStop_state = S(temp_list); break; } temp_list = temp_list->next; } } break; case 't': *pT_max = safe_strtod(optarg, "-t"); break; case 'd': *pDt_max = safe_strtod(optarg, "-d"); break; case 'P': *pP = safe_strtod(optarg, "-P"); break; case 'F': *pDF = safe_strtod(optarg, "-F"); break; case 'v': *pV = safe_strtod(optarg, "-v"); break; case 'x': *pX_step = safe_strtod(optarg, "-x"); break; case 'T': env->T = safe_strtod(optarg, "-T"); break; case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break; case 's': parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings); assert(num_strings >= 2 && num_strings <= 4); state_name = strings[0]; if (state_index(*pState_list, state_name) != -1) { fprintf(stderr, "Invalid option '%s', repeated state name '%s'\n", optarg, state_name); exit(1); } tension_model = index_tension_model(n_tension_models, tension_models, strings[1]); if (num_strings == 4) tension_group = safe_strtoi(strings[2], optarg); else tension_group = 0; if (num_strings >= 3) tension_models[tension_model].params = strings[num_strings-1]; #ifdef DEBUG fprintf(stderr, "%s:\t%s -> state %s, model %s, params %s, group %d\n", __FUNCTION__, optarg, state_name, tension_models[tension_model].name, tension_models[tension_model].params, tension_group); #endif INIT_MODEL(state_name, &tension_models[tension_model], tension_models[tension_model].params, model_params); push(pState_list, create_state(state_name, tension_models[tension_model].handler, tension_models[tension_model].inverse_handler, tension_group, model_params, tension_models[tension_model].destructor, 0), 1); *pState_list = tail(*pState_list); state_name = NULL; free_parsed_list(num_strings, strings); break; case 'N': n = safe_strtoi(optarg, "-N"); #ifdef DEBUG fprintf(stderr, "%s:\tadding %d domains to %s\n", __FUNCTION__, n, S(*pState_list)->name); #endif S(*pState_list)->num_domains += n; break; case 'k': parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings); assert(num_strings >= 3 && num_strings <= 4); initial_state = state_index(*pState_list, strings[0]); final_state = state_index(*pState_list, strings[1]); k_model = index_k_model(n_k_models, k_models, strings[2]); #ifdef DEBUG fprintf(stderr, "%s:\t%s -> %s reaction from state %s (%d) -> state %s (%d)\n", __FUNCTION__, optarg, strings[2], strings[0], initial_state, strings[1], final_state); #endif assert(initial_state != final_state); if (num_strings == 4) k_models[k_model].params = strings[3]; INIT_MODEL(state_name, &k_models[k_model], k_models[k_model].params, model_params); push(pTransition_list, create_transition(list_element(*pState_list, initial_state), list_element(*pState_list, final_state), k_models[k_model].k, model_params, k_models[k_model].destructor), 1); free_parsed_list(num_strings, strings); break; case 'D': print_state_graph(stdout, *pState_list, *pTransition_list); exit(0); break; case 'V': flags = FLAG_OUTPUT_FULL_CURVE; break; case '?': fprintf(stderr, "unrecognized option '%c'\n", optopt); /* fall through to default case */ default: help(prog_name, *pP, *pDF, *pT_max, *pDt_max, *pV, *pX_step, *pStop_state, env, n_tension_models, tension_models, n_k_models, k_models, k_model, *pState_list); exit(1); } } /* check the arguments */ assert(*pState_list != NULL); /* no domains! */ assert(*pP > 0.0); assert(*pP < 1.0); assert(*pDt_max > 0.0); assert(*pV > 0.0); assert(env->T > 0.0); crosslink_groups(*pState_list, *pTransition_list); return; } @ <>= int index_tension_model(int n_models, tension_model_getopt_t *models, char *name) { int i; for (i=0; i>= int index_k_model(int n_models, k_model_getopt_t *models, char *name) { int i; for (i=0; i>= /* requires int num_param_args and char **param_args in the current scope * usage: * INIT_MODEL("folded", folded_model, folded_param_string, folded_params); * defined as a macro, so it can work on both tension_model_getopt_t and k_model_getopt_t types. */ #define INIT_MODEL(role, model, param_string, param_pointer) \ do { \ parse_list_string((param_string), SEP, '{', '}', \ &num_param_args, ¶m_args); \ if (num_param_args != (model)->num_params) { \ fprintf(stderr, \ "%s model %s expected %d params,\n", \ (role), (model)->name, (model)->num_params); \ fprintf(stderr, \ "not the %d params in '%s'\n", \ num_param_args, (param_string)); \ assert(num_param_args == (model)->num_params); \ } \ if ((model)->creator) \ param_pointer = (*(model)->creator)(param_args); \ else \ param_pointer = NULL; \ free_parsed_list(num_param_args, param_args); \ } while (0); @ First we define some safe string parsing functions. Currently these abort on any irregularity, but in the future they could print error messages, etc. [[static]] because the functions are currently defined in each [[*.c]] file that uses them, but perhaps they should be moved to [[utils.h/c]] or some such instead\ldots <>= static int safe_strtoi(const char *s, const char *description) { char *endp = NULL; long int ret; assert(s != NULL); ret = strtol(s, &endp, 10); if (endp[0] != '\0') { //strlen(endp) > 0) { fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n", endp, s, description); assert(1==0); //strlen(endp) == 0); } if (endp == s) { fprintf(stderr, "Error: unparsable empty string '%s' for %s\n", s, description); assert(endp != s); } else if (errno == ERANGE) { fprintf(stderr, "Error: '%s' out of range for %s\n", s, description); assert(errno != ERANGE); } return (int) ret; } static double safe_strtod(const char *s, const char *description) { char *endp = NULL; double ret; assert(s != NULL); ret = strtod(s, &endp); if (endp[0] != '\0') { //strlen(endp) > 0) { fprintf(stderr, "Error: unparsable '%s' while parsing '%s' for %s\n", endp, s, description); assert(1==0); //strlen(endp) == 0); } if (endp == s) { fprintf(stderr, "Error: unparsable empty string '%s' for %s\n", s, description); assert(endp != s); } else if (errno == ERANGE) { fprintf(stderr, "Error: '%s' out of range for %s\n", s, description); assert(errno != ERANGE); } return ret; } @ \phantomsection \appendix \addcontentsline{toc}{section}{Appendicies} \section{sawsim.c details} \subsection{Layout} The general layout of our simulation code is: <<*>>= //#define DEBUG <> <> <> <> <> <
> @ We include [[math.h]], so don't forget to link to the libm with `-lm'. <>= #include /* assert() */ #include /* malloc(), free(), rand(), strto*() */ #include /* fprintf(), stdout */ #include /* strlen, strtok() */ #include /* exp(), M_PI, sqrt() */ #include /* pid_t (returned by getpid()) */ #include /* getpid() (for seeding rand()), getopt() */ #include /* errno, ERANGE (for safe_strto*()) */ #include "global.h" #include "list.h" #include "tension_balance.h" #include "k_model.h" #include "tension_model.h" #include "parse.h" #include "accel_k.h" @ <>= <> <> <> <> <> <> <> <> @ <>= <> <> @ <>= <> <> <> <> <> <> <> <> <> @ <>= <> <> @ <>= void setup(void) { srand(getpid()*time(NULL)); /* seed rand() */ } @ <>= /* in octal b/c of prefixed '0' */ #define FLAG_OUTPUT_FULL_CURVE 01 #define FLAG_OUTPUT_TRANSITION_FORCES 02 @ <>= static unsigned long int flags = 0; @ \subsection{Utilities} \label{sec.utils} <>= #define MAX(a,b) ((a)>(b) ? (a) : (b)) #define MIN(a,b) ((a)<(b) ? (a) : (b)) @ Note that [[STRMATCH]] chokes if one of the strings is [[NULL]]. <>= // Check if two strings match, return 1 if they do static char *temp_string_A; static char *temp_string_B; #define STRMATCH(a,b) (temp_string_A=a, temp_string_B=b, \ strlen(temp_string_A) != strlen(temp_string_B) ? 0 : \ !strncmp(temp_string_A,temp_string_B,strlen(temp_string_B)+1) ) /* +1 to also compare the '\0' */ @ We also define a macro for our [[check]] unit testing <>= #define CHECK_ERR(max_err, expected, received) \ do { \ fail_unless( (received-expected)/expected < max_err, \ "relative error %g >= %g in %s (Expected %g, received %g)", \ (received-expected)/expected, max_err, #received, \ expected, received); \ fail_unless(-(received-expected)/expected < max_err, \ "relative error %g >= %g in %s (Expected %g, received %g)", \ -(received-expected)/expected, max_err, #received, \ expected, received); \ } while(0) @ <>= int happens(double probability) { assert(probability >= 0.0); assert(probability <= 1.0); 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*/ } @ \subsection{Adaptive timesteps} \label{sec.adaptive_dt_implementation} $F(x)$ increases with $x$, possibly exploding (e.g.\ in the worm-like chain model), so basing the timestep on the unfolding probability at the current tension is dangerous, and we need to search for a $dt$ for which $P_N(F(x+v*dt))>= double search_dt(transition_t *transition, list_t *state_list, environment_t *env, double x, double max_prob, double max_dt, double v) { /* Returns a valid timestep dt in seconds for the given transition. * Takes a list of all states, a pointer env to the environmental data, * a starting separation x in m, a max_prob between 0 and 1, * max_dt in s, stretching velocity v in m/s. */ double F, k, P, dtCur, dtU, dtUCur, dtL, dt; /* get upper bound using the current position */ F = find_tension(state_list, env, x, 0, 1); /* BUG. repeated calculation */ //printf("Start with x = %g (F = %g)\n", x, F); k = accel_k(transition->k, F, env, transition->k_params); //printf("x %g\tF %g\tk %g\n", x, F, k); dtU = P1(max_prob, N_INIT(transition)) / k; /* upper bound on dt */ if (dtU > max_dt) { //printf("overshot max_dt\n"); dtU = max_dt; } if (v == 0) /* discrete stepping, so force is temporarily constant */ return dtU; /* set a lower bound on dt too */ dtL = 0.0; /* The dt determined above may produce illegitimate forces or ks. * Reduce the upper bound until we have valid ks. */ dt = dtU; F = find_tension(state_list, env, x+v*dt, 0, 1); while (F == HUGE_VAL) { /* reduce step until we hit a valid force */ dtU /= 2.0; dt = dtU; F = find_tension(state_list, env, x+v*dt, 0, 1); } //printf("Try for dt = %g (F = %g)\n", dt, F); k = accel_k(transition->k, F, env, transition->k_params); /* returned k may be -1.0 */ //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k); while (k == -1.0) { /* reduce step until we hit a valid k */ dtU /= 2.0; dt = dtU; /* hopefully, we can use the max dt, see if it works */ F = find_tension(state_list, env, x+v*dt, 0, 1); //printf("Try for dt = %g (F = %g)\n", dt, F); k = accel_k(transition->k, F, env, transition->k_params); //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k); } assert(dtU > 1e-14); /* timestep to valid k too small */ /* safe timestep back from x+dtU */ dtUCur = P1(max_prob, N_INIT(transition)) / k; if (dtUCur >= dt) return dt; /* dtU is safe. */ /* dtU wasn't safe, lets see what would be. */ while (dtU > 1.1*dtL) { /* until the range is reasonably small */ dt = (dtU + dtL) / 2.0; F = find_tension(state_list, env, x+v*dt, 0, 1); //printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU); k = accel_k(transition->k, F, env, transition->k_params); dtCur = P1(max_prob, N_INIT(transition)) / k; //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\tdtCur = %g\n", x, F, dt, v*dt, k, dtCur); if (dtCur > dt) /* safe timestep back from x+dt covers dt */ dtL = dt; else if (dtCur < dt) { /* unsafe timestep back from x+dt, but... */ dtU = dt; /* ... stepping out only dtCur would be safe */ dtUCur = dtCur; } else break; /* dtCur = dt */ } return MAX(dtUCur, dtL); } @ Determine $dt$ for an array of potentially different folded domains. We apply [[search_dt()]] to each populated domain to find the maximum $dt$ the domain can handle. The final $dt$ is less than the individual domain $dt$s (to satisfy [[search_dt()]], see above). If $v>0$ (continuous pulling), $dt$ is also reduced to satisfy $|F(x+v*dt)-F(x)|>= <> double determine_dt(list_t *state_list, list_t *transition_list, environment_t *env, double F, double x, double max_dF, double max_prob, double dt_max, double v, int transition) { /* Returns the timestep dt in seconds. * Takes the state and transition lists, a pointer to the environment, * the force F(x) in N, an extension x in m, max_dF in N, * max_prob between 0 and 1, dt_max in s, v in m/s, and transition in {0,1}*/ double dt=dt_max, new_dt, F_new; assert(max_dF > 0.0); assert(max_prob > 0.0); assert(max_prob < 1.0); assert(dt_max > 0.0); assert(state_list != NULL); assert(transition_list != NULL); transition_list = head(transition_list); if (v > 0) { /* Ensure |dF| < max_dF */ F_new = find_tension(state_list, env, x+v*dt, transition, 1); while (fabs(F_new - F) > max_dF) { dt /= 2.0; F_new = find_tension(state_list, env, x+v*dt, transition, 1); } } /* Ensure all individual domains pass search_dt() */ while (transition_list != NULL) { if (T(transition_list)->initial_state->num_domains > 0) { new_dt = search_dt(T(transition_list),state_list,env,x,max_prob,dt,v); dt = MIN(dt, new_dt); } transition_list = transition_list->next; } return dt; } @ \subsection{State data} \label{sec.state_data} The domains exist in one of an array of possible states. Each state has a tension model which determines it's force/extension curve (possibly [[null]] for rigid states, see Appendix \ref{sec.null_tension}). Groups are, for example, for two domain states with WLC tensions; one with $p=0.4\U{nm}$, $L=10\U{nm}$ and the other with $p=0.4\U{nm}$, $L=28\U{nm}$. Since the persistence lengths are the same, you would like to add the contour lengths of all the domains in \emph{both} states, and plug that total contour length into a single WLC evaluation (see Section \ref{sec.find_tension}). <>= typedef struct state_struct { char *name; /* name string */ one_dim_fn_t *tension_handler; /* tension handler */ one_dim_fn_t *inverse_tension_handler; /* inverse tension handler */ int tension_group; /* for grouping similar tension states */ void *tension_params; /* pointer to tension parameters */ destroy_data_func_t *destroy_tension; int num_domains; /* number of domains currently in state */ /* cached lookups generated from the above data */ int num_out_trans; /* length of out_trans array */ struct transition_struct **out_trans; /* transitions leaving this state */ int num_grouped_states; /* length of grouped states array */ struct state_struct **grouped_states; /* states grouped with this one (doesn't include self) */ } state_t; /* get the state data for the current list node */ #define S(list) ((state_t *)(list)->d) @ [[tension_params]] is a pointer to the parameters used by the handler function when determining the tension. [[destroy_tension]] points to a function for freeing [[tension_params]]. We it in the state structure so that [[destroy_state]] and will have an easier job cleaning up. [[create_]] and [[destroy_state]] are simple wrappers around [[malloc]] and [[free]]. <>= state_t *create_state(char *name, one_dim_fn_t *tension_handler, one_dim_fn_t *inverse_tension_handler, int tension_group, void *tension_params, destroy_data_func_t *destroy_tension, int num_domains) { state_t *ret = (state_t *)malloc(sizeof(state_t)); assert(ret != NULL); /* make a copy of the name, so we aren't dependent on the original */ ret->name = (char *)malloc(sizeof(char)*(strlen(name)+1)); assert(ret->name != NULL); strcpy(ret->name, name); /* we know ret->name is long enough */ ret->tension_handler = tension_handler; ret->inverse_tension_handler = inverse_tension_handler; ret->tension_group = tension_group; ret->tension_params = tension_params; ret->destroy_tension = destroy_tension; ret->num_domains = num_domains; ret->num_out_trans = 0; ret->out_trans = NULL; ret->num_grouped_states = 0; ret->grouped_states = NULL; #ifdef DEBUG fprintf(stderr, "generated state %s (%p) with tension handler %p, inverse handler %p, params %p, group %d.\n", ret->name, ret, ret->tension_handler, ret->inverse_tension_handler, ret->tension_params, ret->tension_group); #endif return ret; } void destroy_state(state_t *state) { if (state) { #ifdef DEBUG fprintf(stderr, "freeing state %s (%p)\n", state->name, state); #endif if (state->name) free(state->name); if (state->destroy_tension) (*state->destroy_tension)(state->tension_params); if (state->out_trans) free(state->out_trans); if (state->grouped_states) free(state->grouped_states); free(state); } } @ We'll be storing the states in a list (see Appendix \ref{sec.lists}), so we also define a convenience function for cleaning up. <>= void destroy_state_list(list_t *state_list) { state_list = head(state_list); while (state_list != NULL) destroy_state((state_t *) pop(&state_list)); } @ We can also define a convenient way to get a state index from it's name. <>= int state_index(list_t *state_list, char *name) { int i=0; state_list = head(state_list); while (state_list != NULL) { if (STRMATCH(S(state_list)->name, name)) return i; state_list = state_list->next; i++; } return -1; } @ \subsection{Transition data} \label{sec.transition_data} Once you've created states (Sections \ref{sec.main}, \ref{sec.model_selection}, and \ref{sec.state_data}), you need to create transitions between them. <>= typedef struct transition_struct { state_t *initial_state; state_t *final_state; k_func_t *k; /* transition rate model */ void *k_params; /* pointer to k parameters */ destroy_data_func_t *destroy_k; } transition_t; /* get the transition data for the current list node */ #define T(list) ((transition_t *)(list)->d) /* get the number of initial state population for a transition state */ #define N_INIT(transition) ((transition)->initial_state->num_domains) <> <> @ [[k]] is a pointer to the function determining the transition rate for a given tension. [[k_params]] and [[destroy_k]] have similar roles with regards to [[k]] as [[tension_params]] and [[destroy_tension]] have with regards to [[tension_handler]] (see Section \ref{sec.state_data}). [[create_]] and [[destroy_transition]] are simple wrappers around [[malloc]] and [[free]]. <>= transition_t *create_transition(state_t *initial_state, state_t *final_state, k_func_t *k, void *k_params, destroy_data_func_t *destroy_k) { transition_t *ret = (transition_t *)malloc(sizeof(transition_t)); assert(ret != NULL); assert(initial_state != NULL); assert(final_state != NULL); ret->initial_state = initial_state; ret->final_state = final_state; ret->k = k; ret->k_params = k_params; ret->destroy_k = destroy_k; #ifdef DEBUG fprintf(stderr, "generated transition %p from %s to %s, params %p\n", ret, ret->initial_state->name, ret->final_state->name, ret->k_params); #endif return ret; } void destroy_transition(transition_t *transition) { if (transition) { #ifdef DEBUG fprintf(stderr, "freeing transition %p (%s to %s)\n", transition, transition->initial_state->name, transition->final_state->name); #endif if (transition->destroy_k) (*transition->destroy_k)(transition->k_params); free(transition); } } @ We'll be storing the transitions in a list (see Appendix \ref{sec.lists}), so we also define a convenience function for cleaning up. <>= void destroy_transition_list(list_t *transition_list) { transition_list = head(transition_list); while (transition_list != NULL) destroy_transition((transition_t *) pop(&transition_list)); } @ \subsection{Graphviz output} \label{sec.graphviz_output} <>= void print_state_graph(FILE *file, list_t *state_list, list_t *transition_list) { state_list = head(state_list); transition_list = head(transition_list); fprintf(file, "digraph states {\n node [shape=record,color=black,fillcolor=white,style=filled];\n\n"); while (1) { fprintf(file, " n%d [label=\"%s\"]\n", state_index(state_list, S(state_list)->name), S(state_list)->name); if (state_list->next == NULL) break; state_list = state_list->next; } fprintf(file, "\n"); while (transition_list != NULL) { fprintf(file, " n%d -> n%d;\n", state_index(state_list, T(transition_list)->initial_state->name), state_index(state_list, T(transition_list)->final_state->name)); transition_list = transition_list->next; } fprintf(file, "}\n"); } @ \subsection{Domain model and group handling} <>= <> <> @ Scan through a list of states and construct an array of all the active groups. <>= void get_active_groups(list_t *state_list, int *pNum_active_groups, one_dim_fn_t ***pPer_group_handlers, one_dim_fn_t ***pPer_group_inverse_handlers, void ***pPer_group_data /* pointer to array of pointers to tension_handler_data_t */) { void **active_groups=NULL; one_dim_fn_t *handler, *old_handler, *inverse_handler, *old_inverse_handler; one_dim_fn_t **per_group_handlers=NULL, **per_group_inverse_handlers=NULL; void **per_group_data=NULL; /* array of pointers to tension_handler_data_t */ int i, j, num_domains, group, old_group, num_active_groups=0; list_t *active_states_list=NULL; tension_handler_data_t *tdata=NULL; state_t *state; /* get all the active groups in a list */ state_list = head(state_list); #ifdef DEBUG fprintf(stderr, "%s:\t", __FUNCTION__); list_print(stderr, state_list, "states"); #endif while (state_list != NULL) { state = S(state_list); handler = state->tension_handler; inverse_handler = state->inverse_tension_handler; group = state->tension_group; num_domains = state->num_domains; if (list_index(active_states_list, state) == -1) { /* we haven't added this state (or it's associated group) yet */ if (num_domains > 0 && handler != NULL) { /* add the state */ num_active_groups++; push(&active_states_list, state, 1); #ifdef DEBUG fprintf(stderr, "%s:\tactive %s state (%p) with %d domain(s)\n", __FUNCTION__, state->name, state, state->num_domains); #endif for (i=0; i < state->num_grouped_states; i++) { if (state->grouped_states[i]->num_domains > 0) { /* add this related (and active) state now */ assert(state->grouped_states[i]->tension_handler == handler); assert(state->grouped_states[i]->tension_group == group); push(&active_states_list, state->grouped_states[i], 1); #ifdef DEBUG fprintf(stderr, "%s:\tactive grouped %s state (%p) with %d domain(s)\n", __FUNCTION__, state->grouped_states[i]->name, state->grouped_states[i], state->grouped_states[i]->num_domains); #endif } } } /* else empty state or NULL handler */ } /* else state already added */ state_list = state_list->next; } #ifdef DEBUG fprintf(stderr, "%s:\t", __FUNCTION__); list_print(stderr, active_states_list, "active states"); #endif assert(num_active_groups <= list_length(active_states_list)); assert(num_active_groups > 0); /* allocate the arrays we'll be returning */ per_group_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups); assert(per_group_handlers != NULL); per_group_inverse_handlers = (one_dim_fn_t **)malloc(sizeof(one_dim_fn_t *)*num_active_groups); assert(per_group_inverse_handlers != NULL); per_group_data = (void **)malloc(sizeof(void *)*num_active_groups); assert(per_group_data != NULL); #ifdef DEBUG fprintf(stderr, "%s:\tallocating per group data %p:", __FUNCTION__, per_group_data); #endif for (i=0; itension_handler; inverse_handler = state->inverse_tension_handler; group = state->tension_group; num_domains = state->num_domains; if (handler != old_handler || inverse_handler != old_inverse_handler || group != old_group) { j++; /* move to the next group */ tdata = (tension_handler_data_t *) per_group_data[j]; per_group_handlers[j] = handler; per_group_inverse_handlers[j] = inverse_handler; tdata->group_tension_params = NULL; tdata->env = NULL; tdata->persist = NULL; /* we might want to grab any appropriate data from the old persist here... */ } #ifdef DEBUG fprintf(stderr, "%s:\tadding active state %s to group %d (%p:%p:%d) with %d domain(s)\n", __FUNCTION__, state->name, j, handler, inverse_handler, group, num_domains); #endif for (i=0; igroup_tension_params, state->tension_params, 1); } #ifdef DEBUG fprintf(stderr, "%s:\tpushed %d copies of %p onto data %p's param list %p\n", __FUNCTION__, num_domains, state->tension_params, tdata, tdata->group_tension_params); #endif old_handler = handler; old_inverse_handler = inverse_handler; old_group = group; } /* free old groups */ if (*pPer_group_handlers != NULL) free(*pPer_group_handlers); if (*pPer_group_inverse_handlers != NULL) free(*pPer_group_inverse_handlers); if (*pPer_group_data != NULL) { for (i=0; i<*pNum_active_groups; i++) free((*pPer_group_data)[i]); free(*pPer_group_data); } /* set new groups */ #ifdef DEBUG fprintf(stderr, "%s:\t%d groups, tension handler array %p, inverse tension handler array %p, handler data array %p (lead data element %p)\n", __FUNCTION__, num_active_groups, per_group_handlers, per_group_inverse_handlers, per_group_data, per_group_data[0]); #endif *pNum_active_groups = num_active_groups; *pPer_group_handlers = per_group_handlers; *pPer_group_inverse_handlers = per_group_inverse_handlers; *pPer_group_data = per_group_data; } @ <>= void crosslink_groups(list_t *state_groups, list_t *transition_list) { list_t *list, *out_trans=NULL, *associates=NULL; one_dim_fn_t *handler; int i, n, group; state_groups = head(state_groups); while (state_groups != NULL) { /* find transitions leaving this state */ #ifdef DEBUG fprintf(stderr, "finding transitions leaving %s\n", S(state_groups)->name); #endif list = head(transition_list); while (list != NULL) { if (T(list)->initial_state == S(state_groups) && list_index(associates, T(list)->initial_state) != -1) { push(&out_trans, T(list), 1); } list = list->next; } n = list_length(out_trans); S(state_groups)->num_out_trans = n; S(state_groups)->out_trans = (transition_t **)malloc(sizeof(transition_t *)*n); assert(S(state_groups)->out_trans != NULL); for (i=0; iout_trans[i] = (transition_t *)pop(&out_trans); } /* find states grouped with this state */ #ifdef DEBUG fprintf(stderr, "finding states grouped with %s\n", S(state_groups)->name); #endif handler = S(state_groups)->tension_handler; group = S(state_groups)->tension_group; list = head(state_groups); while (list != NULL) { if (S(list)->tension_handler == handler && S(list)->tension_group == group && list != state_groups) { push(&associates, S(list), 1); } list = list->next; } n = list_length(associates); S(state_groups)->num_grouped_states = n; S(state_groups)->grouped_states = (state_t **)malloc(sizeof(state_t *)*n); assert(S(state_groups)->grouped_states != NULL); for (i=0; igrouped_states[i] = (state_t *)pop(&associates); } state_groups = state_groups->next; } } @ \section{String parsing} For handling command line arguments, we need parse delimited strings (e.g.\ [["param1,param2,param3..."]]). The model handling in getopt is set up to handle a fixed number of arguments for each model, so models (like [[kramers_integ]]) that have complicated parameters (location of spline knots) need to take care of parsing those parameters themselves. We implement this parsing in [[parse.c]], define the interface in [[parse.h]], and the the unit testing in [[check_parse.c]]. <>= #ifndef PARSE #define PARSE <> <> <> #endif /* PARSE */ @ <>= NW_SAWSIM_MODS += parse CHECK_BINS += check_parse check_parse_MODS = @ <>= #define SEP ',' /* argument separator character */ @ <>= extern void parse_list_string(char *string, char sep, char deeper, char shallower, int *num, char ***string_array); extern void free_parsed_list(int num, char **string_array); @ [[parse_list_string]] allocates memory, don't forget to free it afterward with [[free_parsed_list]]. It does not alter the original. The string may start off with a [[deeper]] character (i.e. [["{x,y}"]]), and when it does, brace stripping will set leave [[x,y]], where the pointer is one character in on the copied string. However, when we go to free the memory, we need a pointer to the beginning of the string. In order to accommodate this for a string with $N$ argument, allocate a pointer array with $N+1$ elements, let the first $N$ elements point to the separated arguments, and let the last element point to the start of the copied string regardless of braces. <>= /* TODO, split out into parse.hc */ static int next_delim_index(char *string, char sep, char deeper, char shallower) { int i=0, depth = 0; while (string[i] != '\0' && !(string[i] == sep && depth == 0)) { if (string[i] == deeper) {depth++;} else if (string[i] == shallower) {depth--; assert(depth >= 0);} i++; } return i; } void parse_list_string(char *string, char sep, char deeper, char shallower, int *num, char ***string_array) { char *str=NULL, **ret=NULL; int i, j, n; if (string==NULL || strlen(string) == 0) { /* handle the trivial cases */ *num = 0; *string_array = NULL; return; } /* make a copy of the string, so we don't change the original */ str = (char *)malloc(sizeof(char)*(strlen(string)+1)); assert(str != NULL); strcpy(str, string); /* we know str is long enough */ /* count the number of regions, so we can allocate pointers to them */ i=-1; n=0; do { n++; i++; /* move on to next argument */ i += next_delim_index(str+i, sep, deeper, shallower); //fprintf( stderr, "delim at %d (%d) ('%s' -> '%s')\n", i, (int)str[i], str, str+i+1); fflush(stderr); } while (str[i] != '\0'); ret = (char **)malloc(sizeof(char *)*(n+1)); assert(ret != NULL); /* replace the separators with '\0' & assign pointers */ ret[n] = str; /* point to the front of the copied string */ j=0; ret[0] = str; for(i=1; i>= <> <> #include "parse.h" <> @ <>= #include /* assert() */ #include /* NULL */ #include /* fprintf(), stdout *//*!!*/ #include /* strlen() */ #include "parse.h" @ \subsection{Parsing unit tests} Here we check to make sure the various functions work as expected, using \citetalias{sw:check}. <>= <> <> <> <> <> <
> @ <>= #include /* EXIT_SUCCESS and EXIT_FAILURE */ #include /* printf() */ #include /* assert() */ #include /* strlen() */ <> #include "parse.h" @ <>= <> <> @ <>= Suite *test_suite (void) { Suite *s = suite_create ("parse"); <> <> return s; } @ <>= /* START_TEST(test_next_delim_index) { fail_unless(next_delim_index("", ',', '{', '}')==0, NULL); fail_unless(next_delim_index(",arg,{str,ing},test", ',', '{', '}')==0, NULL); fail_unless(next_delim_index("arg,{str,ing},test", ',', '{', '}')==3, NULL); fail_unless(next_delim_index("{str,ing},test", ',', '{', '}')==9, NULL); fail_unless(next_delim_index("test", ',', '{', '}')==4, NULL); } END_TEST */ START_TEST(test_parse_list_null) { int num_param_args; char **param_args; parse_list_string(NULL, SEP, '{', '}', &num_param_args, ¶m_args); fail_unless(num_param_args == 0, NULL); fail_unless(param_args == NULL, NULL); } END_TEST START_TEST(test_parse_list_single_simple) { int num_param_args; char **param_args; parse_list_string("arg", SEP, '{', '}', &num_param_args, ¶m_args); fail_unless(num_param_args == 1, NULL); fail_unless(STRMATCH(param_args[0],"arg"), NULL); } END_TEST START_TEST(test_parse_list_single_compound) { int num_param_args; char **param_args; parse_list_string("{x,y,z}", SEP, '{', '}', &num_param_args, ¶m_args); fail_unless(num_param_args == 1, NULL); fail_unless(STRMATCH(param_args[0],"x,y,z"), "got '%s', expected '%s'", param_args[0], "x,y,z"); } END_TEST START_TEST(test_parse_list_double_simple) { int num_param_args; char **param_args; parse_list_string("abc,def", SEP, '{', '}', &num_param_args, ¶m_args); fail_unless(num_param_args == 2, NULL); fail_unless(STRMATCH(param_args[0],"abc"), NULL); fail_unless(STRMATCH(param_args[1],"def"), NULL); } END_TEST START_TEST(test_parse_list_compound) { int num_param_args; char **param_args; parse_list_string("abc,{def,ghi}", SEP, '{', '}', &num_param_args, ¶m_args); fail_unless(num_param_args == 2, NULL); fail_unless(STRMATCH(param_args[0],"abc"), NULL); fail_unless(STRMATCH(param_args[1],"def,ghi"), NULL); } END_TEST @ <>= TCase *tc_parse_list_string = tcase_create("parse list string"); @ <>= //tcase_add_test(tc_parse_list_string, test_next_delim_index); tcase_add_test(tc_parse_list_string, test_parse_list_null); tcase_add_test(tc_parse_list_string, test_parse_list_single_simple); tcase_add_test(tc_parse_list_string, test_parse_list_single_compound); tcase_add_test(tc_parse_list_string, test_parse_list_double_simple); tcase_add_test(tc_parse_list_string, test_parse_list_compound); suite_add_tcase(s, tc_parse_list_string); @ \section{Unit tests} Here we check to make sure the various functions work as expected, using \citetalias{sw:check}. <>= <> <> <> <> <> <> <> <> <
> @ <>= #include @ <>= @ <>= <> <> <> <> <> @ <>= Suite *test_suite (void) { Suite *s = suite_create ("sawsim"); <> <> <> <> <> <> <> <> /* tcase_add_checked_fixture(tc_strip_address, setup_strip_address, teardown_strip_address); */ return s; } @ <
>= int main(void) { int number_failed; Suite *s = test_suite(); SRunner *sr = srunner_create(s); srunner_run_all(sr, CK_ENV); number_failed = srunner_ntests_failed(sr); srunner_free(sr); return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE; } @ \subsection{Model tests} Check the searching with [[linear_k]]. Check overwhelming force treatment with the heavyside-esque step [[?]]. Hmm, I'm not really sure what I was doing here... <>= double linear_k(double F, environment_t *env, void *params) { double Fnot = *(double *)params; return Fnot+F; } START_TEST(test_determine_dt_linear_k) { state_t folded={0}; transition_t unfold={0}; environment_t env = {0}; double F[]={0,1,2,3,4,5,6}; double dt_max=3.0, Fnot=3.0, x=1.0, max_p=0.1, max_dF=0.1, v=1.0; list_t *state_list=NULL, *transition_list=NULL; int i; env.T = 300.0; folded.tension_handler = &hooke_handler; folded.num_domains = 1; unfold.initial_state = &folded; unfold.k = linear_k; unfold.k_params = &Fnot; push(&state_list, &folded, 1); push(&transition_list, &unfold, 1); for( i=0; i < sizeof(F)/sizeof(double); i++) { //fail_unless(determine_dt(state_list, transition_list, &env, x, max_p, max_dF, dt_max, v)==1, NULL); } } END_TEST @ <>= TCase *tc_determine_dt = tcase_create("Determine dt"); @ <>= tcase_add_test(tc_determine_dt, test_determine_dt_linear_k); suite_add_tcase(s, tc_determine_dt); @ <>= @ <>= @ <>= @ <>= @ <>= @ <>= @ <>= @ <>= @ <>= @ \section{Balancing group extensions} \label{sec.tension_balance} For a given total extension $x$ (of the piezo), the various domain models (WLC, FJC, Hookean springs, \ldots) and groups extend different amounts, and we need to tweak the portion that each extends to balance the tension amongst the active groups. Since the tension balancing is separable from the bulk of the simulation, we break it out into a separate module. The interface is defined in [[tension_balance.h]], the implementation in [[tension_balance.c]], and the unit testing in [[check_tension_balance.c]] <>= #ifndef TENSION_BALANCE_H #define TENSION_BALANCE_H <> <> #endif /* TENSION_BALANCE_H */ @ <>= <> <> <> <> <> <> <> @ <>= NW_SAWSIM_MODS += tension_balance CHECK_BINS += check_tension_balance check_tension_balance_MODS = @ The entire force balancing problem reduces to a solving two nested one-dimensional root-finding problems. First we define the one dimensional tension $F(x_0)$ (where $i=0$ is the index of the first populated group). $x(x_0)$ is also strictly monotonically increasing, so we can solve for $x_0$ such that $\sum_i x_i = x$. [[last_x]] stores the last successful value of $x$, since we expect to be taking small steps ($x \approx$ [[last_x]]). Setting [[last_x = -1]] indicates that the groups have changed. <>= double tension_balance(int num_tension_groups, one_dim_fn_t **tension_handlers, one_dim_fn_t **inverse_tension_handlers, void **params, /* array of pointers to tension_handler_data_t */ double last_x, double x, double *stiffness); @ <>= double tension_balance(int num_tension_groups, one_dim_fn_t **tension_handlers, one_dim_fn_t **inverse_tension_handlers, void **params, /* array of pointers to tension_handler_data_t */ double last_x, double x, double *stiffness) { static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL,NULL}; double F, xo_guess, xo, lb, ub; double min_relx=1e-6, min_rely=1e-6; int max_steps=200, i; assert(num_tension_groups > 0); assert(tension_handlers != NULL); assert(params != NULL); assert(num_tension_groups > 0); if (last_x == -1 || x_xo_data.xi == NULL) { /* new group setup, reset x_xo_data */ /* x_xo_data tension_handler and handler_data arrays freed by find_tension() */ if (x_xo_data.xi != NULL) free(x_xo_data.xi); /* construct new x_xo_data */ x_xo_data.n_groups = num_tension_groups; x_xo_data.tension_handler = tension_handlers; x_xo_data.inverse_tension_handler = inverse_tension_handlers; x_xo_data.handler_data = params; #ifdef DEBUG fprintf(stderr, "%s:\t%d groups, tension handler array %p, handler data array %p (lead data element %p), x_of_xo_data %p\n", __FUNCTION__, x_xo_data.n_groups, x_xo_data.tension_handler, x_xo_data.handler_data, x_xo_data.handler_data[0], &x_xo_data); #endif x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups); assert(x_xo_data.xi != NULL); for (i=0; i last_x) { lb = x_xo_data.xi[0]; ub = x_xo_data.xi[0]+ x-last_x; /* apply all change to x0 */ } else if (x < last_x) { lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */ ub = x_xo_data.xi[0]; } else { /* x == last_x */ #ifdef DEBUG fprintf(stderr, "not moving\n"); #endif lb = x_xo_data.xi[0]; ub = x_xo_data.xi[0]; } } #ifdef DEBUG fprintf(stderr, "%s:\tsolve x(x0) so x = %g with bounds x0 in [%g, %g]\n", __FUNCTION__, x, lb, ub); #endif xo = oneD_solve(x_of_xo, &x_xo_data, x, lb, ub, min_relx, min_rely, max_steps, NULL); F = x_of_xo(xo, &x_xo_data); /* HACK, bonus call to setup x_xo_data.xi */ #ifdef DEBUG if (num_tension_groups > 1) { fprintf(stderr, "%s:\tbalanced for x = %g with ", __FUNCTION__, x); for(i=0; i>= <> @ <>= typedef struct x_of_xo_data_struct { int n_groups; one_dim_fn_t **tension_handler; /* array of fn pointers */ one_dim_fn_t **inverse_tension_handler; /* inverse fn pointers */ void **handler_data; /* array of pointers to tension_handler_data_t */ double *xi; /* array of group extensions */ } x_of_xo_data_t; @ <>= double x_of_xo(double xo, void *pdata) { x_of_xo_data_t *data = (x_of_xo_data_t *)pdata; double F, x=0, xi, xi_guess, lb, ub, slope; int i; double min_relx=1e-6, min_rely=1e-6; int max_steps=200; assert(data->n_groups > 0); data->xi[0] = xo; #ifdef DEBUG fprintf(stderr, "%s:\tget initial tension at %g with %d groups, handler %p, data %p (x_xo_data %p)\n", __FUNCTION__, xo, data->n_groups, data->tension_handler[0], data->handler_data[0], data); fflush(stderr); #endif F = (*data->tension_handler[0])(xo, data->handler_data[0]); #ifdef DEBUG fprintf(stderr, "%s:\tfinding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F); fflush(stderr); #endif x += xo; if (data->xi) data->xi[0] = xo; for (i=1; i < data->n_groups; i++) { if (data->inverse_tension_handler[i] != NULL) { xi = (*data->inverse_tension_handler[i])(F, data->handler_data[i]); } else { /* invert numerically */ if (data->xi[i] == -1) xi_guess = length_scale; /* HACK */ else xi_guess = data->xi[i]; #ifdef DEBUG fprintf(stderr, "%s:\tsearching for proper x[%d] for F[%d](x[%d]) = %g N (handler %p, data %p)\n", __FUNCTION__, i, i, i, F, data->tension_handler[i], data->handler_data[i]); #endif oneD_bracket(data->tension_handler[i], data->handler_data[i], F, xi_guess, &lb, &ub); xi = oneD_solve(data->tension_handler[i], data->handler_data[i], F, lb, ub, min_relx, min_rely, max_steps, NULL); } x += xi; data->xi[i] = xi; #ifdef DEBUG fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F); #endif } #ifdef DEBUG fprintf(stderr, "%s:\tfound x(x0=%g) = %g\n", __FUNCTION__, xo, x); #endif return x; } @ Determine the stiffness of a single tension group by taking a small step [[dx]] back from the position [[x]] for which we want the stiffness. The touchy part is determining a good [[dx]] and ensuring the whole interval is on [[x>=0]]. <>= double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double x, double relx) { double dx, xl, F, Fl, stiffness; assert(x >= 0); assert(relx > 0 && relx < 1); if (x == 0 || relx == 0) { dx = 1e-12; /* HACK, how to get length scale? */ xl = x-dx; if (xl < 0) { xl = 0; x = xl+dx; } } else { dx = x*relx; xl = x-dx; } F = (*tension_handler)(x, handler_data); Fl = (*tension_handler)(xl, handler_data); stiffness = (F-Fl)/dx; assert(stiffness >= 0); return stiffness; } @ Attempt to calculate the chain stiffness by adding group stiffnesses as springs in series. In the case where a group has undefined stiffness (e.g. piston tension, Section \ref{sec.piston_tension}), this algorithm breaks down. In that case, [[tension_balance()]] (\ref{sec.tension_balance}) should use [[full_chain_stiffness()]] which uses the full chain's $F(x)$ rather than that of the individual domains'. We attempt the series approach first because, lacking the numerical inversion inside [[tension_balance()]], it is both faster and more accurate than the full-chain derivative. <>= double chain_stiffness(x_of_xo_data_t *x_xo_data, double relx) { double stiffness, gstiff; int i; /* add group stiffnesses in series */ for (i=0; i < x_xo_data->n_groups; i++) { #ifdef DEBUG fprintf(stderr, "%s:\tgetting stiffness of active state %d of %d for x[%d]=%g, relx=%g\n", __FUNCTION__, i, x_xo_data->n_groups, i, x_xo_data->xi[i], relx); fflush(stderr); #endif gstiff = group_stiffness(x_xo_data->tension_handler[i], x_xo_data->handler_data[i], x_xo_data->xi[i], relx); if (gstiff == 0.0); return 0.0; stiffness += 1.0/gstiff; } stiffness = 1.0 / stiffness; return stiffness; } @ Determine the chain stiffness using only [[tension_balance()]]. This works around difficulties with tension models that have undefined stiffness (see discussion for [[chain_stiffness()]] above). <>= double full_chain_stiffness(int num_tension_groups, one_dim_fn_t **tension_handlers, one_dim_fn_t **inverse_tension_handlers, void **params, /* array of pointers to tension_handler_data_t */ double last_x, double x, double relx, double F /* already evaluated F(x) */) { double dx, xl, Fl, stiffness; assert(x >= 0); assert(relx > 0 && relx < 1); /* Other option for dx that we ignore because of our poor tension_balance() * resolution (i.e. bad slopes if you zoom in too close): * if (last_x != -1.0 && last_x != x) // last_x limited * dx fabs(x-last_x); * else * dx = HUGE_VAL; * if (1==1) { // constant max-value limited * dx_p = 1e-12; * dx = MIN(dx, dx_p); * } * if (x != 0 && relx != 0) { // relx limited * dx_p = x*relx; * dx = MIN(dx, dx_p); * } * TODO, determine which of (min_relx, min_rely, max_steps) is actually * limiting tension_balance. */ dx = 1e-12; /* HACK, how to get length scale? */ xl = x-dx; if (xl >= 0) { Fl = tension_balance(num_tension_groups, tension_handlers, inverse_tension_handlers, params, last_x, xl, NULL); } else { xl = x; Fl = F; x += dx; F = tension_balance(num_tension_groups, tension_handlers, inverse_tension_handlers, params, last_x, x, NULL); } stiffness = (F-Fl)/dx; assert(stiffness >= 0); return stiffness; } @ Solve $f(x) = y$ to a certain precision in $x$ or $y$ in a limited number of steps for monotonically increasing $f(x)$. Simple bisection, so it's robust and converges fairly quickly. You can set the maximum number of steps to take, but if you have enough processor power, it's better to limit your precision with the relative limits [[min_relx]] and [[y]]. These ensure that the width of the bracket is small on the length scale of it's larger side. Note that \emph{both} relative limits must be satisfied for the search to stop. <>= /* equivalent to gsl_func_t */ typedef double one_dim_fn_t(double x, void *params); @ <>= double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub, double min_relx, double min_rely, int max_steps, double *pYx); @ <>= /* TODO, replace with GSL one d root finder http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Root_002dFinding.html */ double oneD_solve(one_dim_fn_t *f, void *params, double y, double lb, double ub, double min_relx, double min_rely, int max_steps, double *pYx) { double yx, ylb, yub, x; int i=0; assert(ub >= lb); ylb = (*f)(lb, params); yub = (*f)(ub, params); /* check some simple cases */ if (ylb == yub) { if (ylb != y) return HUGE_VAL; /* error! f(x)=y not bracketed */ else return lb; /* any x on the range [lb, ub] would work */ } if (ylb == y) { x=lb; yx=ylb; goto end; } if (yub == y) { x=ub; yx=yub; goto end; } #ifdef DEBUG fprintf(stderr, "%s:\tstarting bracket: lb %g, x ~, ub %g\tylb %g, y %g, yub %g\n", __FUNCTION__, lb, ub, ylb, y, yub); #endif assert(ylb < y); assert(yub > y); for (i=0; i y) { ub=x; yub=yx; } else /* yx < y */ { lb=x; ylb=yx; } } end: if (pYx) *pYx = yx; #ifdef DEBUG fprintf(stderr, "%s:\tsolved f(%g) = %g ~= %g\n", __FUNCTION__, x, yx, y); #endif return x; } @ The one dimensional solver needs a bracketed solution, which sometimes we don't have. Generate bracketing $x$ values through bisection or exponential growth. <>= void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub); @ <>= void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub) { double yx, step, x=xguess; int i=0; #ifdef DEBUG fprintf(stderr, "%s:\tbracketing %g (params %p)\n", __FUNCTION__, y, params); #endif yx = (*f)(x, params); #ifdef DEBUG fprintf(stderr, "%s:\tbracketing %g, start at f(%g) = %g guessing x = %g (params %p)\n", __FUNCTION__, y, x, yx, xguess, params); #endif if (yx > y) { assert(x > 0.0); *ub = x; *lb = 0; #ifdef DEBUG fprintf(stderr, "%s:\tbracketted by [0, guess = %g]\n", __FUNCTION__, x); #endif } else { *lb = x; if (x == 0) x = length_scale/2.0; /* HACK */ while (yx < y) { x *= 2.0; yx = (*f)(x, params); #ifdef DEBUG fprintf(stderr, "%s:\tincreasing to f(%g) = %g\n", __FUNCTION__, x, yx); #endif } *ub = x; #ifdef DEBUG fprintf(stderr, "%s:\tbracketted ouput %g by [guess = %g, %g]\n", __FUNCTION__, y, xguess, x); #endif } //fprintf(stdout, "ub %g, %g lb\n", *ub, *lb); } @ \subsection{Balancing implementation} <>= //#define DEBUG <> <> #include "tension_balance.h" double length_scale = 1e-10; /* HACK */ <> <> @ <>= #include /* assert() */ #include /* NULL */ #include /* HUGE_VAL, macro constant, so don't need to link to math */ #include /* fprintf(), stdout */ #include "global.h" @ \subsection{Balancing unit tests} Here we check to make sure the various functions work as expected, using \citetalias{sw:check}. <>= <> <> <> <
> @ <>= #include /* EXIT_SUCCESS and EXIT_FAILURE */ #include /* assert() */ <> #include "global.h" #include "tension_balance.h" @ <>= <> <> @ <>= Suite *test_suite (void) { Suite *s = suite_create ("tension balance"); <> <> return s; } @ <>= <> double hooke(double x, void *pK) { assert(x >= 0); return *((double*)pK) * x; } START_TEST(test_single_function) { double k=5, x=3, last_x=2.0, F; one_dim_fn_t *handlers[] = {&hooke}; one_dim_fn_t *inverse_handlers[] = {NULL}; void *data[] = {&k}; F = tension_balance(1, handlers, inverse_handlers, data, last_x, x, NULL); fail_unless(F = k*x, NULL); } END_TEST @ We can also test balancing two springs with different spring constants. The analytic solution for a general number of springs is given in Appendix \ref{sec.math_hooke}. <>= START_TEST(test_double_hooke) { double k1=5, k2=4, x=3, last_x=-1.0, F, Fe, x1e, x2e; one_dim_fn_t *handlers[] = {&hooke, &hooke}; one_dim_fn_t *inverse_handlers[] = {NULL, NULL}; void *data[] = {&k1, &k2}; F = tension_balance(2, handlers, inverse_handlers, data, last_x, x, NULL); x1e = x*k2/(k1+k2); Fe = k1*x1e; x2e = x1e*k1/k2; //fail_unless(-(F-Fe)/Fe < 1e-6, "relative error %g > 1e-6 (Expected %g, got %g)",-(F-Fe)/Fe, Fe, F); //CHECK_ERR(1e-6, x1e, xi[0]); //CHECK_ERR(1e-6, x2e, xi[1]); CHECK_ERR(1e-6, Fe, F); } END_TEST @ TODO: test [[tension_balance()]] with another balance after the [[x_of_xo_data]] initializing balance we ran above. <>= TCase *tc_tbfunc = tcase_create("tension balance function"); @ <>= tcase_add_test(tc_tbfunc, test_single_function); tcase_add_test(tc_tbfunc, test_double_hooke); suite_add_tcase(s, tc_tbfunc); @ \section{Lists} \label{sec.lists} The globular domains (and cantilever params, etc.) are saved in bi-directional lists. Since the list handling is so general, and separable from the bulk of the simulation, we break it out into a separate module. The interface is defined in [[list.h]], the implementation in [[list.c]], and the unit testing in [[check_list.c]] <>= #ifndef LIST_H #define LIST_H <> <> <> #endif /* LIST_H */ @ <>= <> <> <> <> <> <> <> <> <> <> <> <> @ <>= <> <> <> <> <> <> <> <> <> <> <> <> <> @ <>= NW_SAWSIM_MODS += list CHECK_BINS += check_list check_list_MODS = @ <>= typedef struct list_struct { struct list_struct *next; struct list_struct *prev; void *d; /* data */ } list_t; <> @ [[head]] and [[tail]] return pointers to the head and tail nodes of the list: <>= list_t *head(list_t *list); list_t *tail(list_t *list); @ <>= list_t *head(list_t *list) { if (list == NULL) return list; while (list->prev != NULL) { list = list->prev; } return list; } list_t *tail(list_t *list) { if (list == NULL) return list; while (list->next != NULL) { list = list->next; } return list; } @ <>= int list_length(list_t *list); @ <>= int list_length(list_t *list) { int i; if (list == NULL) return 0; list = head(list); i = 1; while (list->next != NULL) { list = list->next; i += 1; } return i; } @ [[push]] inserts a node relative to the current position in the list without changing the current position. However, if the list we're pushing onto is [[NULL]], the current position isn't defined so we set it to that of the pushed domain. If below is zero, we insert above the current node (push(a>B>c,d,0) -> a>d>B>c), otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c). <>= void push(list_t **pList, void *data, int below); @ <>= void push(list_t **pList, void *data, int below) { list_t *list, *new_node; assert(pList != NULL); list = *pList; new_node = create_node(data); if (list == NULL) *pList = new_node; else if (below==0) { /* insert above */ if (list->prev != NULL) list->prev->next = new_node; new_node->prev = list->prev; list->prev = new_node; new_node->next = list; } else { /* insert below */ if (list->next != NULL) list->next->prev = new_node; new_node->next = list->next; list->next = new_node; new_node->prev = list; } } @ [[pop]] removes the current domain node, moving the current position to the node after it, unless that node is [[NULL]], in which case move the current position to the node before it. <>= void *pop(list_t **pList); @ <>= void *pop(list_t **pList) { list_t *list, *popped; void *data; assert(pList != NULL); list = *pList; assert(list != NULL); /* not an empty list */ popped = list; /* bypass the popped node */ if (list->prev != NULL) list->prev->next = list->next; if (list->next != NULL) { list->next->prev = list->prev; *pList = list->next; } else *pList = list->prev; /* note: list->prev may == NULL. that's ok ;) */ data = popped->d; destroy_node(popped); return data; } @ [[list_destroy]] just pops all the nodes off a list, possibly excuting a clean-up function on their each node's data. If the cleanup function is [[NULL]], no cleanup function is called. The nodes are not popped in any particular order. <>= void list_destroy(list_t **pList, destroy_data_func_t *destroy); @ <>= void list_destroy(list_t **pList, destroy_data_func_t *destroy) { list_t *list; void *data; assert(pList != NULL); list = *pList; *pList = NULL; assert(list != NULL); /* not an empty list */ while (list != NULL) { data = pop(&list); if (destroy != NULL) destroy(data); } } @ [[list_to_array]] converts a list to an array of pointers, preserving order. <>= void list_to_array(list_t *list, int *length, void ***array); @ <>= void list_to_array(list_t *list, int *pLength, void ***pArray) { void **array; int i,length; assert(list != NULL); assert(pLength != NULL); assert(pArray != NULL); length = list_length(list); array = (void **)malloc(sizeof(void *)*length); assert(array != NULL); list = head(list); i=0; while (list != NULL) { array[i++] = list->d; list = list->next; } *pLength = length; *pArray = array; } @ [[move]] moves the current element along the list in either direction. <>= void move(list_t **pList, int downstream); @ <>= void move(list_t **pList, int downstream) { list_t *list, *flipper; void *data; assert(pList != NULL); list = *pList; /* a>B>c>d */ assert(list != NULL); /* not an empty list */ if (downstream == 0) flipper = list->prev; /* flipper is A */ else flipper = list->next; /* flipper is C */ /* ensure there is somewhere to go in the direction we're heading */ assert(flipper != NULL); /* remove the current element */ data = pop(&list); /* data=B, a>C>d */ /* and put it back in in the correct spot */ list = flipper; if (downstream == 0) { push(&list, data, 0); /* b>A>c>d */ list = list->prev; /* B>a>c>d */ } else { push(&list, data, 1); /* a>C>b>d */ list = list->next; /* a>c>B>d */ } *pList = list; } @ [[sort]] sorts a list from smallest to largest according to a user specified comparison. <>= typedef int (comparison_fn_t)(void *A, void *B); @ For example <>= int int_comparison(void *A, void *B) { int a,b; a = *((int *)A); b = *((int *)B); if (a > b) return 1; else if (a == b) return 0; else return -1; } @ <>= void sort(list_t **pList, comparison_fn_t *comp); @ Here we use the bubble sort algorithm. <>= void sort(list_t **pList, comparison_fn_t *comp) { list_t *list; int stable=0; assert(pList != NULL); list = *pList; assert(list != NULL); /* not an empty list */ while (stable == 0) { stable = 1; list = head(list); while (list->next != NULL) { if ((*comp)(list->d, list->next->d) > 0) { stable = 0; move(&list, 1 /* downstream */); } else list = list->next; } } *pList = list; } @ [[list_index]] finds the location of [[data]] in [[list]]. Returns [[-1]] if [[data]] is not in [[list]]. <>= int list_index(list_t *list, void *data); @ <>= int list_index(list_t *list, void *data) { int i=0; list = head(list); while (list != NULL) { if (list->d == data) return i; list = list->next; i++; } return -1; } @ [[list_element]] returns the data bound to the $i^\text{th}$ node of the list. <>= void *list_element(list_t *list, int i); @ <>= void *list_element(list_t *list, int i) { int j=0; list = head(list); while (list != NULL) { if (i == j) return list->d; list = list->next; j++; } assert(0==1); } @ [[unique]] removes duplicates from a sorted list according to a user specified comparison. Don't do this unless you have other pointers any dynamically allocated data in your list, because [[unique]] just drops any non-unique data without freeing it. <>= void unique(list_t **pList, comparison_fn_t *comp); @ <>= void unique(list_t **pList, comparison_fn_t *comp) { list_t *list; assert(pList != NULL); list = *pList; assert(list != NULL); /* not an empty list */ list = head(list); while (list->next != NULL) { if ((*comp)(list->d, list->next->d) == 0) { pop(&list); } else list = list->next; } *pList = list; } @ [[list_print]] prints a list to a [[FILE *]] stream. <>= void list_print(FILE *file, list_t *list, char *list_name); @ <>= void list_print(FILE *file, list_t *list, char *list_name) { fprintf(file, "%s:", (list_name==NULL) ? "list" : list_name); list = head(list); while (list != NULL) { fprintf(file, " %p", list->d); list = list->next; } fprintf(file, "\n"); } @ [[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]]. <>= list_t *create_node(void *data) { list_t *ret = (list_t *)malloc(sizeof(list_t)); assert(ret != NULL); ret->prev = NULL; ret->next = NULL; ret->d = data; return ret; } void destroy_node(list_t *node) { if (node) free(node); } @ The user must free the data pointed to by the node on their own. \subsection{List implementation} <>= <> <> #include "list.h" <> @ <>= #include /* assert() */ #include /* malloc(), free(), rand() */ #include /* fprintf(), stdout, FILE */ #include "global.h" /* destroy_data_func_t */ @ \subsection{List unit tests} Here we check to make sure the various functions work as expected, using \citetalias{sw:check}. <>= <> <> <> <
> @ <>= #include /* EXIT_SUCCESS and EXIT_FAILURE */ #include /* FILE */ <> #include "global.h" #include "list.h" @ <>= <> <> <> @ <>= Suite *test_suite (void) { Suite *s = suite_create ("list"); <> <> <> <> return s; } @ <>= START_TEST(test_push) { list_t *list=NULL; int i, p, e, n=3; for (i=0; id == 0 ); for (i=0; i>= TCase *tc_push = tcase_create("push"); @ <>= tcase_add_test(tc_push, test_push); suite_add_tcase(s, tc_push); @ <>= @ <>= @ <>= @ \section{Function string evaluation} 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). We want the ability to handle fairly general functions, but don't want to reinvent the wheel by writing our own parser/evaluator. The solution is to run a scripting language as a subprocess accessed via pipes. We use \citetalias{sw:python} here, but it would be a simple matter to replace it with another evaluator if so desired. We spawn the subprocess using the standard [[pipe]], [[fork]], [[exec]] idiom. This is one of the more complicated software ideas in \prog, so we'll go into more detail than we have been. 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. We persevere despite the difficulties, because without command-line access to new functions, the saddle-point Kramers' approximation is of very limited utiltiy. 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. Then you can either statically or dynamically link to those hardcoded functions. While much less flexible, this approach would be a fairly simple-to-implement fallback. Because this functionality is independent of the rest of the simulation, we'll split its definition out into its own set of files. The interface is defined in [[string_eval.h]], the implementation in [[string_eval.c]], and the unit testing in [[check_string_eval.c]]. <>= #ifndef STRING_EVAL_H #define STRING_EVAL_H <> <> <> <> #endif /* STRING_EVAL_H */ @ <>= NW_SAWSIM_MODS += string_eval CHECK_BINS += check_string_eval check_string_eval_MODS = @ For an introduction to POSIX process control, see\\ \url{http://www.ibm.com/developerworks/aix/library/au-speakingunix8/} (very simple, but a good intro.), \\ \url{http://www.ibm.com/developerworks/aix/library/au-unixprocess.html} (more detail), and of course, the relavent [[man]] pages. We start our subprocess with [[execvp]], one of the [[exec]] family of functions. [[execvp]] replaces the calling process' program with a new program. The [[p]] in [[execvp]] signifies that the new program will be found by searching the the paths listed in the [[PATH]] environment variable (this may be a security hole if someone messes about with [[PATH]] before you run \prog, but if they had the ability to change [[PATH]], the security of \prog\ is the least of your worries). The new program needs command line arguments, just like it would if you were running it from a shell. The [[v]] in [[execvp]] signifies that these command line arguments will be provided as an array of [[NULL]] terminated strings, with the final array entry being a [[NULL]] pointer. Now that we know how [[execvp]] works, we store it's arguments in some definitions, to make it easy to change the evaluating subprocess to, say, Ruby, or the users personal evaluation language. <>= #define SUBPROCESS "python" //#define SUBPROCESS_ARGS {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL} static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "import sys;sys.ps1='';sys.ps2=''", (char *)NULL}; //static char *SUBPROCESS_ARGS[] = {SUBPROCESS, "-ic", "pass", (char *)NULL}; @ The [[i]] option lets Python know that it should run in interactive mode. In it's standard mode, python reads all of the supplied instructions, and then acts on them in a single pass. In interactive mode, python acts on every instruction as soon as it is recieved. 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. %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. Since the call to [[execvp]] replaces the calling program, we need to split out original program in half useing [[fork]]. The parent half of the program can continue on to run our simulation, while the child half [[exec]]s and turns into Python. [[fork]] returns two copies of the same program, executing the original code. In the child process [[fork]] returns 0, and in the parent it returns the process ID of the child. We use this difference allows to write seperate parent/child code in the section immediately following the [[fork]]. 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. The [[pipe]] function creates an unnamed pipe, and returns the file descriptors for reading and writing into the new pipe. We need two pipes, one for the subprocess's [[stdin]] and one for its [[stdout]]. We store the pipe file descriptors in an array of 4 [[int]]s, and use the following definitions for access. <>= #define PIPE_READ 0 /* the end you read from */ #define PIPE_WRITE 1 /* the end you write to */ #define STDIN 0 /* initial index of stdin pair */ #define STDOUT 2 /* initial index of stdout pair */ @ So [[pfd[STDIN+PIPE_READ]]] is the file descriptor you would read from for the [[stdin]] pipe, and similarly for the other combinations. 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. <>= extern void string_eval_setup(FILE **pIn, FILE **pOut); @ <>= void string_eval_setup(FILE **pIn, FILE **pOut) { pid_t pid; int pfd[4]; /* file descriptors for stdin and stdout */ int rc; assert(pipe(pfd+STDIN) != -1); /* stdin pair (in, out) */ assert(pipe(pfd+STDOUT) != -1); /* stdout pair (in, out) */ pid = fork(); /* split process into two copies */ if (pid == -1) { /* fork error */ perror("fork error"); exit(1); } else if (pid == 0) { /* child */ close(pfd[STDIN+PIPE_WRITE]); /* close stdin pipe input */ close(pfd[STDOUT+PIPE_READ]); /* close stdout pipe output */ dup2(pfd[STDIN+PIPE_READ],0); /* wire stdin pipe output to stdin (closes original stdin) */ dup2(pfd[STDOUT+PIPE_WRITE],1); /* wire stdout pipe input to stdout (closes original stdout) */ execvp(SUBPROCESS, SUBPROCESS_ARGS); /* run subprocess */ perror("exec error"); /* exec shouldn't return */ _exit(1); } else { /* parent */ close(pfd[STDIN+PIPE_READ]); /* close stdin pipe output */ close(pfd[STDOUT+PIPE_WRITE]); /* close stdout pipe input */ *pIn = fdopen(pfd[STDIN+PIPE_WRITE], "w"); /* 'upgrade' our file descriptors to streams */ if ( *pIn == NULL ) { perror("fdopen (in)"); exit(1); } *pOut = fdopen(pfd[STDOUT+PIPE_READ], "r"); if ( *pOut == NULL ) { perror("fdopen (out)"); exit(1); } } } @ To use the evaluating subprocess, we just pipe in our command, and read out the results. For the simple cases we expect here, we restrict ourselves to a single line of returned text. <>= extern void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output); @ <>= void string_eval(FILE *in, FILE *out, char *input, int buflen, char *output) { int rc; rc = fprintf(in, "%s", input); assert(rc == strlen(input)); fflush(in); fflush(out); alarm(1); /* set a one second timeout on the read */ assert( fgets(output, buflen, out) != NULL ); alarm(0); /* cancel the timeout */ //fprintf(stderr, "eval: %s ----> %s", input, output); } @ The [[alarm]] calls set and clear a timeout on the returned output. 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. This protects against invalid input for which a line of output is not printed to [[stdout]]. Other invalid inputs (e.g.\ those generating multiple lines on [[stdout]] or a single line on [[stdout]] and more in [[stderr]]) are silently ignored. If you are getting strange results, check your python code seperately. TODO, better error handling. Cleanup is fairly straightforward, we just close the connecting streams from the parent side. With the stdin pipe close on the parent side, the reading child will recive the broken pipe signal [[SIGPIPE]], and closes. The parent waits to confirm the child closing, recieves the child's exit status, and cleans up to prevent zombies. As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits. <>= extern void string_eval_teardown(FILE **pIn, FILE **pOut); @ <>= void string_eval_teardown(FILE **pIn, FILE **pOut) { pid_t pid=0; int stat_loc; /* redirect Python's stderr */ fprintf(*pIn, "sys.stderr = open('/dev/null', 'w')\n"); fflush(*pIn); /* close pipes */ assert( fclose(*pIn) == 0 ); *pIn = NULL; assert( fclose(*pOut) == 0 ); *pOut = NULL; /* wait for python to exit */ while (pid <= 0) { pid = wait(&stat_loc); if (pid < 0) { perror("pid"); } } /* if (WIFEXITED(stat_loc)) { printf("child exited with status %d\n", WEXITSTATUS(stat_loc)); } else if (WIFSIGNALED(stat_loc)) { printf("child terminated with signal %d\n", WTERMSIG(stat_loc)); } */ } @ The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals. \subsection{String evaluation implementation} <>= <> <> #include "string_eval.h" <> <> @ <>= #include /* assert() */ #include /* NULL */ #include /* fprintf(), stdout, fdopen() */ #include /* strlen() */ #include /* pid_t */ #include /* pipe(), fork(), execvp(), alarm() */ #include /* wait() */ @ <>= <> <> @ <>= <> <> <> @ \subsection{String evaluation unit tests} Here we check to make sure the various functions work as expected, using \citetalias{sw:check}. <>= <> <> <> <> <
> @ <>= #include /* EXIT_SUCCESS and EXIT_FAILURE */ #include /* printf() */ #include /* assert() */ #include /* strlen() */ #include /* SIGKILL */ <> #include "string_eval.h" @ <>= <> <> @ <>= Suite *test_suite (void) { Suite *s = suite_create ("string eval"); <> <> return s; } @ <>= START_TEST(test_setup_teardown) { FILE *in, *out; string_eval_setup(&in, &out); string_eval_teardown(&in, &out); } END_TEST START_TEST(test_invalid_command) { FILE *in, *out; char input[80], output[80]={}; string_eval_setup(&in, &out); sprintf(input, "print undefined_name__traceback_should_be_printed_to_stderr\n"); string_eval(in, out, input, 80, output); string_eval_teardown(&in, &out); } END_TEST START_TEST(test_simple_eval) { FILE *in, *out; char input[80], output[80]={}; string_eval_setup(&in, &out); sprintf(input, "print 3+4*5\n"); string_eval(in, out, input, 80, output); fail_unless(STRMATCH(output,"23\n"), NULL); string_eval_teardown(&in, &out); } END_TEST START_TEST(test_multiple_evals) { FILE *in, *out; char input[80], output[80]={}; string_eval_setup(&in, &out); sprintf(input, "print 3+4*5\n"); string_eval(in, out, input, 80, output); fail_unless(STRMATCH(output,"23\n"), NULL); sprintf(input, "print (3**2 + 4**2)**0.5\n"); string_eval(in, out, input, 80, output); fail_unless(STRMATCH(output,"5.0\n"), NULL); string_eval_teardown(&in, &out); } END_TEST START_TEST(test_eval_with_state) { FILE *in, *out; char input[80], output[80]={}; string_eval_setup(&in, &out); sprintf(input, "print 3+4*5\n"); fprintf(in, "A = 3\n"); sprintf(input, "print A*3\n"); string_eval(in, out, input, 80, output); fail_unless(STRMATCH(output,"9\n"), NULL); string_eval_teardown(&in, &out); } END_TEST @ <>= TCase *tc_string_eval = tcase_create("string_eval"); @ <>= tcase_add_test(tc_string_eval, test_setup_teardown); tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL); tcase_add_test(tc_string_eval, test_simple_eval); tcase_add_test(tc_string_eval, test_multiple_evals); tcase_add_test(tc_string_eval, test_eval_with_state); suite_add_tcase(s, tc_string_eval); @ \section{Accelerating function evaluation} My first version-0.3 code was running very slowly. With the options suggested in the help ([[-v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8]]), the program spent almost all of it's time in the functions [[oneD_solve]], [[wlc]], and [[wlc_handler]], making 0.66 million calls to [[oneD_solve]] and [[9.6]] million calls each to the WLC functions. That's only 15 calls per solution, so the search algorithm seems reasonable. The number of evaluation calls could be drastically reduced, however, by implementing an $x(x_0)$ lookup table. <>= #ifndef ACCEL_K_H #define ACCEL_K_H <> double accel_k(k_func_t *k, double F, environment_t *env, void *params); void free_accels(); #endif /* ACCEL_K_H */ @ <>= NW_SAWSIM_MODS += accel_k CHECK_BINS += check_accel_k check_accel_k_MODS = interp k_model parse string_eval tavl @ <>= <> #include /* assert() */ #include /* realloc(), free(), NULL */ #include "global.h" /* environment_t */ #include "k_model.h" /* k_func_t */ #include "interp.h" /* interp_* */ #include "accel_k.h" typedef struct accel_k_struct { interp_table_t *itable; k_func_t *k; environment_t *env; void *params; } accel_k_t; /* keep an array of all the ks we accelerate, so the caller doesn't have to worry */ static int num_accels = 0; static accel_k_t *accels=NULL; /* Wrap k in the standard f(x) acceleration form */ static double k_wrap(double F, void *params) { accel_k_t *p = (accel_k_t *)params; return (*p->k)(F, p->env, p->params); } static int k_tol(double FA, double kA, double FB, double kB) { assert(FB > FA); if (FB-FA > 1e-12) { //printf("unacceptable tol (x %g > 1 (y %g)\n", fabs(xA-xB), fabs(yA-yB)); return 1; /* unnacceptable */ } else { //printf("acceptable tol\n"); return 0; /* acceptable */ } } static int add_accel_k(k_func_t *k, environment_t *env, void *params) { int i=num_accels; accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels); assert(accels != NULL); accels[i].itable = interp_table_allocate(&k_wrap, &k_tol); accels[i].k = k; accels[i].env = env; accels[i].params = params; return i; } void free_accels() { int i; for (i=0; i>= <> <> <> <> <> <
> @ <>= #include /* EXIT_SUCCESS, EXIT_FAILURE */ #include /* sprintf() */ #include /* assert() */ #include /* exp() */ #include /* errno, ERANGE */ <> #include "global.h" #include "k_model.h" #include "accel_k.h" @ <>= <> <> @ <>= Suite *test_suite (void) { Suite *s = suite_create ("accelerated k model"); <> <> return s; } @ <>= TCase *tc_accel_k = tcase_create("accelerated k model"); @ <>= tcase_add_test(tc_accel_k, test_accel_k_accuracy); suite_add_tcase(s, tc_accel_k); @ <>= START_TEST(test_accel_k_accuracy) { double k, k_accel, F, Fm=0.0, FM=300e-12, dF=0.5e-12; k_func_t *k_func = bell_k; char *param_strings[] = {"3.3e-4", "0.25e-9"}; void *params = string_create_bell_param_t(param_strings); environment_t env = {}; env.T = 300.0; for(F=Fm; F <= FM; F+=dF) { k = k_func(F, &env, params); k_accel = accel_k(k_func, F, &env, params); CHECK_ERR(1e-8, k, k_accel); } destroy_bell_param_t(params); } END_TEST @ \section{Tension models} \label{sec.tension_models} TODO: tension model intro. 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. Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files. 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]]. <>= #ifndef TENSION_MODEL_H #define TENSION_MODEL_H <> <> <> <> <> <> <> <> <> #endif /* TENSION_MODEL_H */ @ <>= NW_SAWSIM_MODS += tension_model CHECK_BINS += check_tension_model check_tension_model_MODS = list tension_balance @ <>= TENSION_MODEL_MODS = tension_model parse list tension_balance TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \ $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) \ $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.h) $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR) notangle -Rtension-model-utils.c $< > $@ $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR) $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%) clean_tension_model_utils : rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils @ The pipe symbol [[|]] marks the prerequisits following it (in this case, the directories) as ‘order-only’ prerequisites. The timestamp on these prerequisits does not effect whether the rules are executed. This is appropriate for directories, where we don't need to recompile after an unrelated has been added to the directory, but only when the source prerequisites change. See the [[make]] documentation for more details (\url{http://www.gnu.org/software/make/manual/html_node/Prerequisite-Types.html}). \subsection{Null} \label{sec.null_tension} For unstretchable domains. <>= {NULL, NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL} @ \subsection{Constant} \label{sec.const_tension} An infinitely stretchable domain providing a constant tension. Obviously this cannot be inverted, so you can't put this domain in series with any other domains. We include it mostly for testing purposes. <>= <> <> @ <>= <> <> <> @ <>= extern double const_tension_handler(double x, void *pdata); @ <>= double const_tension_handler(double x, void *pdata) { tension_handler_data_t *data = (tension_handler_data_t *)pdata; list_t *list = data->group_tension_params; double F; assert(x >= 0.0); list = head(list); assert(list != NULL); /* empty active group?! */ F = ((const_tension_param_t *)list->d)->F; #ifdef DEBUG fprintf(stderr, "%s: tension %g N\n", __FUNCTION__, F); #endif while (list != NULL) { assert(((const_tension_param_t *)list->d)->F == F); list = list->next; } return F; } @ <>= typedef struct const_tension_param_struct { double F; /* tension (force) in N */ } const_tension_param_t; @ <>= extern void *string_create_const_tension_param_t(char **param_strings); extern void destroy_const_tension_param_t(void *p); @ <>= const_tension_param_t *create_const_tension_param_t(double F) { const_tension_param_t *ret = (const_tension_param_t *) malloc(sizeof(const_tension_param_t)); assert(ret != NULL); ret->F = F; return ret; } void *string_create_const_tension_param_t(char **param_strings) { return create_const_tension_param_t( safe_strtod(param_strings[0],__FUNCTION__)); } void destroy_const_tension_param_t(void *p) { if (p) free(p); } @ <>= extern int num_const_tension_params; extern char *const_tension_param_descriptions[]; extern char const_tension_param_string[]; @ <>= int num_const_tension_params = 1; char *const_tension_param_descriptions[] = {"tension F, N"}; char const_tension_param_string[] = "0"; @ <>= {&const_tension_handler, NULL, "const", "an infinitely stretchable domain with constant tension", num_const_tension_params, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t} @ \subsection{Hooke} \label{sec.hooke} <>= <> <> <> <> @ <>= <> <> <> @ The tension of a single spring is given by $F=kx$ for some spring constant $k$. The behavior of a series of springs $k_i$ in series is given by $$ F = \p({\sum_i \frac{1}{k_i}})^{-1} x $$ For a simple proof, see Appendix \ref{sec.math_hooke}. <>= extern double hooke_handler(double x, void *pdata); extern double inverse_hooke_handler(double F, void *pdata); @ First we define a function that computes the effective spring constant (stored in a single [[hooke_param_t]]) for the entire group. <>= static void hooke_param_grouper(tension_handler_data_t *data, hooke_param_t *param) { list_t *list = data->group_tension_params; double k=0.0; list = head(list); assert(list != NULL); /* empty active group?! */ while (list != NULL) { assert( ((hooke_param_t *)list->d)->k > 0 ); k += 1.0/ ((hooke_param_t *)list->d)->k; #ifdef DEBUG fprintf(stderr, "%s: Hooke spring %g N/m in series\n", __FUNCTION__, ((hooke_param_t *)list->d)->k); #endif list = list->next; } k = 1.0 / k; #ifdef DEBUG fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N (from data %p)\n", __FUNCTION__, k, x, k*x, data); #endif param->k = k; } @ Everyone knows Hooke's law: $F=kx$. <>= double hooke_handler(double x, void *pdata) { hooke_param_t param = {0}; assert(x >= 0.0); hooke_param_grouper((tension_handler_data_t *)pdata, ¶m); return param.k*x; } @ Inverting Hooke's law gives $x=F/k$. <>= double inverse_hooke_handler(double F, void *pdata) { hooke_param_t param = {0}; assert(F >= 0.0); hooke_param_grouper((tension_handler_data_t *)pdata, ¶m); return F/param.k; } @ Not much to keep track of for a single effective spring, just the spring constant $k$. <>= typedef struct hooke_param_struct { double k; /* spring constant in N/m */ } hooke_param_t; @ We wrap up our Hooke implementation with some book-keeping functions. <>= extern void *string_create_hooke_param_t(char **param_strings); extern void destroy_hooke_param_t(void *p); @ <>= hooke_param_t *create_hooke_param_t(double k) { hooke_param_t *ret = (hooke_param_t *) malloc(sizeof(hooke_param_t)); assert(ret != NULL); ret->k = k; return ret; } void *string_create_hooke_param_t(char **param_strings) { return create_hooke_param_t(safe_strtod(param_strings[0], __FUNCTION__)); } void destroy_hooke_param_t(void *p) { if (p) free(p); } @ <>= extern int num_hooke_params; extern char *hooke_param_descriptions[]; extern char hooke_param_string[]; @ <>= int num_hooke_params = 1; char *hooke_param_descriptions[] = {"spring constant k, N/m"}; char hooke_param_string[]="0.05"; @ <>= {hooke_handler, inverse_hooke_handler, "hooke", "a Hookean spring (F=kx)", num_hooke_params, hooke_param_descriptions, (char *)hooke_param_string, &string_create_hooke_param_t, &destroy_hooke_param_t} @ \subsection{Worm-like chain} \label{sec.wlc} We can model several unfolded domains with a single worm-like chain. <>= <> <> <> <> @ <>= <> <> <> <> @ <>= <> <> <> @ The combination of all unfolded domains is modeled as a worm like chain, with the force $F_{\text{WLC}}$ approximately given by $$ F_{\text{WLC}} \approx \frac{k_B T}{p} \p[{\frac{1}{4} \p({\frac{1}{(1-x/L)^2} - 1}) +\frac{x}{L}}], $$ where $x$ is the distance between the fixed ends, $k_B$ is Boltzmann's constant, $T$ is the absolute temperature, $p$ is the persistence length, and $L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}. <>= static double wlc(double x, double T, double p, double L) {/* N m K m m */ double xL=0.0; /* uses global k_B */ assert(x >= 0); assert(T > 0); assert(p > 0); assert(L > 0); if (x >= L) return HUGE_VAL; xL = x/L; /* unitless */ //printf("T %g\tp %g\tL %g\tx/L %g\tk %g\n", T, p, L, x/L, // k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL )); return k_B*T/p * ( 0.25*(1.0/((1-xL)*(1-xL))-1) + xL ); } @ Inverting the approximate WLC is fairly straightforward, if a bit tedious. Letting $x_L \equiv x/L$ and $F_T \equiv Fp/k_B T$ we have \begin{align} F_T &\approx \frac{1}{4} \p({\frac{1}{(1-x_L)^2} - 1}) + x_L \\ F_T (1-x_L)^2 &\approx \frac{1}{4} - \frac{1}{4} (1-x_L)^2 + x_L(1-x_L)^2 \\ 0 &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) (1-x_L)^2 + x_L(1-x_L)^2 \\ &\approx \frac{1}{4} - \p({F_T + \frac{1}{4}}) \p({1-2x_L+x_L^2}) + x_L - 2x_L^2 + x_L^3 \\ &\approx \p[{\frac{1}{4} - \p({F_T + \frac{1}{4}})}] + x_L \p[{2F_T + \frac{1}{2} + 1}] + x_L^2 \p[{-\p({F_T + \frac{1}{4}}) - 2}] + x_L^3 \\ &\approx -F_T + x_L \p({2F_T + \frac{3}{2}}) - x_L^2 \p({F_T + \frac{9}{4}}) + x_L^3 \\ &\approx ax_L^3 + bx_L^2 + cx_L + d \end{align} where $a \equiv 1$, $b \equiv -F_T - 9/4$, $c \equiv 2F_T + 3/2$, and $d \equiv -F_T$. % From \citet{wikipedia_cubic_function} the discriminant %\begin{align} % \Delta &= -4b^3d + b^2c^2 - 4ac^3 + 18abcd - 27a^2d^2 \\ % &= -4F_T\p({F_T + \frac{9}{4}})^3 % + \p({F_T + \frac{9}{4}})^2 \cdot \p({2F_T + \frac{3}{2}})^2 % - 4 \p({2F_T + \frac{3}{2}})^3 % + 18F_T \p({F_T + \frac{9}{4}}) \cdot \p({2F_T + \frac{3}{2}}) % - 27F_T^2 \\ % &= -4F_T \p({F_T^3 + \frac{27}{4}F_T^2 +\frac{243}{16}F_T +\frac{729}{64}}) %% a^3 + 3a^2b + 3ab^2 + b^3 % + \p({F_T^2 + \frac{9}{2}F_T + \frac{81}{16}}) % \cdot \p({4F_T^2 + 6F_T + \frac{9}{4}}) \\ % &\qquad - 4 \p({8F_T^3 + 18F_T^2 + \frac{27}{2}F_T + \frac{27}{8}}) %% a^3 + 3a^2b + 3ab^2 + b^3 % + 18F_T \p({2F_T^2 + 6F_T + \frac{27}{8}}) % - 27F_T^2 \\ % &= -4F_T^4 - 27F_T^3 - \frac{243}{4}F_T^2 - \frac{729}{16}F_T % + 4F_T^4 + 24F_T^3 +\frac{99}{2}F_T^2 +\frac{81}{2}F_T +\frac{729}{64}\\ % &\qquad - 32F_T^3 - 72F_T^2 - 54F_T - \frac{27}{2} % + 36F_T^3 + 108F_T^2 + \frac{243}{4}F_T % - 27F_T^2 \\ % &= F_T^4 \p({-4 + 4}) + F_T^3 \p({-27 + 24 - 32 + 36}) % + F_T^2 \p({-\frac{243}{4} + \frac{99}{2} - 72 + 108 -27}) \\ % &\qquad + F_T \p({-\frac{729}{16} + \frac{81}{2} - 54 + \frac{243}{4}}) % \p({\frac{729}{64} - \frac{27}{2}}) \\ % &= F_T^3 - \frac{9}{4}F_T^2 + \frac{27}{16}F_T - \frac{135}{64} %\end{align} We can plug these coefficients into the GSL cubic solver to invert the WLC\citep{galassi05}. The applicable root is always the one which comes before the singularity, which will be the smallest real root. <>= static double inverse_wlc(double F, double T, double p, double L) {/* m N K m m */ double FT = F*p/(k_B*T); /* uses global k_B */ double xL0, xL1, xL2; int num_roots; assert(F >= 0); assert(T > 0); assert(p > 0); assert(L > 0); if (F == HUGE_VAL) return L; num_roots = gsl_poly_solve_cubic(-(FT+2.25),2*FT+1.5,-FT, &xL0, &xL1, &xL2); assert(num_roots > 0); assert(xL0 >= -DOUBLE_PRECISION); assert(xL0 < 1); if (xL0 < 0) xL0 = 0.0; return xL0*L; } @ First we define a function that computes the effective WLC parameters (stored in a single [[wlc_param_t]]) for the entire group. <>= static void wlc_param_grouper(tension_handler_data_t *data, wlc_param_t *param) { list_t *list = data->group_tension_params; double p=0.0, L=0.0; list = head(list); assert(list != NULL); /* empty active group?! */ p = ((wlc_param_t *)list->d)->p; while (list != NULL) { /* enforce consistent p */ assert( ((wlc_param_t *)list->d)->p == p); L += ((wlc_param_t *)list->d)->L; #ifdef DEBUG fprintf(stderr, "%s: WLC section %g m long in series\n", __FUNCTION__, ((wlc_param_t *)list->d)->L); #endif list = list->next; } #ifdef DEBUG fprintf(stderr, "%s: WLC(x=%g m, T=%g K, p=%g m, L=%g m) = %g N\n", __FUNCTION__, x, data->env->T, p, L, wlc(x, data->env->T, p, L)); #endif param->p = p; param->L = L; } @ <>= extern double wlc_handler(double x, void *pdata); @ This model requires that each unfolded domain in the group have the same persistence length. <>= double wlc_handler(double x, void *pdata) { tension_handler_data_t *data = (tension_handler_data_t *)pdata; wlc_param_t param = {0}; assert(x >= 0.0); wlc_param_grouper(data, ¶m); return wlc(x, data->env->T, param.p, param.L); } @ <>= extern double inverse_wlc_handler(double F, void *pdata); @ <>= double inverse_wlc_handler(double F, void *pdata) { tension_handler_data_t *data = (tension_handler_data_t *)pdata; wlc_param_t param = {0}; wlc_param_grouper(data, ¶m); return inverse_wlc(F, data->env->T, param.p, param.L); } @ <>= typedef struct wlc_param_struct { double p; /* WLC persistence length (m) */ double L; /* WLC contour length (m) */ } wlc_param_t; @ <>= extern void *string_create_wlc_param_t(char **param_strings); extern void destroy_wlc_param_t(void *p /* wlc_param_t * */); @ <>= wlc_param_t *create_wlc_param_t(double p, double L) { wlc_param_t *ret = (wlc_param_t *) malloc(sizeof(wlc_param_t)); assert(ret != NULL); ret->p = p; ret->L = L; return ret; } void *string_create_wlc_param_t(char **param_strings) { return create_wlc_param_t(safe_strtod(param_strings[0], __FUNCTION__), safe_strtod(param_strings[1], __FUNCTION__)); } void destroy_wlc_param_t(void *p /* wlc_param_t * */) { if (p) free(p); } @ 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. TODO (now we copy the string before parsing, but still do this for clarity). For example, <>= char string[] = "abc"; string[1] = 'x'; @ works, but <>= char *string = "abc"; string[1] = 'x'; @ segfaults. <>= extern int num_wlc_params; extern char *wlc_param_descriptions[]; extern char wlc_param_string[]; @ <>= int num_wlc_params = 2; char *wlc_param_descriptions[] = {"persistence length p, m", "contour length L, m"}; char wlc_param_string[]="0.39e-9,28.4e-9"; @ <>= {&wlc_handler, &inverse_wlc_handler, "wlc", "a worm-like chain (F=kbT/p{.25[(1-x/L)^-2 -1] + x/L})", num_wlc_params, wlc_param_descriptions, (char *)wlc_param_string, &string_create_wlc_param_t, &destroy_wlc_param_t} @ Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696. \subsection{Freely-jointed chain} \label{sec.fjc} <>= <> <> <> @ <>= <> <> <> @ The freely jointed chain model approximates the protein as a series of rigid links of length $l$ with freely rotating joints. The tension of a chain of $N$ such links is given by $$ F(x) = \frac{k_B T}{l} \mathcal{L}^{-1}\p({\frac{x}{L}}) \;, $$ 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}. <>= static double langevin(double x, void *params) { if (x==0) return 0; return 1.0/tanh(x) - 1/x; } static double inv_langevin(double x) { double lb=0.0, ub=1.0, ret; oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub); ret = oneD_solve(&langevin, NULL, x, lb, ub, 1e-6, 1e-6, 300, NULL); /* TODO, break out oneD_* to own file */ return ret; } static double fjc(double x, double T, double l, int N) { double L = l*(double)N; if (x == 0) return 0; //printf("x %g\tL%g\tx/L %g\n", x, L, x/L); return k_B*T/l * inv_langevin(x/L); } @ <>= extern double fjc_handler(double x, void *pdata); @ <>= double fjc_handler(double x, void *pdata) { tension_handler_data_t *data = (tension_handler_data_t *)pdata; list_t *list = data->group_tension_params; double l; int N=0; list = head(list); assert(list != NULL); /* empty active group?! */ l = ((fjc_param_t *)list->d)->l; while (list != NULL) { /* enforce consistent l */ assert( ((fjc_param_t *)list->d)->l == l); N += ((fjc_param_t *)list->d)->N; #ifdef DEBUG fprintf(stderr, "%s: FJC section %d links long in series\n", __FUNCTION__, ((fjc_param_t *)list->d)->N); #endif list = list->next; } #ifdef DEBUG fprintf(stderr, "%s: FJC(x=%g m, T=%g K, l=%g m, N=%d m) = %g N\n", __FUNCTION__, x, data->env->T, l, N, fjc(x, data->env->T, l, N)); #endif return fjc(x, data->env->T, l, N); } @ <>= typedef struct fjc_param_struct { double l; /* FJC link length (m) */ int N; /* FJC number of links */ } fjc_param_t; @ <>= extern void *string_create_fjc_param_t(char **param_strings); extern void destroy_fjc_param_t(void *p); @ <>= fjc_param_t *create_fjc_param_t(double l, double N) { fjc_param_t *ret = (fjc_param_t *) malloc(sizeof(fjc_param_t)); assert(ret != NULL); ret->l = l; ret->N = N; return ret; } void *string_create_fjc_param_t(char **param_strings) { return create_fjc_param_t(safe_strtod(param_strings[0], __FUNCTION__), safe_strtod(param_strings[1], __FUNCTION__)); } void destroy_fjc_param_t(void *p) { if (p) free(p); } @ <>= extern int num_fjc_params; extern char *fjc_param_descriptions[]; extern char fjc_param_string[]; @ <>= int num_fjc_params = 2; char *fjc_param_descriptions[] = {"link length l, m", "link count N, unitless integer"}; char fjc_param_string[]="0.5e-9,200"; @ <>= {fjc_handler, NULL, "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} @ \subsection{Piston} \label{sec.piston_tension} An infinitely stretchable domain with no tension for extensions less than a particular threshold $L$, and infinite tension for $x>L$. The tension at $x=L$ is undefined, but may be any positive value. The model is called the ``piston'' model because it resembles a frictionless piston sliding in a rigid cylinder of length $L$. $$ F = \begin{cases} 0 & \text{if $xL$}. \end{cases} $$ Note that because of it's infinte stiffness at $x=L$, fully extended domains with this tension model will be identical to completely rigid domains. However, there is a distinction between this model and rigid domains for $x>= <> <> <> <> @ <>= <> <> <> <> @ First we define a function that computes the effective piston parameters (stored in a single [[piston_tension_param_t]]) for the entire group. <>= static void piston_tension_param_grouper(tension_handler_data_t *data, piston_tension_param_t *param) { list_t *list = data->group_tension_params; double L=0; list = head(list); assert(list != NULL); /* empty active group?! */ while (list != NULL) { L += ((piston_tension_param_t *)list->d)->L; list = list->next; } param->L = L; } <>= extern double piston_tension_handler(double x, void *pdata); @ <>= double piston_tension_handler(double x, void *pdata) { tension_handler_data_t *data = (tension_handler_data_t *)pdata; piston_tension_param_t param = {0}; double F; piston_tension_param_grouper(data, ¶m); if (x <= param.L) F = 0; else F = HUGE_VAL; #ifdef DEBUG fprintf(stderr, "%s: x %g, L %g, tension %g N\n", __FUNCTION__, x, param.L, F); #endif return F; } @ We abuse [[x_of_xo()]] and [[oneD_solve()]] with this inverse definition (see Section \ref{sec.tension_balance}), because the $x(F=0)<=L$, but our function returns $x(F=0)=0$. This is fine when the total extension \emph{is} zero, but cheats a bit when there is some total extension. For any non-zero extension, the initial active group produces some tension (we hope. True for all our other tension models). This tension fully extends the piston group (of which there should be only one, since all piston states can get lumped into the same tension group.). If the total extension is $\le$ the target extension, the full extension is accurate, so the inverse was valid after all. If not, [[oneD_solve()]] will halve the extension of the first group, reducing the overall tension. For total extension less than the piston extension, this bisection continues for [[max_steps]], leaving $x_0$ reduced by a factor of $2^\text{max\_steps}$, a very small number. Because of this, the returned tension $F_0(x_0)$ is very small, even though the total extension found by [[x_of_xo()]] is still $>$ the piston length. While this works, it is clearly not the most elegant or robust solution. TODO: think of (and implement) a better procedure. <>= extern double inverse_piston_tension_handler(double x, void *pdata); @ <>= double inverse_piston_tension_handler(double F, void *pdata) { tension_handler_data_t *data = (tension_handler_data_t *)pdata; piston_tension_param_t param = {0}; assert(F >= 0.0); piston_tension_param_grouper(data, ¶m); if (F == 0.0) return 0; return param.L; } @ <>= typedef struct piston_tension_param_struct { double L; /* length of piston in m */ } piston_tension_param_t; @ <>= extern void *string_create_piston_tension_param_t(char **param_strings); extern void destroy_piston_tension_param_t(void *p); @ <>= piston_tension_param_t *create_piston_tension_param_t(double L) { piston_tension_param_t *ret = (piston_tension_param_t *) malloc(sizeof(piston_tension_param_t)); assert(ret != NULL); ret->L = L; return ret; } void *string_create_piston_tension_param_t(char **param_strings) { return create_piston_tension_param_t( safe_strtod(param_strings[0],__FUNCTION__)); } void destroy_piston_tension_param_t(void *p) { if (p) free(p); } @ <>= extern int num_piston_tension_params; extern char *piston_tension_param_descriptions[]; extern char piston_tension_param_string[]; @ <>= int num_piston_tension_params = 1; char *piston_tension_param_descriptions[] = {"length L, m"}; char piston_tension_param_string[] = "0"; @ <>= {&piston_tension_handler, &inverse_piston_tension_handler, "piston", "a domain that slides without tension for xL.", num_piston_tension_params, piston_tension_param_descriptions, piston_tension_param_string, &string_create_piston_tension_param_t, &destroy_piston_tension_param_t} @ \subsection{Tension model implementation} <>= //#define DEBUG <> <> #include "tension_model.h" <> <> <> @ <>= #include /* assert() */ #include /* NULL, strto*() */ #include /* HUGE_VAL, sqrt(), exp() */ #include /* fprintf(), stdout */ #include /* errno, ERANGE */ #include "global.h" #include "list.h" #include "tension_balance.h" /* oneD_*() */ @ <>= <> <> <> <> <> @ <>= <> <> <> <> <> <> @ <>= <> <> <> <> <> <> @ \subsection{Utilities} The tension models can get complicated, and users may want to reassure themselves that this portion of the input physics are functioning properly. The stand-alone program [[tension_model_utils]] demonstrates the tension model interface without the overhead of having to understand a full simulation. [[tension_model_utils]] takes command line model arguments like the full simulation, and prints $F(x)$ data to the screen for the selected model over a range of $x$. <>= <> <> <> <> <> <> @ <>= int main(int argc, char **argv) { <> tension_model_getopt_t *model = NULL; void *params; environment_t env; unsigned int flags; one_dim_fn_t *tension_handler[NUM_TENSION_MODELS] = {0}; tension_handler_data_t tdata; int num_param_args; /* for INIT_MODEL() */ char **param_args; /* for INIT_MODEL() */ double Fmax,Xmax; get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model, &Fmax, &Xmax, &flags); setup(); if (flags & VFLAG) { printf("#initializing model %s with parameters %s\n", model->name, model->params); } INIT_MODEL("utility", model, model->params, params); tdata.group_tension_params = NULL; push(&tdata.group_tension_params, params, 1); tdata.env = &env; tdata.persist = NULL; if (model->handler == NULL) { printf("No tension function for model %s\n", model->name); exit(0); } { int i,N=200; double x=0, F=0; printf("#Distance (m)\tForce (N)\n"); for (i=0; i<=N; i++) { x = Xmax*i/(double)N; F = (*model->handler)(x, &tdata); if (F < 0 || F > Fmax) break; printf("%g\t%g\n", x, F); } if (flags & VFLAG && i <= N) { /* explain exit condition */ if (F < 0) printf("Impossible force %g\n", F); else if (F > Fmax) printf("Reached large force limit %g > %g\n", F, Fmax); } } params = pop(&tdata.group_tension_params); if (model->destructor) (*model->destructor)(params); return 0; } @ <>= #include #include #include /* strlen() */ #include /* assert() */ #include /* errno, ERANGE */ #include "global.h" #include "parse.h" #include "list.h" #include "tension_model.h" @ <>= <> #define VFLAG 1 /* verbose */ <> <> <> @ <>= <> <> <> <> @ <>= void help(char *prog_name, environment_t *env, int n_tension_models, tension_model_getopt_t *tension_models, int tension_model, double Fmax, double Xmax) { int i, j; printf("usage: %s [options]\n", prog_name); printf("Version %s\n\n", VERSION); printf("A utility for understanding the available tension models\n\n"); printf("Environmental options:\n"); printf("-T T\tTemperature T (currently %g K)\n", env->T); printf("-C T\tYou can also set the temperature T in Celsius\n"); printf("Model options:\n"); printf("-m model\tFolded tension model (currently %s)\n", tension_models[tension_model].name); printf("-a args\tFolded tension model argument string (currently %s)\n", tension_models[tension_model].params); printf("F(x) is calculated for a range of x and printed\n"); printf("For example:\n"); printf(" #Distance (m)\tForce (N)\n"); printf(" 123.456\t7.89\n"); printf(" ...\n"); printf("-F\tSet the maximum F value for the standard mode F(x) output (currently %g)\n", Fmax); printf("-X\tSet the maximum x value for the standard mode F(x) output (currently %g)\n", Xmax); printf("-V\tChange output to verbose mode\n"); printf("-h\tPrint this help and exit\n"); printf("\n"); printf("Tension models:\n"); for (i=0; i>= void get_options(int argc, char **argv, environment_t *env, int n_tension_models, tension_model_getopt_t *tension_models, tension_model_getopt_t **model, double *Fmax, double *Xmax, unsigned int *flags) { char *prog_name = NULL; char c, options[] = "T:C:m:a:F:X:Vh"; int tension_model=0, num_strings; extern char *optarg; extern int optind, optopt, opterr; assert(argc > 0); /* setup defaults */ prog_name = argv[0]; env->T = 300.0; /* K */ *Fmax = 1e5; /* N */ *Xmax = 1e-6; /* m */ *flags = 0; *model = tension_models; while ((c=getopt(argc, argv, options)) != -1) { switch(c) { case 'T': env->T = safe_strtod(optarg, "-T"); break; case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break; case 'm': tension_model = index_tension_model(n_tension_models, tension_models, optarg); *model = tension_models+tension_model; break; case 'a': tension_models[tension_model].params = optarg; break; case 'F': *Fmax = safe_strtod(optarg, "-F"); break; case 'X': *Xmax = safe_strtod(optarg, "-X"); break; case 'V': *flags |= VFLAG; break; case '?': fprintf(stderr, "unrecognized option '%c'\n", optopt); /* fall through to default case */ default: help(prog_name, env, n_tension_models, tension_models, tension_model, *Fmax, *Xmax); exit(1); } } return; } @ \subsection{Tension model unit tests} Here we check to make sure the various functions work as expected, using \citetalias{sw:check}. <>= <> <> <> <> <> <
> @ <>= #include /* EXIT_SUCCESS, EXIT_FAILURE */ #include /* sprintf() */ #include /* assert() */ #include /* exp() */ #include /* errno, ERANGE */ <> #include "global.h" #include "list.h" #include "tension_balance.h" /* oneD_*() */ #include "tension_model.h" @ <>= <> <> <> <> <> <> <> @ <>= Suite *test_suite (void) { Suite *s = suite_create ("tension model"); <> <> <> <> <> <> <> <> <> <> return s; } @ \subsubsection{Constant} <>= TCase *tc_const_tension = tcase_create("Constant tension"); @ <>= tcase_add_test(tc_const_tension, test_const_tension_create_destroy); suite_add_tcase(s, tc_const_tension); @ <>= START_TEST(test_const_tension_create_destroy) { char *tension[] = {"1","2.2", "3"}; char *params[] = {tension[0]}; void *p = NULL; int i; for( i=0; i < sizeof(tension)/sizeof(char *); i++) { params[0] = tension[i]; p = string_create_const_tension_param_t(params); destroy_const_tension_param_t(p); } } END_TEST @ TODO: In order to test the constant tension handler itself, we'd have to construct a group. \subsubsection{Hooke} <>= TCase *tc_hooke = tcase_create("Hooke"); @ <>= tcase_add_test(tc_hooke, test_hooke_create_destroy); suite_add_tcase(s, tc_hooke); @ <>= START_TEST(test_hooke_create_destroy) { char *k[] = {"1","2.2", "3"}; char *params[] = {k[0]}; void *p = NULL; int i; for( i=0; i < sizeof(k)/sizeof(char *); i++) { params[0] = k[i]; p = string_create_hooke_param_t(params); destroy_hooke_param_t(p); } } END_TEST @ TODO: In order to test the Hooke tension handler itself, we'd have to construct a group. \subsubsection{Worm-like chain} <>= TCase *tc_wlc = tcase_create("WLC"); @ <>= tcase_add_test(tc_wlc, test_wlc_at_zero); tcase_add_test(tc_wlc, test_wlc_at_half); suite_add_tcase(s, tc_wlc); @ <>= <> START_TEST(test_wlc_at_zero) { double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30; fail_unless(abs(wlc(x, T, p, L)) < lim, \ "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g", x, T, p, L, abs(wlc(x, T, p, L)), lim); } END_TEST START_TEST(test_wlc_at_half) { double T=1.0, L=1.0, p=0.1*k_B, x=0.5; /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1) * = 0.25 * 3 = 3/4 * linear term = x/L = 0.5 * nonlinear + linear = 0.75 + 0.5 = 1.25 * wlc = 10e21*1.25 = 12.5e21 */ fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16, "wlc(%g, %g, %g, %g) = %g != %g", x, T, p, L, wlc(x, T, p, L), 12.5e21); } END_TEST @ \subsubsection{Freely-jointed chain} <>= TCase *tc_fjc = tcase_create("FJC"); @ <>= tcase_add_test(tc_fjc, test_fjc_at_zero); tcase_add_test(tc_fjc, test_fjc_at_half); suite_add_tcase(s, tc_fjc); @ <>= <> START_TEST(test_fjc_at_zero) { int N=10; double T=1.0, l=1.0, p=0.1, x=0.0, lim=1e-30; fail_unless(abs(fjc(x, T, l, N)) < lim, \ "abs(fjc(x=%g,T=%g,l=%g,N=%d)) = %g >= %g", x, T, l, N, abs(fjc(x, T, l, N)), lim); } END_TEST START_TEST(test_fjc_at_half) { int N = 10; double T=1.0/k_B, l=1.0, x=5.0; /* prefactor = k_B T / l = k_B (1.0/k_B) / 1.0 = 1.0 J/nm = 1.0e21 pN * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1) * = 0.25 * 3 = 3/4 * linear term = x/L = 0.5 * nonlinear + linear = 0.75 + 0.5 = 1.25 * fjc = 10e21*1.25 = 12.5e21 */ fail_unless(fjc(x, T, l, N)-12.5e21 < 1e16, "fjc(%g, %g, %g, %d) = %g != %g", x, T, l, N, fjc(x, T, l, N), 12.5e21); } END_TEST @ \subsubsection{Piston} <>= TCase *tc_piston = tcase_create("Piston"); @ <>= tcase_add_test(tc_piston, test_piston_create_destroy); suite_add_tcase(s, tc_piston); @ <>= START_TEST(test_piston_create_destroy) { char *L[] = {"1","2.2", "3"}; char *params[] = {L[0]}; void *p = NULL; int i; for( i=0; i < sizeof(L)/sizeof(char *); i++) { params[0] = L[i]; p = string_create_piston_tension_param_t(params); destroy_piston_tension_param_t(p); } } END_TEST @ TODO: In order to test the piston tension handler itself, we'd have to construct a group. \section{Unfolding rate models} \label{sec.k_models} We have two main choices for determining $k$: Bell's model, which assumes the folded domains exist in a single `folded' state and make instantaneous, stochastic transitions over a free energy barrier; and Kramers' model, which treats the unfolding as a thermalized diffusion process. We deal with these options like we dealt with the different tension models: by bundling all model-specific parameters into structures, with model specific functions for determining $k$. <>= typedef double k_func_t(double F, environment_t *env, void *params); @ Because these models are independent of the rest of the simulation, we'll split their definitions out into their own set of files. 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]]. We go against [[noweb]]'s suggested naming scheme of $\langle \text{\emph{filename}} \rangle$, because \LaTeX\ doesn't like underscores outside math mode. Instead we replace the underscores (`\verb+_+') with hyphens (`\verb+-+'). You could use spaces instead (`\verb+ +'), but then your [[notangle]]s would look like [[notangle -R'k model.h']], which I don't like as much. <>= #ifndef K_MODEL_H #define K_MODEL_H <> <> <> <> <> <> <> <> #endif /* K_MODEL_H */ @ <>= NW_SAWSIM_MODS += k_model CHECK_BINS += check_k_model check_k_model_MODS = parse string_eval @ [[check_k_model]] also depends on the parse module. <>= K_MODEL_MODS = k_model parse string_eval K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \ $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(K_MODEL_MODS:%=$(BUILD_DIR)/%.h) $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR) notangle -Rk-model-utils.c $< > $@ $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR) $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%) clean_k_model_utils : rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils @ \subsection{Null} \label{sec.null_k} For domains that are never folded. <>= @ <>= @ <>= @ <>= {"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL} @ \subsection{Constant rate model} \label{sec.k_const} This is a pretty boring model, but it is usefull for testing the engine. $$ k = k_0\;, $$ where $k_0$ is the constant unfolding rate. <>= <> <> @ <>= <> <> <> @ <>= typedef struct const_k_param_struct { double knot; /* transition rate k_0 (frac population per s) */ } const_k_param_t; @ <>= double const_k(double F, environment_t *env, void *const_k_params); @ <>= double const_k(double F, environment_t *env, void *const_k_params) { /* Returns the rate constant k in frac pop per s. */ const_k_param_t *p = (const_k_param_t *) const_k_params; assert(p != NULL); assert(p->knot > 0); return p->knot; } @ <>= void *string_create_const_k_param_t(char **param_strings); void destroy_const_k_param_t(void *p); @ <>= const_k_param_t *create_const_k_param_t(double knot) { const_k_param_t *ret = (const_k_param_t *) malloc(sizeof(const_k_param_t)); assert(ret != NULL); ret->knot = knot; return ret; } void *string_create_const_k_param_t(char **param_strings) { return create_const_k_param_t(safe_strtod(param_strings[0], __FUNCTION__)); } void destroy_const_k_param_t(void *p /* const_k_param_t * */) { if (p) free(p); } @ <>= extern int num_const_k_params; extern char *const_k_param_descriptions[]; extern char const_k_param_string[]; @ <>= int num_const_k_params = 1; char *const_k_param_descriptions[] = {"unforced unfolding rate knot, % pop per s"}; char const_k_param_string[]="1"; @ <>= {"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} @ \subsection{Bell's model} \label{sec.bell} Quantitatively, Bell's model gives $k$ as $$ k = k_0 \cdot e^{F dx / k_B T} \;, $$ where $k_0$ is the unforced unfolding rate, $F$ is the applied unfolding force, $dx$ is the distance to the transition state, and $k_B T$ is the thermal energy\citep{schlierf06}. <>= <> <> @ <>= <> <> <> @ <>= typedef struct bell_param_struct { double knot; /* transition rate k_0 (frac population per s) */ double dx; /* `distance to transition state' (nm) */ } bell_param_t; @ <>= double bell_k(double F, environment_t *env, void *bell_params); @ <>= double bell_k(double F, environment_t *env, void *bell_params) { /* Returns the rate constant k in frac pop per s. * Takes F in N, T in K, knot in frac pop per s, and dx in m. * uses global k_B in J/K */ bell_param_t *p = (bell_param_t *) bell_params; assert(F >= 0); assert(env->T > 0); assert(p != NULL); assert(p->knot > 0); //assert(p->dx >= 0); /* p < 0 for refolding reactions, etc. */ /* printf("F %g\tT %g\tknot %g\tdx %g\tFdx/kbT %g\tk %g\n", F, env->T, p->knot, p->dx, F*p->dx/(k_B*env->T), p->knot * exp(F*p->dx / (k_B*env->T) )); printf("Fdx %g, kbT %g\n", F*p->dx, k_B*env->T); */ return p->knot * exp(F*p->dx / (k_B*env->T) ); } @ <>= void *string_create_bell_param_t(char **param_strings); void destroy_bell_param_t(void *p); @ <>= bell_param_t *create_bell_param_t(double knot, double dx) { bell_param_t *ret = (bell_param_t *) malloc(sizeof(bell_param_t)); assert(ret != NULL); ret->knot = knot; ret->dx = dx; return ret; } void *string_create_bell_param_t(char **param_strings) { return create_bell_param_t(safe_strtod(param_strings[0], __FUNCTION__), safe_strtod(param_strings[1], __FUNCTION__)); } void destroy_bell_param_t(void *p /* bell_param_t * */) { if (p) free(p); } @ <>= extern int num_bell_params; extern char *bell_param_descriptions[]; extern char bell_param_string[]; @ <>= int num_bell_params = 2; char *bell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"}; char bell_param_string[]="3.3e-4,0.25e-9"; @ <>= {"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} @ Initialized with Titin I27 parameters\citep{carrion-vazquez99a}. % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696 \subsection{Linker stiffness corrected Bell model} \label{sec.kbell} Walton et. al showed that the Bell model constant force approximation doesn't hold for biotin-streptavadin binding, and I extended their results to I27 for the 2009 Biophysical Society Annual Meeting\cite{walton08,king09}. More details to follow when I get done with the conference\ldots We adjust Bell's model to $$ k = k_0 \cdot e^{\frac{F dx - \frac{1}{2} \kappa dx^2}{k_B T}} \;, $$ where $k_0$ is the unforced unfolding rate, $F$ is the applied unfolding force, $\kappa$ is the effective spring constant, $dx$ is the distance to the transition state, and $k_B T$ is the thermal energy\citep{schlierf06}. <>= <> <> @ <>= <> <> <> @ <>= typedef struct kbell_param_struct { double knot; /* transition rate k_0 (frac population per s) */ double dx; /* `distance to transition state' (nm) */ } kbell_param_t; @ <>= double kbell_k(double F, environment_t *env, void *kbell_params); @ <>= double kbell_k(double F, environment_t *env, void *kbell_params) { /* Returns the rate constant k in frac pop per s. * Takes F in N, T in K, knot in frac pop per s, and dx in m. * uses global k_B in J/K */ kbell_param_t *p = (kbell_param_t *) kbell_params; assert(F >= 0); assert(env->T > 0); assert(env->stiffness >= 0); assert(p != NULL); assert(p->knot > 0); assert(p->dx >= 0); return p->knot * exp((F*p->dx - 0.5*env->stiffness*(p->dx)*(p->dx))/(k_B*env->T) ); } @ <>= void *string_create_kbell_param_t(char **param_strings); void destroy_kbell_param_t(void *p); @ <>= kbell_param_t *create_kbell_param_t(double knot, double dx) { kbell_param_t *ret = (kbell_param_t *) malloc(sizeof(kbell_param_t)); assert(ret != NULL); ret->knot = knot; ret->dx = dx; return ret; } void *string_create_kbell_param_t(char **param_strings) { return create_kbell_param_t(safe_strtod(param_strings[0], __FUNCTION__), safe_strtod(param_strings[1], __FUNCTION__)); } void destroy_kbell_param_t(void *p /* kbell_param_t * */) { if (p) free(p); } @ <>= extern int num_kbell_params; extern char *kbell_param_descriptions[]; extern char kbell_param_string[]; @ <>= int num_kbell_params = 2; char *kbell_param_descriptions[] = {"unforced unfolding rate knot, % pop per s", "distance to transition state, m"}; char kbell_param_string[]="3.3e-4,0.25e-9"; @ <>= {"kbell", "thermalized, two-state transitions, corrected for effective linker stiffness", &kbell_k, num_kbell_params, kbell_param_descriptions, (char *)kbell_param_string, &string_create_kbell_param_t, &destroy_kbell_param_t} @ Initialized with Titin I27 parameters\citep{carrion-vazquez99a}. % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696 \subsection{Kramer's model} \label{sec.kramers} Kramer's model gives $k$ as %$$ % \frac{1}{k} = \frac{1}{D} % \int_{x_\text{min}}^{x_\text{max}} % e^{\frac{-U_F(x)}{k_B T}} % \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx' % dx, %$$ %where $D$ is the diffusion constant, %$U_F(x)$ is the free energy along the unfolding coordinate, and %$x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$ % before the minimum and after the maximum, respectively \citep{schlierf06}. $$ k \approx \frac{D}{l_c l_{ts}} e^{-E_b/k_B T}, $$ where $D$ is the diffusion constant, $$ l_a = \sqrt{\frac{2\pi k_B T}{(\partial^2 E(x') / \partial x'^2)_{x_a}}}, $$ are length scales where $x_c(F)$ is the location of the bound state and $x_{ts}(F)$ is the location of the transition state, $E(F,x)$ is the free energy, and $E_b(F) = E(F, x_{ts}(F)) - E(F, x_c(F))$ is the barrier energy (\citet{vanKampen07} Eqn. XIII.2.2, \citet{hanggi90} Eqn.~4.56c, \citet{evans97} Eqn.~3). <>= <> <> @ <>= <> <> <> @ <>= //typedef double kramers_x_func_t(double F, double T, double *E_params, int n); //typedef double kramers_E_func_t(double F, double T, double *E_params, int n, double x); typedef struct kramers_param_struct { double D; /* diffusion rate (um/s) */ char *xc; char *xts; char *ddEdxx; char *E; FILE *in; FILE *out; //kramers_x_func_t *xc; /* function returning metastable x (nm) */ //kramers_x_func_t *xts; /* function returning transition x (nm) */ //kramers_E_func_t *ddEdxx; /* function returning d^2E/dx^2 (J/nm^2) */ //kramers_E_func_t *E; /* function returning E (J) */ //double *E_params; /* parametrize the energy landscape */ //int n_E_params; /* length of E_params array */ } kramers_param_t; @ <>= extern double kramers_E(double F, environment_t *env, void *kramers_params, double x); /* for k_model_utils */ extern double kramers_k(double F, environment_t *env, void *kramers_params); @ <>= double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T) { static char input[160]={0}; static char output[80]={0}; /* setup the environment */ fprintf(in, "F = %g; T = %g\n", F, T); sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */ string_eval(in, out, input, 80, output); return safe_strtod(output, __FUNCTION__); } double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x) { static char input[160]={0}; static char output[80]={0}; /* setup the environment */ fprintf(in, "F = %g;T = %g;x = %g\n", F, T, x); sprintf(input, "print repr(%s)\n", fn); /* repr() to keep full precision */ string_eval(in, out, input, 80, output); return safe_strtod(output, __FUNCTION__); } double kramers_E(double F, environment_t *env, void *kramers_params, double x) { kramers_param_t *p = (kramers_param_t *) kramers_params; return kramers_E_fn(p->in, p->out, p->E, F, env->T, x); } double kramers_k(double F, environment_t *env, void *kramers_params) { /* Returns the rate constant k in frac pop per s. * Takes F in N, T in K, knot in frac pop per s, and dx in m. * uses global k_B in J/K */ kramers_param_t *p = (kramers_param_t *) kramers_params; double kbT, xc, xts, lc, lts, Eb; assert(F >= 0); assert(env->T > 0); assert(p != NULL); assert(p->D > 0); assert(p->xc != NULL); assert(p->xts != NULL); assert(p->ddEdxx != NULL); assert(p->E != NULL); assert(p->in != NULL); assert(p->out != NULL); kbT = k_B*env->T; xc = kramers_x_fn(p->in, p->out, p->xc, F, env->T); if (xc == -1.0) return -1.0; /* error (Too much force) */ xts = kramers_x_fn(p->in, p->out, p->xts, F, env->T); if (xc == -1.0) return -1.0; /* error (Too much force) */ lc = sqrt(2.0*M_PI*kbT / kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xc)); lts = sqrt(-2.0*M_PI*kbT/ kramers_E_fn(p->in, p->out, p->ddEdxx, F, env->T, xts)); Eb = kramers_E_fn(p->in, p->out, p->E, F, env->T, xts) - kramers_E_fn(p->in, p->out, p->E, F, env->T, xc); //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)); return p->D/(lc*lts) * exp(-Eb/kbT); } @ <>= void *string_create_kramers_param_t(char **param_strings); void destroy_kramers_param_t(void *p); @ <>= kramers_param_t *create_kramers_param_t(double D, char *xc_fn, char *xts_fn, char *E_fn, char *ddEdxx_fn, char *global_define_string) // kramers_x_func_t *xc, kramers_x_func_t *xts, // kramers_E_func_t *ddEdxx, kramers_E_func_t *E, // double *E_params, int n_E_params) { kramers_param_t *ret = (kramers_param_t *) malloc(sizeof(kramers_param_t)); assert(ret != NULL); ret->xc = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1)); assert(ret->xc != NULL); ret->xts = (char *) malloc(sizeof(char)*(strlen(xc_fn)+1)); assert(ret->xts != NULL); ret->E = (char *) malloc(sizeof(char)*(strlen(E_fn)+1)); assert(ret->E != NULL); ret->ddEdxx = (char *) malloc(sizeof(char)*(strlen(ddEdxx_fn)+1)); assert(ret->ddEdxx != NULL); ret->D = D; strcpy(ret->xc, xc_fn); strcpy(ret->xts, xts_fn); strcpy(ret->E, E_fn); strcpy(ret->ddEdxx, ddEdxx_fn); //ret->E_params = E_params; //ret->n_E_params = n_E_params; string_eval_setup(&ret->in, &ret->out); fprintf(ret->in, "kB = %g\n", k_B); fprintf(ret->in, "%s\n", global_define_string); return ret; } /* takes an array of 4 functions: {"(1,2),(2,0),..."} */ void *string_create_kramers_param_t(char **param_strings) { return create_kramers_param_t(safe_strtod(param_strings[0], __FUNCTION__), param_strings[2], param_strings[3], param_strings[4], param_strings[5], param_strings[1]); } void destroy_kramers_param_t(void *p /* kramers_param_t * */) { kramers_param_t *pk = (kramers_param_t *)p; if (pk) { string_eval_teardown(&pk->in, &pk->out); //if (pk->E_params) // free(pk->E_params); free(pk); } } @ For our defaults, we'll follow \citet{schlierf06} and study ddFLN4. 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. We follow \citet{shillcock98} and use $$ E(F,E_b,x,x_s) \equiv \frac{27}{4} E_b \p({\frac{x}{x_s}})^2 \p({2-3\frac{x}{x_s}}) -F x $$ where TODO, variable meanings.%\citep{shillcock98}. The first and second derivatives of $E(F,E_b,x,x_s)$ are \begin{align} 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\\ E''(F,E_b,x,x_s)&=\frac{27 E_b}{2 x_s^2}\p({2-9\frac{x}{x_s}}). \end{align} $E' = ax^2 + bx + c$, with $a=\frac{-243 E_b}{4 x_s^3}$, $b=\frac{27 E_b}{x_s^2}$, and $c=-F$. The roots are therefor at \begin{align} x_\pm &= -\frac{b}{2a}\pm\sqrt{\p({\frac{b}{2a}})^2 - \frac{c}{a}} \\ &= \frac{2 x_s}{9}\pm\sqrt{\p({\frac{2 x_s}{9}})^2-\frac{4Fx_s^3}{243E_b}}\\ &= \frac{2 x_s}{9} \p({1 \pm \sqrt{1-\frac{F x_s}{3 E_b}}}) \end{align} <>= extern int num_kramers_params; extern char *kramers_param_descriptions[]; extern char kramers_param_string[]; @ <>= int num_kramers_params = 6; 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"}; 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)}"; @ Which we initialized with parameters for ddFLN4\citep{schlierf06}. There is a bit of funny buisness in the $x_c$ definition to handle high force cases. When $F x_s > 3 E_b$, the root will be complex (negative energy barrier), and the saddle point analysis breaks down. 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{sec.adaptive_dt} reduces it's timestep appropriately. The [[(cond) and (val_1) or (val_2)]] construction is Python hackery to simulate [[C]]'s ternary operator [[(cond) ? (val_1) : (val_2)]]. It works so long as [[val_1]] is not [[false]]. <>= {"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} @ \subsection{Kramer's model (natural cubic splines)} \label{sec.kramers_integ} Alternatively, we could skip the gaussian approximation and evaluate the double-integral form of Kramers' overdamped model, (\citet{hanggi90} Eqn.~4.56b, \citet{vanKampen07} Eqn.~XIII.2.1, \citet{schlierf06} Eqn.~2) $$ \frac{1}{k} = \frac{1}{D} \int_{x_\text{min}}^{x_\text{max}} e^{\frac{U_F(x)}{k_B T}} \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx' dx, $$ where $D$ is the diffusion constant, $U_F(x)$ is the free energy along the unfolding coordinate, and $x_\text{min}$ and $x_\text{max}$ are arbitrarily chosen to be $6 k_B T$ before the minimum and after the maximum, respectively (note: there is a sign difference from Schlierf)\citep{schlierf06,hanggi90}. We can parametrize the energy unforced landscape $U(x)$ by $N$ \emph{knots} (anchor points), interpolating between them with cubic splines. We'll let the \citetalias{galassi05} natural cubic spline implementation do the heavy lifting here. <>= <> <> <> <> <> @ <>= <> <> <> @ <>= typedef double func_t(double x, void *params); typedef struct kramers_integ_param_struct { double D; /* diffusion rate (um/s) */ double x_min; /* integration bounds */ double x_max; func_t *E; /* function returning E (J) */ void *E_params; /* parametrize the energy landscape */ destroy_data_func_t *destroy_E_params; double F; /* for passing into GSL evaluations */ environment_t *env; } kramers_integ_param_t; @ <>= typedef struct spline_param_struct { int n_knots; /* natural cubic spline knots for U(x) */ double *x; double *y; gsl_spline *spline; /* GSL spline parameters */ gsl_interp_accel *acc; /* GSL spline acceleration data */ } spline_param_t; @ We break the evaluation of $k$ into the two integrals, and use the \citetalias{galassi05} to do the heavy lifting. $$ \text{Integral}_A(x) \equiv \int_0^x e^{\frac{-U_F(x')}{k_B T}} dx' $$ <>= double kramers_integ_k_integrandA(double x, void *params) { kramers_integ_param_t *p = (kramers_integ_param_t *)params; double U = (*p->E)(x, p->E_params); double Fx = p->F * x; double kbT = k_B * p->env->T; //fprintf(stderr, "integrandA(%g) = %g\n", x, exp((U-Fx)/kbT)); return exp(-(U-Fx)/kbT); } double kramers_integ_k_integralA(double x, void *params) { kramers_integ_param_t *p = (kramers_integ_param_t *)params; gsl_function f; double result, abserr; size_t neval; /* for qng */ static gsl_integration_workspace *w = NULL; if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */ f.function = &kramers_integ_k_integrandA; f.params = params; //assert(gsl_integration_qng(&f, p->x_min, x, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0); assert(gsl_integration_qag(&f, p->x_min, x, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0); //fprintf(stderr, "integralA = %g\n", result); return result; } @ $$ \text{Integral}_B(x_\text{min}, x_\text{max}) \equiv \int_{x_\text{min}}^{x_\text{max}} e^{\frac{U_F(x)}{k_B T}} \text{Integral}_A(x) dx, $$ <>= double kramers_integ_k_integrandB(double x, void *params) { kramers_integ_param_t *p = (kramers_integ_param_t *)params; double U = (*p->E)(x, p->E_params); double Fx = p->F * x; double kbT = k_B * p->env->T; double intA = kramers_integ_k_integralA(x, params); //fprintf(stderr, "integrandB(%g) = %g\n", x, exp(-(U-Fx)/kbT))*intA; return exp((U-Fx)/kbT)*intA; } double kramers_integ_k_integralB(double F, environment_t *env, void *params) { kramers_integ_param_t *p = (kramers_integ_param_t *)params; gsl_function f; double result, abserr; size_t neval; /* for qng */ static gsl_integration_workspace *w = NULL; if (w == NULL) w = gsl_integration_workspace_alloc(500); /* for qag */ f.function = &kramers_integ_k_integrandB; f.params = params; p->F = F; p->env = env; //assert(gsl_integration_qng(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, &result, &abserr, &neval) == 0); assert(gsl_integration_qag(&f, p->x_min, p->x_max, 0 /* epsabs */, 1e-6, 500, 6/* key */, w, &result, &abserr) == 0); //fprintf(stderr, "integralB = %g\n", result); return result; } @ With the integrals taken care of, evaluating $k$ is simply $$ k = \frac{D}{\text{Integral}_B(x_\text{min}, x_\text{max})} $$ <>= double kramers_integ_k(double F, environment_t *env, void *kramers_params); @ <>= double kramers_integ_k(double F, environment_t *env, void *kramers_params) { /* Returns the rate constant k in frac pop per s. Takes F in N. */ kramers_integ_param_t *p = (kramers_integ_param_t *) kramers_params; return p->D/kramers_integ_k_integralB(F, env, kramers_params); } @ <>= void *string_create_kramers_integ_param_t(char **param_strings); void destroy_kramers_integ_param_t(void *p); @ <>= kramers_integ_param_t *create_kramers_integ_param_t(double D, double xmin, double xmax, func_t *E, void *E_params, destroy_data_func_t *destroy_E_params) { kramers_integ_param_t *ret = (kramers_integ_param_t *) malloc(sizeof(kramers_integ_param_t)); assert(ret != NULL); ret->D = D; ret->x_min = xmin; ret->x_max = xmax; ret->E = E; ret->E_params = E_params; ret->destroy_E_params = destroy_E_params; return ret; } void *string_create_kramers_integ_param_t(char **param_strings) { //printf("create kramers integ. D: %s, knots: %s, x in [%s, %s]\n", // param_strings[0],param_strings[1],param_strings[2],param_strings[3]); void *E_params = string_create_spline_param_t(param_strings+1); return create_kramers_integ_param_t( safe_strtod(param_strings[0], __FUNCTION__), safe_strtod(param_strings[2], __FUNCTION__), safe_strtod(param_strings[3], __FUNCTION__), &spline_eval, E_params, destroy_spline_param_t); } void destroy_kramers_integ_param_t(void *params) { kramers_integ_param_t *p = (kramers_integ_param_t *)params; if (p) { if (p->E_params) (*p->destroy_E_params)(p->E_params); free(p); } } @ Finally we have the GSL spline wrappers: <>= <> <> @ <>= double spline_eval(double x, void *spline_params) { spline_param_t *p = (spline_param_t *)(spline_params); return gsl_spline_eval(p->spline, x, p->acc); } @ <>= spline_param_t *create_spline_param_t(int n_knots, double *x, double *y) { spline_param_t *ret = (spline_param_t *) malloc(sizeof(spline_param_t)); assert(ret != NULL); ret->n_knots = n_knots; ret->x = x; ret->y = y; ret->acc = gsl_interp_accel_alloc(); ret->spline = gsl_spline_alloc(gsl_interp_cspline, n_knots); gsl_spline_init(ret->spline, x, y, n_knots); return ret; } /* takes an (array of strings) of length one: {"(1,2),(2,0),..."} */ void *string_create_spline_param_t(char **param_strings) { char **knot_coords; int i, num_knots; double *x=NULL, *y=NULL; /* break into ordered pairs */ parse_list_string(param_strings[0], SEP, '(', ')', &num_knots, &knot_coords); x = (double *)malloc(sizeof(double)*num_knots); assert(x != NULL); y = (double *)malloc(sizeof(double)*num_knots); assert(x != NULL); for (i=0; i (%g, %g)\n", i, knot_coords[i], x[i], y[i]); } free_parsed_list(num_knots, knot_coords); return create_spline_param_t(num_knots, x, y); } void destroy_spline_param_t(void *params) { spline_param_t *p = (spline_param_t *)params; if (p) { if (p->spline) gsl_spline_free(p->spline); if (p->acc) gsl_interp_accel_free(p->acc); if (p->x) free(p->x); if (p->y) free(p->y); free(p); } } @ <>= extern int num_kramers_integ_params; extern char *kramers_integ_param_descriptions[]; extern char kramers_integ_param_string[]; @ <>= int num_kramers_integ_params = 4; char *kramers_integ_param_descriptions[] = { "diffusion D, m^2 / s", "spline knots, {(x1,y1),(x2,y2),...}, (m,J)", "starting integration bound x_min, m", "ending integration bound x_max, m"}; //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"; //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"; 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"; /* D from Schlierf 2006, Biophys 90, 4, L33, sup. * Anchor points from private communication with Schlierf, Apr 9th, 2008. * other parameters from from Schlierf 2006, Biophys 90, 4, L33, figure 1. */ @ <>= {"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} @ Initialized with Titin I27 parameters\citep{carrion-vazquez99a}. % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696 \subsection{Unfolding model implementation} <>= <> <> #include "k_model.h" <> <> <> @ <>= #include /* assert() */ #include /* NULL, malloc(), strto*() */ #include /* HUGE_VAL, sqrt(), exp() */ #include /* fprintf(), stdout */ #include /* strlen(), strcpy() */ #include /* errno, ERANGE */ #include #include #include "global.h" #include "parse.h" @ <>= <> <> <> <> <> <> @ <>= <> <> <> <> <> <> @ <>= <> <> <> <> <> <> <> @ \subsection{Unfolding model unit tests} Here we check to make sure the various functions work as expected, using \citetalias{sw:check}. <>= <> <> <> <> <> <
> @ <>= #include /* EXIT_SUCCESS, EXIT_FAILURE */ #include /* sprintf() */ #include /* assert() */ #include /* exp() */ #include /* errno, ERANGE */ <> #include "global.h" #include "k_model.h" @ <>= <> <> <> <> @ <>= Suite *test_suite (void) { Suite *s = suite_create ("k model"); <> <> <> <> return s; } @ \subsubsection{Constant} <>= TCase *tc_const_k = tcase_create("Constant k"); @ <>= tcase_add_test(tc_const_k, test_const_k_create_destroy); tcase_add_test(tc_const_k, test_const_k_over_range); suite_add_tcase(s, tc_const_k); @ <>= START_TEST(test_const_k_create_destroy) { char *knot[] = {"1","2","3","4","5","6"}; char *params[] = {knot[0]}; void *p = NULL; int i; for( i=0; i < sizeof(knot)/sizeof(char *); i++) { params[0] = knot[i]; p = string_create_const_k_param_t(params); destroy_const_k_param_t(p); } } END_TEST START_TEST(test_const_k_over_range) { double F, Fm=0.0, FM=10, dF=.1, T=300.0; char *knot[] = {"1","2","3","4","5","6"}; char *params[] = {knot[0]}; void *p = NULL; environment_t env; char param_string[80]; int i; env.T = T; for( i=0; i < sizeof(knot)/sizeof(char *); i++) { params[0] = knot[i]; p = string_create_const_k_param_t(params); for ( F=Fm; F>= TCase *tc_bell = tcase_create("Bell's k"); @ <>= tcase_add_test(tc_bell, test_bell_k_create_destroy); tcase_add_test(tc_bell, test_bell_k_at_zero); tcase_add_test(tc_bell, test_bell_k_at_one); suite_add_tcase(s, tc_bell); @ <>= START_TEST(test_bell_k_create_destroy) { char *knot[] = {"1","2","3","4","5","6"}; char *dx="1"; char *params[] = {knot[0], dx}; void *p = NULL; int i; for( i=0; i < sizeof(knot)/sizeof(char *); i++) { params[0] = knot[i]; p = string_create_bell_param_t(params); destroy_bell_param_t(p); } } END_TEST START_TEST(test_bell_k_at_zero) { double F=0.0, T=300.0; char *knot[] = {"1","2","3","4","5","6"}; char *dx="1"; char *params[] = {knot[0], dx}; void *p = NULL; environment_t env; char param_string[80]; int i; env.T = T; for( i=0; i < sizeof(knot)/sizeof(char *); i++) { params[0] = knot[i]; p = string_create_bell_param_t(params); fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__), "bell_k(%g, %g, &{%s,%s}) = %g != %s", F, T, knot[i], dx, bell_k(F, &env, p), knot[i]); destroy_bell_param_t(p); } } END_TEST START_TEST(test_bell_k_at_one) { double T=300.0; char *knot[] = {"1","2","3","4","5","6"}; char *dx="1"; char *params[] = {knot[0], dx}; double F= k_B*T/safe_strtod(dx, __FUNCTION__); void *p = NULL; environment_t env; char param_string[80]; int i; env.T = T; for( i=0; i < sizeof(knot)/sizeof(char *); i++) { params[0] = knot[i]; p = string_create_bell_param_t(params); CHECK_ERR(1e-6, safe_strtod(knot[i], __FUNCTION__)*exp(1.0), bell_k(F, &env, p)); //fail_unless(bell_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__)*exp(1.0), // "bell_k(%g, %g, &{%s,%s}) = %g != %g", // F, T, knot[i], dx, bell_k(F, &env, p), safe_strtod(knot[i], __FUNCTION__)*exp(1.0)); destroy_bell_param_t(p); } } END_TEST @ <>= @ <>= @ <>= @ <>= <> START_TEST(test_something) { double k=5, x=3, last_x=2.0, F; one_dim_fn_t *handlers[] = {&hooke}; void *data[] = {&k}; double xi[] = {0.0}; int active[] = {1}; int new_active_groups = 1; F = k_model(1, handlers, data, active, new_active_groups, xi, last_x, x); fail_unless(F = k*x, NULL); } END_TEST @ \subsection{Utilities} The unfolding models can get complicated, and are sometimes hard to explain to others. The stand-alone program [[k_model_utils]] demonstrates the unfolding model interface without the overhead of having to understand a full simulation. [[k_model_utils]] can operate in a general mode, where it prints $k(F)$ data to the screen, or special mode, where the operation depends on the specific model selected. <>= <> <> <> <> <> <> @ <>= int main(int argc, char **argv) { <> k_model_getopt_t *model = NULL; void *params; enum mode_t mode; environment_t env; unsigned int flags; int num_param_args; /* for INIT_MODEL() */ char **param_args; /* for INIT_MODEL() */ double Fmax; double special_xmin; double special_xmax; get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model, &Fmax, &special_xmin, &special_xmax, &flags); if (flags & VFLAG) { printf("#initializing model %s with parameters %s\n", model->name, model->params); } INIT_MODEL("utility", model, model->params, params); switch (mode) { case M_K_OF_F : if (model->k == NULL) { printf("No k function for model %s\n", model->name); exit(0); } if (Fmax <= 0) { printf("Fmax = %g <= 0\n", Fmax); exit(1); } { double F, k=1.0; int i, N=200; printf("#F (N)\tk (%% pop. per s)\n"); for (i=0; i<=N; i++) { F = Fmax*i/(double)N; k = (*model->k)(F, &env, params); if (k < 0) break; printf("%g\t%g\n", F, k); } } break; case M_SPECIAL : if (model->k == NULL || model->k == &bell_k) { printf("No special mode for model %s\n", model->name); exit(1); } if (special_xmax <= special_xmin) { printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin); exit(1); } { double Xrng=(special_xmax-special_xmin), x, E; int i, N=500; printf("#x (m)\tE (J)\n"); for (i=0; i<=N; i++) { x = special_xmin + Xrng*i/(double)N; E = multi_model_E(model->k, params, &env, x); printf("%g\t%g\n", x, E); } } break; default : break; } if (model->destructor) (*model->destructor)(params); return 0; } @ <>= double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x) { kramers_integ_param_t *p = (kramers_integ_param_t *)model_params; double E; if (k == kramers_integ_k) E = (*p->E)(x, p->E_params); else E = kramers_E(0, env, model_params, x); return E; } @ <>= #include #include #include /* strlen() */ #include /* assert() */ #include /* errno, ERANGE */ #include "global.h" #include "parse.h" #include "string_eval.h" #include "k_model.h" @ <>= <> #define VFLAG 1 /* verbose */ enum mode_t {M_K_OF_F, M_SPECIAL}; <> <> <> <> <> @ <>= <> <> <> <> @ <>= void help(char *prog_name, environment_t *env, int n_k_models, k_model_getopt_t *k_models, int k_model, double Fmax, double special_xmin, double special_xmax) { int i, j; printf("usage: %s [options]\n", prog_name); printf("Version %s\n\n", VERSION); printf("A utility for understanding the available unfolding models\n\n"); printf("Environmental options:\n"); printf("-T T\tTemperature T (currently %g K)\n", env->T); printf("-C T\tYou can also set the temperature T in Celsius\n"); printf("Model options:\n"); printf("-k model\tTransition rate model (currently %s)\n", k_models[k_model].name); printf("-K args\tTransition rate model argument string (currently %s)\n", k_models[k_model].params); printf("There are two output modes. In standard mode, k(F) is printed\n"); printf("For example:\n"); printf(" #Force (N)\tk (%% pop. per s)\n"); printf(" 123.456\t7.89\n"); printf(" ...\n"); printf("In special mode, the output depends on the model.\n"); printf("For models defining an energy function E(x), that function is printed\n"); printf(" #Distance (m)\tE (J)\n"); printf(" 0.0012\t0.0034\n"); printf(" ...\n"); printf("-m\tChange output to standard mode\n"); printf("-M\tChange output to special mode\n"); printf("-F\tSet the maximum F value for the standard mode k(F) output (currently %g)\n", Fmax); printf("-x\tSet the minimum x value for the special mode E(x) output (currently %g)\n", special_xmin); printf("-X\tSet the maximum x value for the special mode E(x) output (currently %g)\n", special_xmax); printf("-V\tChange output to verbose mode\n"); printf("-h\tPrint this help and exit\n"); printf("\n"); printf("Unfolding rate models:\n"); for (i=0; i>= void get_options(int argc, char **argv, environment_t *env, int n_k_models, k_model_getopt_t *k_models, enum mode_t *mode, k_model_getopt_t **model, double *Fmax, double *special_xmin, double *special_xmax, unsigned int *flags) { char *prog_name = NULL; char c, options[] = "T:C:k:K:mMF:x:X:Vh"; int k_model=0; extern char *optarg; extern int optind, optopt, opterr; assert(argc > 0); /* setup defaults */ prog_name = argv[0]; env->T = 300.0; /* K */ *mode = M_K_OF_F; *flags = 0; *model = k_models; *Fmax = 1e-9; /* N */ *special_xmax = 1e-8; *special_xmin = 0.0; while ((c=getopt(argc, argv, options)) != -1) { switch(c) { case 'T': env->T = safe_strtod(optarg, "-T"); break; case 'C': env->T = safe_strtod(optarg, "-C")+273.15; break; case 'k': k_model = index_k_model(n_k_models, k_models, optarg); *model = k_models+k_model; break; case 'K': k_models[k_model].params = optarg; break; case 'm': *mode = M_K_OF_F; break; case 'M': *mode = M_SPECIAL; break; case 'F': *Fmax = safe_strtod(optarg, "-F"); break; case 'x': *special_xmin = safe_strtod(optarg, "-x"); break; case 'X': *special_xmax = safe_strtod(optarg, "-X"); break; case 'V': *flags |= VFLAG; break; case '?': fprintf(stderr, "unrecognized option '%c'\n", optopt); /* fall through to default case */ default: help(prog_name, env, n_k_models, k_models, k_model, *Fmax, *special_xmin, *special_xmax); exit(1); } } return; } @ \section{\LaTeX\ documentation} The comment blocks in this [[noweb]] file are in \LaTeX\ with \BibTeX\ citations. The comment blocks are extracted (with nicely formatted code blocks), using <>= $(BUILD_DIR)/sawsim.tex : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR) noweave -latex -index -delay $< > $@ $(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR) cp -f $< $@ @ and compiled using <>= $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \ | $(DOC_DIR) $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim" $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim" $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim" $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim" $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim" $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim" mv $(BUILD_DIR)/sawsim.pdf $@ @ The first [[pdflatex]] pass reads through the \TeX\ source and extracts citation information, the \BibTeX\ pass looks through the \BibTeX\ databases for the cited works, and the remaining two [[pdflatex]] passes insert the citations and finalize the various indices and numberings. [[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up. <>= semi-clean_latex : rm -f $(BUILD_DIR)/sawsim.aux $(BUILD_DIR)/sawsim.bbl \ $(BUILD_DIR)/sawsim.blg $(BUILD_DIR)/sawsim.log \ $(BUILD_DIR)/sawsim.out $(BUILD_DIR)/sawsim.toc \ $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib clean_latex : semi-clean_latex rm -f $(DOC_DIR)/sawsim.pdf @ \section{Makefile} [[make]] is a common utility on *nix systems for managing dependencies while building software. In this case, we have several \emph{targets} that we'd like to build: the main simulation program \prog; the testing programs [[check_sawsim]], [[check_tension_balance]], and [[check_list]]; the documentation [[sawsim.pdf]]; and the [[Makefile]] itself. Besides the generated files, there is the utility target [[clean]] for removing all generated files except [[Makefile]]. % and % [[dist]] for generating a distributable tar file. Extract the makefile with `[[notangle -Rmakefile src/sawsim.nw | sed 's/ /\t/' > Makefile]]'. The sed is needed because notangle replaces tabs with eight spaces, but [[make]] needs tabs. <>= # Customize compilation CC = gcc CFLAGS = LDFLAGS = # Decide what directories to put things in SRC_DIR = src BUILD_DIR = build BIN_DIR = bin DOC_DIR = doc # Note: Cannot use, for example, `./build', because make eats the `./' # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...) # Modules (X.c, X.h) defined in the noweb file NW_SAWSIM_MODS = # Check binaries (check_*) defined in the noweb file CHECK_BINS = # Modules defined outside the noweb file FREE_SAWSIM_MODS = interp tavl <> <> <> <> <> <> <> SAWSIM_MODS = $(NW_SAWSIM_MODS) $(FREE_SAWSIM_MODS) # Everything needed for sawsim under one roof. sawsim.c must come first SAWSIM_SRC = $(BUILD_DIR)/sawsim.c $(BUILD_DIR)/global.h \ $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) # Libraries needed to compile sawsim LIBS = gsl gslcblas m CHECK_LIBS = gsl gslcblas m check # The non-check binaries generated BINS = sawsim tension_model_utils k_model_utils sawsim_profile DOCS = sawsim.pdf # Define the major targets all : ./Makefile all_bin all_doc .PHONY: all_bin all_doc all_bin : $(BINS:%=$(BIN_DIR)/%) all_doc : $(DOCS:%=$(DOC_DIR)/%) view : $(DOC_DIR)/sawsim.pdf xpdf $< & profile : $(BIN_DIR)/sawsim_profile | $(BIN_DIR) $(BIN_DIR)/sawsim_profile -v1e-6 -Mhooke -A.05 -U1 -kkramers_integ \ -mnull -Mwlc -A0.39e-9,28e-9 -F8 gprof $(BIN_DIR)/sawsim_profile gmon.out > $@ check : $(CHECK_BINS:%=$(BIN_DIR)/%) $(BIN_DIR)/check_sawsim $(SHELL) -e -c 'for B in $^; do ./$$B; done' clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \ clean_tension_model_utils \ clean_k_model_utils clean_latex clean_check_sawsim rm -f $(BIN_DIR)/sawsim $(BIN_DIR)/sawsim_static \ $(BIN_DIR)/sawsim_profile $(BUILD_DIR)/sawsim.c \ $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \ $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \ $(BUILD_DIR)/global.h ./gmon.out $(SHELL) -e -c "rm -rf $(BUILD_DIR) $(DOC_DIR)" # Various builds of sawsim $(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR) $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%) $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR) $(CC) $(CFLAGS) $(LDFLAGS) -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%) $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR) $(CC) $(CFLAGS) $(LDFLAGS) -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%) # Create the directories $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) : mkdir $@ # Copy over the source external to sawsim # Note: Cannot use, for example, `./build', because make eats the `./' # and then $(BUILD_DIR) doesn't match in $$(subst $(BUILD_DIR), ...) .SECONDEXPANSION: $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\ | $(BUILD_DIR) cp -f $< $@ .SECONDEXPANSION: $(FREE_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $$(subst $(BUILD_DIR),$(SRC_DIR),$$@)\ | $(BUILD_DIR) cp -f $< $@ ## Basic source generated with noweb # The central files sawsim.c and global.h... $(BUILD_DIR)/sawsim.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR) notangle $< > $@ $(BUILD_DIR)/global.h : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR) notangle -Rglobal.h $< > $@ # ... and the modules $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.h) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR) notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@ $(NW_SAWSIM_MODS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR) notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@ $(CHECK_BINS:%=$(BUILD_DIR)/%.c) : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR) notangle -R$(subst _,-,$(@:$(BUILD_DIR)/%=%)) $< > $@ # Note: I use `_' as a space in my file names, but noweb doesn't like # that so we use `-' to name the noweb roots and substitute here. ## Building the unit-test programs # for sawsim.c... $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw notangle -Rchecks $< > $@ $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \ $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \ $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR) $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%) clean_check_sawsim : rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c # ... and the modules .SECONDEXPANSION: $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \ $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).c \ $$(subst $(BIN_DIR)/check_,$(BUILD_DIR)/,$$@).h \ $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\ $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\ $$(BUILD_DIR)/global.h | $(BIN_DIR) $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \ $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\ $(CHECK_LIBS:%=-l%) # TODO: clean up the required modules too $(CHECK_BINS:%=clean_%) : rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c) # Cleaning up the modules .SECONDEXPANSION: $(SAWSIM_MODS:%=clean_%) : rm -f $(@:clean_%=$(BUILD_DIR)/%.c) $(@:clean_%=$(BUILD_DIR)/%.h) <> <> <> Makefile : $(SRC_DIR)/sawsim.nw notangle -Rmakefile $< | sed 's/ /\t/' > $@ @ This is fairly self-explanatory. We compile a dynamically linked version ([[sawsim]]) and a statically linked version ([[sawsim_static]]). 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. \section{Math} \subsection{Hookean springs in series} \label{sec.math_hooke} \begin{theorem} The effective spring constant for $n$ springs $k_i$ in series is given by $$ k_\text{eff} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}. $$ \end{theorem} \begin{proof} \begin{align} F &= k_i x_i &&\forall i \le n \\ x_i &= \frac{k_j}{k_i} x_j &&\forall i,j \le n \\ x &= \sum_i x_i = x_1 \sum_i \frac{k_1}{k_i} \\ F &= k_1 x_1 = k_\text{eff} x \\ k_\text{eff} &= k_1 \frac{x_1}{x} = \frac{k_1}{\sum_i \frac{k_1}{k_i}} = \p({\sum_{i \le n}\frac{1}{k_i}})^{-1}. \end{align} \end{proof} \phantomsection \addcontentsline{toc}{section}{References} \bibliography{sawsim} \end{document}