source: trunk/GSASIIElem.py @ 4605

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

more RMCProfile setup stuff added

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