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 the -P option to sawsim now sets the target probability
+for the per-state transition $P_N$, not the per-domain transisiton
+$P_1$.
+
<<version definition>>=
-#define VERSION "0.8"
+#define VERSION "0.9"
@
\subsection{License}
So the probability of a given domain transitioning in a short time
$dt$ is given by
$$
- P = k \cdot dt.
+ P_1 = k \cdot dt.
$$
-For a group of $N$ identical domains, and therefore identicalunfolding
-probabilities $P$, the probability of $n$ domains unfolding in a given
-timestep is given by the binomial distribution
+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^n(1-P)^(N-n).
+ 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)^N,
+ P(0) = (1-P_1)^N,
$$
so the probability that \emph{at least one} domain unfolds is
$$
- p(n>0) = 1-(1-P)^N.
+ P_N = 1-(1-P_1)^N.
$$
Note that for small $P$,
$$
- p(n>0) = 1-(1-NP+\mathcal{O}(P^2)) = NP - \mathcal{O}(P^2) \approx NP.
+ 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)))
+
+@
+
+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, int num_domains, transition_t *transition)
-{ /* returns 1 or 0, F in N, dt in s, pointer to env. data, number of 'reactant' domains, pointer to transition structure */
+int domain_transitions(double F, double dt, environment_t *env, transition_t *transition)
+{ /* returns 1 or 0, F in N, dt in s, pointer to env. data, pointer to transition structure */
double k;
k = accel_k(transition->k, F, env, transition->k_params);
//(*transition->k)(F, env, domain->k_params);
//printf("k = %g,\tdt = %g,\tk dt = %g\n", k, dt, k*dt);
- return happens(1-pow((1.0-k*dt), num_domains)); /* N dice rolls for prob. k*dt event */
+ return happens(PN(k*dt, N_INIT(transition)));
}
@ [[happens]] is a random decision making function defined in Appendix \ref{sec.utils}.
{ /* 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)->initial_state->num_domains, T(transition_list))) {
+ && 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--;
$F(x)$ increases with $x$, possibly exploding (e.g.\ 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
+for which $P_N(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
* a starting separation x in m, a target_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(state_list, env, x, 0); /* BUG. repeated calculation */
//printf("Start with x = %g (F = %g)\n", x, F);
k = accel_k(transition->k, F, env, transition->k_params);
//printf("x %g\tF %g\tk %g\n", x, F, k);
- dtU = target_prob / k; /* P = k dt, dtU is an upper bound on dt */
+ dtU = P1(target_prob, N_INIT(transition)) / k; /* upper bound on dt */
if (dtU > max_dt) {
//printf("overshot max_dt\n");
dtU = max_dt;
//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(target_prob, N_INIT(transition)) / k;
if (dtUCur >= dt)
return dt; /* dtU is safe. */
F = find_tension(state_list, env, x+v*dt, 0);
//printf("Try for dt = %g (F = %g) (dt bounds %g, %g)\n", dt, F, dtL, dtU);
k = accel_k(transition->k, F, env, transition->k_params);
- dtCur = target_prob / k;
+ dtCur = P1(target_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;
/* 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
--- /dev/null
+* Version 0.8 MOSTLY PASSED
+
+Problems with 50_1e-6_0.05_1e-3_1e-9_300's
+ z = K0/kv: 17106.1 (expected 20000.0)
+Strange banded structure too. Banding most pronounced for smaller
+forces. The banding is due to the double-unfolding problem. Reducing
+P from 1e-3 to 1e-5 (which takes a lot longer to run), gave
+50_1e-6_0.05_1e-3_1e-9_299, which looks much nicer.
+
+Otherwise things look ok.
+
+==> bell_rate_1_1_1_1_1_7.24296e22.compare <==
+z = K0/kv: 1.00068 (expected 1.0)
+b = dx/kBT: 0.999718 (expected 1.00000051031)
+
+==> bell_rate_1_1_1_1_5_7.24296e22.compare <==
+z = K0/kv: 0.985894 (expected 1.0)
+b = dx/kBT: 4.99089 (expected 5.00000255156)
+
+==> bell_rate_1_1_1_5_1_7.24296e22.compare <==
+z = K0/kv: 4.98609 (expected 5.0)
+b = dx/kBT: 0.981462 (expected 1.00000051031)
+
+==> bell_rate_1_1_5_1_1_7.24296e22.compare <==
+z = K0/kv: 0.197084 (expected 0.2)
+b = dx/kBT: 0.999987 (expected 1.00000051031)
+
+==> bell_rate_1_1e-6_0.05_1e-3_1e-9_300.compare <==
+z = K0/kv: 20512.9 (expected 20000.0)
+b = dx/kBT: 2.41336e+11 (expected 241432123206.0)
+
+==> bell_rate_1_5_1_1_1_7.24296e22.compare <==
+z = K0/kv: 0.204091 (expected 0.2)
+b = dx/kBT: 0.982143 (expected 1.00000051031)
+
+==> bell_rate_50_1e-6_0.05_1e-3_1e-9_300.compare <==
+z = K0/kv: 17106.1 (expected 20000.0)
+b = dx/kBT: 2.4366e+11 (expected 241432123206.0)
+
+==> bell_rate_5_1_1_1_1_7.24296e22.compare <==
+z = K0/kv: 1.00114 (expected 1.0)
+b = dx/kBT: 1.00475 (expected 1.00000051031)
+
+* Version 0.9 PASSED
+==> bell_rate_1_1_1_1_1_7.24296e22.compare <==
+z = K0/kv: 1.01076 (expected 1.0)
+b = dx/kBT: 0.984221 (expected 1.00000051031)
+
+==> bell_rate_1_1_1_1_5_7.24296e22.compare <==
+z = K0/kv: 1.0312 (expected 1.0)
+b = dx/kBT: 4.93899 (expected 5.00000255156)
+
+==> bell_rate_1_1_1_5_1_7.24296e22.compare <==
+z = K0/kv: 5.18165 (expected 5.0)
+b = dx/kBT: 0.776549 (expected 1.00000051031)
+
+==> bell_rate_1_1_5_1_1_7.24296e22.compare <==
+z = K0/kv: 0.198328 (expected 0.2)
+b = dx/kBT: 1.00328 (expected 1.00000051031)
+
+==> bell_rate_1_1e-6_0.05_1e-3_1e-9_300.compare <==
+z = K0/kv: 17096.4 (expected 20000.0)
+b = dx/kBT: 2.4364e+11 (expected 241432123206.0)
+
+==> bell_rate_1_5_1_1_1_7.24296e22.compare <==
+z = K0/kv: 0.201196 (expected 0.2)
+b = dx/kBT: 0.984595 (expected 1.00000051031)
+
+==> bell_rate_50_1e-6_0.05_1e-3_1e-9_300.compare <==
+z = K0/kv: 18893 (expected 20000.0)
+b = dx/kBT: 2.42285e+11 (expected 241432123206.0)
+
+==> bell_rate_5_1_1_1_1_7.24296e22.compare <==
+z = K0/kv: 1.01665 (expected 1.0)
+b = dx/kBT: 0.979307 (expected 1.00000051031)
--- /dev/null
+* Version 0.8 PASSED
+
+==> const_rate_1_1_1_1.compare <==
+force scale: -1.01187 (expected -1.0)
+initial value: 991.59 (expected 1000.0)
+
+==> const_rate_1_1_1_2.compare <==
+force scale: -0.504875 (expected -0.5)
+initial value: 995.942 (expected 1000.0)
+
+==> const_rate_1_1_2_1.compare <==
+force scale: -1.95167 (expected -2.0)
+initial value: 510.585 (expected 500.0)
+
+==> const_rate_1_2_1_1.compare <==
+force scale: -2.02266 (expected -2.0)
+initial value: 498.409 (expected 500.0)
+
+==> const_rate_50_1_1_1.compare <==
+force scale: -1.01996 (expected -1.0)
+initial value: 49270 (expected 50000.0)
+
+==> const_rate_50_1_1_2.compare <==
+force scale: -0.507475 (expected -0.5)
+initial value: 49484.4 (expected 50000.0)
+
+==> const_rate_50_1_2_1.compare <==
+force scale: -2.02617 (expected -2.0)
+initial value: 24796.2 (expected 25000.0)
+
+==> const_rate_50_2_1_1.compare <==
+force scale: -2.03676 (expected -2.0)
+initial value: 24674.6 (expected 25000.0)
+
+* Version 0.9 PASSED
+
+==> const_rate_1_1_1_1.compare <==
+force scale: -1.0122 (expected -1.0)
+initial value: 996.249 (expected 1000.0)
+
+==> const_rate_1_1_1_2.compare <==
+force scale: -0.498754 (expected -0.5)
+initial value: 1006.67 (expected 1000.0)
+
+==> const_rate_1_1_2_1.compare <==
+force scale: -1.97399 (expected -2.0)
+initial value: 504.787 (expected 500.0)
+
+==> const_rate_1_2_1_1.compare <==
+force scale: -2.01052 (expected -2.0)
+initial value: 497.477 (expected 500.0)
+
+==> const_rate_50_1_1_1.compare <==
+force scale: -1.00178 (expected -1.0)
+initial value: 49942.3 (expected 50000.0)
+
+==> const_rate_50_1_1_2.compare <==
+force scale: -0.498233 (expected -0.5)
+initial value: 50155.7 (expected 50000.0)
+
+==> const_rate_50_1_2_1.compare <==
+force scale: -2.00883 (expected -2.0)
+initial value: 24920.6 (expected 25000.0)
+
+==> const_rate_50_2_1_1.compare <==
+force scale: -1.99658 (expected -2.0)
+initial value: 25031.5 (expected 25000.0)
projects / sawsim.git / commitdiff