Changeset 2110
- Timestamp:
- Jan 3, 2016 12:40:04 PM (7 years ago)
- Location:
- trunk
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIgrid.py
r2109 r2110 68 68 wxID_RELOADDRAWATOMS,wxID_ATOMSDISAGL,wxID_ATOMMOVE,wxID_MAKEMOLECULE, 69 69 wxID_ASSIGNATMS2RB,wxID_ATOMSPDISAGL, wxID_ISODISP,wxID_ADDHATOM,wxID_UPDATEHATOM, 70 ] = [wx.NewId() for item in range(17)] 70 wxID_WAVEVARY, 71 ] = [wx.NewId() for item in range(18)] 71 72 72 73 [ wxID_DRAWATOMSTYLE, wxID_DRAWATOMLABEL, wxID_DRAWATOMCOLOR, wxID_DRAWATOMRESETCOLOR, … … 1292 1293 self.PostfillDataMenu() 1293 1294 1294 # Phase / Imcommensurate "waves" tab 1295 # Phase / Imcommensurate "waves" tab 1295 1296 self.WavesData = wx.MenuBar() 1296 1297 self.PrefillDataMenu(self.WavesData,helpType='Wave Data', helpLbl='Imcommensurate wave data') 1297 1298 self.WavesData.Append(menu=wx.Menu(title=''),title='Select tab') 1298 self.WavesDataCompute = wx.Menu(title='') 1299 self.WavesDataEdit = wx.Menu(title='') 1300 self.WavesData.Append(menu=self.WavesDataEdit, title='Edit') 1301 self.WavesDataEdit.Append(id=wxID_WAVEVARY, kind=wx.ITEM_NORMAL,text='Global wave vary', 1302 help='Global setting of wave vary flags') 1299 1303 self.PostfillDataMenu() 1300 1304 -
trunk/GSASIImath.py
r2089 r2110 1082 1082 ''' 1083 1083 H: array ops X hklt proj to hkl 1084 HP: array nRefBlk xops X hklt proj to hkl1084 HP: array ops X hklt proj to hkl 1085 1085 Hij: array 2pi^2[a*^2h^2 b*^2k^2 c*^2l^2 a*b*hk a*c*hl b*c*kl] of projected hklm to hkl space 1086 1086 ''' … … 1133 1133 dGdMzSt = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXZt*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2) 1134 1134 dGdMzSx = np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZx*np.cos(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2) 1135 dGdMzS = np.concatenate((np.sum(dGdMzSt,axis=-1),dGdMzSx),axis=-1) 1136 # GSASIIpath.IPyBreak() 1137 return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS],[dGdMzC,dGdMzS] 1138 1139 def ModulationDerv2(H,HP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt): 1140 ''' 1141 H: array refBlk x ops X hklt proj to hkl 1142 HP: array refBlk x ops X hklt proj to hkl 1143 Hij: array 2pi^2[a*^2h^2 b*^2k^2 c*^2l^2 a*b*hk a*c*hl b*c*kl] of projected hklm to hkl space 1144 ''' 1145 1146 Mf = [H.shape[0],]+list(waveShapes[0]) #=[ops,atoms,waves,2] (sin+cos frac mods) 1147 dGdMfC = np.zeros(Mf) 1148 dGdMfS = np.zeros(Mf) 1149 Mx = [H.shape[0],]+list(waveShapes[1]) #=[ops,atoms,waves,6] (sin+cos pos mods) 1150 dGdMxC = np.zeros(Mx) 1151 dGdMxS = np.zeros(Mx) 1152 Mu = [H.shape[0],]+list(waveShapes[2]) #=[ops,atoms,waves,12] (sin+cos Uij mods) 1153 dGdMuC = np.zeros(Mu) 1154 dGdMuS = np.zeros(Mu) 1155 1156 D = twopi*H[:,:,3,nxs]*glTau[nxs,nxs,:] #m*e*tau; refBlk x ops X ngl 1157 HdotX = twopi*np.inner(HP,Xmod) #refBlk x ops x atoms X ngl 1158 HdotXD = HdotX+D[:,:,nxs,:] 1159 if nWaves[2]: 1160 Umod = np.swapaxes((UmodAB),2,4) #atoms x waves x ngl x 3x3 (symmetric so I can do this!) 1161 HuH = np.sum(HP[:,:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1) #refBlk x ops x atoms x waves x ngl 1162 HuH = np.sum(HP[:,:,nxs,nxs,nxs]*np.inner(HP,Umod),axis=-1) #refBlk x ops x atoms x waves x ngl 1163 HbH = np.exp(-np.sum(HuH,axis=-2)) #refBlk x ops x atoms x ngl; sum waves - OK vs Modulation version 1164 part1 = -np.exp(-HuH)*Fmod #refBlk x ops x atoms x waves x ngl 1165 dUdAu = Hij[:,:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[0]),0,4)[nxs,nxs,:,:,:,:] #ops x atoms x waves x ngl x 6sinUij 1166 dUdBu = Hij[:,:,nxs,nxs,nxs,:]*np.rollaxis(G2lat.UijtoU6(SCtauU[1]),0,4)[nxs,nxs,:,:,:,:] #ops x atoms x waves x ngl x 6cosUij 1167 dGdMuCa = np.sum(part1[:,:,:,:,:,nxs]*dUdAu*np.cos(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2) #ops x atoms x waves x 6uij; G-L sum 1168 dGdMuCb = np.sum(part1[:,:,:,:,:,nxs]*dUdBu*np.cos(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2) 1169 dGdMuC = np.concatenate((dGdMuCa,dGdMuCb),axis=-1) #ops x atoms x waves x 12uij 1170 dGdMuSa = np.sum(part1[:,:,:,:,:,nxs]*dUdAu*np.sin(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2) #ops x atoms x waves x 6uij; G-L sum 1171 dGdMuSb = np.sum(part1[:,:,:,:,:,nxs]*dUdBu*np.sin(HdotXD)[:,:,:,nxs,:,nxs]*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2) 1172 dGdMuS = np.concatenate((dGdMuSa,dGdMuSb),axis=-1) #ops x atoms x waves x 12uij 1173 else: 1174 HbH = np.ones_like(HdotX) 1175 dHdXA = twopi*HP[:,:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[0],-1,-2)[nxs,nxs,:,:,:,:] #ops x atoms x sine waves x ngl x xyz 1176 dHdXB = twopi*HP[:,:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[1],-1,-2)[nxs,nxs,:,:,:,:] #ditto - cos waves 1177 # ops x atoms x waves x 2xyz - real part - good 1178 dGdMxCa = -np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXA*np.sin(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2) 1179 dGdMxCb = -np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXB*np.sin(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2) 1180 dGdMxC = np.concatenate((dGdMxCa,dGdMxCb),axis=-1) 1181 # ops x atoms x waves x 2xyz - imag part - good 1182 dGdMxSa = np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXA*np.cos(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2) 1183 dGdMxSb = np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXB*np.cos(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2) 1184 dGdMxS = np.concatenate((dGdMxSa,dGdMxSb),axis=-1) 1185 # ZigZag/Block waves - problems here? 1186 dHdXZt = -twopi*HP[:,:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[2],-1,-2)[nxs,nxs,:,:,:,:] #ops x atoms x ngl x 2(ZigZag/Block Tminmax) 1187 dHdXZx = twopi*HP[:,:,nxs,nxs,:]*np.swapaxes(SCtauX[3],-1,-2)[nxs,nxs,:,:,:] #ops x atoms x ngl x 3(ZigZag/Block XYZmax) 1188 dGdMzCt = -np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXZt*np.sin(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2) 1189 dGdMzCx = -np.sum((Fmod*HbH)[:,:,:,:,nxs]*(dHdXZx*np.sin(HdotXD)[:,:,:,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2) 1190 dGdMzC = np.concatenate((np.sum(dGdMzCt,axis=-1),dGdMzCx),axis=-1) 1191 dGdMzSt = np.sum((Fmod*HbH)[:,:,:,nxs,:,nxs]*(dHdXZt*np.cos(HdotXD)[:,:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,nxs,:,nxs],axis=-2) 1192 dGdMzSx = np.sum((Fmod*HbH)[:,:,:,:,nxs]*(dHdXZx*np.cos(HdotXD)[:,:,:,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2) 1135 1193 dGdMzS = np.concatenate((np.sum(dGdMzSt,axis=-1),dGdMzSx),axis=-1) 1136 1194 # GSASIIpath.IPyBreak() -
trunk/GSASIIphsGUI.py
r2109 r2110 2547 2547 G2frame.bottomSizer = ShowAtomInfo() 2548 2548 mainSizer.Add(G2frame.bottomSizer) 2549 2549 #wxID_WAVEVARY 2550 2550 SetPhaseWindow(G2frame.dataFrame,G2frame.waveData,mainSizer,Scroll) 2551 2552 def OnWaveVary(event): 2553 print 'set vary flags for all waves - TBD' 2551 2554 2552 2555 ################################################################################ … … 6397 6400 if data['General']['Type'] in ['modulated','magnetic']: 6398 6401 FillSelectPageMenu(TabSelectionIdDict, G2frame.dataFrame.WavesData) 6402 G2frame.dataFrame.Bind(wx.EVT_MENU, OnWaveVary, id=G2gd.wxID_WAVEVARY) 6399 6403 # Draw Options 6400 6404 FillSelectPageMenu(TabSelectionIdDict, G2frame.dataFrame.DataDrawOptions) -
trunk/GSASIIstrIO.py
r2070 r2110 2407 2407 for i,name in enumerate(Snames): 2408 2408 ptlbls += '%12s' % (name) 2409 ptstr += '%12. 6f' % (hapData[4][i])2409 ptstr += '%12.1f' % (hapData[4][i]) 2410 2410 if mustrainSig[1][i]: 2411 2411 refine = True 2412 sigstr += '%12. 6f' % (mustrainSig[1][i])2412 sigstr += '%12.1f' % (mustrainSig[1][i]) 2413 2413 else: 2414 2414 sigstr += 12*' ' -
trunk/GSASIIstrMath.py
r2097 r2110 632 632 return np.array(waveTypes),FSSdata,XSSdata,USSdata,MSSdata 633 633 634 def GetSSTauM(SGOps,SSOps,pfx,calcControls,XData):635 636 Natoms = calcControls['Natoms'][pfx]637 maxSSwave = calcControls['maxSSwave'][pfx]638 Smult = np.zeros((Natoms,len(SGOps)))639 TauT = np.zeros((Natoms,len(SGOps)))640 for ix,xyz in enumerate(XData.T):641 for isym,(sop,ssop) in enumerate(zip(SGOps,SSOps)):642 sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(0,sop,ssop,xyz)643 Smult[ix][isym] = sdet644 TauT[ix][isym] = tauT645 return Smult,TauT646 634 #def GetSSTauM(SGOps,SSOps,pfx,calcControls,XData): 635 # 636 # Natoms = calcControls['Natoms'][pfx] 637 # maxSSwave = calcControls['maxSSwave'][pfx] 638 # Smult = np.zeros((Natoms,len(SGOps))) 639 # TauT = np.zeros((Natoms,len(SGOps))) 640 # for ix,xyz in enumerate(XData.T): 641 # for isym,(sop,ssop) in enumerate(zip(SGOps,SSOps)): 642 # sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(0,sop,ssop,xyz) 643 # Smult[ix][isym] = sdet 644 # TauT[ix][isym] = tauT 645 # return Smult,TauT 646 # 647 647 def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict): 648 648 ''' Compute structure factors for all h,k,l for phase … … 742 742 Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T 743 743 Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT) 744 if 'T' in calcControls[hfx+'histType']: 744 if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms 745 745 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr]) 746 746 fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr]) … … 748 748 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) 749 749 fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr]) 750 fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & uniquehkl751 fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl750 fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real 2 x blkSize x nTwin; sum over atoms & uniq hkl 751 fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag 752 752 if SGData['SGInv']: #centrosymmetric; B=0 753 753 fbs[0] *= 0. 754 if 'P' in calcControls[hfx+'histType']: 755 refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) 754 fas[1] *= 0. 755 if 'PWDR' in calcControls[hfx+'histType']: #PWDR: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2 756 refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) 756 757 refl.T[10] = atan2d(fbs[0],fas[0]) #ignore f' & f" 757 else: 758 else: #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2 758 759 if len(TwinLaw) > 1: 759 760 refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2 #FcT from primary twin element … … 761 762 np.sum(TwinFr*np.sum(TwMask[nxs,:,:]*fbs,axis=0)**2,axis=-1) #Fc sum over twins 762 763 refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0]) #ignore f' & f" 763 else: 764 # GSASIIpath.IPyBreak() 765 else: # checked correct!! 764 766 refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2 765 767 refl.T[7] = np.copy(refl.T[9]) 766 768 refl.T[10] = atan2d(fbs[0],fas[0]) #ignore f' & f" 767 # refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"769 # refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f" 768 770 iBeg += blkSize 769 771 # print ' %d sf time %.4f\r'%(nRef,time.time()-time0) … … 772 774 '''Compute structure factor derivatives on single reflections - keep as it works for twins 773 775 but is slower for powders/nontwins 776 input: 777 778 :param dict refDict: where 779 'RefList' list where each ref = h,k,l,it,d,... 780 'FF' dict of form factors - filled in below 781 :param np.array G: reciprocal metric tensor 782 :param str hfx: histogram id string 783 :param str pfx: phase id string 784 :param dict SGData: space group info. dictionary output from SpcGroup 785 :param dict calcControls: 786 :param dict ParmDict: 787 788 :returns: dict dFdvDict: dictionary of derivatives 774 789 ''' 775 790 phfx = pfx.split(':')[0]+hfx … … 914 929 dFdua[iref] = [2.*TwMask[it]*(SA[it]*dfadua[it][0]+SA[it]*dfadua[it][1]+SB[it]*dfbdua[it][0]+SB[it]*dfbdua[it][1]) for it in range(nTwin)] 915 930 dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2 916 917 931 else: #these are good for no twin single crystals 918 932 dFdfr[iref] = (2.*SA*(dfadfr[0]+dfadfr[1])+2.*SB*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(Uniq) … … 964 978 def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict): 965 979 '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only 966 faster than StructureFactorDerv 980 faster than StructureFactorDerv - correct for powders/nontwins!! 981 input: 982 983 :param dict refDict: where 984 'RefList' list where each ref = h,k,l,it,d,... 985 'FF' dict of form factors - filled in below 986 :param np.array G: reciprocal metric tensor 987 :param str hfx: histogram id string 988 :param str pfx: phase id string 989 :param dict SGData: space group info. dictionary output from SpcGroup 990 :param dict calcControls: 991 :param dict parmDict: 992 993 :returns: dict dFdvDict: dictionary of derivatives 967 994 ''' 968 995 phfx = pfx.split(':')[0]+hfx … … 1033 1060 fotp = FPP*Tcorr 1034 1061 if 'T' in calcControls[hfx+'histType']: 1035 fa = np.array([ np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])1036 fb = np.array([ np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])1062 fa = np.array([fot*cosp,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr]) 1063 fb = np.array([fot*sinp,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr]) 1037 1064 else: 1038 fa = np.array([ np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])1039 fb = np.array([ np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])1065 fa = np.array([fot*cosp,-Flack*FPP*sinp*Tcorr]) 1066 fb = np.array([fot*sinp,Flack*FPP*cosp*Tcorr]) 1040 1067 # GSASIIpath.IPyBreak() 1041 1068 fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl array(2,refBlk,nTwins) … … 1108 1135 1109 1136 def StructureFactorDervTw(refDict,G,hfx,pfx,SGData,calcControls,parmDict): 1110 '''Compute structure factor derivatives on blocks of reflections - for twins only 1111 faster than StructureFactorDerv - not working yet 1137 '''Compute structure factor derivatives on single reflections - for twins only 1138 input: 1139 1140 :param dict refDict: where 1141 'RefList' list where each ref = h,k,l,it,d,... 1142 'FF' dict of form factors - filled in below 1143 :param np.array G: reciprocal metric tensor 1144 :param str hfx: histogram id string 1145 :param str pfx: phase id string 1146 :param dict SGData: space group info. dictionary output from SpcGroup 1147 :param dict calcControls: 1148 :param dict parmDict: 1149 1150 :returns: dict dFdvDict: dictionary of derivatives 1112 1151 ''' 1113 1152 phfx = pfx.split(':')[0]+hfx … … 1118 1157 FFtables = calcControls['FFtables'] 1119 1158 BLtables = calcControls['BLtables'] 1120 TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])1121 1159 TwDict = refDict.get('TwDict',{}) 1122 if 'S' in calcControls[hfx+'histType']: 1123 NTL = calcControls[phfx+'NTL'] 1124 NM = calcControls[phfx+'TwinNMN']+1 1125 TwinLaw = calcControls[phfx+'TwinLaw'] 1126 TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))]) 1127 TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1)) 1160 NTL = calcControls[phfx+'NTL'] 1161 NM = calcControls[phfx+'TwinNMN']+1 1162 TwinLaw = calcControls[phfx+'TwinLaw'] 1163 TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))]) 1164 TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1)) 1128 1165 nTwin = len(TwinLaw) 1129 1166 nRef = len(refDict['RefList']) … … 1144 1181 dFdua = np.zeros((nRef,nTwin,mSize,6)) 1145 1182 dFdbab = np.zeros((nRef,nTwin,2)) 1146 dFdfl = np.zeros((nRef,nTwin))1147 1183 dFdtw = np.zeros((nRef,nTwin)) 1148 Flack = 1.01149 if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:1150 Flack = 1.-2.*parmDict[phfx+'Flack']1151 1184 time0 = time.time() 1152 1185 nref = len(refDict['RefList'])/100 1153 #reflection processing begins here - big arrays! 1154 iBeg = 0 1155 blkSize = 32 #no. of reflections in a block - optimized for speed 1156 while iBeg < nRef: 1157 iFin = min(iBeg+blkSize,nRef) 1158 refl = refDict['RefList'][iBeg:iFin] #array(blkSize,nItems) 1159 H = refl.T[:3] 1186 for iref,refl in enumerate(refDict['RefList']): 1187 if 'T' in calcControls[hfx+'histType']: 1188 FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12]) 1189 H = np.array(refl[:3]) 1160 1190 H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(3,nTwins) or (3) 1161 1191 TwMask = np.any(H,axis=-1) 1162 for ir in range(blkSize): 1163 iref = ir+iBeg 1164 if iref in TwDict: 1165 for i in TwDict[iref]: 1166 for n in range(NTL): 1167 H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM]) 1168 TwMask = np.any(H,axis=-1) 1169 SQ = 1./(2.*refl.T[4])**2 # or (sin(theta)/lambda)**2 1192 if iref in TwDict: 1193 for i in TwDict[iref]: 1194 for n in range(NTL): 1195 H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM]) 1196 TwMask = np.any(H,axis=-1) 1197 SQ = 1./(2.*refl[4])**2 # or (sin(theta)/lambda)**2 1170 1198 SQfactor = 8.0*SQ*np.pi**2 1171 if 'T' in calcControls[hfx+'histType']:1172 if 'P' in calcControls[hfx+'histType']:1173 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])1174 else:1175 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])1176 FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)1177 FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)1178 1199 dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor) 1179 Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))1200 Bab = parmDict[phfx+'BabA']*dBabdA 1180 1201 Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) 1181 FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)1202 FF = refDict['FF']['FF'][iref].T[Tindx].T 1182 1203 Uniq = np.inner(H,SGMT) # array(nSGOp,3) or (nTwin,nSGOp,3) 1183 1204 Phi = np.inner(H,SGT) 1184 1205 phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T 1185 sinp = np.sin(phase) #refBlk x nOps x nAtoms1206 sinp = np.sin(phase) 1186 1207 cosp = np.cos(phase) 1187 1208 occ = Mdata*Fdata/len(SGT) 1188 1209 biso = -SQfactor*Uisodata[:,nxs] 1189 Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)* len(TwinLaw),axis=1).T1210 Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1) 1190 1211 HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1) 1191 Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T1192 Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)1193 1212 Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))]) 1194 Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,nTwin,len(SGT),6)) 1195 fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr 1213 Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6))) 1214 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 1215 Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ 1216 fot = (FF+FP-Bab)*Tcorr 1196 1217 fotp = FPP*Tcorr 1197 if 'T' in calcControls[hfx+'histType']: 1198 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr]) 1199 fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr]) 1200 else: 1201 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) 1202 fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr]) 1218 fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-FPP*sinp*Tcorr]) 1219 fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,FPP*cosp*Tcorr]) 1203 1220 # GSASIIpath.IPyBreak() 1204 fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl array(2, refBlk,nTwins)1221 fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl array(2,nTwins) 1205 1222 fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl 1206 fax = np.array([-fot*sinp,-fotp*cosp]) #positions array(2,refBlk,ntwi,nEqv,nAtoms) 1223 if SGData['SGInv']: #centrosymmetric; B=0 1224 fbs[0] *= 0. 1225 fas[1] *= 0. 1226 fax = np.array([-fot*sinp,-fotp*cosp]) #positions array(2,ntwi,nEqv,nAtoms) 1207 1227 fbx = np.array([fot*cosp,-fotp*sinp]) 1208 1228 #sum below is over Uniq 1209 dfadfr = np.sum(fa/occ,axis=-2) #array(2,refBlk,ntwin,nAtom) Fdata != 0 avoids /0. problem 1210 dfadba = np.sum(-cosp*Tcorr,axis=-2) #array(refBlk,nAtom) 1211 dfadfr = np.swapaxes(dfadfr,0,1) 1212 dfadx = np.swapaxes(np.array([np.sum(twopi*Uniq[it,nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[it,:,:,:,:,nxs],axis=-2) for it in range(nTwin)]),0,1) 1213 dfadui = np.swapaxes(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fa,axis=-2),0,1) #array(Ops,refBlk,nTw,nAtoms) 1214 dfadua = np.swapaxes(np.array([np.sum(-Hij[it,nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[it,:,:,:,:,nxs],axis=-2) for it in range(nTwin)]),0,1) 1215 # array(nTwin,2,refBlk,nAtom,3) & array(nTwin,2,refBlk,nAtom,6) 1229 dfadfr = np.sum(fa/occ,axis=-2) #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem 1230 dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1) 1231 dfadui = np.sum(-SQfactor*fa,axis=-2) 1232 dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1233 dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1234 # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6) 1216 1235 if not SGData['SGInv']: 1217 1236 dfbdfr = np.sum(fb/occ,axis=-2) #Fdata != 0 avoids /0. problem 1218 df bdba = np.sum(-sinp*Tcorr,axis=-2)1219 dfbd fr = np.swapaxes(dfbdfr,0,1)1220 dfbd x = np.swapaxes(np.array([np.sum(twopi*Uniq[it,nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[it,:,:,:,:,nxs],axis=-2) for it in range(nTwin)]),0,1)1221 dfbd ui = np.swapaxes(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fb,axis=-2),0,1) #array(Ops,refBlk,nTw,nAtoms)1222 dfbdua = np. swapaxes(np.array([np.sum(-Hij[it,nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[it,:,:,:,:,nxs],axis=-2) for it in range(nTwin)]),0,1)1237 dfadba /= 2. 1238 dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2. 1239 dfbdui = np.sum(-SQfactor*fb,axis=-2) 1240 dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)]) 1241 dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)]) 1223 1242 else: 1224 1243 dfbdfr = np.zeros_like(dfadfr) … … 1227 1246 dfbdua = np.zeros_like(dfadua) 1228 1247 dfbdba = np.zeros_like(dfadba) 1229 dfadfl = 0.01230 dfbdfl = 0.01231 #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!1232 1248 SA = fas[0]+fas[1] 1233 1249 SB = fbs[0]+fbs[1] 1234 dFdfr[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs]*(SA[:,it,nxs]*(np.sum(dfadfr[:,:,it],axis=1))+ \ 1235 SB[:,it,nxs]*(np.sum(dfbdfr[:,:,it],axis=1)))*Mdata/len(SGMT) for it in range(nTwin)]),0,1) 1236 dFdx[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs,nxs]*(SA[:,it,nxs,nxs]*np.sum(dfadx[:,:,it],axis=1)+ \ 1237 SB[:,it,nxs,nxs]*np.sum(dfbdx[:,:,it],axis=1)) for it in range(nTwin)]),0,1) 1238 dFdui[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs]*(SA[:,it,nxs]*np.sum(dfadui[:,:,it],axis=1)+ \ 1239 SB[:,it,nxs]*np.sum(dfbdui[:,:,it],axis=1)) for it in range(nTwin)]),0,1) 1240 dFdua[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs,nxs]*(SA[:,it,nxs,nxs]*np.sum(dfadua[:,:,it],axis=1)+ \ 1241 SB[:,it,nxs,nxs]*np.sum(dfbdua[:,:,it],axis=1)) for it in range(nTwin)]),0,1) 1242 dFdtw[iBeg:iFin] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2 1243 1244 # dFdbab[iBeg:iFin] = 2.*(fas[0,nxs]*np.array([np.sum(dfadba.T*dBabdA,axis=0),np.sum(-dfadba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])+ \ 1245 # fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T 1250 dFdfr[iref] = [2.*TwMask[it]*(SA[it]*dfadfr[0,it]+SA[it]*dfadfr[1,it]+SB[it]*dfbdfr[0,it]+SB[it]*dfbdfr[1,it])*Mdata/len(Uniq[it]) for it in range(nTwin)] 1251 dFdx[iref] = [2.*TwMask[it]*(SA[it]*dfadx[it,0]+SA[it]*dfadx[it,1]+SB[it]*dfbdx[it,0]+SB[it]*dfbdx[it,1]) for it in range(nTwin)] 1252 dFdui[iref] = [2.*TwMask[it]*(SA[it]*dfadui[0,it]+SA[it]*dfadui[1,it]+SB[it]*dfbdui[0,it]+SB[it]*dfbdui[1,it]) for it in range(nTwin)] 1253 dFdua[iref] = [2.*TwMask[it]*(SA[it]*dfadua[it,0]+SA[it]*dfadua[it,1]+SB[it]*dfbdua[it,0]+SB[it]*dfbdua[it,1]) for it in range(nTwin)] 1254 if SGData['SGInv']: #centrosymmetric; B=0 1255 dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2 1256 else: 1257 dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2 1258 dFdbab[iref] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \ 1259 fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T 1246 1260 # GSASIIpath.IPyBreak() 1247 iBeg += blkSize 1248 print ' %d derivative time %.4f\r'%(nRef,time.time()-time0) 1249 #loop over atoms - each dict entry is list of derivatives for all the reflections 1261 print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0) 1262 #loop over atoms - each dict entry is list of derivatives for all the reflections 1250 1263 for i in range(len(Mdata)): #these all OK? 1251 1264 dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0) … … 1260 1273 dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0) 1261 1274 dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0) 1262 dFdvDict[phfx+'Flack'] = 4.*dFdfl.T 1275 dFdvDict[phfx+'BabA'] = dFdbab.T[0] 1276 dFdvDict[phfx+'BabU'] = dFdbab.T[1] 1277 for i in range(nTwin): 1278 dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i] 1279 return dFdvDict 1280 1281 def StructureFactorDervTw2(refDict,G,hfx,pfx,SGData,calcControls,parmDict): 1282 '''Compute structure factor derivatives on blocks of reflections - for twins only 1283 faster than StructureFactorDervTw 1284 input: 1285 1286 :param dict refDict: where 1287 'RefList' list where each ref = h,k,l,it,d,... 1288 'FF' dict of form factors - filled in below 1289 :param np.array G: reciprocal metric tensor 1290 :param str hfx: histogram id string 1291 :param str pfx: phase id string 1292 :param dict SGData: space group info. dictionary output from SpcGroup 1293 :param dict calcControls: 1294 :param dict parmDict: 1295 1296 :returns: dict dFdvDict: dictionary of derivatives 1297 ''' 1298 phfx = pfx.split(':')[0]+hfx 1299 ast = np.sqrt(np.diag(G)) 1300 Mast = twopisq*np.multiply.outer(ast,ast) 1301 SGMT = np.array([ops[0].T for ops in SGData['SGOps']]) 1302 SGT = np.array([ops[1] for ops in SGData['SGOps']]) 1303 FFtables = calcControls['FFtables'] 1304 BLtables = calcControls['BLtables'] 1305 TwDict = refDict.get('TwDict',{}) 1306 NTL = calcControls[phfx+'NTL'] 1307 NM = calcControls[phfx+'TwinNMN']+1 1308 TwinLaw = calcControls[phfx+'TwinLaw'] 1309 TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))]) 1310 TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1)) 1311 nTwin = len(TwinLaw) 1312 nRef = len(refDict['RefList']) 1313 Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict) 1314 mSize = len(Mdata) 1315 FF = np.zeros(len(Tdata)) 1316 if 'NC' in calcControls[hfx+'histType']: 1317 FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam']) 1318 elif 'X' in calcControls[hfx+'histType']: 1319 FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata]) 1320 FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata]) 1321 Uij = np.array(G2lat.U6toUij(Uijdata)) 1322 bij = Mast*Uij.T 1323 dFdvDict = {} 1324 dFdfr = np.zeros((nRef,nTwin,mSize)) 1325 dFdx = np.zeros((nRef,nTwin,mSize,3)) 1326 dFdui = np.zeros((nRef,nTwin,mSize)) 1327 dFdua = np.zeros((nRef,nTwin,mSize,6)) 1328 dFdbab = np.zeros((nRef,nTwin,2)) 1329 dFdtw = np.zeros((nRef,nTwin)) 1330 time0 = time.time() 1331 nref = len(refDict['RefList'])/100 1332 for iref,refl in enumerate(refDict['RefList']): 1333 if 'T' in calcControls[hfx+'histType']: 1334 FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12]) 1335 H = np.array(refl[:3]) 1336 H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(3,nTwins) or (3) 1337 TwMask = np.any(H,axis=-1) 1338 if iref in TwDict: 1339 for i in TwDict[iref]: 1340 for n in range(NTL): 1341 H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM]) 1342 TwMask = np.any(H,axis=-1) 1343 SQ = 1./(2.*refl[4])**2 # or (sin(theta)/lambda)**2 1344 SQfactor = 8.0*SQ*np.pi**2 1345 dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor) 1346 Bab = parmDict[phfx+'BabA']*dBabdA 1347 Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) 1348 FF = refDict['FF']['FF'][iref].T[Tindx].T 1349 Uniq = np.inner(H,SGMT) # array(nSGOp,3) or (nTwin,nSGOp,3) 1350 Phi = np.inner(H,SGT) 1351 phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T 1352 sinp = np.sin(phase) 1353 cosp = np.cos(phase) 1354 occ = Mdata*Fdata/len(SGT) 1355 biso = -SQfactor*Uisodata[:,nxs] 1356 Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1) 1357 HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1) 1358 Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))]) 1359 Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6))) 1360 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 1361 Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ 1362 fot = (FF+FP-Bab)*Tcorr 1363 fotp = FPP*Tcorr 1364 fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-FPP*sinp*Tcorr]) 1365 fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,FPP*cosp*Tcorr]) 1366 # GSASIIpath.IPyBreak() 1367 fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl array(2,nTwins) 1368 fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl 1369 if SGData['SGInv']: #centrosymmetric; B=0 1370 fbs[0] *= 0. 1371 fas[1] *= 0. 1372 fax = np.array([-fot*sinp,-fotp*cosp]) #positions array(2,ntwi,nEqv,nAtoms) 1373 fbx = np.array([fot*cosp,-fotp*sinp]) 1374 #sum below is over Uniq 1375 dfadfr = np.sum(fa/occ,axis=-2) #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem 1376 dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1) 1377 dfadui = np.sum(-SQfactor*fa,axis=-2) 1378 dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1379 dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1380 # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6) 1381 if not SGData['SGInv']: 1382 dfbdfr = np.sum(fb/occ,axis=-2) #Fdata != 0 avoids /0. problem 1383 dfadba /= 2. 1384 dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2. 1385 dfbdui = np.sum(-SQfactor*fb,axis=-2) 1386 dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)]) 1387 dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)]) 1388 else: 1389 dfbdfr = np.zeros_like(dfadfr) 1390 dfbdx = np.zeros_like(dfadx) 1391 dfbdui = np.zeros_like(dfadui) 1392 dfbdua = np.zeros_like(dfadua) 1393 dfbdba = np.zeros_like(dfadba) 1394 SA = fas[0]+fas[1] 1395 SB = fbs[0]+fbs[1] 1396 dFdfr[iref] = [2.*TwMask[it]*(SA[it]*dfadfr[0,it]+SA[it]*dfadfr[1,it]+SB[it]*dfbdfr[0,it]+SB[it]*dfbdfr[1,it])*Mdata/len(Uniq[it]) for it in range(nTwin)] 1397 dFdx[iref] = [2.*TwMask[it]*(SA[it]*dfadx[it,0]+SA[it]*dfadx[it,1]+SB[it]*dfbdx[it,0]+SB[it]*dfbdx[it,1]) for it in range(nTwin)] 1398 dFdui[iref] = [2.*TwMask[it]*(SA[it]*dfadui[0,it]+SA[it]*dfadui[1,it]+SB[it]*dfbdui[0,it]+SB[it]*dfbdui[1,it]) for it in range(nTwin)] 1399 dFdua[iref] = [2.*TwMask[it]*(SA[it]*dfadua[it,0]+SA[it]*dfadua[it,1]+SB[it]*dfbdua[it,0]+SB[it]*dfbdua[it,1]) for it in range(nTwin)] 1400 if SGData['SGInv']: #centrosymmetric; B=0 1401 dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2 1402 else: 1403 dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2 1404 dFdbab[iref] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \ 1405 fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T 1406 # GSASIIpath.IPyBreak() 1407 print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0) 1408 #loop over atoms - each dict entry is list of derivatives for all the reflections 1409 for i in range(len(Mdata)): #these all OK? 1410 dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0) 1411 dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0) 1412 dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0) 1413 dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0) 1414 dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0) 1415 dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0) 1416 dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0) 1417 dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0) 1418 dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0) 1419 dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0) 1420 dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0) 1263 1421 dFdvDict[phfx+'BabA'] = dFdbab.T[0] 1264 1422 dFdvDict[phfx+'BabU'] = dFdbab.T[1] … … 1317 1475 Uij = np.array(G2lat.U6toUij(Uijdata)).T 1318 1476 bij = Mast*Uij 1319 blkSize = 100#no. of reflections in a block1477 blkSize = 32 #no. of reflections in a block 1320 1478 nRef = refDict['RefList'].shape[0] 1321 1479 if not len(refDict['FF']): … … 1340 1498 H = refl.T[:4] #array(blkSize,4) 1341 1499 HP = H[:3]+modQ[:,nxs]*H[3:] #projected hklm to hkl 1342 # H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(blkSize,nTwins,4) or (blkSize,4)1343 # TwMask = np.any(H ,axis=-1)1500 # HP = np.squeeze(np.inner(HP.T,TwinLaw)) #maybe array(blkSize,nTwins,4) or (blkSize,4) 1501 # TwMask = np.any(HP,axis=-1) 1344 1502 # if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i] 1345 1503 # for ir in range(blkSize): … … 1348 1506 # for i in TwDict[iref]: 1349 1507 # for n in range(NTL): 1350 # H [ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])1351 # TwMask = np.any(H ,axis=-1)1508 # HP[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM]) 1509 # TwMask = np.any(HP,axis=-1) 1352 1510 SQ = 1./(2.*refl.T[5])**2 #array(blkSize) 1353 1511 SQfactor = 4.0*SQ*twopisq #ditto prev. … … 1407 1565 1408 1566 def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict): 1409 'Needs a doc string - no twins' 1567 ''' 1568 Compute super structure factor derivatives for all h,k,l,m for phase - no twins 1569 input: 1570 1571 :param dict refDict: where 1572 'RefList' list where each ref = h,k,l,m,it,d,... 1573 'FF' dict of form factors - filled in below 1574 :param int im: = 1 (could be eliminated) 1575 :param np.array G: reciprocal metric tensor 1576 :param str hfx: histogram id string 1577 :param str pfx: phase id string 1578 :param dict SGData: space group info. dictionary output from SpcGroup 1579 :param dict SSGData: super space group info. 1580 :param dict calcControls: 1581 :param dict ParmDict: 1582 1583 :returns: dict dFdvDict: dictionary of derivatives 1584 ''' 1410 1585 phfx = pfx.split(':')[0]+hfx 1411 1586 ast = np.sqrt(np.diag(G)) … … 1568 1743 if not iref%100 : 1569 1744 print ' %d derivative time %.4f\r'%(iref,time.time()-time0), 1745 for i in range(len(Mdata)): #loop over atoms 1746 dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i] 1747 dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i] 1748 dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i] 1749 dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i] 1750 dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i] 1751 dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i] 1752 dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i] 1753 dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i] 1754 dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i] 1755 dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i] 1756 dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i] 1757 for j in range(FSSdata.shape[1]): #loop over waves Fzero & Fwid? 1758 dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i] 1759 dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i] 1760 nx = 0 1761 if waveTypes[i] in ['Block','ZigZag']: 1762 nx = 1 1763 dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i] #ZigZag/Block waves (if any) 1764 dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i] 1765 dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i] 1766 dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i] 1767 dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i] 1768 for j in range(XSSdata.shape[1]-nx): #loop over waves 1769 dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i] 1770 dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i] 1771 dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i] 1772 dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i] 1773 dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i] 1774 dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i] 1775 for j in range(USSdata.shape[1]): #loop over waves 1776 dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i] 1777 dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i] 1778 dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i] 1779 dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i] 1780 dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i] 1781 dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i] 1782 dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i] 1783 dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i] 1784 dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i] 1785 dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i] 1786 dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i] 1787 dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i] 1788 1789 # GSASIIpath.IPyBreak() 1790 dFdvDict[phfx+'BabA'] = dFdbab.T[0] 1791 dFdvDict[phfx+'BabU'] = dFdbab.T[1] 1792 return dFdvDict 1793 1794 def SStructureFactorDerv2(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict): 1795 'Needs a doc string - no twins' 1796 phfx = pfx.split(':')[0]+hfx 1797 ast = np.sqrt(np.diag(G)) 1798 Mast = twopisq*np.multiply.outer(ast,ast) 1799 SGInv = SGData['SGInv'] 1800 SGMT = np.array([ops[0].T for ops in SGData['SGOps']]) 1801 SGT = np.array([ops[1] for ops in SGData['SGOps']]) 1802 SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']]) 1803 SSGT = np.array([ops[1] for ops in SSGData['SSGOps']]) 1804 FFtables = calcControls['FFtables'] 1805 BLtables = calcControls['BLtables'] 1806 nRef = len(refDict['RefList']) 1807 Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict) 1808 mSize = len(Mdata) #no. atoms 1809 waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict) 1810 ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast) 1811 waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast) 1812 modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']]) 1813 FF = np.zeros(len(Tdata)) 1814 if 'NC' in calcControls[hfx+'histType']: 1815 FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam']) 1816 elif 'X' in calcControls[hfx+'histType']: 1817 FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata]) 1818 FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata]) 1819 Uij = np.array(G2lat.U6toUij(Uijdata)).T 1820 bij = Mast*Uij 1821 if not len(refDict['FF']): 1822 if 'N' in calcControls[hfx+'histType']: 1823 dat = G2el.getBLvalues(BLtables) #will need wave here for anom. neutron b's 1824 else: 1825 dat = G2el.getFFvalues(FFtables,0.) 1826 refDict['FF']['El'] = dat.keys() 1827 refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat))) 1828 dFdvDict = {} 1829 dFdfr = np.zeros((nRef,mSize)) 1830 dFdx = np.zeros((nRef,mSize,3)) 1831 dFdui = np.zeros((nRef,mSize)) 1832 dFdua = np.zeros((nRef,mSize,6)) 1833 dFdbab = np.zeros((nRef,2)) 1834 dFdfl = np.zeros((nRef)) 1835 dFdtw = np.zeros((nRef)) 1836 dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2)) 1837 dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6)) 1838 dFdGz = np.zeros((nRef,mSize,5)) 1839 dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12)) 1840 Flack = 1.0 1841 if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict: 1842 Flack = 1.-2.*parmDict[phfx+'Flack'] 1843 time0 = time.time() 1844 iBeg = 0 1845 blkSize = 4 #no. of reflections in a block - optimized for speed 1846 while iBeg < nRef: 1847 iFin = min(iBeg+blkSize,nRef) 1848 refl = refDict['RefList'][iBeg:iFin] #array(blkSize,nItems) 1849 H = refl.T[:4] 1850 HP = H[:3].T+modQ*H.T[:,3:] #projected hklm to hkl 1851 SQ = 1./(2.*refl.T[4+im])**2 # or (sin(theta)/lambda)**2 1852 SQfactor = 8.0*SQ*np.pi**2 1853 if 'T' in calcControls[hfx+'histType']: 1854 if 'P' in calcControls[hfx+'histType']: 1855 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[15]) 1856 else: 1857 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[13]) 1858 FP = np.repeat(FP.T,len(SGT),axis=0) 1859 FPP = np.repeat(FPP.T,len(SGT),axis=0) 1860 dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor) 1861 Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)) 1862 Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) 1863 FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0) 1864 Uniq = np.inner(H.T,SSGMT) 1865 Phi = np.inner(H.T,SSGT) 1866 UniqP = np.inner(HP,SGMT) 1867 if SGInv: #if centro - expand HKL sets 1868 Uniq = np.hstack((Uniq,-Uniq)) 1869 Phi = np.hstack((Phi,-Phi)) 1870 UniqP = np.hstack((UniqP,-UniqP)) 1871 FF = np.vstack((FF,FF)) 1872 Bab = np.concatenate((Bab,Bab)) 1873 phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata+Xdata).T)+Phi[:,:,nxs]) 1874 sinp = np.sin(phase) 1875 cosp = np.cos(phase) 1876 occ = Mdata*Fdata/Uniq.shape[1] 1877 biso = -SQfactor*Uisodata[:,nxs] 1878 Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1],axis=1).T #ops x atoms 1879 HbH = -np.sum(UniqP[:,:,nxs,:3]*np.inner(UniqP[:,:,:3],bij),axis=-1) #ops x atoms 1880 Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in np.reshape(UniqP,(-1,3))]) #atoms x 3x3 1881 Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(iFin-iBeg,-1,6)) #atoms x 6 1882 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) #ops x atoms 1883 # GSASIIpath.IPyBreak() 1884 Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0] #ops x atoms 1885 fot = np.reshape(FF+FP[nxs,:]-Bab[:,nxs],cosp.shape)*Tcorr #ops x atoms 1886 fotp = FPP*Tcorr #ops x atoms 1887 GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms 1888 dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv2(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt) 1889 # GfpuA is 2 x ops x atoms 1890 # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts 1891 fa = np.array([fot*cosp,-Flack*FPP*sinp*Tcorr]) # array(2,nEqv,nAtoms) 1892 fb = np.array([fot*sinp,Flack*FPP*cosp*Tcorr]) #or array(2,nEqv,nAtoms) 1893 fag = fa*GfpuA[0]-fb*GfpuA[1] 1894 fbg = fb*GfpuA[0]+fa*GfpuA[1] 1895 1896 fas = np.sum(np.sum(fag,axis=-1),axis=-1) # 2 x refBlk 1897 fbs = np.sum(np.sum(fbg,axis=-1),axis=-1) 1898 fax = np.array([-fot*sinp,-fotp*cosp]) #positions; 2 x ops x atoms 1899 fbx = np.array([fot*cosp,-fotp*sinp]) 1900 fax = fax*GfpuA[0]-fbx*GfpuA[1] 1901 fbx = fbx*GfpuA[0]+fax*GfpuA[1] 1902 #sum below is over Uniq 1903 dfadfr = np.sum(fag/occ,axis=-2) #Fdata != 0 ever avoids /0. problem 1904 dfbdfr = np.sum(fbg/occ,axis=-2) #Fdata != 0 avoids /0. problem 1905 dfadba = np.sum(-cosp*Tcorr,axis=-2) 1906 dfbdba = np.sum(-sinp*Tcorr,axis=-2) 1907 dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fag,axis=-2) 1908 dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fbg,axis=-2) 1909 dfadx = np.sum(twopi*Uniq[nxs,:,:,nxs,:3]*fax[:,:,:,:,nxs],axis=-3) #2 refBlk x x nAtom x 3xyz; sum nOps 1910 dfbdx = np.sum(twopi*Uniq[nxs,:,:,nxs,:3]*fbx[:,:,:,:,nxs],axis=-3) #2 refBlk x x nAtom x 3xyz; sum nOps 1911 dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fag[:,:,:,:,nxs],axis=-3) #2 x nAtom x 6Uij; sum nOps 1912 dfbdua = np.sum(-Hij[nxs,:,:,nxs,:]*fbg[:,:,:,:,nxs],axis=-3) #2 x nAtom x 6Uij; sum nOps 1913 # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps 1914 dfadGf = np.sum(fa[:,:,:,:,nxs,nxs]*dGdf[0][nxs,:,nxs,:,:,:]-fb[:,:,:,:,nxs,nxs]*dGdf[1][nxs,:,nxs,:,:,:],axis=2) 1915 dfbdGf = np.sum(fb[:,:,:,:,nxs,nxs]*dGdf[0][nxs,:,nxs,:,:,:]+fa[:,:,:,:,nxs,nxs]*dGdf[1][nxs,:,nxs,:,:,:],axis=2) 1916 dfadGx = np.sum(fa[:,:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:,:]-fb[:,:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:,:],axis=2) 1917 dfbdGx = np.sum(fb[:,:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:,:]+fa[:,:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:,:],axis=2) 1918 dfadGz = np.sum(fa[:,:,:,:,nxs]*dGdz[0][nxs,:,:,:,:]-fb[:,:,:,:,nxs]*dGdz[1][nxs,:,:,:,:],axis=2) 1919 dfbdGz = np.sum(fb[:,:,:,:,nxs]*dGdz[0][nxs,:,:,:,:]+fa[:,:,:,:,nxs]*dGdz[1][nxs,:,:,:,:],axis=2) 1920 dfadGu = np.sum(fa[:,:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=2) 1921 dfbdGu = np.sum(fb[:,:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=2) 1922 if not SGData['SGInv']: #Flack derivative 1923 dfadfl = np.sum(np.sum(-FPP*Tcorr*sinp,axis=-1),axis=-1) 1924 dfbdfl = np.sum(np.sum(FPP*Tcorr*cosp,axis=-1),axis=-1) 1925 else: 1926 dfadfl = 1.0 1927 dfbdfl = 1.0 1928 #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3! 1929 SA = fas[0]+fas[1] #float = A+A' (might be array[nTwin]) 1930 SB = fbs[0]+fbs[1] #float = B+B' (might be array[nTwin]) 1931 if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro? 1932 dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl #array(nRef,) 1933 dFdfr[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadfr[0]+fas[1,:,nxs]*dfadfr[1])*Mdata/len(Uniq)+ \ 1934 2.*(fbs[0,:,nxs]*dfbdfr[0]-fbs[1,:,nxs]*dfbdfr[1])*Mdata/len(Uniq) 1935 dFdx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadx[0]+fas[1,:,nxs,nxs]*dfadx[1])+ \ 1936 2.*(fbs[0,:,nxs,nxs]*dfbdx[0]+fbs[1,:,nxs,nxs]*dfbdx[1]) 1937 dFdui[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadui[0]+fas[1,:,nxs]*dfadui[1])+ \ 1938 2.*(fbs[0,:,nxs]*dfbdui[0]-fbs[1,:,nxs]*dfbdui[1]) 1939 dFdua[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadua[0]+fas[1,:,nxs,nxs]*dfadua[1])+ \ 1940 2.*(fbs[0,:,nxs,nxs]*dfbdua[0]+fbs[1,:,nxs,nxs]*dfbdua[1]) 1941 dFdGf[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs,nxs]*dfadGf[0]+fas[1,:,nxs,nxs,nxs]*dfadGf[1])+ \ 1942 2.*(fbs[0,:,nxs,nxs,nxs]*dfbdGf[0]+fbs[1,:,nxs,nxs,nxs]*dfbdGf[1]) 1943 dFdGx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs,nxs]*dfadGx[0]+fas[1,:,nxs,nxs,nxs]*dfadGx[1])+ \ 1944 2.*(fbs[0,:,nxs,nxs,nxs]*dfbdGx[0]-fbs[1,:,nxs,nxs,nxs]*dfbdGx[1]) 1945 dFdGz[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadGz[0]+fas[1,:,nxs,nxs]*dfadGz[1])+ \ 1946 2.*(fbs[0,:,nxs,nxs]*dfbdGz[0]+fbs[1,:,nxs,nxs]*dfbdGz[1]) 1947 dFdGu[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs,nxs]*dfadGu[0]+fas[1,:,nxs,nxs,nxs]*dfadGu[1])+ \ 1948 2.*(fbs[0,:,nxs,nxs,nxs]*dfbdGu[0]+fbs[1,:,nxs,nxs,nxs]*dfbdGu[1]) 1949 else: #OK, I think 1950 dFdfr[iBeg:iFin] = 2.*(SA[:,nxs]*(dfadfr[0]+dfadfr[1])+SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(Uniq) #array(nRef,nAtom) 1951 dFdx[iBeg:iFin] = 2.*(SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1])) #array(nRef,nAtom,3) 1952 dFdui[iBeg:iFin] = 2.*(SA[:,nxs]*(dfadui[0]+dfadui[1])+SB[:,nxs]*(dfbdui[0]+dfbdui[1])) #array(nRef,nAtom) 1953 dFdua[iBeg:iFin] = 2.*(SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1])) #array(nRef,nAtom,6) 1954 dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl #array(nRef,) 1955 1956 dFdGf[iBeg:iFin] = 2.*(SA[:,nxs,nxs,nxs]*dfadGf[0]+SB[:,nxs,nxs,nxs]*dfbdGf[1]) #array(nRef,natom,nwave,2) 1957 dFdGx[iBeg:iFin] = 2.*(SA[:,nxs,nxs,nxs]*dfadGx[0]+SB[:,nxs,nxs,nxs]*dfbdGx[1]) #array(nRef,natom,nwave,6) 1958 dFdGz[iBeg:iFin] = 2.*(SA[:,nxs,nxs]*dfadGz[0]+SB[:,nxs,nxs]*dfbdGz[1]) #array(nRef,natom,5) 1959 dFdGu[iBeg:iFin] = 2.*(SA[:,nxs,nxs,nxs]*dfadGu[0]+SB[:,nxs,nxs,nxs]*dfbdGu[1]) #array(nRef,natom,nwave,12) 1960 # GSASIIpath.IPyBreak() 1961 # dFdbab[iBeg:iFin] = 2.*fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \ 1962 # 2.*fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T 1963 #loop over atoms - each dict entry is list of derivatives for all the reflections 1964 print ' %d derivative time %.4f\r'%(iBeg,time.time()-time0), 1965 iBeg += blkSize 1570 1966 for i in range(len(Mdata)): #loop over atoms 1571 1967 dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i] … … 3162 3558 else: 3163 3559 if len(TwinLaw) > 1: 3164 dFdvDict = StructureFactorDerv (refDict,G,hfx,pfx,SGData,calcControls,parmDict)3165 else: 3560 dFdvDict = StructureFactorDervTw2(refDict,G,hfx,pfx,SGData,calcControls,parmDict) 3561 else: #correct!! 3166 3562 dFdvDict = StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict) 3167 3563 ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase) … … 3362 3758 def errRefine(values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup,dlg=None): 3363 3759 '''Computes the point-by-point discrepancies between every data point in every histogram 3364 and the observed value 3760 and the observed value. Used in the Jacobian, Hessian & numeric least-squares to compute function 3365 3761 3366 3762 :returns: an np array of differences between observed and computed diffraction values.
Note: See TracChangeset
for help on using the changeset viewer.