source: trunk/GSASIIElem.py @ 411

Last change on this file since 411 was 411, checked in by vondreele, 10 years ago

fix very old (!) bug in psvfcj.for
implement neutron scattering lengths in GSASII including the dozen anomalous scattering isotopes
isotope choice is in General
so GSASII will now do neutron CW Rietveld refinements
some cleanup of the constraints GUI
remove varyList from GSASIImapvars.py globals

  • Property svn:keywords set to Date Author Revision URL Id
File size: 16.3 KB
Line 
1"""Element: functions for element types
2   Copyright: 2008, Robert B. Von Dreele (Argonne National Laboratory)
3"""
4########### SVN repository information ###################
5# $Date: 2011-11-11 20:14:35 +0000 (Fri, 11 Nov 2011) $
6# $Author: vondreele $
7# $Revision: 411 $
8# $URL: trunk/GSASIIElem.py $
9# $Id: GSASIIElem.py 411 2011-11-11 20:14:35Z vondreele $
10########### SVN repository information ###################
11
12import wx
13import math
14import sys
15import os.path
16import wx.lib.colourselect as wscs
17import GSASIIpath
18import numpy as np
19
20def GetFormFactorCoeff(El):
21    """Read X-ray form factor coefficients from `atomdata.asc` file
22
23    :param El: element 1-2 character symbol case irrevelant
24    :return: `FormFactors`: list of form factor dictionaries
25   
26    Each X-ray form factor dictionary is:
27   
28    * `Symbol`: 4 character element symbol with valence (e.g. 'NI+2')
29    * `Z`: atomic number
30    * `fa`: 4 A coefficients
31    * `fb`: 4 B coefficients
32    * `fc`: C coefficient
33   
34    """
35    ElS = El.upper()
36    ElS = ElS.rjust(2)
37    filename = os.path.join(os.path.split(sys.argv[0])[0],'atmdata.dat')
38    try:
39        FFdata = open(filename,'Ur')
40    except:
41        wx.MessageBox(message="File atmdata.dat not found in directory %s" % sys.path[0],
42            caption="No atmdata.dat file",style=wx.OK | wx.ICON_EXCLAMATION | wx.STAY_ON_TOP)
43        sys.exit()
44    S = '1'
45    FormFactors = []
46    while S:
47        S = FFdata.readline()
48        if S[3:5] == ElS:
49            if S[5:6] != '_' and S[8] not in ['N','M']:
50                Z=int(S[:2])
51                Symbol = S[3:7].strip()
52                S = S[12:]
53                fa = (float(S[:7]),float(S[14:21]),float(S[28:35]),float(S[42:49]))
54                fb = (float(S[7:14]),float(S[21:28]),float(S[35:42]),float(S[49:56]))
55                FormFac = {'Symbol':Symbol,'Z':Z,'fa':fa,'fb':fb,'fc':float(S[56:63])}
56                FormFactors.append(FormFac)               
57    FFdata.close()
58    return FormFactors
59   
60def GetFFC5(ElSym):
61    '''Get 5 term form factor and Compton scattering data
62    @param ElSym: str(1-2 character element symbol with proper case);
63    @return El: dictionary with 5 term form factor & compton coefficients
64    '''
65    import FormFactors as FF
66    El = {}
67    FF5 = FF.FFac5term[ElSym]
68    El['fa'] = FF5[:5]
69    El['fc'] = FF5[5]
70    El['fb'] = FF5[6:]
71    Cmp5 = FF.Compton[ElSym]
72    El['cmpz'] = Cmp5[0]
73    El['cmpa'] = Cmp5[1:6]
74    El['cmpb'] = Cmp5[6:]
75    return El
76   
77def GetAtomInfo(El):
78   
79    import ElementTable as ET
80    Elements = []
81    for elem in ET.ElTable:
82        Elements.append(elem[0][0])
83    if len(El) in [2,4]:
84        ElS = El.upper()[:2].rjust(2)
85    else:
86        ElS = El.upper()[:1].rjust(2)
87    filename = os.path.join(os.path.split(sys.argv[0])[0],'atmdata.dat')
88    try:
89        FFdata = open(filename,'Ur')
90    except:
91        wx.MessageBox(message="File atmdata.dat not found in directory %s" % sys.path[0],
92            caption="No atmdata.dat file",style=wx.OK | wx.ICON_EXCLAMATION | wx.STAY_ON_TOP)
93        sys.exit()
94    S = '1'
95    AtomInfo = {}
96    Isotopes = {}
97    Mass = []
98    while S:
99        S = FFdata.readline()
100        if S[3:5] == ElS:
101            if S[5:6] == '_':
102                if not Mass:                                 #picks 1st one; natural abundance or 1st isotope
103                    Mass = float(S[10:19])
104                if S[6] in [' ','1','2','3','4','5','6','7','8','9']:                       
105                    isoName = S[6:9]
106                    if isoName == '   ':
107                        isoName = 'Nat. Abund.'              #natural abundance
108                    if S[76:78] in ['LS','BW']:     #special anomalous scattering length info
109                        St = [S[10:19],S[19:25],S[25:31],S[31:38],S[38:44],S[44:50],
110                            S[50:56],S[56:62],S[62:68],S[68:74],]
111                        Vals = []
112                        for item in St:
113                            if item.strip():
114                                Vals.append(float(item.strip()))
115                        Isotopes[isoName.rstrip()] = Vals                       
116                    else:
117                        Isotopes[isoName.rstrip()] = [float(S[10:19]),float(S[19:25])]
118                elif S[5:9] == '_SIZ':
119                    Z=int(S[:2])
120                    Symbol = S[3:5].strip().lower().capitalize()
121                    Drad = float(S[12:22])
122                    Arad = float(S[22:32])
123                    Vdrad = float(S[32:38])
124                    Color = ET.ElTable[Elements.index(Symbol)][6]
125    FFdata.close()
126    AtomInfo={'Symbol':Symbol,'Isotopes':Isotopes,'Mass':Mass,'Z':Z,'Drad':Drad,'Arad':Arad,'Vdrad':Vdrad,'Color':Color}   
127    return AtomInfo
128     
129def GetXsectionCoeff(El):
130    """Read atom orbital scattering cross sections for fprime calculations via Cromer-Lieberman algorithm
131    @param El: 2 character element symbol
132    @return: Orbs: list of orbitals each a dictionary with detailed orbital information used by FPcalc
133    each dictionary is:
134    'OrbName': Orbital name read from file
135    'IfBe' 0/2 depending on orbital
136    'BindEn': binding energy
137    'BB': BindEn/0.02721
138    'XSectIP': 5 cross section inflection points
139    'ElEterm': energy correction term
140    'SEdge': absorption edge for orbital
141    'Nval': 10/11 depending on IfBe
142    'LEner': 10/11 values of log(energy)
143    'LXSect': 10/11 values of log(cross section)
144    """
145    AU = 2.80022e+7
146    C1 = 0.02721
147    ElS = El.upper()
148    ElS = ElS.ljust(2)
149    filename = os.path.join(os.path.split(sys.argv[0])[0],'Xsect.dat')
150    try:
151        xsec = open(filename,'Ur')
152    except:
153        wx.MessageBox(message="File Xsect.dat not found in directory %s" % sys.path[0],
154            caption="No Xsect.dat file",style=wx.OK | wx.ICON_EXCLAMATION |wx.STAY_ON_TOP)
155        sys.exit()
156    S = '1'
157    Orbs = []
158    while S:
159        S = xsec.readline()
160        if S[:2] == ElS:
161            S = S[:-1]+xsec.readline()[:-1]+xsec.readline()
162            OrbName = S[9:14]
163            S = S[14:]
164            IfBe = int(S[0])
165            S = S[1:]
166            val = S.split()
167            BindEn = float(val[0])
168            BB = BindEn/C1
169            Orb = {'OrbName':OrbName,'IfBe':IfBe,'BindEn':BindEn,'BB':BB}
170            Energy = []
171            XSect = []
172            for i in range(11):
173                Energy.append(float(val[2*i+1]))
174                XSect.append(float(val[2*i+2]))
175            XSecIP = []
176            for i in range(5): XSecIP.append(XSect[i+5]/AU)
177            Orb['XSecIP'] = XSecIP
178            if IfBe == 0:
179                Orb['SEdge'] = XSect[10]/AU
180                Nval = 11
181            else:
182                Orb['ElEterm'] = XSect[10]
183                del Energy[10]
184                del XSect[10]
185                Nval = 10
186                Orb['SEdge'] = 0.0
187            Orb['Nval'] = Nval
188            D = dict(zip(Energy,XSect))
189            Energy.sort()
190            X = []
191            for key in Energy:
192                X.append(D[key])
193            XSect = X
194            LEner = []
195            LXSect = []
196            for i in range(Nval):
197                LEner.append(math.log(Energy[i]))
198                if XSect[i] > 0.0:
199                    LXSect.append(math.log(XSect[i]))
200                else:
201                    LXSect.append(0.0)
202            Orb['LEner'] = LEner
203            Orb['LXSect'] = LXSect
204            Orbs.append(Orb)
205    xsec.close()
206    return Orbs
207   
208def GetMagFormFacCoeff(El):
209    """Read magnetic form factor data from atomdata.asc file
210    @param El: 2 character element symbol
211    @return: MagFormFactors: list of all magnetic form factors dictionaries for element El.
212    each dictionary contains:
213    'Symbol':Symbol
214    'Z':Z
215    'mfa': 4 MA coefficients
216    'nfa': 4 NA coefficients
217    'mfb': 4 MB coefficients
218    'nfb': 4 NB coefficients
219    'mfc': MC coefficient
220    'nfc': NC coefficient
221    """
222    ElS = El.upper()
223    ElS = ElS.rjust(2)
224    filename = os.path.join(os.path.split(sys.argv[0])[0],'atmdata.dat')
225    try:
226        FFdata = open(filename,'Ur')
227    except:
228        wx.MessageBox(message="File atmdata.dat not found in directory %s" % sys.path[0],
229            caption="No atmdata.dat file",style=wx.OK | wx.ICON_EXCLAMATION |wx.STAY_ON_TOP)
230        sys.exit()
231    S = '1'
232    MagFormFactors = []
233    while S:
234        S = FFdata.readline()
235        if S[3:5] == ElS:
236            if S[8:9] == 'M':
237                SN = FFdata.readline()               #'N' is assumed to follow 'M' in Atomdata.asc
238                Z=int(S[:2])
239                Symbol = S[3:7]
240                S = S[12:]
241                SN = SN[12:]
242                mfa = (float(S[:7]),float(S[14:21]),float(S[28:35]),float(S[42:49]))
243                mfb = (float(S[7:14]),float(S[21:28]),float(S[35:42]),float(S[49:56]))
244                nfa = (float(SN[:7]),float(SN[14:21]),float(SN[28:35]),float(SN[42:49]))
245                nfb = (float(SN[7:14]),float(SN[21:28]),float(SN[35:42]),float(SN[49:56]))
246                FormFac = {'Symbol':Symbol,'Z':Z,'mfa':mfa,'nfa':nfa,'mfb':mfb,'nfb':nfb,
247                    'mfc':float(S[56:63]),'nfc':float(SN[56:63])}
248                MagFormFactors.append(FormFac)
249    FFdata.close()
250    return MagFormFactors
251
252def ScatFac(El, SQ):
253    """compute value of form factor
254    @param El: element dictionary defined in GetFormFactorCoeff
255    @param SQ: (sin-theta/lambda)**2
256    @return: real part of form factor
257    """
258    fa = np.array(El['fa'])
259    fb = np.array(El['fb'])
260    t = -fb[:,np.newaxis]*SQ
261    return np.sum(fa[:,np.newaxis]*np.exp(t)[:],axis=0)+El['fc']
262       
263def BlenRes(BLdata,wave):
264    FP = []
265    FPP = []
266    Emev = 81.80703/wave**2
267    for BL in BLdata:
268        if len(BL) >= 6:
269            Emev = 81.80703/wave**2
270            G2 = BL[5]**2
271            T = [Emev-BL[4],0,0]
272            D = [T**2+G2,0,0]
273            fp = T/D
274            fpp = 1.0/D
275            if len(BL) == 8:
276                T = Emev-BL[7]
277                D = T**2+G2
278                fp += BL[6]*T/D
279                fpp += BL[6]/D
280            if len(BL) == 10:
281                T = Emev-BL[9]
282                D = T**2+G2
283                fp += BL[8]*T/D
284                fpp += BL[8]/D
285            FP.append(BL[2]*fp)
286            FPP.append(-BL[3]*fpp)
287        else:
288            FP.append(0.0)
289            FPP.append(0.0)
290    return np.array(FP),np.array(FPP)
291   
292def ComptonFac(El,SQ):
293    """compute Compton scattering factor
294    @param El: element dictionary
295    @param SQ: (sin-theta/lambda)**2
296    @return: compton scattering factor
297    """   
298    ca = np.array(El['cmpa'])
299    cb = np.array(El['cmpb'])
300    t = -cb[:,np.newaxis]*SQ       
301    return El['cmpz']-np.sum(ca[:,np.newaxis]*np.exp(t),axis=0) 
302           
303def FPcalc(Orbs, KEv):
304    """Compute real & imaginary resonant X-ray scattering factors
305    @param Orbs: list of orbital dictionaries as defined in GetXsectionCoeff
306    @param KEv: x-ray energy in keV
307    @return: C: (f',f",mu): real, imaginary parts of resonant scattering & atomic absorption coeff.
308    """
309    def Aitken(Orb, LKev):
310        Nval = Orb['Nval']
311        j = Nval-1
312        LEner = Orb['LEner']
313        for i in range(Nval):
314            if LEner[i] <= LKev: j = i
315        if j > Nval-3: j= Nval-3
316        T = [0,0,0,0,0,0]
317        LXSect = Orb['LXSect']
318        for i in range(3):
319           T[i] = LXSect[i+j]
320           T[i+3] = LEner[i+j]-LKev
321        T[1] = (T[0]*T[4]-T[1]*T[3])/(LEner[j+1]-LEner[j])
322        T[2] = (T[0]*T[5]-T[2]*T[3])/(LEner[j+2]-LEner[j])
323        T[2] = (T[1]*T[5]-T[2]*T[4])/(LEner[j+2]-LEner[j+1])
324        C = T[2]
325        return C
326   
327    def DGauss(Orb,CX,RX,ISig):
328        ALG = (0.11846344252810,0.23931433524968,0.284444444444,
329        0.23931433524968,0.11846344252810)
330        XLG = (0.04691007703067,0.23076534494716,0.5,
331        0.76923465505284,0.95308992296933)
332       
333        D = 0.0
334        B2 = Orb['BB']**2
335        R2 = RX**2
336        XSecIP = Orb['XSecIP']
337        for i in range(5):
338            X = XLG[i]
339            X2 = X**2
340            XS = XSecIP[i]
341            if ISig == 0:
342                S = BB*(XS*(B2/X2)-CX*R2)/(R2*X2-B2)
343            elif ISig == 1:
344                S = 0.5*BB*B2*XS/(math.sqrt(X)*(R2*X2-X*B2))
345            elif ISig == 2:
346                T = X*X2*R2-B2/X
347                S = 2.0*BB*(XS*B2/(T*X2**2)-(CX*R2/T))
348            else:
349                S = BB*B2*(XS-Orb['SEdge']*X2)/(R2*X2**2-X2*B2)
350            A = ALG[i]
351            D += A*S
352        return D
353   
354    AU = 2.80022e+7
355    C1 = 0.02721
356    C = 137.0367
357    FP = 0.0
358    FPP = 0.0
359    Mu = 0.0
360    LKev = math.log(KEv)
361    RX = KEv/C1
362    if Orbs:
363        for Orb in Orbs:
364            CX = 0.0
365            BB = Orb['BB']
366            BindEn = Orb['BindEn']
367            if Orb['IfBe'] != 0: ElEterm = Orb['ElEterm']
368            if BindEn <= KEv:
369                CX = math.exp(Aitken(Orb,LKev))
370                Mu += CX
371                CX /= AU
372            Corr = 0.0
373            if Orb['IfBe'] == 0 and BindEn >= KEv:
374                CX = 0.0
375                FPI = DGauss(Orb,CX,RX,3)
376                Corr = 0.5*Orb['SEdge']*BB**2*math.log((RX-BB)/(-RX-BB))/RX
377            else:
378                FPI = DGauss(Orb,CX,RX,Orb['IfBe'])
379                if CX != 0.0: Corr = -0.5*CX*RX*math.log((RX+BB)/(RX-BB))
380            FPI = (FPI+Corr)*C/(2.0*math.pi**2)
381            FPPI = C*CX*RX/(4.0*math.pi)
382            FP += FPI
383            FPP += FPPI
384        FP -= ElEterm
385   
386    return (FP, FPP, Mu)
387   
388
389class PickElement(wx.Dialog):
390    "Makes periodic table widget for picking element - caller maintains element list"
391    Elem=None
392    def _init_ctrls(self, prnt,oneOnly):
393        wx.Dialog.__init__(self, id=-1, name='PickElement',
394              parent=prnt, pos=wx.DefaultPosition, 
395              style=wx.DEFAULT_DIALOG_STYLE, title='Pick Element')
396        import ElementTable as ET
397        self.SetClientSize(wx.Size(770, 250))
398       
399        i=0
400        for E in ET.ElTable:
401            if oneOnly:
402                color=E[4]
403            else:
404                color=E[6]
405            PickElement.ElButton(self,name=E[0],
406               pos=wx.Point(E[1]*40+25,E[2]*24+24),tip=E[3],color=color,oneOnly=oneOnly)
407            i+=1
408
409    def __init__(self, parent,oneOnly=False):
410        self._init_ctrls(parent,oneOnly)
411       
412    def ElButton(self, name, pos, tip, color, oneOnly):
413        Black = wx.Colour(0,0,0)
414        if oneOnly:
415            El = wscs.ColourSelect(label=name[0], parent=self,colour=color,
416                pos=pos, size=wx.Size(40,23), style=wx.RAISED_BORDER)
417#            El.SetLabel(name)
418            El.Bind(wx.EVT_BUTTON, self.OnElButton)
419        else:
420            El = wx.ComboBox(choices=name, parent=self, pos=pos, size=wx.Size(40,23),
421                style=wx.CB_READONLY, value=name[0])
422            El.Bind(wx.EVT_COMBOBOX,self.OnElButton)
423       
424        El.SetBackgroundColour(color)
425        El.SetToolTipString(tip)
426
427    def OnElButton(self, event):
428        El = event.GetEventObject().GetLabel()
429        self.Elem = El
430        self.EndModal(wx.ID_OK)       
431       
432class DeleteElement(wx.Dialog):
433    "Delete element from selected set widget"
434    def _init_ctrls(self, parent,choice):
435        l = len(choice)-1
436        wx.Dialog.__init__(self, id=-1, name='Delete', parent=parent, 
437              pos=wx.DefaultPosition, size=wx.Size(max(128,64+l*24), 87),
438              style=wx.DEFAULT_DIALOG_STYLE, title='Delete Element')
439        self.Show(True)
440        self.SetAutoLayout(True)
441        self.SetHelpText('Select element to delete')
442        self.SetWindowVariant(wx.WINDOW_VARIANT_SMALL)
443
444        i = 0
445        Elem = []
446        for Elem in choice:
447            self.ElButton(id=-1,name=Elem,pos=wx.Point(16+i*24, 16))
448            i+=1
449             
450    def __init__(self, parent,choice):
451        DeleteElement.El = ' '
452        self._init_ctrls(parent,choice)
453
454    def ElButton(self, id, name, pos):
455        White = wx.Colour(255, 255, 255)
456        El = wscs.ColourSelect(label=name, parent=self, colour = White,
457            pos=pos, size=wx.Size(24, 23), style=wx.RAISED_BORDER)
458        El.Bind(wx.EVT_BUTTON, self.OnDeleteButton)
459   
460    def OnDeleteButton(self, event):
461        DeleteElement.El=event.GetEventObject().GetLabel()
462        self.EndModal(wx.ID_OK)
463       
464    def GetDeleteElement(self):
465        return DeleteElement.El
466       
467
Note: See TracBrowser for help on using the repository browser.