Changeset 1887
- Timestamp:
- Jun 12, 2015 3:13:05 PM (8 years ago)
- Location:
- trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIddataGUI.py
r1885 r1887 780 780 781 781 def OnAddTwin(event): 782 twinMat = np.array([[ 1,0,0],[0,1,0],[0,0,1]])783 twinVal = [ 1.0,False]782 twinMat = np.array([[-1,0,0],[0,-1,0],[0,0,-1]]) #inversion by default 783 twinVal = [0.0,False] 784 784 UseList[G2frame.hist]['Twins'].append([twinMat,twinVal]) 785 785 addtwin.SetValue(False) … … 802 802 try: 803 803 val = float(Obj.GetValue()) 804 if 0. > val > 1.: \804 if 0. > val > 1.: 805 805 raise ValueError 806 806 except ValueError: 807 807 val = UseList[G2frame.hist]['Twins'][it][1][0] 808 808 UseList[G2frame.hist]['Twins'][it][1][0] = val 809 Obj.SetValue('%.3f'%(val)) 809 sumTw = 0. 810 for it,twin in enumerate(UseList[G2frame.hist]['Twins']): 811 if it: 812 sumTw += twin[1][0] 813 UseList[G2frame.hist]['Twins'][0][1][0] = 1.-sumTw 814 wx.CallLater(100,RepaintHistogramInfo) 810 815 811 816 def OnTwinRef(event): … … 818 823 it = Indx[Obj.GetId()] 819 824 del UseList[G2frame.hist]['Twins'][it] 825 sumTw = 0. 826 for it,twin in enumerate(UseList[G2frame.hist]['Twins']): 827 if it: 828 sumTw += twin[1][0] 829 UseList[G2frame.hist]['Twins'][0][1][0] = 1.-sumTw 820 830 wx.CallLater(100,RepaintHistogramInfo) 821 831 -
trunk/GSASIIstrIO.py
r1886 r1887 2426 2426 print >>pFile,sigstr 2427 2427 2428 def PrintTwinsAndSig(pfx, hapData,TwinSig):2428 def PrintTwinsAndSig(pfx,twinData,TwinSig): 2429 2429 print >>pFile,'\n Twin Law fractions : ' 2430 2430 ptlbls = ' names :' 2431 2431 ptstr = ' values:' 2432 2432 sigstr = ' sig :' 2433 for it em in hapData:2434 ptlbls += ' %12s'%(item)2435 ptstr += '%12.3f'%( hapData[item][0])2436 if pfx+ itemin TwinSig:2437 sigstr += '%12.3f'%(TwinSig[pfx+ item])2433 for it,item in enumerate(twinData): 2434 ptlbls += ' twin #%d'%(it) 2435 ptstr += '%12.3f'%(item[1][0]) 2436 if pfx+'TwinFr;'+str(it) in TwinSig: 2437 sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr;'+str(it)]) 2438 2438 else: 2439 2439 sigstr += 12*' ' … … 2534 2534 if pfx+item in sigDict: 2535 2535 BabSig[pfx+item] = sigDict[pfx+item] 2536 i tem = 'TwinFr'2537 it = 02538 sumTwFr = 0.2539 while True:2540 try:2541 hapData['TwinFr'][1][0] = parmDict[pfx+'TwinFr;'+str(it)]2542 if pfx+'TwinFr;'+str(it) in sigDict:2543 TwinFrSig[pfx+'TwinFr;'+str(it)] = sigDict[pfx+'TwinFr;'+str(it)]2544 if it:2545 sumTwFr += hapData['TwinFr'][1][0]2546 it += 12547 except KeyError:2548 break2549 # hapData['TwinFr'][1][0] 2536 if 'Twins' in hapData: 2537 it = 0 2538 sumTwFr = 0. 2539 while True: 2540 try: 2541 hapData['Twins'][it][1][0] = parmDict[pfx+'TwinFr;'+str(it)] 2542 if pfx+'TwinFr;'+str(it) in sigDict: 2543 TwinFrSig[pfx+'TwinFr;'+str(it)] = sigDict[pfx+'TwinFr;'+str(it)] 2544 if it: 2545 sumTwFr += hapData['Twins'][it][1][0] 2546 it += 1 2547 except KeyError: 2548 break 2549 hapData['Twins'][0][1][0] = 1.-sumTwFr 2550 2550 2551 2551 if Print: … … 2609 2609 if pfx+'Flack' in ScalExtSig: 2610 2610 print >>pFile,' Flack parameter : %10.3f, sig %10.3f'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack']) 2611 if pfx+'TwinFr;1' in TwinFrSig:2612 PrintTwin FrAndSig(pfx,hapData['TwinFr'],TwinFrSig)2611 if len(TwinFrSig): 2612 PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig) 2613 2613 2614 2614 ################################################################################ -
trunk/GSASIIstrMath.py
r1886 r1887 918 918 TwinFr = np.array([parmDict[phfx+'TwinFr;'+str(i)] for i in range(len(TwinLaw))]) 919 919 if len(TwinLaw) > 1: 920 TwinFr[0] = 1.-np.sum(TwinFr[1:]) 920 TwinFr[0] = 1.-np.sum(TwinFr[1:]) 921 nTwin = len(TwinLaw) 921 922 nRef = len(refDict['RefList']) 922 923 Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict) … … 931 932 bij = Mast*Uij.T 932 933 dFdvDict = {} 933 dFdfr = np.zeros((nRef,mSize)) 934 dFdx = np.zeros((nRef,mSize,3)) 935 dFdui = np.zeros((nRef,mSize)) 936 dFdua = np.zeros((nRef,mSize,6)) 937 dFdbab = np.zeros((nRef,2)) 938 dFdfl = np.zeros(nRef) 934 dFdfr = np.squeeze(np.zeros((nRef,nTwin,mSize))) 935 dFdx = np.squeeze(np.zeros((nRef,nTwin,mSize,3))) 936 dFdui = np.squeeze(np.zeros((nRef,nTwin,mSize))) 937 dFdua = np.squeeze(np.zeros((nRef,nTwin,mSize,6))) 938 dFdbab = np.squeeze(np.zeros((nRef,nTwin,2))) 939 dFdfl = np.squeeze(np.zeros((nRef,nTwin))) 940 dFdtw = np.zeros((nRef,nTwin)) 939 941 Flack = 1.0 940 942 if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType']: … … 958 960 occ = Mdata*Fdata/len(Uniq) 959 961 biso = -SQfactor*Uisodata[:,np.newaxis] 960 Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)* len(TwinLaw),axis=1).T962 Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1).T 961 963 HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1) 962 964 Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))]) 963 Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),( len(TwinLaw),-1,6)))965 Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6))) 964 966 Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T 965 967 Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT) 966 968 fot = (FF+FP-Bab)*occ*Tcorr 967 969 fotp = FPP*occ*Tcorr 968 if 'T' in calcControls[hfx+'histType']: 969 fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) 970 fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr]) 971 else: 972 fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) 973 fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr]) 970 fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) 971 fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr]) 974 972 fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl 975 973 fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl 976 974 fax = np.array([-fot*sinp,-fotp*cosp]) #positions 977 975 fbx = np.array([fot*cosp,-fotp*sinp]) 978 #sum below is over Uniq - twin effects??976 #sum below is over Uniq 979 977 dfadfr = np.sum(fa/occ,axis=-2) #Fdata != 0 ever avoids /0. problem 980 dfad x = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,np.newaxis],axis=-2)978 dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1) 981 979 dfadui = np.sum(-SQfactor*fa,axis=-2) 982 dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,np.newaxis],axis=-2) 983 dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1) 980 if len(TwinLaw) > 1: 981 dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax[it],-2,-1)[:,:,:,np.newaxis],axis=-2) for it in range(nTwin)]) 982 dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa[it],-2,-1)[:,:,:,np.newaxis],axis=-2) for it in range(nTwin)]) 983 else: 984 dfadx = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,np.newaxis],axis=-2) 985 dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,np.newaxis],axis=-2) 984 986 if not SGData['SGInv']: 985 987 dfbdfr = np.sum(fb/occ,axis=-2) #Fdata != 0 ever avoids /0. problem 986 dfbdx = np.sum(twopi*Uniq*np.swapaxes(fbx,-2,-1)[:,:,:,np.newaxis],axis=2)987 dfbdui = np.sum(-SQfactor*fb,axis=-2)988 dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,np.newaxis],axis=2)989 988 dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1) 990 989 dfadfl = np.sum(-fotp[:,np.newaxis]*sinp) 991 990 dfbdfl = np.sum(fotp[:,np.newaxis]*cosp) 991 dfbdui = np.sum(-SQfactor*fb,axis=-2) 992 if len(TwinLaw) > 1: 993 dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx[it],-2,-1)[:,:,:,np.newaxis],axis=2) for it in range(nTwin)]) 994 dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb[it],-2,-1)[:,:,:,np.newaxis],axis=2) for it in range(nTwin)]) 995 else: 996 dfbdx = np.sum(twopi*Uniq*np.swapaxes(fbx,-2,-1)[:,:,:,np.newaxis],axis=2) 997 dfbdua = np.sum(-Hij*np.swapaxes(fb,-2,-1)[:,:,:,np.newaxis],axis=2) 992 998 else: 993 999 dfbdfr = np.zeros_like(dfadfr) … … 1011 1017 SA = fas[0]-fbs[1] 1012 1018 SB = fbs[0]+fas[1] 1013 dFdfr[iref] = 2.*SA*(dfadfr[0]+dfbdfr[1])*Mdata/len(Uniq)+ \ 1014 2.*SB*(dfbdfr[0]+dfadfr[1])*Mdata/len(Uniq) 1015 dFdx[iref] = 2.*SA*(dfadx[0]+dfbdx[1])+2.*SB*(dfbdx[0]+dfadx[1]) 1016 dFdui[iref] = 2.*SA*(dfadui[0]+dfbdui[1])+2.*SB*(dfbdui[0]+dfadui[1]) 1017 dFdua[iref] = 2.*SA*(dfadua[0]+dfbdua[1])+2.*SB*(dfbdua[0]+dfadua[1]) 1018 dFdfl[iref] = -4.*SA*dfadfl-4.*SB*dfbdfl 1019 if nTwin > 1: 1020 dFdfr[iref] = [2.*SA[it]*(dfadfr[it][0]+dfbdfr[it][1])*Mdata/len(Uniq[it])+ \ 1021 2.*SB[it]*(dfbdfr[it][0]+dfadfr[it][1])*Mdata/len(Uniq[it]) for it in range(nTwin)] 1022 dFdx[iref] = [2.*SA[it]*(dfadx[it][0]+dfbdx[it][1])+2.*SB[it]*(dfbdx[it][0]+dfadx[it][1]) for it in range(nTwin)] 1023 dFdui[iref] = [2.*SA[it]*(dfadui[it][0]+dfbdui[it][1])+2.*SB[it]*(dfbdui[it][0]+dfadui[it][1]) for it in range(nTwin)] 1024 dFdua[iref] = [2.*SA[it]*(dfadua[it][0]+dfbdua[it][1])+2.*SB[it]*(dfbdua[it][0]+dfadua[it][1]) for it in range(nTwin)] 1025 dFdfl[iref] = -SA*dfadfl-SB*dfbdfl 1026 dFdtw[iref] = 2.*SA+2.*SB 1027 else: 1028 dFdfr[iref] = 2.*SA*(dfadfr[0]+dfbdfr[1])*Mdata/len(Uniq)+ \ 1029 2.*SB*(dfbdfr[0]+dfadfr[1])*Mdata/len(Uniq) 1030 dFdx[iref] = 2.*SA*(dfadx[0]+dfbdx[1])+2.*SB*(dfbdx[0]+dfadx[1]) 1031 dFdui[iref] = 2.*SA*(dfadui[0]+dfbdui[1])+2.*SB*(dfbdui[0]+dfadui[1]) 1032 dFdua[iref] = 2.*SA*(dfadua[0]+dfbdua[1])+2.*SB*(dfbdua[0]+dfadua[1]) 1033 dFdfl[iref] = -SA*dfadfl-SB*dfbdfl 1019 1034 dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \ 1020 1035 2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T 1021 1036 1022 1037 #loop over atoms - each dict entry is list of derivatives for all the reflections 1023 for i in range(len(Mdata)): 1024 dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i] 1025 dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i] 1026 dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i] 1027 dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i] 1028 dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i] 1029 dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i] 1030 dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i] 1031 dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i] 1032 dFdvDict[pfx+'AU12:'+str(i)] = 0.5*dFdua.T[3][i] 1033 dFdvDict[pfx+'AU13:'+str(i)] = 0.5*dFdua.T[4][i] 1034 dFdvDict[pfx+'AU23:'+str(i)] = 0.5*dFdua.T[5][i] 1038 if nTwin > 1: 1039 for i in range(len(Mdata)): 1040 dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,np.newaxis],axis=0) 1041 dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,np.newaxis],axis=0) 1042 dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,np.newaxis],axis=0) 1043 dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,np.newaxis],axis=0) 1044 dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,np.newaxis],axis=0) 1045 dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,np.newaxis],axis=0) 1046 dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,np.newaxis],axis=0) 1047 dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,np.newaxis],axis=0) 1048 dFdvDict[pfx+'AU12:'+str(i)] = np.sum(0.5*dFdua.T[3][i]*TwinFr[:,np.newaxis],axis=0) 1049 dFdvDict[pfx+'AU13:'+str(i)] = np.sum(0.5*dFdua.T[4][i]*TwinFr[:,np.newaxis],axis=0) 1050 dFdvDict[pfx+'AU23:'+str(i)] = np.sum(0.5*dFdua.T[5][i]*TwinFr[:,np.newaxis],axis=0) 1051 dFdvDict[phfx+'Flack'] = np.sum(dFdfl.T*TwinFr[:,np.newaxis],axis=0) 1052 else: 1053 for i in range(len(Mdata)): 1054 dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i] 1055 dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i] 1056 dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i] 1057 dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i] 1058 dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i] 1059 dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i] 1060 dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i] 1061 dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i] 1062 dFdvDict[pfx+'AU12:'+str(i)] = 0.5*dFdua.T[3][i] 1063 dFdvDict[pfx+'AU13:'+str(i)] = 0.5*dFdua.T[4][i] 1064 dFdvDict[pfx+'AU23:'+str(i)] = 0.5*dFdua.T[5][i] 1065 dFdvDict[phfx+'Flack'] = dFdfl.T 1035 1066 dFdvDict[phfx+'BabA'] = dFdbab.T[0] 1036 1067 dFdvDict[phfx+'BabU'] = dFdbab.T[1] 1037 dFdvDict[phfx+'Flack'] = dFdfl.T 1068 if nTwin > 1: 1069 for i in range(nTwin-1): #skip the base twin element 1070 dFdvDict[phfx+'TwinFr;'+str(i+1)] = dFdtw.T[i+1] 1038 1071 return dFdvDict 1039 1072
Note: See TracChangeset
for help on using the changeset viewer.