source: trunk/GSASIIElem.py @ 2486

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

fix mag form factor lookup - some aren't chemical valences (e.g. Fe+4)
fix mag form factor periodic table - nonmagnetic atoms now not shown
fix powder reflection mark & diff curve placement issues
fix issues with mag structure drawings - fill unit cell, etc.
fix structure transform to make magnetic cells (still needs making of constraints)
fix problem of mustrain & hstrain coeff changes with change in space group
add print to file of covariance matrix - user request
fix magnetic moment site symmetry restrictions
work on mag structure factor calcs.
change EXP file importer to ignore magnetic symmetry from GSAS

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