source: trunk/GSASIIElem.py @ 1360

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

add SRM660b to ImageCalibrants?.py
revamp element stuff to use atmdata.py instead of atmdata.dat
should be immune to old gpx files.

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