1
2
3
4
5
6
7
8
9
10
11
12
13
14
15 __id__ = '$Id: debyefitting.py 6718 2011-08-23 21:33:20Z yshang $'
16
18 '''Class to perform Debye fitting
19
20 This fits the isotropic ADPs (Uiso) of a specified atom in a material
21 according to the Debye model.
22
23 Attributes
24 x -- Temperature values (default none).
25 y -- ADP values of the atom at temperatures specfied by x
26 dy -- Uncertainty in y (default None).
27 mass -- The mass of the atom in atomic mass units (e.g., 12 for carbon)
28 _residual -- The residual function to be optimized.
29 results -- The full results returned from the optimizer (using
30 full_results = True). See the leastsq
31 function from scipy.optimize for the format.
32 thetaD -- The Debye temperature extracted from the ADPs by fitting the
33 Debye model.
34 offset -- The static offset in the ADPs extracted by fitting the Debye
35 model.
36
37 '''
38
39 - def __init__(self, x, y, mass, dy = None):
40 '''Constructor of the fitting object
41
42 x -- A numpy array of temperatures.
43 y -- The thermal parameters over the calculation range.
44 mass -- The mass of the atom in atomic mass units.
45 dy -- The uncertainty in the thermal parameters over the
46 calculation range. If this is None (default), then the
47 uncertainties are set to unity.
48
49 '''
50 self.x = x
51 self.y = y
52 self.dy = dy
53 self.mass = mass
54 self._residual = None
55 self.results = None
56 self.thetaD = None
57 self.offset = None
58 self.yfit = None
59
60 if self.dy is None:
61 import numpy
62 self.dy = numpy.ones_like(self.x)
63 return
64
66 """Make the vector residual function to be optimized."""
67 if self.x is None or self.y is None or self.mass is None:
68 m = "Temperature, ADPs and mass must be specified"
69 raise AttributeError(m)
70
71 def residual(*pars):
72
73 self.thetaD, self.offset = pars[0]
74 self.yfit = debye(self.x, self.mass, self.thetaD)
75 self.yfit += self.offset
76
77 chiv = (self.y - self.yfit) / self.dy
78
79 return chiv
80
81 self._residual = residual
82
83 return
84
86 """Optimize the recipe created above using scipy.
87
88 The FitRecipe we created in makeRecipe has a 'residual' method that we
89 can be minimized using a scipy optimizer. The details are described in
90 the source.
91
92 """
93 from scipy.optimize.minpack import leastsq
94
95 self.results = leastsq(self._residual, [100, 0], full_output = True)
96
97 self.thetaD, self.offset = self.results[0]
98
99 def _calc_error(args):
100 """Get the results and covarience of the fitted parameters"""
101 from numpy import array
102 from math import sqrt
103 p, cov, info, mesg, success = args
104 chisq=sum(info["fvec"]*info["fvec"])
105 dof=len(info["fvec"])-len(p)
106 sigma = array([sqrt(cov[i,i])*sqrt(chisq/dof) for i in range(len(p))])
107 return p, sigma
108
109 params, sig = _calc_error(self.results)
110
111 [self.thetaD, self.offset], [self.cov_thetaD, self.cov_offset] = params, sig
112
113 return
114
116 """Create the recipe, fit it and get the results."""
117 self._makeResidual()
118 self._optimize()
119 return
120
122 """Save the results to a text file"""
123 content = []
124 content.append('# ThetaD %e / %e ' % (self.thetaD, self.cov_thetaD))
125 content.append('# offset %e / %e ' % (self.offset, self.cov_offset))
126 content.append('x\ty\tyfit\tdy')
127 for ii in range(len(self.x)):
128 str = '%f\t%f\t%f\t%f' % (self.x[ii], self.y[ii],self.yfit[ii], self.dy[ii])
129 content.append(str)
130 content = '\n'.join(content)
131
132 f = open(filename, 'w')
133 f.write(content)
134
135 return
136
137
139 """Plot the results contained within a refined FitRecipe."""
140 x = self.x
141 y = self.y
142
143
144
145 import pylab
146 pylab.plot(x, y, 'b.', label = "data")
147 ylabel = "${\Theta}_D = %.2f K, offset = %.2f \AA^2$" %\
148 (self.thetaD, self.offset)
149 pylab.plot(x, self.yfit, 'g-', label = ylabel)
150 pylab.legend(loc = (0.0,0.8))
151 pylab.xlabel("Temperature (K)")
152 pylab.ylabel("$U_{iso} (\AA^2)$")
153 pylab.show()
154 return
155
156
157
158 import numpy
159
161 """A wrapped version of 'adps' that can handle an array of T-values."""
162 y = numpy.array([adps(m, thetaD, x) for x in T])
163 return y
164
165 -def adps(m, thetaD, T):
166 """Calculates atomic displacement factors within the Debye model
167
168 <u^2> = (3h^2/4 pi^2 m kB thetaD)(phi(thetaD/T)/(thetaD/T) + 1/4)
169
170 arguments:
171 m -- float -- mass of the ion in atomic mass units (e.g., C = 12)
172 thetaD -- float -- Debye Temperature
173 T -- float -- temperature.
174
175 return:
176 Uiso -- float -- the thermal factor from the Debye recipe at temp T
177
178 """
179 h = 6.6260755e-34
180 kB = 1.3806503e-23
181 amu = 1.66053886e-27
182
183 def __phi(x):
184 """evaluates the phi integral needed in Debye calculation
185
186 phi(x) = (1/x) int_0^x xi/(exp(xi)-1) dxi
187
188 arguments:
189 x -- float -- value of thetaD (Debye temperature)/T
190
191 returns:
192 phi -- float -- value of the phi function
193
194 """
195 def __debyeKernel(xi):
196 """function needed by debye calculators
197
198 """
199 y = xi/(numpy.exp(xi)-1)
200 return y
201
202 import scipy.integrate
203
204 int = scipy.integrate.quad(__debyeKernel, 0, x)
205 phi = (1/x) * int[0]
206
207 return phi
208
209
210 m = m * amu
211 u2 = (3*h**2 / (4 * numpy.pi**2 *m *kB *thetaD))*\
212 (__phi(thetaD/T)/(thetaD/T) + 1./4.)
213 u2 *= 1e20
214 return u2
215