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 Analytically, 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 :math:`K`, the population follows
46 p(t) = exp(-tK) = exp(-(f-f_0)K/kv) = p(f) \\;.
48 Therefore, a histogram of unfolding vs. force :math:`p(f)` normalized
49 to :math:`p(0)=1` should follow
52 -w N dp/dF = w N pK/kv = w N K/kv exp(-FK/kv)
54 where :math:`w` is the binwidth and :math:`N` is the number of tested
60 from numpy import arange, exp, log
62 from ..histogram import Histogram
63 from ..fit import HistogramModelFitter
64 from ..manager import get_manager
65 from ..sawsim import SawsimRunner
66 from ..sawsim_histogram import sawsim_histogram
69 def probability_distribution(x, params):
70 """Exponential decay decay probability distribution.
72 .. math:: 1/\\tau \\cdot exp(-x/\\tau)
74 p = params # convenient alias
75 p[0] = abs(p[0]) # cannot normalize negative tau.
76 return (1.0/p[0]) * exp(-x/p[0])
79 class ExponentialModelFitter (HistogramModelFitter):
80 """Exponential decay model fitter.
82 def model(self, params):
83 """An exponential decay model.
87 .. math:: y \\propto \\cdot e^{-t/\\tau}
89 self._model_data.counts = (
90 self.info['binwidth']*self.info['N']*probability_distribution(
91 self._model_data.bin_centers, params))
92 return self._model_data
94 def guess_initial_params(self):
95 return [self._data.mean]
97 def guess_scale(self, params):
98 return [self._data.mean]
100 def model_string(self):
101 return 'p(x) = A exp(-t/tau)'
103 def param_string(self, params):
105 for name,param in zip(['tau'], params):
106 pstrings.append('%s=%g' % (name, param))
107 return ', '.join(pstrings)
111 def constant_rate(sawsim_runner, num_domains=1, unfolding_rate=1,
112 spring_constant=1, velocity=1, N=100):
113 loading_rate = float(spring_constant * velocity)
114 tau = loading_rate / unfolding_rate
115 w = 0.2 * tau # calculate bin width (in force)
116 A = w*num_domains*N / tau
118 # A exp(-x/tau) = 0.001
119 # -x/tau = log(0.001/A)
120 # x = -log(0.001/A)*tau = log(A/0.001)*tau = log(1e3*A)*tau
121 theory.bin_edges = arange(start=0, stop=log(1e3*A)*tau, step=w)
122 theory.bin_centers = theory.bin_edges[:-1] + w/2
123 theory.counts = w*num_domains*N*probability_distribution(
124 theory.bin_centers, [tau])
127 max_force_step = w/10.0
128 max_time_step = max_force_step / loading_rate
130 '-d %(max_time_step)g -F %(max_force_step)g -v %(velocity)g '
131 '-s cantilever,hooke,%(spring_constant)g -N1 '
132 '-s folded,null -N%(num_domains)d -s unfolded,null '
133 '-k folded,unfolded,const,%(unfolding_rate)g -q folded '
135 sim = sawsim_histogram(sawsim_runner, param_string, N=N,
136 bin_edges=theory.bin_edges)
138 e = ExponentialModelFitter(sim)
140 sim_tau = abs(params[0])
141 for s,t,n in [(sim_tau, tau, 'tau')]:
142 assert abs(s - t)/w < 3, (
143 'simulation %s = %g != %g = %s (bin width = %g)' % (n,s,t,n,w))
144 return sim.residual(theory)
147 def test(num_domains=[1,50], unfolding_rate=[1,2], spring_constant=[1,2],
148 velocity=[1,2], threshold=0.2):
150 sr = SawsimRunner(manager=m)
152 for nd in num_domains:
153 for ur in unfolding_rate:
154 for sc in spring_constant:
156 residual = constant_rate(
157 sawsim_runner=sr, num_domains=nd,
158 unfolding_rate=ur, spring_constant=sc, velocity=v)
159 assert residual < threshold, residual