source: trunk/GSASIIElem.py @ 2546

Last change on this file since 2546 was 2546, checked in by vondreele, 5 years ago

cleanup of all GSASII main routines - remove unused variables, dead code, etc.

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