Changeset 1979


Ignore:
Timestamp:
Sep 28, 2015 8:59:29 AM (6 years ago)
Author:
vondreele
Message:

fix to SS structure factor now works
start SSderiv & block hkl processing

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIplot.py

    r1978 r1979  
    320320    '''
    321321    from matplotlib.patches import Circle,CirclePolygon
     322    global HKL,HKLF,HKLref
    322323    HKLref = hklRef
    323     global HKL,HKLF,HKLref
    324324   
    325325    def OnSCKeyPress(event):
  • trunk/GSASIIstrMath.py

    r1978 r1979  
    971971    phfx = pfx.split(':')[0]+hfx
    972972    ast = np.sqrt(np.diag(G))
    973     Mast = twopisq*np.multiply.outer(ast,ast)
    974    
    975     SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     973    Mast = twopisq*np.multiply.outer(ast,ast)   
     974#    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
    976975    SGInv = SGData['SGInv']
    977     SGT = np.array([ops[1] for ops in SGData['SGOps']])
     976#    SGT = np.array([ops[1] for ops in SGData['SGOps']])
    978977    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
    979978    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
    980979    FFtables = calcControls['FFtables']
    981980    BLtables = calcControls['BLtables']
     981    Flack = 1.0
     982    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
     983        Flack = 1.-2.*parmDict[phfx+'Flack']
     984    TwinLaw = np.array([[[1,0,0],[0,1,0],[0,0,1]],])
     985    TwDict = refDict.get('TwDict',{})           
     986    if 'S' in calcControls[hfx+'histType']:
     987        NTL = calcControls[phfx+'NTL']
     988        NM = calcControls[phfx+'TwinNMN']+1
     989        TwinLaw = calcControls[phfx+'TwinLaw']
     990        TwinFr = np.array([parmDict[phfx+'TwinFr:'+str(i)] for i in range(len(TwinLaw))])
     991        TwinInv = list(np.where(calcControls[phfx+'TwinInv'],-1,1))
    982992    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
    983993    waveTypes,FSSdata,XSSdata,USSdata,MSSdata = GetAtomSSFXU(pfx,calcControls,parmDict)
    984     SStauM = list(GetSSTauM(SGData['SGOps'],SSGData['SSGOps'],pfx,calcControls,Xdata))
    985     if SGData['SGInv']:
    986         SStauM[0] = np.hstack((SStauM[0],SStauM[0]))
    987         SStauM[1] = np.hstack((SStauM[1],SStauM[1]))
    988994    modQ = np.array([parmDict[pfx+'mV0'],parmDict[pfx+'mV1'],parmDict[pfx+'mV2']])
    989995    FF = np.zeros(len(Tdata))
     
    9951001    Uij = np.array(G2lat.U6toUij(Uijdata))
    9961002    bij = Mast*Uij.T
     1003    blkSize = 100       #no. of reflections in a block
     1004    nRef = refDict['RefList'].shape[0]
    9971005    if not len(refDict['FF']):
    9981006        if 'N' in calcControls[hfx+'histType']:
    9991007            dat = G2el.getBLvalues(BLtables)        #will need wave here for anom. neutron b's
     1008            refDict['FF']['El'] = dat.keys()
     1009            refDict['FF']['FF'] = np.ones((nRef,len(dat)))*dat.values()           
    10001010        else:
    10011011            dat = G2el.getFFvalues(FFtables,0.)       
    1002         refDict['FF']['El'] = dat.keys()
    1003         refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))
     1012            refDict['FF']['El'] = dat.keys()
     1013            refDict['FF']['FF'] = np.ones((nRef,len(dat)))
     1014            for iref,ref in enumerate(refDict['RefList']):
     1015                SQ = 1./(2.*ref[4])**2
     1016                dat = G2el.getFFvalues(FFtables,SQ)
     1017                refDict['FF']['FF'][iref] *= dat.values()
    10041018    time0 = time.time()
    10051019    nref = len(refDict['RefList'])/100   
    1006     for iref,refl in enumerate(refDict['RefList']):
    1007         if 'NT' in calcControls[hfx+'histType']:
    1008             FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl[14+im])
    1009         fbs = np.array([0,0])
    1010         H = refl[:4]
    1011         SQ = 1./(2.*refl[4+im])**2
    1012         SQfactor = 4.0*SQ*twopisq
    1013         Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
    1014         if not np.any(refDict['FF']['FF'][iref]):                #no form factors - 1st time thru StructureFactor
    1015             if 'N' in calcControls[hfx+'histType']:
    1016                 dat = G2el.getBLvalues(BLtables)
    1017                 refDict['FF']['FF'][iref] = dat.values()
    1018             else:       #'X'
    1019                 dat = G2el.getFFvalues(FFtables,SQ)
    1020                 refDict['FF']['FF'][iref] = dat.values()
     1020#reflection processing begins here - big arrays!
     1021    iBeg = 0
     1022    while iBeg < nRef:
     1023        iFin = min(iBeg+blkSize,nRef)
     1024        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
     1025        H = refl.T[:4]                          #array(blkSize,4)
     1026        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(blkSize,nTwins,3) or (blkSize,3)
     1027        TwMask = np.any(H,axis=-1)
     1028        if TwinLaw.shape[0] > 1 and TwDict: #need np.inner(TwinLaw[?],TwDict[iref][i])*TwinInv[i]
     1029            for ir in range(blkSize):
     1030                iref = ir+iBeg
     1031                if iref in TwDict:
     1032                    for i in TwDict[iref]:
     1033                        for n in range(NTL):
     1034                            H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
     1035            TwMask = np.any(H,axis=-1)
     1036        SQ = 1./(2.*refl.T[5])**2               #array(blkSize)
     1037        SQfactor = 4.0*SQ*twopisq               #ditto prev.
     1038        if 'T' in calcControls[hfx+'histType']:
     1039            if 'P' in calcControls[hfx+'histType']:
     1040                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
     1041            else:
     1042                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
     1043            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
     1044            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
     1045        Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT)*len(TwinLaw))
    10211046        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    1022         FF = refDict['FF']['FF'][iref][Tindx]
    1023         HM = np.zeros(4)
    1024         HM[:3] += H[3]*modQ
    1025         HM += H
    1026         Uniq = np.inner(HM[:3],SGMT)
    1027         SSUniq = np.inner(H,SSGMT)
    1028         Phi = np.inner(HM[:3],SGT)
    1029         SSPhi = np.inner(H,SSGT)
     1047        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
     1048        Uniq = np.inner(H,SSGMT)
     1049        Phi = np.inner(H,SSGT)
    10301050        if SGInv:   #if centro - expand HKL sets
    10311051            Uniq = np.vstack((Uniq,-Uniq))
    1032             SSUniq = np.vstack((SSUniq,-SSUniq))
    10331052            Phi = np.hstack((Phi,-Phi))
    1034             SSPhi = np.hstack((SSPhi,-SSPhi))
    10351053#        GSASIIpath.IPyBreak()
    1036         GfpuA = G2mth.Modulation(waveTypes,SSUniq,FSSdata,XSSdata,USSdata,Mast) #2 x sym X atoms
    1037         phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+SSPhi[:,np.newaxis])
     1054        GfpuA = G2mth.Modulation(waveTypes,Uniq,FSSdata,XSSdata,USSdata,Mast) #2 x sym X atoms
     1055        phase = twopi*(np.inner(Uniq[:,:3],(dXdata.T+Xdata.T))+Phi[:,np.newaxis])
    10381056        sinp = np.sin(phase)
    10391057        cosp = np.cos(phase)
    10401058        biso = -SQfactor*Uisodata
    10411059        Tiso = np.where(biso<1.,np.exp(biso),1.0)
    1042         HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
     1060        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq[:,:3]])
    10431061        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
    10441062        Tcorr = Tiso*Tuij*Mdata*Fdata/len(Uniq)
    10451063        fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr])     #2 x sym x atoms
    1046         fb = np.array([FPP*cosp*Tcorr,(FF+FP-Bab)*sinp*Tcorr])      #swapped around - better?
     1064        fb = np.array([FPP*cosp*Tcorr,(FF+FP-Bab)*sinp*Tcorr])      #swapped around - better
    10471065        fa *= GfpuA
    10481066        fb *= GfpuA       
Note: See TracChangeset for help on using the changeset viewer.