From: W. Trevor King <wking@drexel.edu>
Date: Wed, 19 Aug 2009 16:59:49 +0000 (-0400)
Subject: Extended bell and const rate tests to multi-domain chains
X-Git-Tag: v0.9~1
X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=ae0666a896172a7bad990f422a09ac15ecf84df2;p=sawsim.git

Extended bell and const rate tests to multi-domain chains
---

diff --git a/.be/bugs/be888cb1-605e-48dd-ae76-6d418b8373f9/comments/3c6ed1f3-0702-4ce7-a3f0-e171d2403b52/body b/.be/bugs/be888cb1-605e-48dd-ae76-6d418b8373f9/comments/3c6ed1f3-0702-4ce7-a3f0-e171d2403b52/body
new file mode 100644
index 0000000..013fa5d
--- /dev/null
+++ b/.be/bugs/be888cb1-605e-48dd-ae76-6d418b8373f9/comments/3c6ed1f3-0702-4ce7-a3f0-e171d2403b52/body
@@ -0,0 +1 @@
+To answer "How much does the X stretch?"
diff --git a/.be/bugs/be888cb1-605e-48dd-ae76-6d418b8373f9/comments/3c6ed1f3-0702-4ce7-a3f0-e171d2403b52/values b/.be/bugs/be888cb1-605e-48dd-ae76-6d418b8373f9/comments/3c6ed1f3-0702-4ce7-a3f0-e171d2403b52/values
new file mode 100644
index 0000000..dcd9419
--- /dev/null
+++ b/.be/bugs/be888cb1-605e-48dd-ae76-6d418b8373f9/comments/3c6ed1f3-0702-4ce7-a3f0-e171d2403b52/values
@@ -0,0 +1,8 @@
+Author: W. Trevor King <wking@drexel.edu>
+
+
+Content-type: text/plain
+
+
+Date: Sat, 15 Aug 2009 11:36:32 +0000
+
diff --git a/.be/bugs/be888cb1-605e-48dd-ae76-6d418b8373f9/values b/.be/bugs/be888cb1-605e-48dd-ae76-6d418b8373f9/values
new file mode 100644
index 0000000..51a41c9
--- /dev/null
+++ b/.be/bugs/be888cb1-605e-48dd-ae76-6d418b8373f9/values
@@ -0,0 +1,17 @@
+creator: W. Trevor King <wking@drexel.edu>
+
+
+reporter: W. Trevor King <wking@drexel.edu>
+
+
+severity: minor
+
+
+status: open
+
+
+summary: Allow selectable tension-group length output in -V runs.
+
+
+time: Sat, 15 Aug 2009 11:36:03 +0000
+
diff --git a/testing/bell_rate/Makefile b/testing/bell_rate/Makefile
new file mode 100644
index 0000000..ac6fd1b
--- /dev/null
+++ b/testing/bell_rate/Makefile
@@ -0,0 +1,10 @@
+all : bell_rate.png
+
+view : all
+	qiv bell_rate.png
+
+clean :
+	rm -f bell_rate_*.gp bell_rate_*.hist bell_rate_*.d bell_rate_*.png bell_rate_*.compare
+
+bell_rate.png : run_test.sh
+	./$<
diff --git a/testing/bell_rate/bell_rate.sh b/testing/bell_rate/bell_rate.sh
new file mode 100755
index 0000000..2f4e5e6
--- /dev/null
+++ b/testing/bell_rate/bell_rate.sh
@@ -0,0 +1,136 @@
+#!/bin/bash
+#
+# Use a single hooke domain with a constant unfolding rate, so
+# unfolding force is proportional to time, and connect multiple
+# identical Bell unfolding domains.  We expect an peaked death
+# histogram.
+#
+# e.g.
+#   spring constant k=df/dx    \__
+#   velocity        v=dx/dt    /  `--> df/dt = kv --> f = kvt  (+ f(t=0))
+#   unfolding rate  K= a e^(bf)   (frac/s)
+#   let f(t=0) = 0
+#   let P(f) be the population existing at force f (normed so P(t=f=0)=1).
+#   the probability of dying (normed # deaths) between f and f+df is given by
+#        dD = K P(f) * dt = K P(f) df/kv     since dt = df/kv
+#   this is what we plot in the normalized histogram with dD in a bin of width
+#   df.  The continuous analog is the function dD/df = K/kv P(f).
+#
+#   To find P(f), we need to solve
+#        dP/dt = -KP
+#        dP/P = -K dt = -a e^{bf} df/kv = -z e^{bf} df       where z ≡ a/kv
+#        ln(P) = -z/b e^{bf) + c
+#        P = C e^{-z/b e^{bf}}
+#   Applying our f(t=0)=0 boundary condition we see that C = P(t=f=0) e^{z/b}.
+#   Plugging P(f) into our formula for dD/df yields
+#        dD/df = K/kv P(f) = z e^{bf}  P(t=f=0) e^{z/b} e^{-z/b e^{bf}}
+#        dD/df = P(t=f=0) z e^{bf + z/b(1-e^{bf})}
+#
+#   The histogram scaling isn't dD/df, though, it's dD/d(bin)
+#        dD/d(bin) = dD/df df/d(bin)
+#                  = binwidth P(t=f=0) z e^{bf + z/b(1-e^{bf})}
+#
+# Reminders:
+#   spring constant k=df/dx       (N/m)
+#   velocity        v=dx/dt       (m/s)
+#   unfolding rate  K= a e^(bf) = K0 e^(fdx/kBT)  (frac/s)
+# usage: bell_rate.sh <num_domains> <v> <k> <K0> <dx> <T>
+
+if [ $# -ne 6 ]
+then
+    echo "usage: $0 <num_domains> <v> <k> <K0> <dx> <T>"
+    echo ""
+    echo "where"
+    echo "spring constant k=df/dx             (N/m)"
+    echo "velocity        v=dx/dt             (m/s)"
+    echo "unfolding rate  K= K0 e^(fdx/kBT)   (frac/s)"
+    exit 1
+fi
+
+NDOM=$1
+V=$2
+SPRING=$3
+K0=$4
+DX=$5
+T=$6
+DEBUG=0
+
+# since we're on the cluster, we'll go crazy and get really good statistics ;)
+N=10000
+
+SAWSIM=../../bin/sawsim
+DATA=bell_rate_${NDOM}_${V}_${SPRING}_${K0}_${DX}_${T}.d
+HIST=bell_rate_${NDOM}_${V}_${SPRING}_${K0}_${DX}_${T}.hist
+COMPFILE=bell_rate_${NDOM}_${V}_${SPRING}_${K0}_${DX}_${T}.compare
+GPFILE=bell_rate_${NDOM}_${V}_${SPRING}_${K0}_${DX}_${T}.gp
+PNGFILE=bell_rate_${NDOM}_${V}_${SPRING}_${K0}_${DX}_${T}.png
+
+# set up the theory
+# our predicted histogram (see notes above) is
+#   dD/d(bin) = binwidth P0 z e^{bf + z/b(1-e^{bf})}
+# where z ≡ a/kv, and b = Δx/kBT
+KB=1.3806503e-23 # J/K
+B=`python -c "print ($DX/($KB*$T))"`
+Z=`python -c "print (float($K0)/($SPRING*$V))"`
+# find a good df (histogram bin width).
+#    currently, scale based only on e^(bf) term in exponent.
+df=`python -c "print (0.01/$B)"`
+
+THEORY_FN="$df * $N * $NDOM * $Z * exp($B*x + $Z/$B *(1 - exp($B*x)))"
+if [ $DEBUG -eq 1 ]; then
+    echo "b = dx/kBT = $B"
+    echo "z = K0/kv = $Z"
+    echo "Theory --> P(x=f) = $THEORY_FN"
+fi
+
+# run the experiments
+> $DATA
+if [ "$ISACLUSTER" -eq 1 ]
+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 \
+        | grep -v '^#' | cut -f1 >> $DATA
+else
+    # Sawsim >= 0.6
+    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 \
+        | grep -v '^#' | cut -f1 >> $DATA
+fi
+
+# histogram the data
+if [ $DEBUG -eq 1 ]; then
+    echo "cat $DATA | stem_leaf -v-1 -b$df -c | sed 's/://' > $HIST"
+fi
+cat $DATA | stem_leaf -v-1 -b$df -c | sed 's/://' > $HIST
+wc -l "$HIST"
+# fit the data to an exponential decay
+if [ $DEBUG -eq 1 ]; then
+    echo "cat $HIST | ./fit_bell"
+fi
+FIT=`cat $HIST | ./fit_bell` || exit 1
+FITZ=`echo "$FIT" | sed -n 's/a:[[:space:]]*//p'`
+FITB=`echo "$FIT" | sed -n 's/b:[[:space:]]*//p'`
+echo -e "z = K0/kv:\t$FITZ\t(expected $Z)" > $COMPFILE
+echo -e "b = dx/kBT:\t$FITB\t(expected $B)" >> $COMPFILE
+cat $COMPFILE
+
+# generate gnuplot script
+if [ $DEBUG -eq 1 ]; then
+    echo "generate Gnuplot script"
+fi
+echo "set terminal png" > $GPFILE
+echo "set output '$PNGFILE'" >> $GPFILE
+echo "set logscale y" >> $GPFILE
+echo "set xlabel 'Force'" >> $GPFILE
+echo "set ylabel 'Deaths'" >> $GPFILE
+echo "theory(x) = $THEORY_FN" >> $GPFILE
+echo "fit(x) = $df * $N *  $NDOM * $FITZ * exp($FITB*x + $FITZ/$FITB*(1-exp($FITB*x)))" >> $GPFILE
+echo "plot  theory(x), '$HIST' with points, fit(x)" >> $GPFILE
+
+if [ $DEBUG -eq 1 ]; then
+    echo "Generate PNG"
+fi
+gnuplot $GPFILE
diff --git a/testing/bell_rate/params b/testing/bell_rate/params
new file mode 100644
index 0000000..702a353
--- /dev/null
+++ b/testing/bell_rate/params
@@ -0,0 +1,11 @@
+# num domains	velocity	spring constant	unfolding rate	unfolding distance	temperature
+# 7.24296e22 = 1/kB  (in SI units)
+1	1	1	1	1	7.24296e22
+5	1	1	1	1	7.24296e22
+1	5	1	1	1	7.24296e22
+1	1	5	1	1	7.24296e22
+1	1	1	5	1	7.24296e22
+1	1	1	1	5	7.24296e22
+# Now use reasonable protein parameters
+1	1e-6	0.05	1e-3	1e-9	300
+50	1e-6	0.05	1e-3	1e-9	300
diff --git a/testing/bell_rate/run_test.sh b/testing/bell_rate/run_test.sh
index 78421ca..f77a230 100755
--- a/testing/bell_rate/run_test.sh
+++ b/testing/bell_rate/run_test.sh
@@ -1,114 +1,13 @@
 #!/bin/bash
 #
-# Use a single hooke domain with a bell unfolding rate, so unfolding force
-# is proportional to time, and we expect a peaked death histogram.
+# run const_rate.sh with a variety of parameters
 #
-# e.g.
-#   spring constant k=df/dx=1    \__
-#   velocity        v=dx/dt=1    /  `--> df/dt = kv = 1  --> f = t  (+ f0)
-#   unfolding rate  K= a e^(bf)    (frac/s)
-#   let f0 = 0
-#   let P(f) be the population existing at force f (normed so P(0)=0).
-#   the probability of dying (normed # deaths) between f and f+df is given by
-#        dD = K P(f) * dt = K P(f) df/kv     since dt = df/kv
-#   this is what whe plot in the normalized histogram with dD in a bin of width
-#   df.  The continuous analog is the function dD/df = K/kv P(f).
-#
-#        dP/dt = - KP
-#        dP/P = -K dt = -a e^{bf} df/kv = -z e^{bf} df       where z ≡ a/kv
-#        ln(P) = -z/b e^{bf) + c
-#        P = C e^{-z/b e^{bf}}
-#   Applying our f0=0 boundary condition we see that C = P0 e^{z/b}.
-#   Plugging P(f) into our formula for dD/df yields
-#        dD/df = K/kv P(f) = z e^{bf}  P0 e^{z/b} e^{-z/b e^{bf}}
-#        dD/df = P0 z e^{bf + z/b(1-e^{bf})}
-#
-#   The histogram scaling isn't dD/df, though, it's dD/d(bin)
-#        dD/d(bin) = dD/df df/d(bin) = binwidth P0 z e^{bf + z/b(1-e^{bf})}
-#
-#   test with lots of runs of
-#              v- T=1/kB    v- v=1            v- k=1
-#     sawsim -T7.24296e22 -v1 -s folded,hooke,1 -N1 -s unfolded,null \
-#         -k "folded,unfolded,bell,{1,1}" -q folded
-#                                   ^- {a,dx}={1,1}
-#   (kB = 1.3806503 10-23 J/K, we use T=1/kB to get the bell K=e^{-f})
-#
-#   60 runs of that takes ~1 second:
-#     time xtimes 60 sawsim -T7.24296e22 -v1 -s folded,hooke,1 -N1 \
-#         -s unfolded,null -k "folded,unfolded,bell,{1,1}" -q folded \
-#         > /dev/null
-#
-# usage: bell_rate.sh
-
-# since we're on the cluster, we'll go crazy and get really good statistics ;)
-N=10000
-
-SAWSIM=../../bin/sawsim
-DATA=bell_rate.d
-HIST=bell_rate.hist
-GPFILE=bell_rate.gp
-PNGFILE=bell_rate.png
-
-# set up the theory
-#   dD/d(bin) = binwidth P0 z e^{bf + z/b(1-e^{bf})}
-# where z ≡ a/kv, and b = Δx/kBT
-KB=1.3806503e-23 # J/K
-#STYLE="ones"
-STYLE="protein"
-if [ "$STYLE" == "ones" ]
-then
-    df=1e-1  # histogram bin width
-    DX=1
-    T=7.24296e22 # 1/kB
-    A=1
-    K=1
-    V=1
-else
-    df=1e-12
-    DX=1e-9
-    T=300
-    A=1e-3
-    K=0.05
-    V=1e-6
-fi
-B=`echo "print ($DX/($KB*$T))" | python`
-Z=`echo "print (float($A)/($K*$V))" | python`
-THEORY_FN="$df * $N * $Z * exp($B*x + $Z/$B *(1 - exp($B*x)))"
-
-# run the experiments
-> $DATA
-if [ "$ISACLUSTER" -eq 1 ]
-then
-    # Sawsim <= 0.5
-    #qcmd "xtimes $N $SAWSIM -T$T -v$V -mhooke -a$K -kbell -K$A,$DX -F1 | grep -v '^#' >> $DATA"
-    # Sawsim >= 0.6
-    qcmd "xtimes $N $SAWSIM -T$T -v$V -s folded,hooke,$K -N1 -s unfolded,null -k 'folded,unfolded,bell,{$A,$DX}' -q folded | grep -v '^#' | cut -f1 >> $DATA"
-else
-    # Sawsim <= 0.5
-    #xtimes $N $SAWSIM -T$T -v$V -mhooke -a$K -kbell -K$A,$DX -F1 | grep -v '^#' >> $DATA
-    # Sawsim >= 0.6
-    xtimes $N $SAWSIM -T$T -v$V -s folded,hooke,$K -N1 -s unfolded,null -k "folded,unfolded,bell,{$A,$DX}" -q folded | grep -v '^#' | cut -f1 >> $DATA
-fi
-
-# histogram the data
-cat $DATA | stem_leaf -v-1 -b$df | sed 's/://' > $HIST
-
-# fit the data
-FIT=`cat $HIST | ./fit_bell`
-#echo "$FIT"
-FITZ=`echo "$FIT" | sed -n 's/a:[[:space:]]*//p'`
-FITB=`echo "$FIT" | sed -n 's/b:[[:space:]]*//p'`
-echo -e "z:\t$FITZ\t(expected $Z)"
-echo -e "b:\t$FITB\t(expected $B)"
+# usage: run_test.sh
 
-# generate gnuplot script
-echo "set terminal png" > $GPFILE
-echo "set output '$PNGFILE'" >> $GPFILE
-echo "set logscale y" >> $GPFILE
-echo "set xlabel 'Time'" >> $GPFILE
-echo "set ylabel 'Deaths'" >> $GPFILE
-echo "theory(x) = $THEORY_FN" >> $GPFILE
-echo "fit(x) = $df * $N * $FITZ * exp($FITB*x + $FITZ/$FITB*(1-exp($FITB*x)))" >> $GPFILE
-echo "plot  theory(x), '$HIST' with points, fit(x)" >> $GPFILE
+while read LINE
+do
+  echo ./bell_rate.sh $LINE
+  ./bell_rate.sh $LINE
+done < <(cat params | grep -v '^#')
 
-gnuplot $GPFILE
+exit 0
diff --git a/testing/const_rate/Makefile b/testing/const_rate/Makefile
new file mode 100644
index 0000000..b638aa8
--- /dev/null
+++ b/testing/const_rate/Makefile
@@ -0,0 +1,10 @@
+all : const_rate.png
+
+view : all
+	qiv const_rate.png
+
+clean :
+	rm -f const_rate_*.gp const_rate_*.hist const_rate_*.d const_rate_*.png const_rate_*.compare
+
+const_rate.png : run_test.sh
+	./$<
diff --git a/testing/const_rate/const_rate.sh b/testing/const_rate/const_rate.sh
index 56ae14f..baad49c 100755
--- a/testing/const_rate/const_rate.sh
+++ b/testing/const_rate/const_rate.sh
@@ -1,7 +1,9 @@
 #!/bin/bash
 #
-# Use a single hooke domain with a constant unfolding rate, so unfolding force
-# is proportional to time, and we expect an exponential population decay.
+# Use a single hooke domain with a constant unfolding rate, so
+# unfolding force is proportional to time, and connect multiple
+# identical constant unfolding domains.  We expect an exponential
+# population decay.
 #
 # e.g.
 #   spring constant k=df/dx      \__
@@ -13,22 +15,18 @@
 #   so histogram of unfolding vs. force p(F) normalized to p(0)=1 should follow
 #     -binwidth N0 dp/dF = binwidth N0 pK/kv = binwidth N0 K/kv exp(-FK/kv)
 #
-#   test with lots of runs of
-#     sawsim -v1 -mhooke -a1 -kconst -K1 -F1
-#   which is not a problem, because 200 runs of that takes ~1 second:
-#     time xtimes 200 sawsim -v1 -mhooke -a1 -kconst -K1 -F1 > /dev/null
-#
-# usage: const_rate.sh <velocity> <spring_constant> <unfolding_rate>
+# usage: const_rate.sh <num_domains> <velocity> <spring_constant> <unfolding_rate>
 
-if [ $# -ne 3 ]
+if [ $# -ne 4 ]
 then
-    echo "usage: const_rate.sh <velocity> <spring_constant> <unfolding_rate>"
+    echo "usage: $0 <num_domains> <velocity> <spring_constant> <unfolding_rate>"
     exit 1
 fi
 
-V=$1
-K=$2
-RATE=$3
+NDOM=$1
+V=$2
+K=$3
+RATE=$4
 DEBUG=0
 
 # since we're on the cluster, we'll go crazy and get really good statistics ;)
@@ -36,16 +34,16 @@ N=10000
 df=`python -c "print 0.1/$RATE"`
 
 SAWSIM=../../bin/sawsim
-DATA=const_rate_${V}_${K}_$RATE.d
-HIST=const_rate_${V}_${K}_$RATE.hist
-COMPFILE=const_rate_${V}_${K}_$RATE.compare
-GPFILE=const_rate_${V}_${K}_$RATE.gp
-PNGFILE=const_rate_${V}_${K}_$RATE.png
+DATA=const_rate_${NDOM}_${V}_${K}_${RATE}.d
+HIST=const_rate_${NDOM}_${V}_${K}_${RATE}.hist
+COMPFILE=const_rate_${NDOM}_${V}_${K}_${RATE}.compare
+GPFILE=const_rate_${NDOM}_${V}_${K}_${RATE}.gp
+PNGFILE=const_rate_${NDOM}_${V}_${K}_${RATE}.png
 
 # set up the theory
 # our predicted histogram (see notes above) is
 #  hist(F) = binwidth N0 K/kv exp(-FK/kv) = A exp(F/B)
-theoryA=`python -c "print $df*$N*$RATE/float($K*$V)"`
+theoryA=`python -c "print $df*$N*$NDOM*$RATE/float($K*$V)"`
 theoryB=`python -c "print float(-$K*$V)/$RATE"`
 THEORY_FN="$theoryA * exp(x / ($theoryB))"
 
@@ -54,26 +52,24 @@ THEORY_FN="$theoryA * exp(x / ($theoryB))"
 > $DATA
 if [ "$ISACLUSTER" -eq 1 ]
 then
-    # Sawsim <= 0.5
-    #qcmd "xtimes $N $SAWSIM -v$V -mhooke -a$K -kconst -K$RATE -F1 | grep -v '^#' >> $DATA"
     # Sawsim >= 0.6
-    qcmd xtimes $N $SAWSIM -v$V -s folded,hooke,$K -N1 -s unfolded,null \
+    qcmd xtimes $N $SAWSIM -v$V -s cantilever,hooke,$K -N1 \
+	-s folded,null -N$NDOM -s unfolded,null \
         -k folded,unfolded,const,$RATE -q folded \
         | grep -v '^#' | cut -f1 >> $DATA
 else
-    # Sawsim <= 0.5
-    #xtimes $N $SAWSIM -v$V -mhooke -a$K -kconst -K$RATE -F1 | grep -v '^#' >> $DATA
     # Sawsim >= 0.6
-    xtimes $N $SAWSIM -v$V -s folded,hooke,$K -N1 -s unfolded,null \
+    xtimes $N $SAWSIM -v$V -s cantilever,hooke,$K -N1 \
+	-s folded,null -N$NDOM -s unfolded,null \
         -k folded,unfolded,const,$RATE -q folded \
         | grep -v '^#' | cut -f1 >> $DATA
 fi
 
 # histogram the data
 if [ $DEBUG -eq 1 ]; then
-    echo "cat $DATA | stem_leaf -v-1 -b$df | sed 's/://' > $HIST"
+    echo "cat $DATA | stem_leaf -v-1 -b$df -c | sed 's/://' > $HIST"
 fi
-cat $DATA | stem_leaf -v-1 -b$df | sed 's/://' > $HIST
+cat $DATA | stem_leaf -v-1 -b$df -c | sed 's/://' > $HIST
 
 # fit the data to an exponential decay
 if [ $DEBUG -eq 1 ]; then
diff --git a/testing/const_rate/params b/testing/const_rate/params
index 150a6f0..0a589c3 100644
--- a/testing/const_rate/params
+++ b/testing/const_rate/params
@@ -1,5 +1,9 @@
-# velocity	spring constant	unfolding rate
-1	1	1
-2	1	1
-1	2	1
-1	1	2
+# num domains	velocity	spring constant	unfolding rate
+1	1	1	1
+1	2	1	1
+1	1	2	1
+1	1	1	2
+50	1	1	1
+50	2	1	1
+50	1	2	1
+50	1	1	2