1 | """Element: functions for element types |
---|
2 | Copyright: 2008, Robert B. Von Dreele & Brian H. Toby (Argonne National Laboratory) |
---|
3 | """ |
---|
4 | ########### SVN repository information ################### |
---|
5 | # $Date: 2012-02-15 20:36:53 +0000 (Wed, 15 Feb 2012) $ |
---|
6 | # $Author: vondreele $ |
---|
7 | # $Revision: 484 $ |
---|
8 | # $URL: trunk/GSASIIElem.py $ |
---|
9 | # $Id: GSASIIElem.py 484 2012-02-15 20:36:53Z vondreele $ |
---|
10 | ########### SVN repository information ################### |
---|
11 | |
---|
12 | import math |
---|
13 | import os.path |
---|
14 | import GSASIIpath |
---|
15 | import numpy as np |
---|
16 | |
---|
17 | def GetFormFactorCoeff(El): |
---|
18 | """Read X-ray form factor coefficients from `atomdata.asc` file |
---|
19 | |
---|
20 | :param El: element 1-2 character symbol case irrevelant |
---|
21 | :return: `FormFactors`: list of form factor dictionaries |
---|
22 | |
---|
23 | Each X-ray form factor dictionary is: |
---|
24 | |
---|
25 | * `Symbol`: 4 character element symbol with valence (e.g. 'NI+2') |
---|
26 | * `Z`: atomic number |
---|
27 | * `fa`: 4 A coefficients |
---|
28 | * `fb`: 4 B coefficients |
---|
29 | * `fc`: C coefficient |
---|
30 | |
---|
31 | """ |
---|
32 | ElS = El.upper() |
---|
33 | ElS = ElS.rjust(2) |
---|
34 | filename = os.path.join(os.path.split(__file__)[0],'atmdata.dat') |
---|
35 | try: |
---|
36 | FFdata = open(filename,'Ur') |
---|
37 | except: |
---|
38 | print "**** ERROR - File atmdata.dat not found in directory %s" % sys.path[0] |
---|
39 | sys.exit() |
---|
40 | S = '1' |
---|
41 | FormFactors = [] |
---|
42 | while S: |
---|
43 | S = FFdata.readline() |
---|
44 | if S[3:5] == ElS: |
---|
45 | if S[5:6] != '_' and S[8] not in ['N','M']: |
---|
46 | Z=int(S[:2]) |
---|
47 | Symbol = S[3:7].strip() |
---|
48 | S = S[12:] |
---|
49 | fa = (float(S[:7]),float(S[14:21]),float(S[28:35]),float(S[42:49])) |
---|
50 | fb = (float(S[7:14]),float(S[21:28]),float(S[35:42]),float(S[49:56])) |
---|
51 | FormFac = {'Symbol':Symbol,'Z':Z,'fa':fa,'fb':fb,'fc':float(S[56:63])} |
---|
52 | FormFactors.append(FormFac) |
---|
53 | FFdata.close() |
---|
54 | return FormFactors |
---|
55 | |
---|
56 | def GetFFC5(ElSym): |
---|
57 | '''Get 5 term form factor and Compton scattering data |
---|
58 | @param ElSym: str(1-2 character element symbol with proper case); |
---|
59 | @return El: dictionary with 5 term form factor & compton coefficients |
---|
60 | ''' |
---|
61 | import FormFactors as FF |
---|
62 | El = {} |
---|
63 | FF5 = FF.FFac5term[ElSym] |
---|
64 | El['fa'] = FF5[:5] |
---|
65 | El['fc'] = FF5[5] |
---|
66 | El['fb'] = FF5[6:] |
---|
67 | Cmp5 = FF.Compton[ElSym] |
---|
68 | El['cmpz'] = Cmp5[0] |
---|
69 | El['cmpa'] = Cmp5[1:6] |
---|
70 | El['cmpb'] = Cmp5[6:] |
---|
71 | return El |
---|
72 | |
---|
73 | def GetAtomInfo(El): |
---|
74 | |
---|
75 | import ElementTable as ET |
---|
76 | Elements = [] |
---|
77 | for elem in ET.ElTable: |
---|
78 | Elements.append(elem[0][0]) |
---|
79 | if len(El) in [2,4]: |
---|
80 | ElS = El.upper()[:2].rjust(2) |
---|
81 | else: |
---|
82 | ElS = El.upper()[:1].rjust(2) |
---|
83 | filename = os.path.join(os.path.split(__file__)[0],'atmdata.dat') |
---|
84 | try: |
---|
85 | FFdata = open(filename,'Ur') |
---|
86 | except: |
---|
87 | print '**** ERROR - File atmdata.dat not found in directory %s' % sys.path[0] |
---|
88 | sys.exit() |
---|
89 | S = '1' |
---|
90 | AtomInfo = {} |
---|
91 | Isotopes = {} |
---|
92 | Mass = [] |
---|
93 | while S: |
---|
94 | S = FFdata.readline() |
---|
95 | if S[3:5] == ElS: |
---|
96 | if S[5:6] == '_': |
---|
97 | if not Mass: #picks 1st one; natural abundance or 1st isotope |
---|
98 | Mass = float(S[10:19]) |
---|
99 | if S[6] in [' ','1','2','3','4','5','6','7','8','9']: |
---|
100 | isoName = S[6:9] |
---|
101 | if isoName == ' ': |
---|
102 | isoName = 'Nat. Abund.' #natural abundance |
---|
103 | if S[76:78] in ['LS','BW']: #special anomalous scattering length info |
---|
104 | St = [S[10:19],S[19:25],S[25:31],S[31:38],S[38:44],S[44:50], |
---|
105 | S[50:56],S[56:62],S[62:68],S[68:74],] |
---|
106 | Vals = [] |
---|
107 | for item in St: |
---|
108 | if item.strip(): |
---|
109 | Vals.append(float(item.strip())) |
---|
110 | Isotopes[isoName.rstrip()] = Vals |
---|
111 | else: |
---|
112 | Isotopes[isoName.rstrip()] = [float(S[10:19]),float(S[19:25])] |
---|
113 | elif S[5:9] == '_SIZ': |
---|
114 | Z=int(S[:2]) |
---|
115 | Symbol = S[3:5].strip().lower().capitalize() |
---|
116 | Drad = float(S[12:22]) |
---|
117 | Arad = float(S[22:32]) |
---|
118 | Vdrad = float(S[32:38]) |
---|
119 | Color = ET.ElTable[Elements.index(Symbol)][6] |
---|
120 | FFdata.close() |
---|
121 | AtomInfo={'Symbol':Symbol,'Isotopes':Isotopes,'Mass':Mass,'Z':Z,'Drad':Drad,'Arad':Arad,'Vdrad':Vdrad,'Color':Color} |
---|
122 | return AtomInfo |
---|
123 | |
---|
124 | def GetXsectionCoeff(El): |
---|
125 | """Read atom orbital scattering cross sections for fprime calculations via Cromer-Lieberman algorithm |
---|
126 | @param El: 2 character element symbol |
---|
127 | @return: Orbs: list of orbitals each a dictionary with detailed orbital information used by FPcalc |
---|
128 | each dictionary is: |
---|
129 | 'OrbName': Orbital name read from file |
---|
130 | 'IfBe' 0/2 depending on orbital |
---|
131 | 'BindEn': binding energy |
---|
132 | 'BB': BindEn/0.02721 |
---|
133 | 'XSectIP': 5 cross section inflection points |
---|
134 | 'ElEterm': energy correction term |
---|
135 | 'SEdge': absorption edge for orbital |
---|
136 | 'Nval': 10/11 depending on IfBe |
---|
137 | 'LEner': 10/11 values of log(energy) |
---|
138 | 'LXSect': 10/11 values of log(cross section) |
---|
139 | """ |
---|
140 | AU = 2.80022e+7 |
---|
141 | C1 = 0.02721 |
---|
142 | ElS = El.upper() |
---|
143 | ElS = ElS.ljust(2) |
---|
144 | filename = os.path.join(os.path.split(__file__)[0],'Xsect.dat') |
---|
145 | try: |
---|
146 | xsec = open(filename,'Ur') |
---|
147 | except: |
---|
148 | print '**** ERROR - File Xsect.dat not found in directory %s' % sys.path[0] |
---|
149 | sys.exit() |
---|
150 | S = '1' |
---|
151 | Orbs = [] |
---|
152 | while S: |
---|
153 | S = xsec.readline() |
---|
154 | if S[:2] == ElS: |
---|
155 | S = S[:-1]+xsec.readline()[:-1]+xsec.readline() |
---|
156 | OrbName = S[9:14] |
---|
157 | S = S[14:] |
---|
158 | IfBe = int(S[0]) |
---|
159 | S = S[1:] |
---|
160 | val = S.split() |
---|
161 | BindEn = float(val[0]) |
---|
162 | BB = BindEn/C1 |
---|
163 | Orb = {'OrbName':OrbName,'IfBe':IfBe,'BindEn':BindEn,'BB':BB} |
---|
164 | Energy = [] |
---|
165 | XSect = [] |
---|
166 | for i in range(11): |
---|
167 | Energy.append(float(val[2*i+1])) |
---|
168 | XSect.append(float(val[2*i+2])) |
---|
169 | XSecIP = [] |
---|
170 | for i in range(5): XSecIP.append(XSect[i+5]/AU) |
---|
171 | Orb['XSecIP'] = XSecIP |
---|
172 | if IfBe == 0: |
---|
173 | Orb['SEdge'] = XSect[10]/AU |
---|
174 | Nval = 11 |
---|
175 | else: |
---|
176 | Orb['ElEterm'] = XSect[10] |
---|
177 | del Energy[10] |
---|
178 | del XSect[10] |
---|
179 | Nval = 10 |
---|
180 | Orb['SEdge'] = 0.0 |
---|
181 | Orb['Nval'] = Nval |
---|
182 | D = dict(zip(Energy,XSect)) |
---|
183 | Energy.sort() |
---|
184 | X = [] |
---|
185 | for key in Energy: |
---|
186 | X.append(D[key]) |
---|
187 | XSect = X |
---|
188 | LEner = [] |
---|
189 | LXSect = [] |
---|
190 | for i in range(Nval): |
---|
191 | LEner.append(math.log(Energy[i])) |
---|
192 | if XSect[i] > 0.0: |
---|
193 | LXSect.append(math.log(XSect[i])) |
---|
194 | else: |
---|
195 | LXSect.append(0.0) |
---|
196 | Orb['LEner'] = LEner |
---|
197 | Orb['LXSect'] = LXSect |
---|
198 | Orbs.append(Orb) |
---|
199 | xsec.close() |
---|
200 | return Orbs |
---|
201 | |
---|
202 | def GetMagFormFacCoeff(El): |
---|
203 | """Read magnetic form factor data from atomdata.asc file |
---|
204 | @param El: 2 character element symbol |
---|
205 | @return: MagFormFactors: list of all magnetic form factors dictionaries for element El. |
---|
206 | each dictionary contains: |
---|
207 | 'Symbol':Symbol |
---|
208 | 'Z':Z |
---|
209 | 'mfa': 4 MA coefficients |
---|
210 | 'nfa': 4 NA coefficients |
---|
211 | 'mfb': 4 MB coefficients |
---|
212 | 'nfb': 4 NB coefficients |
---|
213 | 'mfc': MC coefficient |
---|
214 | 'nfc': NC coefficient |
---|
215 | """ |
---|
216 | ElS = El.upper() |
---|
217 | ElS = ElS.rjust(2) |
---|
218 | filename = os.path.join(os.path.split(__file__)[0],'atmdata.dat') |
---|
219 | try: |
---|
220 | FFdata = open(filename,'Ur') |
---|
221 | except: |
---|
222 | print '**** ERROR - File atmdata.dat not found in directory %s' % sys.path[0] |
---|
223 | sys.exit() |
---|
224 | S = '1' |
---|
225 | MagFormFactors = [] |
---|
226 | while S: |
---|
227 | S = FFdata.readline() |
---|
228 | if S[3:5] == ElS: |
---|
229 | if S[8:9] == 'M': |
---|
230 | SN = FFdata.readline() #'N' is assumed to follow 'M' in Atomdata.asc |
---|
231 | Z=int(S[:2]) |
---|
232 | Symbol = S[3:7] |
---|
233 | S = S[12:] |
---|
234 | SN = SN[12:] |
---|
235 | mfa = (float(S[:7]),float(S[14:21]),float(S[28:35]),float(S[42:49])) |
---|
236 | mfb = (float(S[7:14]),float(S[21:28]),float(S[35:42]),float(S[49:56])) |
---|
237 | nfa = (float(SN[:7]),float(SN[14:21]),float(SN[28:35]),float(SN[42:49])) |
---|
238 | nfb = (float(SN[7:14]),float(SN[21:28]),float(SN[35:42]),float(SN[49:56])) |
---|
239 | FormFac = {'Symbol':Symbol,'Z':Z,'mfa':mfa,'nfa':nfa,'mfb':mfb,'nfb':nfb, |
---|
240 | 'mfc':float(S[56:63]),'nfc':float(SN[56:63])} |
---|
241 | MagFormFactors.append(FormFac) |
---|
242 | FFdata.close() |
---|
243 | return MagFormFactors |
---|
244 | |
---|
245 | def ScatFac(El, SQ): |
---|
246 | """compute value of form factor |
---|
247 | @param El: element dictionary defined in GetFormFactorCoeff |
---|
248 | @param SQ: (sin-theta/lambda)**2 |
---|
249 | @return: real part of form factor |
---|
250 | """ |
---|
251 | fa = np.array(El['fa']) |
---|
252 | fb = np.array(El['fb']) |
---|
253 | t = -fb[:,np.newaxis]*SQ |
---|
254 | return np.sum(fa[:,np.newaxis]*np.exp(t)[:],axis=0)+El['fc'] |
---|
255 | |
---|
256 | def BlenRes(BLdata,wave): |
---|
257 | FP = [] |
---|
258 | FPP = [] |
---|
259 | Emev = 81.80703/wave**2 |
---|
260 | for BL in BLdata: |
---|
261 | if len(BL) >= 6: |
---|
262 | Emev = 81.80703/wave**2 |
---|
263 | G2 = BL[5]**2 |
---|
264 | T = [Emev-BL[4],0,0] |
---|
265 | D = [T**2+G2,0,0] |
---|
266 | fp = T/D |
---|
267 | fpp = 1.0/D |
---|
268 | if len(BL) == 8: |
---|
269 | T = Emev-BL[7] |
---|
270 | D = T**2+G2 |
---|
271 | fp += BL[6]*T/D |
---|
272 | fpp += BL[6]/D |
---|
273 | if len(BL) == 10: |
---|
274 | T = Emev-BL[9] |
---|
275 | D = T**2+G2 |
---|
276 | fp += BL[8]*T/D |
---|
277 | fpp += BL[8]/D |
---|
278 | FP.append(BL[2]*fp) |
---|
279 | FPP.append(-BL[3]*fpp) |
---|
280 | else: |
---|
281 | FP.append(0.0) |
---|
282 | FPP.append(0.0) |
---|
283 | return np.array(FP),np.array(FPP) |
---|
284 | |
---|
285 | def ComptonFac(El,SQ): |
---|
286 | """compute Compton scattering factor |
---|
287 | @param El: element dictionary |
---|
288 | @param SQ: (sin-theta/lambda)**2 |
---|
289 | @return: compton scattering factor |
---|
290 | """ |
---|
291 | ca = np.array(El['cmpa']) |
---|
292 | cb = np.array(El['cmpb']) |
---|
293 | t = -cb[:,np.newaxis]*SQ |
---|
294 | return El['cmpz']-np.sum(ca[:,np.newaxis]*np.exp(t),axis=0) |
---|
295 | |
---|
296 | def FPcalc(Orbs, KEv): |
---|
297 | """Compute real & imaginary resonant X-ray scattering factors |
---|
298 | @param Orbs: list of orbital dictionaries as defined in GetXsectionCoeff |
---|
299 | @param KEv: x-ray energy in keV |
---|
300 | @return: C: (f',f",mu): real, imaginary parts of resonant scattering & atomic absorption coeff. |
---|
301 | """ |
---|
302 | def Aitken(Orb, LKev): |
---|
303 | Nval = Orb['Nval'] |
---|
304 | j = Nval-1 |
---|
305 | LEner = Orb['LEner'] |
---|
306 | for i in range(Nval): |
---|
307 | if LEner[i] <= LKev: j = i |
---|
308 | if j > Nval-3: j= Nval-3 |
---|
309 | T = [0,0,0,0,0,0] |
---|
310 | LXSect = Orb['LXSect'] |
---|
311 | for i in range(3): |
---|
312 | T[i] = LXSect[i+j] |
---|
313 | T[i+3] = LEner[i+j]-LKev |
---|
314 | T[1] = (T[0]*T[4]-T[1]*T[3])/(LEner[j+1]-LEner[j]) |
---|
315 | T[2] = (T[0]*T[5]-T[2]*T[3])/(LEner[j+2]-LEner[j]) |
---|
316 | T[2] = (T[1]*T[5]-T[2]*T[4])/(LEner[j+2]-LEner[j+1]) |
---|
317 | C = T[2] |
---|
318 | return C |
---|
319 | |
---|
320 | def DGauss(Orb,CX,RX,ISig): |
---|
321 | ALG = (0.11846344252810,0.23931433524968,0.284444444444, |
---|
322 | 0.23931433524968,0.11846344252810) |
---|
323 | XLG = (0.04691007703067,0.23076534494716,0.5, |
---|
324 | 0.76923465505284,0.95308992296933) |
---|
325 | |
---|
326 | D = 0.0 |
---|
327 | B2 = Orb['BB']**2 |
---|
328 | R2 = RX**2 |
---|
329 | XSecIP = Orb['XSecIP'] |
---|
330 | for i in range(5): |
---|
331 | X = XLG[i] |
---|
332 | X2 = X**2 |
---|
333 | XS = XSecIP[i] |
---|
334 | if ISig == 0: |
---|
335 | S = BB*(XS*(B2/X2)-CX*R2)/(R2*X2-B2) |
---|
336 | elif ISig == 1: |
---|
337 | S = 0.5*BB*B2*XS/(math.sqrt(X)*(R2*X2-X*B2)) |
---|
338 | elif ISig == 2: |
---|
339 | T = X*X2*R2-B2/X |
---|
340 | S = 2.0*BB*(XS*B2/(T*X2**2)-(CX*R2/T)) |
---|
341 | else: |
---|
342 | S = BB*B2*(XS-Orb['SEdge']*X2)/(R2*X2**2-X2*B2) |
---|
343 | A = ALG[i] |
---|
344 | D += A*S |
---|
345 | return D |
---|
346 | |
---|
347 | AU = 2.80022e+7 |
---|
348 | C1 = 0.02721 |
---|
349 | C = 137.0367 |
---|
350 | FP = 0.0 |
---|
351 | FPP = 0.0 |
---|
352 | Mu = 0.0 |
---|
353 | LKev = math.log(KEv) |
---|
354 | RX = KEv/C1 |
---|
355 | if Orbs: |
---|
356 | for Orb in Orbs: |
---|
357 | CX = 0.0 |
---|
358 | BB = Orb['BB'] |
---|
359 | BindEn = Orb['BindEn'] |
---|
360 | if Orb['IfBe'] != 0: ElEterm = Orb['ElEterm'] |
---|
361 | if BindEn <= KEv: |
---|
362 | CX = math.exp(Aitken(Orb,LKev)) |
---|
363 | Mu += CX |
---|
364 | CX /= AU |
---|
365 | Corr = 0.0 |
---|
366 | if Orb['IfBe'] == 0 and BindEn >= KEv: |
---|
367 | CX = 0.0 |
---|
368 | FPI = DGauss(Orb,CX,RX,3) |
---|
369 | Corr = 0.5*Orb['SEdge']*BB**2*math.log((RX-BB)/(-RX-BB))/RX |
---|
370 | else: |
---|
371 | FPI = DGauss(Orb,CX,RX,Orb['IfBe']) |
---|
372 | if CX != 0.0: Corr = -0.5*CX*RX*math.log((RX+BB)/(RX-BB)) |
---|
373 | FPI = (FPI+Corr)*C/(2.0*math.pi**2) |
---|
374 | FPPI = C*CX*RX/(4.0*math.pi) |
---|
375 | FP += FPI |
---|
376 | FPP += FPPI |
---|
377 | FP -= ElEterm |
---|
378 | |
---|
379 | return (FP, FPP, Mu) |
---|
380 | |
---|
381 | |
---|