1 # Copyright (C) 2008-2010 W. Trevor King <wking@drexel.edu>
3 # This program is free software: you can redistribute it and/or modify
4 # it under the terms of the GNU General Public License as published by
5 # the Free Software Foundation, either version 3 of the License, or
6 # (at your option) any later version.
8 # This program is distributed in the hope that it will be useful,
9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 # GNU General Public License for more details.
13 # You should have received a copy of the GNU General Public License
14 # along with this program. If not, see <http://www.gnu.org/licenses/>.
16 # The author may be contacted at <wking@drexel.edu> on the Internet, or
17 # write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
18 # Philadelphia PA 19104, USA.
20 """Test a Hookian chain with constant unfolding rate.
22 With the constant velocity experiment and a Hookian domain, the
23 unfolding force is proportional to time, so we expect exponential
24 population decay, even with multiple unfolding domains.
26 For example, with a spring constant
30 and a pulling velocity
34 we have a loading rate
36 .. math:: df/dt = df/dx dx/dt = kv \;,
40 .. math:: f = kvt + f_0 \;.
42 With an unfolding rate constant
44 .. math:: K = 1 \text{frac/s}
46 The population follows
49 dp/dt = Kp = 1 \text{s^{-1}}
50 p(t) = exp(-tK) = exp(-(f-f_0)K/kv) = p(f) \;.
52 Therefore, a histogram of unfolding vs. force :math:`p(f)` normalized
53 to :math:`p(0)=1` should follow
56 -w N dp/dF = w N pK/kv = w N K/kv exp(-FK/kv)
58 where :math:`w` is the binwidth and :math:`N` is the number of tested
62 from numpy import arange, exp, log
64 from ..histogram import Histogram
65 from ..fit import HistogramModelFitter
66 from ..manager import get_manager
67 from ..sawsim import SawsimRunner
68 from ..sawsim_histogram import sawsim_histogram
71 def probability_distribution(x, params):
72 """Exponential decay decay probability distribution.
74 .. math:: 1/\\tau \cdot exp(-x/\\tau)
76 p = params # convenient alias
77 p[0] = abs(p[0]) # cannot normalize negative tau.
78 return (1.0/p[0]) * exp(-x/p[0])
81 class ExponentialModelFitter (HistogramModelFitter):
82 """Exponential decay model fitter.
84 def model(self, params):
85 """An exponential decay model.
89 .. math:: y \\propto \cdot e^{-t/\tau}
91 self._model_data.counts = (
92 self.info['binwidth']*self.info['N']*probability_distribution(
93 self._model_data.bin_centers, params))
94 return self._model_data
96 def guess_initial_params(self):
97 return [self._data.mean]
99 def guess_scale(self, params):
100 return [self._data.mean]
102 def model_string(self):
103 return 'p(x) = A exp(-t/tau)'
105 def param_string(self, params):
107 for name,param in zip(['tau'], params):
108 pstrings.append('%s=%g' % (name, param))
109 return ', '.join(pstrings)
113 def constant_rate(sawsim_runner, num_domains=1, unfolding_rate=1,
114 spring_constant=1, velocity=1, N=400):
115 loading_rate = spring_constant * velocity
116 tau = loading_rate / unfolding_rate
117 w = 0.1 * tau # calculate bin width (in force)
118 A = w*num_domains*N / tau
120 # A exp(-x/tau) = 0.001
121 # -x/tau = log(0.001/A)
122 # x = -log(0.001/A)*tau = log(A/0.001)*tau = log(1e3*A)*tau
123 theory.bin_edges = arange(start=0, stop=log(1e3*A)*tau, step=w)
124 theory.bin_centers = theory.bin_edges[:-1] + w/2
125 theory.counts = w*num_domains*N*probability_distribution(
126 theory.bin_centers, [tau])
130 '-v %(velocity)s -s cantilever,hooke,%(spring_constant)s -N1 '
131 '-s folded,null -N%(num_domains)s -s unfolded,null '
132 '-k folded,unfolded,const,%(unfolding_rate)s -q folded '
134 sim = sawsim_histogram(sawsim_runner, param_string, N=N,
135 bin_edges=theory.bin_edges)
137 e = ExponentialModelFitter(sim)
139 sim_tau = abs(params[0])
140 assert (sim_tau - tau)/tau, 'simulation tau = %g != %g = tau' % (
142 return sim.residual(theory)
145 def test(num_domains=[1,50], unfolding_rate=[1,2], spring_constant=[1,2],
146 velocity=[1,2], threshold=0.3):
148 sr = SawsimRunner(manager=m)
150 for nd in num_domains:
151 for ur in unfolding_rate:
152 for sc in spring_constant:
154 residual = constant_rate(
155 sawsim_runner=sr, num_domains=nd,
156 unfolding_rate=ur, spring_constant=sc, velocity=v)
157 assert residual < threshold, residual