source: trunk/GSASIIElem.py @ 4756

Last change on this file since 4756 was 4659, checked in by toby, 5 years ago

extensive Origin 1->2 updates; merge SetupGeneral? into a single routine in G2Elem; read sym ops from CIF; fix sites sym/mult after cood xform
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 24.6 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: 2020-11-21 20:15:49 +0000 (Sat, 21 Nov 2020) $
10# $Author: toby $
11# $Revision: 4659 $
12# $URL: trunk/GSASIIElem.py $
13# $Id: GSASIIElem.py 4659 2020-11-21 20:15:49Z toby $
14########### SVN repository information ###################
15
16import math
17import sys
18import os.path
19import GSASIIpath
20GSASIIpath.SetVersionNumber("$Revision: 4659 $")
21import numpy as np
22import atmdata
23import GSASIImath as G2mth
24import ElementTable as ET
25
26
27getElSym = lambda sym: sym.split('+')[0].split('-')[0].capitalize()
28def GetFormFactorCoeff(El):
29    """Read X-ray form factor coefficients from `atomdata.py` file
30
31    :param str El: element 1-2 character symbol, case irrevelant
32    :return: `FormFactors`: list of form factor dictionaries
33   
34    Each X-ray form factor dictionary is:
35   
36    * `Symbol`: 4 character element symbol with valence (e.g. 'NI+2')
37    * `Z`: atomic number
38    * `fa`: 4 A coefficients
39    * `fb`: 4 B coefficients
40    * `fc`: C coefficient
41   
42    """
43   
44    Els = El.capitalize().strip()
45    valences = [ky for ky in atmdata.XrayFF.keys() if Els == getElSym(ky)]
46    FormFactors = [atmdata.XrayFF[val] for val in valences]
47    for Sy,FF in zip(valences,FormFactors):
48        FF.update({'Symbol':Sy.upper()})
49    return FormFactors
50   
51def GetFFtable(atomTypes):
52    ''' returns a dictionary of form factor data for atom types found in atomTypes
53
54    :param list atomTypes: list of atom types
55    :return: FFtable, dictionary of form factor data; key is atom type
56
57    '''
58    FFtable = {}
59    for El in atomTypes:
60        FFs = GetFormFactorCoeff(getElSym(El))
61        for item in FFs:
62            if item['Symbol'] == El.upper():
63                FFtable[El] = item
64    return FFtable
65   
66def GetMFtable(atomTypes,Landeg):
67    ''' returns a dictionary of magnetic form factor data for atom types found in atomTypes
68
69    :param list atomTypes: list of atom types
70    :param list Landeg: Lande g factors for atomTypes
71    :return: FFtable, dictionary of form factor data; key is atom type
72
73    '''
74    MFtable = {}
75    for El,gfac in zip(atomTypes,Landeg):
76        MFs = GetMagFormFacCoeff(getElSym(El))
77        for item in MFs:
78            if item['Symbol'] == El.upper():
79                item['gfac'] = gfac
80                MFtable[El] = item
81    return MFtable
82   
83def GetBLtable(General):
84    ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
85
86    :param dict General: dictionary of phase info.; includes AtomTypes & Isotopes
87    :return: BLtable, dictionary of scattering length data; key is atom type
88    '''
89    atomTypes = General['AtomTypes']
90    BLtable = {}
91    isotope = General['Isotope']
92    for El in atomTypes:
93        ElS = getElSym(El)
94        if 'Nat' in isotope[El]:
95            BLtable[El] = [isotope[El],atmdata.AtmBlens[ElS+'_']]
96        else:
97            BLtable[El] = [isotope[El],atmdata.AtmBlens[ElS+'_'+isotope[El]]]
98    return BLtable
99       
100def getFFvalues(FFtables,SQ,ifList=False):
101    'Needs a doc string'
102    if ifList:
103        FFvals = []
104        for El in FFtables:
105            FFvals.append(ScatFac(FFtables[El],SQ)[0])
106    else:
107        FFvals = {}
108        for El in FFtables:
109            FFvals[El] = ScatFac(FFtables[El],SQ)[0]
110    return FFvals
111   
112def getBLvalues(BLtables,ifList=False):
113    'Needs a doc string'
114    if ifList:
115        BLvals = []
116        for El in BLtables:
117            if 'BW-LS' in El:
118                BLvals.append(BLtables[El][1]['BW-LS'][0])
119            else:
120                BLvals.append(BLtables[El][1]['SL'][0])
121    else:
122        BLvals = {}
123        for El in BLtables:
124            if 'BW-LS' in El:
125                BLvals[El] = BLtables[El][1]['BW-LS'][0]
126            else:
127                BLvals[El] = BLtables[El][1]['SL'][0]
128    return BLvals
129       
130def getMFvalues(MFtables,SQ,ifList=False):
131    'Needs a doc string'
132    if ifList:
133        MFvals = []
134        for El in MFtables:
135            MFvals.append(MagScatFac(MFtables[El],SQ)[0])
136    else:
137        MFvals = {}
138        for El in MFtables:
139            MFvals[El] = MagScatFac(MFtables[El],SQ)[0]
140    return MFvals
141   
142def GetFFC5(ElSym):
143    '''Get 5 term form factor and Compton scattering data
144
145    :param ElSym: str(1-2 character element symbol with proper case);
146    :return El: dictionary with 5 term form factor & compton coefficients
147    '''
148    import FormFactors as FF
149    El = {}
150    FF5 = FF.FFac5term[ElSym]
151    El['fa'] = FF5[:5]
152    El['fc'] = FF5[5]
153    El['fb'] = FF5[6:]
154    Cmp5 = FF.Compton[ElSym]
155    El['cmpz'] = Cmp5[0]
156    El['cmpa'] = Cmp5[1:6]
157    El['cmpb'] = Cmp5[6:]
158    return El
159
160def GetBVS(Pair,atSeq,Valences):
161    Els = Pair.strip().split('-')
162    iAt = atSeq.index(Els[0])
163    iVal = Valences[iAt][0]
164    if Els[1] in ['O','F','Cl']:
165        iEls = ['O','F','Cl'].index(Els[1])
166        if iVal in atmdata.BVScoeff:
167            return atmdata.BVScoeff[iVal][iEls]
168        else:
169            return 0.0
170    elif Els[1] in ['Br','I','S','Se','Te','N','P','As','H','D']:
171        iEls = ['Br','I','S','Se','Te','N','P','As','H','D'].index(Els[1])
172        if Els[0] in atmdata.BVSnotOFCl:
173            return atmdata.BVSnotOFCl[Els[0]][iEls]
174        else:
175            return 0.0
176    else:
177        return 0.0
178   
179def CheckElement(El):
180    '''Check if element El is in the periodic table
181
182    :param str El: One or two letter element symbol, capitaliztion ignored
183    :returns: True if the element is found
184
185    '''
186    Elements = []
187    for elem in ET.ElTable:
188        Elements.append(elem[0][0])
189    if El.capitalize() in Elements:
190        return True
191    else:
192        return False 
193
194def FixValence(El):
195    'Returns the element symbol, even when a valence is present'
196    if '+' in El[-1]: #converts An+/- to A+/-n
197        num = El[-2]
198        El = El.split(num)[0]+'+'+num
199    if '+0' in El:
200        El = El.split('+0')[0]
201    if '-' in El[-1]:
202        num = El[-2]
203        El = El.split(num)[0]+'-'+num
204    if '-0' in El:
205        El = El.split('-0')[0]
206    return El
207   
208def GetAtomInfo(El,ifMag=False):
209    'reads element information from atmdata.py'
210    Elem = ET.ElTable
211    if ifMag:
212        Elem = ET.MagElTable
213    Elements = [elem[0][0] for elem in Elem]
214    AtomInfo = {}
215    ElS = getElSym(El)
216    if El not in atmdata.XrayFF and El not in atmdata.MagFF:
217        if ElS not in atmdata.XrayFF:
218            print('Atom type '+El+' not found, using H')
219            ElS = 'H'
220#            return # not sure what this element should be!
221        print('Atom type '+El+' not found, using '+ElS)
222        El = ElS
223    AtomInfo.update(dict(zip(['Drad','Arad','Vdrad','Hbrad'],atmdata.AtmSize[ElS])))
224    AtomInfo['Symbol'] = El
225    AtomInfo['Color'] = ET.ElTable[Elements.index(ElS)][6]
226    AtomInfo['Z'] = atmdata.XrayFF[ElS]['Z']
227    isotopes = [ky for ky in atmdata.AtmBlens.keys() if ElS == ky.split('_')[0]]
228    isotopes.sort()
229    AtomInfo['Mass'] = atmdata.AtmBlens[isotopes[0]]['Mass']    #default to nat. abund.
230    AtomInfo['Isotopes'] = {}
231    for isotope in isotopes:
232        data = atmdata.AtmBlens[isotope]
233        if isotope == ElS+'_':
234            AtomInfo['Isotopes']['Nat. Abund.'] = data
235        else:
236            AtomInfo['Isotopes'][isotope.split('_')[1]] = data
237    AtomInfo['Lande g'] = 2.0
238    return AtomInfo
239   
240def GetElInfo(El,inst):
241    ElemSym = El.strip().capitalize()
242    if 'X' in inst['Type'][0]: 
243        keV = 12.397639/G2mth.getWave(inst)               
244        FpMu = FPcalc(GetXsectionCoeff(ElemSym), keV)
245        ElData = GetFormFactorCoeff(ElemSym)[0]
246        ElData['FormulaNo'] = 0.0
247        ElData.update(GetAtomInfo(ElemSym))
248        ElData.update(dict(zip(['fp','fpp','mu'],FpMu)))
249        ElData.update(GetFFC5(El))
250    else: #'N'eutron
251        ElData = {}
252        ElData.update(GetAtomInfo(ElemSym))
253        ElData['FormulaNo'] = 0.0
254        ElData.update({'mu':0.0,'fp':0.0,'fpp':0.0})
255    return ElData
256       
257def GetXsectionCoeff(El):
258    """Read atom orbital scattering cross sections for fprime calculations via Cromer-Lieberman algorithm
259
260    :param El: 2 character element symbol
261    :return: Orbs: list of orbitals each a dictionary with detailed orbital information used by FPcalc
262
263    each dictionary is:
264
265    * 'OrbName': Orbital name read from file
266    * 'IfBe' 0/2 depending on orbital
267    * 'BindEn': binding energy
268    * 'BB': BindEn/0.02721
269    * 'XSectIP': 5 cross section inflection points
270    * 'ElEterm': energy correction term
271    * 'SEdge': absorption edge for orbital
272    * 'Nval': 10/11 depending on IfBe
273    * 'LEner': 10/11 values of log(energy)
274    * 'LXSect': 10/11 values of log(cross section)
275
276    """
277    AU = 2.80022e+7
278    C1 = 0.02721
279    ElS = El.upper()
280    ElS = ElS.ljust(2)
281    filename = os.path.join(os.path.split(__file__)[0],'Xsect.dat')
282    try:
283        xsec = open(filename,'r')
284    except:
285        print ('**** ERROR - File Xsect.dat not found in directory %s'%os.path.split(filename)[0])
286        sys.exit()
287    S = '1'
288    Orbs = []
289    while S:
290        S = xsec.readline()
291        if S[:2] == ElS:
292            S = S[:-1]+xsec.readline()[:-1]+xsec.readline()
293            OrbName = S[9:14]
294            S = S[14:]
295            IfBe = int(S[0])
296            S = S[1:]
297            val = S.split()
298            BindEn = float(val[0])
299            BB = BindEn/C1
300            Orb = {'OrbName':OrbName,'IfBe':IfBe,'BindEn':BindEn,'BB':BB}
301            Energy = []
302            XSect = []
303            for i in range(11):
304                Energy.append(float(val[2*i+1]))
305                XSect.append(float(val[2*i+2]))
306            XSecIP = []
307            for i in range(5): XSecIP.append(XSect[i+5]/AU)
308            Orb['XSecIP'] = XSecIP
309            if IfBe == 0:
310                Orb['SEdge'] = XSect[10]/AU
311                Nval = 11
312            else:
313                Orb['ElEterm'] = XSect[10]
314                del Energy[10]
315                del XSect[10]
316                Nval = 10
317                Orb['SEdge'] = 0.0
318            Orb['Nval'] = Nval
319            D = dict(zip(Energy,XSect))
320            Energy.sort()
321            X = []
322            for key in Energy:
323                X.append(D[key])
324            XSect = X
325            LEner = []
326            LXSect = []
327            for i in range(Nval):
328                LEner.append(math.log(Energy[i]))
329                if XSect[i] > 0.0:
330                    LXSect.append(math.log(XSect[i]))
331                else:
332                    LXSect.append(0.0)
333            Orb['LEner'] = LEner
334            Orb['LXSect'] = LXSect
335            Orbs.append(Orb)
336    xsec.close()
337    return Orbs
338   
339def GetMagFormFacCoeff(El):
340    """Read magnetic form factor data from atmdata.py
341
342    :param El: 2 character element symbol
343    :return: MagFormFactors: list of all magnetic form factors dictionaries for element El.
344
345    each dictionary contains:
346
347    * 'Symbol':Symbol
348    * 'Z':Z
349    * 'mfa': 4 MA coefficients
350    * 'nfa': 4 NA coefficients
351    * 'mfb': 4 MB coefficients
352    * 'nfb': 4 NB coefficients
353    * 'mfc': MC coefficient
354    * 'nfc': NC coefficient
355   
356    """
357    Els = El.capitalize().strip()
358    MagFormFactors = []
359    mags = [ky for ky in atmdata.MagFF.keys() if Els == getElSym(ky)]
360    for mag in mags:
361        magData = {}
362        data = atmdata.MagFF[mag]
363        magData['Symbol'] = mag.upper()
364        magData['Z'] = atmdata.XrayFF[getElSym(mag)]['Z']
365        magData['mfa'] = [data['M'][i] for i in [0,2,4,6]]
366        magData['mfb'] = [data['M'][i] for i in [1,3,5,7]]
367        magData['mfc'] = data['M'][8]
368        magData['nfa'] = [data['N'][i] for i in [0,2,4,6]]
369        magData['nfb'] = [data['N'][i] for i in [1,3,5,7]]
370        magData['nfc'] = data['N'][8]
371        MagFormFactors.append(magData)
372    return MagFormFactors
373
374def ScatFac(El, SQ):
375    """compute value of form factor
376
377    :param El: element dictionary defined in GetFormFactorCoeff
378    :param SQ: (sin-theta/lambda)**2
379    :return: real part of form factor
380    """
381    fa = np.array(El['fa'])
382    fb = np.array(El['fb'])
383    t = -fb[:,np.newaxis]*SQ
384    return np.sum(fa[:,np.newaxis]*np.exp(t)[:],axis=0)+El['fc']
385       
386def MagScatFac(El, SQ):
387    """compute value of form factor
388
389    :param El: element dictionary defined in GetFormFactorCoeff
390    :param SQ: (sin-theta/lambda)**2
391    :param gfac: Lande g factor (normally = 2.0)
392    :return: real part of form factor
393    """
394    mfa = np.array(El['mfa'])
395    mfb = np.array(El['mfb'])
396    nfa = np.array(El['nfa'])
397    nfb = np.array(El['nfb'])
398    mt = -mfb[:,np.newaxis]*SQ
399    nt = -nfb[:,np.newaxis]*SQ
400    MMF = np.sum(mfa[:,np.newaxis]*np.exp(mt)[:],axis=0)+El['mfc']
401    MMF0 = np.sum(mfa)+El['mfc']
402    NMF = np.sum(nfa[:,np.newaxis]*np.exp(nt)[:],axis=0)+El['nfc']
403    NMF0 = np.sum(nfa)+El['nfc']
404    MF0 = MMF0+(2.0/El['gfac']-1.0)*NMF0
405    return (MMF+(2.0/El['gfac']-1.0)*NMF)/MF0
406       
407def BlenResCW(Els,BLtables,wave):
408    FP = np.zeros(len(Els))
409    FPP = np.zeros(len(Els))
410    for i,El in enumerate(Els):
411        BL = BLtables[El][1]
412        if 'BW-LS' in BL:
413            Re,Im,E0,gam,A,E1,B,E2 = BL['BW-LS'][1:]
414            Emev = 81.80703/wave**2
415            T0 = Emev-E0
416            T1 = Emev-E1
417            T2 = Emev-E2
418            D0 = T0**2+gam**2
419            D1 = T1**2+gam**2
420            D2 = T2**2+gam**2
421            FP[i] = Re*(T0/D0+A*T1/D1+B*T2/D2)+BL['BW-LS'][0]
422            FPP[i] = -Im*(1/D0+A/D1+B/D2)
423        else:
424            FPP[i] = BL['SL'][1]    #for Li, B, etc.
425    return FP,FPP
426   
427def BlenResTOF(Els,BLtables,wave):
428    FP = np.zeros((len(Els),len(wave)))
429    FPP = np.zeros((len(Els),len(wave)))
430    BL = [BLtables[el][1] for el in Els]
431    for i,El in enumerate(Els):
432        if 'BW-LS' in BL[i]:
433            Re,Im,E0,gam,A,E1,B,E2 = BL[i]['BW-LS'][1:]
434            Emev = 81.80703/wave**2
435            T0 = Emev-E0
436            T1 = Emev-E1
437            T2 = Emev-E2
438            D0 = T0**2+gam**2
439            D1 = T1**2+gam**2
440            D2 = T2**2+gam**2
441            FP[i] = Re*(T0/D0+A*T1/D1+B*T2/D2)+BL[i]['BW-LS'][0]
442            FPP[i] = -Im*(1/D0+A/D1+B/D2)
443        else:
444            FPP[i] = np.ones(len(wave))*BL[i]['SL'][1]    #for Li, B, etc.
445    return FP,FPP
446   
447def ComptonFac(El,SQ):
448    """compute Compton scattering factor
449
450    :param El: element dictionary
451    :param SQ: (sin-theta/lambda)**2
452    :return: compton scattering factor
453    """   
454    ca = np.array(El['cmpa'])
455    cb = np.array(El['cmpb'])
456    t = -cb[:,np.newaxis]*SQ       
457    return El['cmpz']-np.sum(ca[:,np.newaxis]*np.exp(t),axis=0) 
458           
459def FPcalc(Orbs, KEv):
460    """Compute real & imaginary resonant X-ray scattering factors
461
462    :param Orbs: list of orbital dictionaries as defined in GetXsectionCoeff
463    :param KEv: x-ray energy in keV
464    :return: C: (f',f",mu): real, imaginary parts of resonant scattering & atomic absorption coeff.
465    """
466    def Aitken(Orb, LKev):
467        Nval = Orb['Nval']
468        j = Nval-1
469        LEner = Orb['LEner']
470        for i in range(Nval):
471            if LEner[i] <= LKev: j = i
472        if j > Nval-3: j= Nval-3
473        T = [0,0,0,0,0,0]
474        LXSect = Orb['LXSect']
475        for i in range(3):
476           T[i] = LXSect[i+j]
477           T[i+3] = LEner[i+j]-LKev
478        T[1] = (T[0]*T[4]-T[1]*T[3])/(LEner[j+1]-LEner[j])
479        T[2] = (T[0]*T[5]-T[2]*T[3])/(LEner[j+2]-LEner[j])
480        T[2] = (T[1]*T[5]-T[2]*T[4])/(LEner[j+2]-LEner[j+1])
481        C = T[2]
482        return C
483   
484    def DGauss(Orb,CX,RX,ISig):
485        ALG = (0.11846344252810,0.23931433524968,0.284444444444,
486        0.23931433524968,0.11846344252810)
487        XLG = (0.04691007703067,0.23076534494716,0.5,
488        0.76923465505284,0.95308992296933)
489       
490        D = 0.0
491        B2 = Orb['BB']**2
492        R2 = RX**2
493        XSecIP = Orb['XSecIP']
494        for i in range(5):
495            X = XLG[i]
496            X2 = X**2
497            XS = XSecIP[i]
498            if ISig == 0:
499                S = BB*(XS*(B2/X2)-CX*R2)/(R2*X2-B2)
500            elif ISig == 1:
501                S = 0.5*BB*B2*XS/(math.sqrt(X)*(R2*X2-X*B2))
502            elif ISig == 2:
503                T = X*X2*R2-B2/X
504                S = 2.0*BB*(XS*B2/(T*X2**2)-(CX*R2/T))
505            else:
506                S = BB*B2*(XS-Orb['SEdge']*X2)/(R2*X2**2-X2*B2)
507            A = ALG[i]
508            D += A*S
509        return D
510   
511    AU = 2.80022e+7
512    C1 = 0.02721
513    C = 137.0367
514    FP = 0.0
515    FPP = 0.0
516    Mu = 0.0
517    LKev = math.log(KEv)
518    RX = KEv/C1
519    if Orbs:
520        for Orb in Orbs:
521            CX = 0.0
522            BB = Orb['BB']
523            BindEn = Orb['BindEn']
524            if Orb['IfBe'] != 0: ElEterm = Orb['ElEterm']
525            if BindEn <= KEv:
526                CX = math.exp(Aitken(Orb,LKev))
527                Mu += CX
528                CX /= AU
529            Corr = 0.0
530            if Orb['IfBe'] == 0 and BindEn >= KEv:
531                CX = 0.0
532                FPI = DGauss(Orb,CX,RX,3)
533                Corr = 0.5*Orb['SEdge']*BB**2*math.log((RX-BB)/(-RX-BB))/RX
534            else:
535                FPI = DGauss(Orb,CX,RX,Orb['IfBe'])
536                if CX != 0.0: Corr = -0.5*CX*RX*math.log((RX+BB)/(RX-BB))
537            FPI = (FPI+Corr)*C/(2.0*math.pi**2)
538            FPPI = C*CX*RX/(4.0*math.pi)
539            FP += FPI
540            FPP += FPPI
541        FP -= ElEterm
542   
543    return (FP, FPP, Mu)
544   
545mapDefault = {'MapType':'','RefList':'','GridStep':0.25,'Show bonds':True,
546                'rho':[],'rhoMax':0.,'mapSize':10.0,'cutOff':50.,'Flip':False}
547
548def SetupGeneral(data, dirname):
549    '''Initialize the General sections of the Phase tree contents
550    Called by SetupGeneral in GSASIIphsGUI and in GSASIIscriptable.SetupGeneral
551    '''
552    generalData = data['General']
553    atomData = data['Atoms']
554    generalData['AtomTypes'] = []
555    generalData['Isotopes'] = {}
556# various patches
557    if 'Isotope' not in generalData:
558        generalData['Isotope'] = {}
559    if 'Data plot type' not in generalData:
560        generalData['Data plot type'] = 'Mustrain'
561    if 'POhkl' not in generalData:
562        generalData['POhkl'] = [0,0,1]
563    if 'Map' not in generalData:
564        generalData['Map'] = mapDefault.copy()
565    if 'Flip' not in generalData:
566        generalData['Flip'] = {'RefList':'','GridStep':0.25,'Norm element':'None',
567            'k-factor':0.1,'k-Max':20.,}
568    if 'testHKL' not in generalData['Flip']:
569        generalData['Flip']['testHKL'] = [[0,0,2],[2,0,0],[1,1,1],[0,2,0],[1,2,3]]
570    if 'doPawley' not in generalData:
571        generalData['doPawley'] = False     #ToDo: change to ''
572    if 'Pawley dmin' not in generalData:
573        generalData['Pawley dmin'] = 1.0
574    if 'Pawley dmax' not in generalData:
575        generalData['Pawley dmax'] = 100.0
576    if 'Pawley neg wt' not in generalData:
577        generalData['Pawley neg wt'] = 0.0
578    if '3Dproj' not in generalData:
579        generalData['3Dproj'] = ''
580    if 'doDysnomia' not in generalData:
581        generalData['doDysnomia'] = False
582    if 'Algolrithm' in generalData.get('MCSA controls',{}) or \
583        'MCSA controls' not in generalData:
584        generalData['MCSA controls'] = {'Data source':'','Annealing':[0.7,0.1,250],
585        'dmin':2.8,'Algorithm':'log','fast parms':[0.8,0.6],'log slope':0.9,
586        'Cycles':1,'Results':[],'newDmin':True}
587    if 'AtomPtrs' not in generalData:
588        generalData['AtomPtrs'] = [3,1,7,9]
589        if generalData['Type'] == 'macromolecular':
590            generalData['AtomPtrs'] = [6,4,10,12]
591        elif generalData['Type'] == 'magnetic':
592            generalData['AtomPtrs'] = [3,1,10,12]
593    if generalData['Modulated']:
594        if 'Super' not in generalData:
595            generalData['Super'] = 1
596            generalData['SuperVec'] = [[0.,0.,0.],False,1]
597            generalData['SSGData'] = {}
598        if '4DmapData' not in generalData:
599            generalData['4DmapData'] = mapDefault.copy()
600            generalData['4DmapData'].update({'MapType':'Fobs'})
601        atomData = data['Atoms']
602        for atom in atomData:
603#                if 'SS1' not in atom:
604#                    atom += [[],[],{'SS1':{'waveType':'Fourier','Sfrac':[],'Spos':[],'Sadp':[],'Smag':[]}}]
605            if isinstance(atom[-1],dict) and 'waveType' in atom[-1]['SS1']:
606                waveType = atom[-1]['SS1']['waveType']
607                for parm in ['Sfrac','Spos','Sadp','Smag']:
608                    if len(atom[-1]['SS1'][parm]):
609                        wType = 'Fourier'
610                        if parm == 'Sfrac':
611                            if 'Crenel' in waveType:
612                                wType = 'Crenel'
613                        elif parm == 'Spos':
614                            if not 'Crenel' in waveType:
615                                wType = waveType
616                        atom[-1]['SS1'][parm] = [wType,]+list(atom[-1]['SS1'][parm])
617                del atom[-1]['SS1']['waveType']
618    else:
619        generalData['Super'] = 0
620    if 'Modulated' not in generalData:
621        generalData['Modulated'] = False
622    if 'HydIds' not in generalData:
623        generalData['HydIds'] = {}
624    if generalData['Type'] == 'magnetic':
625        if 'SGGray' not in generalData['SGData']:
626            generalData['SGData']['SGGray'] = False
627    if 'Resolution' in generalData['Map']:
628        generalData['Map']['GridStep'] = generalData['Map']['Resolution']/2.
629        generalData['Flip']['GridStep'] = generalData['Flip']['Resolution']/2.
630        del generalData['Map']['Resolution']
631        del generalData['Flip']['Resolution']
632    if 'Compare' not in generalData:
633        generalData['Compare'] = {'Oatoms':'','Tatoms':'',
634                'Tilts':{'Otilts':[],'Ttilts':[]},
635            'Bonds':{'Obonds':[],'Tbonds':[]},'Vects':{'Ovec':[],'Tvec':[]},
636            'dVects':{'Ovec':[],'Tvec':[]},'Sampling':1.0}
637    if 'Sampling' not in generalData['Compare']:
638        generalData['Compare']['Sampling'] = 1.0
639
640# end of patches
641    cx,ct,cs,cia = generalData['AtomPtrs']
642    generalData['NoAtoms'] = {}
643    generalData['BondRadii'] = []
644    generalData['AngleRadii'] = []
645    generalData['vdWRadii'] = []
646    generalData['AtomMass'] = []
647    generalData['Color'] = []
648    if generalData['Type'] == 'magnetic':
649        generalData['MagDmin'] = generalData.get('MagDmin',1.0)
650        landeg = generalData.get('Lande g',[])
651    generalData['Mydir'] = dirname
652    badList = {}
653    for iat,atom in enumerate(atomData):
654        atom[ct] = atom[ct].lower().capitalize()              #force to standard form
655        if generalData['AtomTypes'].count(atom[ct]):
656            generalData['NoAtoms'][atom[ct]] += atom[cx+3]*float(atom[cs+1])
657        elif atom[ct] != 'UNK':
658            Info = GetAtomInfo(atom[ct])
659            if not Info:
660                if atom[ct] not in badList:
661                    badList[atom[ct]] = 0
662                badList[atom[ct]] += 1
663                atom[ct] = 'UNK'
664                continue
665            atom[ct] = Info['Symbol'] # N.B. symbol might be changed by GetAtomInfo
666            generalData['AtomTypes'].append(atom[ct])
667            generalData['Z'] = Info['Z']
668            generalData['Isotopes'][atom[ct]] = Info['Isotopes']
669            generalData['BondRadii'].append(Info['Drad'])
670            generalData['AngleRadii'].append(Info['Arad'])
671            generalData['vdWRadii'].append(Info['Vdrad'])
672            if atom[ct] in generalData['Isotope']:
673                if generalData['Isotope'][atom[ct]] not in generalData['Isotopes'][atom[ct]]:
674                    isotope = list(generalData['Isotopes'][atom[ct]].keys())[-1]
675                    generalData['Isotope'][atom[ct]] = isotope
676                generalData['AtomMass'].append(Info['Isotopes'][generalData['Isotope'][atom[ct]]]['Mass'])
677            else:
678                generalData['Isotope'][atom[ct]] = 'Nat. Abund.'
679                if 'Nat. Abund.' not in generalData['Isotopes'][atom[ct]]:
680                    isotope = list(generalData['Isotopes'][atom[ct]].keys())[-1]
681                    generalData['Isotope'][atom[ct]] = isotope
682                generalData['AtomMass'].append(Info['Mass'])
683            generalData['NoAtoms'][atom[ct]] = atom[cx+3]*float(atom[cs+1])
684            generalData['Color'].append(Info['Color'])
685            if generalData['Type'] == 'magnetic':
686                if len(landeg) < len(generalData['AtomTypes']):
687                    landeg.append(2.0)
688    if generalData['Type'] == 'magnetic':
689        generalData['Lande g'] = landeg[:len(generalData['AtomTypes'])]
690
691    F000X = 0.
692    F000N = 0.
693    for i,elem in enumerate(generalData['AtomTypes']):
694        F000X += generalData['NoAtoms'][elem]*generalData['Z']
695        isotope = generalData['Isotope'][elem]
696        F000N += generalData['NoAtoms'][elem]*generalData['Isotopes'][elem][isotope]['SL'][0]
697    generalData['F000X'] = F000X
698    generalData['F000N'] = F000N
699    generalData['Mass'] = G2mth.getMass(generalData)
700 
701    if badList:
702        msg = 'Warning: element symbol(s) not found:'
703        for key in badList:
704            msg += '\n\t' + key
705            if badList[key] > 1:
706                msg += ' (' + str(badList[key]) + ' times)'
707        #wx.MessageBox(msg,caption='Element symbol error')
708        raise ValueError("Phase error:\n" + msg)       
Note: See TracBrowser for help on using the repository browser.