Changeset 2070
- Timestamp:
- Nov 25, 2015 2:19:24 PM (7 years ago)
- Location:
- trunk
- Files:
-
- 5 edited
-
GSASIImath.py (modified) (9 diffs)
-
GSASIIphsGUI.py (modified) (1 diff)
-
GSASIIplot.py (modified) (2 diffs)
-
GSASIIstrIO.py (modified) (1 diff)
-
GSASIIstrMath.py (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIImath.py
r2067 r2070 957 957 nx = 0 958 958 if 'ZigZag' in waveTypes[iatm]: 959 n t= 1959 nx = 1 960 960 Tmm = Ax[iatm][0][:2] 961 961 XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]]) … … 963 963 slopeUp = 2.*XYZmax/DT 964 964 slopeDn = 2.*XYZmax/(1.-DT) 965 XmodZ[iatm][0] = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-XYZmax+slopeUp*((t-Tmm[0])%1.),XYZmax-slopeDn*((t-Tmm[1])%1.)) for t in glTau]).T965 XmodZ[iatm][0] += np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-XYZmax+slopeUp*((t-Tmm[0])%1.),XYZmax-slopeDn*((t-Tmm[1])%1.)) for t in glTau]).T 966 966 elif 'Block' in waveTypes[iatm]: 967 n t= 1967 nx = 1 968 968 Tmm = Ax[iatm][0][:2] 969 969 XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]]) 970 XmodZ[iatm][0] = np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-XYZmax,XYZmax) for t in glTau]).T970 XmodZ[iatm][0] += np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-XYZmax,XYZmax) for t in glTau]).T 971 971 tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau #Xwaves x 32 972 XmodA[iatm] = Ax[iatm,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X 3 X 32973 XmodB[iatm] = Bx[iatm,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto972 XmodA[iatm][nx:] = Ax[iatm,nx:,:,nxs]*np.sin(twopi*tauX[nxs,:,nxs,:]) #atoms X waves X 3 X 32 973 XmodB[iatm][nx:] = Bx[iatm,nx:,:,nxs]*np.cos(twopi*tauX[nxs,:,nxs,:]) #ditto 974 974 Xmod = np.sum(XmodA+XmodB+XmodZ,axis=1) #atoms X 3 X 32; sum waves 975 975 Xmod = np.swapaxes(Xmod,1,2) … … 981 981 else: 982 982 Umod = 1.0 983 # GSASIIpath.IPyBreak() 983 984 return nWaves,Fmod,Xmod,Umod,glTau,glWt 984 985 … … 1023 1024 Af = np.array(FSSdata[0]).T #sin frac mods x waves x atoms 1024 1025 Bf = np.array(FSSdata[1]).T #cos frac mods... 1025 Ax = np.array(XSSdata[:3]).T #...cos pos mods 1026 Ax = np.array(XSSdata[:3]).T #...cos pos mods x waves x atoms 1026 1027 Bx = np.array(XSSdata[3:]).T #...cos pos mods 1027 1028 Au = Mast*np.array(G2lat.U6toUij(USSdata[:6])).T #atoms x waves x sin Uij mods 1028 1029 Bu = Mast*np.array(G2lat.U6toUij(USSdata[6:])).T #...cos Uij mods 1029 1030 nWaves = [Af.shape[1],Ax.shape[1],Au.shape[1]] 1030 nf = 01031 1031 nx = 0 1032 # if 'Fourier' in waveTypes: 1033 # nf = 0 1034 # nx = 0 1035 # XmodZ = np.zeros((Ax.shape[0],Ax.shape[1],3,32)) 1036 # FmodZ = np.zeros((Af.shape[0],Af.shape[1],32)) 1037 # if 'Crenel' in waveTypes: 1038 # nC = np.where('Crenel' in waveTypes) 1039 # nf = 1 1040 # #FmodZ = 0 replace 1041 # else: 1042 # nx = 1 1043 # if 'Sawtooth' in waveTypes: 1044 # nS = np.where('Sawtooth' in waveTypes) 1045 # #XmodZ = 0 replace 1046 tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau #Xwaves x 32 1047 StauX = np.ones_like(Ax)[:,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:] #atoms X waves X 3(xyz) X 32 1048 CtauX = np.ones_like(Bx)[:,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:] #ditto 1032 StauX = np.zeros((Ax.shape[0],Ax.shape[1],3,32)) #atoms x waves x 3 x 32 1033 CtauX = np.zeros((Ax.shape[0],Ax.shape[1],3,32)) 1034 ZtauXx = np.zeros((Ax.shape[0],3,32)) #atoms x XYZmax x 32 1035 ZtauXt = np.zeros((Ax.shape[0],2,32)) #atoms x Tminmax x 32 1036 for iatm in range(Ax.shape[0]): 1037 nx = 0 1038 if 'ZigZag' in waveTypes[iatm]: 1039 nx = 1 1040 # Tmm = Ax[iatm][0][:2] 1041 # XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]]) 1042 # DT = Tmm[1]-Tmm[0] 1043 # slopeUp = 2.*XYZmax/DT 1044 # slopeDn = 2.*XYZmax/(1.-DT) 1045 # XmodZ[iatm][0] += np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-XYZmax+slopeUp*((t-Tmm[0])%1.),XYZmax-slopeDn*((t-Tmm[1])%1.)) for t in glTau]).T 1046 elif 'Block' in waveTypes[iatm]: 1047 nx = 1 1048 Tmm = Ax[iatm][0][:2] 1049 XYZmax = np.array([Ax[iatm][0][2],Bx[iatm][0][0],Bx[iatm][0][1]]) 1050 ZtauXt[iatm] = np.ones(2)[:,nxs] 1051 ZtauXx[iatm] += np.array([np.where(Tmm[0] < t%1. <= Tmm[1],-np.ones(3),np.ones(3)) for t in glTau]).T 1052 tauX = np.arange(1.,nWaves[1]+1-nx)[:,nxs]*glTau #Xwaves x 32 1053 StauX[iatm][nx:] = np.ones_like(Ax)[iatm,nx:,:,nxs]*np.sin(twopi*tauX)[nxs,:,nxs,:] #atoms X waves X 3(xyz) X 32 1054 CtauX[iatm][nx:] = np.ones_like(Bx)[iatm,nx:,:,nxs]*np.cos(twopi*tauX)[nxs,:,nxs,:] #ditto 1055 # GSASIIpath.IPyBreak() 1049 1056 if nWaves[0]: 1050 1057 tauF = np.arange(1.,nWaves[0]+1-nf)[:,nxs]*glTau #Fwaves x 32 1051 StauF = np.ones_like(Af)[:, nf:,nxs]*np.sin(twopi*tauF)[nxs,:,:] #also dFmod/dAf1052 CtauF = np.ones_like(Bf)[:, nf:,nxs]*np.cos(twopi*tauF)[nxs,:,:] #also dFmod/dBf1058 StauF = np.ones_like(Af)[:,:,nxs]*np.sin(twopi*tauF)[nxs,:,:] #also dFmod/dAf 1059 CtauF = np.ones_like(Bf)[:,:,nxs]*np.cos(twopi*tauF)[nxs,:,:] #also dFmod/dBf 1053 1060 else: 1054 1061 StauF = 1.0 … … 1068 1075 UmodA = 0. 1069 1076 UmodB = 0. 1070 return waveShapes,[StauF,CtauF],[StauX,CtauX ],[StauU,CtauU],UmodA+UmodB1077 return waveShapes,[StauF,CtauF],[StauX,CtauX,ZtauXt,ZtauXx],[StauU,CtauU],UmodA+UmodB 1071 1078 1072 1079 def ModulationDerv(H,HP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt): … … 1116 1123 dGdMxSb = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(dHdXB*np.cos(HdotXD)[:,:,nxs,:,nxs])*glWt[nxs,nxs,nxs,:,nxs],axis=-2) 1117 1124 dGdMxS = np.concatenate((dGdMxSa,dGdMxSb),axis=-1) 1125 # ZigZag/Block waves 1126 dHdXZt = twopi*np.ones(HP.shape[0])[:,nxs,nxs,nxs]*np.swapaxes(SCtauX[2],-1,-2)[nxs,:,:,:] #??ops x atoms x 32 x 3(ZigZag/Block Tminmax) 1127 dHdXZx = twopi*HP[:,nxs,nxs,:]*np.swapaxes(SCtauX[3],-1,-2)[nxs,:,:,:] #ops x atoms x 32 x 3(ZigZag/Block XYZmax) 1128 dGdMzCt = -np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZt*np.sin(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2) 1129 dGdMzCx = -np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZx*np.sin(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2) 1130 dGdMzC = np.concatenate((dGdMzCt,dGdMzCx),axis=-1) 1131 dGdMzSt = np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZt*np.cos(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2) 1132 dGdMzSx = np.sum((Fmod*HbH)[:,:,:,nxs]*(dHdXZx*np.cos(HdotXD)[:,:,:,nxs])*glWt[nxs,nxs,:,nxs],axis=-2) 1133 dGdMzS = np.concatenate((dGdMzSt,dGdMzSx),axis=-1) 1118 1134 # ops x atoms x waves x 2xyz - imaginary part - good 1119 return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS] 1120 1121 def posFourier(tau,psin,pcos,smul): 1122 A = np.array([ps[:,np.newaxis]*np.sin(2*np.pi*(i+1)*tau) for i,ps in enumerate(psin)]) #*smul 1135 # GSASIIpath.IPyBreak() 1136 return [dGdMfC,dGdMfS],[dGdMxC,dGdMxS],[dGdMuC,dGdMuS],[dGdMzC,dGdMzS] 1137 1138 def posFourier(tau,psin,pcos): 1139 A = np.array([ps[:,np.newaxis]*np.sin(2*np.pi*(i+1)*tau) for i,ps in enumerate(psin)]) 1123 1140 B = np.array([pc[:,np.newaxis]*np.cos(2*np.pi*(i+1)*tau) for i,pc in enumerate(pcos)]) 1124 1141 return np.sum(A,axis=0)+np.sum(B,axis=0) … … 1175 1192 sop,ssop,icent = G2spc.OpsfromStringOps(opr,SGData,SSGData) 1176 1193 sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(tau,sop,ssop,atxyz) 1177 smul = sdet # n-fold vs m operator on wave1178 1194 tauT *= icent #invert wave on -1 1179 1195 wave = np.zeros(3) … … 1206 1222 ccof.append(spos[0][3:]) 1207 1223 if len(scof): 1208 wave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof) ,smul),axis=1)1224 wave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof)),axis=1) 1209 1225 if len(Sadp): 1210 1226 scof = [] … … 1213 1229 scof.append(sadp[0][:6]) 1214 1230 ccof.append(sadp[0][6:]) 1215 uwave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof) ,smul),axis=1)1231 uwave += np.sum(posFourier(tauT,np.array(scof),np.array(ccof)),axis=1) 1216 1232 if atom[cia] == 'A': 1217 1233 X,U = G2spc.ApplyStringOps(opr,SGData,atxyz+wave,atuij+uwave) -
trunk/GSASIIphsGUI.py
r2064 r2070 2417 2417 waveSizer.Add(waveHead) 2418 2418 if len(waveBlk): 2419 n Four= 02419 nx = 0 2420 2420 for iwave,wave in enumerate(waveBlk): 2421 if waveType == 'Fourier':2422 nFour += 12423 2421 if not iwave: 2424 CSI = G2spc.GetSSfxuinel(waveType,nFour,xyz,SGData,SSGData) 2422 if waveType in ['ZigZag','Block']: 2423 nx = 1 2424 CSI = G2spc.GetSSfxuinel(waveType,1,xyz,SGData,SSGData) 2425 2425 else: 2426 CSI = G2spc.GetSSfxuinel('Fourier', nFour,xyz,SGData,SSGData)2426 CSI = G2spc.GetSSfxuinel('Fourier',iwave+1-nx,xyz,SGData,SSGData) 2427 2427 waveName = 'Fourier' 2428 2428 if Stype == 'Sfrac': -
trunk/GSASIIplot.py
r2065 r2070 3121 3121 scof.append(spos[0][:3]) 3122 3122 ccof.append(spos[0][3:]) 3123 wave += G2mth.posFourier(tau,np.array(scof),np.array(ccof) ,1)3123 wave += G2mth.posFourier(tau,np.array(scof),np.array(ccof)) 3124 3124 if mapData['Flip']: 3125 3125 Title = 'Charge flip' … … 3150 3150 Plot.set_xlabel('t') 3151 3151 Plot.set_ylabel(r'$\mathsf{\Delta}$%s'%(Ax)) 3152 Slab = np.hstack((slab,slab,slab)) 3152 Slab = np.hstack((slab,slab,slab)) 3153 3153 acolor = mpl.cm.get_cmap('RdYlGn') 3154 3154 if 'delt' in MapType: -
trunk/GSASIIstrIO.py
r2064 r2070 1168 1168 for Stype in ['Sfrac','Spos','Sadp','Smag']: 1169 1169 Waves = AtomSS[Stype] 1170 nx = 0 1170 1171 for iw,wave in enumerate(Waves): 1171 1172 if not iw: 1172 CSI = G2spc.GetSSfxuinel(waveType,iw+1,at[cx:cx+3],SGData,SSGData) 1173 if waveType in ['ZigZag','Block']: 1174 nx = 1 1175 CSI = G2spc.GetSSfxuinel(waveType,1,at[cx:cx+3],SGData,SSGData) 1173 1176 else: 1174 CSI = G2spc.GetSSfxuinel('Fourier',iw+1 ,at[cx:cx+3],SGData,SSGData)1177 CSI = G2spc.GetSSfxuinel('Fourier',iw+1-nx,at[cx:cx+3],SGData,SSGData) 1175 1178 uId,uCoef = CSI[Stype] 1176 1179 stiw = str(i)+':'+str(iw) -
trunk/GSASIIstrMath.py
r2067 r2070 1164 1164 dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2)) 1165 1165 dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6)) 1166 dFdGz = np.zeros((nRef,mSize,5)) 1166 1167 dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12)) 1167 1168 Flack = 1.0 … … 1210 1211 fotp = FPP*Tcorr #ops x atoms 1211 1212 GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms 1212 dGdf,dGdx,dGdu = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt)1213 dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt) 1213 1214 # GfpuA is 2 x ops x atoms 1214 # derivs are: ops x atoms x waves x 2,6, or 12parms as [real,imag] parts1215 # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts 1215 1216 fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms) 1216 1217 fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr]) #or array(2,nEqv,nAtoms) … … 1241 1242 dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1) 1242 1243 dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1) 1244 dfadGz = np.sum(fa[:,it,:,0,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:]-fb[:,it,:,0,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:],axis=1) 1245 dfbdGz = np.sum(fb[:,it,:,0,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:]+fa[:,it,:,0,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:],axis=1) 1243 1246 dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1) 1244 1247 dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1) … … 1253 1256 dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1) 1254 1257 dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1) 1258 dfadGz = np.sum(fa[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]-fb[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1) 1259 dfbdGz = np.sum(fb[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]+fa[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1) 1255 1260 dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1) 1256 1261 dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1) … … 1298 1303 dFdGf[iref] = [2.*TwMask[it]*(SA[it]*dfadGf[1]+SB[it]*dfbdGf[1]) for it in range(nTwin)] 1299 1304 dFdGx[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)] 1305 dFdGz[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)] 1300 1306 dFdGu[iref] = [2.*TwMask[it]*(SA[it]*dfadGu[1]+SB[it]*dfbdGu[1]) for it in range(nTwin)] 1301 1307 else: #these are good for no twin single crystals … … 1308 1314 dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1]) #array(nRef,natom,nwave,2) 1309 1315 dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1]) #array(nRef,natom,nwave,6) 1316 dFdGz[iref] = 2.*(SA*dfadGz[0]+SB*dfbdGz[1]) #array(nRef,natom,5) 1310 1317 dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1]) #array(nRef,natom,nwave,12) 1311 1318 dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \ … … 1330 1337 dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i] 1331 1338 dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i] 1332 for j in range(XSSdata.shape[1]): #loop over waves 'Tzero','Xslope','Yslope','Zslope'?1339 for j in range(XSSdata.shape[1]): #loop over waves 1333 1340 dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j)] = dFdGx.T[0][j][i] 1334 1341 dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j)] = dFdGx.T[1][j][i] … … 1337 1344 dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j)] = dFdGx.T[4][j][i] 1338 1345 dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j)] = dFdGx.T[5][j][i] 1346 dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i] 1347 dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[0][i] 1348 dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[0][i] 1349 dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[0][i] 1350 dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[0][i] 1339 1351 for j in range(USSdata.shape[1]): #loop over waves 1340 1352 dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i]
Note: See TracChangeset
for help on using the changeset viewer.
