source: trunk/GSASIIElem.py @ 1365

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

fix atmdata.py for Rh (& other anomalous b isotopes)
fix powder imports to include old style PRCF records in GSAS iparm files
format of penalty fxn. prints
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 13.8 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-28 16:57:03 +0000 (Wed, 28 May 2014) $
10# $Author: vondreele $
11# $Revision: 1365 $
12# $URL: trunk/GSASIIElem.py $
13# $Id: GSASIIElem.py 1365 2014-05-28 16:57:03Z vondreele $
14########### SVN repository information ###################
15
16import math
17import os.path
18import GSASIIpath
19GSASIIpath.SetVersionNumber("$Revision: 1365 $")
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    isotopes = General['Isotopes']
71    isotope = General['Isotope']
72    for El in atomTypes:
73        ElS = getElSym(El)
74        if 'Nat' in isotope[El]:
75            BLtable[El] = [isotope[El],atmdata.AtmBlens[ElS+'_']]
76        else:
77            BLtable[El] = [isotope[El],atmdata.AtmBlens[ElS+'_'+isotope[El]]]
78    return BLtable
79       
80def getFFvalues(FFtables,SQ,ifList=False):
81    'Needs a doc string'
82    if ifList:
83        FFvals = []
84        for El in FFtables:
85            FFvals.append(ScatFac(FFtables[El],SQ)[0])
86    else:
87        FFvals = {}
88        for El in FFtables:
89            FFvals[El] = ScatFac(FFtables[El],SQ)[0]
90    return FFvals
91   
92def getBLvalues(BLtables,ifList=False):
93    'Needs a doc string'
94    if ifList:
95        BLvals = []
96        for El in BLtables:
97            BLvals.append(BLtables[El][1]['SL'][0])
98    else:
99        BLvals = {}
100        for El in BLtables:
101            BLvals[El] = BLtables[El][1]['SL'][0]
102    return BLvals
103       
104def GetFFC5(ElSym):
105    '''Get 5 term form factor and Compton scattering data
106
107    :param ElSym: str(1-2 character element symbol with proper case);
108    :return El: dictionary with 5 term form factor & compton coefficients
109    '''
110    import FormFactors as FF
111    El = {}
112    FF5 = FF.FFac5term[ElSym]
113    El['fa'] = FF5[:5]
114    El['fc'] = FF5[5]
115    El['fb'] = FF5[6:]
116    Cmp5 = FF.Compton[ElSym]
117    El['cmpz'] = Cmp5[0]
118    El['cmpa'] = Cmp5[1:6]
119    El['cmpb'] = Cmp5[6:]
120    return El
121   
122def CheckElement(El):
123    '''Check if element El is in the periodic table
124
125    :param str El: One or two letter element symbol, capitaliztion ignored
126    :returns: True if the element is found
127
128    '''
129    import ElementTable as ET
130    Elements = []
131    for elem in ET.ElTable:
132        Elements.append(elem[0][0])
133    if El.capitalize() in Elements:
134        return True
135    else:
136        return False 
137
138def FixValence(El):
139    'Returns the element symbol, even when a valence is present'
140    if '+' in El[-1]: #converts An+/- to A+/-n
141        num = El[-2]
142        El = El.split(num)[0]+'+'+num
143    if '+0' in El:
144        El = El.split('+0')[0]
145    if '-' in El[-1]:
146        num = El[-2]
147        El = El.split(num)[0]+'-'+num
148    if '-0' in El:
149        El = El.split('-0')[0]
150    return El
151   
152def GetAtomInfo(El):
153    'reads element information from atmdata.py'
154    import ElementTable as ET
155    Elements = [elem[0][0] for elem in ET.ElTable]
156    AtomInfo = {}
157    ElS = getElSym(El)
158    AtomInfo.update(dict(zip(['Drad','Arad','Vdrad','Hbrad'],atmdata.AtmSize[ElS])))
159    AtomInfo['Symbol'] = El
160    AtomInfo['Color'] = ET.ElTable[Elements.index(ElS)][6]
161    AtomInfo['Z'] = atmdata.XrayFF[El]['Z']
162    isotopes = [ky for ky in atmdata.AtmBlens.keys() if ElS == ky.split('_')[0]]
163    isotopes.sort()
164    AtomInfo['Mass'] = atmdata.AtmBlens[isotopes[0]]['Mass']    #default to nat. abund.
165    AtomInfo['Isotopes'] = {}
166    for isotope in isotopes:
167        data = atmdata.AtmBlens[isotope]
168        if isotope == ElS+'_':
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    MagFormFactors = []
275    mags = [ky for ky in atmdata.MagFF.keys() if El == getElSym(ky)]
276    for mag in mags:
277        magData = {}
278        data = atmdata.MagFF[mag]
279        magData['Symbol'] = mag
280        magData['Z'] = atmdata.XrayFF[getElSym(mags)]['Z']
281        magData['mfa'] = [data['M'][i] for i in [0,2,4,6]]
282        magdata['mfb'] = [data['M'][i] for i in [1,3,5,7]]
283        magdata['mfc'] = data['M'][8]
284        magData['nfa'] = [data['N'][i] for i in [0,2,4,6]]
285        magdata['nfb'] = [data['N'][i] for i in [1,3,5,7]]
286        magdata['nfc'] = data['N'][8]
287        magdata['g-fac'] = data['N'][9]
288        MagFormFactors.append(magdata)
289    return MagFormFactors
290
291def ScatFac(El, SQ):
292    """compute value of form factor
293
294    :param El: element dictionary defined in GetFormFactorCoeff
295    :param SQ: (sin-theta/lambda)**2
296    :return: real part of form factor
297    """
298    fa = np.array(El['fa'])
299    fb = np.array(El['fb'])
300    t = -fb[:,np.newaxis]*SQ
301    return np.sum(fa[:,np.newaxis]*np.exp(t)[:],axis=0)+El['fc']
302       
303def BlenRes(Elist,BLtables,wave):
304    FP = np.zeros(len(Elist))
305    FPP = np.zeros(len(Elist))
306    Emev = 81.80703/wave**2
307    for i,El in enumerate(Elist):
308        BL = BLtables[El]
309        if len(BL) >= 6:
310            Emev = 81.80703/wave**2
311            G2 = BL[5]**2
312            T = [Emev-BL[4],0,0]
313            D = [T**2+G2,0,0]
314            fp = T/D
315            fpp = 1.0/D
316            if len(BL) == 8:
317                T = Emev-BL[7]
318                D = T**2+G2
319                fp += BL[6]*T/D
320                fpp += BL[6]/D
321            if len(BL) == 10:
322                T = Emev-BL[9]
323                D = T**2+G2
324                fp += BL[8]*T/D
325                fpp += BL[8]/D
326            FP[i] = (BL[2]*fp)
327            FPP[i] = (-BL[3]*fpp)
328        else:
329            FP[i] = 0.0
330            FPP[i] = 0.0
331    return FP,FPP
332   
333#def BlenRes(BLdata,wave):
334#    FP = []
335#    FPP = []
336#    Emev = 81.80703/wave**2
337#    for BL in BLdata:
338#        if len(BL) >= 6:
339#            G2 = BL[5]**2
340#            T = [Emev-BL[4],0,0]
341#            D = [T**2+G2,0,0]
342#            fp = T/D
343#            fpp = 1.0/D
344#            if len(BL) == 8:
345#                T = Emev-BL[7]
346#                D = T**2+G2
347#                fp += BL[6]*T/D
348#                fpp += BL[6]/D
349#            if len(BL) == 10:
350#                T = Emev-BL[9]
351#                D = T**2+G2
352#                fp += BL[8]*T/D
353#                fpp += BL[8]/D
354#            FP.append(BL[2]*fp)
355#            FPP.append(-BL[3]*fpp)
356#        else:
357#            FP.append(0.0)
358#            FPP.append(0.0)
359#    return np.array(FP),np.array(FPP)
360   
361def ComptonFac(El,SQ):
362    """compute Compton scattering factor
363
364    :param El: element dictionary
365    :param SQ: (sin-theta/lambda)**2
366    :return: compton scattering factor
367    """   
368    ca = np.array(El['cmpa'])
369    cb = np.array(El['cmpb'])
370    t = -cb[:,np.newaxis]*SQ       
371    return El['cmpz']-np.sum(ca[:,np.newaxis]*np.exp(t),axis=0) 
372           
373def FPcalc(Orbs, KEv):
374    """Compute real & imaginary resonant X-ray scattering factors
375
376    :param Orbs: list of orbital dictionaries as defined in GetXsectionCoeff
377    :param KEv: x-ray energy in keV
378    :return: C: (f',f",mu): real, imaginary parts of resonant scattering & atomic absorption coeff.
379    """
380    def Aitken(Orb, LKev):
381        Nval = Orb['Nval']
382        j = Nval-1
383        LEner = Orb['LEner']
384        for i in range(Nval):
385            if LEner[i] <= LKev: j = i
386        if j > Nval-3: j= Nval-3
387        T = [0,0,0,0,0,0]
388        LXSect = Orb['LXSect']
389        for i in range(3):
390           T[i] = LXSect[i+j]
391           T[i+3] = LEner[i+j]-LKev
392        T[1] = (T[0]*T[4]-T[1]*T[3])/(LEner[j+1]-LEner[j])
393        T[2] = (T[0]*T[5]-T[2]*T[3])/(LEner[j+2]-LEner[j])
394        T[2] = (T[1]*T[5]-T[2]*T[4])/(LEner[j+2]-LEner[j+1])
395        C = T[2]
396        return C
397   
398    def DGauss(Orb,CX,RX,ISig):
399        ALG = (0.11846344252810,0.23931433524968,0.284444444444,
400        0.23931433524968,0.11846344252810)
401        XLG = (0.04691007703067,0.23076534494716,0.5,
402        0.76923465505284,0.95308992296933)
403       
404        D = 0.0
405        B2 = Orb['BB']**2
406        R2 = RX**2
407        XSecIP = Orb['XSecIP']
408        for i in range(5):
409            X = XLG[i]
410            X2 = X**2
411            XS = XSecIP[i]
412            if ISig == 0:
413                S = BB*(XS*(B2/X2)-CX*R2)/(R2*X2-B2)
414            elif ISig == 1:
415                S = 0.5*BB*B2*XS/(math.sqrt(X)*(R2*X2-X*B2))
416            elif ISig == 2:
417                T = X*X2*R2-B2/X
418                S = 2.0*BB*(XS*B2/(T*X2**2)-(CX*R2/T))
419            else:
420                S = BB*B2*(XS-Orb['SEdge']*X2)/(R2*X2**2-X2*B2)
421            A = ALG[i]
422            D += A*S
423        return D
424   
425    AU = 2.80022e+7
426    C1 = 0.02721
427    C = 137.0367
428    FP = 0.0
429    FPP = 0.0
430    Mu = 0.0
431    LKev = math.log(KEv)
432    RX = KEv/C1
433    if Orbs:
434        for Orb in Orbs:
435            CX = 0.0
436            BB = Orb['BB']
437            BindEn = Orb['BindEn']
438            if Orb['IfBe'] != 0: ElEterm = Orb['ElEterm']
439            if BindEn <= KEv:
440                CX = math.exp(Aitken(Orb,LKev))
441                Mu += CX
442                CX /= AU
443            Corr = 0.0
444            if Orb['IfBe'] == 0 and BindEn >= KEv:
445                CX = 0.0
446                FPI = DGauss(Orb,CX,RX,3)
447                Corr = 0.5*Orb['SEdge']*BB**2*math.log((RX-BB)/(-RX-BB))/RX
448            else:
449                FPI = DGauss(Orb,CX,RX,Orb['IfBe'])
450                if CX != 0.0: Corr = -0.5*CX*RX*math.log((RX+BB)/(RX-BB))
451            FPI = (FPI+Corr)*C/(2.0*math.pi**2)
452            FPPI = C*CX*RX/(4.0*math.pi)
453            FP += FPI
454            FPP += FPPI
455        FP -= ElEterm
456   
457    return (FP, FPP, Mu)
458   
459
Note: See TracBrowser for help on using the repository browser.