source: trunk/GSASIIElem.py @ 1364

Last change on this file since 1364 was 1364, checked in by vondreele, 7 years ago

read TOPAZ structure factor files

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 13.8 KB
Line 
1# -*- coding: utf-8 -*-
2"""
3*GSASIIElem: functions for element types*
4-----------------------------------------
5
6"""
7# Copyright: 2008, Robert B. Von Dreele & Brian H. Toby (Argonne National Laboratory)
8########### SVN repository information ###################
9# $Date: 2014-05-27 17:14:59 +0000 (Tue, 27 May 2014) $
10# $Author: vondreele $
11# $Revision: 1364 $
12# $URL: trunk/GSASIIElem.py $
13# $Id: GSASIIElem.py 1364 2014-05-27 17:14:59Z vondreele $
14########### SVN repository information ###################
15
16import math
17import os.path
18import GSASIIpath
19GSASIIpath.SetVersionNumber("$Revision: 1364 $")
20import numpy as np
21import atmdata
22
23getElSym = lambda sym: sym.split('+')[0].split('-')[0]
24def GetFormFactorCoeff(El):
25    """Read X-ray form factor coefficients from `atomdata.py` file
26
27    :param str El: element 1-2 character symbol, case irrevelant
28    :return: `FormFactors`: list of form factor dictionaries
29   
30    Each X-ray form factor dictionary is:
31   
32    * `Symbol`: 4 character element symbol with valence (e.g. 'NI+2')
33    * `Z`: atomic number
34    * `fa`: 4 A coefficients
35    * `fb`: 4 B coefficients
36    * `fc`: C coefficient
37   
38    """
39   
40    Els = El.capitalize().strip()
41    valences = [ky for ky in atmdata.XrayFF.keys() if Els == getElSym(ky)]
42    FormFactors = [atmdata.XrayFF[val] for val in valences]
43    for Sy,FF in zip(valences,FormFactors):
44        FF.update({'Symbol':Sy.upper()})
45    return FormFactors
46   
47def GetFFtable(atomTypes):
48    ''' returns a dictionary of form factor data for atom types found in atomTypes
49
50    :param list atomTypes: list of atom types
51    :return: FFtable, dictionary of form factor data; key is atom type
52
53    '''
54    FFtable = {}
55    for El in atomTypes:
56        FFs = GetFormFactorCoeff(getElSym(El))
57        for item in FFs:
58            if item['Symbol'] == El.upper():
59                FFtable[El] = item
60    return FFtable
61   
62def GetBLtable(General):
63    ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
64
65    :param dict General: dictionary of phase info.; includes AtomTypes & Isotopes
66    :return: BLtable, dictionary of scattering length data; key is atom type
67    '''
68    atomTypes = General['AtomTypes']
69    BLtable = {}
70    isotopes = General['Isotopes']
71    isotope = General['Isotope']
72    for El in atomTypes:
73        ElS = getElSym(El)
74        if 'Nat' in isotope[El]:
75            BLtable[El] = [isotope[El],atmdata.AtmBlens[ElS+'_']]
76        else:
77            BLtable[El] = [isotope[El],atmdata.AtmBlens[ElS+'_'+isotope[El]]]
78    return BLtable
79       
80def getFFvalues(FFtables,SQ,ifList=False):
81    'Needs a doc string'
82    if ifList:
83        FFvals = []
84        for El in FFtables:
85            FFvals.append(ScatFac(FFtables[El],SQ)[0])
86    else:
87        FFvals = {}
88        for El in FFtables:
89            FFvals[El] = ScatFac(FFtables[El],SQ)[0]
90    return FFvals
91   
92def getBLvalues(BLtables,ifList=False):
93    'Needs a doc string'
94    if ifList:
95        BLvals = []
96        for El in BLtables:
97            BLvals.append(BLtables[El][1]['SL'][0])
98    else:
99        BLvals = {}
100        for El in BLtables:
101            BLvals[El] = BLtables[El][1]['SL'][0]
102    return BLvals
103       
104def GetFFC5(ElSym):
105    '''Get 5 term form factor and Compton scattering data
106
107    :param ElSym: str(1-2 character element symbol with proper case);
108    :return El: dictionary with 5 term form factor & compton coefficients
109    '''
110    import FormFactors as FF
111    El = {}
112    FF5 = FF.FFac5term[ElSym]
113    El['fa'] = FF5[:5]
114    El['fc'] = FF5[5]
115    El['fb'] = FF5[6:]
116    Cmp5 = FF.Compton[ElSym]
117    El['cmpz'] = Cmp5[0]
118    El['cmpa'] = Cmp5[1:6]
119    El['cmpb'] = Cmp5[6:]
120    return El
121   
122def CheckElement(El):
123    '''Check if element El is in the periodic table
124
125    :param str El: One or two letter element symbol, capitaliztion ignored
126    :returns: True if the element is found
127
128    '''
129    import ElementTable as ET
130    Elements = []
131    for elem in ET.ElTable:
132        Elements.append(elem[0][0])
133    if El.capitalize() in Elements:
134        return True
135    else:
136        return False 
137
138def FixValence(El):
139    'Returns the element symbol, even when a valence is present'
140    if '+' in El[-1]: #converts An+/- to A+/-n
141        num = El[-2]
142        El = El.split(num)[0]+'+'+num
143    if '+0' in El:
144        El = El.split('+0')[0]
145    if '-' in El[-1]:
146        num = El[-2]
147        El = El.split(num)[0]+'-'+num
148    if '-0' in El:
149        El = El.split('-0')[0]
150    return El
151   
152def GetAtomInfo(El):
153    'reads element information from atmdata.py'
154    import ElementTable as ET
155    Elements = [elem[0][0] for elem in ET.ElTable]
156    AtomInfo = {}
157    ElS = getElSym(El)
158    AtomInfo.update(dict(zip(['Drad','Arad','Vdrad','Hbrad'],atmdata.AtmSize[ElS])))
159    AtomInfo['Symbol'] = El
160    AtomInfo['Color'] = ET.ElTable[Elements.index(ElS)][6]
161    AtomInfo['Z'] = atmdata.XrayFF[El]['Z']
162    isotopes = [ky for ky in atmdata.AtmBlens.keys() if ElS == ky.split('_')[0]]
163    isotopes.sort()
164    AtomInfo['Mass'] = atmdata.AtmBlens[isotopes[0]]['Mass']    #default to nat. abund.
165    AtomInfo['Isotopes'] = {}
166    for isotope in isotopes:
167        data = atmdata.AtmBlens[isotope]
168        if isotope == ElS+'_':
169            if data['SL']:
170                AtomInfo['Isotopes']['Nat. Abund.'] = data
171        else:
172            AtomInfo['Isotopes'][isotope.split('_')[1]] = data
173    return AtomInfo
174       
175def GetXsectionCoeff(El):
176    """Read atom orbital scattering cross sections for fprime calculations via Cromer-Lieberman algorithm
177
178    :param El: 2 character element symbol
179    :return: Orbs: list of orbitals each a dictionary with detailed orbital information used by FPcalc
180
181    each dictionary is:
182
183    * 'OrbName': Orbital name read from file
184    * 'IfBe' 0/2 depending on orbital
185    * 'BindEn': binding energy
186    * 'BB': BindEn/0.02721
187    * 'XSectIP': 5 cross section inflection points
188    * 'ElEterm': energy correction term
189    * 'SEdge': absorption edge for orbital
190    * 'Nval': 10/11 depending on IfBe
191    * 'LEner': 10/11 values of log(energy)
192    * 'LXSect': 10/11 values of log(cross section)
193
194    """
195    AU = 2.80022e+7
196    C1 = 0.02721
197    ElS = El.upper()
198    ElS = ElS.ljust(2)
199    filename = os.path.join(os.path.split(__file__)[0],'Xsect.dat')
200    try:
201        xsec = open(filename,'Ur')
202    except:
203        print '**** ERROR - File Xsect.dat not found in directory %s' % sys.path[0]
204        sys.exit()
205    S = '1'
206    Orbs = []
207    while S:
208        S = xsec.readline()
209        if S[:2] == ElS:
210            S = S[:-1]+xsec.readline()[:-1]+xsec.readline()
211            OrbName = S[9:14]
212            S = S[14:]
213            IfBe = int(S[0])
214            S = S[1:]
215            val = S.split()
216            BindEn = float(val[0])
217            BB = BindEn/C1
218            Orb = {'OrbName':OrbName,'IfBe':IfBe,'BindEn':BindEn,'BB':BB}
219            Energy = []
220            XSect = []
221            for i in range(11):
222                Energy.append(float(val[2*i+1]))
223                XSect.append(float(val[2*i+2]))
224            XSecIP = []
225            for i in range(5): XSecIP.append(XSect[i+5]/AU)
226            Orb['XSecIP'] = XSecIP
227            if IfBe == 0:
228                Orb['SEdge'] = XSect[10]/AU
229                Nval = 11
230            else:
231                Orb['ElEterm'] = XSect[10]
232                del Energy[10]
233                del XSect[10]
234                Nval = 10
235                Orb['SEdge'] = 0.0
236            Orb['Nval'] = Nval
237            D = dict(zip(Energy,XSect))
238            Energy.sort()
239            X = []
240            for key in Energy:
241                X.append(D[key])
242            XSect = X
243            LEner = []
244            LXSect = []
245            for i in range(Nval):
246                LEner.append(math.log(Energy[i]))
247                if XSect[i] > 0.0:
248                    LXSect.append(math.log(XSect[i]))
249                else:
250                    LXSect.append(0.0)
251            Orb['LEner'] = LEner
252            Orb['LXSect'] = LXSect
253            Orbs.append(Orb)
254    xsec.close()
255    return Orbs
256   
257def GetMagFormFacCoeff(El):
258    """Read magnetic form factor data from atomdata.asc file
259
260    :param El: 2 character element symbol
261    :return: MagFormFactors: list of all magnetic form factors dictionaries for element El.
262
263    each dictionary contains:
264
265    * 'Symbol':Symbol
266    * 'Z':Z
267    * 'mfa': 4 MA coefficients
268    * 'nfa': 4 NA coefficients
269    * 'mfb': 4 MB coefficients
270    * 'nfb': 4 NB coefficients
271    * 'mfc': MC coefficient
272    * 'nfc': NC coefficient
273   
274    """
275    MagFormFactors = []
276    mags = [ky for ky in atmdata.MagFF.keys() if El == getElSym(ky)]
277    for mag in mags:
278        magData = {}
279        data = atmdata.MagFF[mag]
280        magData['Symbol'] = mag
281        magData['Z'] = atmdata.XrayFF[getElSym(mags)]['Z']
282        magData['mfa'] = [data['M'][i] for i in [0,2,4,6]]
283        magdata['mfb'] = [data['M'][i] for i in [1,3,5,7]]
284        magdata['mfc'] = data['M'][8]
285        magData['nfa'] = [data['N'][i] for i in [0,2,4,6]]
286        magdata['nfb'] = [data['N'][i] for i in [1,3,5,7]]
287        magdata['nfc'] = data['N'][8]
288        magdata['g-fac'] = data['N'][9]
289        MagFormFactors.append(magdata)
290    return MagFormFactors
291
292def ScatFac(El, SQ):
293    """compute value of form factor
294
295    :param El: element dictionary defined in GetFormFactorCoeff
296    :param SQ: (sin-theta/lambda)**2
297    :return: real part of form factor
298    """
299    fa = np.array(El['fa'])
300    fb = np.array(El['fb'])
301    t = -fb[:,np.newaxis]*SQ
302    return np.sum(fa[:,np.newaxis]*np.exp(t)[:],axis=0)+El['fc']
303       
304def BlenRes(Elist,BLtables,wave):
305    FP = np.zeros(len(Elist))
306    FPP = np.zeros(len(Elist))
307    Emev = 81.80703/wave**2
308    for i,El in enumerate(Elist):
309        BL = BLtables[El]
310        if len(BL) >= 6:
311            Emev = 81.80703/wave**2
312            G2 = BL[5]**2
313            T = [Emev-BL[4],0,0]
314            D = [T**2+G2,0,0]
315            fp = T/D
316            fpp = 1.0/D
317            if len(BL) == 8:
318                T = Emev-BL[7]
319                D = T**2+G2
320                fp += BL[6]*T/D
321                fpp += BL[6]/D
322            if len(BL) == 10:
323                T = Emev-BL[9]
324                D = T**2+G2
325                fp += BL[8]*T/D
326                fpp += BL[8]/D
327            FP[i] = (BL[2]*fp)
328            FPP[i] = (-BL[3]*fpp)
329        else:
330            FP[i] = 0.0
331            FPP[i] = 0.0
332    return FP,FPP
333   
334#def BlenRes(BLdata,wave):
335#    FP = []
336#    FPP = []
337#    Emev = 81.80703/wave**2
338#    for BL in BLdata:
339#        if len(BL) >= 6:
340#            G2 = BL[5]**2
341#            T = [Emev-BL[4],0,0]
342#            D = [T**2+G2,0,0]
343#            fp = T/D
344#            fpp = 1.0/D
345#            if len(BL) == 8:
346#                T = Emev-BL[7]
347#                D = T**2+G2
348#                fp += BL[6]*T/D
349#                fpp += BL[6]/D
350#            if len(BL) == 10:
351#                T = Emev-BL[9]
352#                D = T**2+G2
353#                fp += BL[8]*T/D
354#                fpp += BL[8]/D
355#            FP.append(BL[2]*fp)
356#            FPP.append(-BL[3]*fpp)
357#        else:
358#            FP.append(0.0)
359#            FPP.append(0.0)
360#    return np.array(FP),np.array(FPP)
361   
362def ComptonFac(El,SQ):
363    """compute Compton scattering factor
364
365    :param El: element dictionary
366    :param SQ: (sin-theta/lambda)**2
367    :return: compton scattering factor
368    """   
369    ca = np.array(El['cmpa'])
370    cb = np.array(El['cmpb'])
371    t = -cb[:,np.newaxis]*SQ       
372    return El['cmpz']-np.sum(ca[:,np.newaxis]*np.exp(t),axis=0) 
373           
374def FPcalc(Orbs, KEv):
375    """Compute real & imaginary resonant X-ray scattering factors
376
377    :param Orbs: list of orbital dictionaries as defined in GetXsectionCoeff
378    :param KEv: x-ray energy in keV
379    :return: C: (f',f",mu): real, imaginary parts of resonant scattering & atomic absorption coeff.
380    """
381    def Aitken(Orb, LKev):
382        Nval = Orb['Nval']
383        j = Nval-1
384        LEner = Orb['LEner']
385        for i in range(Nval):
386            if LEner[i] <= LKev: j = i
387        if j > Nval-3: j= Nval-3
388        T = [0,0,0,0,0,0]
389        LXSect = Orb['LXSect']
390        for i in range(3):
391           T[i] = LXSect[i+j]
392           T[i+3] = LEner[i+j]-LKev
393        T[1] = (T[0]*T[4]-T[1]*T[3])/(LEner[j+1]-LEner[j])
394        T[2] = (T[0]*T[5]-T[2]*T[3])/(LEner[j+2]-LEner[j])
395        T[2] = (T[1]*T[5]-T[2]*T[4])/(LEner[j+2]-LEner[j+1])
396        C = T[2]
397        return C
398   
399    def DGauss(Orb,CX,RX,ISig):
400        ALG = (0.11846344252810,0.23931433524968,0.284444444444,
401        0.23931433524968,0.11846344252810)
402        XLG = (0.04691007703067,0.23076534494716,0.5,
403        0.76923465505284,0.95308992296933)
404       
405        D = 0.0
406        B2 = Orb['BB']**2
407        R2 = RX**2
408        XSecIP = Orb['XSecIP']
409        for i in range(5):
410            X = XLG[i]
411            X2 = X**2
412            XS = XSecIP[i]
413            if ISig == 0:
414                S = BB*(XS*(B2/X2)-CX*R2)/(R2*X2-B2)
415            elif ISig == 1:
416                S = 0.5*BB*B2*XS/(math.sqrt(X)*(R2*X2-X*B2))
417            elif ISig == 2:
418                T = X*X2*R2-B2/X
419                S = 2.0*BB*(XS*B2/(T*X2**2)-(CX*R2/T))
420            else:
421                S = BB*B2*(XS-Orb['SEdge']*X2)/(R2*X2**2-X2*B2)
422            A = ALG[i]
423            D += A*S
424        return D
425   
426    AU = 2.80022e+7
427    C1 = 0.02721
428    C = 137.0367
429    FP = 0.0
430    FPP = 0.0
431    Mu = 0.0
432    LKev = math.log(KEv)
433    RX = KEv/C1
434    if Orbs:
435        for Orb in Orbs:
436            CX = 0.0
437            BB = Orb['BB']
438            BindEn = Orb['BindEn']
439            if Orb['IfBe'] != 0: ElEterm = Orb['ElEterm']
440            if BindEn <= KEv:
441                CX = math.exp(Aitken(Orb,LKev))
442                Mu += CX
443                CX /= AU
444            Corr = 0.0
445            if Orb['IfBe'] == 0 and BindEn >= KEv:
446                CX = 0.0
447                FPI = DGauss(Orb,CX,RX,3)
448                Corr = 0.5*Orb['SEdge']*BB**2*math.log((RX-BB)/(-RX-BB))/RX
449            else:
450                FPI = DGauss(Orb,CX,RX,Orb['IfBe'])
451                if CX != 0.0: Corr = -0.5*CX*RX*math.log((RX+BB)/(RX-BB))
452            FPI = (FPI+Corr)*C/(2.0*math.pi**2)
453            FPPI = C*CX*RX/(4.0*math.pi)
454            FP += FPI
455            FPP += FPPI
456        FP -= ElEterm
457   
458    return (FP, FPP, Mu)
459   
460
Note: See TracBrowser for help on using the repository browser.