Changeset 2097
- Timestamp:
- Dec 18, 2015 10:58:57 AM (7 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIphsGUI.py
r2077 r2097 1392 1392 name = atomData[r][c] 1393 1393 if len(name) in [2,4]: 1394 atomData[r][c-1] = name[:2]+' (%d)'%(r+1)1394 atomData[r][c-1] = name[:2]+'%d'%(r+1) 1395 1395 else: 1396 atomData[r][c-1] = name[:1]+' (%d)'%(r+1)1396 atomData[r][c-1] = name[:1]+'%d'%(r+1) 1397 1397 PE.Destroy() 1398 1398 SetupGeneral() … … 2053 2053 else: 2054 2054 data['Atoms'] = result 2055 Atoms.ClearSelection() 2056 data['Drawing']['Atoms'] = [] 2055 2057 OnReloadDrawAtoms(event) 2056 2058 FillAtomsGrid(Atoms) … … 2064 2066 def OnDistAnglePrt(event): 2065 2067 'save distances and angles to a file' 2066 fp = file(os.path.abspath(os.path.splitext(G2frame.GSASprojectfile 2067 )[0]+'.disagl'),'w') 2068 fp = file(os.path.abspath(os.path.splitext(G2frame.GSASprojectfile)[0]+'.disagl'),'w') 2068 2069 OnDistAngle(event,fp=fp) 2069 2070 fp.close() -
trunk/GSASIIstrMath.py
r2096 r2097 648 648 ''' Compute structure factors for all h,k,l for phase 649 649 puts the result, F^2, in each ref[8] in refList 650 operates on blocks of 100 reflections for speed 650 651 input: 651 652 652 653 :param dict refDict: where 653 'RefList' list where each ref = h,k,l, m,d,...654 'RefList' list where each ref = h,k,l,it,d,... 654 655 'FF' dict of form factors - filed in below 655 656 :param np.array G: reciprocal metric tensor … … 769 770 770 771 def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict): 771 '''Compute structure factor derivatives on blocks ofreflections - keep as it works for twins772 '''Compute structure factor derivatives on single reflections - keep as it works for twins 772 773 but is slower for powders/nontwins 773 774 ''' … … 962 963 963 964 def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict): 964 '''Compute structure factor derivatives on blocks of reflections - works for powders/nontwins965 '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only 965 966 faster than StructureFactorDerv 967 ''' 968 phfx = pfx.split(':')[0]+hfx 969 ast = np.sqrt(np.diag(G)) 970 Mast = twopisq*np.multiply.outer(ast,ast) 971 SGMT = np.array([ops[0].T for ops in SGData['SGOps']]) 972 SGT = np.array([ops[1] for ops in SGData['SGOps']]) 973 FFtables = calcControls['FFtables'] 974 BLtables = calcControls['BLtables'] 975 nRef = len(refDict['RefList']) 976 Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict) 977 mSize = len(Mdata) 978 FF = np.zeros(len(Tdata)) 979 if 'NC' in calcControls[hfx+'histType']: 980 FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam']) 981 elif 'X' in calcControls[hfx+'histType']: 982 FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata]) 983 FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata]) 984 Uij = np.array(G2lat.U6toUij(Uijdata)) 985 bij = Mast*Uij.T 986 dFdvDict = {} 987 dFdfr = np.zeros((nRef,mSize)) 988 dFdx = np.zeros((nRef,mSize,3)) 989 dFdui = np.zeros((nRef,mSize)) 990 dFdua = np.zeros((nRef,mSize,6)) 991 dFdbab = np.zeros((nRef,2)) 992 dFdfl = np.zeros((nRef)) 993 dFdtw = np.zeros((nRef)) 994 Flack = 1.0 995 if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict: 996 Flack = 1.-2.*parmDict[phfx+'Flack'] 997 time0 = time.time() 998 nref = len(refDict['RefList'])/100 999 #reflection processing begins here - big arrays! 1000 iBeg = 0 1001 blkSize = 32 #no. of reflections in a block - optimized for speed 1002 while iBeg < nRef: 1003 iFin = min(iBeg+blkSize,nRef) 1004 refl = refDict['RefList'][iBeg:iFin] #array(blkSize,nItems) 1005 H = refl.T[:3] 1006 SQ = 1./(2.*refl.T[4])**2 # or (sin(theta)/lambda)**2 1007 SQfactor = 8.0*SQ*np.pi**2 1008 if 'T' in calcControls[hfx+'histType']: 1009 if 'P' in calcControls[hfx+'histType']: 1010 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14]) 1011 else: 1012 FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12]) 1013 FP = np.repeat(FP.T,len(SGT),axis=0) 1014 FPP = np.repeat(FPP.T,len(SGT),axis=0) 1015 dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor) 1016 Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)) 1017 Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) 1018 FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0) 1019 Uniq = np.inner(H.T,SGMT) # array(nSGOp,3) or (nTwin,nSGOp,3) 1020 Phi = np.inner(H.T,SGT) 1021 phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T 1022 sinp = np.sin(phase) #refBlk x nOps x nAtoms 1023 cosp = np.cos(phase) 1024 occ = Mdata*Fdata/len(SGT) 1025 biso = -SQfactor*Uisodata[:,nxs] 1026 Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T 1027 HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1) 1028 Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T 1029 Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT) 1030 Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))]) 1031 Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,len(SGT),6)) 1032 fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr 1033 fotp = FPP*Tcorr 1034 if 'T' in calcControls[hfx+'histType']: 1035 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr]) 1036 fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr]) 1037 else: 1038 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) 1039 fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr]) 1040 # GSASIIpath.IPyBreak() 1041 fas = np.sum(np.sum(fa,axis=-1),axis=-1) #real sum over atoms & unique hkl array(2,refBlk,nTwins) 1042 fbs = np.sum(np.sum(fb,axis=-1),axis=-1) #imag sum over atoms & uniq hkl 1043 fax = np.array([-fot*sinp,-fotp*cosp]) #positions array(2,refBlk,ntwi,nEqv,nAtoms) 1044 fbx = np.array([fot*cosp,-fotp*sinp]) 1045 #sum below is over Uniq 1046 dfadfr = np.sum(fa/occ,axis=-2) #array(2,refBlk,ntwin,nAtom) Fdata != 0 avoids /0. problem 1047 dfadba = np.sum(-cosp*Tcorr,axis=-2) #array(refBlk,nAtom) 1048 dfadx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[:,:,:,:,nxs],axis=-2) 1049 dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fa,axis=-2) #array(Ops,refBlk,nTw,nAtoms) 1050 dfadua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[:,:,:,:,nxs],axis=-2) 1051 # array(2,refBlk,nAtom,3) & array(2,refBlk,nAtom,6) 1052 if not SGData['SGInv']: 1053 dfbdfr = np.sum(fb/occ,axis=-2) #Fdata != 0 avoids /0. problem 1054 dfbdba = np.sum(-sinp*Tcorr,axis=-2) 1055 dfadfl = np.sum(np.sum(-FPP*Tcorr*sinp,axis=-1),axis=-1) 1056 dfbdfl = np.sum(np.sum(FPP*Tcorr*cosp,axis=-1),axis=-1) 1057 dfbdx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[:,:,:,:,nxs],axis=-2) 1058 dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fb,axis=-2) 1059 dfbdua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[:,:,:,:,nxs],axis=-2) 1060 else: 1061 dfbdfr = np.zeros_like(dfadfr) 1062 dfbdx = np.zeros_like(dfadx) 1063 dfbdui = np.zeros_like(dfadui) 1064 dfbdua = np.zeros_like(dfadua) 1065 dfbdba = np.zeros_like(dfadba) 1066 dfadfl = 0.0 1067 dfbdfl = 0.0 1068 #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3! 1069 SA = fas[0]+fas[1] 1070 SB = fbs[0]+fbs[1] 1071 if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro 1072 dFdfr[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadfr[0]+fas[1,:,nxs]*dfadfr[1])*Mdata/len(SGMT)+ \ 1073 2.*(fbs[0,:,nxs]*dfbdfr[0]-fbs[1,:,nxs]*dfbdfr[1])*Mdata/len(SGMT) 1074 dFdx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadx[0]+fas[1,:,nxs,nxs]*dfadx[1])+ \ 1075 2.*(fbs[0,:,nxs,nxs]*dfbdx[0]+fbs[1,:,nxs,nxs]*dfbdx[1]) 1076 dFdui[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadui[0]+fas[1,:,nxs]*dfadui[1])+ \ 1077 2.*(fbs[0,:,nxs]*dfbdui[0]-fbs[1,:,nxs]*dfbdui[1]) 1078 dFdua[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadua[0]+fas[1,:,nxs,nxs]*dfadua[1])+ \ 1079 2.*(fbs[0,:,nxs,nxs]*dfbdua[0]+fbs[1,:,nxs,nxs]*dfbdua[1]) 1080 else: 1081 dFdfr[iBeg:iFin] = (2.*SA[:,nxs]*(dfadfr[0]+dfadfr[1])+2.*SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(SGMT) 1082 dFdx[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+2.*SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1]) 1083 dFdui[iBeg:iFin] = 2.*SA[:,nxs]*(dfadui[0]+dfadui[1])+2.*SB[:,nxs]*(dfbdui[0]+dfbdui[1]) 1084 dFdua[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+2.*SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1]) 1085 dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl #array(nRef,) 1086 dFdbab[iBeg:iFin] = 2.*(fas[0,nxs]*np.array([np.sum(dfadba.T*dBabdA,axis=0),np.sum(-dfadba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])+ \ 1087 fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T 1088 # GSASIIpath.IPyBreak() 1089 iBeg += blkSize 1090 print ' %d derivative time %.4f\r'%(nRef,time.time()-time0) 1091 #loop over atoms - each dict entry is list of derivatives for all the reflections 1092 for i in range(len(Mdata)): 1093 dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i] 1094 dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i] 1095 dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i] 1096 dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i] 1097 dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i] 1098 dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i] 1099 dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i] 1100 dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i] 1101 dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i] 1102 dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i] 1103 dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i] 1104 dFdvDict[phfx+'Flack'] = 4.*dFdfl.T 1105 dFdvDict[phfx+'BabA'] = dFdbab.T[0] 1106 dFdvDict[phfx+'BabU'] = dFdbab.T[1] 1107 return dFdvDict 1108 1109 def StructureFactorDervTw(refDict,G,hfx,pfx,SGData,calcControls,parmDict): 1110 '''Compute structure factor derivatives on blocks of reflections - for twins only 1111 faster than StructureFactorDerv - not working yet 966 1112 ''' 967 1113 phfx = pfx.split(':')[0]+hfx … … 993 1139 bij = Mast*Uij.T 994 1140 dFdvDict = {} 995 if nTwin > 1: 996 dFdfr = np.zeros((nRef,nTwin,mSize)) 997 dFdx = np.zeros((nRef,nTwin,mSize,3)) 998 dFdui = np.zeros((nRef,nTwin,mSize)) 999 dFdua = np.zeros((nRef,nTwin,mSize,6)) 1000 dFdbab = np.zeros((nRef,nTwin,2)) 1001 dFdfl = np.zeros((nRef,nTwin)) 1002 dFdtw = np.zeros((nRef,nTwin)) 1003 else: 1004 dFdfr = np.zeros((nRef,mSize)) 1005 dFdx = np.zeros((nRef,mSize,3)) 1006 dFdui = np.zeros((nRef,mSize)) 1007 dFdua = np.zeros((nRef,mSize,6)) 1008 dFdbab = np.zeros((nRef,2)) 1009 dFdfl = np.zeros((nRef)) 1010 dFdtw = np.zeros((nRef)) 1141 dFdfr = np.zeros((nRef,nTwin,mSize)) 1142 dFdx = np.zeros((nRef,nTwin,mSize,3)) 1143 dFdui = np.zeros((nRef,nTwin,mSize)) 1144 dFdua = np.zeros((nRef,nTwin,mSize,6)) 1145 dFdbab = np.zeros((nRef,nTwin,2)) 1146 dFdfl = np.zeros((nRef,nTwin)) 1147 dFdtw = np.zeros((nRef,nTwin)) 1011 1148 Flack = 1.0 1012 1149 if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict: … … 1023 1160 H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(3,nTwins) or (3) 1024 1161 TwMask = np.any(H,axis=-1) 1025 if TwinLaw.shape[0] > 1 and TwDict: 1026 for ir in range(blkSize): 1027 iref = ir+iBeg 1028 if iref in TwDict: 1029 for i in TwDict[iref]: 1030 for n in range(NTL): 1031 H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM]) 1032 TwMask = np.any(H,axis=-1) 1162 for ir in range(blkSize): 1163 iref = ir+iBeg 1164 if iref in TwDict: 1165 for i in TwDict[iref]: 1166 for n in range(NTL): 1167 H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM]) 1168 TwMask = np.any(H,axis=-1) 1033 1169 SQ = 1./(2.*refl.T[4])**2 # or (sin(theta)/lambda)**2 1034 1170 SQfactor = 8.0*SQ*np.pi**2 … … 1056 1192 Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT) 1057 1193 Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))]) 1058 if nTwin > 1: 1059 Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,nTwin,len(SGT),6)) 1060 else: 1061 Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,len(SGT),6)) 1194 Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,nTwin,len(SGT),6)) 1062 1195 fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr 1063 1196 fotp = FPP*Tcorr … … 1076 1209 dfadfr = np.sum(fa/occ,axis=-2) #array(2,refBlk,ntwin,nAtom) Fdata != 0 avoids /0. problem 1077 1210 dfadba = np.sum(-cosp*Tcorr,axis=-2) #array(refBlk,nAtom) 1078 if nTwin > 1: 1079 dfadfr = np.swapaxes(dfadfr,0,1) 1080 dfadx = np.swapaxes(np.array([np.sum(twopi*Uniq[it,nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[it,:,:,:,:,nxs],axis=-2) for it in range(nTwin)]),0,1) 1081 dfadui = np.swapaxes(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fa,axis=-2),0,1) #array(Ops,refBlk,nTw,nAtoms) 1082 dfadua = np.swapaxes(np.array([np.sum(-Hij[it,nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[it,:,:,:,:,nxs],axis=-2) for it in range(nTwin)]),0,1) 1083 # array(nTwin,2,refBlk,nAtom,3) & array(nTwin,2,refBlk,nAtom,6) 1084 else: 1085 dfadx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[:,:,:,:,nxs],axis=-2) 1086 dfadui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fa,axis=-2) #array(Ops,refBlk,nTw,nAtoms) 1087 dfadua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[:,:,:,:,nxs],axis=-2) 1088 # array(2,refBlk,nAtom,3) & array(2,refBlk,nAtom,6) 1211 dfadfr = np.swapaxes(dfadfr,0,1) 1212 dfadx = np.swapaxes(np.array([np.sum(twopi*Uniq[it,nxs,:,nxs,:,:]*np.swapaxes(fax,-2,-1)[it,:,:,:,:,nxs],axis=-2) for it in range(nTwin)]),0,1) 1213 dfadui = np.swapaxes(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fa,axis=-2),0,1) #array(Ops,refBlk,nTw,nAtoms) 1214 dfadua = np.swapaxes(np.array([np.sum(-Hij[it,nxs,:,nxs,:,:]*np.swapaxes(fa,-2,-1)[it,:,:,:,:,nxs],axis=-2) for it in range(nTwin)]),0,1) 1215 # array(nTwin,2,refBlk,nAtom,3) & array(nTwin,2,refBlk,nAtom,6) 1089 1216 if not SGData['SGInv']: 1090 1217 dfbdfr = np.sum(fb/occ,axis=-2) #Fdata != 0 avoids /0. problem 1091 1218 dfbdba = np.sum(-sinp*Tcorr,axis=-2) 1092 if len(TwinLaw) > 1: 1093 dfbdfr = np.swapaxes(dfbdfr,0,1) 1094 dfbdx = np.swapaxes(np.array([np.sum(twopi*Uniq[it,nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[it,:,:,:,:,nxs],axis=-2) for it in range(nTwin)]),0,1) 1095 dfbdui = np.swapaxes(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fb,axis=-2),0,1) #array(Ops,refBlk,nTw,nAtoms) 1096 dfbdua = np.swapaxes(np.array([np.sum(-Hij[it,nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[it,:,:,:,:,nxs],axis=-2) for it in range(nTwin)]),0,1) 1097 else: 1098 dfadfl = np.sum(np.sum(-FPP*Tcorr*sinp,axis=-1),axis=-1) 1099 dfbdfl = np.sum(np.sum(FPP*Tcorr*cosp,axis=-1),axis=-1) 1100 dfbdx = np.sum(twopi*Uniq[nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[:,:,:,:,nxs],axis=-2) 1101 dfbdui = np.sum(-SQfactor[nxs,:,nxs,nxs]*fb,axis=-2) 1102 dfbdua = np.sum(-Hij[nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[:,:,:,:,nxs],axis=-2) 1219 dfbdfr = np.swapaxes(dfbdfr,0,1) 1220 dfbdx = np.swapaxes(np.array([np.sum(twopi*Uniq[it,nxs,:,nxs,:,:]*np.swapaxes(fbx,-2,-1)[it,:,:,:,:,nxs],axis=-2) for it in range(nTwin)]),0,1) 1221 dfbdui = np.swapaxes(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fb,axis=-2),0,1) #array(Ops,refBlk,nTw,nAtoms) 1222 dfbdua = np.swapaxes(np.array([np.sum(-Hij[it,nxs,:,nxs,:,:]*np.swapaxes(fb,-2,-1)[it,:,:,:,:,nxs],axis=-2) for it in range(nTwin)]),0,1) 1103 1223 else: 1104 1224 dfbdfr = np.zeros_like(dfadfr) … … 1112 1232 SA = fas[0]+fas[1] 1113 1233 SB = fbs[0]+fbs[1] 1114 if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro 1115 dFdfr[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadfr[0]+fas[1,:,nxs]*dfadfr[1])*Mdata/len(SGMT)+ \ 1116 2.*(fbs[0,:,nxs]*dfbdfr[0]-fbs[1,:,nxs]*dfbdfr[1])*Mdata/len(SGMT) 1117 dFdx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadx[0]+fas[1,:,nxs,nxs]*dfadx[1])+ \ 1118 2.*(fbs[0,:,nxs,nxs]*dfbdx[0]+fbs[1,:,nxs,nxs]*dfbdx[1]) 1119 dFdui[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadui[0]+fas[1,:,nxs]*dfadui[1])+ \ 1120 2.*(fbs[0,:,nxs]*dfbdui[0]-fbs[1,:,nxs]*dfbdui[1]) 1121 dFdua[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadua[0]+fas[1,:,nxs,nxs]*dfadua[1])+ \ 1122 2.*(fbs[0,:,nxs,nxs]*dfbdua[0]+fbs[1,:,nxs,nxs]*dfbdua[1]) 1123 else: 1124 if nTwin > 1: 1125 dFdfr[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs]*(SA[:,it,nxs]*(np.sum(dfadfr[:,:,it],axis=1))+ \ 1126 SB[:,it,nxs]*(np.sum(dfbdfr[:,:,it],axis=1)))*Mdata/len(SGMT) for it in range(nTwin)]),0,1) 1127 dFdx[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs,nxs]*(SA[:,it,nxs,nxs]*np.sum(dfadx[:,:,it],axis=1)+ \ 1128 SB[:,it,nxs,nxs]*np.sum(dfbdx[:,:,it],axis=1)) for it in range(nTwin)]),0,1) 1129 dFdui[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs]*(SA[:,it,nxs]*np.sum(dfadui[:,:,it],axis=1)+ \ 1130 SB[:,it,nxs]*np.sum(dfbdui[:,:,it],axis=1)) for it in range(nTwin)]),0,1) 1131 dFdua[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs,nxs]*(SA[:,it,nxs,nxs]*np.sum(dfadua[:,:,it],axis=1)+ \ 1132 SB[:,it,nxs,nxs]*np.sum(dfbdua[:,:,it],axis=1)) for it in range(nTwin)]),0,1) 1133 dFdtw[iBeg:iFin] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2 1234 dFdfr[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs]*(SA[:,it,nxs]*(np.sum(dfadfr[:,:,it],axis=1))+ \ 1235 SB[:,it,nxs]*(np.sum(dfbdfr[:,:,it],axis=1)))*Mdata/len(SGMT) for it in range(nTwin)]),0,1) 1236 dFdx[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs,nxs]*(SA[:,it,nxs,nxs]*np.sum(dfadx[:,:,it],axis=1)+ \ 1237 SB[:,it,nxs,nxs]*np.sum(dfbdx[:,:,it],axis=1)) for it in range(nTwin)]),0,1) 1238 dFdui[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs]*(SA[:,it,nxs]*np.sum(dfadui[:,:,it],axis=1)+ \ 1239 SB[:,it,nxs]*np.sum(dfbdui[:,:,it],axis=1)) for it in range(nTwin)]),0,1) 1240 dFdua[iBeg:iFin] = np.swapaxes(np.array([2.*TwMask[:,it,nxs,nxs]*(SA[:,it,nxs,nxs]*np.sum(dfadua[:,:,it],axis=1)+ \ 1241 SB[:,it,nxs,nxs]*np.sum(dfbdua[:,:,it],axis=1)) for it in range(nTwin)]),0,1) 1242 dFdtw[iBeg:iFin] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2 1134 1243 1135 else: #these are good for no twin single crystals 1136 dFdfr[iBeg:iFin] = (2.*SA[:,nxs]*(dfadfr[0]+dfadfr[1])+2.*SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(SGMT) 1137 dFdx[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadx[0]+dfadx[1])+2.*SB[:,nxs,nxs]*(dfbdx[0]+dfbdx[1]) 1138 dFdui[iBeg:iFin] = 2.*SA[:,nxs]*(dfadui[0]+dfadui[1])+2.*SB[:,nxs]*(dfbdui[0]+dfbdui[1]) 1139 dFdua[iBeg:iFin] = 2.*SA[:,nxs,nxs]*(dfadua[0]+dfadua[1])+2.*SB[:,nxs,nxs]*(dfbdua[0]+dfbdua[1]) 1140 dFdfl[iBeg:iFin] = -SA*dfadfl-SB*dfbdfl #array(nRef,) 1141 dFdbab[iBeg:iFin] = 2.*(fas[0,nxs]*np.array([np.sum(dfadba.T*dBabdA,axis=0),np.sum(-dfadba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])+ \ 1142 fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T 1244 # dFdbab[iBeg:iFin] = 2.*(fas[0,nxs]*np.array([np.sum(dfadba.T*dBabdA,axis=0),np.sum(-dfadba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])+ \ 1245 # fbs[0,nxs]*np.array([np.sum(dfbdba.T*dBabdA,axis=0),np.sum(-dfbdba.T*parmDict[phfx+'BabA']*SQfactor*dBabdA,axis=0)])).T 1143 1246 # GSASIIpath.IPyBreak() 1144 1247 iBeg += blkSize 1145 1248 print ' %d derivative time %.4f\r'%(nRef,time.time()-time0) 1146 1249 #loop over atoms - each dict entry is list of derivatives for all the reflections 1147 if nTwin > 1: 1148 for i in range(len(Mdata)): #these all OK? 1149 dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0) 1150 dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0) 1151 dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0) 1152 dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0) 1153 dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0) 1154 dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0) 1155 dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0) 1156 dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0) 1157 dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0) 1158 dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0) 1159 dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0) 1160 else: 1161 for i in range(len(Mdata)): 1162 dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i] 1163 dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i] 1164 dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i] 1165 dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i] 1166 dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i] 1167 dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i] 1168 dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i] 1169 dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i] 1170 dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i] 1171 dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i] 1172 dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i] 1173 dFdvDict[phfx+'Flack'] = 4.*dFdfl.T 1250 for i in range(len(Mdata)): #these all OK? 1251 dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0) 1252 dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0) 1253 dFdvDict[pfx+'dAy:'+str(i)] = np.sum(dFdx.T[1][i]*TwinFr[:,nxs],axis=0) 1254 dFdvDict[pfx+'dAz:'+str(i)] = np.sum(dFdx.T[2][i]*TwinFr[:,nxs],axis=0) 1255 dFdvDict[pfx+'AUiso:'+str(i)] = np.sum(dFdui.T[i]*TwinFr[:,nxs],axis=0) 1256 dFdvDict[pfx+'AU11:'+str(i)] = np.sum(dFdua.T[0][i]*TwinFr[:,nxs],axis=0) 1257 dFdvDict[pfx+'AU22:'+str(i)] = np.sum(dFdua.T[1][i]*TwinFr[:,nxs],axis=0) 1258 dFdvDict[pfx+'AU33:'+str(i)] = np.sum(dFdua.T[2][i]*TwinFr[:,nxs],axis=0) 1259 dFdvDict[pfx+'AU12:'+str(i)] = 2.*np.sum(dFdua.T[3][i]*TwinFr[:,nxs],axis=0) 1260 dFdvDict[pfx+'AU13:'+str(i)] = 2.*np.sum(dFdua.T[4][i]*TwinFr[:,nxs],axis=0) 1261 dFdvDict[pfx+'AU23:'+str(i)] = 2.*np.sum(dFdua.T[5][i]*TwinFr[:,nxs],axis=0) 1262 dFdvDict[phfx+'Flack'] = 4.*dFdfl.T 1174 1263 dFdvDict[phfx+'BabA'] = dFdbab.T[0] 1175 1264 dFdvDict[phfx+'BabU'] = dFdbab.T[1] 1176 if nTwin > 1: 1177 for i in range(nTwin): 1178 dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i] 1265 for i in range(nTwin): 1266 dFdvDict[phfx+'TwinFr:'+str(i)] = dFdtw.T[i] 1179 1267 return dFdvDict 1180 1268 … … 1183 1271 Compute super structure factors for all h,k,l,m for phase 1184 1272 puts the result, F^2, in each ref[8+im] in refList 1273 works on blocks of 100 reflections for speed 1185 1274 input: 1186 1275 1187 1276 :param dict refDict: where 1188 'RefList' list where each ref = h,k,l,m, d,...1277 'RefList' list where each ref = h,k,l,m,it,d,... 1189 1278 'FF' dict of form factors - filed in below 1190 1279 :param np.array G: reciprocal metric tensor … … 1318 1407 1319 1408 def SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict): 1409 'Needs a doc string - no twins' 1410 phfx = pfx.split(':')[0]+hfx 1411 ast = np.sqrt(np.diag(G)) 1412 Mast = twopisq*np.multiply.outer(ast,ast) 1413 SGInv = SGData['SGInv'] 1414 SGMT = np.array([ops[0].T for ops in SGData['SGOps']]) 1415 SGT = np.array([ops[1] for ops in SGData['SGOps']]) 1416 SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']]) 1417 SSGT = np.array([ops[1] for ops in SSGData['SSGOps']]) 1418 FFtables = calcControls['FFtables'] 1419 BLtables = calcControls['BLtables'] 1420 nRef = len(refDict['RefList']) 1421 Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict) 1422 mSize = len(Mdata) #no. atoms 1423 waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict) 1424 ngl,nWaves,Fmod,Xmod,Umod,glTau,glWt = G2mth.makeWaves(waveTypes,FSSdata,XSSdata,USSdata,Mast) 1425 waveShapes,SCtauF,SCtauX,SCtauU,UmodAB = G2mth.makeWavesDerv(ngl,waveTypes,FSSdata,XSSdata,USSdata,Mast) 1426 modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']]) 1427 FF = np.zeros(len(Tdata)) 1428 if 'NC' in calcControls[hfx+'histType']: 1429 FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam']) 1430 elif 'X' in calcControls[hfx+'histType']: 1431 FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata]) 1432 FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata]) 1433 Uij = np.array(G2lat.U6toUij(Uijdata)).T 1434 bij = Mast*Uij 1435 if not len(refDict['FF']): 1436 if 'N' in calcControls[hfx+'histType']: 1437 dat = G2el.getBLvalues(BLtables) #will need wave here for anom. neutron b's 1438 else: 1439 dat = G2el.getFFvalues(FFtables,0.) 1440 refDict['FF']['El'] = dat.keys() 1441 refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat))) 1442 dFdvDict = {} 1443 dFdfr = np.zeros((nRef,mSize)) 1444 dFdx = np.zeros((nRef,mSize,3)) 1445 dFdui = np.zeros((nRef,mSize)) 1446 dFdua = np.zeros((nRef,mSize,6)) 1447 dFdbab = np.zeros((nRef,2)) 1448 dFdfl = np.zeros((nRef)) 1449 dFdtw = np.zeros((nRef)) 1450 dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2)) 1451 dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6)) 1452 dFdGz = np.zeros((nRef,mSize,5)) 1453 dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12)) 1454 Flack = 1.0 1455 if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict: 1456 Flack = 1.-2.*parmDict[phfx+'Flack'] 1457 time0 = time.time() 1458 nRef = len(refDict['RefList'])/100 1459 for iref,refl in enumerate(refDict['RefList']): 1460 if 'T' in calcControls[hfx+'histType']: 1461 FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12+im]) 1462 H = np.array(refl[:4]) 1463 HP = H[:3]+modQ*H[3:] #projected hklm to hkl 1464 SQ = 1./(2.*refl[4+im])**2 # or (sin(theta)/lambda)**2 1465 SQfactor = 8.0*SQ*np.pi**2 1466 dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor) 1467 Bab = parmDict[phfx+'BabA']*dBabdA 1468 Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata]) 1469 FF = refDict['FF']['FF'][iref].T[Tindx] 1470 Uniq = np.inner(H,SSGMT) 1471 Phi = np.inner(H,SSGT) 1472 UniqP = np.inner(HP,SGMT) 1473 if SGInv: #if centro - expand HKL sets 1474 Uniq = np.vstack((Uniq,-Uniq)) 1475 Phi = np.hstack((Phi,-Phi)) 1476 UniqP = np.vstack((UniqP,-UniqP)) 1477 phase = twopi*(np.inner(Uniq[:,:3],(dXdata+Xdata).T)+Phi[:,nxs]) 1478 sinp = np.sin(phase) 1479 cosp = np.cos(phase) 1480 occ = Mdata*Fdata/Uniq.shape[0] 1481 biso = -SQfactor*Uisodata[:,nxs] 1482 Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),Uniq.shape[0],axis=1).T #ops x atoms 1483 HbH = -np.sum(UniqP[:,nxs,:3]*np.inner(UniqP[:,:3],bij),axis=-1) #ops x atoms 1484 Hij = np.array([Mast*np.multiply.outer(U[:3],U[:3]) for U in UniqP]) #atoms x 3x3 1485 Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij]) #atoms x 6 1486 Tuij = np.where(HbH<1.,np.exp(HbH),1.0) #ops x atoms 1487 Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/Uniq.shape[0] #ops x atoms 1488 fot = (FF+FP-Bab)*Tcorr #ops x atoms 1489 fotp = FPP*Tcorr #ops x atoms 1490 GfpuA = G2mth.Modulation(Uniq,UniqP,nWaves,Fmod,Xmod,Umod,glTau,glWt) #2 x sym X atoms 1491 dGdf,dGdx,dGdu,dGdz = G2mth.ModulationDerv(Uniq,UniqP,Hij,nWaves,waveShapes,Fmod,Xmod,UmodAB,SCtauF,SCtauX,SCtauU,glTau,glWt) 1492 # GfpuA is 2 x ops x atoms 1493 # derivs are: ops x atoms x waves x 2,6,12, or 5 parms as [real,imag] parts 1494 fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr]) # array(2,nEqv,nAtoms) 1495 fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr]) #or array(2,nEqv,nAtoms) 1496 fag = fa*GfpuA[0]-fb*GfpuA[1] 1497 fbg = fb*GfpuA[0]+fa*GfpuA[1] 1498 1499 fas = np.sum(np.sum(fag,axis=1),axis=1) # 2 x twin 1500 fbs = np.sum(np.sum(fbg,axis=1),axis=1) 1501 fax = np.array([-fot*sinp,-fotp*cosp]) #positions; 2 x ops x atoms 1502 fbx = np.array([fot*cosp,-fotp*sinp]) 1503 fax = fax*GfpuA[0]-fbx*GfpuA[1] 1504 fbx = fbx*GfpuA[0]+fax*GfpuA[1] 1505 #sum below is over Uniq 1506 dfadfr = np.sum(fag/occ,axis=1) #Fdata != 0 ever avoids /0. problem 1507 dfbdfr = np.sum(fbg/occ,axis=1) #Fdata != 0 avoids /0. problem 1508 dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1) 1509 dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1) 1510 dfadui = np.sum(-SQfactor*fag,axis=1) 1511 dfbdui = np.sum(-SQfactor*fbg,axis=1) 1512 dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2) #2 x nAtom x 3xyz; sum nOps 1513 dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2) 1514 dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2) #2 x nAtom x 6Uij; sum nOps 1515 dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2) #these are correct also for twins above 1516 # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps 1517 dfadGf = np.sum(fa[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1) 1518 dfbdGf = np.sum(fb[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1) 1519 dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1) 1520 dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1) 1521 dfadGz = np.sum(fa[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]-fb[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1) 1522 dfbdGz = np.sum(fb[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]+fa[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1) 1523 dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1) 1524 dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1) 1525 # GSASIIpath.IPyBreak() 1526 if not SGData['SGInv']: #Flack derivative 1527 dfadfl = np.sum(-FPP*Tcorr*sinp) 1528 dfbdfl = np.sum(FPP*Tcorr*cosp) 1529 else: 1530 dfadfl = 1.0 1531 dfbdfl = 1.0 1532 #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3! 1533 SA = fas[0]+fas[1] #float = A+A' (might be array[nTwin]) 1534 SB = fbs[0]+fbs[1] #float = B+B' (might be array[nTwin]) 1535 if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro? 1536 dFdfl[iref] = -SA*dfadfl-SB*dfbdfl #array(nRef,) 1537 dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+ \ 1538 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq) 1539 dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+ \ 1540 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1]) 1541 dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+ \ 1542 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1]) 1543 dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+ \ 1544 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1]) 1545 dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+ \ 1546 2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1]) 1547 dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+ \ 1548 2.*(fbs[0]*dfbdGx[0]-fbs[1]*dfbdGx[1]) 1549 dFdGz[iref] = 2.*(fas[0]*dfadGz[0]+fas[1]*dfadGz[1])+ \ 1550 2.*(fbs[0]*dfbdGz[0]+fbs[1]*dfbdGz[1]) 1551 dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+ \ 1552 2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1]) 1553 else: #OK, I think 1554 dFdfr[iref] = 2.*(SA*dfadfr[0]+SA*dfadfr[1]+SB*dfbdfr[0]+SB*dfbdfr[1])*Mdata/len(Uniq) #array(nRef,nAtom) 1555 dFdx[iref] = 2.*(SA*dfadx[0]+SA*dfadx[1]+SB*dfbdx[0]+SB*dfbdx[1]) #array(nRef,nAtom,3) 1556 dFdui[iref] = 2.*(SA*dfadui[0]+SA*dfadui[1]+SB*dfbdui[0]+SB*dfbdui[1]) #array(nRef,nAtom) 1557 dFdua[iref] = 2.*(SA*dfadua[0]+SA*dfadua[1]+SB*dfbdua[0]+SB*dfbdua[1]) #array(nRef,nAtom,6) 1558 dFdfl[iref] = -SA*dfadfl-SB*dfbdfl #array(nRef,) 1559 1560 dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1]) #array(nRef,natom,nwave,2) 1561 dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1]) #array(nRef,natom,nwave,6) 1562 dFdGz[iref] = 2.*(SA*dfadGz[0]+SB*dfbdGz[1]) #array(nRef,natom,5) 1563 dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1]) #array(nRef,natom,nwave,12) 1564 # GSASIIpath.IPyBreak() 1565 dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \ 1566 2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T 1567 #loop over atoms - each dict entry is list of derivatives for all the reflections 1568 if not iref%100 : 1569 print ' %d derivative time %.4f\r'%(iref,time.time()-time0), 1570 for i in range(len(Mdata)): #loop over atoms 1571 dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i] 1572 dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i] 1573 dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i] 1574 dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i] 1575 dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i] 1576 dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i] 1577 dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i] 1578 dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i] 1579 dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i] 1580 dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i] 1581 dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i] 1582 for j in range(FSSdata.shape[1]): #loop over waves Fzero & Fwid? 1583 dFdvDict[pfx+'Fsin:'+str(i)+':'+str(j)] = dFdGf.T[0][j][i] 1584 dFdvDict[pfx+'Fcos:'+str(i)+':'+str(j)] = dFdGf.T[1][j][i] 1585 nx = 0 1586 if waveTypes[i] in ['Block','ZigZag']: 1587 nx = 1 1588 dFdvDict[pfx+'Tmin:'+str(i)+':0'] = dFdGz.T[0][i] #ZigZag/Block waves (if any) 1589 dFdvDict[pfx+'Tmax:'+str(i)+':0'] = dFdGz.T[1][i] 1590 dFdvDict[pfx+'Xmax:'+str(i)+':0'] = dFdGz.T[2][i] 1591 dFdvDict[pfx+'Ymax:'+str(i)+':0'] = dFdGz.T[3][i] 1592 dFdvDict[pfx+'Zmax:'+str(i)+':0'] = dFdGz.T[4][i] 1593 for j in range(XSSdata.shape[1]-nx): #loop over waves 1594 dFdvDict[pfx+'Xsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[0][j][i] 1595 dFdvDict[pfx+'Ysin:'+str(i)+':'+str(j+nx)] = dFdGx.T[1][j][i] 1596 dFdvDict[pfx+'Zsin:'+str(i)+':'+str(j+nx)] = dFdGx.T[2][j][i] 1597 dFdvDict[pfx+'Xcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[3][j][i] 1598 dFdvDict[pfx+'Ycos:'+str(i)+':'+str(j+nx)] = dFdGx.T[4][j][i] 1599 dFdvDict[pfx+'Zcos:'+str(i)+':'+str(j+nx)] = dFdGx.T[5][j][i] 1600 for j in range(USSdata.shape[1]): #loop over waves 1601 dFdvDict[pfx+'U11sin:'+str(i)+':'+str(j)] = dFdGu.T[0][j][i] 1602 dFdvDict[pfx+'U22sin:'+str(i)+':'+str(j)] = dFdGu.T[1][j][i] 1603 dFdvDict[pfx+'U33sin:'+str(i)+':'+str(j)] = dFdGu.T[2][j][i] 1604 dFdvDict[pfx+'U12sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[3][j][i] 1605 dFdvDict[pfx+'U13sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[4][j][i] 1606 dFdvDict[pfx+'U23sin:'+str(i)+':'+str(j)] = 2.*dFdGu.T[5][j][i] 1607 dFdvDict[pfx+'U11cos:'+str(i)+':'+str(j)] = dFdGu.T[6][j][i] 1608 dFdvDict[pfx+'U22cos:'+str(i)+':'+str(j)] = dFdGu.T[7][j][i] 1609 dFdvDict[pfx+'U33cos:'+str(i)+':'+str(j)] = dFdGu.T[8][j][i] 1610 dFdvDict[pfx+'U12cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[9][j][i] 1611 dFdvDict[pfx+'U13cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[10][j][i] 1612 dFdvDict[pfx+'U23cos:'+str(i)+':'+str(j)] = 2.*dFdGu.T[11][j][i] 1613 1614 # GSASIIpath.IPyBreak() 1615 dFdvDict[phfx+'BabA'] = dFdbab.T[0] 1616 dFdvDict[phfx+'BabU'] = dFdbab.T[1] 1617 return dFdvDict 1618 1619 def SStructureFactorDervTw(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict): 1320 1620 'Needs a doc string' 1321 1621 phfx = pfx.split(':')[0]+hfx … … 1361 1661 refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat))) 1362 1662 dFdvDict = {} 1363 if nTwin > 1: 1364 dFdfr = np.zeros((nRef,nTwin,mSize)) 1365 dFdx = np.zeros((nRef,nTwin,mSize,3)) 1366 dFdui = np.zeros((nRef,nTwin,mSize)) 1367 dFdua = np.zeros((nRef,nTwin,mSize,6)) 1368 dFdbab = np.zeros((nRef,nTwin,2)) 1369 dFdfl = np.zeros((nRef,nTwin)) 1370 dFdtw = np.zeros((nRef,nTwin)) 1371 dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1])) 1372 dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3)) 1373 dFdGz = np.zeros((nRef,nTwin,mSize,5)) 1374 dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6)) 1375 else: 1376 dFdfr = np.zeros((nRef,mSize)) 1377 dFdx = np.zeros((nRef,mSize,3)) 1378 dFdui = np.zeros((nRef,mSize)) 1379 dFdua = np.zeros((nRef,mSize,6)) 1380 dFdbab = np.zeros((nRef,2)) 1381 dFdfl = np.zeros((nRef)) 1382 dFdtw = np.zeros((nRef)) 1383 dFdGf = np.zeros((nRef,mSize,FSSdata.shape[1],2)) 1384 dFdGx = np.zeros((nRef,mSize,XSSdata.shape[1],6)) 1385 dFdGz = np.zeros((nRef,mSize,5)) 1386 dFdGu = np.zeros((nRef,mSize,USSdata.shape[1],12)) 1663 dFdfr = np.zeros((nRef,nTwin,mSize)) 1664 dFdx = np.zeros((nRef,nTwin,mSize,3)) 1665 dFdui = np.zeros((nRef,nTwin,mSize)) 1666 dFdua = np.zeros((nRef,nTwin,mSize,6)) 1667 dFdbab = np.zeros((nRef,nTwin,2)) 1668 dFdfl = np.zeros((nRef,nTwin)) 1669 dFdtw = np.zeros((nRef,nTwin)) 1670 dFdGf = np.zeros((nRef,nTwin,mSize,FSSdata.shape[1])) 1671 dFdGx = np.zeros((nRef,nTwin,mSize,XSSdata.shape[1],3)) 1672 dFdGz = np.zeros((nRef,nTwin,mSize,5)) 1673 dFdGu = np.zeros((nRef,nTwin,mSize,USSdata.shape[1],6)) 1387 1674 Flack = 1.0 1388 1675 if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict: … … 1395 1682 H = np.array(refl[:4]) 1396 1683 HP = H[:3]+modQ*H[3:] #projected hklm to hkl 1397 # H = np.squeeze(np.inner(H.T,TwinLaw)) #maybe array(4,nTwins) or (4)1398 #TwMask = np.any(H,axis=-1)1399 #if TwinLaw.shape[0] > 1 and TwDict:1400 #if iref in TwDict:1401 #for i in TwDict[iref]:1402 #for n in range(NTL):1403 #H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])1404 #TwMask = np.any(H,axis=-1)1684 H = np.inner(H.T,TwinLaw) #maybe array(4,nTwins) or (4) 1685 TwMask = np.any(H,axis=-1) 1686 if TwinLaw.shape[0] > 1 and TwDict: 1687 if iref in TwDict: 1688 for i in TwDict[iref]: 1689 for n in range(NTL): 1690 H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM]) 1691 TwMask = np.any(H,axis=-1) 1405 1692 SQ = 1./(2.*refl[4+im])**2 # or (sin(theta)/lambda)**2 1406 1693 SQfactor = 8.0*SQ*np.pi**2 … … 1451 1738 dfadui = np.sum(-SQfactor*fag,axis=1) 1452 1739 dfbdui = np.sum(-SQfactor*fbg,axis=1) 1453 if nTwin > 1: 1454 dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1455 dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1456 dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1457 dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1458 # array(2,nTwin,nAtom,3) & array(2,nTwin,nAtom,6) & array(2,nTwin,nAtom,12) 1459 dfadGf = np.sum(fa[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1) 1460 dfbdGf = np.sum(fb[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1) 1461 dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1) 1462 dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1) 1463 dfadGz = np.sum(fa[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]-fb[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1) 1464 dfbdGz = np.sum(fb[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]+fa[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1) 1465 dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1) 1466 dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1) 1467 else: 1468 dfadx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fax,-2,-1)[:,:,:,nxs],axis=-2) #2 x nAtom x 3xyz; sum nOps 1469 dfbdx = np.sum(twopi*Uniq[:,:3]*np.swapaxes(fbx,-2,-1)[:,:,:,nxs],axis=-2) 1470 dfadua = np.sum(-Hij*np.swapaxes(fag,-2,-1)[:,:,:,nxs],axis=-2) #2 x nAtom x 6Uij; sum nOps 1471 dfbdua = np.sum(-Hij*np.swapaxes(fbg,-2,-1)[:,:,:,nxs],axis=-2) #these are correct also for twins above 1472 # array(2,nAtom,nWave,2) & array(2,nAtom,nWave,6) & array(2,nAtom,nWave,12); sum on nOps 1473 dfadGf = np.sum(fa[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1) 1474 dfbdGf = np.sum(fb[:,:,:,nxs,nxs]*dGdf[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdf[1][nxs,:,:,:,:],axis=1) 1475 dfadGx = np.sum(fa[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1) 1476 dfbdGx = np.sum(fb[:,:,:,nxs,nxs]*dGdx[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdx[1][nxs,:,:,:,:],axis=1) 1477 dfadGz = np.sum(fa[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]-fb[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1) 1478 dfbdGz = np.sum(fb[:,:,0,nxs,nxs]*dGdz[0][nxs,:,:,:]+fa[:,:,0,nxs,nxs]*dGdz[1][nxs,:,:,:],axis=1) 1479 dfadGu = np.sum(fa[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]-fb[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1) 1480 dfbdGu = np.sum(fb[:,:,:,nxs,nxs]*dGdu[0][nxs,:,:,:,:]+fa[:,:,:,nxs,nxs]*dGdu[1][nxs,:,:,:,:],axis=1) 1740 dfadx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1741 dfbdx = np.array([np.sum(twopi*Uniq[it,:,:3]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1742 dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fag,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1743 dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fbg,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)]) 1744 # array(2,nTwin,nAtom,3) & array(2,nTwin,nAtom,6) & array(2,nTwin,nAtom,12) 1745 dfadGf = np.sum(fa[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1) 1746 dfbdGf = np.sum(fb[:,it,:,:,nxs,nxs]*dGdf[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdf[1][nxs,nxs,:,:,:,:],axis=1) 1747 dfadGx = np.sum(fa[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1) 1748 dfbdGx = np.sum(fb[:,it,:,:,nxs,nxs]*dGdx[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdx[1][nxs,nxs,:,:,:,:],axis=1) 1749 dfadGz = np.sum(fa[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]-fb[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1) 1750 dfbdGz = np.sum(fb[:,it,:,0,nxs,nxs]*dGdz[0][nxs,nxs,:,:,:]+fa[:,it,:,0,nxs,nxs]*dGdz[1][nxs,nxs,:,:,:],axis=1) 1751 dfadGu = np.sum(fa[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]-fb[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1) 1752 dfbdGu = np.sum(fb[:,it,:,:,nxs,nxs]*dGdu[0][nxs,nxs,:,:,:,:]+fa[:,it,:,:,nxs,nxs]*dGdu[1][nxs,nxs,:,:,:,:],axis=1) 1481 1753 # GSASIIpath.IPyBreak() 1482 1754 if not SGData['SGInv'] and len(TwinLaw) == 1: #Flack derivative … … 1489 1761 SA = fas[0]+fas[1] #float = A+A' (might be array[nTwin]) 1490 1762 SB = fbs[0]+fbs[1] #float = B+B' (might be array[nTwin]) 1491 if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro 1492 dFdfl[iref] = -SA*dfadfl-SB*dfbdfl #array(nRef,) 1493 dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+ \ 1494 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq) 1495 dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+ \ 1496 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1]) 1497 dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+ \ 1498 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1]) 1499 dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+ \ 1500 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1]) 1501 dFdGf[iref] = 2.*(fas[0]*dfadGf[0]+fas[1]*dfadGf[1])+ \ 1502 2.*(fbs[0]*dfbdGf[0]+fbs[1]*dfbdGf[1]) 1503 dFdGx[iref] = 2.*(fas[0]*dfadGx[0]+fas[1]*dfadGx[1])+ \ 1504 2.*(fbs[0]*dfbdGx[0]-fbs[1]*dfbdGx[1]) 1505 dFdGz[iref] = 2.*(fas[0]*dfadGz[0]+fas[1]*dfadGz[1])+ \ 1506 2.*(fbs[0]*dfbdGz[0]+fbs[1]*dfbdGz[1]) 1507 dFdGu[iref] = 2.*(fas[0]*dfadGu[0]+fas[1]*dfadGu[1])+ \ 1508 2.*(fbs[0]*dfbdGu[0]+fbs[1]*dfbdGu[1]) 1509 else: 1510 if nTwin > 1: 1511 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)] 1512 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)] 1513 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)] 1514 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)] 1515 dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2 1763 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)] 1764 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)] 1765 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)] 1766 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)] 1767 dFdtw[iref] = np.sum(TwMask*fas,axis=0)**2+np.sum(TwMask*fbs,axis=0)**2 1516 1768 1517 dFdGf[iref] = [2.*TwMask[it]*(SA[it]*dfadGf[1]+SB[it]*dfbdGf[1]) for it in range(nTwin)] 1518 dFdGx[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)] 1519 dFdGz[iref] = [2.*TwMask[it]*(SA[it]*dfadGz[1]+SB[it]*dfbdGz[1]) for it in range(nTwin)] 1520 dFdGu[iref] = [2.*TwMask[it]*(SA[it]*dfadGu[1]+SB[it]*dfbdGu[1]) for it in range(nTwin)] 1521 else: #these are good for no twin single crystals 1522 dFdfr[iref] = 2.*(SA*dfadfr[0]+SA*dfadfr[1]+SB*dfbdfr[0]+SB*dfbdfr[1])*Mdata/len(Uniq) #array(nRef,nAtom) 1523 dFdx[iref] = 2.*(SA*dfadx[0]+SA*dfadx[1]+SB*dfbdx[0]+SB*dfbdx[1]) #array(nRef,nAtom,3) 1524 dFdui[iref] = 2.*(SA*dfadui[0]+SA*dfadui[1]+SB*dfbdui[0]+SB*dfbdui[1]) #array(nRef,nAtom) 1525 dFdua[iref] = 2.*(SA*dfadua[0]+SA*dfadua[1]+SB*dfbdua[0]+SB*dfbdua[1]) #array(nRef,nAtom,6) 1526 dFdfl[iref] = -SA*dfadfl-SB*dfbdfl #array(nRef,) 1527 1528 dFdGf[iref] = 2.*(SA*dfadGf[0]+SB*dfbdGf[1]) #array(nRef,natom,nwave,2) 1529 dFdGx[iref] = 2.*(SA*dfadGx[0]+SB*dfbdGx[1]) #array(nRef,natom,nwave,6) 1530 dFdGz[iref] = 2.*(SA*dfadGz[0]+SB*dfbdGz[1]) #array(nRef,natom,5) 1531 dFdGu[iref] = 2.*(SA*dfadGu[0]+SB*dfbdGu[1]) #array(nRef,natom,nwave,12) 1769 dFdGf[iref] = [2.*TwMask[it]*(SA[it]*dfadGf[1]+SB[it]*dfbdGf[1]) for it in range(nTwin)] 1770 dFdGx[iref] = [2.*TwMask[it]*(SA[it]*dfadGx[1]+SB[it]*dfbdGx[1]) for it in range(nTwin)] 1771 dFdGz[iref] = [2.*TwMask[it]*(SA[it]*dfadGz[1]+SB[it]*dfbdGz[1]) for it in range(nTwin)] 1772 dFdGu[iref] = [2.*TwMask[it]*(SA[it]*dfadGu[1]+SB[it]*dfbdGu[1]) for it in range(nTwin)] 1532 1773 # GSASIIpath.IPyBreak() 1533 1774 dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \ … … 2914 3155 if parmDict[phfx+'Scale'] < 0.: 2915 3156 parmDict[phfx+'Scale'] = .001 2916 if im: 2917 dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict) 3157 if im: # split to nontwin/twin versions 3158 if len(TwinLaw) > 1: 3159 dFdvDict = SStructureFactorDervTw(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict) 3160 else: 3161 dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict) 2918 3162 else: 2919 3163 if len(TwinLaw) > 1:
Note: See TracChangeset
for help on using the changeset viewer.