source: trunk/GSASIIElem.py

Last change on this file was 5572, checked in by toby, 5 months ago

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