testing/kramers_rate/run_test.sh

   1 #!/bin/bash
   2 # -*- coding: utf-8 -*-
   3 #
   4 # Test the Kramers integral k(F) implementation on the cusp potential.
   5 #
   6 # The test is escape from a cusp potential from Kramers' 1940 paper.
   7 # He solves the cusp potential in the text immediately after Eq. 17 & Fig. 2
   8 # We can translate Kramers' notation to Hänggi's more modern notation with
   9 #  Kramers    Hängii    Meaning               Why I think so
  10 #  2πω        ω         Curvature at minimum  Hänggi's footnote to Eq. 1.5
  11 #  η          γ         Viscosity/damping     Compare H. Eq.1.5  with K. Eq. 25
  12 #  Q          E_b       Energy barrier height K. Fig. 1 (p 291)
  13 #  T          k_B·T≡1/β Absolute temperature  K: “k_B = 1” (text after Eq. 5)
  14 #  r          k         Reaction rate         seems more reasonable than k ;)
  15 # With these translations, Kramers' Eq. is
  16 #  k ≅ ω²/2πγ √(πβE_b) e^(-βE_b)
  17 # for E(x) = ½ω²(|x|-Δx)².
  18 # The barrier is at x_b=0 with minimums at x_ac=±Δx.
  19 # The barrier height E_b = E(0) = ½ω²Δx²
  20 #
  21 # As a force is applied to tip the potential,
  22 #   E(F,x) = ½ω²(|x|-Δx)²-Fx.
  23 # Focusing only on the area left of the barrier (x<0), we have
  24 #   E_L(F,x) = ½ω²(-x-Δx)²-Fx =  ½ω²(x+Δx)²-Fx
  25 #   E_L'(F,x) = ω²(x+Δx)-F
  26 #   E_L"(F,x) = ω²,
  27 # just like the F=0 case.  Kramers' solution still holds then, but with a
  28 # different barrier height.  The minimum moves to
  29 #   x_a = F/ω² - Δx.
  30 # so
  31 #   Eb(F) = E(F,0) - E(F,x_a) = E(0,0) - E(F,x_a)
  32 #         = ½ω²Δx² - ½ω²[(F/ω² - Δx)+Δx]²-F(F/ω² - Δx)
  33 #         = ½ω²Δx² - ½F² - F²/ω² - FΔx
  34 #         = ½ω²Δx² - FΔx - (½+1/ω²)F²
  35 #            ----+-----       ^- reduction due to decreasing barrier distance
  36 #                ^- Bell model
  37 #
  38 # Usage: kramers_rate.sh
  39
  40 # Since we're on the cluster, we'll go crazy and get really good statistics ;)
  41 FMAX=1
  42 K_UTILS=../../bin/k_model_utils
  43 DATA=kramers_rate.d
  44 HIST=kramers_rate.hist
  45 GPFILE=kramers_rate.gp
  46 PNGFILE=kramers_rate.png
  47
  48 # set up the physics
  49 #   E(x) = ½ω²(|x|-Δx)²
  50 OMEGA=3
  51 GAMMA=100
  52 KB=1.3806503e-23 # m²·kg/s²·K
  53 T=`python -c "print (1 / ($KB))"`  # pick units where k_B = 1
  54 DX=1
  55 EB=`python -c "print (0.5 * (($OMEGA)*($DX))**2)"`
  56
  57 # set up the theory
  58 E_FN="0.5*($OMEGA)**2 * (abs(x)-($DX))**2"
  59 BETA=`python -c "print (1.0 / (($KB)*($T)))"`
  60 #  k ≅ ω²/2πγ √(πβE_b) e^(-βE_b)
  61 BETAEB=`python -c "print (($BETA)*($EB))"`
  62 THEORY_K="BETAEB(F) = ($BETA)*(($EB) - (0.5+1/($OMEGA)**2)*F**2 - F*($DX))
  63 K(F) = ($OMEGA)**2/(2*pi*($GAMMA)) * sqrt(pi*BETAEB(F)) * exp(-BETAEB(F))"
  64
  65 # Get the diffusion coefficient from the Einstein-Smoluchowski relation
  66 #   D = µ_p·k_B·T = k_B·T/γ
  67 # where µ_p=1/γ is the particle mobility
  68 # See, for example,
  69 #  http://en.wikipedia.org/wiki/Einstein_relation_(kinetic_theory)
  70 D=`python -c "print (1/(($BETA)*($GAMMA)))"`
  71 # pick safe integration bounds
  72 XMIN=`python -c "print (-2*($DX))"`
  73 XMAX="$DX"
  74 # set up the cubic-spline approximation to a sharp cusp… get close anyway…
  75 SPLINE_PTS=`python -c "from sys import stdout
  76 dx = (($XMAX)-($XMIN))/100.0
  77 stdout.write('{')
  78 for i in range(100) :
  79     x = ($XMIN)+i*dx
  80     stdout.write('(%g,%g),' % (x, ($E_FN)))
  81 stdout.write('(%g,%g)}\n' % (($XMAX), ($E_FN)))"`
  82
  83 # Run the simulation
  84 echo "$K_UTILS -kkramers_integ -K \"$D,$SPLINE_PTS,$XMIN,$XMAX\" -x$XMIN -X$XMAX -T$T -F$FMAX > $DATA"
  85 $K_UTILS -kkramers_integ -K "$D,$SPLINE_PTS,$XMIN,$XMAX" -x$XMIN -X$XMAX -T$T -F$FMAX > $DATA
  86
  87 # fit the data
  88 FIT=`cat $DATA | ./fit_kramers`
  89 FITGAMMA=`echo "$FIT" | sed -n 's/gamma:[[:space:]]*//p'`
  90 FITOMEGA=`echo "$FIT" | sed -n 's/omega:[[:space:]]*//p'`
  91 FITT=`echo "$FIT" | sed -n 's/temp:[[:space:]]*//p'`
  92 FITDX=`echo "$FIT" | sed -n 's/delta x:[[:space:]]*//p'`
  93 echo -e "gamma:\t$FITGAMMA\t(expected $GAMMA)"
  94 echo -e "omega:\t$FITOMEGA\t(expected $OMEGA)"
  95 echo -e "temp:\t$FITT\t(expected $T)"
  96 echo -e "delta x:\t$FITDX\t(expected $DX)"
  97
  98 FITEB=`python -c "print (0.5 * (($FITOMEGA)*($FITDX))**2)"`
  99 FITE_FN="0.5*($FITOMEGA)**2 * (abs(x)-($FITDX))**2"
 100 FITBETA=`python -c "print (1.0 / (($KB)*($FITT)))"`
 101 FITBETAEB=`python -c "print (($FITBETA)*($FITEB))"`
 102 FIT_K="FITBETAEB(F) = ($FITBETA)*(($FITEB) - (0.5+1/($FITOMEGA)**2)*F**2 - F*($FITDX))
 103 fitK(F) = ($FITOMEGA)**2/(2*pi*($FITGAMMA)) * sqrt(pi*FITBETAEB(F)) * exp(-FITBETAEB(F))"
 104
 105 # generate gnuplot script
 106 echo "set terminal png" > $GPFILE
 107 echo "set output '$PNGFILE'" >> $GPFILE
 108 echo "set logscale y" >> $GPFILE
 109 echo "set xlabel 'Time'" >> $GPFILE
 110 echo "set ylabel 'Deaths'" >> $GPFILE
 111 echo "$THEORY_K" >> $GPFILE
 112 echo "$FIT_K" >> $GPFILE
 113 echo "plot K(x), fitK(x), '$DATA' with points" >> $GPFILE
 114
 115 gnuplot $GPFILE