testing/const_rate/const_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 constant unfolding domains.  We expect an exponential
   6 # population decay.
   7 #
   8 # e.g.
   9 #   spring constant k=df/dx      \__
  10 #   velocity        v=dx/dt      /  `--> df/dt = kv --> f = kvt  (+ f0)
  11 #   unfolding rate  K=1 frac/s
  12 #
  13 #   Population as a function of force/time follows
  14 #     p(t) = exp(-tK) = exp(-FK/kv) = p(F)
  15 #   so histogram of unfolding vs. force p(F) normalized to p(0)=1 should follow
  16 #     -binwidth N0 dp/dF = binwidth N0 pK/kv = binwidth N0 K/kv exp(-FK/kv)
  17 #
  18 # usage: const_rate.sh <num_domains> <velocity> <spring_constant> <unfolding_rate>
  19
  20 if [ $# -ne 4 ]
  21 then
  22     echo "usage: $0 <num_domains> <velocity> <spring_constant> <unfolding_rate>"
  23     exit 1
  24 fi
  25
  26 NDOM=$1
  27 V=$2
  28 K=$3
  29 RATE=$4
  30 DEBUG=0
  31
  32 # since we're on the cluster, we'll go crazy and get really good statistics ;)
  33 N=10000
  34 df=`python -c "print 0.1/$RATE"`
  35
  36 SAWSIM=../../bin/sawsim
  37 DATA=const_rate_${NDOM}_${V}_${K}_${RATE}.d
  38 HIST=const_rate_${NDOM}_${V}_${K}_${RATE}.hist
  39 COMPFILE=const_rate_${NDOM}_${V}_${K}_${RATE}.compare
  40 GPFILE=const_rate_${NDOM}_${V}_${K}_${RATE}.gp
  41 PNGFILE=const_rate_${NDOM}_${V}_${K}_${RATE}.png
  42
  43 # set up the theory
  44 # our predicted histogram (see notes above) is
  45 #  hist(F) = binwidth N0 K/kv exp(-FK/kv) = A exp(F/B)
  46 theoryA=`python -c "print $df*$N*$NDOM*$RATE/float($K*$V)"`
  47 theoryB=`python -c "print float(-$K*$V)/$RATE"`
  48 THEORY_FN="$theoryA * exp(x / ($theoryB))"
  49
  50
  51 # run the experiments
  52 > $DATA
  53 if [ "$ISACLUSTER" -eq 1 ]
  54 then
  55     # Sawsim >= 0.6
  56     qcmd xtimes $N $SAWSIM -v$V -s cantilever,hooke,$K -N1 \
  57         -s folded,null -N$NDOM -s unfolded,null \
  58         -k folded,unfolded,const,$RATE -q folded \
  59         | grep -v '^#' | cut -f1 >> $DATA
  60 else
  61     # Sawsim >= 0.6
  62     xtimes $N $SAWSIM -v$V -s cantilever,hooke,$K -N1 \
  63         -s folded,null -N$NDOM -s unfolded,null \
  64         -k folded,unfolded,const,$RATE -q folded \
  65         | grep -v '^#' | cut -f1 >> $DATA
  66 fi
  67
  68 # histogram the data
  69 if [ $DEBUG -eq 1 ]; then
  70     echo "cat $DATA | stem_leaf -v-1 -b$df -c | sed 's/://' > $HIST"
  71 fi
  72 cat $DATA | stem_leaf -v-1 -b$df -c | sed 's/://' > $HIST
  73
  74 # fit the data to an exponential decay
  75 if [ $DEBUG -eq 1 ]; then
  76     echo "cat $HIST | ./fit_exponential"
  77 fi
  78 FIT=`cat $HIST | ./fit_exponential` || exit 1
  79 B=`echo "$FIT" | sed -n 's/time constant:[[:space:]]*//p'`
  80 A=`echo "$FIT" | sed -n 's/initial value:[[:space:]]*//p'`
  81 echo -e "force scale:\t$B\t(expected $theoryB)" > $COMPFILE
  82 echo -e "initial value:\t$A\t(expected $theoryA)" >> $COMPFILE
  83 cat $COMPFILE
  84
  85 # generate gnuplot script
  86 if [ $DEBUG -eq 1 ]; then
  87     echo "generate Gnuplot script"
  88 fi
  89 echo "set terminal png" > $GPFILE
  90 echo "set output '$PNGFILE'" >> $GPFILE
  91 echo "set logscale y" >> $GPFILE
  92 echo "set xlabel 'Force'" >> $GPFILE
  93 echo "set ylabel 'Deaths'" >> $GPFILE
  94 echo "theory(x) = $THEORY_FN" >> $GPFILE
  95 echo "fit(x) = $A * exp(x / $B)" >> $GPFILE
  96 echo "plot  theory(x), '$HIST' with points, fit(x)" >> $GPFILE
  97
  98 if [ $DEBUG -eq 1 ]; then
  99     echo "Generate PNG"
 100 fi
 101 gnuplot $GPFILE