source: trunk/GSASIIElem.py @ 3191

Last change on this file since 3191 was 3191, checked in by vondreele, 6 years ago

fix integer divides in responding to unit cell changes on Unit Cells window
show StarFile? errors from cif importer
add missing nonstandard space group symbols for face centered orthorhombics
begin implementation of cif import magnetic phases into 2 phase result
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 16.3 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: 2017-12-13 14:26:51 +0000 (Wed, 13 Dec 2017) $
10# $Author: vondreele $
11# $Revision: 3191 $
12# $URL: trunk/GSASIIElem.py $
13# $Id: GSASIIElem.py 3191 2017-12-13 14:26:51Z vondreele $
14########### SVN repository information ###################
15
16import math
17import sys
18import os.path
19import GSASIIpath
20GSASIIpath.SetVersionNumber("$Revision: 3191 $")
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 H')
200            ElS = 'H'
201#            return # not sure what this element should be!
202        print('Atom type '+El+' not found, using '+ElS)
203        El = ElS
204    AtomInfo.update(dict(zip(['Drad','Arad','Vdrad','Hbrad'],atmdata.AtmSize[ElS])))
205    AtomInfo['Symbol'] = El
206    AtomInfo['Color'] = ET.ElTable[Elements.index(ElS)][6]
207    AtomInfo['Z'] = atmdata.XrayFF[ElS]['Z']
208    isotopes = [ky for ky in atmdata.AtmBlens.keys() if ElS == ky.split('_')[0]]
209    isotopes.sort()
210    AtomInfo['Mass'] = atmdata.AtmBlens[isotopes[0]]['Mass']    #default to nat. abund.
211    AtomInfo['Isotopes'] = {}
212    for isotope in isotopes:
213        data = atmdata.AtmBlens[isotope]
214        if isotope == ElS+'_':
215            AtomInfo['Isotopes']['Nat. Abund.'] = data
216        else:
217            AtomInfo['Isotopes'][isotope.split('_')[1]] = data
218    AtomInfo['Lande g'] = 2.0
219    return AtomInfo
220   
221def GetElInfo(El,inst):
222    ElemSym = El.strip().capitalize()
223    if 'X' in inst['Type'][0]: 
224        keV = 12.397639/G2mth.getWave(inst)               
225        FpMu = FPcalc(GetXsectionCoeff(ElemSym), keV)
226        ElData = GetFormFactorCoeff(ElemSym)[0]
227        ElData['FormulaNo'] = 0.0
228        ElData.update(GetAtomInfo(ElemSym))
229        ElData.update(dict(zip(['fp','fpp','mu'],FpMu)))
230        ElData.update(GetFFC5(El))
231    else: #'N'eutron
232        ElData = {}
233        ElData.update(GetAtomInfo(ElemSym))
234        ElData['FormulaNo'] = 0.0
235        ElData.update({'mu':0.0,'fp':0.0,'fpp':0.0})
236    return ElData
237       
238def GetXsectionCoeff(El):
239    """Read atom orbital scattering cross sections for fprime calculations via Cromer-Lieberman algorithm
240
241    :param El: 2 character element symbol
242    :return: Orbs: list of orbitals each a dictionary with detailed orbital information used by FPcalc
243
244    each dictionary is:
245
246    * 'OrbName': Orbital name read from file
247    * 'IfBe' 0/2 depending on orbital
248    * 'BindEn': binding energy
249    * 'BB': BindEn/0.02721
250    * 'XSectIP': 5 cross section inflection points
251    * 'ElEterm': energy correction term
252    * 'SEdge': absorption edge for orbital
253    * 'Nval': 10/11 depending on IfBe
254    * 'LEner': 10/11 values of log(energy)
255    * 'LXSect': 10/11 values of log(cross section)
256
257    """
258    AU = 2.80022e+7
259    C1 = 0.02721
260    ElS = El.upper()
261    ElS = ElS.ljust(2)
262    filename = os.path.join(os.path.split(__file__)[0],'Xsect.dat')
263    try:
264        xsec = open(filename,'Ur')
265    except:
266        print ('**** ERROR - File Xsect.dat not found in directory %s'%os.path.split(filename)[0])
267        sys.exit()
268    S = '1'
269    Orbs = []
270    while S:
271        S = xsec.readline()
272        if S[:2] == ElS:
273            S = S[:-1]+xsec.readline()[:-1]+xsec.readline()
274            OrbName = S[9:14]
275            S = S[14:]
276            IfBe = int(S[0])
277            S = S[1:]
278            val = S.split()
279            BindEn = float(val[0])
280            BB = BindEn/C1
281            Orb = {'OrbName':OrbName,'IfBe':IfBe,'BindEn':BindEn,'BB':BB}
282            Energy = []
283            XSect = []
284            for i in range(11):
285                Energy.append(float(val[2*i+1]))
286                XSect.append(float(val[2*i+2]))
287            XSecIP = []
288            for i in range(5): XSecIP.append(XSect[i+5]/AU)
289            Orb['XSecIP'] = XSecIP
290            if IfBe == 0:
291                Orb['SEdge'] = XSect[10]/AU
292                Nval = 11
293            else:
294                Orb['ElEterm'] = XSect[10]
295                del Energy[10]
296                del XSect[10]
297                Nval = 10
298                Orb['SEdge'] = 0.0
299            Orb['Nval'] = Nval
300            D = dict(zip(Energy,XSect))
301            Energy.sort()
302            X = []
303            for key in Energy:
304                X.append(D[key])
305            XSect = X
306            LEner = []
307            LXSect = []
308            for i in range(Nval):
309                LEner.append(math.log(Energy[i]))
310                if XSect[i] > 0.0:
311                    LXSect.append(math.log(XSect[i]))
312                else:
313                    LXSect.append(0.0)
314            Orb['LEner'] = LEner
315            Orb['LXSect'] = LXSect
316            Orbs.append(Orb)
317    xsec.close()
318    return Orbs
319   
320def GetMagFormFacCoeff(El):
321    """Read magnetic form factor data from atmdata.py
322
323    :param El: 2 character element symbol
324    :return: MagFormFactors: list of all magnetic form factors dictionaries for element El.
325
326    each dictionary contains:
327
328    * 'Symbol':Symbol
329    * 'Z':Z
330    * 'mfa': 4 MA coefficients
331    * 'nfa': 4 NA coefficients
332    * 'mfb': 4 MB coefficients
333    * 'nfb': 4 NB coefficients
334    * 'mfc': MC coefficient
335    * 'nfc': NC coefficient
336   
337    """
338    Els = El.capitalize().strip()
339    MagFormFactors = []
340    mags = [ky for ky in atmdata.MagFF.keys() if Els == getElSym(ky)]
341    for mag in mags:
342        magData = {}
343        data = atmdata.MagFF[mag]
344        magData['Symbol'] = mag.upper()
345        magData['Z'] = atmdata.XrayFF[getElSym(mag)]['Z']
346        magData['mfa'] = [data['M'][i] for i in [0,2,4,6]]
347        magData['mfb'] = [data['M'][i] for i in [1,3,5,7]]
348        magData['mfc'] = data['M'][8]
349        magData['nfa'] = [data['N'][i] for i in [0,2,4,6]]
350        magData['nfb'] = [data['N'][i] for i in [1,3,5,7]]
351        magData['nfc'] = data['N'][8]
352        MagFormFactors.append(magData)
353    return MagFormFactors
354
355def ScatFac(El, SQ):
356    """compute value of form factor
357
358    :param El: element dictionary defined in GetFormFactorCoeff
359    :param SQ: (sin-theta/lambda)**2
360    :return: real part of form factor
361    """
362    fa = np.array(El['fa'])
363    fb = np.array(El['fb'])
364    t = -fb[:,np.newaxis]*SQ
365    return np.sum(fa[:,np.newaxis]*np.exp(t)[:],axis=0)+El['fc']
366       
367def MagScatFac(El, SQ):
368    """compute value of form factor
369
370    :param El: element dictionary defined in GetFormFactorCoeff
371    :param SQ: (sin-theta/lambda)**2
372    :param gfac: Lande g factor (normally = 2.0)
373    :return: real part of form factor
374    """
375    mfa = np.array(El['mfa'])
376    mfb = np.array(El['mfb'])
377    nfa = np.array(El['nfa'])
378    nfb = np.array(El['nfb'])
379    mt = -mfb[:,np.newaxis]*SQ
380    nt = -nfb[:,np.newaxis]*SQ
381    MMF = np.sum(mfa[:,np.newaxis]*np.exp(mt)[:],axis=0)+El['mfc']
382    NMF = np.sum(nfa[:,np.newaxis]*np.exp(nt)[:],axis=0)+El['nfc']
383    return MMF+(2.0/El['gfac']-1.0)*NMF
384       
385def BlenResCW(Els,BLtables,wave):
386    FP = np.zeros(len(Els))
387    FPP = np.zeros(len(Els))
388    for i,El in enumerate(Els):
389        BL = BLtables[El][1]
390        if 'BW-LS' in BL:
391            Re,Im,E0,gam,A,E1,B,E2 = BL['BW-LS'][1:]
392            Emev = 81.80703/wave**2
393            T0 = Emev-E0
394            T1 = Emev-E1
395            T2 = Emev-E2
396            D0 = T0**2+gam**2
397            D1 = T1**2+gam**2
398            D2 = T2**2+gam**2
399            FP[i] = Re*(T0/D0+A*T1/D1+B*T2/D2)+BL[i]['BW-LS'][0]
400            FPP[i] = -Im*(1/D0+A/D1+B/D2)
401        else:
402            FPP[i] = BL['SL'][1]    #for Li, B, etc.
403    return FP,FPP
404   
405def BlenResTOF(Els,BLtables,wave):
406    FP = np.zeros((len(Els),len(wave)))
407    FPP = np.zeros((len(Els),len(wave)))
408    BL = [BLtables[el][1] for el in Els]
409    for i,El in enumerate(Els):
410        if 'BW-LS' in BL[i]:
411            Re,Im,E0,gam,A,E1,B,E2 = BL[i]['BW-LS'][1:]
412            Emev = 81.80703/wave**2
413            T0 = Emev-E0
414            T1 = Emev-E1
415            T2 = Emev-E2
416            D0 = T0**2+gam**2
417            D1 = T1**2+gam**2
418            D2 = T2**2+gam**2
419            FP[i] = Re*(T0/D0+A*T1/D1+B*T2/D2)+BL[i]['BW-LS'][0]
420            FPP[i] = -Im*(1/D0+A/D1+B/D2)
421        else:
422            FPP[i] = np.ones(len(wave))*BL[i]['SL'][1]    #for Li, B, etc.
423    return FP,FPP
424   
425def ComptonFac(El,SQ):
426    """compute Compton scattering factor
427
428    :param El: element dictionary
429    :param SQ: (sin-theta/lambda)**2
430    :return: compton scattering factor
431    """   
432    ca = np.array(El['cmpa'])
433    cb = np.array(El['cmpb'])
434    t = -cb[:,np.newaxis]*SQ       
435    return El['cmpz']-np.sum(ca[:,np.newaxis]*np.exp(t),axis=0) 
436           
437def FPcalc(Orbs, KEv):
438    """Compute real & imaginary resonant X-ray scattering factors
439
440    :param Orbs: list of orbital dictionaries as defined in GetXsectionCoeff
441    :param KEv: x-ray energy in keV
442    :return: C: (f',f",mu): real, imaginary parts of resonant scattering & atomic absorption coeff.
443    """
444    def Aitken(Orb, LKev):
445        Nval = Orb['Nval']
446        j = Nval-1
447        LEner = Orb['LEner']
448        for i in range(Nval):
449            if LEner[i] <= LKev: j = i
450        if j > Nval-3: j= Nval-3
451        T = [0,0,0,0,0,0]
452        LXSect = Orb['LXSect']
453        for i in range(3):
454           T[i] = LXSect[i+j]
455           T[i+3] = LEner[i+j]-LKev
456        T[1] = (T[0]*T[4]-T[1]*T[3])/(LEner[j+1]-LEner[j])
457        T[2] = (T[0]*T[5]-T[2]*T[3])/(LEner[j+2]-LEner[j])
458        T[2] = (T[1]*T[5]-T[2]*T[4])/(LEner[j+2]-LEner[j+1])
459        C = T[2]
460        return C
461   
462    def DGauss(Orb,CX,RX,ISig):
463        ALG = (0.11846344252810,0.23931433524968,0.284444444444,
464        0.23931433524968,0.11846344252810)
465        XLG = (0.04691007703067,0.23076534494716,0.5,
466        0.76923465505284,0.95308992296933)
467       
468        D = 0.0
469        B2 = Orb['BB']**2
470        R2 = RX**2
471        XSecIP = Orb['XSecIP']
472        for i in range(5):
473            X = XLG[i]
474            X2 = X**2
475            XS = XSecIP[i]
476            if ISig == 0:
477                S = BB*(XS*(B2/X2)-CX*R2)/(R2*X2-B2)
478            elif ISig == 1:
479                S = 0.5*BB*B2*XS/(math.sqrt(X)*(R2*X2-X*B2))
480            elif ISig == 2:
481                T = X*X2*R2-B2/X
482                S = 2.0*BB*(XS*B2/(T*X2**2)-(CX*R2/T))
483            else:
484                S = BB*B2*(XS-Orb['SEdge']*X2)/(R2*X2**2-X2*B2)
485            A = ALG[i]
486            D += A*S
487        return D
488   
489    AU = 2.80022e+7
490    C1 = 0.02721
491    C = 137.0367
492    FP = 0.0
493    FPP = 0.0
494    Mu = 0.0
495    LKev = math.log(KEv)
496    RX = KEv/C1
497    if Orbs:
498        for Orb in Orbs:
499            CX = 0.0
500            BB = Orb['BB']
501            BindEn = Orb['BindEn']
502            if Orb['IfBe'] != 0: ElEterm = Orb['ElEterm']
503            if BindEn <= KEv:
504                CX = math.exp(Aitken(Orb,LKev))
505                Mu += CX
506                CX /= AU
507            Corr = 0.0
508            if Orb['IfBe'] == 0 and BindEn >= KEv:
509                CX = 0.0
510                FPI = DGauss(Orb,CX,RX,3)
511                Corr = 0.5*Orb['SEdge']*BB**2*math.log((RX-BB)/(-RX-BB))/RX
512            else:
513                FPI = DGauss(Orb,CX,RX,Orb['IfBe'])
514                if CX != 0.0: Corr = -0.5*CX*RX*math.log((RX+BB)/(RX-BB))
515            FPI = (FPI+Corr)*C/(2.0*math.pi**2)
516            FPPI = C*CX*RX/(4.0*math.pi)
517            FP += FPI
518            FPP += FPPI
519        FP -= ElEterm
520   
521    return (FP, FPP, Mu)
522   
523
Note: See TracBrowser for help on using the repository browser.