Changeset 1864


Ignore:
Timestamp:
May 22, 2015 11:08:59 AM (7 years ago)
Author:
vondreele
Message:

revise single crystal derivatives - now same as for old GSAS
Scale & extinction parms - good derivs for Uiso, F & F2 refinement
X,Y,Z 7 Uij approx. correct. Still some variation from numeric vals.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIstrMath.py

    r1862 r1864  
    938938        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
    939939       
    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
    942942        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
    943943        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
     
    948948        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
    949949        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)]).T
    956950        if not SGData['SGInv']:
    957951            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
     
    960954            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
    961955            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
    967982        #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)):
    969984        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
    970985        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
     
    24772492                w = 1.0/ref[6+im]
    24782493                if ref[3+im] > 0:
    2479                     wdf[iref] = 2.0*Fo*w*(Fo-Fc)
     2494                    wdf[iref] = 2.0*Fc*w*(Fo-Fc)
    24802495                    for j,var in enumerate(varylist):
    24812496                        if var in dFdvDict:
     
    24852500                            depDerivDict[var][iref] = w*dFdvDict[var][iref]*parmDict[phfx+'Scale']*ref[11+im]
    24862501                    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
    24882503                    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!
    24912506                        if phfx+item in varylist and phfx+item in dervDict:
    24922507                            dMdvh[varylist.index(phfx+item)][iref] = w*dervDict[phfx+item]/ref[11+im]
     
    27542769                        Fo = np.sqrt(ref[5+im])
    27552770                        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)?
    27572772                        if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]:    #skip sp.gp. absences (mul=0)
    27582773                            ref[3+im] = abs(ref[3+im])      #mark as allowed
Note: See TracChangeset for help on using the changeset viewer.