Changeset 942


Ignore:
Timestamp:
Jun 5, 2013 9:26:36 AM (10 years ago)
Author:
vondreele
Message:

mods for MC/SA:
moved scat fac routines from GSASIIstrIO.py & GSASIIstrMath.py to GSASIIElem.py

Location:
trunk
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIElem.py

    r939 r942  
    5959    return FormFactors
    6060   
     61def GetFFtable(atomTypes):
     62    ''' returns a dictionary of form factor data for atom types found in atomTypes
     63
     64    :param list atomTypes: list of atom types
     65    :return: FFtable, dictionary of form factor data; key is atom type
     66
     67    '''
     68    FFtable = {}
     69    for El in atomTypes:
     70        FFs = GetFormFactorCoeff(El.split('+')[0].split('-')[0])
     71        for item in FFs:
     72            if item['Symbol'] == El.upper():
     73                FFtable[El] = item
     74    return FFtable
     75   
     76def GetBLtable(General):
     77    ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
     78
     79    :param dict General: dictionary of phase info.; includes AtomTypes & Isotopes
     80    :return: BLtable, dictionary of scattering length data; key is atom type
     81    '''
     82    atomTypes = General['AtomTypes']
     83    BLtable = {}
     84    isotopes = General['Isotopes']
     85    isotope = General['Isotope']
     86    for El in atomTypes:
     87        BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]
     88    return BLtable
     89       
     90def getFFvalues(FFtables,SQ):
     91    'Needs a doc string'
     92    FFvals = {}
     93    for El in FFtables:
     94        FFvals[El] = ScatFac(FFtables[El],SQ)[0]
     95    return FFvals
     96   
     97def getBLvalues(BLtables):
     98    'Needs a doc string'
     99    BLvals = {}
     100    for El in BLtables:
     101        BLvals[El] = BLtables[El][1][1]
     102    return BLvals
     103       
    61104def GetFFC5(ElSym):
    62105    '''Get 5 term form factor and Compton scattering data

  • trunk/GSASIIIO.py

    r941 r942  
    5757        inst[item] = list(inst[item])
    5858    return inst
    59 
    6059
    6160def FileDlgFixExt(dlg,file):

  • trunk/GSASIImath.py

    r940 r942  
    19001900
    19011901
    1902 def mcsaSearch(data,reflType,reflData,covData,pgbar):
     1902def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
    19031903    gamFW = lambda s,g: math.exp(math.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.)
    19041904   
    … …  
    19121912            upper.append(limits[1])
    19131913           
    1914     def getAtomparms(item,pfx,parmDict,varyList):
     1914    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
    19151915        parmDict[pfx+'Atype'] = item['atType']
    1916         pstr = ['x','y','z']
     1916        aTypes |= set([item['atType'],])
     1917        pstr = ['Ax','Ay','Az']
     1918        XYZ = [0,0,0]
    19171919        for i in range(3):
    19181920            name = pfx+pstr[i]
    19191921            parmDict[name] = item['Pos'][0][i]
     1922            XYZ[i] = parmDict[name]
    19201923            if item['Pos'][1][i]:
    19211924                varyList += [name,]
    … …  
    19231926                lower.append(limits[0])
    19241927                upper.append(limits[1])
     1928        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
    19251929           
    1926     def getRBparms(item,pfx,parmDict,varyList):
    1927         parmDict['MolCent'] = item['MolCent']
    1928         parmDict['RBId'] = item['RBId']
     1930    def getRBparms(item,mfx,aTypes,RBdata,parmDict,varyList):
     1931        parmDict[mfx+'MolCent'] = item['MolCent']
     1932        parmDict[mfx+'RBId'] = item['RBId']
    19291933        pstr = ['Px','Py','Pz']
    19301934        ostr = ['Qa','Qi','Qj','Qk']
    19311935        for i in range(3):
    1932             name = pfx+pstr[i]
     1936            name = mfx+pstr[i]
    19331937            parmDict[name] = item['Pos'][0][i]
    19341938            if item['Pos'][1][i]:
    … …  
    19381942                upper.append(limits[1])
    19391943        for i in range(4):
    1940             name = pfx+ostr[i]
     1944            name = mfx+ostr[i]
    19411945            parmDict[name] = item['Ori'][0][i]
    19421946            if item['Ovar'] == 'AV' and i:
    … …  
    19521956        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
    19531957            for i in range(len(item['Tor'][0])):
    1954                 name = pfx+'Tor'+str(i)
     1958                name = mfx+'Tor'+str(i)
    19551959                parmDict[name] = item['Tor'][0][i]
    19561960                if item['Tor'][1][i]:
    … …  
    19591963                    lower.append(limits[0])
    19601964                    upper.append(limits[1])
    1961 
    1962     def getFFvalues(FFtables,SQ):
    1963         FFvals = {}
    1964         for El in FFtables:
    1965             FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0]
    1966         return FFvals
    1967        
    1968     def getBLvalues(BLtables):
    1969         BLvals = {}
    1970         for El in BLtables:
    1971             BLvals[El] = BLtables[El][1][1]
    1972         return BLvals
    1973            
    1974 #    def mcsaSfCalc(refList,G,SGData,calcControls,parmDict):
     1965        aTypes |= set(RBdata[item['Type']][item['RBId']]['rbTypes'])
     1966
     1967    def GetAtomFX(pfx,Natoms,parmDict):
     1968        'Needs a doc string'
     1969        Tdata = Natoms*[' ',]
     1970        Mdata = np.zeros(Natoms)
     1971        Fdata = np.zeros(Natoms)
     1972        Xdata = np.zeros((3,Natoms))
     1973        keys = {'Atype:':Tdata,'Amul:':Mdata,
     1974            'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2]}
     1975        for iatm in range(Natoms):
     1976            for key in keys:
     1977                parm = pfx+key+str(iatm)
     1978                if parm in parmDict:
     1979                    keys[key][iatm] = parmDict[parm]
     1980        return Tdata,Mdata,Xdata
     1981   
     1982#    def mcsaSfCalc(refList,G,SGData,parmDict):
    19751983#        ''' Compute structure factors for all h,k,l for phase
    19761984#        input:
    … …  
    19781986#            G:      reciprocal metric tensor
    19791987#            SGData: space group info. dictionary output from SpcGroup
    1980 #            calcControls:
    19811988#            ParmDict:
    19821989#        puts result F^2 in each ref[8] in refList
    … …  
    19861993#        ast = np.sqrt(np.diag(G))
    19871994#        Mast = twopisq*np.multiply.outer(ast,ast)
    1988 #        FFtables = calcControls['FFtables']
    1989 #        BLtables = calcControls['BLtables']
    1990 #        Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
     1995#        Tdata,Mdata,Fdata,Xdata = GetAtomFX(pfx,calcControls,parmDict)
    19911996#        FF = np.zeros(len(Tdata))
    19921997#        if 'N' in parmDict[hfx+'Type']:
    … …  
    19952000#            FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
    19962001#            FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
    1997 #        maxPos = len(SGData['SGOps'])
    19982002#        Uij = np.array(G2lat.U6toUij(Uijdata))
    19992003#        bij = Mast*Uij.T
    … …  
    20132017#            Uniq = refl[11]
    20142018#            phi = refl[12]
    2015 #            phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
     2019#            phase = twopi*(np.inner(Uniq,(Xdata.T))+phi[:,np.newaxis])
    20162020#            sinp = np.sin(phase)
    20172021#            cosp = np.cos(phase)
    20182022#            occ = Mdata*Fdata/len(Uniq)
    2019 #            biso = -SQfactor*Uisodata
    2020 #            Tiso = np.where(biso<1.,np.exp(biso),1.0)
    2021 #            HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
    2022 #            Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
    2023 #            Tcorr = Tiso*Tuij
    2024 #            fa = np.array([(FF+FP-Bab)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
     2023#            fa = np.array([(FF+FP-Bab)*occ*cosp,-FPP*occ*sinp])
    20252024#            fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
    20262025#            if not SGData['SGInv']:
    … …  
    20362035    sq4pi = np.sqrt(4*np.pi)
    20372036    generalData = data['General']
     2037    SGData = generalData['SGData']
    20382038    pId = data['pId']
    2039     pfx = str(pId)+'::'
    2040     Atoms = data['Atoms']                       #if any
     2039    fixAtoms = data['Atoms']                       #if any
     2040    cx,ct,cs = generalData['AtomPtrs'][:3]
     2041    aTypes = set([])
     2042    parmDict = {}
     2043    varyList = []
     2044    atNo = 0
     2045    for atm in fixAtoms:
     2046        pfx = ':'+str(atNo)+':'
     2047        parmDict[pfx+'Atype'] = atm[ct]
     2048        aTypes |= set([atm[ct],])
     2049        pstr = ['Ax','Ay','Az']
     2050        parmDict[pfx+'Amul'] = atm[cs+1]
     2051        for i in range(3):
     2052            name = pfx+pstr[i]
     2053            parmDict[name] = atm[cx+i]
     2054        atNo += 1       
    20412055    mcsaControls = generalData['MCSA controls']
    20422056    reflName = mcsaControls['Data source']
    20432057    phaseName = generalData['Name']
    20442058    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
    2045     parmDict = {}
    2046     varyList = []
    20472059    upper = []
    20482060    lower = []
    … …  
    20522064            getMDparms(item,mfx,parmDict,varyList)
    20532065        elif item['Type'] == 'Atom':
    2054             getAtomparms(item,mfx,parmDict,varyList)
     2066            pfx = mfx+str(atNo)+':'
     2067            getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList)
     2068            atNo += 1
    20552069        elif item['Type'] in ['Residue','Vector']:
    2056             getRBparms(item,mfx,parmDict,varyList)
     2070            pfx = mfx+':'
     2071            getRBparms(item,pfx,aTypes,RBdata,parmDict,varyList)
     2072    FFtables = G2el.GetFFtable(aTypes)
     2073    refs = []
    20572074    if 'PWDR' in reflName:
    2058         refs = []
    20592075        for ref in reflData:
    20602076            h,k,l,m,d,pos,sig,gam,Fsq = ref[:9]
    20612077            if d >= mcsaControls['dmin']:
    2062                 FWHM = gamFW(sig,gam)/2.35482        #sqrt(8ln2) --> sig from FWHM
    2063                 Fsq *= m
     2078                sig = gamFW(sig,gam)/sq8ln2        #--> sig from FWHM
     2079                SQ = 0.25/d**2
     2080                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
     2081                FFs = G2el.getFFvalues(FFtables,SQ)
     2082                refs.append([h,k,l,m,f*m,pos,sig,FFs,Uniq,phi])
     2083        nRef = len(refs)
     2084        rcov = np.zeros((nRef,nRef))
     2085        for iref,refI in enumerate(refs):
     2086            rcov[iref][iref] = 1./(sq4pi*refI[6])
     2087            for jref,refJ in enumerate(refs[:iref]):
     2088                t1 = refI[6]**2+refJ[6]**2
     2089                t2 = (refJ[5]-refI[5])**2/(2.*t1)
     2090                if t2 > 10.:
     2091                    rcov[iref][jref] = 0.
     2092                else:
     2093                    rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2))
     2094        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
     2095        Rdiag = np.sqrt(np.diag(rcov))
     2096        Rnorm = np.outer(Rdiag,Rdiag)
     2097        rcov /= Rnorm
    20642098    elif 'Pawley' in reflName:
    2065         refs = []
    20662099        covMatrix = covData['covMatrix']
    2067         varyList = covData['varyList']
     2100        vList = covData['varyList']
    20682101        for iref,refI in enumerate(reflData):
    20692102            h,k,l,m,d,v,f,s = refI
    20702103            if d >= mcsaControls['dmin'] and v:       #skip unrefined ones
    2071                 nameI = pfx+'PWLref:'+str(iref)
    2072                 if nameI in covData['varyList']:
    2073                     refs.append([h,k,l,f*m,0.,0.])
     2104                SQ = 0.25/d**2
     2105                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
     2106                FFs = G2el.getFFvalues(FFtables,SQ)
     2107                refs.append([h,k,l,m,f*m,0.,0.,FFs,Uniq,phi])
    20742108        nRef = len(refs)
    2075         print nRef       
    2076         covTerms = np.zeros((nRef,nRef))
    2077         print covTerms.shape
    2078        
     2109        rcov = np.zeros((nRef,nRef))       
    20792110        Iref = 0
    20802111        for iref,refI in enumerate(reflData):
    … …  
    20822113                nameI = pfx+'PWLref:'+str(iref)
    20832114                if nameI in covData['varyList']:
    2084                     Iindx = varyList.index(nameI)
    2085                     covTerms[Iref][Iref] = covMatrix[Iindx][Iindx]
     2115                    Iindx = vList.index(nameI)
     2116                    rcov[Iref][Iref] = covMatrix[Iindx][Iindx]
    20862117                    Jref = 0
    20872118                    for jref,refJ in enumerate(reflData[:iref]):
    … …  
    20892120                            nameJ = pfx+'PWLref:'+str(jref)
    20902121                            try:
    2091                                 covTerms[Iref][Jref] = covMatrix[varyList.index(nameI)][varyList.index(nameJ)]
     2122                                rcov[Iref][Jref] = covMatrix[vList.index(nameI)][vList.index(nameJ)]
    20922123                            except ValueError:
    2093                                 covTerms[Iref][Jref] = covTerms[Iref][Jref-1]
     2124                                rcov[Iref][Jref] = rcov[Iref][Jref-1]
    20942125                            Jref += 1
     2126                else:
     2127                    rcov[Iref] = rcov[Iref-1]
     2128                    rcov[Iref][Iref] = rcov[Iref-1][Iref-1]
    20952129                Iref += 1
    2096         print covTerms
     2130        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
     2131        Rdiag = np.sqrt(np.diag(rcov))
     2132        Rnorm = np.outer(Rdiag,Rdiag)
     2133        rcov /= Rnorm
     2134    elif 'HKLF' in reflName:
     2135        for ref in reflData:
     2136            [h,k,l,m,d],f = ref[:5],ref[6]
     2137            if d >= mcsaControls['dmin']:
     2138                SQ = 0.25/d**2
     2139                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
     2140                FFs = G2el.getFFvalues(FFtables,SQ)
     2141                refs.append([h,k,l,m,f*m,0.,0.,FFs,Uniq,phi])
     2142        rcov = np.identity(len(refs))
     2143       
     2144    for parm in parmDict:
     2145        print parm,parmDict[parm]
    20972146                       
    2098                    
    2099 #    covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
    2100 #        'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
    2101     elif 'HKLF' in reflName:
    2102         refs = []
    2103         for ref in reflData:
    2104             h,k,l,m,d,Fsq = ref[:5],ref[6]
    2105             if d >= mcsaControls['dmin']:
    2106                 Fsq *= m
    2107                 refs.append([h,k,l,m,Fsq])
    2108         rcov = np.identity(len(refs))   
    2109                        
    2110        
     2147#    XYZ,aTypes = UpdateMCSAxyz(Bmat,MCSA)       
    21112148           
    21122149#    generalData['MCSA controls'] = {'Data source':'','Annealing':[50.,0.001,50,1.e-6],

  • trunk/GSASIIphsGUI.py

    r941 r942  
    43704370            print '**** ERROR - No data defined for MC/SA run'
    43714371            return
     4372        RBdata = G2frame.PatternTree.GetItemPyData(   
     4373            G2gd.GetPatternTreeItemId(G2frame,G2frame.root,'Rigid bodies'))
    43724374        MCSAmodels = MCSAdata['Models']
    43734375        if not len(MCSAmodels):
    … …  
    43824384        pgbar.SetSize(Size)
    43834385        try:
    4384             MCSAdata['Results'] = G2mth.mcsaSearch(data,reflType,reflData,covData,pgbar)
     4386            MCSAdata['Results'] = G2mth.mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar)
    43854387        finally:
    43864388            pgbar.Destroy()
    … …  
    49414943        page = event.GetSelection()
    49424944        text = G2frame.dataDisplay.GetPageText(page)
     4945        Atoms.ClearGrid()
    49434946        if text == 'Atoms':
    49444947            G2gd.SetDataMenuBar(G2frame,G2frame.dataFrame.AtomsMenu)

  • trunk/GSASIIstrIO.py

    r939 r942  
    483483    print >>pFile,' Initial shift factor: ','%.3f'%(Controls['shift factor'])
    484484   
    485 def GetFFtable(General):
    486     ''' returns a dictionary of form factor data for atom types found in General
    487 
    488     :param dict General: dictionary of phase info.; includes AtomTypes
    489     :return: FFtable, dictionary of form factor data; key is atom type
    490 
    491     '''
    492     atomTypes = General['AtomTypes']
    493     FFtable = {}
    494     for El in atomTypes:
    495         FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
    496         for item in FFs:
    497             if item['Symbol'] == El.upper():
    498                 FFtable[El] = item
    499     return FFtable
    500    
    501 def GetBLtable(General):
    502     ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
    503 
    504     :param dict General: dictionary of phase info.; includes AtomTypes & Isotopes
    505     :return: BLtable, dictionary of scattering length data; key is atom type
    506     '''
    507     atomTypes = General['AtomTypes']
    508     BLtable = {}
    509     isotopes = General['Isotopes']
    510     isotope = General['Isotope']
    511     for El in atomTypes:
    512         BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]
    513     return BLtable
    514        
    515485def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
    516486    'needs a doc string'
    … …  
    900870        pId = PhaseData[name]['pId']
    901871        pfx = str(pId)+'::'
    902         FFtable = GetFFtable(General)
    903         BLtable = GetBLtable(General)
     872        FFtable = G2el.GetFFtable(General['AtomTypes'])
     873        BLtable = G2el.GetBLtable(General)
    904874        FFtables.update(FFtable)
    905875        BLtables.update(BLtable)

  • trunk/GSASIIstrMath.py

    r939 r942  
    516516    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
    517517   
    518 def getFFvalues(FFtables,SQ):
    519     'Needs a doc string'
    520     FFvals = {}
    521     for El in FFtables:
    522         FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0]
    523     return FFvals
    524    
    525 def getBLvalues(BLtables):
    526     'Needs a doc string'
    527     BLvals = {}
    528     for El in BLtables:
    529         BLvals[El] = BLtables[El][1][1]
    530     return BLvals
    531        
    532518def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
    533519    ''' Compute structure factors for all h,k,l for phase
    … …  
    557543        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
    558544        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
    559     maxPos = len(SGData['SGOps'])
    560545    Uij = np.array(G2lat.U6toUij(Uijdata))
    561546    bij = Mast*Uij.T
    … …  
    568553        if not len(refl[-1]):                #no form factors
    569554            if 'N' in parmDict[hfx+'Type']:
    570                 refl[-1] = getBLvalues(BLtables)
     555                refl[-1] = G2el.getBLvalues(BLtables)
    571556            else:       #'X'
    572                 refl[-1] = getFFvalues(FFtables,SQ)
     557                refl[-1] = G2el.getFFvalues(FFtables,SQ)
    573558        for i,El in enumerate(Tdata):
    574559            FF[i] = refl[-1][El]           
    … …  
    611596        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
    612597        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
    613     maxPos = len(SGData['SGOps'])       
    614598    Uij = np.array(G2lat.U6toUij(Uijdata))
    615599    bij = Mast*Uij.T
Note: See TracChangeset for help on using the changeset viewer.