Changeset 2111
- Timestamp:
- Jan 4, 2016 1:36:49 PM (6 years ago)
- Location:
- trunk
- Files:
-
- 4 edited
-
GSASIIindex.py (modified) (1 diff)
-
GSASIIphsGUI.py (modified) (1 diff)
-
GSASIIpwdGUI.py (modified) (2 diffs)
-
GSASIIstrMath.py (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIindex.py
r2090 r2111 313 313 for v,r in zip(ssopt['ModVec'],Vref): 314 314 if r: 315 ranges += [slice(0.0 5,1.98,.05),]315 ranges += [slice(0.02,0.98,.02),] 316 316 values += [v,] 317 317 dmin = getDmin(peaks)-0.005 -
trunk/GSASIIphsGUI.py
r2110 r2111 1691 1691 if nextName: 1692 1692 nextneigh = G2mth.FindNeighbors(data,nextName,AtNames,notName=neigh[0]) 1693 neigh[1][1].append(nextneigh[1][1][0]) 1693 if nextneigh[0]: 1694 neigh[1][1].append(nextneigh[1][1][0]) 1694 1695 neigh[2] = max(0,nH) #set expected no. H's needed 1695 1696 if len(neigh[1][0]): -
trunk/GSASIIpwdGUI.py
r2109 r2111 2470 2470 Id = Indx[ObjId] 2471 2471 try: 2472 value = min( 1.0,max(-1.0,float(Obj.GetValue())))2472 value = min(0.98,max(-0.98,float(Obj.GetValue()))) 2473 2473 except ValueError: 2474 2474 value = ssopt['ModVec'][Id] … … 2481 2481 ObjId = Obj.GetId() 2482 2482 Id,valObj = Indx[ObjId] 2483 move = Obj.GetValue()*0.0 012483 move = Obj.GetValue()*0.01 2484 2484 Obj.SetValue(0) 2485 value = min( 1.0,max(-1.0,float(valObj.GetValue())+move))2485 value = min(0.98,max(-0.98,float(valObj.GetValue())+move)) 2486 2486 valObj.SetValue('%.4f'%(value)) 2487 2487 ssopt['ModVec'][Id] = value -
trunk/GSASIIstrMath.py
r2110 r2111 771 771 # print ' %d sf time %.4f\r'%(nRef,time.time()-time0) 772 772 773 def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):774 '''Compute structure factor derivatives on single reflections - keep as it works for twins775 but is slower for powders/nontwins776 input:777 778 :param dict refDict: where779 'RefList' list where each ref = h,k,l,it,d,...780 'FF' dict of form factors - filled in below781 :param np.array G: reciprocal metric tensor782 :param str hfx: histogram id string783 :param str pfx: phase id string784 :param dict SGData: space group info. dictionary output from SpcGroup785 :param dict calcControls:786 :param dict ParmDict:787 788 :returns: dict dFdvDict: dictionary of derivatives789 '''790 phfx = pfx.split(':')[0]+hfx791 ast = np.sqrt(np.diag(G))792 Mast = twopisq*np.multiply.outer(ast,ast)793 SGMT = np.array([ops[0].T for ops in SGData['SGOps']])794 SGT = np.array([ops[1] for ops in SGData['SGOps']])795 FFtables = calcControls['FFtables']796 BLtables = calcControls['BLtables']797 TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])798 TwDict = refDict.get('TwDict',{})799 if 'S' in calcControls[hfx+'histType']:800 NTL = calcControls[phfx+'NTL']801 NM = calcControls[phfx+'TwinNMN']+1802 TwinLaw = calcControls[phfx+'TwinLaw']803 TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])804 TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))805 nTwin = len(TwinLaw)806 nRef = len(refDict['RefList'])807 Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)808 mSize = len(Mdata)809 FF = np.zeros(len(Tdata))810 if 'NC' in calcControls[hfx+'histType']:811 FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])812 elif 'X' in calcControls[hfx+'histType']:813 FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])814 FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])815 Uij = np.array(G2lat.U6toUij(Uijdata))816 bij = Mast*Uij.T817 dFdvDict = {}818 if nTwin > 1:819 dFdfr = np.zeros((nRef,nTwin,mSize))820 dFdx = np.zeros((nRef,nTwin,mSize,3))821 dFdui = np.zeros((nRef,nTwin,mSize))822 dFdua = np.zeros((nRef,nTwin,mSize,6))823 dFdbab = np.zeros((nRef,nTwin,2))824 dFdfl = np.zeros((nRef,nTwin))825 dFdtw = np.zeros((nRef,nTwin))826 else:827 dFdfr = np.zeros((nRef,mSize))828 dFdx = np.zeros((nRef,mSize,3))829 dFdui = np.zeros((nRef,mSize))830 dFdua = np.zeros((nRef,mSize,6))831 dFdbab = np.zeros((nRef,2))832 dFdfl = np.zeros((nRef))833 dFdtw = np.zeros((nRef))834 Flack = 1.0835 if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:836 Flack = 1.-2.*parmDict[phfx+'Flack']837 time0 = time.time()838 nref = len(refDict['RefList'])/100839 for iref,refl in enumerate(refDict['RefList']):840 if 'T' in calcControls[hfx+'histType']:841 FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])842 H = np.array(refl[:3])843 H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(3,nTwins) or (3)844 TwMask = np.any(H,axis=-1)845 if TwinLaw.shape[0] > 1 and TwDict:846 if iref in TwDict:847 for i in TwDict[iref]:848 for n in range(NTL):849 H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])850 TwMask = np.any(H,axis=-1)851 SQ = 1./(2.*refl[4])**2 # or (sin(theta)/lambda)**2852 SQfactor = 8.0*SQ*np.pi**2853 dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)854 Bab = parmDict[phfx+'BabA']*dBabdA855 Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])856 FF = refDict['FF']['FF'][iref].T[Tindx].T857 Uniq = np.inner(H,SGMT) # array(nSGOp,3) or (nTwin,nSGOp,3)858 Phi = np.inner(H,SGT)859 phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T860 sinp = np.sin(phase)861 cosp = np.cos(phase)862 occ = Mdata*Fdata/len(SGT)863 biso = -SQfactor*Uisodata[:,nxs]864 Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)865 HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)866 Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])867 Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))868 Tuij = np.where(HbH<1.,np.exp(HbH),1.0)869 Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ870 fot = (FF+FP-Bab)*Tcorr871 fotp = FPP*Tcorr872 fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])873 fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])874 # GSASIIpath.IPyBreak()875 fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl array(2,nTwins)876 fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl877 fax = np.array([-fot*sinp,-fotp*cosp]) #positions array(2,ntwi,nEqv,nAtoms)878 fbx = np.array([fot*cosp,-fotp*sinp])879 #sum below is over Uniq880 dfadfr = np.sum(fa/occ,axis=-2) #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem881 dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)882 dfadui = np.sum(-SQfactor*fa,axis=-2)883 if nTwin > 1:884 dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])885 dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])886 # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)887 else:888 dfadx = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2)889 dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,nxs],axis=-2)890 # array(2,nAtom,3) & array(2,nAtom,6)891 if not SGData['SGInv']:892 dfbdfr = np.sum(fb/occ,axis=-2) #Fdata != 0 avoids /0. problem893 dfadba /= 2.894 dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.895 dfbdui = np.sum(-SQfactor*fb,axis=-2)896 if len(TwinLaw) > 1:897 dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])898 dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])899 else:900 dfadfl = np.sum(-FPP*Tcorr*sinp)901 dfbdfl = np.sum(FPP*Tcorr*cosp)902 dfbdx = np.sum(twopi*Uniq*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=2)903 dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,nxs],axis=2)904 else:905 dfbdfr = np.zeros_like(dfadfr)906 dfbdx = np.zeros_like(dfadx)907 dfbdui = np.zeros_like(dfadui)908 dfbdua = np.zeros_like(dfadua)909 dfbdba = np.zeros_like(dfadba)910 dfadfl = 0.0911 dfbdfl = 0.0912 #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3!913 SA = fas[0]+fas[1]914 SB = fbs[0]+fbs[1]915 if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro916 dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+ \917 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)918 dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+ \919 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])920 dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+ \921 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])922 dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+ \923 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])924 else:925 if nTwin > 1:926 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)]927 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)]928 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)]929 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)]930 dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2931 else: #these are good for no twin single crystals932 dFdfr[iref] = (2.*SA*(dfadfr[0]+dfadfr[1])+2.*SB*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(Uniq)933 dFdx[iref] = 2.*SA*(dfadx[0]+dfadx[1])+2.*SB*(dfbdx[0]+dfbdx[1])934 dFdui[iref] = 2.*SA*(dfadui[0]+dfadui[1])+2.*SB*(dfbdui[0]+dfbdui[1])935 dFdua[iref] = 2.*SA*(dfadua[0]+dfadua[1])+2.*SB*(dfbdua[0]+dfbdua[1])936 dFdfl[iref] = -SA*dfadfl-SB*dfbdfl #array(nRef,)937 dFdbab[iref] = fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \938 fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T939 # GSASIIpath.IPyBreak()940 if not iref%100 :941 print ' %d derivative time %.4f\r'%(iref,time.time()-time0),942 print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0)943 #loop over atoms - each dict entry is list of derivatives for all the reflections944 if nTwin > 1:945 for i in range(len(Mdata)): #these all OK?946 dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)947 dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)948 dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)949 dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)950 dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)951 dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)952 dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)953 dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)954 dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)955 dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)956 dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)957 else:958 for i in range(len(Mdata)):959 dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]960 dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]961 dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]962 dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]963 dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]964 dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]965 dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]966 dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]967 dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]968 dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]969 dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]970 dFdvDict[phfx+'Flack'] = 4.*dFdfl.T971 dFdvDict[phfx+'BabA'] = dFdbab.T[0]972 dFdvDict[phfx+'BabU'] = dFdbab.T[1]973 if nTwin > 1:974 for i in range(nTwin):975 dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]976 return dFdvDict773 #def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict): 774 # '''Compute structure factor derivatives on single reflections - keep as it works for twins 775 # but is slower for powders/nontwins 776 # input: 777 # 778 # :param dict refDict: where 779 # 'RefList' list where each ref = h,k,l,it,d,... 780 # 'FF' dict of form factors - filled in below 781 # :param np.array G: reciprocal metric tensor 782 # :param str hfx: histogram id string 783 # :param str pfx: phase id string 784 # :param dict SGData: space group info. dictionary output from SpcGroup 785 # :param dict calcControls: 786 # :param dict ParmDict: 787 # 788 # :returns: dict dFdvDict: dictionary of derivatives 789 # ''' 790 # phfx = pfx.split(':')[0]+hfx 791 # ast = np.sqrt(np.diag(G)) 792 # Mast = twopisq*np.multiply.outer(ast,ast) 793 # SGMT = np.array([ops[0].T for ops in SGData['SGOps']]) 794 # SGT = np.array([ops[1] for ops in SGData['SGOps']]) 795 # FFtables = calcControls['FFtables'] 796 # BLtables = calcControls['BLtables'] 797 # TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],]) 798 # TwDict = refDict.get('TwDict',{}) 799 # if 'S' in calcControls[hfx+'histType']: 800 # NTL = calcControls[phfx+'NTL'] 801 # NM = calcControls[phfx+'TwinNMN']+1 802 # TwinLaw = calcControls[phfx+'TwinLaw'] 803 # TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))]) 804 # TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1)) 805 # nTwin = len(TwinLaw) 806 # nRef = len(refDict['RefList']) 807 # Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict) 808 # mSize = len(Mdata) 809 # FF = np.zeros(len(Tdata)) 810 # if 'NC' in calcControls[hfx+'histType']: 811 # FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam']) 812 # elif 'X' in calcControls[hfx+'histType']: 813 # FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata]) 814 # FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata]) 815 # Uij = np.array(G2lat.U6toUij(Uijdata)) 816 # bij = Mast*Uij.T 817 # dFdvDict = {} 818 # if nTwin > 1: 819 # dFdfr = np.zeros((nRef,nTwin,mSize)) 820 # dFdx = np.zeros((nRef,nTwin,mSize,3)) 821 # dFdui = np.zeros((nRef,nTwin,mSize)) 822 # dFdua = np.zeros((nRef,nTwin,mSize,6)) 823 # dFdbab = np.zeros((nRef,nTwin,2)) 824 # dFdfl = np.zeros((nRef,nTwin)) 825 # dFdtw = np.zeros((nRef,nTwin)) 826 # else: 827 # dFdfr = np.zeros((nRef,mSize)) 828 # dFdx = np.zeros((nRef,mSize,3)) 829 # dFdui = np.zeros((nRef,mSize)) 830 # dFdua = np.zeros((nRef,mSize,6)) 831 # dFdbab = np.zeros((nRef,2)) 832 # dFdfl = np.zeros((nRef)) 833 # dFdtw = np.zeros((nRef)) 834 # Flack = 1.0 835 # if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict: 836 # Flack = 1.-2.*parmDict[phfx+'Flack'] 837 # time0 = time.time() 838 # nref = len(refDict['RefList'])/100 839 # for iref,refl in enumerate(refDict['RefList']): 840 # if 'T' in calcControls[hfx+'histType']: 841 # FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12]) 842 # H = np.array(refl[:3]) 843 # H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(3,nTwins) or (3) 844 # TwMask = np.any(H,axis=-1) 845 # if TwinLaw.shape[0] > 1 and TwDict: 846 # if iref in TwDict: 847 # for i in TwDict[iref]: 848 # for n in range(NTL): 849 # H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM]) 850 # TwMask = np.any(H,axis=-1) 851 # SQ = 1./(2.*refl[4])**2 # or (sin(theta)/lambda)**2 852 # SQfactor = 8.0*SQ*np.pi**2 853 # dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor) 854 # Bab = parmDict[phfx+'BabA']*dBabdA 855 # Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) 856 # FF = refDict['FF']['FF'][iref].T[Tindx].T 857 # Uniq = np.inner(H,SGMT) # array(nSGOp,3) or (nTwin,nSGOp,3) 858 # Phi = np.inner(H,SGT) 859 # phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T 860 # sinp = np.sin(phase) 861 # cosp = np.cos(phase) 862 # occ = Mdata*Fdata/len(SGT) 863 # biso = -SQfactor*Uisodata[:,nxs] 864 # Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1) 865 # HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1) 866 # Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))]) 867 # Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6))) 868 # Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 869 # Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ 870 # fot = (FF+FP-Bab)*Tcorr 871 # fotp = FPP*Tcorr 872 # fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) 873 # fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr]) 874 ## GSASIIpath.IPyBreak() 875 # fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl array(2,nTwins) 876 # fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl 877 # fax = np.array([-fot*sinp,-fotp*cosp]) #positions array(2,ntwi,nEqv,nAtoms) 878 # fbx = np.array([fot*cosp,-fotp*sinp]) 879 # #sum below is over Uniq 880 # dfadfr = np.sum(fa/occ,axis=-2) #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem 881 # dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1) 882 # dfadui = np.sum(-SQfactor*fa,axis=-2) 883 # if nTwin > 1: 884 # dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 885 # dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 886 # # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6) 887 # else: 888 # dfadx = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2) 889 # dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,nxs],axis=-2) 890 # # array(2,nAtom,3) & array(2,nAtom,6) 891 # if not SGData['SGInv']: 892 # dfbdfr = np.sum(fb/occ,axis=-2) #Fdata != 0 avoids /0. problem 893 # dfadba /= 2. 894 # dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2. 895 # dfbdui = np.sum(-SQfactor*fb,axis=-2) 896 # if len(TwinLaw) > 1: 897 # dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)]) 898 # dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)]) 899 # else: 900 # dfadfl = np.sum(-FPP*Tcorr*sinp) 901 # dfbdfl = np.sum(FPP*Tcorr*cosp) 902 # dfbdx = np.sum(twopi*Uniq*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=2) 903 # dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,nxs],axis=2) 904 # else: 905 # dfbdfr = np.zeros_like(dfadfr) 906 # dfbdx = np.zeros_like(dfadx) 907 # dfbdui = np.zeros_like(dfadui) 908 # dfbdua = np.zeros_like(dfadua) 909 # dfbdba = np.zeros_like(dfadba) 910 # dfadfl = 0.0 911 # dfbdfl = 0.0 912 # #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3! 913 # SA = fas[0]+fas[1] 914 # SB = fbs[0]+fbs[1] 915 # if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro 916 # dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+ \ 917 # 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq) 918 # dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+ \ 919 # 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1]) 920 # dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+ \ 921 # 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1]) 922 # dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+ \ 923 # 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1]) 924 # else: 925 # if nTwin > 1: 926 # 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)] 927 # 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)] 928 # 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)] 929 # 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)] 930 # dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2 931 # else: #these are good for no twin single crystals 932 # dFdfr[iref] = (2.*SA*(dfadfr[0]+dfadfr[1])+2.*SB*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(Uniq) 933 # dFdx[iref] = 2.*SA*(dfadx[0]+dfadx[1])+2.*SB*(dfbdx[0]+dfbdx[1]) 934 # dFdui[iref] = 2.*SA*(dfadui[0]+dfadui[1])+2.*SB*(dfbdui[0]+dfbdui[1]) 935 # dFdua[iref] = 2.*SA*(dfadua[0]+dfadua[1])+2.*SB*(dfbdua[0]+dfbdua[1]) 936 # dFdfl[iref] = -SA*dfadfl-SB*dfbdfl #array(nRef,) 937 # dFdbab[iref] = fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \ 938 # fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T 939 ## GSASIIpath.IPyBreak() 940 # if not iref%100 : 941 # print ' %d derivative time %.4f\r'%(iref,time.time()-time0), 942 # print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0) 943 # #loop over atoms - each dict entry is list of derivatives for all the reflections 944 # if nTwin > 1: 945 # for i in range(len(Mdata)): #these all OK? 946 # dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0) 947 # dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0) 948 # dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0) 949 # dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0) 950 # dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0) 951 # dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0) 952 # dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0) 953 # dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0) 954 # dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0) 955 # dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0) 956 # dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0) 957 # else: 958 # for i in range(len(Mdata)): 959 # dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i] 960 # dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i] 961 # dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i] 962 # dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i] 963 # dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i] 964 # dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i] 965 # dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i] 966 # dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i] 967 # dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i] 968 # dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i] 969 # dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i] 970 # dFdvDict[phfx+'Flack'] = 4.*dFdfl.T 971 # dFdvDict[phfx+'BabA'] = dFdbab.T[0] 972 # dFdvDict[phfx+'BabU'] = dFdbab.T[1] 973 # if nTwin > 1: 974 # for i in range(nTwin): 975 # dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i] 976 # return dFdvDict 977 977 978 978 def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict): … … 1134 1134 return dFdvDict 1135 1135 1136 def StructureFactorDervTw(refDict,G,hfx,pfx,SGData,calcControls,parmDict):1137 '''Compute structure factor derivatives on single reflections - for twins only1138 input:1139 1140 :param dict refDict: where1141 'RefList' list where each ref = h,k,l,it,d,...1142 'FF' dict of form factors - filled in below1143 :param np.array G: reciprocal metric tensor1144 :param str hfx: histogram id string1145 :param str pfx: phase id string1146 :param dict SGData: space group info. dictionary output from SpcGroup1147 :param dict calcControls:1148 :param dict parmDict:1149 1150 :returns: dict dFdvDict: dictionary of derivatives1151 '''1152 phfx = pfx.split(':')[0]+hfx1153 ast = np.sqrt(np.diag(G))1154 Mast = twopisq*np.multiply.outer(ast,ast)1155 SGMT = np.array([ops[0].T for ops in SGData['SGOps']])1156 SGT = np.array([ops[1] for ops in SGData['SGOps']])1157 FFtables = calcControls['FFtables']1158 BLtables = calcControls['BLtables']1159 TwDict = refDict.get('TwDict',{})1160 NTL = calcControls[phfx+'NTL']1161 NM = calcControls[phfx+'TwinNMN']+11162 TwinLaw = calcControls[phfx+'TwinLaw']1163 TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])1164 TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))1165 nTwin = len(TwinLaw)1166 nRef = len(refDict['RefList'])1167 Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)1168 mSize = len(Mdata)1169 FF = np.zeros(len(Tdata))1170 if 'NC' in calcControls[hfx+'histType']:1171 FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])1172 elif 'X' in calcControls[hfx+'histType']:1173 FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])1174 FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])1175 Uij = np.array(G2lat.U6toUij(Uijdata))1176 bij = Mast*Uij.T1177 dFdvDict = {}1178 dFdfr = np.zeros((nRef,nTwin,mSize))1179 dFdx = np.zeros((nRef,nTwin,mSize,3))1180 dFdui = np.zeros((nRef,nTwin,mSize))1181 dFdua = np.zeros((nRef,nTwin,mSize,6))1182 dFdbab = np.zeros((nRef,nTwin,2))1183 dFdtw = np.zeros((nRef,nTwin))1184 time0 = time.time()1185 nref = len(refDict['RefList'])/1001186 for iref,refl in enumerate(refDict['RefList']):1187 if 'T' in calcControls[hfx+'histType']:1188 FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])1189 H = np.array(refl[:3])1190 H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(3,nTwins) or (3)1191 TwMask = np.any(H,axis=-1)1192 if iref in TwDict:1193 for i in TwDict[iref]:1194 for n in range(NTL):1195 H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])1196 TwMask = np.any(H,axis=-1)1197 SQ = 1./(2.*refl[4])**2 # or (sin(theta)/lambda)**21198 SQfactor = 8.0*SQ*np.pi**21199 dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)1200 Bab = parmDict[phfx+'BabA']*dBabdA1201 Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])1202 FF = refDict['FF']['FF'][iref].T[Tindx].T1203 Uniq = np.inner(H,SGMT) # array(nSGOp,3) or (nTwin,nSGOp,3)1204 Phi = np.inner(H,SGT)1205 phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T1206 sinp = np.sin(phase)1207 cosp = np.cos(phase)1208 occ = Mdata*Fdata/len(SGT)1209 biso = -SQfactor*Uisodata[:,nxs]1210 Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1)1211 HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)1212 Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])1213 Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))1214 Tuij = np.where(HbH<1.,np.exp(HbH),1.0)1215 Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ1216 fot = (FF+FP-Bab)*Tcorr1217 fotp = FPP*Tcorr1218 fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-FPP*sinp*Tcorr])1219 fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,FPP*cosp*Tcorr])1220 # GSASIIpath.IPyBreak()1221 fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl array(2,nTwins)1222 fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl1223 if SGData['SGInv']: #centrosymmetric; B=01224 fbs[0] *= 0.1225 fas[1] *= 0.1226 fax = np.array([-fot*sinp,-fotp*cosp]) #positions array(2,ntwi,nEqv,nAtoms)1227 fbx = np.array([fot*cosp,-fotp*sinp])1228 #sum below is over Uniq1229 dfadfr = np.sum(fa/occ,axis=-2) #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem1230 dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)1231 dfadui = np.sum(-SQfactor*fa,axis=-2)1232 dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])1233 dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])1234 # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)1235 if not SGData['SGInv']:1236 dfbdfr = np.sum(fb/occ,axis=-2) #Fdata != 0 avoids /0. problem1237 dfadba /= 2.1238 dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.1239 dfbdui = np.sum(-SQfactor*fb,axis=-2)1240 dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])1241 dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])1242 else:1243 dfbdfr = np.zeros_like(dfadfr)1244 dfbdx = np.zeros_like(dfadx)1245 dfbdui = np.zeros_like(dfadui)1246 dfbdua = np.zeros_like(dfadua)1247 dfbdba = np.zeros_like(dfadba)1248 SA = fas[0]+fas[1]1249 SB = fbs[0]+fbs[1]1250 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)]1251 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)]1252 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)]1253 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)]1254 if SGData['SGInv']: #centrosymmetric; B=01255 dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**21256 else:1257 dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**21258 dFdbab[iref] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \1259 fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T1260 # GSASIIpath.IPyBreak()1261 print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0)1262 #loop over atoms - each dict entry is list of derivatives for all the reflections1263 for i in range(len(Mdata)): #these all OK?1264 dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)1265 dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)1266 dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0)1267 dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0)1268 dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0)1269 dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0)1270 dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0)1271 dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0)1272 dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0)1273 dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0)1274 dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0)1275 dFdvDict[phfx+'BabA'] = dFdbab.T[0]1276 dFdvDict[phfx+'BabU'] = dFdbab.T[1]1277 for i in range(nTwin):1278 dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i]1279 return dFdvDict1136 #def StructureFactorDervTw(refDict,G,hfx,pfx,SGData,calcControls,parmDict): 1137 # '''Compute structure factor derivatives on single reflections - for twins only 1138 # input: 1139 # 1140 # :param dict refDict: where 1141 # 'RefList' list where each ref = h,k,l,it,d,... 1142 # 'FF' dict of form factors - filled in below 1143 # :param np.array G: reciprocal metric tensor 1144 # :param str hfx: histogram id string 1145 # :param str pfx: phase id string 1146 # :param dict SGData: space group info. dictionary output from SpcGroup 1147 # :param dict calcControls: 1148 # :param dict parmDict: 1149 # 1150 # :returns: dict dFdvDict: dictionary of derivatives 1151 # ''' 1152 # phfx = pfx.split(':')[0]+hfx 1153 # ast = np.sqrt(np.diag(G)) 1154 # Mast = twopisq*np.multiply.outer(ast,ast) 1155 # SGMT = np.array([ops[0].T for ops in SGData['SGOps']]) 1156 # SGT = np.array([ops[1] for ops in SGData['SGOps']]) 1157 # FFtables = calcControls['FFtables'] 1158 # BLtables = calcControls['BLtables'] 1159 # TwDict = refDict.get('TwDict',{}) 1160 # NTL = calcControls[phfx+'NTL'] 1161 # NM = calcControls[phfx+'TwinNMN']+1 1162 # TwinLaw = calcControls[phfx+'TwinLaw'] 1163 # TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))]) 1164 # TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1)) 1165 # nTwin = len(TwinLaw) 1166 # nRef = len(refDict['RefList']) 1167 # Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict) 1168 # mSize = len(Mdata) 1169 # FF = np.zeros(len(Tdata)) 1170 # if 'NC' in calcControls[hfx+'histType']: 1171 # FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam']) 1172 # elif 'X' in calcControls[hfx+'histType']: 1173 # FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata]) 1174 # FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata]) 1175 # Uij = np.array(G2lat.U6toUij(Uijdata)) 1176 # bij = Mast*Uij.T 1177 # dFdvDict = {} 1178 # dFdfr = np.zeros((nRef,nTwin,mSize)) 1179 # dFdx = np.zeros((nRef,nTwin,mSize,3)) 1180 # dFdui = np.zeros((nRef,nTwin,mSize)) 1181 # dFdua = np.zeros((nRef,nTwin,mSize,6)) 1182 # dFdbab = np.zeros((nRef,nTwin,2)) 1183 # dFdtw = np.zeros((nRef,nTwin)) 1184 # time0 = time.time() 1185 # nref = len(refDict['RefList'])/100 1186 # for iref,refl in enumerate(refDict['RefList']): 1187 # if 'T' in calcControls[hfx+'histType']: 1188 # FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12]) 1189 # H = np.array(refl[:3]) 1190 # H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(3,nTwins) or (3) 1191 # TwMask = np.any(H,axis=-1) 1192 # if iref in TwDict: 1193 # for i in TwDict[iref]: 1194 # for n in range(NTL): 1195 # H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM]) 1196 # TwMask = np.any(H,axis=-1) 1197 # SQ = 1./(2.*refl[4])**2 # or (sin(theta)/lambda)**2 1198 # SQfactor = 8.0*SQ*np.pi**2 1199 # dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor) 1200 # Bab = parmDict[phfx+'BabA']*dBabdA 1201 # Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) 1202 # FF = refDict['FF']['FF'][iref].T[Tindx].T 1203 # Uniq = np.inner(H,SGMT) # array(nSGOp,3) or (nTwin,nSGOp,3) 1204 # Phi = np.inner(H,SGT) 1205 # phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T 1206 # sinp = np.sin(phase) 1207 # cosp = np.cos(phase) 1208 # occ = Mdata*Fdata/len(SGT) 1209 # biso = -SQfactor*Uisodata[:,nxs] 1210 # Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1) 1211 # HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1) 1212 # Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))]) 1213 # Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6))) 1214 # Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 1215 # Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ 1216 # fot = (FF+FP-Bab)*Tcorr 1217 # fotp = FPP*Tcorr 1218 # fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-FPP*sinp*Tcorr]) 1219 # fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,FPP*cosp*Tcorr]) 1220 ## GSASIIpath.IPyBreak() 1221 # fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl array(2,nTwins) 1222 # fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl 1223 # if SGData['SGInv']: #centrosymmetric; B=0 1224 # fbs[0] *= 0. 1225 # fas[1] *= 0. 1226 # fax = np.array([-fot*sinp,-fotp*cosp]) #positions array(2,ntwi,nEqv,nAtoms) 1227 # fbx = np.array([fot*cosp,-fotp*sinp]) 1228 # #sum below is over Uniq 1229 # dfadfr = np.sum(fa/occ,axis=-2) #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem 1230 # dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1) 1231 # dfadui = np.sum(-SQfactor*fa,axis=-2) 1232 # dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1233 # dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1234 # # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6) 1235 # if not SGData['SGInv']: 1236 # dfbdfr = np.sum(fb/occ,axis=-2) #Fdata != 0 avoids /0. problem 1237 # dfadba /= 2. 1238 # dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2. 1239 # dfbdui = np.sum(-SQfactor*fb,axis=-2) 1240 # dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)]) 1241 # dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)]) 1242 # else: 1243 # dfbdfr = np.zeros_like(dfadfr) 1244 # dfbdx = np.zeros_like(dfadx) 1245 # dfbdui = np.zeros_like(dfadui) 1246 # dfbdua = np.zeros_like(dfadua) 1247 # dfbdba = np.zeros_like(dfadba) 1248 # SA = fas[0]+fas[1] 1249 # SB = fbs[0]+fbs[1] 1250 # 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)] 1251 # 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)] 1252 # 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)] 1253 # 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)] 1254 # if SGData['SGInv']: #centrosymmetric; B=0 1255 # dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2 1256 # else: 1257 # dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2 1258 # dFdbab[iref] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \ 1259 # fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T 1260 ## GSASIIpath.IPyBreak() 1261 # print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0) 1262 # #loop over atoms - each dict entry is list of derivatives for all the reflections 1263 # for i in range(len(Mdata)): #these all OK? 1264 # dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0) 1265 # dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0) 1266 # dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0) 1267 # dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0) 1268 # dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0) 1269 # dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0) 1270 # dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0) 1271 # dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0) 1272 # dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0) 1273 # dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0) 1274 # dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0) 1275 # dFdvDict[phfx+'BabA'] = dFdbab.T[0] 1276 # dFdvDict[phfx+'BabU'] = dFdbab.T[1] 1277 # for i in range(nTwin): 1278 # dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i] 1279 # return dFdvDict 1280 1280 1281 1281 def StructureFactorDervTw2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
Note: See TracChangeset
for help on using the changeset viewer.
