source: trunk/GSASIIElem.py @ 1344

Last change on this file since 1344 was 1295, checked in by vondreele, 11 years ago

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