Changeset 1864
- Timestamp:
- May 22, 2015 11:08:59 AM (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIstrMath.py
r1862 r1864 938 938 fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp]) 939 939 940 fas = np.sum(np.sum(fa,axis=1),axis=1) 941 fbs = np.sum(np.sum(fb,axis=1),axis=1) 940 fas = np.sum(np.sum(fa,axis=1),axis=1) #real sum over atoms & unique hkl 941 fbs = np.sum(np.sum(fb,axis=1),axis=1) #imag sum over atoms & uniq hkl 942 942 fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp]) #positions 943 943 fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp]) … … 948 948 dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2) 949 949 dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1) 950 #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for al2O3!951 dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)952 dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])953 dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])954 dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])955 dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T956 950 if not SGData['SGInv']: 957 951 dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2) #Fdata != 0 ever avoids /0. problem … … 960 954 dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2) 961 955 dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1) 962 dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq) 963 dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1]) 964 dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1]) 965 dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1]) 966 dFdbab[iref] += 2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T 956 else: 957 dfbdfr = np.zeros_like(dfadfr) 958 dfbdx = np.zeros_like(dfadx) 959 dfbdui = np.zeros_like(dfadui) 960 dfbdua = np.zeros_like(dfadua) 961 dfbdba = np.zeros_like(dfadba) 962 #NB: the above have been checked against PA(1:10,1:2) in strfctr.for for Al2O3! 963 if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro 964 dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)+ \ 965 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq) 966 dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])+ \ 967 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1]) 968 dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])+ \ 969 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1]) 970 dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])+ \ 971 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1]) 972 else: 973 SA = fas[0]-fbs[1] 974 SB = fbs[0]+fas[1] 975 dFdfr[iref] = 2.*SA*(dfadfr[0]+dfbdfr[1])*Mdata/len(Uniq)+ \ 976 2.*SB*(dfbdfr[0]+dfadfr[1])*Mdata/len(Uniq) 977 dFdx[iref] = 2.*SA*(dfadx[0]+dfbdx[1])+2.*SB*(dfbdx[0]+dfadx[1]) 978 dFdui[iref] = 2.*SA*(dfadui[0]+dfbdui[1])+2.*SB*(dfbdui[0]+dfadui[1]) 979 dFdua[iref] = 2.*SA*(dfadua[0]+dfbdua[1])+2.*SB*(dfbdua[0]+dfadua[1]) 980 dFdbab[iref] = 2.*fas[0]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \ 981 2.*fbs[0]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T 967 982 #loop over atoms - each dict entry is list of derivatives for all the reflections 968 for i in range(len(Mdata)): 983 for i in range(len(Mdata)): 969 984 dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i] 970 985 dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i] … … 2477 2492 w = 1.0/ref[6+im] 2478 2493 if ref[3+im] > 0: 2479 wdf[iref] = 2.0*F o*w*(Fo-Fc)2494 wdf[iref] = 2.0*Fc*w*(Fo-Fc) 2480 2495 for j,var in enumerate(varylist): 2481 2496 if var in dFdvDict: … … 2485 2500 depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im] 2486 2501 if phfx+'Scale' in varylist: 2487 dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9+im]*ref[11+im] 2502 dMdvh[varylist.index(phfx+'Scale')][iref] = w*ref[9+im]*ref[11+im] #OK 2488 2503 elif phfx+'Scale' in dependentVars: 2489 depDerivDict[phfx+'Scale'][iref] = w*ref[9+im]*ref[11+im] 2490 for item in ['Ep','Es','Eg']: 2504 depDerivDict[phfx+'Scale'][iref] = w*ref[9+im]*ref[11+im] #OK 2505 for item in ['Ep','Es','Eg']: #OK! 2491 2506 if phfx+item in varylist and phfx+item in dervDict: 2492 2507 dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im] … … 2754 2769 Fo = np.sqrt(ref[5+im]) 2755 2770 Fc = np.sqrt(ref[7+im]) 2756 w = 2.0* (Fo/ref[6+im])**2# 1/sig(F)?2771 w = 2.0*Fo/ref[6+im] # 1/sig(F)? 2757 2772 if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]: #skip sp.gp. absences (mul=0) 2758 2773 ref[3+im] = abs(ref[3+im]) #mark as allowed
Note: See TracChangeset
for help on using the changeset viewer.