From 1b3b40b2ff15fcc9d9486afd218240d3c6706af5 Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Fri, 21 Aug 2009 08:57:06 -0400 Subject: [PATCH] -P option now sets target P_N, not P_1. Bumped to version 0.9. This cleans up the trouble I had before with testing/bell_rate/bell_rate.sh 50 1e-6 0.05 1e-3 1e-9 300 where the large number of domains was making P_N much larger than P_1. This lead to more "ignored double unfoldings", creating a banded structure in the simulated histogram. The new '-P sets target P_N' configuration makes it more straightforward to use -P for it's intended purpose, avoiding "ignored double unfoldings", without forcing unnecessarily long timesteps once the population of the limiting state is reduced. Note that the second purpose of the -P option, avoiding significant force changes in a single timestep, is either unaffected (for states with population 1) or improved (for states with population > 1), but this should perhaps be separately controllable through it's own "--dF-max" option, or similar... --- src/sawsim.nw | 66 ++++++++++++++++++++++-------- testing/bell_rate/Makefile | 2 +- testing/bell_rate/NOTES | 75 ++++++++++++++++++++++++++++++++++ testing/bell_rate/bell_rate.sh | 2 +- testing/const_rate/NOTES | 67 ++++++++++++++++++++++++++++++ 5 files changed, 192 insertions(+), 20 deletions(-) create mode 100644 testing/bell_rate/NOTES create mode 100644 testing/const_rate/NOTES diff --git a/src/sawsim.nw b/src/sawsim.nw index 850e7a9..9917130 100644 --- a/src/sawsim.nw +++ b/src/sawsim.nw @@ -158,8 +158,12 @@ 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 the -P option to sawsim now sets the target probability +for the per-state transition $P_N$, not the per-domain transisiton +$P_1$. + <>= -#define VERSION "0.8" +#define VERSION "0.9" @ \subsection{License} @@ -546,38 +550,57 @@ $$ 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 +<>= +/* find multi-domain transition rate from the single domain rate */ +#define PN(P,N) (1.0-pow(1.0-(P), (N))) + +@ + +We can also define a macro to find the approprate $dt$ to achieve a +target $P_N$ with a given $k$ and $N$. +\begin{align} + P_N &= 1-(1-P_1)^N \\ + 1-P_1 &= (1-P_N)^{1/N} \\ + P_1 &= 1-(1-P_N)^{1/N} +\end{align} +<>= +#define P1(PN,N) (1.0-pow(1.0-(PN), 1.0/(N))) +@ + We take some time to discuss the meaning of $p(n>1)$ (i.e. multi-unfolding timesteps) in Section \ref{sec.timescales}. <>= -int domain_transitions(double F, double dt, environment_t *env, 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}. @@ -587,7 +610,7 @@ int random_transitions(list_t *transition_list, double tension, { /* 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--; @@ -1299,7 +1322,7 @@ int happens(double probability) $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 @@ -1315,14 +1338,14 @@ double search_dt(transition_t *transition, * 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; @@ -1355,7 +1378,8 @@ double search_dt(transition_t *transition, //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. */ @@ -1365,7 +1389,7 @@ double search_dt(transition_t *transition, 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; @@ -1549,6 +1573,12 @@ typedef struct transition_struct { /* get the transition data for the current list node */ #define T(list) ((transition_t *)(list)->d) + +/* get the number of initial state population for a transition state */ +#define N_INIT(transition) ((transition)->initial_state->num_domains) + +<> +<> @ [[k]] is a pointer to the function determining the transition rate for a given tension. [[k_params]] and [[destroy_k]] have similar roles with regards to [[k]] as [[tension_params]] and diff --git a/testing/bell_rate/Makefile b/testing/bell_rate/Makefile index ac6fd1b..557ae35 100644 --- a/testing/bell_rate/Makefile +++ b/testing/bell_rate/Makefile @@ -4,7 +4,7 @@ view : all qiv bell_rate.png clean : - rm -f bell_rate_*.gp bell_rate_*.hist bell_rate_*.d bell_rate_*.png bell_rate_*.compare + rm -f bell_rate_*.gp bell_rate_*.hist bell_rate_*.d bell_rate_*.png bell_rate_*.compare STDIN* bell_rate.png : run_test.sh ./$< diff --git a/testing/bell_rate/NOTES b/testing/bell_rate/NOTES new file mode 100644 index 0000000..87375aa --- /dev/null +++ b/testing/bell_rate/NOTES @@ -0,0 +1,75 @@ +* 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) diff --git a/testing/bell_rate/bell_rate.sh b/testing/bell_rate/bell_rate.sh index 2f4e5e6..6aee491 100755 --- a/testing/bell_rate/bell_rate.sh +++ b/testing/bell_rate/bell_rate.sh @@ -90,7 +90,7 @@ then # Sawsim >= 0.6 qcmd -w14:00:00 xtimes $N $SAWSIM -T$T -v$V -s cantilever,hooke,$SPRING -N1 \ -s folded,null -N$NDOM -s unfolded,null \ - -k "'folded,unfolded,bell,{$K0,$DX}'" -q folded -P 1e-5 \ + -k "'folded,unfolded,bell,{$K0,$DX}'" -q folded \ | grep -v '^#' | cut -f1 >> $DATA else # Sawsim >= 0.6 diff --git a/testing/const_rate/NOTES b/testing/const_rate/NOTES new file mode 100644 index 0000000..7a066b5 --- /dev/null +++ b/testing/const_rate/NOTES @@ -0,0 +1,67 @@ +* 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) -- 2.26.2