Changeset 942
- Timestamp:
- Jun 5, 2013 9:26:36 AM (10 years ago)
- Location:
- trunk
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIElem.py
r939 r942 59 59 return FormFactors 60 60 61 def GetFFtable(atomTypes): 62 ''' returns a dictionary of form factor data for atom types found in atomTypes 63 64 :param list atomTypes: list of atom types 65 :return: FFtable, dictionary of form factor data; key is atom type 66 67 ''' 68 FFtable = {} 69 for El in atomTypes: 70 FFs = GetFormFactorCoeff(El.split('+')[0].split('-')[0]) 71 for item in FFs: 72 if item['Symbol'] == El.upper(): 73 FFtable[El] = item 74 return FFtable 75 76 def GetBLtable(General): 77 ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General 78 79 :param dict General: dictionary of phase info.; includes AtomTypes & Isotopes 80 :return: BLtable, dictionary of scattering length data; key is atom type 81 ''' 82 atomTypes = General['AtomTypes'] 83 BLtable = {} 84 isotopes = General['Isotopes'] 85 isotope = General['Isotope'] 86 for El in atomTypes: 87 BLtable[El] = [isotope[El],isotopes[El][isotope[El]]] 88 return BLtable 89 90 def getFFvalues(FFtables,SQ): 91 'Needs a doc string' 92 FFvals = {} 93 for El in FFtables: 94 FFvals[El] = ScatFac(FFtables[El],SQ)[0] 95 return FFvals 96 97 def getBLvalues(BLtables): 98 'Needs a doc string' 99 BLvals = {} 100 for El in BLtables: 101 BLvals[El] = BLtables[El][1][1] 102 return BLvals 103 61 104 def GetFFC5(ElSym): 62 105 '''Get 5 term form factor and Compton scattering data -
trunk/GSASIIIO.py
r941 r942 57 57 inst[item] = list(inst[item]) 58 58 return inst 59 60 59 61 60 def FileDlgFixExt(dlg,file): -
trunk/GSASIImath.py
r940 r942 1900 1900 1901 1901 1902 def mcsaSearch(data, reflType,reflData,covData,pgbar):1902 def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar): 1903 1903 gamFW = lambda s,g: math.exp(math.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.) 1904 1904 … … 1912 1912 upper.append(limits[1]) 1913 1913 1914 def getAtomparms(item,pfx, parmDict,varyList):1914 def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList): 1915 1915 parmDict[pfx+'Atype'] = item['atType'] 1916 pstr = ['x','y','z'] 1916 aTypes |= set([item['atType'],]) 1917 pstr = ['Ax','Ay','Az'] 1918 XYZ = [0,0,0] 1917 1919 for i in range(3): 1918 1920 name = pfx+pstr[i] 1919 1921 parmDict[name] = item['Pos'][0][i] 1922 XYZ[i] = parmDict[name] 1920 1923 if item['Pos'][1][i]: 1921 1924 varyList += [name,] … … 1923 1926 lower.append(limits[0]) 1924 1927 upper.append(limits[1]) 1928 parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData)) 1925 1929 1926 def getRBparms(item, pfx,parmDict,varyList):1927 parmDict[ 'MolCent'] = item['MolCent']1928 parmDict[ 'RBId'] = item['RBId']1930 def getRBparms(item,mfx,aTypes,RBdata,parmDict,varyList): 1931 parmDict[mfx+'MolCent'] = item['MolCent'] 1932 parmDict[mfx+'RBId'] = item['RBId'] 1929 1933 pstr = ['Px','Py','Pz'] 1930 1934 ostr = ['Qa','Qi','Qj','Qk'] 1931 1935 for i in range(3): 1932 name = pfx+pstr[i]1936 name = mfx+pstr[i] 1933 1937 parmDict[name] = item['Pos'][0][i] 1934 1938 if item['Pos'][1][i]: … … 1938 1942 upper.append(limits[1]) 1939 1943 for i in range(4): 1940 name = pfx+ostr[i]1944 name = mfx+ostr[i] 1941 1945 parmDict[name] = item['Ori'][0][i] 1942 1946 if item['Ovar'] == 'AV' and i: … … 1952 1956 if 'Tor' in item: #'Tor' not there for 'Vector' RBs 1953 1957 for i in range(len(item['Tor'][0])): 1954 name = pfx+'Tor'+str(i)1958 name = mfx+'Tor'+str(i) 1955 1959 parmDict[name] = item['Tor'][0][i] 1956 1960 if item['Tor'][1][i]: … … 1959 1963 lower.append(limits[0]) 1960 1964 upper.append(limits[1]) 1961 1962 def getFFvalues(FFtables,SQ): 1963 FFvals = {} 1964 for El in FFtables: 1965 FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0] 1966 return FFvals 1967 1968 def getBLvalues(BLtables): 1969 BLvals = {} 1970 for El in BLtables: 1971 BLvals[El] = BLtables[El][1][1] 1972 return BLvals 1973 1974 # def mcsaSfCalc(refList,G,SGData,calcControls,parmDict): 1965 aTypes |= set(RBdata[item['Type']][item['RBId']]['rbTypes']) 1966 1967 def GetAtomFX(pfx,Natoms,parmDict): 1968 'Needs a doc string' 1969 Tdata = Natoms*[' ',] 1970 Mdata = np.zeros(Natoms) 1971 Fdata = np.zeros(Natoms) 1972 Xdata = np.zeros((3,Natoms)) 1973 keys = {'Atype:':Tdata,'Amul:':Mdata, 1974 'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2]} 1975 for iatm in range(Natoms): 1976 for key in keys: 1977 parm = pfx+key+str(iatm) 1978 if parm in parmDict: 1979 keys[key][iatm] = parmDict[parm] 1980 return Tdata,Mdata,Xdata 1981 1982 # def mcsaSfCalc(refList,G,SGData,parmDict): 1975 1983 # ''' Compute structure factors for all h,k,l for phase 1976 1984 # input: … … 1978 1986 # G: reciprocal metric tensor 1979 1987 # SGData: space group info. dictionary output from SpcGroup 1980 # calcControls:1981 1988 # ParmDict: 1982 1989 # puts result F^2 in each ref[8] in refList … … 1986 1993 # ast = np.sqrt(np.diag(G)) 1987 1994 # Mast = twopisq*np.multiply.outer(ast,ast) 1988 # FFtables = calcControls['FFtables'] 1989 # BLtables = calcControls['BLtables'] 1990 # Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict) 1995 # Tdata,Mdata,Fdata,Xdata = GetAtomFX(pfx,calcControls,parmDict) 1991 1996 # FF = np.zeros(len(Tdata)) 1992 1997 # if 'N' in parmDict[hfx+'Type']: … … 1995 2000 # FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata]) 1996 2001 # FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata]) 1997 # maxPos = len(SGData['SGOps'])1998 2002 # Uij = np.array(G2lat.U6toUij(Uijdata)) 1999 2003 # bij = Mast*Uij.T … … 2013 2017 # Uniq = refl[11] 2014 2018 # phi = refl[12] 2015 # phase = twopi*(np.inner(Uniq,( dXdata.T+Xdata.T))+phi[:,np.newaxis])2019 # phase = twopi*(np.inner(Uniq,(Xdata.T))+phi[:,np.newaxis]) 2016 2020 # sinp = np.sin(phase) 2017 2021 # cosp = np.cos(phase) 2018 2022 # occ = Mdata*Fdata/len(Uniq) 2019 # biso = -SQfactor*Uisodata 2020 # Tiso = np.where(biso<1.,np.exp(biso),1.0) 2021 # HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq]) 2022 # Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 2023 # Tcorr = Tiso*Tuij 2024 # fa = np.array([(FF+FP-Bab)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr]) 2023 # fa = np.array([(FF+FP-Bab)*occ*cosp,-FPP*occ*sinp]) 2025 2024 # fas = np.sum(np.sum(fa,axis=1),axis=1) #real 2026 2025 # if not SGData['SGInv']: … … 2036 2035 sq4pi = np.sqrt(4*np.pi) 2037 2036 generalData = data['General'] 2037 SGData = generalData['SGData'] 2038 2038 pId = data['pId'] 2039 pfx = str(pId)+'::' 2040 Atoms = data['Atoms'] #if any 2039 fixAtoms = data['Atoms'] #if any 2040 cx,ct,cs = generalData['AtomPtrs'][:3] 2041 aTypes = set([]) 2042 parmDict = {} 2043 varyList = [] 2044 atNo = 0 2045 for atm in fixAtoms: 2046 pfx = ':'+str(atNo)+':' 2047 parmDict[pfx+'Atype'] = atm[ct] 2048 aTypes |= set([atm[ct],]) 2049 pstr = ['Ax','Ay','Az'] 2050 parmDict[pfx+'Amul'] = atm[cs+1] 2051 for i in range(3): 2052 name = pfx+pstr[i] 2053 parmDict[name] = atm[cx+i] 2054 atNo += 1 2041 2055 mcsaControls = generalData['MCSA controls'] 2042 2056 reflName = mcsaControls['Data source'] 2043 2057 phaseName = generalData['Name'] 2044 2058 MCSAObjs = data['MCSA']['Models'] #list of MCSA models 2045 parmDict = {}2046 varyList = []2047 2059 upper = [] 2048 2060 lower = [] … … 2052 2064 getMDparms(item,mfx,parmDict,varyList) 2053 2065 elif item['Type'] == 'Atom': 2054 getAtomparms(item,mfx,parmDict,varyList) 2066 pfx = mfx+str(atNo)+':' 2067 getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList) 2068 atNo += 1 2055 2069 elif item['Type'] in ['Residue','Vector']: 2056 getRBparms(item,mfx,parmDict,varyList) 2070 pfx = mfx+':' 2071 getRBparms(item,pfx,aTypes,RBdata,parmDict,varyList) 2072 FFtables = G2el.GetFFtable(aTypes) 2073 refs = [] 2057 2074 if 'PWDR' in reflName: 2058 refs = []2059 2075 for ref in reflData: 2060 2076 h,k,l,m,d,pos,sig,gam,Fsq = ref[:9] 2061 2077 if d >= mcsaControls['dmin']: 2062 FWHM = gamFW(sig,gam)/2.35482 #sqrt(8ln2) --> sig from FWHM 2063 Fsq *= m 2078 sig = gamFW(sig,gam)/sq8ln2 #--> sig from FWHM 2079 SQ = 0.25/d**2 2080 Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:] 2081 FFs = G2el.getFFvalues(FFtables,SQ) 2082 refs.append([h,k,l,m,f*m,pos,sig,FFs,Uniq,phi]) 2083 nRef = len(refs) 2084 rcov = np.zeros((nRef,nRef)) 2085 for iref,refI in enumerate(refs): 2086 rcov[iref][iref] = 1./(sq4pi*refI[6]) 2087 for jref,refJ in enumerate(refs[:iref]): 2088 t1 = refI[6]**2+refJ[6]**2 2089 t2 = (refJ[5]-refI[5])**2/(2.*t1) 2090 if t2 > 10.: 2091 rcov[iref][jref] = 0. 2092 else: 2093 rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2)) 2094 rcov += (rcov.T-np.diagflat(np.diagonal(rcov))) 2095 Rdiag = np.sqrt(np.diag(rcov)) 2096 Rnorm = np.outer(Rdiag,Rdiag) 2097 rcov /= Rnorm 2064 2098 elif 'Pawley' in reflName: 2065 refs = []2066 2099 covMatrix = covData['covMatrix'] 2067 v aryList = covData['varyList']2100 vList = covData['varyList'] 2068 2101 for iref,refI in enumerate(reflData): 2069 2102 h,k,l,m,d,v,f,s = refI 2070 2103 if d >= mcsaControls['dmin'] and v: #skip unrefined ones 2071 nameI = pfx+'PWLref:'+str(iref) 2072 if nameI in covData['varyList']: 2073 refs.append([h,k,l,f*m,0.,0.]) 2104 SQ = 0.25/d**2 2105 Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:] 2106 FFs = G2el.getFFvalues(FFtables,SQ) 2107 refs.append([h,k,l,m,f*m,0.,0.,FFs,Uniq,phi]) 2074 2108 nRef = len(refs) 2075 print nRef 2076 covTerms = np.zeros((nRef,nRef)) 2077 print covTerms.shape 2078 2109 rcov = np.zeros((nRef,nRef)) 2079 2110 Iref = 0 2080 2111 for iref,refI in enumerate(reflData): … … 2082 2113 nameI = pfx+'PWLref:'+str(iref) 2083 2114 if nameI in covData['varyList']: 2084 Iindx = v aryList.index(nameI)2085 covTerms[Iref][Iref] = covMatrix[Iindx][Iindx]2115 Iindx = vList.index(nameI) 2116 rcov[Iref][Iref] = covMatrix[Iindx][Iindx] 2086 2117 Jref = 0 2087 2118 for jref,refJ in enumerate(reflData[:iref]): … … 2089 2120 nameJ = pfx+'PWLref:'+str(jref) 2090 2121 try: 2091 covTerms[Iref][Jref] = covMatrix[varyList.index(nameI)][varyList.index(nameJ)]2122 rcov[Iref][Jref] = covMatrix[vList.index(nameI)][vList.index(nameJ)] 2092 2123 except ValueError: 2093 covTerms[Iref][Jref] = covTerms[Iref][Jref-1]2124 rcov[Iref][Jref] = rcov[Iref][Jref-1] 2094 2125 Jref += 1 2126 else: 2127 rcov[Iref] = rcov[Iref-1] 2128 rcov[Iref][Iref] = rcov[Iref-1][Iref-1] 2095 2129 Iref += 1 2096 print covTerms 2130 rcov += (rcov.T-np.diagflat(np.diagonal(rcov))) 2131 Rdiag = np.sqrt(np.diag(rcov)) 2132 Rnorm = np.outer(Rdiag,Rdiag) 2133 rcov /= Rnorm 2134 elif 'HKLF' in reflName: 2135 for ref in reflData: 2136 [h,k,l,m,d],f = ref[:5],ref[6] 2137 if d >= mcsaControls['dmin']: 2138 SQ = 0.25/d**2 2139 Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:] 2140 FFs = G2el.getFFvalues(FFtables,SQ) 2141 refs.append([h,k,l,m,f*m,0.,0.,FFs,Uniq,phi]) 2142 rcov = np.identity(len(refs)) 2143 2144 for parm in parmDict: 2145 print parm,parmDict[parm] 2097 2146 2098 2099 # covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals, 2100 # 'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict} 2101 elif 'HKLF' in reflName: 2102 refs = [] 2103 for ref in reflData: 2104 h,k,l,m,d,Fsq = ref[:5],ref[6] 2105 if d >= mcsaControls['dmin']: 2106 Fsq *= m 2107 refs.append([h,k,l,m,Fsq]) 2108 rcov = np.identity(len(refs)) 2109 2110 2147 # XYZ,aTypes = UpdateMCSAxyz(Bmat,MCSA) 2111 2148 2112 2149 # generalData['MCSA controls'] = {'Data source':'','Annealing':[50.,0.001,50,1.e-6], -
trunk/GSASIIphsGUI.py
r941 r942 4370 4370 print '**** ERROR - No data defined for MC/SA run' 4371 4371 return 4372 RBdata = G2frame.PatternTree.GetItemPyData( 4373 G2gd.GetPatternTreeItemId(G2frame,G2frame.root,'Rigid bodies')) 4372 4374 MCSAmodels = MCSAdata['Models'] 4373 4375 if not len(MCSAmodels): … … 4382 4384 pgbar.SetSize(Size) 4383 4385 try: 4384 MCSAdata['Results'] = G2mth.mcsaSearch(data, reflType,reflData,covData,pgbar)4386 MCSAdata['Results'] = G2mth.mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar) 4385 4387 finally: 4386 4388 pgbar.Destroy() … … 4941 4943 page = event.GetSelection() 4942 4944 text = G2frame.dataDisplay.GetPageText(page) 4945 Atoms.ClearGrid() 4943 4946 if text == 'Atoms': 4944 4947 G2gd.SetDataMenuBar(G2frame,G2frame.dataFrame.AtomsMenu) -
trunk/GSASIIstrIO.py
r939 r942 483 483 print >>pFile,' Initial shift factor: ','%.3f'%(Controls['shift factor']) 484 484 485 def GetFFtable(General):486 ''' returns a dictionary of form factor data for atom types found in General487 488 :param dict General: dictionary of phase info.; includes AtomTypes489 :return: FFtable, dictionary of form factor data; key is atom type490 491 '''492 atomTypes = General['AtomTypes']493 FFtable = {}494 for El in atomTypes:495 FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])496 for item in FFs:497 if item['Symbol'] == El.upper():498 FFtable[El] = item499 return FFtable500 501 def GetBLtable(General):502 ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General503 504 :param dict General: dictionary of phase info.; includes AtomTypes & Isotopes505 :return: BLtable, dictionary of scattering length data; key is atom type506 '''507 atomTypes = General['AtomTypes']508 BLtable = {}509 isotopes = General['Isotopes']510 isotope = General['Isotope']511 for El in atomTypes:512 BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]513 return BLtable514 515 485 def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary): 516 486 'needs a doc string' … … 900 870 pId = PhaseData[name]['pId'] 901 871 pfx = str(pId)+'::' 902 FFtable = G etFFtable(General)903 BLtable = G etBLtable(General)872 FFtable = G2el.GetFFtable(General['AtomTypes']) 873 BLtable = G2el.GetBLtable(General) 904 874 FFtables.update(FFtable) 905 875 BLtables.update(BLtable) -
trunk/GSASIIstrMath.py
r939 r942 516 516 return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata 517 517 518 def getFFvalues(FFtables,SQ):519 'Needs a doc string'520 FFvals = {}521 for El in FFtables:522 FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0]523 return FFvals524 525 def getBLvalues(BLtables):526 'Needs a doc string'527 BLvals = {}528 for El in BLtables:529 BLvals[El] = BLtables[El][1][1]530 return BLvals531 532 518 def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict): 533 519 ''' Compute structure factors for all h,k,l for phase … … 557 543 FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata]) 558 544 FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata]) 559 maxPos = len(SGData['SGOps'])560 545 Uij = np.array(G2lat.U6toUij(Uijdata)) 561 546 bij = Mast*Uij.T … … 568 553 if not len(refl[-1]): #no form factors 569 554 if 'N' in parmDict[hfx+'Type']: 570 refl[-1] = getBLvalues(BLtables)555 refl[-1] = G2el.getBLvalues(BLtables) 571 556 else: #'X' 572 refl[-1] = getFFvalues(FFtables,SQ)557 refl[-1] = G2el.getFFvalues(FFtables,SQ) 573 558 for i,El in enumerate(Tdata): 574 559 FF[i] = refl[-1][El] … … 611 596 FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata]) 612 597 FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata]) 613 maxPos = len(SGData['SGOps'])614 598 Uij = np.array(G2lat.U6toUij(Uijdata)) 615 599 bij = Mast*Uij.T
Note: See TracChangeset
for help on using the changeset viewer.