Changeset 2520
- Timestamp:
- Nov 12, 2016 11:17:34 AM (7 years ago)
- Location:
- trunk
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIconstrGUI.py
r2512 r2520 1184 1184 NB: A = [G11,G22,G33,2*G12,2*G13,2*G23] 1185 1185 ''' 1186 1187 def SetUniqAj(pId,iA,Aname,SGLaue): 1188 if SGLaue in ['4/m','4/mmm'] and iA in [0,1]: 1189 parm = '%d::%s'%(pId,'A0') 1190 elif SGLaue in ['m3','m3m'] and iA in [0,1,2]: 1191 parm = '%d::%s'%(pId,'A0') 1192 elif SGLaue in ['6/m','6/mmm','3m1', '31m', '3'] and iA in [0,1,3]: 1193 parm = '%d::%s'%(pId,'A0') 1194 elif SGLaue in ['3R', '3mR']: 1195 if ia in [0,1,2]: 1196 parm = '%d::%s'%(pId,'A0') 1197 else: 1198 parm = '%d::%s'%(pId,'A3') 1199 else: 1200 parm = '%d::%s'%(pId,Aname) 1201 return parm 1202 1186 1203 Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree() 1187 1204 UseList = newPhase['Histograms'] 1188 1205 detTrans = np.abs(nl.det(Trans)) 1189 1206 invTrans = nl.inv(Trans) 1190 oCell = oldPhase['General']['Cell'][1:7] 1191 nCell = newPhase['General']['Cell'][1:7] 1192 oGmat = G2lat.cell2Gmat(oldPhase['General']['Cell'][1:7])[0] 1193 nGmat = G2lat.cell2Gmat(newPhase['General']['Cell'][1:7])[0] 1194 Gtrans = np.inner(nGmat,nl.inv(oGmat)) 1195 print 'oA',G2lat.cell2A(oldPhase['General']['Cell'][1:7]) 1196 print 'nA',G2lat.cell2A(newPhase['General']['Cell'][1:7]) 1197 print 'oGmat',oGmat 1198 print 'nGmat',nGmat 1199 print 'Gtrans',Gtrans 1207 # print 'invTrans',invTrans 1200 1208 nAcof = G2lat.cell2A(newPhase['General']['Cell'][1:7]) 1201 1209 … … 1206 1214 cx,ct,cs,cia = newPhase['General']['AtomPtrs'] 1207 1215 nAtoms = newPhase['Atoms'] 1216 oAtoms = oldPhase['Atoms'] 1208 1217 oSGData = oldPhase['General']['SGData'] 1209 1218 nSGData = newPhase['General']['SGData'] 1219 oAcof = G2lat.cell2A(oldPhase['General']['Cell'][1:7]) 1220 nAcof = G2lat.cell2A(newPhase['General']['Cell'][1:7]) 1221 oGmat = G2lat.cell2Gmat(oldPhase['General']['Cell'][1:7])[1] 1222 nGmat = G2lat.cell2Gmat(newPhase['General']['Cell'][1:7])[1] 1210 1223 item = G2gd.GetPatternTreeItemId(G2frame,G2frame.root,'Constraints') 1211 1224 constraints = G2frame.PatternTree.GetItemPyData(item) 1212 1225 # GSASIIpath.IPyBreak() 1226 parmDict = {} 1227 varyList = [] 1213 1228 for ia,code in enumerate(atCodes): 1214 1229 atom = nAtoms[ia] … … 1243 1258 #how do I do Uij's for most Trans? 1244 1259 Anames = [['A0','A3','A4'],['A3','A1','A5'],['A4','A5','A2']] 1245 Aids = [[0,0,'A0'],[1,1,'A1'],[2,2,'A2'],[0,1,'A3'],[0,2,'A4'],[1,2,'A5']] 1260 As = ['A0','A1','A2','A3','A4','A5'] 1261 Aids = [[0,0,'A0',-1],[1,1,'A1',-1],[2,2,'A2',-1],[0,1,'A3',2],[0,2,'A4',1],[1,2,'A5',0]] 1246 1262 Axes = ['a','b','c'] 1247 1263 Holds = [] 1248 #how do I invoke Laue symmetry matches for Laue symm change? 1249 for iA,Aid in enumerate(Aids): 1250 IndpCon = [1.0,G2obj.G2VarObj('%d::%s'%(npId,Aid[2]))] 1264 for iA,Aid in enumerate(Aids): 1265 parm = SetUniqAj(opId,iA,Aid[2],oSGData['SGLaue']) 1266 parmDict[parm] = oAcof[iA] 1267 varyList.append(parm) 1268 IndpCon = [1.0,G2obj.G2VarObj(parm)] 1251 1269 DepCons = [] 1252 1270 for iat in range(3): 1253 if nSGData['SGLaue'] in ['-1','2/m']: 1254 if (abs(nAcof[iA]) < 1.e-8) and (abs( Gtrans[iat][Aid[0]]) < 1.e-8):1271 if nSGData['SGLaue'] in ['-1','2/m']: #set holds 1272 if (abs(nAcof[iA]) < 1.e-8) and (abs(Trans[Aid[0],Aid[1]]) < 1.e-8): 1255 1273 if Axes[iat] != oSGData['SGUniq'] and oSGData['SGLaue'] != nSGData['SGLaue']: 1256 1274 HoldObj = G2obj.G2VarObj('%d::%s'%(npId,Aid[2])) … … 1258 1276 constraints['Phase'].append([[0.0,HoldObj],None,None,'h']) 1259 1277 Holds.append(HoldObj) 1260 print constraints['Phase'][-1] 1261 if abs(Gtrans[iat][Aid[0]]) > 1.e-8 and abs(nAcof[iA]) > 1.e-8: 1262 print iat,Aid,Gtrans[iat][Aid[1]],'%d::%s'%(opId,Anames[iat][Aid[1]]) 1263 DepCons.append([Gtrans[iat][Aid[1]],G2obj.G2VarObj('%d::%s'%(opId,Anames[iat][Aid[1]]))]) 1278 continue 1279 # print iA,Aid,iat,invTrans[iat][Aid[0]],invTrans[Aid[1]][iat],Anames[Aid[0]][Aid[1]],parm 1280 if abs(invTrans[iat,Aid[1]]) > 1.e-8 and abs(nAcof[iA]) > 1.e-8: 1281 parm = SetUniqAj(npId,iA,Anames[Aid[0]][Aid[1]],nSGData['SGLaue']) 1282 parmDict[parm] = nAcof[As.index(Aid[2])] 1283 if not parm in varyList: 1284 varyList.append(parm) 1285 DepCons.append([invTrans[Aid[0],Aid[0]]*invTrans[Aid[1],Aid[1]],G2obj.G2VarObj(parm)]) 1264 1286 if len(DepCons) == 1: 1265 1287 constraints['Phase'].append([IndpCon,DepCons[0],None,None,'e']) … … 1268 1290 Dep[0] *= -1 1269 1291 constraints['Phase'].append([IndpCon]+DepCons+[0.0,None,'c']) 1292 # constDict,fixedList,ignored = G2stIO.ProcessConstraints(constraints['Phase']) 1293 # groups,parmlist = G2mv.GroupConstraints(constDict) 1294 # G2mv.GenerateConstraints(groups,parmlist,varyList,constDict,fixedList,parmDict) 1295 # print 'old',parmDict 1296 # G2mv.Dict2Map(parmDict,varyList) 1297 # print 'new',parmDict 1270 1298 for hId,hist in enumerate(UseList): #HAP - seems OK 1271 1299 ohapkey = '%d:%d:'%(opId,hId) -
trunk/GSASIIgrid.py
r2516 r2520 700 700 label=' NB: Nonmagnetic atoms will be deleted from new phase'),0,WACV) 701 701 constr = wx.CheckBox(self.panel,label=' Make constraints between phases?') 702 mainSizer.Add(wx.StaticText(self.panel, \ 703 label=' Constraints not correct for non-diagonal transforms'),0,WACV) 702 704 constr.SetValue(self.ifConstr) 703 705 constr.Bind(wx.EVT_CHECKBOX,OnConstr) -
trunk/GSASIIstrIO.py
r2495 r2520 1556 1556 Vol = 1/np.sqrt(rVsq) 1557 1557 sigVol = Vol**3*np.sqrt(srcvlsq)/2. #ok - checks with GSAS 1558 1558 1559 R123 = Ax[0]*Ax[1]*Ax[2] 1559 1560 dsasdg = np.zeros((3,6)) … … 1585 1586 dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0] 1586 1587 dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr 1588 # GSASIIpath.IPyBreak() 1587 1589 sigMat = np.inner(dadg,np.inner(dadg,vcov)) 1588 1590 var = np.diag(sigMat) -
trunk/GSASIIstrMath.py
r2512 r2520 636 636 return np.array(waveTypes),FSSdata,XSSdata,USSdata,MSSdata 637 637 638 #def GetSSTauM(SGOps,SSOps,pfx,calcControls,XData):639 #640 # Natoms = calcControls['Natoms'][pfx]641 # maxSSwave = calcControls['maxSSwave'][pfx]642 # Smult = np.zeros((Natoms,len(SGOps)))643 # TauT = np.zeros((Natoms,len(SGOps)))644 # for ix,xyz in enumerate(XData.T):645 # for isym,(sop,ssop) in enumerate(zip(SGOps,SSOps)):646 # sdet,ssdet,dtau,dT,tauT = G2spc.getTauT(0,sop,ssop,xyz)647 # Smult[ix][isym] = sdet648 # TauT[ix][isym] = tauT649 # return Smult,TauT650 #651 638 def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict): 652 639 ''' Compute structure factors for all h,k,l for phase … … 824 811 # print ' %d sf time %.4f\r'%(nRef,time.time()-time0) 825 812 826 #def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):827 # '''Compute structure factor derivatives on single reflections - keep as it works for twins828 # but is slower for powders/nontwins829 # input:830 #831 # :param dict refDict: where832 # 'RefList' list where each ref = h,k,l,it,d,...833 # 'FF' dict of form factors - filled in below834 # :param np.array G: reciprocal metric tensor835 # :param str hfx: histogram id string836 # :param str pfx: phase id string837 # :param dict SGData: space group info. dictionary output from SpcGroup838 # :param dict calcControls:839 # :param dict ParmDict:840 #841 # :returns: dict dFdvDict: dictionary of derivatives842 # '''843 # phfx = pfx.split(':')[0]+hfx844 # ast = np.sqrt(np.diag(G))845 # Mast = twopisq*np.multiply.outer(ast,ast)846 # SGMT = np.array([ops[0].T for ops in SGData['SGOps']])847 # SGT = np.array([ops[1] for ops in SGData['SGOps']])848 # FFtables = calcControls['FFtables']849 # BLtables = calcControls['BLtables']850 # TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])851 # TwDict = refDict.get('TwDict',{})852 # if 'S' in calcControls[hfx+'histType']:853 # NTL = calcControls[phfx+'NTL']854 # NM = calcControls[phfx+'TwinNMN']+1855 # TwinLaw = calcControls[phfx+'TwinLaw']856 # TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])857 # TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))858 # nTwin = len(TwinLaw)859 # nRef = len(refDict['RefList'])860 # Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \861 # GetAtomFXU(pfx,calcControls,parmDict)862 # mSize = len(Mdata)863 # FF = np.zeros(len(Tdata))864 # if 'NC' in calcControls[hfx+'histType']:865 # FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])866 # elif 'X' in calcControls[hfx+'histType']:867 # FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])868 # FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])869 # Uij = np.array(G2lat.U6toUij(Uijdata))870 # bij = Mast*Uij.T871 # dFdvDict = {}872 # if nTwin > 1:873 # dFdfr = np.zeros((nRef,nTwin,mSize))874 # dFdx = np.zeros((nRef,nTwin,mSize,3))875 # dFdui = np.zeros((nRef,nTwin,mSize))876 # dFdua = np.zeros((nRef,nTwin,mSize,6))877 # dFdbab = np.zeros((nRef,nTwin,2))878 # dFdfl = np.zeros((nRef,nTwin))879 # dFdtw = np.zeros((nRef,nTwin))880 # else:881 # dFdfr = np.zeros((nRef,mSize))882 # dFdx = np.zeros((nRef,mSize,3))883 # dFdui = np.zeros((nRef,mSize))884 # dFdua = np.zeros((nRef,mSize,6))885 # dFdbab = np.zeros((nRef,2))886 # dFdfl = np.zeros((nRef))887 # dFdtw = np.zeros((nRef))888 # Flack = 1.0889 # if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:890 # Flack = 1.-2.*parmDict[phfx+'Flack']891 # time0 = time.time()892 # nref = len(refDict['RefList'])/100893 # for iref,refl in enumerate(refDict['RefList']):894 # if 'T' in calcControls[hfx+'histType']:895 # FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])896 # H = np.array(refl[:3])897 # H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(3,nTwins) or (3)898 # TwMask = np.any(H,axis=-1)899 # if TwinLaw.shape[0] > 1 and TwDict:900 # if iref in TwDict:901 # for i in TwDict[iref]:902 # for n in range(NTL):903 # H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])904 # TwMask = np.any(H,axis=-1)905 # SQ = 1./(2.*refl[4])**2 # or (sin(theta)/lambda)**2906 # SQfactor = 8.0*SQ*np.pi**2907 # dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)908 # Bab = parmDict[phfx+'BabA']*dBabdA909 # Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])910 # FF = refDict['FF']['FF'][iref].T[Tindx].T911 # Uniq = np.inner(H,SGMT) # array(nSGOp,3) or (nTwin,nSGOp,3)912 # Phi = np.inner(H,SGT)913 # phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T914 # sinp = np.sin(phase)915 # cosp = np.cos(phase)916 # occ = Mdata*Fdata/len(SGT)917 # biso = -SQfactor*Uisodata[:,nxs]918 # Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)919 # HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)920 # Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])921 # Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))922 # Tuij = np.where(HbH<1.,np.exp(HbH),1.0)923 # Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ924 # fot = (FF+FP-Bab)*Tcorr925 # fotp = FPP*Tcorr926 # fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])927 # fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])928 ## GSASIIpath.IPyBreak()929 # fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl array(2,nTwins)930 # fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl931 # fax = np.array([-fot*sinp,-fotp*cosp]) #positions array(2,ntwi,nEqv,nAtoms)932 # fbx = np.array([fot*cosp,-fotp*sinp])933 # #sum below is over Uniq934 # dfadfr = np.sum(fa/occ,axis=-2) #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem935 # dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)936 # dfadui = np.sum(-SQfactor*fa,axis=-2)937 # if nTwin > 1:938 # dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])939 # dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])940 # # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)941 # else:942 # dfadx = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)943 # dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,nxs],axis=-2)944 # # array(2,nAtom,3) & array(2,nAtom,6)945 # if not SGData['SGInv']:946 # dfbdfr = np.sum(fb/occ,axis=-2) #Fdata != 0 avoids /0. problem947 # dfadba /= 2.948 # dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.949 # dfbdui = np.sum(-SQfactor*fb,axis=-2)950 # if len(TwinLaw) > 1:951 # dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])952 # dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])953 # else:954 # dfadfl = np.sum(-FPP*Tcorr*sinp)955 # dfbdfl = np.sum(FPP*Tcorr*cosp)956 # dfbdx = np.sum(twopi*Uniq*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=2)957 # dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,nxs],axis=2)958 # else:959 # dfbdfr = np.zeros_like(dfadfr)960 # dfbdx = np.zeros_like(dfadx)961 # dfbdui = np.zeros_like(dfadui)962 # dfbdua = np.zeros_like(dfadua)963 # dfbdba = np.zeros_like(dfadba)964 # dfadfl = 0.0965 # dfbdfl = 0.0966 # #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!967 # SA = fas[0]+fas[1]968 # SB = fbs[0]+fbs[1]969 # if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro970 # dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+ \971 # 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)972 # dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+ \973 # 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])974 # dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+ \975 # 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])976 # dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+ \977 # 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])978 # else:979 # if nTwin > 1:980 # 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)]981 # 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)]982 # dFdui[iref] = [2.*TwMask[it]*(SA[it]*dfadui[it][0]+SA[it]*dfadui[it][1]+SB[it]*dfbdui[it][0]+SB[it]*dfbdui[it][1]) for it in range(nTwin)]983 # 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)]984 # dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2985 # else: #these are good for no twin single crystals986 # dFdfr[iref] = (2.*SA*(dfadfr[0]+dfadfr[1])+2.*SB*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(Uniq)987 # dFdx[iref] = 2.*SA*(dfadx[0]+dfadx[1])+2.*SB*(dfbdx[0]+dfbdx[1])988 # dFdui[iref] = 2.*SA*(dfadui[0]+dfadui[1])+2.*SB*(dfbdui[0]+dfbdui[1])989 # dFdua[iref] = 2.*SA*(dfadua[0]+dfadua[1])+2.*SB*(dfbdua[0]+dfbdua[1])990 # dFdfl[iref] = -SA*dfadfl-SB*dfbdfl #array(nRef,)991 # dFdbab[iref] = fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \992 # fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T993 ## GSASIIpath.IPyBreak()994 # if not iref%100 :995 # print ' %d derivative time %.4f\r'%(iref,time.time()-time0),996 # print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0)997 # #loop over atoms - each dict entry is list of derivatives for all the reflections998 # if nTwin > 1:999 # for i in range(len(Mdata)): #these all OK?1000 # dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)1001 # dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)1002 # dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)1003 # dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)1004 # dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)1005 # dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)1006 # dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)1007 # dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)1008 # dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)1009 # dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)1010 # dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)1011 # else:1012 # for i in range(len(Mdata)):1013 # dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]1014 # dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]1015 # dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]1016 # dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]1017 # dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]1018 # dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]1019 # dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]1020 # dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]1021 # dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]1022 # dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]1023 # dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]1024 # dFdvDict[phfx+'Flack'] = 4.*dFdfl.T1025 # dFdvDict[phfx+'BabA'] = dFdbab.T[0]1026 # dFdvDict[phfx+'BabU'] = dFdbab.T[1]1027 # if nTwin > 1:1028 # for i in range(nTwin):1029 # dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]1030 # return dFdvDict1031 1032 813 def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict): 1033 814 '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only … … 1194 975 def StructureFactorDervMag(refDict,G,hfx,pfx,SGData,calcControls,parmDict): 1195 976 '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only 1196 faster than StructureFactorDerv - correct for powders/nontwins!!1197 977 input: 1198 978 … … 1232 1012 Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata #flip vectors according to spin flip 1233 1013 Gdata = np.inner(Amat,Gdata.T) #convert back to cart. space MXYZ, Natoms, NOps*Inv*Ncen 1014 dGdM = SGData['MagMom'][nxs,:,nxs]*dGdM 1234 1015 Gdata = np.swapaxes(Gdata,1,2) # put Natoms last - Mxyz,Nops,Natms 1235 1016 # GSASIIpath.IPyBreak() … … 1244 1025 dFdfr = np.zeros((nRef,mSize)) 1245 1026 dFdx = np.zeros((nRef,mSize,3)) 1246 dFdMx = np.zeros(( nRef,mSize,3))1027 dFdMx = np.zeros((3,nRef,mSize)) 1247 1028 dFdui = np.zeros((nRef,mSize)) 1248 1029 dFdua = np.zeros((nRef,mSize,6)) … … 1266 1047 HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1) 1267 1048 Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T 1268 # GSASIIpath.IPyBreak()1269 1049 Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))]) 1270 1050 Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,len(SGT),6)) … … 1289 1069 eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0) 1290 1070 Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Gdata[:,nxs,:,:] #Mxyz,Nref,Nop,Natm = BPM in magstrfc.for OK 1291 #there is still something wrong with the next three lines = dF/dmag not corect yet1292 1071 dqdm = np.array([np.outer(hm,hm)-np.eye(3) for hm in HM.T]).T #Mxyz,Mxyz,Nref (3x3 matrix) 1293 dqmx = np.sum(dqdm[:,:,:,nxs,nxs]*dGdm[:,nxs,nxs,:,:],axis=0) #matrix * vector = vector 1294 dmx = Q*dGdM[:,nxs,:,:]+dqmx #*Mag canceled out of dqmx term 1072 dqmx = np.sum(dqdm[:,:,:,nxs,nxs]*dGdm[:,nxs,nxs,:1,:],axis=0) #matrix * vector = vector 1073 dmx = Q*dGdM[:,nxs,:1,:]+dqmx #*Mag canceled out of dqmx term 1074 # GSASIIpath.IPyBreak() 1295 1075 # 1296 1076 fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:] #Mxyz,Nref,Nop,Natm … … 1303 1083 dfadfr = np.sum(fam/occ,axis=2) #array(Mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem deriv OK 1304 1084 dfadx = np.sum(-twopi*Uniq[nxs,:,:,nxs,:]*famx[:,:,:,:,nxs],axis=2) #deriv OK 1305 dfadmx = np.sum( TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*dmx,axis=2)1085 dfadmx = np.sum(dmx*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:],axis=2) 1306 1086 dfadui = np.sum(-SQfactor[:,nxs,nxs]*fam,axis=2) #array(Ops,refBlk,nAtoms) deriv OK 1307 1087 dfadua = np.sum(-Hij[nxs,:,:,nxs,:]*fam[:,:,:,:,nxs],axis=2) #deriv OK? not U12 & U23 in sarc … … 1309 1089 dfbdfr = np.sum(fbm/occ,axis=2) #array(mxyz,refBlk,nAtom) Fdata != 0 avoids /0. problem 1310 1090 dfbdx = np.sum(-twopi*Uniq[nxs,:,:,nxs,:]*fbmx[:,:,:,:,nxs],axis=2) 1311 dfbdmx = np.sum( TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*dmx,axis=2)1091 dfbdmx = np.sum(dmx*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:],axis=2) 1312 1092 dfbdui = np.sum(-SQfactor[:,nxs,nxs]*fbm,axis=2) #array(Ops,refBlk,nAtoms) 1313 1093 dfbdua = np.sum(-Hij[nxs,:,:,nxs,:]*fbm[:,:,:,:,nxs],axis=2) … … 1315 1095 dFdfr[iBeg:iFin] = 2.*np.sum((fams[:,:,nxs]*dfadfr+fbms[:,:,nxs]*dfbdfr)*Mdata/(2*Nops*Ncen),axis=0) 1316 1096 dFdx[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadx+fbms[:,:,nxs,nxs]*dfbdx,axis=0) 1317 # GSASIIpath.IPyBreak() 1318 dFdMx[iBeg:iFin] = np.reshape(2.*fams[:,:,nxs]*dfadmx+fbms[:,:,nxs]*dfbdmx,(iFin-iBeg,-1,3)) 1097 dFdMx[:,iBeg:iFin,:] = 2.*(fams[:,:,nxs]*dfadmx+fbms[:,:,nxs]*dfbdmx) 1319 1098 dFdui[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs]*dfadui+fbms[:,:,nxs]*dfbdui,axis=0) 1320 1099 dFdua[iBeg:iFin] = 2.*np.sum(fams[:,:,nxs,nxs]*dfadua+fbms[:,:,nxs,nxs]*dfbdua,axis=0) 1100 # GSASIIpath.IPyBreak() 1321 1101 iBeg += blkSize 1322 1102 print ' %d derivative time %.4f\r'%(nRef,time.time()-time0) … … 1327 1107 dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i] 1328 1108 dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i] 1329 dFdvDict[pfx+'AMx:'+str(i)] = dFdMx .T[0][i]1330 dFdvDict[pfx+'AMy:'+str(i)] = dFdMx .T[1][i]1331 dFdvDict[pfx+'AMz:'+str(i)] = dFdMx .T[2][i]1109 dFdvDict[pfx+'AMx:'+str(i)] = dFdMx[0,:,i] 1110 dFdvDict[pfx+'AMy:'+str(i)] = dFdMx[1,:,i] 1111 dFdvDict[pfx+'AMz:'+str(i)] = dFdMx[2,:,i] 1332 1112 dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i] 1333 1113 dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i] … … 1339 1119 # GSASIIpath.IPyBreak() 1340 1120 return dFdvDict 1341 1342 #def StructureFactorDervTw(refDict,G,hfx,pfx,SGData,calcControls,parmDict): 1343 # '''Compute structure factor derivatives on single reflections - for twins only 1344 # input: 1345 # 1346 # :param dict refDict: where 1347 # 'RefList' list where each ref = h,k,l,it,d,... 1348 # 'FF' dict of form factors - filled in below 1349 # :param np.array G: reciprocal metric tensor 1350 # :param str hfx: histogram id string 1351 # :param str pfx: phase id string 1352 # :param dict SGData: space group info. dictionary output from SpcGroup 1353 # :param dict calcControls: 1354 # :param dict parmDict: 1355 # 1356 # :returns: dict dFdvDict: dictionary of derivatives 1357 # ''' 1358 # phfx = pfx.split(':')[0]+hfx 1359 # ast = np.sqrt(np.diag(G)) 1360 # Mast = twopisq*np.multiply.outer(ast,ast) 1361 # SGMT = np.array([ops[0].T for ops in SGData['SGOps']]) 1362 # SGT = np.array([ops[1] for ops in SGData['SGOps']]) 1363 # FFtables = calcControls['FFtables'] 1364 # BLtables = calcControls['BLtables'] 1365 # TwDict = refDict.get('TwDict',{}) 1366 # NTL = calcControls[phfx+'NTL'] 1367 # NM = calcControls[phfx+'TwinNMN']+1 1368 # TwinLaw = calcControls[phfx+'TwinLaw'] 1369 # TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))]) 1370 # TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1)) 1371 # nTwin = len(TwinLaw) 1372 # nRef = len(refDict['RefList']) 1373 # Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \ 1374 # GetAtomFXU(pfx,calcControls,parmDict) 1375 # mSize = len(Mdata) 1376 # FF = np.zeros(len(Tdata)) 1377 # if 'NC' in calcControls[hfx+'histType']: 1378 # FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam']) 1379 # elif 'X' in calcControls[hfx+'histType']: 1380 # FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata]) 1381 # FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata]) 1382 # Uij = np.array(G2lat.U6toUij(Uijdata)) 1383 # bij = Mast*Uij.T 1384 # dFdvDict = {} 1385 # dFdfr = np.zeros((nRef,nTwin,mSize)) 1386 # dFdx = np.zeros((nRef,nTwin,mSize,3)) 1387 # dFdui = np.zeros((nRef,nTwin,mSize)) 1388 # dFdua = np.zeros((nRef,nTwin,mSize,6)) 1389 # dFdbab = np.zeros((nRef,nTwin,2)) 1390 # dFdtw = np.zeros((nRef,nTwin)) 1391 # time0 = time.time() 1392 # nref = len(refDict['RefList'])/100 1393 # for iref,refl in enumerate(refDict['RefList']): 1394 # if 'T' in calcControls[hfx+'histType']: 1395 # FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12]) 1396 # H = np.array(refl[:3]) 1397 # H = np.inner(H.T,TwinLaw) #maybe array(3,nTwins) or (3) 1398 # TwMask = np.any(H,axis=-1) 1399 # if iref in TwDict: 1400 # for i in TwDict[iref]: 1401 # for n in range(NTL): 1402 # H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM]) 1403 # TwMask = np.any(H,axis=-1) 1404 # SQ = 1./(2.*refl[4])**2 # or (sin(theta)/lambda)**2 1405 # SQfactor = 8.0*SQ*np.pi**2 1406 # dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor) 1407 # Bab = parmDict[phfx+'BabA']*dBabdA 1408 # Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) 1409 # FF = refDict['FF']['FF'][iref].T[Tindx].T 1410 # Uniq = np.inner(H,SGMT) # array(nSGOp,3) or (nTwin,nSGOp,3) 1411 # Phi = np.inner(H,SGT) 1412 # phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T 1413 # sinp = np.sin(phase) 1414 # cosp = np.cos(phase) 1415 # occ = Mdata*Fdata/len(SGT) 1416 # biso = -SQfactor*Uisodata[:,nxs] 1417 # Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1) 1418 # HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1) 1419 # Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))]) 1420 # Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)) 1421 # Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 1422 # Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ 1423 # fot = (FF+FP-Bab)*Tcorr 1424 # fotp = FPP*Tcorr 1425 # fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-FPP*sinp*Tcorr]) 1426 # fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,FPP*cosp*Tcorr]) 1427 ## GSASIIpath.IPyBreak() 1428 # fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl array(2,nTwins) 1429 # fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl 1430 # if SGData['SGInv']: #centrosymmetric; B=0 1431 # fbs[0] *= 0. 1432 # fas[1] *= 0. 1433 # fax = np.array([-fot*sinp,-fotp*cosp]) #positions array(2,ntwi,nEqv,nAtoms) 1434 # fbx = np.array([fot*cosp,-fotp*sinp]) 1435 # #sum below is over Uniq 1436 # dfadfr = np.sum(fa/occ,axis=-2) #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem 1437 # dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1) 1438 # dfadui = np.sum(-SQfactor*fa,axis=-2) 1439 # dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1440 # dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1441 # # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6) 1442 # if not SGData['SGInv']: 1443 # dfbdfr = np.sum(fb/occ,axis=-2) #Fdata != 0 avoids /0. problem 1444 # dfadba /= 2. 1445 # dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2. 1446 # dfbdui = np.sum(-SQfactor*fb,axis=-2) 1447 # dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)]) 1448 # dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)]) 1449 # else: 1450 # dfbdfr = np.zeros_like(dfadfr) 1451 # dfbdx = np.zeros_like(dfadx) 1452 # dfbdui = np.zeros_like(dfadui) 1453 # dfbdua = np.zeros_like(dfadua) 1454 # dfbdba = np.zeros_like(dfadba) 1455 # SA = fas[0]+fas[1] 1456 # SB = fbs[0]+fbs[1] 1457 # 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)] 1458 # 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)] 1459 # 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)] 1460 # 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)] 1461 # if SGData['SGInv']: #centrosymmetric; B=0 1462 # dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2 1463 # else: 1464 # dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2 1465 # dFdbab[iref] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \ 1466 # fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T 1467 ## GSASIIpath.IPyBreak() 1468 # print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0) 1469 # #loop over atoms - each dict entry is list of derivatives for all the reflections 1470 # for i in range(len(Mdata)): #these all OK? 1471 # dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0) 1472 # dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0) 1473 # dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0) 1474 # dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0) 1475 # dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0) 1476 # dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0) 1477 # dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0) 1478 # dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0) 1479 # dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0) 1480 # dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0) 1481 # dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0) 1482 # dFdvDict[phfx+'BabA'] = dFdbab.T[0] 1483 # dFdvDict[phfx+'BabU'] = dFdbab.T[1] 1484 # for i in range(nTwin): 1485 # dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i] 1486 # return dFdvDict 1487 1121 1488 1122 def StructureFactorDervTw2(refDict,G,hfx,pfx,SGData,calcControls,parmDict): 1489 1123 '''Compute structure factor derivatives on blocks of reflections - for twins only
Note: See TracChangeset
for help on using the changeset viewer.