Changeset 489 for trunk/GSASIIstruct.py


Ignore:
Timestamp:
Feb 22, 2012 4:08:40 PM (10 years ago)
Author:
vondreele
Message:

fix to speed up scat fac calculations & skip texture calcs if order=0

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIstruct.py

    r484 r489  
    379379def ShowControls(Controls):
    380380    print ' Least squares controls:'
    381     print ' Derivative type: ',Controls['deriv type']
    382     print ' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
     381    print ' Refinement type: ',Controls['deriv type']
     382    if 'Hessian' in Controls['deriv type']:
     383        print ' Maximum number of cycles:',Controls['max cyc']
     384    else:
     385        print ' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
    383386    print ' Initial shift factor: ','%.3f'%(Controls['shift factor'])
    384387   
     
    616619
    617620               
    618             if 'SH Texture' in General:
    619                 textureData = General['SH Texture']
     621            textureData = General['SH Texture']
     622            if textureData['Order']:
     623                phaseDict[pfx+'SHorder'] = textureData['Order']
    620624                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
    621                 phaseDict[pfx+'SHorder'] = textureData['Order']
    622625                for name in ['omega','chi','phi']:
    623626                    phaseDict[pfx+'SH '+name] = textureData['Sample '+name][1]
     
    640643                    ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
    641644                    '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
    642                 if 'SH Texture' in General:
    643                     PrintTexture(textureData)
     645                PrintTexture(textureData)
    644646                   
    645647        elif PawleyRef:
     
    916918            PrintAtomsAndSig(General,Atoms,atomsSig)
    917919       
    918         if 'SH Texture' in General:
    919             textureData = General['SH Texture']   
    920             if textureData['Order']:
    921                 SHtextureSig = {}
    922                 for name in ['omega','chi','phi']:
    923                     aname = pfx+'SH '+name
    924                     textureData['Sample '+name][1] = parmDict[aname]
    925                     if aname in sigDict:
    926                         SHtextureSig['Sample '+name] = sigDict[aname]
    927                 for name in textureData['SH Coeff'][1]:
    928                     aname = pfx+name
    929                     textureData['SH Coeff'][1][name] = parmDict[aname]
    930                     if aname in sigDict:
    931                         SHtextureSig[name] = sigDict[aname]
    932                 PrintSHtextureAndSig(textureData,SHtextureSig)
     920        textureData = General['SH Texture']   
     921        if textureData['Order']:
     922            SHtextureSig = {}
     923            for name in ['omega','chi','phi']:
     924                aname = pfx+'SH '+name
     925                textureData['Sample '+name][1] = parmDict[aname]
     926                if aname in sigDict:
     927                    SHtextureSig['Sample '+name] = sigDict[aname]
     928            for name in textureData['SH Coeff'][1]:
     929                aname = pfx+name
     930                textureData['SH Coeff'][1][name] = parmDict[aname]
     931                if aname in sigDict:
     932                    SHtextureSig[name] = sigDict[aname]
     933            PrintSHtextureAndSig(textureData,SHtextureSig)
    933934
    934935def GetHistogramPhaseData(Phases,Histograms,Print=True):
     
    11061107                    pos = 2.0*asind(wave/(2.0*d))+Zero
    11071108                    if limits[0] < pos < limits[1]:
    1108                         refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0])
     1109                        refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0,{}])
     1110                        #last item should contain form factor values by atom type
    11091111                else:
    11101112                    raise ValueError
     
    16221624                PrintBackgroundSig(Background,backSig)
    16231625
    1624 def GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict):
     1626def GetAtomFXU(pfx,calcControls,parmDict):
    16251627    Natoms = calcControls['Natoms'][pfx]
    16261628    Tdata = Natoms*[' ',]
     
    16441646            if parm in parmDict:
    16451647                keys[key][iatm] = parmDict[parm]
    1646         FFdata.append(FFtables[Tdata[iatm]])
    1647         BLdata.append(BLtables[Tdata[iatm]][1])
    1648     return FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
    1649    
     1648    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
     1649   
     1650def getFFvalues(FFtables,SQ):
     1651    FFvals = {}
     1652    for El in FFtables:
     1653        FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0]
     1654    return FFvals
     1655   
     1656def getBLvalues(BLtables):
     1657    BLvals = {}
     1658    for El in BLtables:
     1659        BLvals[El] = BLtables[El][1]
     1660    return BLvals
     1661       
    16501662def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
    16511663    ''' Compute structure factors for all h,k,l for phase
     
    16651677    FFtables = calcControls['FFtables']
    16661678    BLtables = calcControls['BLtables']
    1667     FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict)
     1679    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1680    FF = np.zeros(len(Tdata))
    16681681    if 'N' in parmDict[hfx+'Type']:
    1669         FP,FPP = G2el.BlenRes(BLdata,parmDict[hfx+'Lam'])
     1682        FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
    16701683    else:
    1671         FP = np.array([El[hfx+'FP'] for El in FFdata])
    1672         FPP = np.array([El[hfx+'FPP'] for El in FFdata])
     1684        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     1685        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
    16731686    maxPos = len(SGData['SGOps'])
    16741687    Uij = np.array(G2lat.U6toUij(Uijdata))
     
    16781691        H = refl[:3]
    16791692        SQ = 1./(2.*refl[4])**2
    1680         if 'N' in parmDict[hfx+'Type']:
    1681             FF = np.array([El[1] for El in BLdata])
    1682         else:       #'X'
    1683             FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
     1693        if not len(refl[-1]):                #no form factors
     1694            if 'N' in parmDict[hfx+'Type']:
     1695                refl[-1] = getBLvalues(BLtables)
     1696            else:       #'X'
     1697                refl[-1] = getFFvalues(FFtables,SQ)
     1698        for i,El in enumerate(Tdata):           
     1699            FF[i] = refl[-1][El]           
    16841700        SQfactor = 4.0*SQ*twopisq
    16851701        Uniq = refl[11]
     
    17121728    FFtables = calcControls['FFtables']
    17131729    BLtables = calcControls['BLtables']
    1714     FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict)
     1730    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1731    FF = np.zeros(len(Tdata))
    17151732    if 'N' in parmDict[hfx+'Type']:
    17161733        FP = 0.
    17171734        FPP = 0.
    17181735    else:
    1719         FP = np.array([El[hfx+'FP'] for El in FFdata])
    1720         FPP = np.array([El[hfx+'FPP'] for El in FFdata])
     1736        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
     1737        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
    17211738    maxPos = len(SGData['SGOps'])       
    17221739    Uij = np.array(G2lat.U6toUij(Uijdata))
     
    17291746    for iref,refl in enumerate(refList):
    17301747        H = np.array(refl[:3])
    1731         SQ = 1./(2.*refl[4])**2          # or (sin(theta)/lambda)**2
    1732         if 'N' in parmDict[hfx+'Type']:
    1733             FF = np.array([El[1] for El in BLdata])
    1734         else:       #'X'
    1735             FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
     1748        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
     1749        for i,El in enumerate(Tdata):           
     1750            FF[i] = refl[-1][El]           
    17361751        SQfactor = 8.0*SQ*np.pi**2
    17371752        Uniq = refl[11]
     
    19161931   
    19171932def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
     1933    POcorr = 1.0
    19181934    if calcControls[phfx+'poType'] == 'MD':
    19191935        MD = parmDict[phfx+'MD']
    1920         MDAxis = calcControls[phfx+'MDAxis']
    1921         sumMD = 0
    1922         for H in refl[11]:           
    1923             cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
    1924             A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
    1925             sumMD += A**3
    1926         POcorr = sumMD/len(refl[11])
     1936        if MD != 1.0:
     1937            MDAxis = calcControls[phfx+'MDAxis']
     1938            sumMD = 0
     1939            for H in refl[11]:           
     1940                cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
     1941                A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
     1942                sumMD += A**3
     1943            POcorr = sumMD/len(refl[11])
    19271944    else:   #spherical harmonics
    1928         POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
     1945        if calcControls[pfx+'SHord']:
     1946            POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
    19291947    return POcorr
    19301948   
    19311949def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
     1950    POcorr = 1.0
    19321951    POderv = {}
    19331952    if calcControls[phfx+'poType'] == 'MD':
     
    19441963        POderv[phfx+'MD'] = sumdMD/len(refl[11])
    19451964    else:   #spherical harmonics
    1946         POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
     1965        if calcControls[pfx+'SHord']:
     1966            POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
    19471967    return POcorr,POderv
    19481968   
     
    21702190                FFtables[El][hfx+'FPP'] = FPP               
    21712191           
     2192def GetFobsSq(Histograms,Phases,parmDict,calcControls):
     2193    histoList = Histograms.keys()
     2194    histoList.sort()
     2195    for histogram in histoList:
     2196        if 'PWDR' in histogram[:4]:
     2197            Histogram = Histograms[histogram]
     2198            hId = Histogram['hId']
     2199            hfx = ':%d:'%(hId)
     2200            Limits = calcControls[hfx+'Limits']
     2201            shl = max(parmDict[hfx+'SH/L'],0.002)
     2202            Ka2 = False
     2203            kRatio = 0.0
     2204            if hfx+'Lam1' in parmDict.keys():
     2205                Ka2 = True
     2206                lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
     2207                kRatio = parmDict[hfx+'I(L2)/I(L1)']
     2208            x,y,w,yc,yb,yd = Histogram['Data']
     2209            ymb = np.array(y-yb)
     2210            ycmb = np.array(yc-yb)
     2211            ratio = np.where(ycmb!=0.,ymb/ycmb,0.0)           
     2212            xB = np.searchsorted(x,Limits[0])
     2213            xF = np.searchsorted(x,Limits[1])
     2214            refLists = Histogram['Reflection Lists']
     2215            for phase in refLists:
     2216                Phase = Phases[phase]
     2217                pId = Phase['pId']
     2218                phfx = '%d:%d:'%(pId,hId)
     2219                refList = refLists[phase]
     2220                sumFo = 0.0
     2221                sumdF = 0.0
     2222                sumFosq = 0.0
     2223                sumdFsq = 0.0
     2224                for refl in refList:
     2225                    if 'C' in calcControls[hfx+'histType']:
     2226                        yp = np.zeros_like(yb)
     2227                        Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
     2228                        iBeg = np.searchsorted(x[xB:xF],refl[5]-fmin)
     2229                        iFin = np.searchsorted(x[xB:xF],refl[5]+fmax)
     2230                        iFin2 = iFin
     2231                        yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here                           
     2232                        if Ka2:
     2233                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
     2234                            Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
     2235                            iBeg2 = np.searchsorted(x,pos2-fmin)
     2236                            iFin2 = np.searchsorted(x,pos2+fmax)
     2237                            yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
     2238                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0))
     2239                    elif 'T' in calcControls[hfx+'histType']:
     2240                        print 'TOF Undefined at present'
     2241                        raise Exception    #no TOF yet
     2242                    Fo = np.sqrt(np.abs(refl[8]))
     2243                    Fc = np.sqrt(np.abs(refl[9]))
     2244                    sumFo += Fo
     2245                    sumFosq += refl[8]**2
     2246                    sumdF += np.abs(Fo-Fc)
     2247                    sumdFsq += (refl[8]-refl[9])**2
     2248                Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
     2249                Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
     2250                Histogram[phfx+'Nref'] = len(refList)
     2251               
    21722252def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
    21732253   
     
    22562336    return yc,yb
    22572337   
    2258 def GetFobsSq(Histograms,Phases,parmDict,calcControls):
    2259     histoList = Histograms.keys()
    2260     histoList.sort()
    2261     for histogram in histoList:
    2262         if 'PWDR' in histogram[:4]:
    2263             Histogram = Histograms[histogram]
    2264             hId = Histogram['hId']
    2265             hfx = ':%d:'%(hId)
    2266             Limits = calcControls[hfx+'Limits']
    2267             shl = max(parmDict[hfx+'SH/L'],0.002)
    2268             Ka2 = False
    2269             kRatio = 0.0
    2270             if hfx+'Lam1' in parmDict.keys():
    2271                 Ka2 = True
    2272                 lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
    2273                 kRatio = parmDict[hfx+'I(L2)/I(L1)']
    2274             x,y,w,yc,yb,yd = Histogram['Data']
    2275             ymb = np.array(y-yb)
    2276             ycmb = np.array(yc-yb)
    2277             ratio = np.where(ycmb!=0.,ymb/ycmb,0.0)           
    2278             xB = np.searchsorted(x,Limits[0])
    2279             xF = np.searchsorted(x,Limits[1])
    2280             refLists = Histogram['Reflection Lists']
    2281             for phase in refLists:
    2282                 Phase = Phases[phase]
    2283                 pId = Phase['pId']
    2284                 phfx = '%d:%d:'%(pId,hId)
    2285                 refList = refLists[phase]
    2286                 sumFo = 0.0
    2287                 sumdF = 0.0
    2288                 sumFosq = 0.0
    2289                 sumdFsq = 0.0
    2290                 for refl in refList:
    2291                     if 'C' in calcControls[hfx+'histType']:
    2292                         yp = np.zeros_like(yb)
    2293                         Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
    2294                         iBeg = np.searchsorted(x[xB:xF],refl[5]-fmin)
    2295                         iFin = np.searchsorted(x[xB:xF],refl[5]+fmax)
    2296                         iFin2 = iFin
    2297                         yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here                           
    2298                         if Ka2:
    2299                             pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
    2300                             Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
    2301                             iBeg2 = np.searchsorted(x,pos2-fmin)
    2302                             iFin2 = np.searchsorted(x,pos2+fmax)
    2303                             yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
    2304                         refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0))
    2305                     elif 'T' in calcControls[hfx+'histType']:
    2306                         print 'TOF Undefined at present'
    2307                         raise Exception    #no TOF yet
    2308                     Fo = np.sqrt(np.abs(refl[8]))
    2309                     Fc = np.sqrt(np.abs(refl[9]))
    2310                     sumFo += Fo
    2311                     sumFosq += refl[8]**2
    2312                     sumdF += np.abs(Fo-Fc)
    2313                     sumdFsq += (refl[8]-refl[9])**2
    2314                 Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
    2315                 Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
    2316                 Histogram[phfx+'Nref'] = len(refList)
    2317                
    23182338def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
    23192339   
Note: See TracChangeset for help on using the changeset viewer.