Changeset 1980
- Timestamp:
- Sep 28, 2015 3:54:31 PM (7 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIImath.py
r1978 r1980 926 926 ################################################################################ 927 927 928 def Modulation (waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast):928 def Modulation2(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast): 929 929 930 930 nxs = np.newaxis … … 969 969 UmodB = Bu[:,:,:,:,nxs]*np.cos(twopi*tauU[nxs,:,nxs,nxs,:]) #ditto 970 970 Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3) #atoms x 3x3 x 32; sum waves 971 HbH = np.exp( np.swapaxes(np.array([-np.inner(h,np.inner(Umod,h)) for h in H[:,:3]]),1,2)).T #Overhauser corr.?971 HbH = np.exp(-np.sum(H[:,:,nxs,nxs,:3]*np.inner(H[:,:,:3],Umod),axis=-1)) # refBlk x ops x atoms x 32 add Overhauser corr.? 972 972 else: 973 973 HbH = 1.0 974 D = H[:, 3][:,nxs]*glTau[nxs,:] #m*tau;ops X 32975 HdotX = np.inner( np.swapaxes(Xmod,1,2),H[:,:3])+D.T[nxs,:,:] #atoms X 32 X ops976 sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt [nxs,:,nxs],axis=1) #atoms X ops; sum for G-L integration977 cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt [nxs,:,nxs],axis=1) #ditto974 D = H[:,:,3][:,:,nxs]*glTau[nxs,nxs,:] #m*tau;refBlk x ops X 32 975 HdotX = np.inner(H[:,:,:3],np.swapaxes(Xmod,1,2))+D[:,:,nxs,:] #refBlk X ops x atoms X 32 976 sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt,axis=-1) #refBlk X ops x atoms; sum for G-L integration 977 cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt,axis=-1) #ditto 978 978 # GSASIIpath.IPyBreak() 979 return cosHA .T,sinHA.T #ops X atoms979 return cosHA,sinHA #each refBlk x ops X atoms 980 980 981 981 Ax = np.array(XSSdata[:3]).T #atoms x waves x sin pos mods … … 985 985 Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T #atoms x waves x sin Uij mods 986 986 Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T #...cos Uij mods 987 GpA = np.array(expModInt(SSUniq,Af,Bf,Ax,Bx,Au,Bu)) 988 return GpA # SGops x atoms 989 987 GpA = np.array(expModInt(Uniq,Af,Bf,Ax,Bx,Au,Bu)) 988 return GpA # 2 x refBlk x SGops x atoms 989 990 #def Modulation(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast): 991 # 992 # nxs = np.newaxis 993 # glTau,glWt = pwd.pygauleg(0.,1.,32) #get Gauss-Legendre intervals & weights 994 # 995 # def expModInt(H,Af,Bf,Ax,Bx,Au,Bu): 996 # ''' 997 # H: array ops X hklt 998 # Ax: array atoms X waves X xyz 999 # Bx: array atoms X waves X xyz 1000 # S: array ops 1001 # ''' 1002 # if 'Fourier' in waveTypes: 1003 # nf = 0 1004 # nx = 0 1005 # XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,32)) 1006 # FmodZ = np.zeros((Af.shape[0],Af.shape[1],32)) 1007 # if 'Crenel' in waveTypes: 1008 # nC = np.where('Crenel' in waveTypes) 1009 # nf = 1 1010 # #FmodZ = 0 replace 1011 # else: 1012 # nx = 1 1013 # if 'Sawtooth' in waveTypes: 1014 # nS = np.where('Sawtooth' in waveTypes) 1015 # #XmodZ = 0 replace 1016 # if Af.shape[1]: 1017 # tauF = np.arange(1.,Af.shape[1]+1-nf)[:,nxs]*glTau #Fwaves x 32 1018 # FmodA = Af[:,nf:,nxs]*np.sin(twopi*tauF[nxs,:,:]) #atoms X Fwaves X 32 1019 # FmodB = Bf[:,nf:,nxs]*np.cos(twopi*tauF[nxs,:,:]) 1020 # Fmod = np.sum(FmodA+FmodB+FmodC,axis=1) #atoms X 32; sum waves 1021 # else: 1022 # Fmod = 1.0 1023 # if Ax.shape[1]: 1024 # tauX = np.arange(1.,Ax.shape[1]+1-nx)[:,nxs]*glTau #Xwaves x 32 1025 # XmodA = Ax[:,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X pos X 32 1026 # XmodB = Bx[:,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto 1027 # Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1) #atoms X pos X 32; sum waves 1028 # if Au.shape[1]: 1029 # tauU = np.arange(1.,Au.shape[1]+1)[:,nxs]*glTau #Uwaves x 32 1030 # UmodA = Au[:,:,:,:,nxs]*np.sin(twopi*tauU[nxs,:,nxs,nxs,:]) #atoms x waves x 3x3 x 32 1031 # UmodB = Bu[:,:,:,:,nxs]*np.cos(twopi*tauU[nxs,:,nxs,nxs,:]) #ditto 1032 # Umod = np.swapaxes(np.sum(UmodA+UmodB,axis=1),1,3) #atoms x 3x3 x 32; sum waves 1033 # HbH = np.exp(np.swapaxes(np.array([-np.inner(h,np.inner(Umod,h)) for h in H[:,:3]]),1,2)).T #Overhauser corr.? 1034 # else: 1035 # HbH = 1.0 1036 # D = H[:,3][:,nxs]*glTau[nxs,:] #m*tau; ops X 32 1037 # HdotX = np.inner(np.swapaxes(Xmod,1,2),H[:,:3])+D.T[nxs,:,:] #atoms X 32 X ops 1038 # sinHA = np.sum(Fmod*HbH*np.sin(twopi*HdotX)*glWt[nxs,:,nxs],axis=1) #atoms X ops; sum for G-L integration 1039 # cosHA = np.sum(Fmod*HbH*np.cos(twopi*HdotX)*glWt[nxs,:,nxs],axis=1) #ditto 1040 ## GSASIIpath.IPyBreak() 1041 # return cosHA.T,sinHA.T #ops X atoms 1042 # 1043 # Ax = np.array(XSSdata[:3]).T #atoms x waves x sin pos mods 1044 # Bx = np.array(XSSdata[3:]).T #...cos pos mods 1045 # Af = np.array(FSSdata[0]).T #sin frac mods x waves x atoms 1046 # Bf = np.array(FSSdata[1]).T #cos frac mods... 1047 # Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T #atoms x waves x sin Uij mods 1048 # Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T #...cos Uij mods 1049 # GpA = np.array(expModInt(SSUniq,Af,Bf,Ax,Bx,Au,Bu)) 1050 # return GpA # 2 x SGops x atoms 1051 # 990 1052 def ModulationDerv(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast): 991 1053 -
trunk/GSASIIstrMath.py
r1979 r1980 953 953 return dFdvDict 954 954 955 def SStructureFactor (refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):955 def SStructureFactor2(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict): 956 956 ''' 957 957 Compute super structure factors for all h,k,l,m for phase … … 968 968 :param dict ParmDict: 969 969 970 ''' 970 ''' 971 nxs = np.newaxis 971 972 phfx = pfx.split(':')[0]+hfx 972 973 ast = np.sqrt(np.diag(G)) … … 1013 1014 refDict['FF']['FF'] = np.ones((nRef,len(dat))) 1014 1015 for iref,ref in enumerate(refDict['RefList']): 1015 SQ = 1./(2.*ref[ 4])**21016 SQ = 1./(2.*ref[5])**2 1016 1017 dat = G2el.getFFvalues(FFtables,SQ) 1017 1018 refDict['FF']['FF'][iref] *= dat.values() 1018 1019 time0 = time.time() 1019 nref = len(refDict['RefList'])/1001020 1020 #reflection processing begins here - big arrays! 1021 1021 iBeg = 0 … … 1024 1024 refl = refDict['RefList'][iBeg:iFin] #array(blkSize,nItems) 1025 1025 H = refl.T[:4] #array(blkSize,4) 1026 H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(blkSize,nTwins,3) or (blkSize,3)1027 TwMask = np.any(H,axis=-1)1028 if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]1029 for ir in range(blkSize):1030 iref = ir+iBeg1031 if iref in TwDict:1032 for i in TwDict[iref]:1033 for n in range(NTL):1034 H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])1035 TwMask = np.any(H,axis=-1)1026 # H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(blkSize,nTwins,4) or (blkSize,4) 1027 # TwMask = np.any(H,axis=-1) 1028 # if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i] 1029 # for ir in range(blkSize): 1030 # iref = ir+iBeg 1031 # if iref in TwDict: 1032 # for i in TwDict[iref]: 1033 # for n in range(NTL): 1034 # H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM]) 1035 # TwMask = np.any(H,axis=-1) 1036 1036 SQ = 1./(2.*refl.T[5])**2 #array(blkSize) 1037 1037 SQfactor = 4.0*SQ*twopisq #ditto prev. 1038 Uniq = np.inner(H.T,SSGMT) 1039 Phi = np.inner(H.T,SSGT) 1040 if SGInv: #if centro - expand HKL sets 1041 Uniq = np.hstack((Uniq,-Uniq)) 1042 Phi = np.hstack((Phi,-Phi)) 1038 1043 if 'T' in calcControls[hfx+'histType']: 1039 1044 if 'P' in calcControls[hfx+'histType']: … … 1041 1046 else: 1042 1047 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12]) 1043 FP = np.repeat(FP.T, len(SGT)*len(TwinLaw),axis=0)1044 FPP = np.repeat(FPP.T, len(SGT)*len(TwinLaw),axis=0)1045 Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor), len(SGT)*len(TwinLaw))1048 FP = np.repeat(FP.T,Uniq.shape[1]*len(TwinLaw),axis=0) 1049 FPP = np.repeat(FPP.T,Uniq.shape[1]*len(TwinLaw),axis=0) 1050 Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),Uniq.shape[1]*len(TwinLaw)) 1046 1051 Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) 1047 FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0) 1048 Uniq = np.inner(H,SSGMT) 1049 Phi = np.inner(H,SSGT) 1050 if SGInv: #if centro - expand HKL sets 1051 Uniq = np.vstack((Uniq,-Uniq)) 1052 Phi = np.hstack((Phi,-Phi)) 1053 # GSASIIpath.IPyBreak() 1054 GfpuA = G2mth.Modulation(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast) #2 x sym X atoms 1055 phase = twopi*(np.inner(Uniq[:,:3],(dXdata.T+Xdata.T))+Phi[:,np.newaxis]) 1052 FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,Uniq.shape[1]*len(TwinLaw),axis=0) 1053 phase = twopi*(np.inner(Uniq[:,:,:3],(dXdata.T+Xdata.T))+Phi[:,:,nxs]) 1056 1054 sinp = np.sin(phase) 1057 1055 cosp = np.cos(phase) 1058 biso = -SQfactor*Uisodata 1059 Tiso = np. where(biso<1.,np.exp(biso),1.0)1060 HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq[:,:3]])1056 biso = -SQfactor*Uisodata[:,nxs] 1057 Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[1]*len(TwinLaw),axis=1).T 1058 HbH = -np.sum(Uniq[:,:,nxs,:3]*np.inner(Uniq[:,:,:3],bij),axis=-1) 1061 1059 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 1062 Tcorr = Tiso*Tuij*Mdata*Fdata/len(Uniq) 1063 fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr]) #2 x sym x atoms 1064 fb = np.array([FPP*cosp*Tcorr,(FF+FP-Bab)*sinp*Tcorr]) #swapped around - better 1060 Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[1] 1061 if 'T' in calcControls[hfx+'histType']: 1062 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr]) 1063 fb = np.array([np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr]) 1064 else: 1065 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) 1066 fb = np.array([Flack*FPP*cosp*Tcorr,np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr]) 1067 GfpuA = G2mth.Modulation2(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast) #2 x refBlk x sym X atoms 1065 1068 fa *= GfpuA 1066 1069 fb *= GfpuA 1067 fas = np.sum(np.sum(fa,axis=1),axis=1) 1068 fbs = np.sum(np.sum(fb,axis=1),axis=1) 1069 fasq = fas**2 1070 fbsq = fbs**2 #imaginary 1071 refl[9+im] = np.sum(fasq)+np.sum(fbsq) 1072 refl[7+im] = np.sum(fasq)+np.sum(fbsq) 1073 refl[10+im] = atan2d(fbs[0],fas[0]) 1074 if not iref%nref: print 'ref no. %d time %.4f\r'%(iref,time.time()-time0), 1075 print '\nref no. %d time %.4f\r'%(iref,time.time()-time0) 1076 1070 # GSASIIpath.IPyBreak() 1071 fas = np.sum(np.sum(fa,axis=-1),axis=-1) 1072 fbs = np.sum(np.sum(fb,axis=-1),axis=-1) 1073 if 'P' in calcControls[hfx+'histType']: 1074 refl.T[10] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) 1075 refl.T[11] = atan2d(fbs[0],fas[0]) #ignore f' & f" 1076 else: 1077 if len(TwinLaw) > 1: 1078 refl.T[10] = np.sum(fas[:,:,0]**2,axis=0)+np.sum(fbs[:,:,0]**2,axis=0) #FcT from primary twin element 1079 refl.T[8] = np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fas,axis=0)**2,axis=-1)+ \ 1080 np.sum(TwinFr*np.sum(TwMask[np.newaxis,:,:]*fbs,axis=0)**2,axis=-1) #Fc sum over twins 1081 refl.T[11] = atan2d(fbs[0].T[0],fas[0].T[0]) #ignore f' & f" 1082 else: 1083 refl.T[10] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2 1084 refl.T[8] = np.copy(refl.T[10]) 1085 refl.T[11] = atan2d(fbs[0],fas[0]) #ignore f' & f" 1086 iBeg += blkSize 1087 print 'nRef %d time %.4f\r'%(nRef,time.time()-time0) 1088 1089 #def SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict): 1090 # ''' 1091 # Compute super structure factors for all h,k,l,m for phase 1092 # puts the result, F^2, in each ref[8+im] in refList 1093 # input: 1094 # 1095 # :param dict refDict: where 1096 # 'RefList' list where each ref = h,k,l,m,d,... 1097 # 'FF' dict of form factors - filed in below 1098 # :param np.array G: reciprocal metric tensor 1099 # :param str pfx: phase id string 1100 # :param dict SGData: space group info. dictionary output from SpcGroup 1101 # :param dict calcControls: 1102 # :param dict ParmDict: 1103 # 1104 # ''' 1105 # phfx = pfx.split(':')[0]+hfx 1106 # ast = np.sqrt(np.diag(G)) 1107 # Mast = twopisq*np.multiply.outer(ast,ast) 1108 # 1109 # SGInv = SGData['SGInv'] 1110 # SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']]) 1111 # SSGT = np.array([ops[1] for ops in SSGData['SSGOps']]) 1112 # FFtables = calcControls['FFtables'] 1113 # BLtables = calcControls['BLtables'] 1114 # Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict) 1115 # waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict) 1116 # FF = np.zeros(len(Tdata)) 1117 # if 'NC' in calcControls[hfx+'histType']: 1118 # FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam']) 1119 # else: 1120 # FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata]) 1121 # FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata]) 1122 # Uij = np.array(G2lat.U6toUij(Uijdata)) 1123 # bij = Mast*Uij.T 1124 # if not len(refDict['FF']): 1125 # if 'N' in calcControls[hfx+'histType']: 1126 # dat = G2el.getBLvalues(BLtables) #will need wave here for anom. neutron b's 1127 # else: 1128 # dat = G2el.getFFvalues(FFtables,0.) 1129 # refDict['FF']['El'] = dat.keys() 1130 # refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat))) 1131 # time0 = time.time() 1132 # nref = len(refDict['RefList'])/100 1133 # for iref,refl in enumerate(refDict['RefList']): 1134 # if 'NT' in calcControls[hfx+'histType']: 1135 # FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14+im]) 1136 # fbs = np.array([0,0]) 1137 # H = refl[:4] 1138 # SQ = 1./(2.*refl[4+im])**2 1139 # SQfactor = 4.0*SQ*twopisq 1140 # Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor) 1141 # if not np.any(refDict['FF']['FF'][iref]): #no form factors - 1st time thru StructureFactor 1142 # if 'N' in calcControls[hfx+'histType']: 1143 # dat = G2el.getBLvalues(BLtables) 1144 # refDict['FF']['FF'][iref] = dat.values() 1145 # else: #'X' 1146 # dat = G2el.getFFvalues(FFtables,SQ) 1147 # refDict['FF']['FF'][iref] = dat.values() 1148 # Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) 1149 # FF = refDict['FF']['FF'][iref][Tindx] 1150 # SSUniq = np.inner(H,SSGMT) 1151 # SSPhi = np.inner(H,SSGT) 1152 # if SGInv: #if centro - expand HKL sets 1153 # SSUniq = np.vstack((SSUniq,-SSUniq)) 1154 # SSPhi = np.hstack((SSPhi,-SSPhi)) 1155 # phase = twopi*(np.inner(SSUniq[:,:3],(dXdata.T+Xdata.T))+SSPhi[:,np.newaxis]) 1156 # sinp = np.sin(phase) 1157 # cosp = np.cos(phase) 1158 # biso = -SQfactor*Uisodata 1159 # Tiso = np.where(biso<1.,np.exp(biso),1.0) 1160 # HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in SSUniq[:,:3]]) 1161 # Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 1162 # Tcorr = Tiso*Tuij*Mdata*Fdata/len(SSUniq) 1163 # fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr]) #2 x sym x atoms 1164 # fb = np.array([FPP*cosp*Tcorr,(FF+FP-Bab)*sinp*Tcorr]) #swapped around - better? 1165 # GfpuA = G2mth.Modulation(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast) #2 x sym X atoms 1166 # GSASIIpath.IPyBreak() 1167 # fa *= GfpuA 1168 # fb *= GfpuA 1169 # fas = np.sum(np.sum(fa,axis=1),axis=1) 1170 # fbs = np.sum(np.sum(fb,axis=1),axis=1) 1171 # fasq = fas**2 1172 # fbsq = fbs**2 #imaginary 1173 # refl[9+im] = np.sum(fasq)+np.sum(fbsq) 1174 # refl[7+im] = np.sum(fasq)+np.sum(fbsq) 1175 # refl[10+im] = atan2d(fbs[0],fas[0]) 1176 # if not iref%nref: print 'ref no. %d time %.4f\r'%(iref,time.time()-time0), 1177 # print '\nref no. %d time %.4f\r'%(iref,time.time()-time0) 1178 # 1077 1179 def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict): 1078 1180 'Needs a doc string' … … 2797 2899 time0 = time.time() 2798 2900 if im: 2799 SStructureFactor (refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)2901 SStructureFactor2(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict) 2800 2902 else: 2801 2903 StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
Note: See TracChangeset
for help on using the changeset viewer.