source: trunk/GSASIIElem.py @ 2473

Last change on this file since 2473 was 2473, checked in by vondreele, 5 years ago

work on magnetic structures - import from EXP, plotting, LS refine I/O, mag. form factors, etc.

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