Changeset 1992
- Timestamp:
- Oct 7, 2015 3:35:11 PM (7 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIImath.py
r1990 r1992 1062 1062 dGdBx = np.sum((Fmod*HbH)[:,:,nxs,:,nxs]*(twopi*dHdXB*np.cos(twopi*HdotXB)-np.sin(twopi*HdotXB))*glWt[nxs,nxs,nxs,:,nxs],axis=-2) 1063 1063 # GSASIIpath.IPyBreak() 1064 return np.array([cosHA,sinHA]),dGdAf,dGdBf,dGdAx,dGdBx,dGdAu,dGdBu #ops X atoms 1064 return np.array([cosHA,sinHA]),np.concatenate((dGdAf,dGdBf),-1), \ 1065 np.concatenate((dGdAx,dGdBx),-1),np.concatenate((dGdAu,dGdBu),-1) #ops X atoms 1065 1066 1066 1067 def posFourier(tau,psin,pcos,smul): -
trunk/GSASIIstrMath.py
r1990 r1992 1086 1086 print 'nRef %d time %.4f\r'%(nRef,time.time()-time0) 1087 1087 1088 #def SStructureFactor(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict):1089 # '''1090 # Compute super structure factors for all h,k,l,m for phase1091 # puts the result, F^2, in each ref[8+im] in refList1092 # input:1093 #1094 # :param dict refDict: where1095 # 'RefList' list where each ref = h,k,l,m,d,...1096 # 'FF' dict of form factors - filed in below1097 # :param np.array G: reciprocal metric tensor1098 # :param str pfx: phase id string1099 # :param dict SGData: space group info. dictionary output from SpcGroup1100 # :param dict calcControls:1101 # :param dict ParmDict:1102 #1103 # '''1104 # phfx = pfx.split(':')[0]+hfx1105 # ast = np.sqrt(np.diag(G))1106 # Mast = twopisq*np.multiply.outer(ast,ast)1107 #1108 # SGInv = SGData['SGInv']1109 # SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])1110 # SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])1111 # FFtables = calcControls['FFtables']1112 # BLtables = calcControls['BLtables']1113 # Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)1114 # waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)1115 # FF = np.zeros(len(Tdata))1116 # if 'NC' in calcControls[hfx+'histType']:1117 # FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])1118 # else:1119 # FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])1120 # FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])1121 # Uij = np.array(G2lat.U6toUij(Uijdata))1122 # bij = Mast*Uij.T1123 # if not len(refDict['FF']):1124 # if 'N' in calcControls[hfx+'histType']:1125 # dat = G2el.getBLvalues(BLtables) #will need wave here for anom. neutron b's1126 # else:1127 # dat = G2el.getFFvalues(FFtables,0.)1128 # refDict['FF']['El'] = dat.keys()1129 # refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))1130 # time0 = time.time()1131 # nref = len(refDict['RefList'])/1001132 # for iref,refl in enumerate(refDict['RefList']):1133 # if 'NT' in calcControls[hfx+'histType']:1134 # FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14+im])1135 # fbs = np.array([0,0])1136 # H = refl[:4]1137 # SQ = 1./(2.*refl[4+im])**21138 # SQfactor = 4.0*SQ*twopisq1139 # Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)1140 # if not np.any(refDict['FF']['FF'][iref]): #no form factors - 1st time thru StructureFactor1141 # if 'N' in calcControls[hfx+'histType']:1142 # dat = G2el.getBLvalues(BLtables)1143 # refDict['FF']['FF'][iref] = dat.values()1144 # else: #'X'1145 # dat = G2el.getFFvalues(FFtables,SQ)1146 # refDict['FF']['FF'][iref] = dat.values()1147 # Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])1148 # FF = refDict['FF']['FF'][iref][Tindx]1149 # SSUniq = np.inner(H,SSGMT)1150 # SSPhi = np.inner(H,SSGT)1151 # if SGInv: #if centro - expand HKL sets1152 # SSUniq = np.vstack((SSUniq,-SSUniq))1153 # SSPhi = np.hstack((SSPhi,-SSPhi))1154 # phase = twopi*(np.inner(SSUniq[:,:3],(dXdata.T+Xdata.T))+SSPhi[:,np.newaxis])1155 # sinp = np.sin(phase)1156 # cosp = np.cos(phase)1157 # biso = -SQfactor*Uisodata1158 # Tiso = np.where(biso<1.,np.exp(biso),1.0)1159 # HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in SSUniq[:,:3]])1160 # Tuij = np.where(HbH<1.,np.exp(HbH),1.0)1161 # Tcorr = Tiso*Tuij*Mdata*Fdata/len(SSUniq)1162 # fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr]) #2 x sym x atoms1163 # fb = np.array([FPP*cosp*Tcorr,(FF+FP-Bab)*sinp*Tcorr]) #swapped around - better?1164 # GfpuA = G2mth.Modulation(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast) #2 x sym X atoms1165 # GSASIIpath.IPyBreak()1166 # fa *= GfpuA1167 # fb *= GfpuA1168 # fas = np.sum(np.sum(fa,axis=1),axis=1)1169 # fbs = np.sum(np.sum(fb,axis=1),axis=1)1170 # fasq = fas**21171 # fbsq = fbs**2 #imaginary1172 # refl[9+im] = np.sum(fasq)+np.sum(fbsq)1173 # refl[7+im] = np.sum(fasq)+np.sum(fbsq)1174 # refl[10+im] = atan2d(fbs[0],fas[0])1175 # if not iref%nref: print 'ref no. %d time %.4f\r'%(iref,time.time()-time0),1176 # print '\nref no. %d time %.4f\r'%(iref,time.time()-time0)1177 #1178 1088 def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict): 1179 1089 'Needs a doc string' … … 1198 1108 nRef = len(refDict['RefList']) 1199 1109 Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict) 1200 mSize = len(Mdata) 1110 mSize = len(Mdata) #no. atoms 1201 1111 waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict) 1202 1112 FF = np.zeros(len(Tdata)) … … 1224 1134 dFdfl = np.zeros((nRef,nTwin)) 1225 1135 dFdtw = np.zeros((nRef,nTwin)) 1136 dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1],2)) 1137 dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],6)) 1138 dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],12)) 1226 1139 else: 1227 1140 dFdfr = np.zeros((nRef,mSize)) … … 1232 1145 dFdfl = np.zeros((nRef)) 1233 1146 dFdtw = np.zeros((nRef)) 1147 dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2)) 1148 dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6)) 1149 dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12)) 1234 1150 Flack = 1.0 1235 1151 if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict: … … 1267 1183 Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0]*len(TwinLaw),axis=1).T #ops x atoms 1268 1184 HbH = -np.sum(Uniq[:,nxs,:3]*np.inner(Uniq[:,:3],bij),axis=-1) #ops x atoms 1269 # GSASIIpath.IPyBreak()1270 1185 Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in Uniq]) #atoms x 3x3 1271 1186 Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6))) 1272 1187 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) #ops x atoms 1273 1188 Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0] #ops x atoms 1274 fot = (FF+FP-Bab)*occ*Tcorr #ops x atoms 1275 fotp = FPP*occ*Tcorr #ops x atoms 1276 GfpuA,dGdAf,dGdBf,dGdAx,dGdBx,dGdAu,dGdBu = \ 1277 G2mth.ModulationDerv(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast) 1189 fot = (FF+FP-Bab)*Tcorr #ops x atoms 1190 fotp = FPP*Tcorr #ops x atoms 1191 GfpuA,dGdf,dGdx,dGdu = G2mth.ModulationDerv(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast) 1278 1192 fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms) 1279 1193 fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr]) … … 1290 1204 dfadfr = np.sum(fag/occ,axis=1) #Fdata != 0 ever avoids /0. problem 1291 1205 dfbdfr = np.sum(fbg/occ,axis=1) #Fdata != 0 avoids /0. problem 1292 dfadba = np.sum(-cosp* (occ*Tcorr)[:,nxs],axis=1)1293 dfbdba = np.sum(-sinp* (occ*Tcorr)[:,nxs],axis=1)1206 dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1) 1207 dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1) 1294 1208 dfadui = np.sum(-SQfactor*fag,axis=1) 1295 1209 dfbdui = np.sum(-SQfactor*fbg,axis=1) … … 1307 1221 # array(2,nAtom,3) & array(2,nAtom,6) 1308 1222 #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for al2O3! 1223 # GSASIIpath.IPyBreak() 1309 1224 if not SGData['SGInv'] and len(TwinLaw) == 1: #Flack derivative 1310 1225 dfadfl = np.sum(-FPP*Tcorr*sinp)
Note: See TracChangeset
for help on using the changeset viewer.