testing/bell_rate/run_test.sh

   1 #!/bin/bash
   2 #
   3 # Use a single hooke domain with a bell unfolding rate, so unfolding force
   4 # is proportional to time, and we expect a peaked death histogram.
   5 #
   6 # e.g.
   7 #   spring constant k=df/dx=1    \__
   8 #   velocity        v=dx/dt=1    /  `--> df/dt = kv = 1  --> f = t  (+ f0)
   9 #   unfolding rate  K= a e^(bf)    (frac/s)
  10 #   let f0 = 0
  11 #   let P(f) be the population existing at force f (normed so P(0)=0).
  12 #   the probability of dying (normed # deaths) between f and f+df is given by
  13 #        dD = K P(f) * dt = K P(f) df/kv     since dt = df/kv
  14 #   this is what whe plot in the normalized histogram with dD in a bin of width
  15 #   df.  The continuous analog is the function dD/df = K/kv P(f).
  16 #
  17 #        dP/dt = - KP
  18 #        dP/P = -K dt = -a e^{bf} df/kv = -z e^{bf} df       where z ≡ a/kv
  19 #        ln(P) = -z/b e^{bf) + c
  20 #        P = C e^{-z/b e^{bf}}
  21 #   Applying our f0=0 boundary condition we see that C = P0 e^{z/b}.
  22 #   Plugging P(f) into our formula for dD/df yields
  23 #        dD/df = K/kv P(f) = z e^{bf}  P0 e^{z/b} e^{-z/b e^{bf}}
  24 #        dD/df = P0 z e^{bf + z/b(1-e^{bf})}
  25 #
  26 #   The histogram scaling isn't dD/df, though, it's dD/d(bin)
  27 #        dD/d(bin) = dD/df df/d(bin) = binwidth P0 z e^{bf + z/b(1-e^{bf})}
  28 #
  29 #   test with lots of runs of
  30 #              v- T=1/kB    v- v=1            v- k=1
  31 #     sawsim -T7.24296e22 -v1 -s folded,hooke,1 -N1 -s unfolded,null \
  32 #         -k "folded,unfolded,bell,{1,1}" -q folded
  33 #                                   ^- {a,dx}={1,1}
  34 #   (kB = 1.3806503 10-23 J/K, we use T=1/kB to get the bell K=e^{-f})
  35 #
  36 #   60 runs of that takes ~1 second:
  37 #     time xtimes 60 sawsim -T7.24296e22 -v1 -s folded,hooke,1 -N1 \
  38 #         -s unfolded,null -k "folded,unfolded,bell,{1,1}" -q folded \
  39 #         > /dev/null
  40 #
  41 # usage: bell_rate.sh
  42
  43 # since we're on the cluster, we'll go crazy and get really good statistics ;)
  44 N=10000
  45
  46 SAWSIM=../../bin/sawsim
  47 DATA=bell_rate.d
  48 HIST=bell_rate.hist
  49 GPFILE=bell_rate.gp
  50 PNGFILE=bell_rate.png
  51
  52 # set up the theory
  53 #   dD/d(bin) = binwidth P0 z e^{bf + z/b(1-e^{bf})}
  54 # where z ≡ a/kv, and b = Δx/kBT
  55 KB=1.3806503e-23 # J/K
  56 #STYLE="ones"
  57 STYLE="protein"
  58 if [ "$STYLE" == "ones" ]
  59 then
  60     df=1e-1  # histogram bin width
  61     DX=1
  62     T=7.24296e22 # 1/kB
  63     A=1
  64     K=1
  65     V=1
  66 else
  67     df=1e-12
  68     DX=1e-9
  69     T=300
  70     A=1e-3
  71     K=0.05
  72     V=1e-6
  73 fi
  74 B=`echo "print ($DX/($KB*$T))" | python`
  75 Z=`echo "print (float($A)/($K*$V))" | python`
  76 THEORY_FN="$df * $N * $Z * exp($B*x + $Z/$B *(1 - exp($B*x)))"
  77
  78 # run the experiments
  79 > $DATA
  80 if [ "$ISACLUSTER" -eq 1 ]
  81 then
  82     # Sawsim <= 0.5
  83     #qcmd "xtimes $N $SAWSIM -T$T -v$V -mhooke -a$K -kbell -K$A,$DX -F1 | grep -v '^#' >> $DATA"
  84     # Sawsim >= 0.6
  85     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"
  86 else
  87     # Sawsim <= 0.5
  88     #xtimes $N $SAWSIM -T$T -v$V -mhooke -a$K -kbell -K$A,$DX -F1 | grep -v '^#' >> $DATA
  89     # Sawsim >= 0.6
  90     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
  91 fi
  92
  93 # histogram the data
  94 cat $DATA | stem_leaf -v-1 -b$df | sed 's/://' > $HIST
  95
  96 # fit the data
  97 FIT=`cat $HIST | ./fit_bell`
  98 #echo "$FIT"
  99 FITZ=`echo "$FIT" | sed -n 's/a:[[:space:]]*//p'`
 100 FITB=`echo "$FIT" | sed -n 's/b:[[:space:]]*//p'`
 101 echo -e "z:\t$FITZ\t(expected $Z)"
 102 echo -e "b:\t$FITB\t(expected $B)"
 103
 104 # generate gnuplot script
 105 echo "set terminal png" > $GPFILE
 106 echo "set output '$PNGFILE'" >> $GPFILE
 107 echo "set logscale y" >> $GPFILE
 108 echo "set xlabel 'Time'" >> $GPFILE
 109 echo "set ylabel 'Deaths'" >> $GPFILE
 110 echo "theory(x) = $THEORY_FN" >> $GPFILE
 111 echo "fit(x) = $df * $N * $FITZ * exp($FITB*x + $FITZ/$FITB*(1-exp($FITB*x)))" >> $GPFILE
 112 echo "plot  theory(x), '$HIST' with points, fit(x)" >> $GPFILE
 113
 114 gnuplot $GPFILE