+% 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 <http://www.gnu.org/licenses/>.
+%
+% The author may be contacted at <wking@drexel.edu> 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,
+% 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[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{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
\today \\
\end{centering}
\bigskip
+
+\tableofcontents
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%% End LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\subsection{Related work}
-TODO References
+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}
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.
+
<<version definition>>=
-#define VERSION "0.4"
-@
+#define VERSION "0.10"
+@
\subsection{License}
<<license>>=
- sawsim - program for simulating single molecule mechanical unfolding.
- Copyright (C) 2008, William Trevor King
+ 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
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.
+ 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, write to the Free Software
- Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
- 02111-1307, USA.
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
The author may be contacted at <wking@drexel.edu> on the Internet, or
write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
Philadelphia PA 19104, USA.
-@
+@
<<license comment>>=
/*
<<license>>
*/
-@
+@
\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}).
-The unfolding system is stored as a chain of \emph{domains}. Domains
-can be folded, globular protein domains, unfolded protein linkers, AFM
-cantilevers, or other stretchable system components.
+\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
\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 whether the domain is folded or unfolded. The transition rate from
-the folded state to the unfolded state is also handled generally with
-function pointers, with implementations for
+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), is applied to each folded
-domain, and used to calculate the probability of that subdomain
-unfolding in the next timestep $dt$. Then the time is stepped
-forward, and the process is repeated.
-
+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.
<<main program>>=
int main(int argc, char **argv)
{
<<tension model getopt array definition>>
<<k model getopt array definition>>
- list_t *domain_list=NULL; /* variables initialized in get_options() */
+
+ /* variables initialized in get_options() */
+ list_t *state_list=NULL, *transition_list=NULL;
environment_t env={0};
- double P, dt_max, v, xstep;
- int num_folded=0, unfolding;
- double x=0, dt, F; /* variables used in the simulation loop */
- get_options(argc, argv, &P, &dt_max, &v, &xstep, &env, NUM_TENSION_MODELS,
- tension_models, NUM_K_MODELS, k_models, &domain_list, &num_folded);
+ 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("#Position (m)\tForce (N)\n");
- if (flags & FLAG_OUTPUT_UNFOLDING_FORCES) printf("#Unfolding Force (N)\n");
- while (num_folded > 0) {
- F = find_tension(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, unfolding);
- if (xstep == 0)
- dt = determine_dt(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, v);
+
+ 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(NUM_TENSION_MODELS, tension_handlers, domain_list, &env, x, P, dt_max, 0);
- unfolding=random_unfoldings(domain_list, F, dt, &env, &num_folded);
- if (flags & FLAG_OUTPUT_FULL_CURVE) printf("%g\t%g\n", x, F);
- if (xstep == 0) {
+ 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 */
- fprintf(stderr, "step completed\n");
- x += xstep;
- dt_max = xstep / v;
+ x += x_step;
+ dt_max = x_step / v;
} else { /* still working on this step */
- fprintf(stderr, "partial step %g of %g\n", dt, dt_max);
dt_max -= dt;
}
}
}
- destroy_domain_list(domain_list);
+ 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_unfoldings]].
+[[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_unfoldings]] in Section \ref{sec.unfolding_rate}.
+[[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
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
+tensions (e.g.\ temperature, pH) are stored in a single structure to
facilitate extension to more complicated models in the future.
-<<environment definition>>=
+<<environment definition>>=
typedef struct environment_struct {
double T;
+ double stiffness;
} environment_t;
-@
+@
<<global.h>>=
#ifndef GLOBAL_H
#define GLOBAL_H
+
+#define DOUBLE_PRECISION 1e-12
+
<<environment definition>>
<<one dimensional function definition>>
<<create/destroy definitions>>
<<model globals>>
#endif /* GLOBAL_H */
-@
+@
Where the [[GLOBAL_H]] defines are protecting [[global.h]] from being
included multiple times.
<<find tension>>
<<determine dt>>
<<happens>>
-<<does domain unfold>>
-<<randomly unfold domains>>
-@
+<<does domain transition>>
+<<randomly transition domains>>
+@
\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 models and groups. The models are used to
-determine a domain's tension (Hookean spring, WLC, \ldots). Groups
-will represent collections of domains which the model should treat as
-a single entity. A domain may belong to different groups or models
-depending on its state. For example, a domain might be modeled as a
-freely-jointed chain when it iss folded and as a worm-like chain class
-when it is unfolded. 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. This allows for non-linear group models such as th
-worm-like or freely-jointed chains. Each model has a tension handler
-function, which gives the tension $F$ needed to stretch a given group
-of domains a distance $x$.
-
-To understand the purpose of groups, consider a chain of proteins
-where two folded proteins $A$ and $B$ are modeled as WLCs with
-persistence length $p_f$ and two unfolded proteins $C$ and $D$ are
-modeled as WLCs with persistence length $p_u$. The proteins are
-attached to a cantilever $E$ qof spring constant $κ$. The string
-would get split into three lists:
-\begin{center}
- \begin{tabular}{llll}
- Model & Group & List & Tension \\
- WLC & 0 & $[A,B]$ & $F_\text{WLC}(x, p_f, L_A+L_B)$ \\
- WLC & 1 & $[C,D]$ & $F_\text{WLC}(x, p_u, L_C+L_D)$ \\
- Hooke & 0 & $[E]$ & $F_\text{Hooke}(x, κ)$ \\
- \end{tabular}
-\end{center}
-Note that group numbers only matter \emph{within} model classes, since
-grouping domains with seperate models doesn't make sense.
+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.
<<find tension definitions>>=
-#define NUM_TENSION_MODELS 5
+#define NUM_TENSION_MODELS 6
-@
+@
<<tension handler array global declaration>>=
extern one_dim_fn_t *tension_handlers[];
@
-<<tension handler array global>>=
+<<tension handler array global>>=
one_dim_fn_t *tension_handlers[] = {
NULL,
&const_tension_handler,
&fjc_handler,
};
-@
+@
<<tension model global declarations>>=
<<tension handler array global declaration>>
-@
+@
<<tension handler types>>=
typedef struct tension_handler_data_struct {
- /* int num_domains; */
- /* domain_t *domains; */
- list_t *group;
+ 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 groups within each model
-class. [[find_tension]] is basically a wrapper to feed proper input
-into the [[tension_balance]] function. [[unfolding]] should be set to
-0 if there were no unfoldings since the previous call to
-[[find_tension]].
+\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}).
<<find tension>>=
-double find_tension(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list, environment_t *env, double x, int unfolding)
-{
+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 void **per_group_params = 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;
+ double F, *pStiffness;
int i;
#ifdef DEBUG
- fprintf(stderr, "finding tension at distance %g%s\n", x, unfolding ? " after an unfolding" : "");
- list_print(stderr, domain_list, "domain list");
+ fprintf(stderr, "%s:\tfinding tension at distance %g%s\n", __FUNCTION__, x, transition ? " after a transition" : "");
#endif
- if (unfolding != 0 || num_active_groups == 0) { /* setup statics */
- /* free old statics */
- if (per_group_handlers != NULL)
- free(per_group_handlers);
- if (per_group_params != NULL) {
- for (i=0; i < num_active_groups; i++) {
- assert(per_group_params[i] != NULL);
- free(per_group_params[i]);
- }
- free(per_group_params);
- }
-
- /* get new statics */
- get_active_groups(domain_list, num_tension_handlers, tension_handlers, &num_active_groups, &per_group_handlers, &per_group_params);
+ 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; i<num_active_groups; i++) {
- tdata = (tension_handler_data_t *) per_group_params[i];
+ tdata = (tension_handler_data_t *) per_group_data[i];
#ifdef DEBUG
- fprintf(stderr, "active group %d of %d", i, num_active_groups);
- list_print(stderr, tdata->group, " parameters");
+ fprintf(stderr, "%s:\tactive group %d of %d (data %p, param list %p)", __FUNCTION__, i, num_active_groups, tdata, tdata->group_tension_params);
+ list_print(stderr, tdata->group_tension_params, " parameters");
#endif
tdata->env = env;
/* tdata->persist continues unchanged */
}
last_x = -1.0;
- } /* else roll with the current statics */
+ } /* 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);
- F = tension_balance(num_active_groups, per_group_handlers, (void **)per_group_params, last_x, x);
last_x = x;
+
return F;
}
@ For the moment, we will restrict our group boundaries to a single
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{app.tension_balance}
+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 it's data argument.
+ 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. See the Section \ref{app.model_selection} for a
-discussion on the command line framework. See the worm-like chain
+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{Unfolding rate}
-\label{sec.unfolding_rate}
+\subsection{Transition rate}
+\label{sec.transition_rate}
-Each folded domain is modeled as unimolecular, first order reaction
+Each state transition is modeled as unimolecular, first order reaction
$$
- \text{Folded} \xrightarrow{k} % in tex: X atop Y
- \text{Unfolded}
+ \text{State 1} \xrightarrow{k} \text{State 2}
$$
With the rate constant $k$ defined by
$$
- \frac{d[\text{Folded}]}{dt} = -k [\text{Folded}].
+ \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
$$
-So the probability of a given protein unfolding in a short time $dt$
-is given by
+ P_1 = k \cdot dt.
$$
- P = 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
+<<PN definition>>=
+/* find multi-domain transition rate from the single domain rate */
+#define PN(P,N) (1.0-pow(1.0-(P), (N)))
-<<does domain unfold>>=
-int domain_unfolds(double F, double dt, environment_t *env, domain_t *domain)
-{ /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to a folded domain */
+@
+
+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}
+<<P1 definition>>=
+#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}.
+
+<<does domain transition>>=
+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(domain->k, F, env, domain->k_params);
- //(*domain->k)(F, env, domain->k_params);
+ 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(k*dt); /* dice roll for prob. k*dt event */
-}
-@ [[happens]] is a random decision making function defined in Appendix \ref{app.utils}.
-
-Only one domain can unfold in each timestep, because the timescale of
-a domain unfolding $dt_u$ is assumed to be less than the simulation
-timestep $dt$, so a domain will completely unfold in a single
-timestep. We adapt our timesteps to keep the probability of a single
-domain unfolding low, so the probability of two domains unfolding in
-the same timestep is negligible. This approach breaks down as the
-adaptive timestep scheme approaches $dt \sim dt_u$, but $dt_u \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.
-<<randomly unfold domains>>=
-int random_unfoldings(list_t *domain_list, double tension,
- double dt, environment_t *env,
- int *num_folded_domains)
-{ /* returns 1 if there was an unfolding and 0 otherwise */
- while (domain_list != NULL) {
- if (D(domain_list)->state == DS_FOLDED
- && domain_unfolds(tension, dt, env, D(domain_list))) {
- if (flags & FLAG_OUTPUT_UNFOLDING_FORCES)
- fprintf(stdout, "%g\n", tension);
- D(domain_list)->state = DS_UNFOLDED;
- (*num_folded_domains)--;
- return 1; /* our one domain has unfolded, stop looking for others */
+ return happens(PN(k*dt, N_INIT(transition)));
+}
+@ [[happens]] is a random decision making function defined in Appendix \ref{sec.utils}.
+
+<<randomly transition domains>>=
+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 */
}
- domain_list = domain_list->next;
+ 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 unfolding in the next
-timestep is small. If we don't adapt our timestep, we risk breaking
-our approximation that $F(x' \in [x, x+v \cdot dt]) = F(x)$ or that
-only one domain unfolds in a given timestep. Because $F(x)$ is
-monotonically increasing, excessively large timesteps will lead to
-erroneously large unfolding forces. 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.
+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{app.adaptive_dt}.
+\ref{sec.adaptive_dt_implementation}.
\section{Models}
-
+
TODO: model intro.
The models provide the physics of the simulation, but the simulation
simulation itself, so we remove the actual model implementation to the
appendices.
-The tension models are defined in Appendix \ref{sec.tension_models}
+The tension models are defined in Appendix \ref{sec.tension_models},
and unfolding models are defined in Appendix \ref{sec.k_models}.
<<model globals>>=
#define k_B 1.3806503e-23 /* J/K */
-@
+
+@
\section{Command line interface}
<<help>>
<<index tension model>>
<<index k model>>
-<<generate domain>>
<<get options>>
-@
+@
\subsection{Model selection}
-\label{app.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
<<get options definitions>>=
<<model getopt definitions>>
-@
+@
<<create/destroy definitions>>=
typedef void *create_data_func_t(char **param_strings);
typedef void destroy_data_func_t(void *param_struct);
-@
+@
<<model getopt definitions>>=
<<tension model getopt definitions>>
<<k model getopt definitions>>
-@
+@
+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 <state-name>,<model>[,<group>] -S <model-params> [-N <initial-population]]] \\
+\end{center}
+e.g.
+\begin{center}
+[[-s unfolded,wlc,1 -S 1e-8,4e-10 -N 5]]
+\end{center}
+This creates the state named [[unfolded]], using the WLC model and one
+as a group number (defaults to 0, see Section \ref{sec.find_tension}).
+The tension parameters are then set to [[1e-8,4e-10]], which in the
+case of the WLC are the contour and persistence lengths respectively.
+Finally we populate the state by adding five domains to the state.
+The name of the state is importent for identifying states later.
+
+Once the domains have been created and populated, you can added
+transition rates following
+\begin{center}
+ [[-k <initial-state>,<final-state>,<model> -K <model-params>]] \\
+\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}
model.
<<tension model getopt definitions>>=
typedef struct tension_model_getopt_struct {
- one_dim_fn_t *handler; /* fn implementing the model on a group */
+ 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 */
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 array definition>>=
tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
<<constant tension model getopt>>,
<<hooke tension model getopt>>,
<<worm-like chain tension model getopt>>,
- <<freely-jointed chain tension model getopt>>
+ <<freely-jointed chain tension model getopt>>,
+ <<piston tension model getopt>>
};
-@
+@
\subsubsection{Unfolding rate}
<<k model getopt definitions>>=
-#define NUM_K_MODELS 5
+#define NUM_K_MODELS 6
typedef struct k_model_getopt_struct {
char *name;
create_data_func_t *creator;
destroy_data_func_t *destructor;
} k_model_getopt_t;
-@
+@
<<k model getopt array definition>>=
k_model_getopt_t k_models[NUM_K_MODELS] = {
<<null k model getopt>>,
<<const k model getopt>>,
<<bell k model getopt>>,
+ <<kbell k model getopt>>,
<<kramers k model getopt>>,
<<kramers integ k model getopt>>
};
-@
+@
\subsection{help}
<<help>>=
-void help(char *prog_name, double P, double dt_max, double v, double xstep,
- environment_t *env,
+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 folded_tension_model, int unfolded_tension_model,
int n_k_models, k_model_getopt_t *k_models,
- int k_model)
+ 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("-P P\tTarget probability for dt (currently %g)\n", P);
- printf("-t dt\tMaximum allowed timestep dt (currently %g)\n", dt_max);
+ 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", xstep);
- printf("\tWhen dx = 0, the pulling is continuous\n");
- printf("\tWhen dx > 0, the pulling occurs in discrete steps\n");
+ 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("Model options:\n");
- printf("The domains exist in either the folded or unfolded state\n");
- printf("The following options change the default behavior in each state.\n");
- printf("-m model[,group]\tFolded tension model (currently %s)\n",
- tension_models[folded_tension_model].name);
- printf("-a args\tFolded tension model argument string (currently %s)\n",
- tension_models[folded_tension_model].params);
- printf("-M model[,group]\tUnfolded tension model (currently %s)\n",
- tension_models[unfolded_tension_model].name);
- printf("-A args\tUnfolded tension model argument string (currently %s)\n",
- tension_models[unfolded_tension_model].params);
- printf("The following options change the unfolding rate.\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("Domain creation options:\n");
- printf("Once you've set up the appropriate default models, you need to add the domains\n");
- printf("-F n\tAdd n folded domains with the current models\n");
- printf("-U n\tAdd n unfolded domains with the current models\n");
- printf("Output mode options:\n");
- printf("There are two output modes. In standard mode, only the unfolding\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(" #Unfolding Force (N)\n");
- printf(" 123.456e-12\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(" #Position (m)\tForce (N)\n");
- printf(" 0.001\t0.0023\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("-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<n_tension_models; i++) {
printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
printf("\t\t\t%s\n", k_models[i].param_descriptions[j]);
printf("\t\tdefaults: %s\n", k_models[i].params);
}
+
printf("Examples:\n");
- printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded\n");
- printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K3.3e-4,.25e-9 -mnull -Mwlc -A0.39e-9,28e-9 -F8\n", prog_name);
- printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded\n");
- printf(" $ %s -v1e-6 -Mhooke -A.05 -U1 -kbell -K5e-5,.225e-9 -mwlc -a1e-8,4e-10 -Mwlc,1 -A0.4e-9,28.1e-9 -F8\n", prog_name);
+ printf("1 spring and 8 folded domains of titin I27 only stretchable when unfolded. Stop once there are no more folded domains.\n");
+ printf(" $ %s -v1e-6 -s cantilever,hooke,0.05 -N1 -s folded,null -N8 -s 'unfolded,wlc,{0.39e-9,28e-9}' -k 'folded,unfolded,bell,{3.3e-4,0.25e-9}' -q folded\n", prog_name);
+ printf("1 spring and 8 folded domains of ubiquitin with different WLC parameters when folded or unfolded. Stop after 0.25 seconds\n");
+ printf(" $ %s -v1e-6 -s cantilever,hooke,0.05 -N1 -s 'folded,wlc,{1e-8,4e-10}' -N8 -s 'unfolded,wlc,1,{0.4e-9,28.1e-9}' -k 'folded,unfolded,bell,{5e-5,0.225e-9}' -t0.25\n", prog_name);
}
-@
+@
\subsection{Get options}
+\label{sec.get_options}
<<get options>>=
-void get_options(int argc, char **argv,
- double *pP, double *pDt_max, double *pV, double *pXstep,
- environment_t *env,
+<<safe strto*>>
+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 **pDomain_list, int *pNum_folded)
+ list_t **pState_list, list_t **pTransition_list)
{
char *prog_name = NULL;
- char c, options[] = "P:t:v:x:T:C:m:a:M:A:k:K:F:U:Vh";
- int ftension_model=0, utension_model=0, k_model=0;
- char *ftension_params=NULL, *utension_params=NULL;
- int i, n, ftension_group=0, utension_group=0;
+ 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);
+ assert(argc > 0);
#ifdef DEBUG
for (i=0; i<n_tension_models; i++) {
- fprintf(stderr, "tension model %i %s:\t model's handler %p\n",
- i, tension_models[i].name, tension_models[i].handler);
+ fprintf(stderr, "tension model %i %s:\t model's handler %p\t inverse handler %p\n",
+ i, tension_models[i].name, tension_models[i].handler, tension_models[i].inverse_handler);
assert(tension_models[i].handler == tension_handlers[i]);
}
#endif
/* setup defaults */
- flags = FLAG_OUTPUT_UNFOLDING_FORCES;
+ flags = FLAG_OUTPUT_TRANSITION_FORCES;
prog_name = argv[0];
- *pP = 1e-3; /* % pop per s */
- *pDt_max = 0.001; /* s */
- *pV = 1e-6; /* m/s */
- *pXstep = 0.0; /* m */
- env->T = 300.0; /* K */
+ *pStop_state = NULL;
+ *pT_max = -1; /* s, -1 removes this limit */
+ *pDt_max = 0.001; /* s */
+ *pP = 1e-3; /* % pop per s */
+ *pDF = 1e-12; /* N */
+ *pV = 1e-6; /* m/s */
+ *pX_step = 0.0; /* m */
+ env->T = 300.0; /* K */
+ *pState_list = NULL;
+ *pTransition_list = NULL;
while ((c=getopt(argc, argv, options)) != -1) {
switch(c) {
- case 'P': *pP = atof(optarg); break;
- case 't': *pDt_max = atof(optarg); break;
- case 'v': *pV = atof(optarg); break;
- case 'x': *pXstep = atof(optarg); break;
- case 'T': env->T = atof(optarg); break;
- case 'C': env->T = atof(optarg)+273.15; break;
- case 'm':
- parse_list_string(optarg, SEP, '{', '}', &num_strings, &strings);
- assert(num_strings == 1 || num_strings == 2);
- if (num_strings == 1)
- ftension_group = 0;
- else
- ftension_group = atoi(strings[1]);
- ftension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
- ftension_params = tension_models[ftension_model].params;
- free_parsed_list(num_strings, strings);
- break;
- case 'a':
- ftension_params = optarg;
+ 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 'M':
+ 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 == 1 || num_strings == 2);
- if (num_strings == 1)
- utension_group = 0;
+ 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
- utension_group = atoi(strings[1]);
- utension_model = index_tension_model(n_tension_models, tension_models, strings[0]);
- utension_params = tension_models[utension_model].params;
+ 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 'A':
- utension_params = optarg;
+ 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':
- k_model = index_k_model(n_k_models, k_models, optarg);
- break;
- case 'K':
- k_models[k_model].params = optarg;
- break;
- case 'F':
- n = atoi(optarg);
- assert(n > 0);
- for (i=0; i<n; i++) {
- push(pDomain_list, generate_domain(DS_FOLDED,
- tension_models+ftension_model,
- ftension_group,
- ftension_params,
- tension_models+utension_model,
- utension_group,
- utension_params,
- k_models+k_model), 1);
- }
- *pNum_folded += n;
+ 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 'U':
- n = atoi(optarg);
- assert(n > 0);
- for (i=0; i<n; i++) {
- push(pDomain_list, generate_domain(DS_UNFOLDED,
- tension_models+ftension_model,
- ftension_group,
- ftension_params,
- tension_models+utension_model,
- utension_group,
- utension_params,
- k_models+k_model), 1);
- }
+ case 'D':
+ print_state_graph(stdout, *pState_list, *pTransition_list);
+ exit(0);
break;
case 'V':
flags = FLAG_OUTPUT_FULL_CURVE;
fprintf(stderr, "unrecognized option '%c'\n", optopt);
/* fall through to default case */
default:
- help(prog_name, *pP, *pDt_max, *pV, *pXstep, env, n_tension_models, tension_models, ftension_model, utension_model, n_k_models, k_models, k_model);
+ 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(*pDomain_list != NULL); /* no domains! */
+ 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;
}
-@
+@
<<index tension model>>=
int index_tension_model(int n_models, tension_model_getopt_t *models, char *name)
}
return i;
}
-@
+@
<<index k model>>=
int index_k_model(int n_models, k_model_getopt_t *models, char *name)
}
return i;
}
-@
+@
<<initialize model definition>>=
/* requires int num_param_args and char **param_args in the current scope
*/
#define INIT_MODEL(role, model, param_string, param_pointer) \
do { \
- parse_list_string(param_string, SEP, '{', '}', \
+ parse_list_string((param_string), SEP, '{', '}', \
&num_param_args, ¶m_args); \
- if (num_param_args != model->num_params) { \
+ if (num_param_args != (model)->num_params) { \
fprintf(stderr, \
"%s model %s expected %d params,\n", \
- role, model->name, model->num_params); \
+ (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); \
+ num_param_args, (param_string)); \
+ assert(num_param_args == (model)->num_params); \
} \
- if (model->creator) \
- param_pointer = (*model->creator)(param_args); \
+ if ((model)->creator) \
+ param_pointer = (*(model)->creator)(param_args); \
else \
param_pointer = NULL; \
free_parsed_list(num_param_args, param_args); \
} while (0);
-@
-
-<<generate domain>>=
-void *generate_domain(enum domain_state_t state,
- tension_model_getopt_t *folded_model,
- int folded_group,
- char *folded_param_string,
- tension_model_getopt_t *unfolded_model,
- int unfolded_group,
- char *unfolded_param_string,
- k_model_getopt_t *k_model)
-{
- void *folded_params, *unfolded_params, *k_params;
- int num_param_args; /* for INIT_MODEL() */
- char **param_args; /* for INIT_MODEL() */
+@
-#ifdef DEBUG
- fprintf(stderr, "generating %s ", state==DS_FOLDED ? "folded" : "unfolded");
- fprintf(stderr, "domain (k: \"%s\", \"%s\", f: \"%s\" (%p) group %d \"%s\", u: \"%s\" (%p) group %d \"%s\")\n",
- k_model->name, k_model->params,
- folded_model->name, folded_model->handler, folded_group, folded_param_string,
- unfolded_model->name, unfolded_model->handler, unfolded_group, unfolded_param_string);
-#endif
- INIT_MODEL("folded", folded_model, folded_param_string, folded_params);
- INIT_MODEL("unfolded", unfolded_model, unfolded_param_string, unfolded_params);
- INIT_MODEL("k", k_model, k_model->params, k_params);
- return create_domain(state,
- k_model->k, k_params, k_model->destructor,
- folded_model->handler, folded_group, folded_params, folded_model->destructor,
- unfolded_model->handler, unfolded_group, unfolded_params, unfolded_model->destructor
- );
+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
+
+<<safe strto*>>=
+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
The general layout of our simulation code is:
<<*>>=
+//#define DEBUG
<<license comment>>
<<includes>>
<<definitions>>
<<globals>>
<<functions>>
<<main program>>
-@
+@
We include [[math.h]], so don't forget to link to the libm with `-lm'.
<<includes>>=
#include <assert.h> /* assert() */
-#include <stdlib.h> /* malloc(), free(), rand() */
+#include <stdlib.h> /* malloc(), free(), rand(), strto*() */
#include <stdio.h> /* fprintf(), stdout */
#include <string.h> /* strlen, strtok() */
#include <math.h> /* exp(), M_PI, sqrt() */
#include <sys/types.h> /* pid_t (returned by getpid()) */
#include <unistd.h> /* getpid() (for seeding rand()), getopt() */
+#include <errno.h> /* errno, ERANGE (for safe_strto*()) */
#include "global.h"
#include "list.h"
#include "tension_balance.h"
#include "tension_model.h"
#include "parse.h"
#include "accel_k.h"
-@
+
+@
<<definitions>>=
<<version definition>>
<<string comparison definition and globals>>
<<initialize model definition>>
<<get options definitions>>
-<<domain definitions>>
-@
+<<state definitions>>
+<<transition definitions>>
+@
<<globals>>=
<<flag globals>>
<<model globals>>
-@
+@
<<functions>>=
-<<create/destroy domain>>
-<<destroy domain list>>
+<<create/destroy state>>
+<<destroy state list>>
+<<state index>>
+<<create/destroy transition>>
+<<destroy transition list>>
+<<print state graph>>
<<group functions>>
<<simulation functions>>
<<boilerplate functions>>
-@
+@
<<boilerplate functions>>=
<<setup>>
<<get option functions>>
-@
+@
<<setup>>=
void setup(void)
{
srand(getpid()*time(NULL)); /* seed rand() */
}
-@
+@
<<flag definitions>>=
/* in octal b/c of prefixed '0' */
-#define FLAG_OUTPUT_FULL_CURVE 01
-#define FLAG_OUTPUT_UNFOLDING_FORCES 02
+#define FLAG_OUTPUT_FULL_CURVE 01
+#define FLAG_OUTPUT_TRANSITION_FORCES 02
@
<<flag globals>>=
@
\subsection{Utilities}
-\label{app.utils}
+\label{sec.utils}
<<max/min definitions>>=
#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]].
<<string comparison definition and globals>>=
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
<<check relative error>>=
-(received-expected)/expected, max_err, #received, \
expected, received); \
} while(0)
-@
+
+@
<<happens>>=
int happens(double probability)
\subsection{Adaptive timesteps}
-\label{app.adaptive_dt}
-
-$F(x)$ increases with $x$, possibly exploding, as in the worm-like chain model,
-so basing the timestep on the the unfolding probability at the current tension
-is dangerous, and we need to search for a $dt$ for which
-$P(F(x+v*dt)) < P_\text{target}$.
-There are two cases to consider.
-In the most common, no domains have unfolded since the last step,
-and we expect the next step to be slightly shorter than the previous one.
-In the less common, domains did unfold in the last step,
-and we expect the next step to be considerably longer than the previous one.
+\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))<P_\text{max}$ (See Section
+\ref{sec.transition_rate} for a discussion of $P_N$). There are two
+cases to consider. In the most common, no domains have unfolded since
+the last step, and we expect the next step to be slightly shorter than
+the previous one. In the less common, domains did unfold in the last
+step, and we expect the next step to be considerably longer than the
+previous one.
<<search dt>>=
-double search_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers,
- list_t *domain_list,
+double search_dt(transition_t *transition,
+ list_t *state_list,
environment_t *env, double x,
- double target_prob, double max_dt, double v)
-{ /* Returns the timestep dt in seconds for the current folded domain.
- * Takes a list of tension handlers, the list of domains,
- * a pointer env to the environmental data, a starting separation x in m,
- * a target_prob between 0 and 1,
+ 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, dtCur, dtU, dtUCur, dtL, dt;
+ double F, k, P, dtCur, dtU, dtUCur, dtL, dt;
/* get upper bound using the current position */
- F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x, 0); /* BUG. repeated calculation */
+ 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(D(domain_list)->k, F, env, D(domain_list)->k_params);
+ k = accel_k(transition->k, F, env, transition->k_params);
//printf("x %g\tF %g\tk %g\n", x, F, k);
- dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */
+ dtU = P1(max_prob, N_INIT(transition)) / k; /* upper bound on dt */
if (dtU > max_dt) {
- //printf("overshot max_dt\n");
+ //printf("overshot max_dt\n");
dtU = max_dt;
}
- if (v == 0) /* discrete stepping, so force is temporatily constant */
+ if (v == 0) /* discrete stepping, so force is temporarily constant */
return dtU;
/* set a lower bound on dt too */
/* The dt determined above may produce illegitimate forces or ks.
* Reduce the upper bound until we have valid ks. */
dt = dtU;
- F = find_tension(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
+ 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(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
+ F = find_tension(state_list, env, x+v*dt, 0, 1);
}
//printf("Try for dt = %g (F = %g)\n", dt, F);
- k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
+ 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);
+ //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(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 0);
+ F = find_tension(state_list, env, x+v*dt, 0, 1);
//printf("Try for dt = %g (F = %g)\n", dt, F);
- k = accel_k(D(domain_list)->k, F, env, D(domain_list)->k_params);
- //printf("x %g\tF %g\tdt %g\tv dt %g\tk %g\n", x, F, dt, v*dt, k);
+ 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 */
- dtUCur = target_prob / k; /* safe timestep back from x+dtU */
+ /* 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(num_tension_handlers, tension_handlers, domain_list, env, x+v*dt, 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(D(domain_list)->k, F, env, D(domain_list)->k_params);
- dtCur = target_prob / k;
+ 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 */
+ dtU = dt; /* ... stepping out only dtCur would be safe */
dtUCur = dtCur;
} else
break; /* dtCur = dt */
}
return MAX(dtUCur, dtL);
}
-@
-To determine $dt$ for an array of potentially different folded domains,
-we need to find the maximum $dt$ that satisfies $k dt < P$ for all domains.
+@
+
+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)|<dF_\text{max}$, which limits errors due to our
+transition rate assumption that the force does not change appreciably
+over a single timestep.
<<determine dt>>=
<<search dt>>
-double determine_dt(int num_tension_handlers, one_dim_fn_t **tension_handlers, list_t *domain_list,
- environment_t *env, double x,
- double target_prob, double dt_max, double v)
+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 list of folded domains, target_prob between 0 and 1,
- * F in N, and T in K. */
- double dt=dt_max, new_dt;
- assert(target_prob > 0.0); assert(target_prob < 1.0);
+ * 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);
-
- while (domain_list != NULL) {
- if (D(domain_list)->state == DS_FOLDED) {
- new_dt = search_dt(num_tension_handlers, tension_handlers, domain_list, env, x, target_prob, dt, v);
+ 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);
}
- domain_list = domain_list->next;
+ transition_list = transition_list->next;
}
return dt;
}
-@
-\subsection{Domain data}
-
-Currently domains exist in two states, folded and unfolded, and the
-only allowed transitions are folded $\rightarrow$ unfolded. Of
-course, it wouldn't be too complicated to extent this to a multi-state
-system, with an array containing the domains group for each possible
-state, and a matrix of transition-rate-calculating functions.
-However, at this point such generality seems unnecessary at this
-point.
-<<domain definitions>>=
-enum domain_state_t {DS_FOLDED,
- DS_UNFOLDED
-};
+@
-typedef struct domain_struct {
- enum domain_state_t state;
- one_dim_fn_t *folded_handler;
- int folded_group;
- one_dim_fn_t *unfolded_handler;
- int unfolded_group;
- k_func_t *k; /* function returning unfolding rate */
- void *folded_params; /* pointer to folded parameters */
- void *unfolded_params; /* pointer to unfolded parameters */
- void *k_params; /* pointer to k parameters */
- destroy_data_func_t *destroy_folded;
- destroy_data_func_t *destroy_unfolded;
- destroy_data_func_t *destroy_k;
-} domain_t;
-
-/* get the domain data for the current list node */
-#define D(list) ((domain_t *)(list)->d)
-/* get the tension params for the current list node */
-#define D_TP(list) (((domain_t *)(list)->d)->state == DS_FOLDED \
- ? ((domain_t *)(list)->d)->folded_params \
- : ((domain_t *)(list)->d)->unfolded_params)
-@
-[[k]] is a pointer to the function determining the unfolding rate for a given tension.
-[[folded_params]] is a pointer to the parameters used by the function pointed to by [[k]].
-[[unfolded_params]] is a pointer to the parameters used by the group-appropriate handler function when determining the tension.
-The [[destroy_*]] pointers point to functions for freeing the memory [[*_params]].
-We store them with the domain data so that [[destroy_domain]] doesn't have to know which type of domain it's cleaning up after.
-
-[[create_]] and [[destroy_domain]] are simple wrappers around [[malloc]] and [[free]].
-<<create/destroy domain>>=
-domain_t *create_domain(enum domain_state_t state,
- k_func_t *k,
- void *k_params,
- destroy_data_func_t *destroy_k,
- one_dim_fn_t *folded_handler,
- int folded_group,
- void *folded_params,
- destroy_data_func_t *destroy_folded,
- one_dim_fn_t *unfolded_handler,
- int unfolded_group,
- void *unfolded_params,
- destroy_data_func_t *destroy_unfolded)
-{
- domain_t *ret = (domain_t *)malloc(sizeof(domain_t));
+\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}).
+<<state definitions>>=
+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]].
+<<create/destroy state>>=
+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);
- if (state == DS_FOLDED) {
- assert(k != NULL); /* the pointer points somewhere valid */
- assert(*k != NULL); /* and there is something useful there */
- } else
- assert(state == DS_UNFOLDED);
- ret->state = state;
- ret->folded_handler = folded_handler;
- ret->folded_group = folded_group;
- ret->unfolded_handler = unfolded_handler;
- ret->unfolded_group = unfolded_group;
- ret->k = k;
- ret->folded_params = folded_params;
- ret->unfolded_params = unfolded_params;
- ret->k_params = k_params;
- ret->destroy_folded = destroy_folded;
- ret->destroy_unfolded = destroy_unfolded;
- ret->destroy_k = destroy_k;
+ /* 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 domain %p. state %d\n", ret, ret->state);
+ 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_domain(domain_t *domain)
+void destroy_state(state_t *state)
{
- if (domain) {
- //printf("domain %p & %p\n", *domain, domain);
- if (domain->destroy_folded)
- (*domain->destroy_folded)(domain->folded_params);
- if (domain->destroy_unfolded)
- (*domain->destroy_unfolded)(domain->unfolded_params);
- if (domain->destroy_k)
- (*domain->destroy_k)(domain->k_params);
- free(domain);
+ 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);
}
}
+
@
-<<destroy domain list>>=
-void destroy_domain_list(list_t *domain_list)
+We'll be storing the states in a list (see Appendix \ref{sec.lists}),
+so we also define a convenience function for cleaning up.
+<<destroy state list>>=
+void destroy_state_list(list_t *state_list)
{
- domain_list = head(domain_list);
- while (domain_list != NULL)
- destroy_domain((domain_t *) pop(&domain_list));
+ state_list = head(state_list);
+ while (state_list != NULL)
+ destroy_state((state_t *) pop(&state_list));
}
-@
-\subsection{Domain model and group handling}
-
-<<group functions>>=
-<<get tension handler>>
-<<get group>>
-<<int comparison function>>
-<<find possible groups>>
-<<is group active>>
-<<get group list>>
-<<get active groups>>
-@
-
-<<get tension handler>>=
-one_dim_fn_t *get_tension_handler(domain_t *domain)
-{
- if (domain->state == DS_FOLDED)
- return domain->folded_handler;
- else {
- if (domain->state != DS_UNFOLDED)
- fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
- assert(domain->state == DS_UNFOLDED);
- return domain->unfolded_handler;
- }
-}
-@
+@
-<<get group>>=
-int get_group(domain_t *domain)
+We can also define a convenient way to get a state index from it's name.
+<<state index>>=
+int state_index(list_t *state_list, char *name)
{
- if (domain->state == DS_FOLDED)
- return domain->folded_group;
- else {
- if (domain->state != DS_UNFOLDED)
- fprintf(stderr, "oops, domain %p state %d\n", domain, domain->state);
- assert(domain->state == DS_UNFOLDED);
- return domain->unfolded_group;
+ 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;
}
-@
-
-We already know all possible tension classes, since they are all
-hardcoded into \prog. However, there may be any number of possible
-groups. We define a function [[find_possible_groups]] to search for
-possible (not neccessarily active) groups. It can search for
-subgroups of a single [[handler]], or by setting [[handler]] to
-[[NULL]], subgroups of \emph{any} handler. It returns a list of group
-integers.
-<<find possible groups>>=
-list_t *find_possible_groups(list_t *list, one_dim_fn_t *handler) {
- list_t *ret = NULL;
- list = head(list);
- while (list != NULL) {
- if (handler == NULL || get_tension_handler(D(list)) == handler) {
- push(&ret, &D(list)->folded_group, 1);
- push(&ret, &D(list)->unfolded_group, 1);
- }
- list = list->next;
- }
- if (ret == NULL)
- return ret; /* no groups with this handler, no processing needed */
+@
- /* sort the ret list, and remove duplicates */
- sort(&ret, &int_comparison);
- unique(&ret, &int_comparison);
+\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.
+<<transition definitions>>=
+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)
+
+<<PN definition>>
+<<P1 definition>>
+@ [[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]].
+<<create/destroy transition>>=
+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, "handler %p possible groups:", handler);
- list = head(ret);
- while (list != NULL) {
- fprintf(stderr, " %d", *((int *)list->d));
- list = list->next;
- }
- fprintf(stderr, "\n");
+ 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;
}
-@
-Search a [[list]] of domains, and determine whether a particular model
-class and group number combination is active.
-<<is group active>>=
-int is_group_active(list_t *list, one_dim_fn_t *handler, int group)
+void destroy_transition(transition_t *transition)
{
- list = head(list);
- while (list != NULL) {
- if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group)
- return 1;
- list = list->next;
+ 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);
}
- return 0;
}
-@
-Search a [[list]] of domains, and return all domains matching a
-particular model class and group number.
-<<get group list>>=
-list_t *get_group_list(list_t *list, one_dim_fn_t *handler, int group)
+@
+
+We'll be storing the transitions in a list (see Appendix
+\ref{sec.lists}), so we also define a convenience function for
+cleaning up.
+<<destroy transition list>>=
+void destroy_transition_list(list_t *transition_list)
{
- list_t *ret = NULL;
- list = head(list);
- while (list != NULL) {
- if (get_tension_handler(D(list)) == handler && get_group(D(list)) == group) {
- push(&ret, D_TP(list), 1); /* add a pointer to the appropriate tension parameters to our new list. */
-#ifdef DEBUG
- fprintf(stderr, "push domain %p %s tension parameters %p onto active group %p %d list %p\n", D(list), D(list)->state == DS_FOLDED ? "folded" : "unfolded", D_TP(list), handler, group, ret);
-#endif
- }
- list = list->next;
+ transition_list = head(transition_list);
+ while (transition_list != NULL)
+ destroy_transition((transition_t *) pop(&transition_list));
+}
+
+@
+
+\subsection{Graphviz output}
+\label{sec.graphviz_output}
+
+<<print state graph>>=
+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;
}
- return ret;
+ 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");
}
-@
-Because all the node data in lists returned by [[get_group_list]] is
-also in the main domain list, you shouldn't destroy the node data
-popped off when destroying the group lists. It will all get cleaned
-up when the main domain list is destroyed.
+@
+
+\subsection{Domain model and group handling}
+
+<<group functions>>=
+<<crosslink groups>>
+<<get active groups>>
+@
-Put all this together to scan through a list of domains and construct
-an array of all the active groups.
+Scan through a list of states and construct an array of all the
+active groups.
<<get active groups>>=
-void get_active_groups(list_t *list,
- int num_tension_handlers, one_dim_fn_t **tension_handlers,
- int *pNum_active_groups, one_dim_fn_t ***pPer_group_handlers, void ***pActive_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, **per_group_handlers=NULL;
- int i, num_active_groups, *group;
- list_t *possible_group_numbers, *active_group_list, *active_groups_list=NULL, *per_group_handlers_list=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 */
- for (i=0; i<num_tension_handlers; i++) {
- handler = tension_handlers[i];
- if (handler == NULL) continue; /* tensionless handler */
- possible_group_numbers = head(find_possible_groups(list, handler));
- while (possible_group_numbers != NULL) {
- group = (int *)pop(&possible_group_numbers);
- if (is_group_active(list, handler, *group) == 1) {
- active_group_list = get_group_list(list, handler, *group);
- push(&active_groups_list, active_group_list, 1);
- push(&per_group_handlers_list, handler, 1);
+ state_list = head(state_list);
#ifdef DEBUG
- fprintf(stderr, "active group %p for %p %d headed by tension parameters %p\n", active_group_list, handler, *group, head(active_group_list)->d);
- list_print(stderr, active_group_list, "active group");
+ 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
- /* allocate the array we'll be returning */
- num_active_groups = list_length(active_groups_list);
- active_groups = (void **) malloc(sizeof(void *)*num_active_groups);
- assert(active_groups != NULL);
+ 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; i<num_active_groups; i++) { /* we might want to recycle any group data structs allocated in the previous round */
+ per_group_data[i] = malloc(sizeof(tension_handler_data_t));
+ assert(per_group_data[i] != NULL);
+#ifdef DEBUG
+ fprintf(stderr, " %d,%p", i, per_group_data[i]);
+#endif
+ }
+#ifdef DEBUG
+ fprintf(stderr, "\n");
+#endif
- /* populate the array we'll be returning */
- active_groups_list = head(active_groups_list);
- for (i=0; i<num_active_groups; i++) {
- handler = pop(&per_group_handlers_list);
- per_group_handlers[i] = handler;
-
- active_groups[i] = (void *) malloc(sizeof(tension_handler_data_t));
- assert(active_groups[i] != NULL);
- active_group_list = pop(&active_groups_list);
- ((tension_handler_data_t *)active_groups[i])->group = active_group_list;
- ((tension_handler_data_t *)active_groups[i])->env = NULL;
- ((tension_handler_data_t *)active_groups[i])->persist = NULL;
+ /* populate the arrays we'll be returning */
+ active_states_list = head(active_states_list);
+ old_handler = NULL;
+ old_inverse_handler = NULL;
+ j = -1; /* j is the current active _group_ index */
+ while (active_states_list != NULL) {
+ state = (state_t *) pop(&active_states_list);
+ handler = state->tension_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; i<num_domains; i++) { /* add num_domains copies of the tension params */
+ push(&tdata->group_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;
}
- assert(active_groups_list == NULL);
- assert(per_group_handlers_list == NULL);
+ /* 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;
- *pActive_groups = active_groups;
+ *pPer_group_inverse_handlers = per_group_inverse_handlers;
+ *pPer_group_data = per_group_data;
+}
+@
+
+<<crosslink groups>>=
+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; i<n; i++) {
+ S(state_groups)->out_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; i<n; i++) {
+ S(state_groups)->grouped_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..."]]).
+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.
<<parse definitions>>
<<parse declarations>>
#endif /* PARSE */
-@
+@
<<parse module makefile lines>>=
NW_SAWSIM_MODS += parse
CHECK_BINS += check_parse
-check_parse_MODS =
-@
+check_parse_MODS =
+@
<<parse definitions>>=
#define SEP ',' /* argument separator character */
-@
+@
<<parse declarations>>=
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.
<<parse includes>>
#include "parse.h"
<<parse delimited list string functions>>
-@
+@
<<parse includes>>=
#include <assert.h> /* assert() */
#include <stdio.h> /* fprintf(), stdout *//*!!*/
#include <string.h> /* strlen() */
#include "parse.h"
-@
+@
\subsection{Parsing unit tests}
Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
<<check-parse.c>>=
+<<license comment>>
<<parse check includes>>
<<string comparison definition and globals>>
<<check relative error>>
<<parse test suite>>
<<main check program>>
-@
+@
<<parse check includes>>=
-#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
-#include <stdio.h> /* printf() */
-#include <assert.h> /* assert() */
-#include <string.h> /* strlen() */
+#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
+#include <stdio.h> /* printf() */
+#include <assert.h> /* assert() */
+#include <string.h> /* strlen() */
<<check includes>>
#include "parse.h"
-@
+@
<<parse test suite>>=
<<parse list string tests>>
<<parse suite definition>>
-@
+@
<<parse suite definition>>=
Suite *test_suite (void)
{
- Suite *s = suite_create ("k model");
+ Suite *s = suite_create ("parse");
<<parse list string test case defs>>
- <<parse list string test case add>>
+ <<parse list string test case adds>>
return s;
}
-@
+@
<<parse list string tests>>=
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
@
<<parse list string test case defs>>=
TCase *tc_parse_list_string = tcase_create("parse list string");
@
-<<parse list string test case add>>=
+<<parse list string test case adds>>=
//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);
@
<<functions>>
<<test suite>>
<<main check program>>
-@
+@
<<check includes>>=
#include <check.h>
-@
+@
<<check globals>>=
-@
+@
<<test suite>>=
-<<F tests>>
<<determine dt tests>>
<<happens tests>>
-<<does domain unfold tests>>
+<<does domain transition tests>>
<<randomly unfold domains tests>>
<<suite definition>>
-@
+@
<<suite definition>>=
Suite *test_suite (void)
{
Suite *s = suite_create ("sawsim");
- <<F test case defs>>
<<determine dt test case defs>>
<<happens test case defs>>
- <<does domain unfold test case defs>>
+ <<does domain transition test case defs>>
<<randomly unfold domains test case defs>>
- <<F test case add>>
- <<determine dt test case add>>
- <<happens test case add>>
- <<does domain unfold test case add>>
- <<randomly unfold domains test case add>>
+ <<determine dt test case adds>>
+ <<happens test case adds>>
+ <<does domain transition test case adds>>
+ <<randomly unfold domains test case adds>>
/*
tcase_add_checked_fixture(tc_strip_address,
*/
return s;
}
-@
+@
<<main check program>>=
int main(void)
srunner_free(sr);
return (number_failed == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
}
-@
-
-\subsection{F tests}
-<<F tests>>=
-<<worm-like chain tests>>
-@
-<<F test case defs>>=
-<<worm-like chain test case def>>
-@
-<<F test case add>>=
-<<worm-like chain test case add>>
-@
-
-<<worm-like chain tests>>=
-extern double wlc(double x, double T, double p, double L);
-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
-@
-<<worm-like chain test case def>>=
-TCase *tc_wlc = tcase_create("WLC");
-@
+@
-<<worm-like chain test case add>>=
-tcase_add_test(tc_wlc, test_wlc_at_zero);
-tcase_add_test(tc_wlc, test_wlc_at_half);
-suite_add_tcase(s, tc_wlc);
-@
\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...
<<determine dt tests>>=
double linear_k(double F, environment_t *env, void *params)
{
START_TEST(test_determine_dt_linear_k)
{
- environment_t env;
- double dt_max=3.0, Fnot=3.0;
+ state_t folded={0};
+ transition_t unfold={0};
+ environment_t env = {0};
double F[]={0,1,2,3,4,5,6};
- domain_t dom; /* use both parts at once for folded/unfolded */
+ 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;
-/*
- dom->next = dom->prev = NULL;
- dom->k_func_t = linear_k;
- dom->folded_params = &Fnot;
- dom->unfolded_params = !!!!!!!!!
- dom->destroy_folded = dom->destroy_unfolded = NULL;
+ 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(folded, unfolded, targ_p, dt_max, v, x, T)==1, NULL);
+ //fail_unless(determine_dt(state_list, transition_list, &env, x, max_p, max_dF, dt_max, v)==1, NULL);
}
-*/
}
END_TEST
@
<<determine dt test case defs>>=
TCase *tc_determine_dt = tcase_create("Determine dt");
@
-<<determine dt test case add>>=
+<<determine dt test case adds>>=
tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
suite_add_tcase(s, tc_determine_dt);
@
@
<<happens test case defs>>=
@
-<<happens test case add>>=
+<<happens test case adds>>=
@
-<<does domain unfold tests>>=
+<<does domain transition tests>>=
@
-<<does domain unfold test case defs>>=
+<<does domain transition test case defs>>=
@
-<<does domain unfold test case add>>=
+<<does domain transition test case adds>>=
@
<<randomly unfold domains tests>>=
@
<<randomly unfold domains test case defs>>=
@
-<<randomly unfold domains test case add>>=
+<<randomly unfold domains test case adds>>=
@
\section{Balancing group extensions}
-\label{app.tension_balance}
+\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
<<license comment>>
<<tension balance function declaration>>
#endif /* TENSION_BALANCE_H */
-@
+@
<<tension balance functions>>=
<<one dimensional bracket>>
<<one dimensional solve>>
<<x of x0>>
+<<group stiffness function>>
+<<chain stiffness function>>
+<<full chain stiffness function>>
<<tension balance function>>
-@
+@
<<tension balance module makefile lines>>=
NW_SAWSIM_MODS += tension_balance
CHECK_BINS += check_tension_balance
-check_tension_balance_MODS =
-@
+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
-non-empty group). $x(x_0)$ is also strictly monotonically increasing,
+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]]
<<tension balance function declaration>>=
double tension_balance(int num_tension_groups,
one_dim_fn_t **tension_handlers,
- void **params,
+ one_dim_fn_t **inverse_tension_handlers,
+ void **params, /* array of pointers to tension_handler_data_t */
double last_x,
- double x);
-@
+ double x,
+ double *stiffness);
+@
<<tension balance function>>=
double tension_balance(int num_tension_groups,
one_dim_fn_t **tension_handlers,
- void **params,
+ one_dim_fn_t **inverse_tension_handlers,
+ void **params, /* array of pointers to tension_handler_data_t */
double last_x,
- double x)
+ double x,
+ double *stiffness)
{
- static x_of_xo_data_t x_xo_data = {0,NULL,NULL,NULL};
+ 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(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<num_tension_groups; i++)
+ x_xo_data.xi[i] = last_x;
+ }
+
if (num_tension_groups == 1) { /* only one group, no balancing required */
xo = x;
+ x_xo_data.xi[0] = xo;
+#ifdef DEBUG
+ fprintf(stderr, "%s:\tonly one active group, no balancing required\n", __FUNCTION__);
+#endif
} else {
- //fprintf(stderr, "balancing for x = %g with ", x);
- //for(i=0; i<num_tension_groups; i++) fprintf(stderr, " %d: %g\t", i, x_xo_data.xi[i]);
- //fprintf(stderr, "\n");
- if (last_x == -1) { /* 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.handler_data = params;
- x_xo_data.xi = (double *) malloc(sizeof(double)*num_tension_groups);
- for (i=0; i<num_tension_groups; i++)
- x_xo_data.xi[i] = -1.0;
- }
+#ifdef DEBUG
+ fprintf(stderr, "%s:\tbalancing for x = %g with", __FUNCTION__, x);
+ for(i=0; i<num_tension_groups; i++) fprintf(stderr, "\tx[%d]=%g", i, x_xo_data.xi[i]);
+ fprintf(stderr, "\n");
+#endif
if (last_x == -1 || 1) { /* take a rough guess and attempt boundary conditions. */
if (x == 0) xo_guess = length_scale;
else xo_guess = x/num_tension_groups;
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");
- lb= x_xo_data.xi[0];
- ub= x_xo_data.xi[0];
+#endif
+ lb = x_xo_data.xi[0];
+ ub = x_xo_data.xi[0];
}
}
#ifdef DEBUG
}
#endif
}
+
F = (*tension_handlers[0])(xo, params[0]);
+
+ if (stiffness != NULL) {
+ *stiffness = chain_stiffness(&x_xo_data, min_relx);
+ if (*stiffness == 0.0) { /* re-calculate based on full chain */
+ *stiffness = full_chain_stiffness(num_tension_groups, tension_handlers,
+ inverse_tension_handlers, params,
+ last_x, x, min_relx, F);
+ }
+ }
return F;
}
-@
+
+@
<<tension balance internal definitions>>=
<<x of x0 definitions>>
-@
+@
<<x of x0 definitions>>=
typedef struct x_of_xo_data_struct {
int n_groups;
- one_dim_fn_t **tension_handler; /* array of fn pointers */
- void **handler_data; /* array of void* pointers */
- double *xi;
+ 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;
-@
+@
<<x of x0>>=
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;
+ 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: finding x(x0=%g). F0(x0) = %g\n", __FUNCTION__, xo, F);
+ 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->xi[i] == -1) xi_guess = length_scale; /* HACK */
- else xi_guess = data->xi[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: searching for proper x[%d] for F[%d](x[%d]) = %g N\n", __FUNCTION__, i, i, i, F);
+ 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);
+ 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: found F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
+ fprintf(stderr, "%s:\tfound F[%d](x[%d]=%g) = %g N\n", __FUNCTION__, i, i, xi, F);
#endif
- data->xi[i] = xi;
- x += xi;
- if (data->xi) data->xi[i] = xi;
}
#ifdef DEBUG
- fprintf(stderr, "%s: found x(x0=%g) = %g\n", __FUNCTION__, xo, x);
+ 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]].
+<<group stiffness function>>=
+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.
+<<chain stiffness function>>=
+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).
+<<full chain stiffness function>>=
+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
<<one dimensional function definition>>=
/* equivalent to gsl_func_t */
typedef double one_dim_fn_t(double x, void *params);
-@
+@
<<one dimensional solve declaration>>=
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 min_relx, double min_rely, int max_steps,
double *pYx);
-@
+@
<<one dimensional solve>>=
/* 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 min_relx, double min_rely, int max_steps,
double *pYx)
{
double yx, ylb, yub, x;
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 bounded */
+ 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; }
#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.
<<one dimensional bracket declaration>>=
void oneD_bracket(one_dim_fn_t *f, void *params, double y, double xguess, double *lb, double *ub);
-@
+@
<<one dimensional bracket>>=
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);
}
//fprintf(stdout, "ub %g, %g lb\n", *ub, *lb);
}
-@
+@
\subsection{Balancing implementation}
<<tension balance internal definitions>>
<<tension balance functions>>
-@
+@
<<tension balance includes>>=
#include <assert.h> /* assert() */
#include <math.h> /* HUGE_VAL, macro constant, so don't need to link to math */
#include <stdio.h> /* 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}.
<<check-tension-balance.c>>=
+<<license comment>>
<<tension balance check includes>>
<<tension balance test suite>>
<<main check program>>
-@
+@
<<tension balance check includes>>=
#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
<<check includes>>
#include "global.h"
#include "tension_balance.h"
-@
+@
<<tension balance test suite>>=
<<tension balance function tests>>
<<tension balance suite definition>>
-@
+@
<<tension balance suite definition>>=
Suite *test_suite (void)
<<tension balance function test case adds>>
return s;
}
-@
+@
<<tension balance function tests>>=
<<check relative error>>
{
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, data, last_x, x);
+ 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{app.math_hooke}.
+The analytic solution for a general number of springs is given in Appendix \ref{sec.math_hooke}.
<<tension balance function tests>>=
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, data, last_x, x);
+ F = tension_balance(2, handlers, inverse_handlers, data, last_x, x, NULL);
x1e = x*k2/(k1+k2);
Fe = k1*x1e;
x2e = x1e*k1/k2;
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.
<<tension balance function test case defs>>=
TCase *tc_tbfunc = tcase_create("tension balance function");
-@
+@
<<tension balance function test case adds>>=
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.
<<list definitions>>
<<list declarations>>
#endif /* LIST_H */
-@
+@
<<list declarations>>=
<<head and tail declarations>>
<<list to array declaration>>
<<move declaration>>
<<sort declaration>>
+<<list index declaration>>
+<<list element declaration>>
<<unique declaration>>
<<list print declaration>>
-@
+@
<<list functions>>=
<<create and destroy node>>
<<list to array>>
<<move>>
<<sort>>
+<<list index>>
+<<list element>>
<<unique>>
<<list print>>
-@
+@
<<list module makefile lines>>=
NW_SAWSIM_MODS += list
CHECK_BINS += check_list
-check_list_MODS =
-@
+check_list_MODS =
+@
<<list definitions>>=
typedef struct list_struct {
} list_t;
<<comparison function typedef>>
-@
+@
[[head]] and [[tail]] return pointers to the head and tail nodes of the list:
<<head and tail declarations>>=
list_t *head(list_t *list);
list_t *tail(list_t *list);
-@
+@
<<head and tail>>=
list_t *head(list_t *list)
{
}
return list;
}
-@
+@
<<list length declaration>>=
int list_length(list_t *list);
-@
+@
<<list length>>=
int list_length(list_t *list)
{
}
return i;
}
-@
+@
[[push]] inserts a node relative to the current position in the list
without changing the current position.
otherwise we insert below the current node (push(a>B>c,d,1) -> a>B>d>c).
<<push declaration>>=
void push(list_t **pList, void *data, int below);
-@
+@
<<push>>=
void push(list_t **pList, void *data, int below)
-{
+{
list_t *list, *new_node;
assert(pList != NULL);
list = *pList;
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.
<<pop declaration>>=
void *pop(list_t **pList);
-@
+@
<<pop>>=
void *pop(list_t **pList)
{
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.
<<list destroy declaration>>=
void list_destroy(list_t **pList, destroy_data_func_t *destroy);
-@
+@
<<list destroy>>=
void list_destroy(list_t **pList, destroy_data_func_t *destroy)
{
destroy(data);
}
}
-@
+@
[[list_to_array]] converts a list to an array of pointers, preserving order.
<<list to array declaration>>=
void list_to_array(list_t *list, int *length, void ***array);
-@
+@
<<list to array>>=
void list_to_array(list_t *list, int *pLength, void ***pArray)
{
*pLength = length;
*pArray = array;
}
-@
+@
[[move]] moves the current element along the list in either direction.
<<move declaration>>=
void move(list_t **pList, int downstream);
-@
+@
<<move>>=
void move(list_t **pList, int downstream)
{
list = flipper;
if (downstream == 0) {
push(&list, data, 0); /* b>A>c>d */
- list = list->prev; /* 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.
<<comparison function typedef>>=
typedef int (comparison_fn_t)(void *A, void *B);
-@
+@
For example
<<int comparison function>>=
int int_comparison(void *A, void *B)
else if (a == b) return 0;
else return -1;
}
-@
+
+@
<<sort declaration>>=
void sort(list_t **pList, comparison_fn_t *comp);
-@
+@
Here we use the bubble sort algorithm.
<<sort>>=
void sort(list_t **pList, comparison_fn_t *comp)
}
*pList = list;
}
-@
+
+@
+
+[[list_index]] finds the location of [[data]] in [[list]]. Returns
+[[-1]] if [[data]] is not in [[list]].
+<<list index declaration>>=
+int list_index(list_t *list, void *data);
+@
+
+<<list index>>=
+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.
+<<list element declaration>>=
+void *list_element(list_t *list, int i);
+@
+<<list element>>=
+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
drops any non-unique data without freeing it.
<<unique declaration>>=
void unique(list_t **pList, comparison_fn_t *comp);
-@
+@
<<unique>>=
void unique(list_t **pList, comparison_fn_t *comp)
{
}
*pList = list;
}
-@
+@
[[list_print]] prints a list to a [[FILE *]] stream.
<<list print declaration>>=
void list_print(FILE *file, list_t *list, char *list_name);
-@
+@
<<list print>>=
void list_print(FILE *file, list_t *list, char *list_name)
{
}
fprintf(file, "\n");
}
-@
+@
[[create_]] and [[destroy_node]] are simple wrappers around [[malloc]] and [[free]].
<<create and destroy node>>=
<<list includes>>
#include "list.h"
<<list functions>>
-@
+@
<<list includes>>=
#include <assert.h> /* assert() */
#include <stdlib.h> /* malloc(), free(), rand() */
#include <stdio.h> /* 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}.
<<check-list.c>>=
+<<license comment>>
<<list check includes>>
<<list test suite>>
<<main check program>>
-@
+@
<<list check includes>>=
#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
<<check includes>>
#include "global.h"
#include "list.h"
-@
+@
<<list test suite>>=
<<push tests>>
<<pop tests>>
<<list suite definition>>
-@
+@
<<list suite definition>>=
Suite *test_suite (void)
<<pop test case adds>>
return s;
}
-@
+@
<<push tests>>=
START_TEST(test_push)
}
}
END_TEST
-@
+@
<<push test case defs>>=
TCase *tc_push = tcase_create("push");
-@
+@
<<push test case adds>>=
tcase_add_test(tc_push, test_push);
suite_add_tcase(s, tc_push);
-@
+@
<<pop tests>>=
-@
+@
<<pop test case defs>>=
-@
+@
<<pop test case adds>>=
-@
+@
\section{Function string evaluation}
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.
+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.
<<string eval function declaration>>
<<string eval teardown declaration>>
#endif /* STRING_EVAL_H */
-@
+@
<<string eval module makefile lines>>=
NW_SAWSIM_MODS += string_eval
CHECK_BINS += check_string_eval
-check_string_eval_MODS =
-@
+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.), \\
//#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 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]].
#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.
<<string eval setup declaration>>=
extern void string_eval_setup(FILE **pIn, FILE **pOut);
-@
+@
<<string eval setup definition>>=
void string_eval_setup(FILE **pIn, FILE **pOut)
{
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 */
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.
+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.
As an added touch, we redirect Python's [[stderr]] before closing the pipes, otherwise it prints a blank line when it exits.
<<string eval teardown declaration>>=
extern void string_eval_teardown(FILE **pIn, FILE **pOut);
-@
+@
<<string eval teardown definition>>=
void string_eval_teardown(FILE **pIn, FILE **pOut)
}
*/
}
-@
+@
The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
\subsection{String evaluation implementation}
#include "string_eval.h"
<<string eval internal definitions>>
<<string eval functions>>
-@
+@
<<string eval includes>>=
#include <assert.h> /* assert() */
#include <sys/types.h> /* pid_t */
#include <unistd.h> /* pipe(), fork(), execvp(), alarm() */
#include <wait.h> /* wait() */
-@
+@
<<string eval internal definitions>>=
<<string eval subprocess definitions>>
<<string eval pipe definitions>>
-@
+@
<<string eval functions>>=
<<string eval setup definition>>
<<string eval function definition>>
<<string eval teardown definition>>
-@
+@
\subsection{String evaluation unit tests}
Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
<<check-string-eval.c>>=
+<<license comment>>
<<string eval check includes>>
<<string comparison definition and globals>>
<<string eval test suite>>
<<main check program>>
-@
+@
<<string eval check includes>>=
-#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
-#include <stdio.h> /* printf() */
-#include <assert.h> /* assert() */
-#include <string.h> /* strlen() */
-#include <signal.h> /* SIGKILL */
+#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE */
+#include <stdio.h> /* printf() */
+#include <assert.h> /* assert() */
+#include <string.h> /* strlen() */
+#include <signal.h> /* SIGKILL */
<<check includes>>
#include "string_eval.h"
-@
+@
<<string eval test suite>>=
<<string eval tests>>
<<string eval suite definition>>
-@
+@
<<string eval suite definition>>=
Suite *test_suite (void)
Suite *s = suite_create ("string eval");
<<string eval test case defs>>
- <<string eval test case add>>
+ <<string eval test case adds>>
return s;
}
-@
+@
<<string eval tests>>=
START_TEST(test_setup_teardown)
FILE *in, *out;
char input[80], output[80]={};
string_eval_setup(&in, &out);
- sprintf(input, "print ABCDefg\n");
- string_eval(in, out, input, 80, output);
+ 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
char input[80], output[80]={};
string_eval_setup(&in, &out);
sprintf(input, "print 3+4*5\n");
- string_eval(in, out, input, 80, output);
+ string_eval(in, out, input, 80, output);
fail_unless(STRMATCH(output,"23\n"), NULL);
string_eval_teardown(&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);
+ 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);
+ string_eval(in, out, input, 80, output);
fail_unless(STRMATCH(output,"5.0\n"), NULL);
string_eval_teardown(&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);
+ string_eval(in, out, input, 80, output);
fail_unless(STRMATCH(output,"9\n"), NULL);
string_eval_teardown(&in, &out);
}
<<string eval test case defs>>=
TCase *tc_string_eval = tcase_create("string_eval");
@
-<<string eval test case add>>=
+<<string eval test case adds>>=
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);
\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]]),
+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.
<<accel-k.h>>=
#ifndef ACCEL_K_H
#define ACCEL_K_H
+<<license comment>>
double accel_k(k_func_t *k, double F, environment_t *env, void *params);
void free_accels();
#endif /* ACCEL_K_H */
-@
+@
<<accel k module makefile lines>>=
NW_SAWSIM_MODS += accel_k
-#CHECK_BINS += check_accel_k
-check_accel_k_MODS =
-@
+CHECK_BINS += check_accel_k
+check_accel_k_MODS = interp k_model parse string_eval tavl
+@
<<accel-k.c>>=
+<<license comment>>
#include <assert.h> /* assert() */
#include <stdlib.h> /* realloc(), free(), NULL */
#include "global.h" /* environment_t */
{
int i=num_accels;
accels = (accel_k_t *)realloc(accels, sizeof(accel_k_t) * ++num_accels);
- assert(accels != NULL);
+ assert(accels != NULL);
accels[i].itable = interp_table_allocate(&k_wrap, &k_tol);
accels[i].k = k;
accels[i].env = env;
* WARNING: we're assuming param and env are the same. */
return interp_table_eval(accels[i].itable, accels+i, F);
}
- }
+ }
/* set up a new acceleration environment */
i = add_accel_k(k, env, params);
return interp_table_eval(accels[i].itable, accels+i, F);
}
-@
+@
+
+\subsection{Unit tests}
+
+Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
+<<check-accel-k.c>>=
+<<license comment>>
+<<accel k check includes>>
+<<check relative error>>
+<<model globals>>
+<<accel k test suite>>
+<<main check program>>
+@
+
+<<accel k check includes>>=
+#include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
+#include <stdio.h> /* sprintf() */
+#include <assert.h> /* assert() */
+#include <math.h> /* exp() */
+#include <errno.h> /* errno, ERANGE */
+<<check includes>>
+#include "global.h"
+#include "k_model.h"
+#include "accel_k.h"
+
+@
+
+<<accel k test suite>>=
+<<accel k suite tests>>
+<<accel k suite definition>>
+@
+
+<<accel k suite definition>>=
+Suite *test_suite (void)
+{
+ Suite *s = suite_create ("accelerated k model");
+ <<accel k test case defs>>
+
+ <<accel k test case adds>>
+ return s;
+}
+@
+
+<<accel k test case defs>>=
+TCase *tc_accel_k = tcase_create("accelerated k model");
+@
+
+<<accel k test case adds>>=
+tcase_add_test(tc_accel_k, test_accel_k_accuracy);
+suite_add_tcase(s, tc_accel_k);
+@
+
+<<accel k suite tests>>=
+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}
<<hooke tension model declarations>>
<<worm-like chain tension model declarations>>
<<freely-jointed chain tension model declarations>>
+<<piston tension model declarations>>
<<find tension definitions>>
<<tension model global declarations>>
#endif /* TENSION_MODEL_H */
-@
+@
<<tension model module makefile lines>>=
NW_SAWSIM_MODS += tension_model
-#CHECK_BINS += check_tension_model
-check_tension_model_MODS =
-@
+CHECK_BINS += check_tension_model
+check_tension_model_MODS = list tension_balance
+@
<<tension model utils makefile lines>>=
TENSION_MODEL_MODS = tension_model parse list tension_balance
TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.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)
- gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
+ $(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
For unstretchable domains.
<<null tension model getopt>>=
-{NULL, "null", "an unstretchable domain", 0, NULL, NULL, NULL, NULL}
-@
+{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.
+
<<constant tension functions>>=
<<constant tension function>>
<<constant tension parameter create/destroy functions>>
-@
+@
<<constant tension model declarations>>=
<<constant tension function declaration>>
<<constant tension parameter create/destroy function declarations>>
<<constant tension model global declarations>>
-@
-
-An infinitely stretchable domain providing a constant tension.
+@
<<constant tension function declaration>>=
extern double const_tension_handler(double x, void *pdata);
-@
+@
<<constant tension function>>=
double const_tension_handler(double x, void *pdata)
{
tension_handler_data_t *data = (tension_handler_data_t *)pdata;
- list_t *list = data->group;
+ list_t *list = data->group_tension_params;
double F;
- assert (x >= 0.0);
+ assert(x >= 0.0);
list = head(list);
assert(list != NULL); /* empty active group?! */
F = ((const_tension_param_t *)list->d)->F;
return F;
}
-@
+@
<<constant tension structure>>=
typedef struct const_tension_param_struct {
double F; /* tension (force) in N */
} const_tension_param_t;
-@
+@
<<constant tension parameter create/destroy function declarations>>=
extern void *string_create_const_tension_param_t(char **param_strings);
extern void destroy_const_tension_param_t(void *p);
-@
+@
<<constant tension parameter create/destroy functions>>=
const_tension_param_t *create_const_tension_param_t(double F)
void *string_create_const_tension_param_t(char **param_strings)
{
- return create_const_tension_param_t(atof(param_strings[0]));
+ return create_const_tension_param_t(
+ safe_strtod(param_strings[0],__FUNCTION__));
}
void destroy_const_tension_param_t(void *p)
extern int num_const_tension_params;
extern char *const_tension_param_descriptions[];
extern char const_tension_param_string[];
-@
+@
<<constant tension model globals>>=
int num_const_tension_params = 1;
char *const_tension_param_descriptions[] = {"tension F, N"};
char const_tension_param_string[] = "0";
-@
+@
<<constant tension model getopt>>=
-{&const_tension_handler, "const", "an infinitely stretchable domain with constant tension", 1, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t}
-@
+{&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}
<<hooke functions>>=
-<<hooke function>>
+<<internal hooke functions>>
+<<hooke handler>>
+<<inverse hooke handler>>
<<hooke parameter create/destroy functions>>
-@
+@
<<hooke tension model declarations>>=
<<hooke tension function declaration>>
<<hooke tension parameter create/destroy function declarations>>
<<hooke tension model global declarations>>
-@
+@
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{app.math_hooke}.
+For a simple proof, see Appendix \ref{sec.math_hooke}.
<<hooke tension function declaration>>=
extern double hooke_handler(double x, void *pdata);
-@
+extern double inverse_hooke_handler(double F, void *pdata);
+@
-<<hooke function>>=
-double hooke_handler(double x, void *pdata)
-{
- tension_handler_data_t *data = (tension_handler_data_t *)pdata;
- list_t *list = data->group;
+First we define a function that computes the effective spring constant
+(stored in a single [[hooke_param_t]]) for the entire group.
+<<internal hooke functions>>=
+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;
- assert (x >= 0.0);
list = head(list);
assert(list != NULL); /* empty active group?! */
while (list != NULL) {
}
k = 1.0 / k;
#ifdef DEBUG
- fprintf(stderr, "%s: effective Hooke spring %g N/m at %g m -> F = %g N\n",
- __FUNCTION__, k, x, k*x);
+ 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
- return k*x;
+ param->k = k;
+}
+
+@
+
+Everyone knows Hooke's law: $F=kx$.
+<<hooke handler>>=
+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$.
+<<inverse hooke handler>>=
+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$.
<<hooke structure>>=
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.
<<hooke tension parameter create/destroy function declarations>>=
extern void *string_create_hooke_param_t(char **param_strings);
extern void destroy_hooke_param_t(void *p);
-@
+
+@
<<hooke parameter create/destroy functions>>=
hooke_param_t *create_hooke_param_t(double k)
void *string_create_hooke_param_t(char **param_strings)
{
- return create_hooke_param_t(atof(param_strings[0]));
+ return create_hooke_param_t(safe_strtod(param_strings[0], __FUNCTION__));
}
void destroy_hooke_param_t(void *p)
extern int num_hooke_params;
extern char *hooke_param_descriptions[];
extern char hooke_param_string[];
-@
+@
<<hooke tension model globals>>=
int num_hooke_params = 1;
char *hooke_param_descriptions[] = {"spring constant k, N/m"};
char hooke_param_string[]="0.05";
-@
+@
<<hooke tension model getopt>>=
-{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}
-@
+{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.
<<worm-like chain functions>>=
-<<worm-like chain function>>
+<<internal worm-like chain functions>>
<<worm-like chain handler>>
+<<inverse worm-like chain handler>>
<<worm-like chain create/destroy functions>>
-@
+@
<<worm-like chain tension model declarations>>=
+
<<worm-like chain handler declaration>>
+<<inverse worm-like chain handler declaration>>
<<worm-like chain create/destroy function declarations>>
<<worm-like chain tension model global declarations>>
-@
+@
+
+<<internal worm-like chain functions>>=
+<<worm-like chain function>>
+<<inverse worm-like chain function>>
+<<worm-like chain parameter grouper>>
+@
The combination of all unfolded domains is modeled as a worm like chain,
with the force $F_{\text{WLC}}$ approximately given by
$p$ is the persistence length, and
$L$ is the total contour length of all unfolded domains\citep{bustamante94,marko95,hatfield99}.
<<worm-like chain function>>=
-double wlc(double x, double T, double p, double L)
-{/* N m K m m */
+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);
// 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 );
}
-@
-This model assumes that each unfolded domain has the same persistence length.
-<<worm-like chain handler declaration>>=
-extern double wlc_handler(double x, void *pdata);
-@
+@
-<<worm-like chain handler>>=
-double wlc_handler(double x, void *pdata)
-{
- tension_handler_data_t *data = (tension_handler_data_t *)pdata;
- list_t *list = data->group;
- double p, L=0.0;
+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.
+<<inverse worm-like chain function>>=
+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.
+<<worm-like chain parameter grouper>>=
+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?! */
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
- return wlc(x, data->env->T, p, L);
+ param->p = p;
+ param->L = L;
+}
+
+@
+
+<<worm-like chain handler declaration>>=
+extern double wlc_handler(double x, void *pdata);
+@
+
+This model requires that each unfolded domain in the group have the
+same persistence length.
+<<worm-like chain handler>>=
+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);
}
-@
+
+@
+
+<<inverse worm-like chain handler declaration>>=
+extern double inverse_wlc_handler(double F, void *pdata);
+@
+
+<<inverse worm-like chain handler>>=
+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);
+}
+
+@
<<worm-like chain structure>>=
typedef struct wlc_param_struct {
double p; /* WLC persistence length (m) */
double L; /* WLC contour length (m) */
} wlc_param_t;
-@
+@
<<worm-like chain create/destroy function declarations>>=
extern void *string_create_wlc_param_t(char **param_strings);
extern void destroy_wlc_param_t(void *p /* wlc_param_t * */);
-@
+@
<<worm-like chain create/destroy functions>>=
wlc_param_t *create_wlc_param_t(double p, double L)
void *string_create_wlc_param_t(char **param_strings)
{
- return create_wlc_param_t(atof(param_strings[0]), atof(param_strings[1]));
+ 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.
extern int num_wlc_params;
extern char *wlc_param_descriptions[];
extern char wlc_param_string[];
-@
+@
<<worm-like chain tension model globals>>=
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";
-@
+@
<<worm-like chain tension model getopt>>=
-{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}
-@
+{&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}
<<freely-jointed chain function>>
<<freely-jointed chain handler>>
<<freely-jointed chain create/destroy functions>>
-@
+@
<<freely-jointed chain tension model declarations>>=
<<freely-jointed chain handler declaration>>
<<freely-jointed chain create/destroy function declarations>>
<<freely-jointed chain tension model global declarations>>
-@
+@
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
$$
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}.
<<freely-jointed chain function>>=
-double langevin(double x, void *params)
+static double langevin(double x, void *params)
{
if (x==0) return 0;
return 1.0/tanh(x) - 1/x;
}
-<<one dimensional bracket declaration>>
-<<one dimensional solve declaration>>
-double inv_langevin(double x)
+
+static double inv_langevin(double x)
{
double lb=0.0, ub=1.0, ret;
oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
return ret;
}
-double fjc(double x, double T, double l, int N)
+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);
}
-@
+
+@
<<freely-jointed chain handler declaration>>=
extern double fjc_handler(double x, void *pdata);
-@
+@
<<freely-jointed chain handler>>=
double fjc_handler(double x, void *pdata)
{
tension_handler_data_t *data = (tension_handler_data_t *)pdata;
- list_t *list = data->group;
+ list_t *list = data->group_tension_params;
double l;
int N=0;
#endif
return fjc(x, data->env->T, l, N);
}
-@
+@
<<freely-jointed chain structure>>=
typedef struct fjc_param_struct {
double l; /* FJC link length (m) */
int N; /* FJC number of links */
} fjc_param_t;
-@
+@
<<freely-jointed chain create/destroy function declarations>>=
extern void *string_create_fjc_param_t(char **param_strings);
extern void destroy_fjc_param_t(void *p);
-@
+@
<<freely-jointed chain create/destroy functions>>=
fjc_param_t *create_fjc_param_t(double l, double N)
void *string_create_fjc_param_t(char **param_strings)
{
- return create_fjc_param_t(atof(param_strings[0]), atof(param_strings[1]));
+ 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);
}
-@
+@
<<freely-jointed chain tension model global declarations>>=
extern int num_fjc_params;
extern char *fjc_param_descriptions[];
extern char fjc_param_string[];
-@
+@
<<freely-jointed chain tension model globals>>=
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";
-@
+@
<<freely-jointed chain tension model getopt>>=
-{fjc_handler, "fjc", "a freely-jointed chain [F=kbT/l L^-1(x/L)]", num_fjc_params, fjc_param_descriptions, (char *)fjc_param_string, &string_create_fjc_param_t, &destroy_fjc_param_t}
-@
+{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 $x<L$}, \\
+ \infty & \text{if $x>L$}.
+ \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<L$, because any reactions that occur at $F=0$
+(e.g. refolding) will have more time to occur.
+
+Because the tension at $x=L$ is undefined, the user must make sure
+domains with this tension model are never the initial active group,
+because it would break [[tension_balance()]] in [[find_tension()]]
+(see Section \ref{sec.tension_balance}).
+
+<<piston tension functions>>=
+<<piston tension parameter grouper>>
+<<piston tension handler>>
+<<inverse piston tension handler>>
+<<piston tension parameter create/destroy functions>>
+@
+
+<<piston tension model declarations>>=
+<<piston tension handler declaration>>
+<<inverse piston tension handler declaration>>
+<<piston tension parameter create/destroy function declarations>>
+<<piston tension model global declarations>>
+
+@
+
+First we define a function that computes the effective piston parameters
+(stored in a single [[piston_tension_param_t]]) for the entire group.
+<<piston tension parameter grouper>>=
+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;
+}
+
+<<piston tension handler declaration>>=
+extern double piston_tension_handler(double x, void *pdata);
+@
+<<piston tension handler>>=
+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.
+
+<<inverse piston tension handler declaration>>=
+extern double inverse_piston_tension_handler(double x, void *pdata);
+@
+<<inverse piston tension handler>>=
+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;
+}
+
+@
+
+<<piston tension structure>>=
+typedef struct piston_tension_param_struct {
+ double L; /* length of piston in m */
+} piston_tension_param_t;
+
+@
+
+<<piston tension parameter create/destroy function declarations>>=
+extern void *string_create_piston_tension_param_t(char **param_strings);
+extern void destroy_piston_tension_param_t(void *p);
+
+@
+
+<<piston tension parameter create/destroy functions>>=
+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);
+}
+
+@
+
+<<piston tension model global declarations>>=
+extern int num_piston_tension_params;
+extern char *piston_tension_param_descriptions[];
+extern char piston_tension_param_string[];
+
+@
+
+<<piston tension model globals>>=
+int num_piston_tension_params = 1;
+char *piston_tension_param_descriptions[] = {"length L, m"};
+char piston_tension_param_string[] = "0";
+
+@
+
+<<piston tension model getopt>>=
+{&piston_tension_handler, &inverse_piston_tension_handler, "piston", "a domain that slides without tension for x<L, but has infinite tension for x>L.", 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}
<<tension model internal definitions>>
<<tension model globals>>
<<tension model functions>>
-@
+@
<<tension model includes>>=
#include <assert.h> /* assert() */
-#include <stdlib.h> /* NULL */
+#include <stdlib.h> /* NULL, strto*() */
#include <math.h> /* HUGE_VAL, sqrt(), exp() */
#include <stdio.h> /* fprintf(), stdout */
+#include <errno.h> /* errno, ERANGE */
#include "global.h"
#include "list.h"
#include "tension_balance.h" /* oneD_*() */
-@
+@
<<tension model internal definitions>>=
<<constant tension structure>>
<<hooke structure>>
<<worm-like chain structure>>
<<freely-jointed chain structure>>
-@
+<<piston tension structure>>
+@
<<tension model globals>>=
<<tension handler array global>>
<<hooke tension model globals>>
<<worm-like chain tension model globals>>
<<freely-jointed chain tension model globals>>
-@
+<<piston tension model globals>>
+@
<<tension model functions>>=
+<<safe strto*>>
<<constant tension functions>>
<<hooke functions>>
<<worm-like chain functions>>
<<freely-jointed chain functions>>
-@
-
+<<piston tension functions>>
+@
\subsection{Utilities}
<<tension model utility getopt functions>>
<<setup>>
<<tension model utility main>>
-@
+@
<<tension model utility main>>=
int main(int argc, char **argv)
tension_handler_data_t tdata;
int num_param_args; /* for INIT_MODEL() */
char **param_args; /* for INIT_MODEL() */
- get_options(argc, argv, &env, NUM_TENSION_MODELS, tension_models, &model, &flags);
+ 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 = NULL;
- push(&tdata.group, params, 1);
+ tdata.group_tension_params = NULL;
+ push(&tdata.group_tension_params, params, 1);
tdata.env = &env;
tdata.persist = NULL;
if (model->handler == NULL) {
exit(0);
}
{
- double dx=1e-10, x=0, F=0;
- printf("#F (N)\tk (%% pop. per s)\n");
- while (F >= 0 && F < 1e5 && x < 1e-6) {
+ 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);
- x += dx;
+ }
+ 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);
+ params = pop(&tdata.group_tension_params);
if (model->destructor)
(*model->destructor)(params);
return 0;
}
-@
+@
+
+<<tension model utility includes>>=
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h> /* strlen() */
+#include <assert.h> /* assert() */
+#include <errno.h> /* errno, ERANGE */
+#include "global.h"
+#include "parse.h"
+#include "list.h"
+#include "tension_model.h"
+@
+
+<<tension model utility definitions>>=
+<<version definition>>
+#define VFLAG 1 /* verbose */
+<<string comparison definition and globals>>
+<<tension model getopt definitions>>
+<<initialize model definition>>
+
+@
+
+<<tension model utility getopt functions>>=
+<<safe strto*>>
+<<index tension model>>
+<<tension model utility help>>
+<<tension model utility get options>>
+@
+
+<<tension model utility help>>=
+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<n_tension_models; i++) {
+ printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
+ for (j=0; j < tension_models[i].num_params ; j++)
+ printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
+ printf("\t\tdefaults: %s\n", tension_models[i].params);
+ }
+}
+@
+
+<<tension model utility get options>>=
+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}.
+<<check-tension-model.c>>=
+<<license comment>>
+<<tension model check includes>>
+<<check relative error>>
+<<model globals>>
+<<tension model test suite>>
+<<main check program>>
+@
+
+<<tension model check includes>>=
+#include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
+#include <stdio.h> /* sprintf() */
+#include <assert.h> /* assert() */
+#include <math.h> /* exp() */
+#include <errno.h> /* errno, ERANGE */
+<<check includes>>
+#include "global.h"
+#include "list.h"
+#include "tension_balance.h" /* oneD_*() */
+#include "tension_model.h"
+@
+
+<<tension model test suite>>=
+<<safe strto*>>
+<<const tension tests>>
+<<hooke tests>>
+<<worm-like chain tests>>
+<<freely-jointed chain tests>>
+<<piston tests>>
+<<tension model suite definition>>
+@
+
+<<tension model suite definition>>=
+Suite *test_suite (void)
+{
+ Suite *s = suite_create ("tension model");
+ <<const tension test case defs>>
+ <<hooke test case defs>>
+ <<worm-like chain test case defs>>
+ <<freely-jointed chain test case defs>>
+ <<piston test case defs>>
+
+ <<const tension test case adds>>
+ <<hooke test case adds>>
+ <<worm-like chain test case adds>>
+ <<freely-jointed chain test case adds>>
+ <<piston test case adds>>
+ return s;
+}
+@
+
+\subsubsection{Constant}
+
+<<const tension test case defs>>=
+TCase *tc_const_tension = tcase_create("Constant tension");
+@
+
+<<const tension test case adds>>=
+tcase_add_test(tc_const_tension, test_const_tension_create_destroy);
+suite_add_tcase(s, tc_const_tension);
+@
+
+<<const tension tests>>=
+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}
+
+<<hooke test case defs>>=
+TCase *tc_hooke = tcase_create("Hooke");
+@
+
+<<hooke test case adds>>=
+tcase_add_test(tc_hooke, test_hooke_create_destroy);
+suite_add_tcase(s, tc_hooke);
+
+@
+
+<<hooke tests>>=
+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}
+
+<<worm-like chain test case defs>>=
+TCase *tc_wlc = tcase_create("WLC");
+@
+
+<<worm-like chain test case adds>>=
+tcase_add_test(tc_wlc, test_wlc_at_zero);
+tcase_add_test(tc_wlc, test_wlc_at_half);
+suite_add_tcase(s, tc_wlc);
+
+@
+
+<<worm-like chain tests>>=
+<<worm-like chain function>>
+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
-<<tension model utility includes>>=
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h> /* strlen() */
-#include <assert.h> /* assert() */
-#include "global.h"
-#include "parse.h"
-#include "list.h"
-#include "tension_model.h"
-@
+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
-<<tension model utility definitions>>=
-<<version definition>>
-#define VFLAG 1 /* verbose */
-<<string comparison definition and globals>>
-<<tension model getopt definitions>>
-<<initialize model definition>>
+@
-@
-<<tension model utility getopt functions>>=
-<<index tension model>>
-<<tension model utility help>>
-<<tension model utility get options>>
-@
+\subsubsection{Freely-jointed chain}
-<<tension model utility help>>=
-void help(char *prog_name,
- environment_t *env,
- int n_tension_models, tension_model_getopt_t *tension_models,
- int tension_model)
+<<freely-jointed chain test case defs>>=
+TCase *tc_fjc = tcase_create("FJC");
+@
+
+<<freely-jointed chain test case adds>>=
+tcase_add_test(tc_fjc, test_fjc_at_zero);
+tcase_add_test(tc_fjc, test_fjc_at_half);
+suite_add_tcase(s, tc_fjc);
+
+@
+
+<<freely-jointed chain tests>>=
+<<freely-jointed chain function>>
+START_TEST(test_fjc_at_zero)
{
- 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 (x)\tForce (N)\n");
- printf(" 123.456\t7.89\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<n_tension_models; i++) {
- printf("\t%s\t%s\n", tension_models[i].name, tension_models[i].description);
- for (j=0; j < tension_models[i].num_params ; j++)
- printf("\t\t\t%s\n", tension_models[i].param_descriptions[j]);
- printf("\t\tdefaults: %s\n", tension_models[i].params);
- }
+ 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
-<<tension model utility get options>>=
-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,
- unsigned int *flags)
+START_TEST(test_fjc_at_half)
{
- char *prog_name = NULL;
- char c, options[] = "T:C:m:a:Vh";
- int tension_model=0, num_strings;
- extern char *optarg;
- extern int optind, optopt, opterr;
+ 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
- assert (argc > 0);
+@
- /* setup defaults */
- prog_name = argv[0];
- env->T = 300.0; /* K */
- *flags = 0;
- *model = tension_models;
+\subsubsection{Piston}
- while ((c=getopt(argc, argv, options)) != -1) {
- switch(c) {
- case 'T': env->T = atof(optarg); break;
- case 'C': env->T = atof(optarg)+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 '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);
- exit(1);
- }
+<<piston test case defs>>=
+TCase *tc_piston = tcase_create("Piston");
+@
+
+<<piston test case adds>>=
+tcase_add_test(tc_piston, test_piston_create_destroy);
+suite_add_tcase(s, tc_piston);
+
+@
+
+<<piston tests>>=
+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);
}
- return;
}
-@
+END_TEST
+
+@ TODO: In order to test the piston tension handler itself, we'd
+have to construct a group.
\section{Unfolding rate models}
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
+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$.
<<k func definition>>=
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]].
<<null k declarations>>
<<const k declarations>>
<<bell k declarations>>
+<<kbell k declarations>>
<<kramers k declarations>>
<<kramers integ k declarations>>
#endif /* K_MODEL_H */
-@
+@
<<k model module makefile lines>>=
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 utils makefile lines>>=
$(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
notangle -Rk-model-utils.c $< > $@
$(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
- gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
+ $(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 k declarations>>=
-@
+@
<<null k functions>>=
-@
+@
<<null k globals>>=
-@
+@
<<null k model getopt>>=
{"null", "a permanently unfolded domain", NULL, 0, NULL, NULL, NULL, NULL}
-@
+@
\subsection{Constant rate model}
\label{sec.k_const}
<<const k functions>>=
<<const k function>>
<<const k structure create/destroy functions>>
-@
+@
<<const k declarations>>=
<<const k function declaration>>
<<const k structure create/destroy function declarations>>
<<const k global declarations>>
-@
+@
<<const k structure>>=
typedef struct const_k_param_struct {
double knot; /* transition rate k_0 (frac population per s) */
} const_k_param_t;
-@
+@
<<const k function declaration>>=
double const_k(double F, environment_t *env, void *const_k_params);
-@
+@
<<const k function>>=
double const_k(double F, environment_t *env, void *const_k_params)
{ /* Returns the rate constant k in frac pop per s. */
assert(p->knot > 0);
return p->knot;
}
-@
+@
<<const k structure create/destroy function declarations>>=
void *string_create_const_k_param_t(char **param_strings);
void destroy_const_k_param_t(void *p);
-@
+@
<<const k structure create/destroy functions>>=
const_k_param_t *create_const_k_param_t(double knot)
void *string_create_const_k_param_t(char **param_strings)
{
- return create_const_k_param_t(atof(param_strings[0]));
+ return create_const_k_param_t(safe_strtod(param_strings[0], __FUNCTION__));
}
void destroy_const_k_param_t(void *p /* const_k_param_t * */)
<<const k model getopt>>=
{"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
+Quantitatively, Bell's model gives $k$ as
$$
k = k_0 \cdot e^{F dx / k_B T} \;,
$$
<<bell k functions>>=
<<bell k function>>
<<bell k structure create/destroy functions>>
-@
+@
<<bell k declarations>>=
<<bell k function declaration>>
<<bell k structure create/destroy function declarations>>
<<bell k global declarations>>
-@
+@
<<bell k structure>>=
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;
-@
+@
<<bell k function declaration>>=
double bell_k(double F, environment_t *env, void *bell_params);
-@
+@
<<bell k function>>=
double bell_k(double F, environment_t *env, void *bell_params)
{ /* Returns the rate constant k in frac pop per s.
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);
+ 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),
*/
return p->knot * exp(F*p->dx / (k_B*env->T) );
}
-@
+@
<<bell k structure create/destroy function declarations>>=
void *string_create_bell_param_t(char **param_strings);
void destroy_bell_param_t(void *p);
-@
+@
<<bell k structure create/destroy functions>>=
bell_param_t *create_bell_param_t(double knot, double dx)
void *string_create_bell_param_t(char **param_strings)
{
- return create_bell_param_t(atof(param_strings[0]), atof(param_strings[1]));
+ 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 * */)
<<bell k model getopt>>=
{"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-vazquez99b}.
+@
+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}.
+
+<<kbell k functions>>=
+<<kbell k function>>
+<<kbell k structure create/destroy functions>>
+@
+
+<<kbell k declarations>>=
+<<kbell k function declaration>>
+<<kbell k structure create/destroy function declarations>>
+<<kbell k global declarations>>
+
+@
+
+<<kbell k structure>>=
+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;
+@
+
+<<kbell k function declaration>>=
+double kbell_k(double F, environment_t *env, void *kbell_params);
+@
+<<kbell k function>>=
+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) );
+}
+@
+
+<<kbell k structure create/destroy function declarations>>=
+void *string_create_kbell_param_t(char **param_strings);
+void destroy_kbell_param_t(void *p);
+@
+
+<<kbell k structure create/destroy functions>>=
+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);
+}
+@
+
+<<kbell k global declarations>>=
+extern int num_kbell_params;
+extern char *kbell_param_descriptions[];
+extern char kbell_param_string[];
+@
+
+<<kbell k globals>>=
+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 k model getopt>>=
+{"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}
<<kramers k functions>>=
<<kramers k function>>
<<kramers k structure create/destroy functions>>
-@
+@
<<kramers k declarations>>=
<<kramers k function declaration>>
<<kramers k structure create/destroy function declarations>>
<<kramers k global declarations>>
-@
+@
<<kramers k structure>>=
//typedef double kramers_x_func_t(double F, double T, double *E_params, int n);
//double *E_params; /* parametrize the energy landscape */
//int n_E_params; /* length of E_params array */
} kramers_param_t;
-@
+@
<<kramers k function declaration>>=
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);
-@
+@
<<kramers k function>>=
double kramers_x_fn(FILE *in, FILE *out, char *fn, double F, double T)
{
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 atof(output);
+ return safe_strtod(output, __FUNCTION__);
}
double kramers_E_fn(FILE *in, FILE *out, char *fn, double F, double T, double x)
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 atof(output);
+ return safe_strtod(output, __FUNCTION__);
}
double kramers_E(double F, environment_t *env, void *kramers_params, double x)
//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);
}
-@
+@
<<kramers k structure create/destroy function declarations>>=
void *string_create_kramers_param_t(char **param_strings);
void destroy_kramers_param_t(void *p);
-@
+@
<<kramers k structure create/destroy functions>>=
kramers_param_t *create_kramers_param_t(double D,
/* takes an array of 4 functions: {"(1,2),(2,0),..."} */
void *string_create_kramers_param_t(char **param_strings)
{
- return create_kramers_param_t(atof(param_strings[0]),
+ return create_kramers_param_t(safe_strtod(param_strings[0], __FUNCTION__),
param_strings[2],
param_strings[3],
param_strings[4],
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
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{app.adaptive_dt} reduces it's timestep appropriately.
+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 k model getopt>>=
{"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}
<<kramers integ B k functions>>
<<kramers integ k function>>
<<kramers integ k structure create/destroy functions>>
-@
+@
<<kramers integ k declarations>>=
<<kramers integ k function declaration>>
<<kramers integ k structure create/destroy function declarations>>
<<kramers integ k global declarations>>
-@
+@
<<kramers integ k structure>>=
typedef double func_t(double x, void *params);
double F; /* for passing into GSL evaluations */
environment_t *env;
} kramers_integ_param_t;
-@
+@
<<spline param structure>>=
typedef struct spline_param_struct {
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.
$$
//fprintf(stderr, "integralA = %g\n", result);
return result;
}
-@
+@
$$
\text{Integral}_B(x_\text{min}, x_\text{max}) \equiv
//fprintf(stderr, "integralB = %g\n", result);
return result;
}
-@
+@
With the integrals taken care of, evaluating $k$ is simply
$$
$$
<<kramers integ k function declaration>>=
double kramers_integ_k(double F, environment_t *env, void *kramers_params);
-@
+@
<<kramers integ k function>>=
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);
}
-@
+@
<<kramers integ k structure create/destroy function declarations>>=
void *string_create_kramers_integ_param_t(char **param_strings);
void destroy_kramers_integ_param_t(void *p);
-@
+@
<<kramers integ k structure create/destroy functions>>=
kramers_integ_param_t *create_kramers_integ_param_t(double D,
//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(atof(param_strings[0]),
- atof(param_strings[2]),
- atof(param_strings[3]),
- &spline_eval, E_params,
- destroy_spline_param_t);
+ 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)
<<spline functions>>=
<<spline function>>
<<create/destroy spline>>
-@
+@
<<spline function>>=
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);
}
-@
+@
<<create/destroy spline>>=
spline_param_t *create_spline_param_t(int n_knots, double *x, double *y)
free(p);
}
}
-@
+@
<<kramers integ k global declarations>>=
extern int num_kramers_integ_params;
<<kramers integ k model getopt>>=
{"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-vazquez99b}.
+@
+Initialized with Titin I27 parameters\citep{carrion-vazquez99a}.
% Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
\subsection{Unfolding model implementation}
<<k model internal definitions>>
<<k model globals>>
<<k model functions>>
-@
+@
<<k model includes>>=
-#include <assert.h> /* assert() */
-#include <stdlib.h> /* NULL, malloc() */
-#include <math.h> /* HUGE_VAL, sqrt(), exp() */
-#include <stdio.h> /* fprintf(), stdout */
-#include <string.h> /* strlen(), strcpy() */
+#include <assert.h> /* assert() */
+#include <stdlib.h> /* NULL, malloc(), strto*() */
+#include <math.h> /* HUGE_VAL, sqrt(), exp() */
+#include <stdio.h> /* fprintf(), stdout */
+#include <string.h> /* strlen(), strcpy() */
+#include <errno.h> /* errno, ERANGE */
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include "global.h"
#include "parse.h"
-@
+@
<<k model internal definitions>>=
<<const k structure>>
<<bell k structure>>
+<<kbell k structure>>
<<kramers k structure>>
<<kramers integ k structure>>
<<spline param structure>>
-@
+@
<<k model globals>>=
<<null k globals>>
<<const k globals>>
<<bell k globals>>
+<<kbell k globals>>
<<kramers k globals>>
<<kramers integ k globals>>
-@
+@
<<k model functions>>=
+<<safe strto*>>
<<null k functions>>
<<const k functions>>
<<bell k functions>>
+<<kbell k functions>>
<<kramers k functions>>
<<kramers integ k functions>>
-@
+@
\subsection{Unfolding model unit tests}
Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
<<check-k-model.c>>=
+<<license comment>>
<<k model check includes>>
<<check relative error>>
<<model globals>>
<<k model test suite>>
<<main check program>>
-@
+@
<<k model check includes>>=
-#include <stdlib.h> /* EXIT_SUCCESS and EXIT_FAILURE, atof() */
-#include <stdio.h> /* sprintf() */
-#include <assert.h> /* assert() */
-#include <math.h> /* exp() */
+#include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
+#include <stdio.h> /* sprintf() */
+#include <assert.h> /* assert() */
+#include <math.h> /* exp() */
+#include <errno.h> /* errno, ERANGE */
<<check includes>>
#include "global.h"
#include "k_model.h"
-@
+@
<<k model test suite>>=
+<<safe strto*>>
<<const k tests>>
<<bell k tests>>
<<k model suite definition>>
-@
+@
<<k model suite definition>>=
Suite *test_suite (void)
<<bell k test case adds>>
return s;
}
-@
+@
\subsubsection{Constant}
<<const k test case defs>>=
TCase *tc_const_k = tcase_create("Constant k");
-@
+@
<<const k test case adds>>=
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);
-@
+@
<<const k tests>>=
int i;
for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
params[0] = knot[i];
- p = string_create_const_k_param_t(params);
+ p = string_create_const_k_param_t(params);
destroy_const_k_param_t(p);
}
}
env.T = T;
for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
params[0] = knot[i];
- p = string_create_const_k_param_t(params);
+ p = string_create_const_k_param_t(params);
for ( F=Fm; F<FM; F+=dF ) {
- fail_unless(const_k(F, &env, p)==atof(knot[i]),
+ fail_unless(const_k(F, &env, p)==safe_strtod(knot[i], __FUNCTION__),
"const_k(%g, %g, &{%s}) = %g != %s",
F, T, knot[i], const_k(F, &env, p), knot[i]);
}
}
}
END_TEST
-@
+@
\subsubsection{Bell}
<<bell k test case defs>>=
TCase *tc_bell = tcase_create("Bell's k");
-@
+@
<<bell k test case adds>>=
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);
-@
+@
<<bell k tests>>=
START_TEST(test_bell_k_create_destroy)
int i;
for( i=0; i < sizeof(knot)/sizeof(char *); i++) {
params[0] = knot[i];
- p = string_create_bell_param_t(params);
+ p = string_create_bell_param_t(params);
destroy_bell_param_t(p);
}
}
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)==atof(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);
char *knot[] = {"1","2","3","4","5","6"};
char *dx="1";
char *params[] = {knot[0], dx};
- double F= k_B*T/atof(dx);
+ double F= k_B*T/safe_strtod(dx, __FUNCTION__);
void *p = NULL;
environment_t env;
char param_string[80];
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, atof(knot[i])*exp(1.0), bell_k(F, &env, p));
- //fail_unless(bell_k(F, &env, p)==atof(knot[i])*exp(1.0),
+ 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), atof(knot[i])*exp(1.0));
+ // 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
-@
+@
<<kramers k tests>>=
-@
+@
-<<kramers k test case def>>=
-@
+<<kramers k test case defs>>=
+@
-<<kramers k test case add>>=
-@
+<<kramers k test case adds>>=
+@
<<k model function tests>>=
<<check relative error>>
fail_unless(F = k*x, NULL);
}
END_TEST
-@
+@
\subsection{Utilities}
<<k model utility getopt functions>>
<<k model utility multi model E>>
<<k model utility main>>
-@
+@
<<k model utility main>>=
int main(int argc, char **argv)
(*model->destructor)(params);
return 0;
}
-@
+@
<<k model utility multi model E>>=
double multi_model_E(k_func_t *k, void *model_params, environment_t *env, double x)
double E;
if (k == kramers_integ_k)
E = (*p->E)(x, p->E_params);
- else
+ else
E = kramers_E(0, env, model_params, x);
return E;
}
-@
+@
<<k model utility includes>>=
#include <stdlib.h>
#include <stdio.h>
-#include <string.h> /* strlen() */
-#include <assert.h> /* assert() */
+#include <string.h> /* strlen() */
+#include <assert.h> /* assert() */
+#include <errno.h> /* errno, ERANGE */
#include "global.h"
#include "parse.h"
#include "string_eval.h"
#include "k_model.h"
-@
+@
<<k model utility definitions>>=
<<version definition>>
<<kramers k structure>>
<<kramers integ k structure>>
-@
+@
<<k model utility getopt functions>>=
+<<safe strto*>>
<<index k model>>
<<k model utility help>>
<<k model utility get options>>
-@
+@
<<k model utility help>>=
void help(char *prog_name,
environment_t *env,
int n_k_models, k_model_getopt_t *k_models,
- int k_model)
+ int k_model, double Fmax, double special_xmin, double special_xmax)
{
int i, j;
printf("usage: %s [options]\n", prog_name);
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(" #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(" #Position (m)\tE (J)\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\n");
- printf("-x\tSet the minimum x value for the special mode E(x) output\n");
- printf("-X\tSet the maximum x value for the special mode E(x) output\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("\t\tdefaults: %s\n", k_models[i].params);
}
}
-@
+@
<<k model utility get options>>=
void get_options(int argc, char **argv, environment_t *env,
extern char *optarg;
extern int optind, optopt, opterr;
- assert (argc > 0);
+ assert(argc > 0);
/* setup defaults */
prog_name = argv[0];
- env->T = 300.0; /* K */
+ env->T = 300.0; /* K */
*mode = M_K_OF_F;
*flags = 0;
*model = k_models;
- *Fmax = 1e-9;
+ *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 = atof(optarg); break;
- case 'C': env->T = atof(optarg)+273.15; break;
+ 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;
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 = atof(optarg); break;
- case 'x': *special_xmin = atof(optarg); break;
- case 'X': *special_xmax = atof(optarg); break;
- case 'V': *flags |= VFLAG; 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);
+ help(prog_name, env, n_k_models, k_models, k_model, *Fmax, *special_xmin, *special_xmax);
exit(1);
}
}
return;
}
-@
+@
\section{\LaTeX\ documentation}
noweave -latex -index -delay $< > $@
$(BUILD_DIR)/sawsim.bib : $(SRC_DIR)/sawsim.bib | $(BUILD_DIR)
cp -f $< $@
-@
+@
and compiled using
<<latex makefile lines>>=
$(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
$(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 $@
@
[[bibtex]] and [[pdflatex]] produce quite a few intermediate files, so provide a means to clean them up.
<<latex makefile lines>>=
-semi-clean_latex :
+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.tex \
- $(BUILD_DIR)/sawsim.bib
+ $(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}
The sed is needed because notangle replaces tabs with eight spaces,
but [[make]] needs tabs.
<<makefile>>=
+# 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), ...)
+# 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_BINS =
+NW_SAWSIM_MODS =
+# Check binaries (check_*) defined in the noweb file
+CHECK_BINS =
# Modules defined outside the noweb file
FREE_SAWSIM_MODS = interp tavl
$(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 "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
+ $(SHELL) -e -c "rmdir $(BUILD_DIR) $(DOC_DIR)"
# Various builds of sawsim
$(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
- gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
+ $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
$(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
- gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
+ $(CC) $(CFLAGS) $(LDFLAGS) -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
$(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
- gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
+ $(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), ...)
+# 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)
$(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)
- gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
+ $(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
$$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
$$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
$$(BUILD_DIR)/global.h | $(BIN_DIR)
- gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
+ $(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 to
+# TODO: clean up the required modules too
$(CHECK_BINS:%=clean_%) :
rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)
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
\section{Math}
\subsection{Hookean springs in series}
-\label{app.math_hooke}
+\label{sec.math_hooke}
\begin{theorem}
The effective spring constant for $n$ springs $k_i$ in series is given by
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}}
+ 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}
\bibliography{sawsim}
\end{document}
-
projects / sawsim.git / blobdiff