source: trunk/GSASIIElem.py @ 2563

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

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