testing/bell_rate/bell_rate.sh

   1 #!/bin/bash
   2 #
   3 # Use a single hooke domain with a constant unfolding rate, so
   4 # unfolding force is proportional to time, and connect multiple
   5 # identical Bell unfolding domains.  We expect an peaked death
   6 # histogram.
   7 #
   8 # e.g.
   9 #   spring constant k=df/dx    \__
  10 #   velocity        v=dx/dt    /  `--> df/dt = kv --> f = kvt  (+ f(t=0))
  11 #   unfolding rate  K= a e^(bf)   (frac/s)
  12 #   let f(t=0) = 0
  13 #   let P(f) be the population existing at force f (normed so P(t=f=0)=1).
  14 #   the probability of dying (normed # deaths) between f and f+df is given by
  15 #        dD = K P(f) * dt = K P(f) df/kv     since dt = df/kv
  16 #   this is what we plot in the normalized histogram with dD in a bin of width
  17 #   df.  The continuous analog is the function dD/df = K/kv P(f).
  18 #
  19 #   To find P(f), we need to solve
  20 #        dP/dt = -KP
  21 #        dP/P = -K dt = -a e^{bf} df/kv = -z e^{bf} df       where z ≡ a/kv
  22 #        ln(P) = -z/b e^{bf) + c
  23 #        P = C e^{-z/b e^{bf}}
  24 #   Applying our f(t=0)=0 boundary condition we see that C = P(t=f=0) e^{z/b}.
  25 #   Plugging P(f) into our formula for dD/df yields
  26 #        dD/df = K/kv P(f) = z e^{bf}  P(t=f=0) e^{z/b} e^{-z/b e^{bf}}
  27 #        dD/df = P(t=f=0) z e^{bf + z/b(1-e^{bf})}
  28 #
  29 #   The histogram scaling isn't dD/df, though, it's dD/d(bin)
  30 #        dD/d(bin) = dD/df df/d(bin)
  31 #                  = binwidth P(t=f=0) z e^{bf + z/b(1-e^{bf})}
  32 #
  33 # Reminders:
  34 #   spring constant k=df/dx       (N/m)
  35 #   velocity        v=dx/dt       (m/s)
  36 #   unfolding rate  K= a e^(bf) = K0 e^(fdx/kBT)  (frac/s)
  37 # usage: bell_rate.sh <num_domains> <v> <k> <K0> <dx> <T>
  38
  39 if [ $# -ne 6 ]
  40 then
  41     echo "usage: $0 <num_domains> <v> <k> <K0> <dx> <T>"
  42     echo ""
  43     echo "where"
  44     echo "spring constant k=df/dx             (N/m)"
  45     echo "velocity        v=dx/dt             (m/s)"
  46     echo "unfolding rate  K= K0 e^(fdx/kBT)   (frac/s)"
  47     exit 1
  48 fi
  49
  50 NDOM=$1
  51 V=$2
  52 SPRING=$3
  53 K0=$4
  54 DX=$5
  55 T=$6
  56 DEBUG=0
  57
  58 # since we're on the cluster, we'll go crazy and get really good statistics ;)
  59 N=10000
  60
  61 SAWSIM=../../bin/sawsim
  62 DATA=bell_rate_${NDOM}_${V}_${SPRING}_${K0}_${DX}_${T}.d
  63 HIST=bell_rate_${NDOM}_${V}_${SPRING}_${K0}_${DX}_${T}.hist
  64 COMPFILE=bell_rate_${NDOM}_${V}_${SPRING}_${K0}_${DX}_${T}.compare
  65 GPFILE=bell_rate_${NDOM}_${V}_${SPRING}_${K0}_${DX}_${T}.gp
  66 PNGFILE=bell_rate_${NDOM}_${V}_${SPRING}_${K0}_${DX}_${T}.png
  67
  68 # set up the theory
  69 # our predicted histogram (see notes above) is
  70 #   dD/d(bin) = binwidth P0 z e^{bf + z/b(1-e^{bf})}
  71 # where z ≡ a/kv, and b = Δx/kBT
  72 KB=1.3806503e-23 # J/K
  73 B=`python -c "print ($DX/($KB*$T))"`
  74 Z=`python -c "print (float($K0)/($SPRING*$V))"`
  75 # find a good df (histogram bin width).
  76 #    currently, scale based only on e^(bf) term in exponent.
  77 df=`python -c "print (0.01/$B)"`
  78
  79 THEORY_FN="$df * $N * $NDOM * $Z * exp($B*x + $Z/$B *(1 - exp($B*x)))"
  80 if [ $DEBUG -eq 1 ]; then
  81     echo "b = dx/kBT = $B"
  82     echo "z = K0/kv = $Z"
  83     echo "Theory --> P(x=f) = $THEORY_FN"
  84 fi
  85
  86 # run the experiments
  87 > $DATA
  88 if [ "$ISACLUSTER" -eq 1 ]
  89 then
  90     # Sawsim >= 0.6
  91     qcmd -w14:00:00 xtimes $N $SAWSIM -T$T -v$V -s cantilever,hooke,$SPRING -N1 \
  92         -s folded,null -N$NDOM -s unfolded,null \
  93         -k "'folded,unfolded,bell,{$K0,$DX}'" -q folded \
  94         | grep -v '^#' | cut -f1 >> $DATA
  95 else
  96     # Sawsim >= 0.6
  97     xtimes $N $SAWSIM -T$T -v$V -s cantilever,hooke,$SPRING -N1 \
  98         -s folded,null -N$NDOM -s unfolded,null \
  99         -k "folded,unfolded,bell,{$K0,$DX}" -q folded \
 100         | grep -v '^#' | cut -f1 >> $DATA
 101 fi
 102
 103 # histogram the data
 104 if [ $DEBUG -eq 1 ]; then
 105     echo "cat $DATA | stem_leaf -v-1 -b$df -c | sed 's/://' > $HIST"
 106 fi
 107 cat $DATA | stem_leaf -v-1 -b$df -c | sed 's/://' > $HIST
 108 wc -l "$HIST"
 109 # fit the data to an exponential decay
 110 if [ $DEBUG -eq 1 ]; then
 111     echo "cat $HIST | ./fit_bell"
 112 fi
 113 FIT=`cat $HIST | ./fit_bell` || exit 1
 114 FITZ=`echo "$FIT" | sed -n 's/a:[[:space:]]*//p'`
 115 FITB=`echo "$FIT" | sed -n 's/b:[[:space:]]*//p'`
 116 echo -e "z = K0/kv:\t$FITZ\t(expected $Z)" > $COMPFILE
 117 echo -e "b = dx/kBT:\t$FITB\t(expected $B)" >> $COMPFILE
 118 cat $COMPFILE
 119
 120 # generate gnuplot script
 121 if [ $DEBUG -eq 1 ]; then
 122     echo "generate Gnuplot script"
 123 fi
 124 echo "set terminal png" > $GPFILE
 125 echo "set output '$PNGFILE'" >> $GPFILE
 126 echo "set logscale y" >> $GPFILE
 127 echo "set xlabel 'Force'" >> $GPFILE
 128 echo "set ylabel 'Deaths'" >> $GPFILE
 129 echo "theory(x) = $THEORY_FN" >> $GPFILE
 130 echo "fit(x) = $df * $N *  $NDOM * $FITZ * exp($FITB*x + $FITZ/$FITB*(1-exp($FITB*x)))" >> $GPFILE
 131 echo "plot  theory(x), '$HIST' with points, fit(x)" >> $GPFILE
 132
 133 if [ $DEBUG -eq 1 ]; then
 134     echo "Generate PNG"
 135 fi
 136 gnuplot $GPFILE