source: trunk/GSASIIElem.py @ 327

Last change on this file since 327 was 327, checked in by vondreele, 12 years ago

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