source: trunk/GSASIIElem.py @ 155

Last change on this file since 155 was 155, checked in by vondreele, 11 years ago

changes for structure drawing

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