Ignore:
Timestamp:
Jun 11, 2015 4:20:44 PM (7 years ago)
Author:
vondreele
Message:

refactor StructureFactor2 & StructureFactorDerv? for twinning - deriv not complete

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIstrMath.py

    r1884 r1886  
    818818    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType']:
    819819        Flack = 1.-2.*parmDict[phfx+'Flack']
     820    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
     821    if 'S' in calcControls[hfx+'histType']:
     822        TwinLaw = calcControls[phfx+'TwinLaw']
     823        TwinFr = np.array([parmDict[phfx+'TwinFr;'+str(i)] for i in range(len(TwinLaw))])
     824        if len(TwinLaw) > 1:
     825            TwinFr[0] = 1.-np.sum(TwinFr[1:])       
    820826    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
    821827    FF = np.zeros(len(Tdata))
     
    846852    while iBeg < nRef:
    847853        iFin = min(iBeg+blkSize,nRef)
    848         refl = refDict['RefList'][iBeg:iFin]
    849         H = refl.T[:3]
    850         SQ = 1./(2.*refl.T[4])**2
    851         SQfactor = 4.0*SQ*twopisq
     854        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
     855        H = refl.T[:3]                          #array(blkSize,3)
     856        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,3,nTwins) or (blkSize,3)
     857        SQ = 1./(2.*refl.T[4])**2               #array(blkSize)
     858        SQfactor = 4.0*SQ*twopisq               #ditto prev.
    852859        if 'T' in calcControls[hfx+'histType']:
    853860            if 'P' in calcControls[hfx+'histType']:
     
    855862            else:
    856863                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
    857             FP = np.repeat(FP.T,len(SGT),axis=0)
    858             FPP = np.repeat(FPP.T,len(SGT),axis=0)
    859         Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT))
     864            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
     865            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
     866        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
    860867        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    861         FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT),axis=0)
    862         Uniq = np.reshape(np.inner(H.T,SGMT),(-1,3))
    863         Phi = np.inner(H.T,SGT).flatten()
    864         phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T)+Phi[:,np.newaxis])
     868        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
     869        Uniq = np.inner(H,SGMT)
     870        Phi = np.inner(H,SGT)
     871        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
    865872        sinp = np.sin(phase)
    866873        cosp = np.cos(phase)
    867874        biso = -SQfactor*Uisodata[:,np.newaxis]
    868         Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
    869         HbH = -np.sum(Uniq.T*np.inner(bij,Uniq),axis=1)
     875        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
     876        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
    870877        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
    871         Tcorr = Tiso*Tuij*Mdata*Fdata/len(SGMT)
    872         fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    873         fa = np.reshape(fa,(2,len(refl),len(SGT),len(Mdata)))   #real A,-b
    874         fas = np.sum(np.sum(fa,axis=2),axis=2)        #real sum over atoms & unique hkl
    875         fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
    876         fb = np.reshape(fb,(2,len(refl),len(SGT),len(Mdata)))   #imag -B,+a       
    877         fbs = np.sum(np.sum(fb,axis=2),axis=2)  #imag sum over atoms & uniq hkl
     878        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
     879        if 'T' in calcControls[hfx+'histType']:
     880            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
     881            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
     882        else:
     883            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
     884            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
     885        fas = np.sum(np.sum(fa,axis=-1),axis=-1)       #real sum over atoms & unique hkl
     886        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag sum over atoms & uniq hkl
    878887        if SGData['SGInv']: #centrosymmetric; B=0
    879888            fbs[0] *= 0.
    880889        if 'P' in calcControls[hfx+'histType']:
    881890            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0)
     891            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    882892        else:
    883             refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
    884         refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"?
     893            if len(TwinLaw) > 1:
     894                refl.T[9] = np.sum(fas[:,:,0]**2,axis=0)+np.sum(fbs[:,:,0]**2,axis=0)   #FcT from primary twin element
     895                refl.T[7] = np.sum(TwinFr*np.sum(fas,axis=0)**2,axis=-1)+   \
     896                    np.sum(TwinFr*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
     897#                what goes in refl.T[8]? (FoT)           
     898                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f"
     899            else:
     900                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
     901                refl.T[7] = np.copy(refl.T[9])               
     902                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    885903#        refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
    886904        iBeg += blkSize
     
    895913    FFtables = calcControls['FFtables']
    896914    BLtables = calcControls['BLtables']
     915    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
     916    if 'S' in calcControls[hfx+'histType']:
     917        TwinLaw = calcControls[phfx+'TwinLaw']
     918        TwinFr = np.array([parmDict[phfx+'TwinFr;'+str(i)] for i in range(len(TwinLaw))])
     919        if len(TwinLaw) > 1:
     920            TwinFr[0] = 1.-np.sum(TwinFr[1:])       
    897921    nRef = len(refDict['RefList'])
    898922    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     
    920944            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
    921945        H = np.array(refl[:3])
     946        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,3,nTwins) or (blkSize,3)
    922947        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
    923948        SQfactor = 8.0*SQ*np.pi**2
     
    925950        Bab = parmDict[phfx+'BabA']*dBabdA
    926951        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    927         FF = refDict['FF']['FF'][iref].T[Tindx]
     952        FF = refDict['FF']['FF'][iref].T[Tindx].T
    928953        Uniq = np.inner(H,SGMT)
    929954        Phi = np.inner(H,SGT)
    930         phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+Phi[np.newaxis,:])
     955        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
    931956        sinp = np.sin(phase)
    932957        cosp = np.cos(phase)
    933958        occ = Mdata*Fdata/len(Uniq)
    934         biso = -SQfactor*Uisodata
    935         Tiso = np.where(biso<1.,np.exp(biso),1.0)
    936         HbH = -np.inner(H,np.inner(bij,H))
    937         Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
    938         Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
    939         Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
    940         Tcorr = Tiso*Tuij
     959        biso = -SQfactor*Uisodata[:,np.newaxis]
     960        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
     961        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
     962        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)))
     964        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
     965        Tcorr = np.reshape(Tiso,Tuij.shape)*Tuij*Mdata*Fdata/len(SGMT)
    941966        fot = (FF+FP-Bab)*occ*Tcorr
    942         fotp = FPP*occ*Tcorr
    943         fa = np.array([fot[:,np.newaxis]*cosp,-Flack*fotp[:,np.newaxis]*sinp])       #non positions
    944         fb = np.array([fot[:,np.newaxis]*sinp,Flack*fotp[:,np.newaxis]*cosp])
    945        
    946         fas = np.sum(np.sum(fa,axis=1),axis=1)      #real sum over atoms & unique hkl
    947         fbs = np.sum(np.sum(fb,axis=1),axis=1)      #imag sum over atoms & uniq hkl
    948         fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*cosp])   #positions
    949         fbx = np.array([fot[:,np.newaxis]*cosp,-fotp[:,np.newaxis]*sinp])
    950         #sum below is over Uniq
    951         dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
    952         dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
    953         dfadui = np.sum(-SQfactor*fa,axis=2)
    954         dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
     967        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])
     974        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl
     975        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
     976        fax = np.array([-fot*sinp,-fotp*cosp])   #positions
     977        fbx = np.array([fot*cosp,-fotp*sinp])
     978        #sum below is over Uniq - twin effects??
     979        dfadfr = np.sum(fa/occ,axis=-2)        #Fdata != 0 ever avoids /0. problem
     980        dfadx = np.sum(twopi*Uniq*np.swapaxes(fax,-2,-1)[:,:,:,np.newaxis],axis=-2)
     981        dfadui = np.sum(-SQfactor*fa,axis=-2)
     982        dfadua = np.sum(-Hij*np.swapaxes(fa,-2,-1)[:,:,:,np.newaxis],axis=-2)
    955983        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
    956984        if not SGData['SGInv']:
    957             dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #Fdata != 0 ever avoids /0. problem
    958             dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)           
    959             dfbdui = np.sum(-SQfactor*fb,axis=2)
    960             dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
     985            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)
    961989            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
    962990            dfadfl = np.sum(-fotp[:,np.newaxis]*sinp)
     
    27472775                        ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
    27482776                        w = 1.0/ref[6+im]   # 1/sig(F^2)
    2749                         ref[7+im] = parmDict[phfx+'Scale']*ref[9+im]*ref[11+im]  #correct Fc^2 for extinction
     2777                        ref[7+im] *= parmDict[phfx+'Scale']*ref[11+im]  #correct Fc^2 for extinction
    27502778                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
    27512779                        if UserRejectHKL(ref,im,calcControls['UsrReject']) and ref[3+im]:    #skip sp.gp. absences (mul=0)
     
    27692797                    if ref[5+im] > 0.:
    27702798                        ref[11+im] = SCExtinction(ref,im,phfx,hfx,pfx,calcControls,parmDict,varylist)[0]
    2771                         ref[7+im] = parmDict[phfx+'Scale']*ref[9+im]*ref[11+im]    #correct Fc^2 for extinction
     2799                        ref[7+im] *= parmDict[phfx+'Scale']*ref[11+im]    #correct Fc^2 for extinction
    27722800                        ref[8+im] = ref[5+im]/(parmDict[phfx+'Scale']*ref[11+im])
    27732801                        Fo = np.sqrt(ref[5+im])
Note: See TracChangeset for help on using the changeset viewer.