Changeset 1125


Ignore:
Timestamp:
Oct 28, 2013 3:32:05 PM (10 years ago)
Author:
vondreele
Message:

StructureFactor2 does SF calculation for entire block of reflection - see substantial speed gains in some cases.
The old StructureFactor? routine is left in as a reference - modified form factor stuff

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASII.py

    r1121 r1125  
    14101410        self.setPoly = False
    14111411        self.setFrame = False
     1412        self.setSpot = False
     1413        self.setArc = False
     1414        self.setRing = False
    14121415        arg = sys.argv
    14131416        if len(arg) > 1:

  • trunk/GSASIIElem.py

    r967 r1125  
    100100    return FFvals
    101101   
    102 def getBLvalues(BLtables):
     102def getBLvalues(BLtables,ifList=False):
    103103    'Needs a doc string'
    104     BLvals = {}
    105     for El in BLtables:
    106         BLvals[El] = BLtables[El][1][1]
     104    if ifList:
     105        BLvals = []
     106        for El in BLtables:
     107            BLvals.append(BLtables[El][1][1])
     108    else:
     109        BLvals = {}
     110        for El in BLtables:
     111            BLvals[El] = BLtables[El][1][1]
    107112    return BLvals
    108113       

  • trunk/GSASIIstrIO.py

    r1110 r1125  
    297297#patch
    298298                if isinstance(datum[1][1],list):
    299                     RefData = {'RefList':[],'FF':[]}
     299                    RefData = {'RefList':[],'FF':{}}
    300300                    for ref in datum[1][1]:
    301301                        RefData['RefList'].append(ref[:11]+[ref[13],])
    302                         RefData['FF'].append(ref[14])
    303302                    RefData['RefList'] = np.array(RefData['RefList'])
    304303                    datum[1][1] = RefData
    305304#end patch
     305                datum[1][1]['FF'] = {}
    306306                HKLFdata['Data'] = datum[1][1]
    307307                HKLFdata[data[1][0]] = data[1][1]       #Instrument parameters
    … …  
    17561756                    Uniq = []
    17571757                    Phi = []
    1758                     FF = []
    17591758                    for h,k,l,d in HKLd:
    17601759                        ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
    … …  
    17681767                                Uniq.append(uniq)
    17691768                                Phi.append(phi)
    1770                                 FF.append({})
    17711769                        else:
    17721770                            raise ValueError
    1773                     Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':FF}
     1771                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{}}
    17741772            elif 'HKLF' in histogram:
    17751773                inst = Histogram['Instrument Parameters'][0]

  • trunk/GSASIIstrMain.py

    r1077 r1125  
    260260        parmDict = {}
    261261        Histo = {histogram:Histograms[histogram],}
    262         hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,False)
     262        hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,Print=False)
    263263        calcControls.update(controlDict)
    264264        histVary,histDict,controlDict = G2stIO.GetHistogramData(Histo,False)

  • trunk/GSASIIstrMath.py

    r1119 r1125  
    561561    Uij = np.array(G2lat.U6toUij(Uijdata))
    562562    bij = Mast*Uij.T
     563    if not len(refDict['FF']):
     564        if 'N' in calcControls[hfx+'histType']:
     565            dat = G2el.getBLvalues(BLtables)
     566        else:
     567            dat = G2el.getFFvalues(FFtables,0.)       
     568        refDict['FF']['El'] = dat.keys()
     569        refDict['FF']['FF'] = np.zeros((len(refDict['RefList']),len(dat)))   
    563570    for iref,refl in enumerate(refDict['RefList']):
    564571        fbs = np.array([0,0])
    … …  
    567574        SQfactor = 4.0*SQ*twopisq
    568575        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
    569         if not len(refDict['FF'][iref]):                #no form factors
     576        if not np.any(refDict['FF']['FF'][iref]):                #no form factors - 1st time thru StructureFactor
    570577            if 'N' in calcControls[hfx+'histType']:
    571                 refDict['FF'][iref] = G2el.getBLvalues(BLtables)
     578                dat = G2el.getBLvalues(BLtables)
     579                refDict['FF']['FF'][iref] = dat.values()
    572580            else:       #'X'
    573                 refDict['FF'][iref] = G2el.getFFvalues(FFtables,SQ)
    574         FF = [refDict['FF'][iref][El] for El in Tdata]         
     581                dat = G2el.getFFvalues(FFtables,SQ)
     582                refDict['FF']['FF'][iref] = dat.values()
     583        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     584        FF = refDict['FF']['FF'][iref][Tindx]
    575585        Uniq = np.inner(H,SGMT)
    576586        Phi = np.inner(H,SGT)
    … …  
    578588        sinp = np.sin(phase)
    579589        cosp = np.cos(phase)
    580         occ = Mdata*Fdata/len(Uniq)
    581590        biso = -SQfactor*Uisodata
    582591        Tiso = np.where(biso<1.,np.exp(biso),1.0)
    583592        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
    584593        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
    585         Tcorr = Tiso*Tuij
    586         fa = np.array([(FF+FP-Bab)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
     594        Tcorr = Tiso*Tuij*Mdata*Fdata/len(Uniq)
     595        fa = np.array([(FF+FP-Bab)*cosp*Tcorr,-FPP*sinp*Tcorr])
    587596        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
    588597        if not SGData['SGInv']:
    589             fb = np.array([(FF+FP-Bab)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
     598            fb = np.array([(FF+FP-Bab)*sinp*Tcorr,FPP*cosp*Tcorr])
    590599            fbs = np.sum(np.sum(fb,axis=1),axis=1)
    591600        fasq = fas**2
    … …  
    593602        refl[9] = np.sum(fasq)+np.sum(fbsq)
    594603        refl[10] = atan2d(fbs[0],fas[0])
     604   
     605def StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
     606    ''' Compute structure factors for all h,k,l for phase
     607    puts the result, F^2, in each ref[8] in refList
     608    input:
     609   
     610    :param dict refDict: where
     611        'RefList' list where each ref = h,k,l,m,d,...
     612        'FF' dict of form factors - filed in below
     613    :param np.array G:      reciprocal metric tensor
     614    :param str pfx:    phase id string
     615    :param dict SGData: space group info. dictionary output from SpcGroup
     616    :param dict calcControls:
     617    :param dict ParmDict:
     618
     619    '''       
     620    twopi = 2.0*np.pi
     621    twopisq = 2.0*np.pi**2
     622    phfx = pfx.split(':')[0]+hfx
     623    ast = np.sqrt(np.diag(G))
     624    Mast = twopisq*np.multiply.outer(ast,ast)
     625    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
     626    SGT = np.array([ops[1] for ops in SGData['SGOps']])
     627    FFtables = calcControls['FFtables']
     628    BLtables = calcControls['BLtables']
     629    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     630    FF = np.zeros(len(Tdata))
     631    if 'N' in calcControls[hfx+'histType']:
     632        FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
     633    else:
     634        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     635        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
     636    Uij = np.array(G2lat.U6toUij(Uijdata))
     637    bij = Mast*Uij.T
     638    if not len(refDict['FF']):                #no form factors - 1st time thru StructureFactor
     639        if 'N' in calcControls[hfx+'histType']:
     640            dat = G2el.getBLvalues(BLtables)
     641            refDict['FF']['El'] = dat.keys()
     642            refDict['FF']['FF'] = np.ones((len(refl),len(dat)))*dat.values()           
     643        else:       #'X'
     644            dat = G2el.getFFvalues(FFtables,0.)
     645            refDict['FF']['El'] = dat.keys()
     646            refDict['FF']['FF'] = np.ones((len(refDict['RefList']),len(dat)))
     647        refDict['ifNew'] = True           
     648#reflection processing begins here - big arrays!
     649    refl = refDict['RefList']
     650    H = refl.T[:3]
     651    SQ = 1./(2.*refl.T[4])**2
     652    SQfactor = 4.0*SQ*twopisq
     653    Bab = np.repeat(parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor),len(SGT))
     654    if refDict['ifNew']:                #no form factors - 1st time thru StructureFactor
     655        for iref in range(len(refl)):
     656            if 'X' in calcControls[hfx+'histType']:
     657                dat = G2el.getFFvalues(FFtables,SQ[iref])
     658                refDict['FF']['FF'][iref] *= dat.values()
     659        refDict['ifNew'] = False
     660    Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     661    FF = np.repeat(refDict['FF']['FF'].T[Tindx].T,len(SGT),axis=0)
     662    Uniq = np.reshape(np.inner(H.T,SGMT),(-1,3))
     663    Phi = np.inner(H.T,SGT).flatten()
     664    phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T)+Phi[:,np.newaxis])
     665    sinp = np.sin(phase)
     666    cosp = np.cos(phase)
     667    biso = -SQfactor*Uisodata[:,np.newaxis]
     668    Tiso = np.repeat(np.where(biso<1.,np.exp(biso),1.0),len(SGT),axis=1).T
     669    HbH = -np.sum(Uniq.T*np.inner(bij,Uniq),axis=1)
     670    Tuij = np.where(HbH<1.,np.exp(HbH),1.0).T
     671    Tcorr = Tiso*Tuij*Mdata*Fdata/len(SGMT)
     672    fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-FPP*sinp*Tcorr])
     673    fa = np.reshape(fa,(2,len(refl),len(SGT),len(Mdata)))
     674    fas = np.sum(np.sum(fa,axis=2),axis=2)        #real
     675    fbs = np.zeros_like(fas)
     676    if not SGData['SGInv']:
     677        fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,FPP*cosp*Tcorr])
     678        fb = np.reshape(fb,(2,len(refl),len(SGT),len(Mdata)))
     679        fbs = np.sum(np.sum(fb,axis=2),axis=2)
     680    fasq = fas**2
     681    fbsq = fbs**2        #imaginary
     682    refl.T[9] = np.sum(fasq,axis=0)+np.sum(fbsq,axis=0)
     683    refl.T[10] = atan2d(fbs[0],fas[0])
    595684   
    596685def StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
    … …  
    628717        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
    629718        Bab = parmDict[phfx+'BabA']*dBabdA
    630         FF = [refDict['FF'][iref][El] for El in Tdata]         
     719        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
     720        FF = refDict['FF']['FF'][iref].T[Tindx]
     721#        FF = [refDict['FF'][iref][El] for El in Tdata]         
    631722        Uniq = np.inner(H,SGMT)
    632723        Phi = np.inner(H,SGT)
    … …  
    13081399        if not Phase['General'].get('doPawley'):
    13091400            time0 = time.time()
    1310             StructureFactor(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     1401            StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    13111402            print 'sf calc time: %.3fs'%(time.time()-time0)
    13121403        time0 = time.time()
    … …  
    17891880            refDict = Histogram['Data']
    17901881            time0 = time.time()
    1791             dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)  #accurate for powders!
     1882            dFdvDict = StructureFactorDerv(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     1883            print 'sf-deriv time: %.3f'%(time.time()-time0)
    17921884            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
    17931885            dMdvh = np.zeros((len(varylist),len(refDict['RefList'])))
    … …  
    19432035            refDict = Histogram['Data']
    19442036            time0 = time.time()
    1945             StructureFactor(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    1946 #            print 'sf-calc time: %.3f'%(time.time()-time0)
     2037            StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
     2038            print 'sf-calc time: %.3f'%(time.time()-time0)
    19472039            df = np.zeros(len(refDict['RefList']))
    19482040            sumwYo = 0
Note: See TracChangeset for help on using the changeset viewer.