- Timestamp:
- Oct 15, 2013 12:25:24 PM (12 years ago)
- Location:
- trunk
- Files:
-
- 11 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified trunk/GSASII.py ¶
r1105 r1106 557 557 Id = self.PatternTree.AppendItem(parent=self.root, 558 558 text='HKLF '+HistName) 559 self.PatternTree.SetItemPyData(Id,[{'wtFactor':1.0,'Dummy':False},rd.Ref List])559 self.PatternTree.SetItemPyData(Id,[{'wtFactor':1.0,'Dummy':False},rd.RefDict]) 560 560 Sub = self.PatternTree.AppendItem(Id,text='Instrument Parameters') 561 561 self.PatternTree.SetItemPyData(Sub,rd.Parameters) … … 2183 2183 for item in self.Refine: item.Enable(True) 2184 2184 for item in self.SeqRefine: item.Enable(True) 2185 item, cookie = self.PatternTree.GetNextChild(self.root, cookie) 2185 item, cookie = self.PatternTree.GetNextChild(self.root, cookie) 2186 2186 if Id: 2187 2187 self.PatternTree.SelectItem(Id) -
TabularUnified trunk/GSASIIIO.py ¶
r1102 r1106 1582 1582 self.InitParameters() 1583 1583 self.InitControls() 1584 self.Ref List = []1584 self.RefDict = {'RefList':[],'Uniq':[],'Phi':[],'FF':[]} 1585 1585 1586 1586 def InitControls(self): … … 1621 1621 HKLmin = [None,None,None] 1622 1622 Fo2max = None 1623 for refl in self.Ref List:1623 for refl in self.RefDict['RefList']: 1624 1624 HKL = refl[:3] 1625 1625 if Fo2max is None: … … 1853 1853 if hnum is None: return True 1854 1854 self.histnam = [choices[hnum]] 1855 1855 1856 def loadParmDict(self): 1856 1857 '''Load the GSAS-II refinable parameters from the tree into a dict (self.parmDict). Update -
TabularUnified trunk/GSASIIgrid.py ¶
r1087 r1106 3069 3069 newdata = np.linspace(start,end,N,True) 3070 3070 if len(newdata) < 2: return # too small a range - reject 3071 data[1][0] = newdata 3072 data[1][1] = np.zeros_like(newdata) 3073 data[1][2] = np.ones_like(newdata) 3074 data[1][3] = np.zeros_like(newdata) 3075 data[1][4] = np.zeros_like(newdata) 3076 data[1][5] = np.zeros_like(newdata) 3071 data[1] = [newdata,np.zeros_like(newdata),np.ones_like(newdata), 3072 np.zeros_like(newdata),np.zeros_like(newdata),np.zeros_like(newdata)] 3077 3073 Tmin = newdata[0] 3078 3074 Tmax = newdata[-1] … … 3093 3089 3094 3090 data = G2frame.PatternTree.GetItemPyData(item) 3091 #patches 3095 3092 if 'wtFactor' not in data[0]: 3096 3093 data[0] = {'wtFactor':1.0} 3094 if isinstance(data[1],list) and kind == 'HKLF': 3095 RefData = {'RefList':[],'Uniq':[],'Phi':[],'FF':[]} 3096 for ref in data[1]: 3097 RefData['RefList'].append(ref[:11]+[ref[13],]) 3098 RefData['Uniq'].append(ref[11]) 3099 RefData['Phi'].append(ref[12]) 3100 RefData['FF'].append(ref[14]) 3101 data[1] = RefData 3102 G2frame.PatternTree.SetItemPyData(item,data) 3103 #end patches 3097 3104 if G2frame.dataDisplay: 3098 3105 G2frame.dataDisplay.Destroy() -
TabularUnified trunk/GSASIImath.py ¶
r1097 r1106 104 104 if Print: 105 105 print ' Hessian refinement on %d variables:'%(n) 106 Lam = np.zeros((n,n)) 106 107 while icycle < maxcyc: 107 108 time0 = time.time() … … 145 146 break 146 147 icycle += 1 147 M = func(x0,*args) 148 nfev += 1 149 Yvec,Amat = Hess(x0,*args) 150 Amatlam = Amat*(One+Lam)/Anorm #scale Amat to Marquardt array 148 else: #after last cycle or if zero cycles 149 M = func(x0,*args) 150 nfev += 1 151 Yvec,Amat = Hess(x0,*args) 152 Adiag = np.sqrt(np.diag(Amat)) 153 Anorm = np.outer(Adiag,Adiag) 154 Amat /= Anorm 155 Amatlam = Amat*(One+Lam)/Anorm #scale Amat to Marquardt array 151 156 try: 152 157 Bmat = nl.inv(Amatlam)*(One+Lam)/Anorm #rescale Bmat to Marquardt array … … 1337 1342 Hmax[2] = ((Hmax[2]+1)/4)*4 1338 1343 1339 def OmitMap(data,reflD ata):1344 def OmitMap(data,reflDict): 1340 1345 '''default doc string 1341 1346 … … 1358 1363 Fhkl = np.zeros(shape=2*Hmax,dtype='c16') 1359 1364 time0 = time.time() 1360 for ref in reflData:1365 for iref,ref in enumerate(reflDict['RefList']): 1361 1366 if ref[4] >= dmin: 1362 1367 Fosq,Fcsq,ph = ref[8:11] 1363 for i,hkl in enumerate(ref [11]):1368 for i,hkl in enumerate(reflDict['Uniq']): #uses uniq 1364 1369 hkl = np.asarray(hkl,dtype='i') 1365 dp = 360.*ref [12][i]1370 dp = 360.*reflDict['Phi'][i] #and phi 1366 1371 a = cosd(ph+dp) 1367 1372 b = sind(ph+dp) … … 1381 1386 return mapData 1382 1387 1383 def FourierMap(data,reflD ata):1388 def FourierMap(data,reflDict): 1384 1389 '''default doc string 1385 1390 … … 1403 1408 # Fhkl[0,0,0] = generalData['F000X'] 1404 1409 time0 = time.time() 1405 for ref in reflData:1410 for iref,ref in enumerate(reflDict['RefList']): 1406 1411 if ref[4] >= dmin: 1407 1412 Fosq,Fcsq,ph = ref[8:11] 1408 for i,hkl in enumerate(ref [11]):1413 for i,hkl in enumerate(reflDict['Uniq'][iref]): #uses uniq 1409 1414 hkl = np.asarray(hkl,dtype='i') 1410 dp = 360.*ref [12][i]1415 dp = 360.*reflDict['Phi'][iref][i] #and phi 1411 1416 a = cosd(ph+dp) 1412 1417 b = sind(ph+dp) … … 1545 1550 return DX 1546 1551 1547 def ChargeFlip(data,reflD ata,pgbar):1552 def ChargeFlip(data,reflDict,pgbar): 1548 1553 '''default doc string 1549 1554 … … 1572 1577 Ehkl = np.zeros(shape=2*Hmax,dtype='c16') #2X64bits per complex no. 1573 1578 time0 = time.time() 1574 for ref in reflData:1579 for iref,ref in enumerate(reflDict['RefList']): 1575 1580 dsp = ref[4] 1576 1581 if dsp >= dmin: … … 1585 1590 ph = ref[10] 1586 1591 ph = rn.uniform(0.,360.) 1587 for i,hkl in enumerate(ref [11]):1592 for i,hkl in enumerate(reflDict['Uniq'][iref]): #uses uniq 1588 1593 hkl = np.asarray(hkl,dtype='i') 1589 dp = 360.*ref [12][i]1594 dp = 360.*reflDict['Phi'][iref][i] #and phi 1590 1595 a = cosd(ph+dp) 1591 1596 b = sind(ph+dp) -
TabularUnified trunk/GSASIIphsGUI.py ¶
r1100 r1106 3113 3113 Cell = generalData['Cell'][1:7] 3114 3114 G,g = G2lat.cell2Gmat(Cell) 3115 for ref in reflData:3115 for iref,ref in enumerate(reflData['RefList']): 3116 3116 H = ref[:3] 3117 3117 ref[4] = np.sqrt(1./G2lat.calc_rDsq2(H,G)) 3118 iabsnt,ref[3],ref[11],ref[12] = G2spc.GenHKLf(H,SGData) 3118 iabsnt,ref[3],Uniq,phi = G2spc.GenHKLf(H,SGData) 3119 reflData['Uniq'][iref] = Uniq 3120 reflData['Phi'][iref] = phi 3119 3121 G2frame.PatternTree.SetItemPyData(Id,[refDict,reflData]) 3120 3122 … … 4416 4418 PatternId = G2gd.GetPatternTreeItemId(G2frame,G2frame.root, reflName) 4417 4419 reflSets = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists')) 4418 reflData = reflSets[phaseName] 4420 try: #patch for old reflection data 4421 reflData = reflSets[phaseName]['RefList'] 4422 except TypeError: 4423 reflData = reflSets[phaseName] 4419 4424 reflType = 'PWDR' 4420 4425 elif 'HKLF' in reflName: 4421 4426 PatternId = G2gd.GetPatternTreeItemId(G2frame,G2frame.root, reflName) 4422 reflData = G2frame.PatternTree.GetItemPyData(PatternId)[1] 4427 try: 4428 reflData = G2frame.PatternTree.GetItemPyData(PatternId)[1]['RefList'] 4429 except TypeError: 4430 reflData = G2frame.PatternTree.GetItemPyData(PatternId)[1] 4423 4431 reflType = 'HKLF' 4424 4432 elif reflName == 'Pawley reflections': … … 4883 4891 reflSets = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists')) 4884 4892 reflData = reflSets[phaseName] 4893 if isinstance(reflDict,list): #patch for old reflection data 4894 RefData = {'RefList':[],'Uniq':[],'Phi':[],'FF':[]} 4895 for ref in reflDict: 4896 RefData['RefList'].append(ref[:11]+[ref[13],]) 4897 RefData['Uniq'].append(ref[11]) 4898 RefData['Phi'].append(ref[12]) 4899 RefData['FF'].append(ref[14]) 4900 reflDict = RefData 4885 4901 elif 'HKLF' in reflName: 4886 4902 PatternId = G2gd.GetPatternTreeItemId(G2frame,G2frame.root, reflName) … … 4967 4983 PatternId = G2gd.GetPatternTreeItemId(G2frame,G2frame.root, reflName) 4968 4984 reflSets = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists')) 4969 reflData = reflSets[phaseName] 4985 reflDict = reflSets[phaseName] 4986 if isinstance(reflDict,list): #patch for old reflection data 4987 RefData = {'RefList':[],'Uniq':[],'Phi':[],'FF':[]} 4988 for ref in reflDict: 4989 RefData['RefList'].append(ref[:11]+[ref[13],]) 4990 RefData['Uniq'].append(ref[11]) 4991 RefData['Phi'].append(ref[12]) 4992 RefData['FF'].append(ref[14]) 4993 reflDict = RefData 4970 4994 elif 'HKLF' in reflName: 4971 4995 PatternId = G2gd.GetPatternTreeItemId(G2frame,G2frame.root, reflName) 4972 reflD ata= G2frame.PatternTree.GetItemPyData(PatternId)[1]4996 reflDict = G2frame.PatternTree.GetItemPyData(PatternId)[1] 4973 4997 else: 4974 4998 print '**** ERROR - No data defined for charge flipping' … … 4982 5006 pgbar.SetSize(Size) 4983 5007 try: 4984 mapData.update(G2mth.ChargeFlip(data,reflD ata,pgbar))5008 mapData.update(G2mth.ChargeFlip(data,reflDict,pgbar)) 4985 5009 finally: 4986 5010 pgbar.Destroy() -
TabularUnified trunk/GSASIIplot.py ¶
r1099 r1106 299 299 300 300 Plot.set_aspect(aspect='equal') 301 HKLref = self.PatternTree.GetItemPyData(self.Sngl)[1] 301 HKLref = self.PatternTree.GetItemPyData(self.Sngl)[1]['RefList'] 302 302 Data = self.PatternTree.GetItemPyData( 303 303 G2gd.GetPatternTreeItemId(self,self.Sngl, 'HKL Plot Controls')) … … 715 715 HKL = [] 716 716 if Phases: 717 for peak in Phases[G2frame.RefList]: 718 HKL.append(peak[:6]) 717 try: 718 for peak in Phases[G2frame.RefList]['RefList']: 719 HKL.append(peak[:6]) 720 except TypeError: 721 for peak in Phases[G2frame.RefList]: 722 HKL.append(peak[:6]) 719 723 HKL = np.array(HKL) 720 724 else: … … 875 879 Phases = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists')) 876 880 for pId,phase in enumerate(Phases): 877 peaks = Phases[phase] 881 try: #patch for old style reflection lists 882 peaks = Phases[phase]['RefList'] 883 except TypeError: 884 peaks = Phases[phase] 878 885 if not peaks: 879 886 continue -
TabularUnified trunk/GSASIIpwdGUI.py ¶
r1028 r1106 2016 2016 if HKLF: 2017 2017 G2gd.SetDataMenuBar(G2frame) 2018 ref List = [refl[:11] for refl in data[1]]2018 refs = data[1]['RefList'] 2019 2019 else: 2020 2020 G2gd.SetDataMenuBar(G2frame,G2frame.dataFrame.ReflMenu) … … 2025 2025 if len(data) > 1: 2026 2026 G2frame.dataFrame.SelectPhase.Enable(True) 2027 refList = np.array([refl[:11] for refl in data[G2frame.RefList]]) 2028 Icorr = np.array([refl[13] for refl in data[G2frame.RefList]]) 2029 I100 = refList.T[8]*Icorr 2027 try: #patch for old reflection lists 2028 refList = np.array(data[G2frame.RefList]['RefList']) 2029 I100 = refList.T[8]*refList.T[11] 2030 except TypeError: 2031 refList = np.array([refl[:11] for refl in data[G2frame.RefList]]) 2032 I100 = refList.T[8]*np.array([refl[13] for refl in data[G2frame.RefList]]) 2030 2033 Imax = np.max(I100) 2031 2034 if Imax: 2032 2035 I100 *= 100.0/Imax 2033 ref List = np.vstack((refList.T,I100)).T2036 refs = np.vstack((refList.T[:11],I100)).T 2034 2037 for i in range(len(refList)): rowLabels.append(str(i)) 2035 2038 if HKLF: … … 2040 2043 2*[wg.GRID_VALUE_FLOAT+':10,2',]+[wg.GRID_VALUE_FLOAT+':10,3',]+ \ 2041 2044 [wg.GRID_VALUE_FLOAT+':10,2',] 2042 G2frame.PeakTable = G2gd.Table(ref List,rowLabels=rowLabels,colLabels=colLabels,types=Types)2045 G2frame.PeakTable = G2gd.Table(refs,rowLabels=rowLabels,colLabels=colLabels,types=Types) 2043 2046 G2frame.dataFrame.SetLabel('Reflection List for '+phaseName) 2044 2047 G2frame.dataDisplay = G2gd.GSGrid(parent=G2frame.dataFrame) -
TabularUnified trunk/GSASIIstrIO.py ¶
r1078 r1106 295 295 HKLFdata = {} 296 296 HKLFdata.update(datum[1][0]) #weight factor 297 #patch 298 if isinstance(datum[1][1],list): 299 RefData = {'RefList':[],'Uniq':[],'Phi':[],'FF':[]} 300 for ref in datum[1][1]: 301 RefData['RefList'].append(ref[:11]+[ref[13],]) 302 RefData['Uniq'].append(ref[11]) 303 RefData['Phi'].append(ref[12]) 304 RefData['FF'].append(ref[14]) 305 datum[1][1] = RefData 306 #end patch 297 307 HKLFdata['Data'] = datum[1][1] 298 308 HKLFdata[data[1][0]] = data[1][1] #Instrument parameters … … 436 446 histogram = Histograms[datum[0]] 437 447 # print 'found ',datum[0] 438 data[0][1][1] = list(histogram['Data'])448 data[0][1][1] = histogram['Data'] 439 449 data[0][1][0].update(histogram['Residuals']) 440 450 for datus in data[1:]: … … 1743 1753 PrintBabinet(hapData['Babinet']) 1744 1754 HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A)) 1745 refList = [] 1746 for h,k,l,d in HKLd: 1747 ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData) 1748 mul *= 2 # for powder overlap of Friedel pairs 1749 if ext: 1750 continue 1751 if 'C' in inst['Type'][0]: 1752 pos = 2.0*asind(wave/(2.0*d))+Zero 1753 if limits[0] < pos < limits[1]: 1754 refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0,{}]) 1755 #last item should contain form factor values by atom type 1756 else: 1757 raise ValueError 1758 if resetRefList: Histogram['Reflection Lists'][phase] = refList 1755 if resetRefList: 1756 refList = [] 1757 Uniq = [] 1758 Phi = [] 1759 FF = [] 1760 for h,k,l,d in HKLd: 1761 ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData) 1762 mul *= 2 # for powder overlap of Friedel pairs 1763 if ext: 1764 continue 1765 if 'C' in inst['Type'][0]: 1766 pos = 2.0*asind(wave/(2.0*d))+Zero 1767 if limits[0] < pos < limits[1]: 1768 refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,0.0]) 1769 Uniq.append(uniq) 1770 Phi.append(phi) 1771 FF.append({}) 1772 else: 1773 raise ValueError 1774 Histogram['Reflection Lists'][phase] = {'RefList':refList,'Uniq':Uniq,'Phi':Phi,'FF':FF} 1759 1775 elif 'HKLF' in histogram: 1760 1776 inst = Histogram['Instrument Parameters'][0] -
TabularUnified trunk/GSASIIstrMath.py ¶
r1103 r1106 528 528 return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata 529 529 530 def StructureFactor(ref List,G,hfx,pfx,SGData,calcControls,parmDict):530 def StructureFactor(refDict,G,hfx,pfx,SGData,calcControls,parmDict): 531 531 ''' Compute structure factors for all h,k,l for phase 532 532 puts the result, F^2, in each ref[8] in refList 533 533 input: 534 534 535 :param list refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv] 535 :param dict refDict: where 536 'RefList' list where each ref = h,k,l,m,d,... 537 'Uniq' list of [equiv h,k,l] 538 'Phi' list of phase[equiv] 539 'FF' dict of form factors - filed in below 536 540 :param np.array G: reciprocal metric tensor 537 541 :param str pfx: phase id string … … 557 561 Uij = np.array(G2lat.U6toUij(Uijdata)) 558 562 bij = Mast*Uij.T 559 for refl in refList:563 for iref,refl in enumerate(refDict['RefList']): 560 564 fbs = np.array([0,0]) 561 565 H = refl[:3] … … 563 567 SQfactor = 4.0*SQ*twopisq 564 568 Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor) 565 if not len(ref l[-1]): #no form factors569 if not len(refDict['FF'][iref]): #no form factors 566 570 if 'N' in calcControls[hfx+'histType']: 567 ref l[-1] = G2el.getBLvalues(BLtables)571 refDict['FF'][iref] = G2el.getBLvalues(BLtables) 568 572 else: #'X' 569 ref l[-1] = G2el.getFFvalues(FFtables,SQ)573 refDict['FF'][iref] = G2el.getFFvalues(FFtables,SQ) 570 574 for i,El in enumerate(Tdata): 571 FF[i] = refl[-1][El] 572 Uniq = refl[11] 573 phi = refl[12] 574 phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis]) 575 FF[i] = refDict['FF'][iref][El] 576 phase = twopi*(np.inner(refDict['Uniq'][iref],(dXdata.T+Xdata.T))+refDict['Phi'][iref][:,np.newaxis]) 575 577 sinp = np.sin(phase) 576 578 cosp = np.cos(phase) 577 occ = Mdata*Fdata/len( Uniq)579 occ = Mdata*Fdata/len(refDict['Uniq'][iref]) 578 580 biso = -SQfactor*Uisodata 579 581 Tiso = np.where(biso<1.,np.exp(biso),1.0) 580 HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])582 HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in refDict['Uniq'][iref]]) 581 583 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 582 584 Tcorr = Tiso*Tuij … … 591 593 refl[10] = atan2d(fbs[0],fas[0]) 592 594 593 def StructureFactor2(refList,G,hfx,pfx,SGData,calcControls,parmDict): 594 ''' Compute structure factors for all h,k,l for phase 595 puts the result, F^2, in each ref[8] in refList - want to do this for blocks of reflections 596 & not one by one. 597 input: 598 599 :param list refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv] 600 :param np.array G: reciprocal metric tensor 601 :param str pfx: phase id string 602 :param dict SGData: space group info. dictionary output from SpcGroup 603 :param dict calcControls: 604 :param dict ParmDict: 605 606 ''' 607 twopi = 2.0*np.pi 608 twopisq = 2.0*np.pi**2 609 phfx = pfx.split(':')[0]+hfx 610 ast = np.sqrt(np.diag(G)) 611 Mast = twopisq*np.multiply.outer(ast,ast) 612 FFtables = calcControls['FFtables'] 613 BLtables = calcControls['BLtables'] 614 Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict) 615 FF = np.zeros(len(Tdata)) 616 if 'N' in calcControls[hfx+'histType']: 617 FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam']) 618 else: 619 FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata]) 620 FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata]) 621 Uij = np.array(G2lat.U6toUij(Uijdata)) 622 bij = Mast*Uij.T 623 for refl in refList: 624 fbs = np.array([0,0]) 625 H = refl[:3] 626 SQ = 1./(2.*refl[4])**2 627 SQfactor = 4.0*SQ*twopisq 628 Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor) 629 if not len(refl[-1]): #no form factors 630 if 'N' in calcControls[hfx+'histType']: 631 refl[-1] = G2el.getBLvalues(BLtables) 632 else: #'X' 633 refl[-1] = G2el.getFFvalues(FFtables,SQ) 634 for i,El in enumerate(Tdata): 635 FF[i] = refl[-1][El] 636 Uniq = refl[11] 637 phi = refl[12] 638 phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis]) 639 sinp = np.sin(phase) 640 cosp = np.cos(phase) 641 occ = Mdata*Fdata/len(Uniq) 642 biso = -SQfactor*Uisodata 643 Tiso = np.where(biso<1.,np.exp(biso),1.0) 644 HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq]) 645 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 646 Tcorr = Tiso*Tuij 647 fa = np.array([(FF+FP-Bab)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr]) 648 fas = np.sum(np.sum(fa,axis=1),axis=1) #real 649 if not SGData['SGInv']: 650 fb = np.array([(FF+FP-Bab)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr]) 651 fbs = np.sum(np.sum(fb,axis=1),axis=1) 652 fasq = fas**2 653 fbsq = fbs**2 #imaginary 654 refl[9] = np.sum(fasq)+np.sum(fbsq) 655 refl[10] = atan2d(fbs[0],fas[0]) 656 657 def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict): 595 def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict): 658 596 'Needs a doc string' 659 597 twopi = 2.0*np.pi … … 675 613 bij = Mast*Uij.T 676 614 dFdvDict = {} 677 dFdfr = np.zeros((len(ref List),len(Mdata)))678 dFdx = np.zeros((len(ref List),len(Mdata),3))679 dFdui = np.zeros((len(ref List),len(Mdata)))680 dFdua = np.zeros((len(ref List),len(Mdata),6))681 dFdbab = np.zeros((len(ref List),2))682 for iref,refl in enumerate(ref List):615 dFdfr = np.zeros((len(refDict['RefList']),len(Mdata))) 616 dFdx = np.zeros((len(refDict['RefList']),len(Mdata),3)) 617 dFdui = np.zeros((len(refDict['RefList']),len(Mdata))) 618 dFdua = np.zeros((len(refDict['RefList']),len(Mdata),6)) 619 dFdbab = np.zeros((len(refDict['RefList']),2)) 620 for iref,refl in enumerate(refDict['RefList']): 683 621 H = np.array(refl[:3]) 684 622 SQ = 1./(2.*refl[4])**2 # or (sin(theta)/lambda)**2 … … 687 625 Bab = parmDict[phfx+'BabA']*dBabdA 688 626 for i,El in enumerate(Tdata): 689 FF[i] = refl[-1][El] 690 Uniq = refl[11] 691 phi = refl[12] 692 phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:]) 627 FF[i] = refDict['FF'][iref][El] 628 phase = twopi*(np.inner((dXdata.T+Xdata.T),refDict['Uniq'][iref])+refDict['Phi'][iref][np.newaxis,:]) 693 629 sinp = np.sin(phase) 694 630 cosp = np.cos(phase) 695 occ = Mdata*Fdata/len( Uniq)631 occ = Mdata*Fdata/len(refDict['Uniq'][iref]) 696 632 biso = -SQfactor*Uisodata 697 633 Tiso = np.where(biso<1.,np.exp(biso),1.0) 698 634 HbH = -np.inner(H,np.inner(bij,H)) 699 Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])635 Hij = np.array([Mast*np.multiply.outer(U,U) for U in refDict['Uniq'][iref]]) 700 636 Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij]) 701 637 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) … … 712 648 #sum below is over Uniq 713 649 dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2) 714 dfadx = np.sum(twopi* Uniq*fax[:,:,:,np.newaxis],axis=2)650 dfadx = np.sum(twopi*refDict['Uniq'][iref]*fax[:,:,:,np.newaxis],axis=2) 715 651 dfadui = np.sum(-SQfactor*fa,axis=2) 716 652 dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2) 717 653 dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1) 718 654 #NB: the above have been checked against PA(1:10,1:2) in strfctr.for 719 dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len( Uniq)655 dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(refDict['Uniq'][iref]) 720 656 dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1]) 721 657 dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1]) … … 724 660 if not SGData['SGInv']: 725 661 dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2) #problem here if occ=0 for some atom 726 dfbdx = np.sum(twopi* Uniq*fbx[:,:,:,np.newaxis],axis=2)662 dfbdx = np.sum(twopi*refDict['Uniq'][iref]*fbx[:,:,:,np.newaxis],axis=2) 727 663 dfbdui = np.sum(-SQfactor*fb,axis=2) 728 664 dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2) 729 665 dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1) 730 dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len( Uniq)666 dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(refDict['Uniq'][iref]) 731 667 dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1]) 732 668 dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1]) … … 750 686 return dFdvDict 751 687 752 def StructureFactorDerv2(refList,G,hfx,pfx,SGData,calcControls,parmDict):753 '''Needs a doc string - want to do this for blocks of reflections754 & not one by one.'''755 twopi = 2.0*np.pi756 twopisq = 2.0*np.pi**2757 phfx = pfx.split(':')[0]+hfx758 ast = np.sqrt(np.diag(G))759 Mast = twopisq*np.multiply.outer(ast,ast)760 FFtables = calcControls['FFtables']761 BLtables = calcControls['BLtables']762 Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)763 FF = np.zeros(len(Tdata))764 if 'N' in calcControls[hfx+'histType']:765 FP = 0.766 FPP = 0.767 else:768 FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])769 FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])770 Uij = np.array(G2lat.U6toUij(Uijdata))771 bij = Mast*Uij.T772 dFdvDict = {}773 dFdfr = np.zeros((len(refList),len(Mdata)))774 dFdx = np.zeros((len(refList),len(Mdata),3))775 dFdui = np.zeros((len(refList),len(Mdata)))776 dFdua = np.zeros((len(refList),len(Mdata),6))777 dFdbab = np.zeros((len(refList),2))778 for iref,refl in enumerate(refList):779 H = np.array(refl[:3])780 SQ = 1./(2.*refl[4])**2 # or (sin(theta)/lambda)**2781 SQfactor = 8.0*SQ*np.pi**2782 dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)783 Bab = parmDict[phfx+'BabA']*dBabdA784 for i,El in enumerate(Tdata):785 FF[i] = refl[-1][El]786 Uniq = refl[11]787 phi = refl[12]788 phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])789 sinp = np.sin(phase)790 cosp = np.cos(phase)791 occ = Mdata*Fdata/len(Uniq)792 biso = -SQfactor*Uisodata793 Tiso = np.where(biso<1.,np.exp(biso),1.0)794 HbH = -np.inner(H,np.inner(bij,H))795 Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])796 Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])797 Tuij = np.where(HbH<1.,np.exp(HbH),1.0)798 Tcorr = Tiso*Tuij799 fot = (FF+FP-Bab)*occ*Tcorr800 fotp = FPP*occ*Tcorr801 fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp]) #non positions802 fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])803 804 fas = np.sum(np.sum(fa,axis=1),axis=1)805 fbs = np.sum(np.sum(fb,axis=1),axis=1)806 fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp]) #positions807 fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])808 #sum below is over Uniq809 dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)810 dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)811 dfadui = np.sum(-SQfactor*fa,axis=2)812 dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)813 dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)814 #NB: the above have been checked against PA(1:10,1:2) in strfctr.for815 dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)816 dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])817 dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])818 dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])819 dFdbab[iref] = np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T820 if not SGData['SGInv']:821 dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2) #problem here if occ=0 for some atom822 dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)823 dfbdui = np.sum(-SQfactor*fb,axis=2)824 dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)825 dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)826 dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)827 dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])828 dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])829 dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])830 dFdbab[iref] += np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T831 #loop over atoms - each dict entry is list of derivatives for all the reflections832 for i in range(len(Mdata)):833 dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]834 dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]835 dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]836 dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]837 dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]838 dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]839 dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]840 dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]841 dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]842 dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]843 dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]844 dFdvDict[pfx+'BabA'] = dFdbab.T[0]845 dFdvDict[pfx+'BabU'] = dFdbab.T[1]846 return dFdvDict847 848 688 def SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varyList): 849 689 ''' Single crystal extinction function; puts correction in ref[13] and returns 850 690 corrections needed for derivatives 851 691 ''' 852 ref[1 3] = 1.0692 ref[11] = 1.0 853 693 dervCor = 1.0 854 694 dervDict = {} … … 910 750 dervDict[phfx+'Eg'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2 911 751 912 ref[1 3] = 1./extCor752 ref[11] = 1./extCor 913 753 return dervCor,dervDict 914 915 754 916 755 def Dict2Values(parmdict, varylist): … … 1043 882 return odfCor,dFdODF 1044 883 1045 def GetPrefOri(refl, G,g,phfx,hfx,SGData,calcControls,parmDict):884 def GetPrefOri(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict): 1046 885 'Needs a doc string' 1047 886 POcorr = 1.0 … … 1051 890 MDAxis = calcControls[phfx+'MDAxis'] 1052 891 sumMD = 0 1053 for H in refl[11]:892 for H in uniq: 1054 893 cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G) 1055 894 A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD) 1056 895 sumMD += A**3 1057 POcorr = sumMD/len( refl[11])896 POcorr = sumMD/len(uniq) 1058 897 else: #spherical harmonics 1059 898 if calcControls[phfx+'SHord']: … … 1061 900 return POcorr 1062 901 1063 def GetPrefOriDerv(refl, G,g,phfx,hfx,SGData,calcControls,parmDict):902 def GetPrefOriDerv(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict): 1064 903 'Needs a doc string' 1065 904 POcorr = 1.0 … … 1070 909 sumMD = 0 1071 910 sumdMD = 0 1072 for H in refl[11]:911 for H in uniq: 1073 912 cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G) 1074 913 A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD) 1075 914 sumMD += A**3 1076 915 sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2) 1077 POcorr = sumMD/len( refl[11])1078 POderv[phfx+'MD'] = sumdMD/len( refl[11])916 POcorr = sumMD/len(uniq) 917 POderv[phfx+'MD'] = sumdMD/len(uniq) 1079 918 else: #spherical harmonics 1080 919 if calcControls[phfx+'SHord']: … … 1096 935 return 0.0 1097 936 1098 def GetIntensityCorr(refl, G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):937 def GetIntensityCorr(refl,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict): 1099 938 'Needs a doc string' 1100 939 Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3] #scale*multiplicity 1101 940 if 'X' in parmDict[hfx+'Type']: 1102 941 Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0] 1103 Icorr *= GetPrefOri(refl, G,g,phfx,hfx,SGData,calcControls,parmDict)942 Icorr *= GetPrefOri(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict) 1104 943 if pfx+'SHorder' in parmDict: 1105 944 Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict) 1106 945 Icorr *= GetAbsorb(refl,hfx,calcControls,parmDict) 1107 refl[1 3] = Icorr1108 1109 def GetIntensityDerv(refl, G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):946 refl[11] = Icorr 947 948 def GetIntensityDerv(refl,uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict): 1110 949 'Needs a doc string' 1111 950 dIdsh = 1./parmDict[hfx+'Scale'] … … 1116 955 else: #'N' 1117 956 dIdPola = 0.0 1118 POcorr,dIdPO = GetPrefOriDerv(refl, G,g,phfx,hfx,SGData,calcControls,parmDict)957 POcorr,dIdPO = GetPrefOriDerv(refl,uniq,G,g,phfx,hfx,SGData,calcControls,parmDict) 1119 958 for iPO in dIdPO: 1120 959 dIdPO[iPO] /= POcorr … … 1372 1211 pId = Phase['pId'] 1373 1212 phfx = '%d:%d:'%(pId,hId) 1374 ref List = refLists[phase]1213 refDict = refLists[phase] 1375 1214 sumFo = 0.0 1376 1215 sumdF = 0.0 1377 1216 sumFosq = 0.0 1378 1217 sumdFsq = 0.0 1379 for refl in ref List:1218 for refl in refDict['RefList']: 1380 1219 if 'C' in calcControls[hfx+'histType']: 1381 1220 yp = np.zeros_like(yb) … … 1384 1223 iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF)) 1385 1224 iFin2 = iFin 1386 yp[iBeg:iFin] = refl[1 3]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here1225 yp[iBeg:iFin] = refl[11]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here 1387 1226 if Ka2: 1388 1227 pos2 = refl[5]+lamRatio*tand(refl[5]/2.0) # + 360/pi * Dlam/lam * tan(th) … … 1394 1233 elif not iBeg2-iFin2: #peak above high limit - done 1395 1234 break 1396 yp[iBeg2:iFin2] += refl[1 3]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) #and here1397 refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[1 3]*(1.+kRatio)),0.0))1235 yp[iBeg2:iFin2] += refl[11]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) #and here 1236 refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[11]*(1.+kRatio)),0.0)) 1398 1237 elif 'T' in calcControls[hfx+'histType']: 1399 1238 print 'TOF Undefined at present' … … 1407 1246 Histogram['Residuals'][phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.) 1408 1247 Histogram['Residuals'][phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.) 1409 Histogram['Residuals'][phfx+'Nref'] = len(ref List)1248 Histogram['Residuals'][phfx+'Nref'] = len(refDict['RefList']) 1410 1249 Histogram['Residuals']['hId'] = hId 1411 1250 elif 'HKLF' in histogram[:4]: … … 1450 1289 raise ValueError 1451 1290 for phase in Histogram['Reflection Lists']: 1452 ref List = Histogram['Reflection Lists'][phase]1291 refDict = Histogram['Reflection Lists'][phase] 1453 1292 Phase = Phases[phase] 1454 1293 pId = Phase['pId'] … … 1463 1302 if not Phase['General'].get('doPawley'): 1464 1303 time0 = time.time() 1465 StructureFactor(ref List,G,hfx,pfx,SGData,calcControls,parmDict)1304 StructureFactor(refDict,G,hfx,pfx,SGData,calcControls,parmDict) 1466 1305 print 'sf calc time: %.3fs'%(time.time()-time0) 1467 1306 time0 = time.time() 1468 for refl in refList: 1307 Uniq = refDict['Uniq'] 1308 for iref,refl in enumerate(refDict['RefList']): 1469 1309 if 'C' in calcControls[hfx+'histType']: 1470 1310 h,k,l = refl[:3] … … 1473 1313 refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift 1474 1314 refl[6:8] = GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict) #peak sig & gam 1475 GetIntensityCorr(refl, G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) #puts corrections in refl[13]1476 refl[1 3] *= Vst*Lorenz1315 GetIntensityCorr(refl,Uniq[iref],G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) #puts corrections in refl[11] 1316 refl[11] *= Vst*Lorenz 1477 1317 if Phase['General'].get('doPawley'): 1478 1318 try: … … 1489 1329 elif not iBeg-iFin: #peak above high limit - done 1490 1330 break 1491 yc[iBeg:iFin] += refl[1 3]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin])) #>90% of time spent here1331 yc[iBeg:iFin] += refl[11]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin])) #>90% of time spent here 1492 1332 if Ka2: 1493 1333 pos2 = refl[5]+lamRatio*tand(refl[5]/2.0) # + 360/pi * Dlam/lam * tan(th) … … 1499 1339 elif not iBeg-iFin: #peak above high limit - done 1500 1340 return yc,yb 1501 yc[iBeg:iFin] += refl[1 3]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin])) #and here1341 yc[iBeg:iFin] += refl[11]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin])) #and here 1502 1342 elif 'T' in calcControls[hfx+'histType']: 1503 1343 print 'TOF Undefined at present' … … 1576 1416 raise ValueError 1577 1417 for phase in Histogram['Reflection Lists']: 1578 ref List = Histogram['Reflection Lists'][phase]1418 refDict = Histogram['Reflection Lists'][phase] 1579 1419 Phase = Phases[phase] 1580 1420 SGData = Phase['General']['SGData'] … … 1587 1427 if not Phase['General'].get('doPawley'): 1588 1428 time0 = time.time() 1589 dFdvDict = StructureFactorDerv(ref List,G,hfx,pfx,SGData,calcControls,parmDict)1429 dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict) 1590 1430 print 'sf-derv time %.3fs'%(time.time()-time0) 1591 1431 ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase) 1592 1432 time0 = time.time() 1593 for iref,refl in enumerate(refList): 1433 Uniq = refDict['Uniq'] 1434 for iref,refl in enumerate(refDict['RefList']): 1594 1435 if 'C' in calcControls[hfx+'histType']: #CW powder 1595 1436 h,k,l = refl[:3] 1596 dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl, G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)1437 dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl,Uniq[iref],G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) 1597 1438 Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl) 1598 1439 iBeg = np.searchsorted(x,refl[5]-fmin) … … 1609 1450 dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,ma.getdata(x[iBeg:iFin])) 1610 1451 for i in range(5): 1611 dMdpk[i] += 100.*cw[iBeg:iFin]*refl[1 3]*refl[9]*dMdipk[i]1452 dMdpk[i] += 100.*cw[iBeg:iFin]*refl[11]*refl[9]*dMdipk[i] 1612 1453 dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])} 1613 1454 if Ka2: … … 1620 1461 dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,ma.getdata(x[iBeg2:iFin2])) 1621 1462 for i in range(5): 1622 dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[1 3]*refl[9]*kRatio*dMdipk2[i]1623 dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[1 3]*dMdipk2[0]1463 dMdpk2[i] = 100.*cw[iBeg2:iFin2]*refl[11]*refl[9]*kRatio*dMdipk2[i] 1464 dMdpk2[5] = 100.*cw[iBeg2:iFin2]*refl[11]*dMdipk2[0] 1624 1465 dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]} 1625 1466 if Phase['General'].get('doPawley'): … … 1801 1642 A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] 1802 1643 G,g = G2lat.A2Gmat(A) #recip & real metric tensors 1803 ref List = Histogram['Data']1804 dFdvDict = StructureFactorDerv(ref List,G,hfx,pfx,SGData,calcControls,parmDict)1644 refDict = Histogram['Data'] 1645 dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict) 1805 1646 ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase) 1806 dMdvh = np.zeros((len(varylist),len(ref List)))1647 dMdvh = np.zeros((len(varylist),len(refDict['RefList']))) 1807 1648 dependentVars = G2mv.GetDependentVars() 1808 1649 depDerivDict = {} 1809 1650 for j in dependentVars: 1810 depDerivDict[j] = np.zeros(shape=(len(ref List)))1651 depDerivDict[j] = np.zeros(shape=(len(refDict['RefList']))) 1811 1652 if calcControls['F**2']: 1812 for iref,ref in enumerate(ref List):1653 for iref,ref in enumerate(refDict['RefList']): 1813 1654 if ref[6] > 0: 1814 1655 dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13] … … 1836 1677 depDerivDict[phfx+item][iref] = w*dervCor*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale'] 1837 1678 else: 1838 for iref,ref in enumerate(ref List):1679 for iref,ref in enumerate(refDict['RefList']): 1839 1680 if ref[5] > 0.: 1840 1681 dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13] … … 1938 1779 A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] 1939 1780 G,g = G2lat.A2Gmat(A) #recip & real metric tensors 1940 ref List = Histogram['Data']1781 refDict = Histogram['Data'] 1941 1782 time0 = time.time() 1942 dFdvDict = StructureFactorDerv(ref List,G,hfx,pfx,SGData,calcControls,parmDict) #accurate for powders!1783 dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict) #accurate for powders! 1943 1784 ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase) 1944 dMdvh = np.zeros((len(varylist),len(ref List)))1785 dMdvh = np.zeros((len(varylist),len(refDict['RefList']))) 1945 1786 dependentVars = G2mv.GetDependentVars() 1946 1787 depDerivDict = {} 1947 1788 for j in dependentVars: 1948 depDerivDict[j] = np.zeros(shape=(len(ref List)))1949 wdf = np.zeros(len(ref List))1789 depDerivDict[j] = np.zeros(shape=(len(refDict['RefList']))) 1790 wdf = np.zeros(len(refDict['RefList'])) 1950 1791 if calcControls['F**2']: 1951 for iref,ref in enumerate(ref List):1792 for iref,ref in enumerate(refDict['RefList']): 1952 1793 if ref[6] > 0: 1953 1794 dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13] … … 1976 1817 depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*dervCor 1977 1818 else: 1978 for iref,ref in enumerate(ref List):1819 for iref,ref in enumerate(refDict['RefList']): 1979 1820 if ref[5] > 0.: 1980 1821 dervCor,dervDict = SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13] … … 2092 1933 A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] 2093 1934 G,g = G2lat.A2Gmat(A) #recip & real metric tensors 2094 ref List = Histogram['Data']1935 refDict = Histogram['Data'] 2095 1936 time0 = time.time() 2096 StructureFactor(ref List,G,hfx,pfx,SGData,calcControls,parmDict)1937 StructureFactor(refDict,G,hfx,pfx,SGData,calcControls,parmDict) 2097 1938 # print 'sf-calc time: %.3f'%(time.time()-time0) 2098 df = np.zeros(len(ref List))1939 df = np.zeros(len(refDict['RefList'])) 2099 1940 sumwYo = 0 2100 1941 sumFo = 0 … … 2104 1945 nobs = 0 2105 1946 if calcControls['F**2']: 2106 for i,ref in enumerate(ref List):1947 for i,ref in enumerate(refDict['RefList']): 2107 1948 if ref[6] > 0: 2108 1949 SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13] 2109 1950 w = 1.0/ref[6] 2110 1951 ref[7] = parmDict[phfx+'Scale']*ref[9] 2111 ref[7] *= ref[1 3] #correct Fc^2 for extinction1952 ref[7] *= ref[11] #correct Fc^2 for extinction 2112 1953 ref[8] = ref[5]/parmDict[phfx+'Scale'] 2113 1954 if w*ref[5] >= calcControls['minF/sig']: … … 2122 1963 sumwYo += (w*ref[5])**2 2123 1964 else: 2124 for i,ref in enumerate(ref List):1965 for i,ref in enumerate(refDict['RefList']): 2125 1966 if ref[5] > 0.: 2126 1967 SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varylist) #puts correction in refl[13] 2127 1968 ref[7] = parmDict[phfx+'Scale']*ref[9] 2128 ref[7] *= ref[1 3] #correct Fc^2 for extinction1969 ref[7] *= ref[11] #correct Fc^2 for extinction 2129 1970 Fo = np.sqrt(ref[5]) 2130 1971 Fc = np.sqrt(ref[7]) -
TabularUnified trunk/imports/G2sfact.py ¶
r1077 r1106 51 51 sigFo = float(sigFo) 52 52 # h,k,l,m,dsp,Fo2,sig,Fc2,Fot2,Fct2,phase,... 53 self.RefList.append([h,k,l,0,0,Fo**2,2.*Fo*sigFo,0,Fo**2,0,0,[],[],0,{}]) 53 self.RefDict['RefList'].append([h,k,l,0,0,Fo**2,2.*Fo*sigFo,0,Fo**2,0,0,0]) 54 self.RefDict['Uniq'].append([]) 55 self.RefDict['Phi'].append([]) 56 self.RefDict['FF'].append({}) 54 57 self.UpdateControls(Type='Fosq',FcalcPresent=False) # set Fobs type & if Fcalc values are loaded 55 58 self.UpdateParameters(Type='SXC',Wave=None) # histogram type … … 97 100 sigFo = float(sigFo) 98 101 # h,k,l,m,dsp,Fo2,sig,Fc2,Fot2,Fct2,phase,... 99 self.RefList.append([h,k,l,0,0,Fo,sigFo,0,Fo,0,0,[],[],0,{}]) 102 self.RefDict['RefList'].append([h,k,l,0,0,Fo,sigFo,0,Fo,0,0,0]) 103 self.RefDict['Uniq'].append([]) 104 self.RefDict['Phi'].append([]) 105 self.RefDict['FF'].append({}) 100 106 self.UpdateControls(Type='Fosq',FcalcPresent=False) # set Fobs type & if Fcalc values are loaded 101 107 self.UpdateParameters(Type='SXC',Wave=None) # histogram type -
TabularUnified trunk/imports/G2sfact_CIF.py ¶
r1077 r1106 155 155 #h,k,l,m,dsp,Fo2,sig,Fc2,Fot2,Fct2,phase,... 156 156 ref = HKL+[0,0,0,0,0,0,0,0,[],[],0,{}] 157 # ref = HKL+[0,0,0,0,0,0,0,0,0] 157 158 if '_refln_f_squared_meas' in itemkeys: 158 159 try: … … 212 213 except: 213 214 pass 214 self.RefList.append(ref) 215 self.RefDict['RefList'].append(ref) 216 self.RefDict['Uniq'].append([]) 217 self.RefDict['Phi'].append([]) 218 self.RefDict['FF'].append({}) 215 219 self.UpdateControls(Type='Fosq',FcalcPresent=FcalcPresent) # set Fobs type & if Fcalc values are loaded 216 220 if blk.get('_diffrn_radiation_probe'):
Note: See TracChangeset
for help on using the changeset viewer.