Changeset 2089
- Timestamp:
- Dec 11, 2015 1:34:43 PM (7 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIImath.py
r2084 r2089 1034 1034 StauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl)) #atoms x waves x 3 x ngl 1035 1035 CtauX = np.zeros((Ax.shape[0],Ax.shape[1],3,ngl)) 1036 ZtauXt = np.zeros((Ax.shape[0],2, ngl)) #atoms x Tminmaxx ngl1036 ZtauXt = np.zeros((Ax.shape[0],2,3,ngl)) #atoms x Tminmax x 3 x ngl 1037 1037 ZtauXx = np.zeros((Ax.shape[0],3,ngl)) #atoms x XYZmax x ngl 1038 1038 for iatm in range(Ax.shape[0]): … … 1042 1042 Tmm = Ax[iatm][0][:2] 1043 1043 XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]]) 1044 for i in range(2): 1045 Tmm[i] += dT 1046 Zp = posZigZag(glTau,Tmm,XYZmax).T 1047 Tmm[i] -= 2.*dT 1048 Zm = posZigZag(glTau,Tmm,XYZmax).T 1049 Tmm[i] += dT 1050 ZtauXt[iatm][i] = (Zp-Zm)[i]/(2.*dT) 1051 # for i in range(3): 1052 # XYZmax[i] += dX 1053 # Zp = posZigZag(glTau,Tmm,XYZmax).T 1054 # XYZmax[i] -= 2.*dX 1055 # Zm = posZigZag(glTau,Tmm,XYZmax).T 1056 # XYZmax[i] += dX 1057 # ZtauXx[iatm][i] = (Zp-Zm)[i]/(2.*dX) 1058 dAdX = posZigZagDerv(glTau,Tmm,XYZmax) 1059 ZtauXx[iatm] = dAdX 1044 ZtauXt[iatm],ZtauXx[iatm] = posZigZagDerv(glTau,Tmm,XYZmax) 1060 1045 elif 'Block' in waveTypes[iatm]: 1061 1046 nx = 1 1062 1047 Tmm = Ax[iatm][0][:2] 1063 1048 XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]]) 1064 for i in range(2): 1065 Tmm[i] += dT 1066 Zp = posBlock(glTau,Tmm,XYZmax).T 1067 Tmm[i] -= 2.*dT 1068 Zm = posBlock(glTau,Tmm,XYZmax).T 1069 Tmm[i] += dT 1070 ZtauXt[iatm][i] = (Zp-Zm)[i]/(2.*dT) 1071 # for i in range(3): 1072 # XYZmax[i] += dX 1073 # Zp = posBlock(glTau,Tmm,XYZmax).T 1074 # XYZmax[i] -= 2.*dX 1075 # Zm = posBlock(glTau,Tmm,XYZmax).T 1076 # XYZmax[i] += dX 1077 ZtauXx[iatm] = posBlockDerv(glTau,Tmm,XYZmax) 1049 ZtauXt[iatm],ZtauXx[iatm] = posBlockDerv(glTau,Tmm,XYZmax) 1078 1050 tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau #Xwaves x ngl 1079 1051 if nx: … … 1145 1117 dHdXA = twopi*HP[:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[0],-1,-2)[nxs,:,:,:,:] #ops x atoms x sine waves x ngl x xyz 1146 1118 dHdXB = twopi*HP[:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[1],-1,-2)[nxs,:,:,:,:] #ditto - cos waves 1119 # ops x atoms x waves x 2xyz - real part - good 1147 1120 dGdMxCa = -np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXA*np.sin(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2) 1148 1121 dGdMxCb = -np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXB*np.sin(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2) 1149 1122 dGdMxC = np.concatenate((dGdMxCa,dGdMxCb),axis=-1) 1150 # ops x atoms x waves x 2xyz - realpart - good1123 # ops x atoms x waves x 2xyz - imag part - good 1151 1124 dGdMxSa = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXA*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2) 1152 1125 dGdMxSb = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXB*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2) 1153 1126 dGdMxS = np.concatenate((dGdMxSa,dGdMxSb),axis=-1) 1154 # ZigZag/Block waves 1155 dHdXZt = twopi*np.ones(HP.shape[0])[:,nxs,nxs,nxs]*np.swapaxes(SCtauX[2],-1,-2)[nxs,:,:,:] #??ops x atoms x ngl x 3(ZigZag/Block Tminmax)1127 # ZigZag/Block waves - problems here? 1128 dHdXZt = -twopi*HP[:,nxs,nxs,nxs,:]*np.swapaxes(SCtauX[2],-1,-2)[nxs,:,:,:,:] #ops x atoms x ngl x 2(ZigZag/Block Tminmax) 1156 1129 dHdXZx = twopi*HP[:,nxs,nxs,:]*np.swapaxes(SCtauX[3],-1,-2)[nxs,:,:,:] #ops x atoms x ngl x 3(ZigZag/Block XYZmax) 1157 dGdMzCt = -np.sum((Fmod*HbH)[:,:, :,nxs]*(dHdXZt*np.sin(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2)1130 dGdMzCt = -np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXZt*np.sin(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2) 1158 1131 dGdMzCx = -np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZx*np.sin(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2) 1159 dGdMzC = np.concatenate(( dGdMzCt,dGdMzCx),axis=-1)1160 dGdMzSt = np.sum((Fmod*HbH)[:,:, :,nxs]*(dHdXZt*np.cos(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2)1132 dGdMzC = np.concatenate((np.sum(dGdMzCt,axis=-1),dGdMzCx),axis=-1) 1133 dGdMzSt = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXZt*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2) 1161 1134 dGdMzSx = np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZx*np.cos(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2) 1162 dGdMzS = np.concatenate(( dGdMzSt,dGdMzSx),axis=-1)1135 dGdMzS = np.concatenate((np.sum(dGdMzSt,axis=-1),dGdMzSx),axis=-1) 1163 1136 # GSASIIpath.IPyBreak() 1164 # ops x atoms x waves x 2xyz - imaginary part - good1165 1137 return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS],[dGdMzC,dGdMzS] 1166 1138 … … 1170 1142 return np.sum(A,axis=0)+np.sum(B,axis=0) 1171 1143 1172 def posZigZag( tau,Tmm,XYZmax):1144 def posZigZag(T,Tmm,Xmax): 1173 1145 DT = Tmm[1]-Tmm[0] 1174 slopeUp = 2.*XYZmax/DT1175 slopeDn = 2.*XYZmax/(1.-DT)1176 A = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-X YZmax+slopeUp*((t-Tmm[0])%1.),XYZmax-slopeDn*((t-Tmm[1])%1.)) for t in tau])1146 Su = 2.*Xmax/DT 1147 Sd = 2.*Xmax/(1.-DT) 1148 A = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-Xmax+Su*((t-Tmm[0])%1.),Xmax-Sd*((t-Tmm[1])%1.)) for t in T]) 1177 1149 return A 1178 1150 1179 def posZigZagDerv( tau,Tmm,XYZmax):1151 def posZigZagDerv(T,Tmm,Xmax): 1180 1152 DT = Tmm[1]-Tmm[0] 1181 slopeUp = 2.*XYZmax/DT 1182 slopeDn = 2.*XYZmax/(1.-DT) 1183 dAdX = np.ones(3)[:,nxs]*np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-1+2*(t-Tmm[0])%1./DT,1-2*(t-Tmm[1])%1./DT) for t in tau]) 1184 return dAdX 1185 1186 1187 def posBlock(tau,Tmm,XYZmax): 1188 A = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-XYZmax,XYZmax) for t in tau]) 1153 Su = 2.*Xmax/DT 1154 Sd = 2.*Xmax/(1.-DT) 1155 dAdT = np.zeros((2,3,len(T))) 1156 dAdT[0] = np.array([np.where(Tmm[0] < t <= Tmm[1],Su*(t-Tmm[0]-1)/DT,-Sd*(t-Tmm[1])/(1.-DT)) for t in T]).T 1157 dAdT[1] = np.array([np.where(Tmm[0] < t <= Tmm[1],-Su*(t-Tmm[0])/DT,Sd*(t-Tmm[1])/(1.-DT)) for t in T]).T 1158 dAdX = np.ones(3)[:,nxs]*np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-1.+2.*(t-Tmm[0])/DT,1.-2.*(t-Tmm[1])%1./DT) for t in T]) 1159 return dAdT,dAdX 1160 1161 def posBlock(T,Tmm,Xmax): 1162 A = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-Xmax,Xmax) for t in T]) 1189 1163 return A 1190 1164 1191 def posBlockDerv(tau,Tmm,XYZmax): 1192 dAdX = np.ones(3)[:,nxs]*np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-1.,1.) for t in tau]) 1193 return dAdX 1165 def posBlockDerv(T,Tmm,Xmax): 1166 dAdT = np.zeros((2,3,len(T))) 1167 ind = np.searchsorted(T,Tmm) 1168 dAdT[0,:,ind[0]] = -Xmax/len(T) 1169 dAdT[1,:,ind[1]] = Xmax/len(T) 1170 dAdX = np.ones(3)[:,nxs]*np.array([np.where(Tmm[0] < t <= Tmm[1],-1.,1.) for t in T]) #OK 1171 return dAdT,dAdX 1194 1172 1195 1173 def fracCrenel(tau,Toff,Twid): -
trunk/GSASIIstrMath.py
r2078 r2089 888 888 dfbdfl = 0.0 889 889 #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3! 890 SA = fas[0]+fas[1] 891 SB = fbs[0]+fbs[1] 890 892 if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro 891 893 dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+ \ … … 898 900 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1]) 899 901 else: 900 SA = fas[0]+fas[1]901 SB = fbs[0]+fbs[1]902 902 if nTwin > 1: 903 903 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)] … … 1258 1258 dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1) 1259 1259 dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1) 1260 # GSASIIpath.IPyBreak() 1260 1261 if not SGData['SGInv'] and len(TwinLaw) == 1: #Flack derivative 1261 1262 dfadfl = np.sum(-FPP*Tcorr*sinp) … … 1268 1269 SB = fbs[0]+fbs[1] #float = B+B' (might be array[nTwin]) 1269 1270 if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro 1270 dFdfr[iref] = 2.*(SA*dfadfr[0]+SA*dfadfr[1]+SB*dfbdfr[0]+SB*dfbdfr[1])*Mdata/len(Uniq) #array(nRef,nAtom)1271 dFdx[iref] = 2.*(SA*dfadx[0]+SA*dfadx[1]+SB*dfbdx[0]+SB*dfbdx[1]) #array(nRef,nAtom,3)1272 dFdui[iref] = 2.*(SA*dfadui[0]+SA*dfadui[1]+SB*dfbdui[0]+SB*dfbdui[1]) #array(nRef,nAtom)1273 dFdua[iref] = 2.*(SA*dfadua[0]+SA*dfadua[1]+SB*dfbdua[0]+SB*dfbdua[1]) #array(nRef,nAtom,6)1274 1271 dFdfl[iref] = -SA*dfadfl-SB*dfbdfl #array(nRef,) 1275 1276 dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1]) #array(nRef,natom,nwave,2) 1277 dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1]) #array(nRef,natom,nwave,6) 1278 dFdGz[iref] = 2.*(SA*dfadGz[0]+SB*dfbdGz[1]) #array(nRef,natom,5) 1279 dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1]) #array(nRef,natom,nwave,12) 1280 # dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+ \ 1281 # 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq) 1282 # dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+ \ 1283 # 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1]) 1284 # dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+ \ 1285 # 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1]) 1286 # dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+ \ 1287 # 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1]) 1288 # dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+ \ 1289 # 2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1]) 1290 # dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+ \ 1291 # 2.*(fbs[0]*dfbdGx[0]+fbs[1]*dfbdGx[1]) 1292 # dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+ \ 1293 # 2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1]) 1272 dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+ \ 1273 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq) 1274 dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+ \ 1275 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1]) 1276 dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+ \ 1277 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1]) 1278 dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+ \ 1279 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1]) 1280 dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+ \ 1281 2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1]) 1282 dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+ \ 1283 2.*(fbs[0]*dfbdGx[0]-fbs[1]*dfbdGx[1]) 1284 dFdGz[iref] = 2.*(fas[0]*dfadGz[0]+fas[1]*dfadGz[1])+ \ 1285 2.*(fbs[0]*dfbdGz[0]+fbs[1]*dfbdGz[1]) 1286 dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+ \ 1287 2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1]) 1294 1288 else: 1295 1289 if nTwin > 1: … … 1315 1309 dFdGz[iref] = 2.*(SA*dfadGz[0]+SB*dfbdGz[1]) #array(nRef,natom,5) 1316 1310 dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1]) #array(nRef,natom,nwave,12) 1311 # GSASIIpath.IPyBreak() 1317 1312 dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \ 1318 1313 2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T … … 1364 1359 dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i] 1365 1360 1366 # GSASIIpath.IPyBreak()1361 # GSASIIpath.IPyBreak() 1367 1362 dFdvDict[phfx+'BabA'] = dFdbab.T[0] 1368 1363 dFdvDict[phfx+'BabU'] = dFdbab.T[1] … … 1993 1988 vec = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']]) 1994 1989 dstsq = G2lat.calc_rDsqSS(np.array([h,k,l,m]),A,vec) 1990 h,k,l = [h+m*vec[0],k+m*vec[1],l+m*vec[2]] #do proj of hklm to hkl so dPdA & dPdV come out right 1995 1991 else: 1996 1992 m = 0 … … 2023 2019 dpdDB = 1./dsp 2024 2020 dpdV = np.array([2.*h*A[0]+k*A[3]+l*A[4],2*k*A[1]+h*A[3]+l*A[5], 2025 2*l*A[2]+h*A[4]+k*A[5]])*m* *parmDict[hfx+'difC']*dsp**3/2.2021 2*l*A[2]+h*A[4]+k*A[5]])*m*parmDict[hfx+'difC']*dsp**3/2. 2026 2022 return dpdA,dpdZ,dpdDC,dpdDA,dpdDB,dpdV 2027 2023
Note: See TracChangeset
for help on using the changeset viewer.