source: trunk/GSASIIElem.py @ 1391

Last change on this file since 1391 was 1391, checked in by vondreele, 9 years ago

remove stray prints
Breit-Wigner & Lynn & Seeger neutron anomalous form factors for CW data
Add 'Ext' to HKLF reflection lists
make sure wave & tbar are in TOF HKLF reflection lists
add '_refln_f_sigma' & '_refln_f_squared_sigma' to possible cif reflection columns
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 13.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: 2014-06-20 18:43:42 +0000 (Fri, 20 Jun 2014) $
10# $Author: vondreele $
11# $Revision: 1391 $
12# $URL: trunk/GSASIIElem.py $
13# $Id: GSASIIElem.py 1391 2014-06-20 18:43:42Z vondreele $
14########### SVN repository information ###################
15
16import math
17import os.path
18import GSASIIpath
19GSASIIpath.SetVersionNumber("$Revision: 1391 $")
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    isotope = General['Isotope']
71    for El in atomTypes:
72        ElS = getElSym(El)
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            if 'BW-LS' in El:
97                BLvals.append(BLtables[El][1]['BW-LS'][0])
98            else:
99                BLvals.append(BLtables[El][1]['SL'][0])
100    else:
101        BLvals = {}
102        for El in BLtables:
103            if 'BW-LS' in El:
104                BLvals[El] = BLtables[El][1]['BW-LS'][0]
105            else:
106                BLvals[El] = BLtables[El][1]['SL'][0]
107    return BLvals
108       
109def GetFFC5(ElSym):
110    '''Get 5 term form factor and Compton scattering data
111
112    :param ElSym: str(1-2 character element symbol with proper case);
113    :return El: dictionary with 5 term form factor & compton coefficients
114    '''
115    import FormFactors as FF
116    El = {}
117    FF5 = FF.FFac5term[ElSym]
118    El['fa'] = FF5[:5]
119    El['fc'] = FF5[5]
120    El['fb'] = FF5[6:]
121    Cmp5 = FF.Compton[ElSym]
122    El['cmpz'] = Cmp5[0]
123    El['cmpa'] = Cmp5[1:6]
124    El['cmpb'] = Cmp5[6:]
125    return El
126   
127def CheckElement(El):
128    '''Check if element El is in the periodic table
129
130    :param str El: One or two letter element symbol, capitaliztion ignored
131    :returns: True if the element is found
132
133    '''
134    import ElementTable as ET
135    Elements = []
136    for elem in ET.ElTable:
137        Elements.append(elem[0][0])
138    if El.capitalize() in Elements:
139        return True
140    else:
141        return False 
142
143def FixValence(El):
144    'Returns the element symbol, even when a valence is present'
145    if '+' in El[-1]: #converts An+/- to A+/-n
146        num = El[-2]
147        El = El.split(num)[0]+'+'+num
148    if '+0' in El:
149        El = El.split('+0')[0]
150    if '-' in El[-1]:
151        num = El[-2]
152        El = El.split(num)[0]+'-'+num
153    if '-0' in El:
154        El = El.split('-0')[0]
155    return El
156   
157def GetAtomInfo(El):
158    'reads element information from atmdata.py'
159    import ElementTable as ET
160    Elements = [elem[0][0] for elem in ET.ElTable]
161    AtomInfo = {}
162    ElS = getElSym(El)
163    AtomInfo.update(dict(zip(['Drad','Arad','Vdrad','Hbrad'],atmdata.AtmSize[ElS])))
164    AtomInfo['Symbol'] = El
165    AtomInfo['Color'] = ET.ElTable[Elements.index(ElS)][6]
166    AtomInfo['Z'] = atmdata.XrayFF[El]['Z']
167    isotopes = [ky for ky in atmdata.AtmBlens.keys() if ElS == ky.split('_')[0]]
168    isotopes.sort()
169    AtomInfo['Mass'] = atmdata.AtmBlens[isotopes[0]]['Mass']    #default to nat. abund.
170    AtomInfo['Isotopes'] = {}
171    for isotope in isotopes:
172        data = atmdata.AtmBlens[isotope]
173        if isotope == ElS+'_':
174            AtomInfo['Isotopes']['Nat. Abund.'] = data
175        else:
176            AtomInfo['Isotopes'][isotope.split('_')[1]] = data
177    return AtomInfo
178       
179def GetXsectionCoeff(El):
180    """Read atom orbital scattering cross sections for fprime calculations via Cromer-Lieberman algorithm
181
182    :param El: 2 character element symbol
183    :return: Orbs: list of orbitals each a dictionary with detailed orbital information used by FPcalc
184
185    each dictionary is:
186
187    * 'OrbName': Orbital name read from file
188    * 'IfBe' 0/2 depending on orbital
189    * 'BindEn': binding energy
190    * 'BB': BindEn/0.02721
191    * 'XSectIP': 5 cross section inflection points
192    * 'ElEterm': energy correction term
193    * 'SEdge': absorption edge for orbital
194    * 'Nval': 10/11 depending on IfBe
195    * 'LEner': 10/11 values of log(energy)
196    * 'LXSect': 10/11 values of log(cross section)
197
198    """
199    AU = 2.80022e+7
200    C1 = 0.02721
201    ElS = El.upper()
202    ElS = ElS.ljust(2)
203    filename = os.path.join(os.path.split(__file__)[0],'Xsect.dat')
204    try:
205        xsec = open(filename,'Ur')
206    except:
207        print '**** ERROR - File Xsect.dat not found in directory %s' % sys.path[0]
208        sys.exit()
209    S = '1'
210    Orbs = []
211    while S:
212        S = xsec.readline()
213        if S[:2] == ElS:
214            S = S[:-1]+xsec.readline()[:-1]+xsec.readline()
215            OrbName = S[9:14]
216            S = S[14:]
217            IfBe = int(S[0])
218            S = S[1:]
219            val = S.split()
220            BindEn = float(val[0])
221            BB = BindEn/C1
222            Orb = {'OrbName':OrbName,'IfBe':IfBe,'BindEn':BindEn,'BB':BB}
223            Energy = []
224            XSect = []
225            for i in range(11):
226                Energy.append(float(val[2*i+1]))
227                XSect.append(float(val[2*i+2]))
228            XSecIP = []
229            for i in range(5): XSecIP.append(XSect[i+5]/AU)
230            Orb['XSecIP'] = XSecIP
231            if IfBe == 0:
232                Orb['SEdge'] = XSect[10]/AU
233                Nval = 11
234            else:
235                Orb['ElEterm'] = XSect[10]
236                del Energy[10]
237                del XSect[10]
238                Nval = 10
239                Orb['SEdge'] = 0.0
240            Orb['Nval'] = Nval
241            D = dict(zip(Energy,XSect))
242            Energy.sort()
243            X = []
244            for key in Energy:
245                X.append(D[key])
246            XSect = X
247            LEner = []
248            LXSect = []
249            for i in range(Nval):
250                LEner.append(math.log(Energy[i]))
251                if XSect[i] > 0.0:
252                    LXSect.append(math.log(XSect[i]))
253                else:
254                    LXSect.append(0.0)
255            Orb['LEner'] = LEner
256            Orb['LXSect'] = LXSect
257            Orbs.append(Orb)
258    xsec.close()
259    return Orbs
260   
261def GetMagFormFacCoeff(El):
262    """Read magnetic form factor data from atomdata.asc file
263
264    :param El: 2 character element symbol
265    :return: MagFormFactors: list of all magnetic form factors dictionaries for element El.
266
267    each dictionary contains:
268
269    * 'Symbol':Symbol
270    * 'Z':Z
271    * 'mfa': 4 MA coefficients
272    * 'nfa': 4 NA coefficients
273    * 'mfb': 4 MB coefficients
274    * 'nfb': 4 NB coefficients
275    * 'mfc': MC coefficient
276    * 'nfc': NC coefficient
277   
278    """
279    MagFormFactors = []
280    mags = [ky for ky in atmdata.MagFF.keys() if El == getElSym(ky)]
281    for mag in mags:
282        magData = {}
283        data = atmdata.MagFF[mag]
284        magData['Symbol'] = mag
285        magData['Z'] = atmdata.XrayFF[getElSym(mags)]['Z']
286        magData['mfa'] = [data['M'][i] for i in [0,2,4,6]]
287        magdata['mfb'] = [data['M'][i] for i in [1,3,5,7]]
288        magdata['mfc'] = data['M'][8]
289        magData['nfa'] = [data['N'][i] for i in [0,2,4,6]]
290        magdata['nfb'] = [data['N'][i] for i in [1,3,5,7]]
291        magdata['nfc'] = data['N'][8]
292        magdata['g-fac'] = data['N'][9]
293        MagFormFactors.append(magdata)
294    return MagFormFactors
295
296def ScatFac(El, SQ):
297    """compute value of form factor
298
299    :param El: element dictionary defined in GetFormFactorCoeff
300    :param SQ: (sin-theta/lambda)**2
301    :return: real part of form factor
302    """
303    fa = np.array(El['fa'])
304    fb = np.array(El['fb'])
305    t = -fb[:,np.newaxis]*SQ
306    return np.sum(fa[:,np.newaxis]*np.exp(t)[:],axis=0)+El['fc']
307       
308def BlenResCW(Els,BLtables,wave):
309    FP = np.zeros(len(Els))
310    FPP = np.zeros(len(Els))
311    for i,El in enumerate(Els):
312        BL = BLtables[El][1]
313        if 'BW-LS' in BL:
314            Re,Im,E0,gam,A,E1,B,E2 = BL['BW-LS'][1:]
315            Emev = 81.80703/wave**2
316            T0 = Emev-E0
317            T1 = Emev-E1
318            T2 = Emev-E2
319            D0 = T0**2+gam**2
320            D1 = T1**2+gam**2
321            D2 = T2**2+gam**2
322            FP[i] = Re*(T0/D0+A*T1/D1+B*T2/D2)
323            FPP[i] = Im*(1/D0+A/D1+B/D2)
324        else:
325            FPP[i] = BL['SL'][1]    #for Li, B, etc.
326    return FP,FPP
327   
328def BlenResTOF(El,BLtables,wave):
329#    FP = np.zeros(len(wave))
330    FPP = np.zeros(len(wave))
331    BL = BLtables[El][1]
332    print BL
333    if 'BW-LS' in BL:
334        Re,Im,E0,gam,A,E1,B,E2 = BL['BW-LS'][1:]
335        Emev = 81.80703/wave**2
336        T0 = Emev-E0
337        T1 = Emev-E1
338        T2 = Emev-E2
339        D0 = T0**2+gam**2
340        D1 = T1**2+gam**2
341        D2 = T2**2+gam**2
342        FP = Re*(T0/D0+A*T1/D1+B*T2/D2)
343        FPP = Im*(1/D0+A/D1+B/D2)
344    else:
345        FPP = np.ones(len(wave))*BL['SL'][1]    #for Li, B, etc.
346    return FP,FPP
347   
348def ComptonFac(El,SQ):
349    """compute Compton scattering factor
350
351    :param El: element dictionary
352    :param SQ: (sin-theta/lambda)**2
353    :return: compton scattering factor
354    """   
355    ca = np.array(El['cmpa'])
356    cb = np.array(El['cmpb'])
357    t = -cb[:,np.newaxis]*SQ       
358    return El['cmpz']-np.sum(ca[:,np.newaxis]*np.exp(t),axis=0) 
359           
360def FPcalc(Orbs, KEv):
361    """Compute real & imaginary resonant X-ray scattering factors
362
363    :param Orbs: list of orbital dictionaries as defined in GetXsectionCoeff
364    :param KEv: x-ray energy in keV
365    :return: C: (f',f",mu): real, imaginary parts of resonant scattering & atomic absorption coeff.
366    """
367    def Aitken(Orb, LKev):
368        Nval = Orb['Nval']
369        j = Nval-1
370        LEner = Orb['LEner']
371        for i in range(Nval):
372            if LEner[i] <= LKev: j = i
373        if j > Nval-3: j= Nval-3
374        T = [0,0,0,0,0,0]
375        LXSect = Orb['LXSect']
376        for i in range(3):
377           T[i] = LXSect[i+j]
378           T[i+3] = LEner[i+j]-LKev
379        T[1] = (T[0]*T[4]-T[1]*T[3])/(LEner[j+1]-LEner[j])
380        T[2] = (T[0]*T[5]-T[2]*T[3])/(LEner[j+2]-LEner[j])
381        T[2] = (T[1]*T[5]-T[2]*T[4])/(LEner[j+2]-LEner[j+1])
382        C = T[2]
383        return C
384   
385    def DGauss(Orb,CX,RX,ISig):
386        ALG = (0.11846344252810,0.23931433524968,0.284444444444,
387        0.23931433524968,0.11846344252810)
388        XLG = (0.04691007703067,0.23076534494716,0.5,
389        0.76923465505284,0.95308992296933)
390       
391        D = 0.0
392        B2 = Orb['BB']**2
393        R2 = RX**2
394        XSecIP = Orb['XSecIP']
395        for i in range(5):
396            X = XLG[i]
397            X2 = X**2
398            XS = XSecIP[i]
399            if ISig == 0:
400                S = BB*(XS*(B2/X2)-CX*R2)/(R2*X2-B2)
401            elif ISig == 1:
402                S = 0.5*BB*B2*XS/(math.sqrt(X)*(R2*X2-X*B2))
403            elif ISig == 2:
404                T = X*X2*R2-B2/X
405                S = 2.0*BB*(XS*B2/(T*X2**2)-(CX*R2/T))
406            else:
407                S = BB*B2*(XS-Orb['SEdge']*X2)/(R2*X2**2-X2*B2)
408            A = ALG[i]
409            D += A*S
410        return D
411   
412    AU = 2.80022e+7
413    C1 = 0.02721
414    C = 137.0367
415    FP = 0.0
416    FPP = 0.0
417    Mu = 0.0
418    LKev = math.log(KEv)
419    RX = KEv/C1
420    if Orbs:
421        for Orb in Orbs:
422            CX = 0.0
423            BB = Orb['BB']
424            BindEn = Orb['BindEn']
425            if Orb['IfBe'] != 0: ElEterm = Orb['ElEterm']
426            if BindEn <= KEv:
427                CX = math.exp(Aitken(Orb,LKev))
428                Mu += CX
429                CX /= AU
430            Corr = 0.0
431            if Orb['IfBe'] == 0 and BindEn >= KEv:
432                CX = 0.0
433                FPI = DGauss(Orb,CX,RX,3)
434                Corr = 0.5*Orb['SEdge']*BB**2*math.log((RX-BB)/(-RX-BB))/RX
435            else:
436                FPI = DGauss(Orb,CX,RX,Orb['IfBe'])
437                if CX != 0.0: Corr = -0.5*CX*RX*math.log((RX+BB)/(RX-BB))
438            FPI = (FPI+Corr)*C/(2.0*math.pi**2)
439            FPPI = C*CX*RX/(4.0*math.pi)
440            FP += FPI
441            FPP += FPPI
442        FP -= ElEterm
443   
444    return (FP, FPP, Mu)
445   
446
Note: See TracBrowser for help on using the repository browser.