Ignore:
Timestamp:
Mar 23, 2018 1:41:13 PM (4 years ago)
Author:
vondreele
Message:

implement new delt/sig subplot for PWDR & SASD plots - nvoked with 'w' key
make new MagStructureFactor2 routine - removed magnetism stuff from StructureFactor2

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIstrMath.py

    r3297 r3320  
    666666    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
    667667    SGT = np.array([ops[1] for ops in SGData['SGOps']])
    668     Ncen = len(SGData['SGCen'])
    669     Nops = len(SGMT)*Ncen*(1+SGData['SGInv'])
    670668    FFtables = calcControls['FFtables']
    671669    BLtables = calcControls['BLtables']
    672     MFtables = calcControls['MFtables']
    673670    Amat,Bmat = G2lat.Gmat2AB(G)
    674671    Flack = 1.0
     
    687684    if not Xdata.size:          #no atoms in phase!
    688685        return
    689     if parmDict[pfx+'isMag']:       #TODO: fix the math - mag moments now along crystal axes
    690         Mag = np.sqrt(np.sum(Gdata**2,axis=0))      #magnitude of moments for uniq atoms
    691         Gdata = np.where(Mag>0.,Gdata/Mag,0.)       #normalze mag. moments
    692         Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
    693         if SGData['SGInv'] and not SGData['SGFixed']:
    694             Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
    695         Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
    696         Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip * det(opM)
    697         Mag = np.tile(Mag[:,nxs],len(SGMT)*Ncen).T
    698         if SGData['SGInv'] and not SGData['SGFixed']:
    699             Mag = np.repeat(Mag,2,axis=0)                  #Mag same shape as Gdata
    700686    if 'NC' in calcControls[hfx+'histType']:
    701687        FP,FPP = G2el.BlenResCW(Tdata,BLtables,parmDict[hfx+'Lam'])
     
    712698        refDict['FF']['El'] = list(dat.keys())
    713699        refDict['FF']['FF'] = np.ones((nRef,len(dat)))*list(dat.values())
    714         refDict['FF']['MF'] = np.zeros((nRef,len(dat)))
    715         for iel,El in enumerate(refDict['FF']['El']):
    716             if El in MFtables:
    717                 refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
    718700    else:       #'X'
    719701        dat = G2el.getFFvalues(FFtables,0.)
     
    759741        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    760742        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
    761         if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:       #TODO: math here??
    762             MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
    763             TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
    764             if SGData['SGInv'] and not SGData['SGFixed']:
    765                 mphase = np.hstack((phase,-phase))
    766             else:
    767                 mphase = phase
    768             mphase = np.array([mphase+twopi*np.inner(cen,H)[:,nxs,nxs] for cen in SGData['SGCen']])
    769             mphase = np.concatenate(mphase,axis=1)              #Nref,full Nop,Natm
    770             sinm = np.sin(mphase)                               #ditto - match magstrfc.for
    771             cosm = np.cos(mphase)                               #ditto
    772             HM = np.inner(Bmat.T,H)                             #put into cartesian space
    773             HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
    774             eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0)
    775             Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Gdata[:,nxs,:,:] #xyz,Nref,Nop,Natm = BPM in magstrfc.for OK
    776             fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
    777             fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
    778             fams = np.sum(np.sum(fam,axis=-1),axis=-1)                          #xyz,Nref
    779             fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)                          #ditto
    780             refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
    781             refl.T[7] = np.copy(refl.T[9])               
    782             refl.T[10] = 0.0    #atan2d(fbs[0],fas[0]) - what is phase for mag refl?
     743        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
     744        if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
     745            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
     746            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
    783747        else:
    784             Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
    785             if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
    786                 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(Flack*FPP,sinp.shape)*sinp*Tcorr])
    787                 fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(Flack*FPP,cosp.shape)*cosp*Tcorr])
    788             else:
    789                 fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
    790                 fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
    791             fas = np.sum(np.sum(fa,axis=-1),axis=-1)  #real 2 x blkSize x nTwin; sum over atoms & uniq hkl
    792             fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag
    793             if SGData['SGInv']: #centrosymmetric; B=0
    794                 fbs[0] *= 0.
    795                 fas[1] *= 0.
    796             if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
    797                 refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
     748            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-Flack*FPP*sinp*Tcorr])
     749            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,Flack*FPP*cosp*Tcorr])
     750        fas = np.sum(np.sum(fa,axis=-1),axis=-1)  #real 2 x blkSize x nTwin; sum over atoms & uniq hkl
     751        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)  #imag
     752        if SGData['SGInv']: #centrosymmetric; B=0
     753            fbs[0] *= 0.
     754            fas[1] *= 0.
     755        if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
     756            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
     757            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
     758        else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
     759            if len(TwinLaw) > 1:
     760                refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
     761                refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
     762                    np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
     763                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
     764            else:   # checked correct!!
     765                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
     766                refl.T[7] = np.copy(refl.T[9])               
    798767                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    799                 if 'N' in calcControls[hfx+'histType'] and parmDict[pfx+'isMag']:
    800                     refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
    801             else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
    802                 if len(TwinLaw) > 1:
    803                     refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
    804                     refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
    805                         np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
    806                     refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
    807                 else:   # checked correct!!
    808                     refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
    809                     refl.T[7] = np.copy(refl.T[9])               
    810                     refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    811768#                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
    812769        iBeg += blkSize
    813770#    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
    814771   
     772def MagStructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
     773    ''' Compute neutron magnetic structure factors for all h,k,l for phase
     774    puts the result, F^2, in each ref[8] in refList
     775    operates on blocks of 100 reflections for speed
     776    input:
     777   
     778    :param dict refDict: where
     779        'RefList' list where each ref = h,k,l,it,d,...
     780        'FF' dict of form factors - filed in below
     781    :param np.array G:      reciprocal metric tensor
     782    :param str pfx:    phase id string
     783    :param dict SGData: space group info. dictionary output from SpcGroup
     784    :param dict calcControls:
     785    :param dict ParmDict:
     786
     787    '''       
     788#    phfx = pfx.split(':')[0]+hfx
     789    ast = np.sqrt(np.diag(G))
     790    Mast = twopisq*np.multiply.outer(ast,ast)
     791    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     792    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     793    Ncen = len(SGData['SGCen'])
     794    Nops = len(SGMT)*Ncen*(1+SGData['SGInv'])
     795    MFtables = calcControls['MFtables']
     796    Amat,Bmat = G2lat.Gmat2AB(G)
     797    TwinLaw = np.ones(1)
     798#    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
     799#    TwDict = refDict.get('TwDict',{})           
     800#    if 'S' in calcControls[hfx+'histType']:
     801#        NTL = calcControls[phfx+'NTL']
     802#        NM = calcControls[phfx+'TwinNMN']+1
     803#        TwinLaw = calcControls[phfx+'TwinLaw']
     804#        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
     805#        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
     806    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata,Gdata = \
     807        GetAtomFXU(pfx,calcControls,parmDict)
     808    if not Xdata.size:          #no atoms in phase!
     809        return
     810    Mag = np.sqrt(np.sum(Gdata**2,axis=0))      #magnitude of moments for uniq atoms
     811    Gdata = np.where(Mag>0.,Gdata/Mag,0.)       #normalze mag. moments
     812    Gdata = np.inner(Bmat,Gdata.T)              #convert to crystal space
     813    Gdata = np.inner(Gdata.T,SGMT).T            #apply sym. ops.
     814    if SGData['SGInv'] and not SGData['SGFixed']:
     815        Gdata = np.hstack((Gdata,-Gdata))       #inversion if any
     816    Gdata = np.hstack([Gdata for icen in range(Ncen)])        #dup over cell centering
     817    Gdata = SGData['MagMom'][nxs,:,nxs]*Gdata   #flip vectors according to spin flip * det(opM)
     818    Gdata = np.inner(Amat,Gdata.T)              #convert back to cart. space MXYZ, Natoms, NOps*Inv*Ncen
     819    Gdata = np.swapaxes(Gdata,1,2)              # put Natoms last
     820    Mag = np.tile(Mag[:,nxs],len(SGMT)*Ncen).T
     821    if SGData['SGInv'] and not SGData['SGFixed']:
     822        Mag = np.repeat(Mag,2,axis=0)                  #Mag same shape as Gdata
     823    Uij = np.array(G2lat.U6toUij(Uijdata))
     824    bij = Mast*Uij.T
     825    blkSize = 100       #no. of reflections in a block - size seems optimal
     826    nRef = refDict['RefList'].shape[0]
     827    SQ = 1./(2.*refDict['RefList'].T[4])**2
     828    refDict['FF']['El'] = list(MFtables.keys())
     829    refDict['FF']['MF'] = np.zeros((nRef,len(MFtables)))
     830    for iel,El in enumerate(refDict['FF']['El']):
     831        refDict['FF']['MF'].T[iel] = G2el.MagScatFac(MFtables[El],SQ)
     832#reflection processing begins here - big arrays!
     833    iBeg = 0
     834    while iBeg < nRef:
     835        iFin = min(iBeg+blkSize,nRef)
     836        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
     837        H = refl.T[:3]                          #array(blkSize,3)
     838#        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,3) or (blkSize,3)
     839#        TwMask = np.any(H,axis=-1)
     840#        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
     841#            for ir in range(blkSize):
     842#                iref = ir+iBeg
     843#                if iref in TwDict:
     844#                    for i in TwDict[iref]:
     845#                        for n in range(NTL):
     846#                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
     847#            TwMask = np.any(H,axis=-1)
     848        SQ = 1./(2.*refl.T[4])**2               #array(blkSize)
     849        SQfactor = 4.0*SQ*twopisq               #ditto prev.
     850        Uniq = np.inner(H.T,SGMT)
     851        Phi = np.inner(H.T,SGT)
     852        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
     853        biso = -SQfactor*Uisodata[:,nxs]
     854        Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT)*len(TwinLaw),axis=1).T
     855        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
     856        Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
     857        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     858        MF = refDict['FF']['MF'][iBeg:iFin].T[Tindx].T   #Nref,Natm
     859        TMcorr = 0.539*(np.reshape(Tiso,Tuij.shape)*Tuij)[:,0,:]*Fdata*Mdata*MF/(2*Nops)     #Nref,Natm
     860        if SGData['SGInv'] and not SGData['SGFixed']:
     861            mphase = np.hstack((phase,-phase))
     862        else:
     863            mphase = phase
     864        mphase = np.array([mphase+twopi*np.inner(cen,H.T)[:,nxs,nxs] for cen in SGData['SGCen']])
     865        mphase = np.concatenate(mphase,axis=1)              #Nref,full Nop,Natm
     866        sinm = np.sin(mphase)                               #ditto - match magstrfc.for
     867        cosm = np.cos(mphase)                               #ditto
     868        HM = np.inner(Bmat.T,H.T)                             #put into cartesian space
     869        HM = HM/np.sqrt(np.sum(HM**2,axis=0))               #Gdata = MAGS & HM = UVEC in magstrfc.for both OK
     870        eDotK = np.sum(HM[:,:,nxs,nxs]*Gdata[:,nxs,:,:],axis=0)
     871        Q = HM[:,:,nxs,nxs]*eDotK[nxs,:,:,:]-Gdata[:,nxs,:,:] #xyz,Nref,Nop,Natm = BPM in magstrfc.for OK
     872        fam = Q*TMcorr[nxs,:,nxs,:]*cosm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
     873        fbm = Q*TMcorr[nxs,:,nxs,:]*sinm[nxs,:,:,:]*Mag[nxs,nxs,:,:]    #ditto
     874        fams = np.sum(np.sum(fam,axis=-1),axis=-1)                          #xyz,Nref
     875        fbms = np.sum(np.sum(fbm,axis=-1),axis=-1)                          #ditto
     876        refl.T[9] = np.sum(fams**2,axis=0)+np.sum(fbms**2,axis=0)
     877        refl.T[7] = np.copy(refl.T[9])               
     878        refl.T[10] = 0.0    #atan2d(fbs[0],fas[0]) - what is phase for mag refl?
     879#        if 'P' in calcControls[hfx+'histType']:     #PXC, PNC & PNT: F^2 = A[0]^2 + A[1]^2 + B[0]^2 + B[1]^2
     880#            refl.T[9] = np.sum(fas**2,axis=0)+np.sum(fbs**2,axis=0) #add fam**2 & fbm**2 here   
     881#            refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
     882#        else:                                       #HKLF: F^2 = (A[0]+A[1])^2 + (B[0]+B[1])^2
     883#            if len(TwinLaw) > 1:
     884#                refl.T[9] = np.sum(fas[:,:,0],axis=0)**2+np.sum(fbs[:,:,0],axis=0)**2   #FcT from primary twin element
     885#                refl.T[7] = np.sum(TwinFr*TwMask*np.sum(fas,axis=0)**2,axis=-1)+   \
     886#                    np.sum(TwinFr*TwMask*np.sum(fbs,axis=0)**2,axis=-1)                        #Fc sum over twins
     887#                refl.T[10] = atan2d(fbs[0].T[0],fas[0].T[0])  #ignore f' & f" & use primary twin
     888#            else:   # checked correct!!
     889#                refl.T[9] = np.sum(fas,axis=0)**2+np.sum(fbs,axis=0)**2
     890#                refl.T[7] = np.copy(refl.T[9])               
     891#                refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
     892##                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
     893        iBeg += blkSize
     894#    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
     895
    815896def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    816897    '''Compute structure factor derivatives on blocks of reflections - for powders/nontwins only
     
    31833264            if im:
    31843265                SStructureFactor(refDict,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
     3266            elif parmDict[pfx+'isMag'] and 'N' in calcControls[hfx+'histType']:
     3267                MagStructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    31853268            else:
    31863269                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
Note: See TracChangeset for help on using the changeset viewer.