Changeset 1919
- Timestamp:
- Jul 6, 2015 8:41:52 AM (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIstrMath.py
r1918 r1919 966 966 Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) 967 967 FF = refDict['FF']['FF'][iref].T[Tindx].T 968 Uniq = np.inner(H,SGMT) # array( 3,nSGOp) or968 Uniq = np.inner(H,SGMT) # array(nSGOp,3) or (nTwin,nSGOp,3) 969 969 Phi = np.inner(H,SGT) 970 970 phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T 971 971 sinp = np.sin(phase) 972 972 cosp = np.cos(phase) 973 occ = Mdata*Fdata/len( Uniq)973 occ = Mdata*Fdata/len(SGT) 974 974 biso = -SQfactor*Uisodata[:,np.newaxis] 975 Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1) .T975 Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*nTwin,axis=1) 976 976 HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1) 977 977 Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))]) 978 978 Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6))) 979 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) .T980 Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*occ979 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) 980 Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ 981 981 fot = (FF+FP-Bab)*Tcorr 982 982 fotp = FPP*Tcorr 983 fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) 983 fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nTwin,nEqv,nAtoms) 984 984 fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr]) 985 fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl 985 fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl array(2,nTwins) 986 986 fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl 987 fax = np.array([-fot*sinp,-fotp*cosp]) #positions 987 fax = np.array([-fot*sinp,-fotp*cosp]) #positions array(2,ntwi,nEqv,nAtoms) 988 988 fbx = np.array([fot*cosp,-fotp*sinp]) 989 989 #sum below is over Uniq 990 dfadfr = np.sum(fa/occ,axis=-2) # Fdata != 0 avoids /0. problem990 dfadfr = np.sum(fa/occ,axis=-2) #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem 991 991 dfadba = np.sum(-cosp*Tcorr[:,np.newaxis],axis=1) 992 992 dfadui = np.sum(-SQfactor*fa,axis=-2) 993 if len(TwinLaw)> 1:993 if nTwin > 1: 994 994 dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)]) 995 995 dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,np.newaxis],axis=-2) for it in range(nTwin)]) 996 # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6) 996 997 else: 997 998 dfadx = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,np.newaxis],axis=-2) 998 999 dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,np.newaxis],axis=-2) 1000 # array(2,nAtom,3) & array(2,nAtom,6) 999 1001 if not SGData['SGInv']: 1000 1002 dfbdfr = np.sum(fb/occ,axis=-2) #Fdata != 0 avoids /0. problem … … 1031 1033 SB = fbs[0]+fas[1] 1032 1034 if nTwin > 1: 1033 dFdfr[iref] = TwMask[:,np.newaxis]*[SA[it]*(dfadfr[0][it]+dfbdfr[1][it])*Mdata/len(Uniq[it])+ \1034 SB[it]*(dfbdfr[0][it]+dfadfr[1][it])*Mdata/len(Uniq[it]) for it in range(nTwin)]1035 dFdx[iref] = TwMask[:,np.newaxis,np.newaxis]*[SA[it]*(dfadx[it][0]+dfbdx[it][1])+SB[it]*(dfbdx[it][0]+dfadx[it][1]) for it in range(nTwin)]1036 dFdui[iref] = TwMask[:,np.newaxis]*[SA[it]*(dfadui[0][it]+dfbdui[1][it])+SB[it]*(dfbdui[0][it]+dfadui[1][it]) for it in range(nTwin)]1037 dFdua[iref] = TwMask[:,np.newaxis,np.newaxis]*[SA[it]*(dfadua[it][0]+dfbdua[it][1])+SB[it]*(dfbdua[it][0]+dfadua[it][1]) for it in range(nTwin)]1038 dFdtw[iref] = TwMask*np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**21035 dFdfr[iref] = [2.*TwMask[it]*SA[it]*(dfadfr[0][it]+dfbdfr[1][it])*Mdata/len(Uniq[it])+ \ 1036 2.*TwMask[it]*SB[it]*(dfbdfr[0][it]+dfadfr[1][it])*Mdata/len(Uniq[it]) for it in range(nTwin)] 1037 dFdx[iref] = [2.*TwMask[it]*SA[it]*(dfadx[it][0]+dfbdx[it][1])+2.*TwMask[it]*SB[it]*(dfbdx[it][0]+dfadx[it][1]) for it in range(nTwin)] 1038 dFdui[iref] = [2.*TwMask[it]*SA[it]*(dfadui[0][it]+dfbdui[1][it])+2.*TwMask[it]*SB[it]*(dfbdui[0][it]+dfadui[1][it]) for it in range(nTwin)] 1039 dFdua[iref] = [2.*TwMask[it]*SA[it]*(dfadua[it][0]+dfbdua[it][1])+2.*TwMask[it]*SB[it]*(dfbdua[it][0]+dfadua[it][1]) for it in range(nTwin)] 1040 dFdtw[iref] = 2.*TwMask*np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2 1039 1041 else: #these are good for no twin single crystals 1040 1042 dFdfr[iref] = 2.*SA*(dfadfr[0]+dfbdfr[1])*Mdata/len(Uniq)+ \ 1041 2.*SB*(dfbdfr[0]+dfadfr[1])*Mdata/len(Uniq) 1042 dFdx[iref] = 2.*SA*(dfadx[0]+dfbdx[1])+2.*SB*(dfbdx[0]+dfadx[1]) 1043 dFdui[iref] = 2.*SA*(dfadui[0]+dfbdui[1])+2.*SB*(dfbdui[0]+dfadui[1]) 1044 dFdua[iref] = 2.*SA*(dfadua[0]+dfbdua[1])+2.*SB*(dfbdua[0]+dfadua[1]) 1045 dFdfl[iref] = -SA*dfadfl-SB*dfbdfl 1043 2.*SB*(dfbdfr[0]+dfadfr[1])*Mdata/len(Uniq) #array(nRef,nAtom) 1044 dFdx[iref] = 2.*SA*(dfadx[0]+dfbdx[1])+2.*SB*(dfbdx[0]+dfadx[1]) #array(nRef,nAtom,3) 1045 dFdui[iref] = 2.*SA*(dfadui[0]+dfbdui[1])+2.*SB*(dfbdui[0]+dfadui[1]) #array(nRef,nAtom) 1046 dFdua[iref] = 2.*SA*(dfadua[0]+dfbdua[1])+2.*SB*(dfbdua[0]+dfadua[1]) #array(nRef,nAtom,6) 1047 dFdfl[iref] = -SA*dfadfl-SB*dfbdfl #array(nRef,) 1046 1048 dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \ 1047 1049 2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T 1048 1049 1050 #loop over atoms - each dict entry is list of derivatives for all the reflections 1050 1051 if nTwin > 1: 1051 1052 for i in range(len(Mdata)): #these all OK? 1052 1053 dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,np.newaxis],axis=0) 1053 dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,np.newaxis],axis=0) /2.1054 dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,np.newaxis],axis=0) /2.1055 dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,np.newaxis],axis=0) /2.1056 dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,np.newaxis],axis=0) /2.1057 dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,np.newaxis],axis=0) /2.1058 dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,np.newaxis],axis=0) /2.1059 dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,np.newaxis],axis=0) /2.1060 dFdvDict[pfx+'AU12:'+str(i)] = np.sum(0.5*dFdua.T[3][i]*TwinFr[:,np.newaxis],axis=0) /2.1061 dFdvDict[pfx+'AU13:'+str(i)] = np.sum(0.5*dFdua.T[4][i]*TwinFr[:,np.newaxis],axis=0) /2.1062 dFdvDict[pfx+'AU23:'+str(i)] = np.sum(0.5*dFdua.T[5][i]*TwinFr[:,np.newaxis],axis=0) /2.1054 dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,np.newaxis],axis=0) 1055 dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,np.newaxis],axis=0) 1056 dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,np.newaxis],axis=0) 1057 dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,np.newaxis],axis=0) 1058 dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,np.newaxis],axis=0) 1059 dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,np.newaxis],axis=0) 1060 dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,np.newaxis],axis=0) 1061 dFdvDict[pfx+'AU12:'+str(i)] = np.sum(0.5*dFdua.T[3][i]*TwinFr[:,np.newaxis],axis=0) 1062 dFdvDict[pfx+'AU13:'+str(i)] = np.sum(0.5*dFdua.T[4][i]*TwinFr[:,np.newaxis],axis=0) 1063 dFdvDict[pfx+'AU23:'+str(i)] = np.sum(0.5*dFdua.T[5][i]*TwinFr[:,np.newaxis],axis=0) 1063 1064 else: 1064 1065 for i in range(len(Mdata)):
Note: See TracChangeset
for help on using the changeset viewer.