Changeset 1592


Ignore:
Timestamp:
Dec 2, 2014 2:16:00 PM (9 years ago)
Author:
vondreele
Message:

finish modulation vector determination - use so.brute - slow but works
add peak indexing option to not scale M20 by 1+X20 - might be useful

Location:
trunk
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIindex.py

    r1587 r1592  
    2929import GSASIIlattice as G2lat
    3030import GSASIIpwd as G2pwd
     31import GSASIIspc as G2spc
     32import GSASIImath as G2mth
    3133import scipy.optimize as so
    3234
    … …  
    171173    return [a,b,c,alp,bet,gam]
    172174   
    173 def calc_M20(peaks,HKL):
     175def calc_M20(peaks,HKL,ifX20=True):
    174176    'needs a doc string'
    175177    diff = 0
    … …  
    197199    else:
    198200        M20 = 0
    199     M20 /= (1.+X20)
     201    if ifX20:
     202        M20 /= (1.+X20)
    200203    return M20,X20
    201204   
    … …  
    257260    return X
    258261   
    259 def findMV(peaks,HKL,ssopt):
    260 #    import basinhopping as bh
    261     print ssopt
    262     return ssopt['ModVec']
     262def findMV(peaks,controls,ssopt,Inst):
     263       
     264    def Val2Vec(vec,Vref,values):
     265        Vec = []
     266        i = 0
     267        for j,r in enumerate(Vref):
     268            if r:
     269                if values.size > 1:
     270                    Vec.append(values[i])
     271                else:
     272                    Vec.append(values)                   
     273                i += 1
     274            else:
     275                Vec.append(vec[j])
     276        return np.array(Vec) 
     277     
     278    def ZSSfunc(values,peaks,dmin,Inst,SGData,SSGData,vec,Vref,maxH,A,wave,Z):
     279        Vec = Val2Vec(vec,Vref,values)
     280        HKL =  G2pwd.getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A)
     281        Peaks = np.array(IndexSSPeaks(peaks,HKL)[1]).T
     282        Qo = 1./Peaks[-2]**2
     283        Qc = G2lat.calc_rDsqZSS(Peaks[4:8],A,Vec,Z,Peaks[0],wave)
     284        return np.sum((Qo-Qc)**2)
     285
     286    if 'C' in Inst['Type'][0]:
     287        wave = G2mth.getWave(Inst)
     288    else:
     289        difC = Inst['difC'][1]
     290    SGData = G2spc.SpcGroup(controls[13])[1]
     291    SSGData = G2spc.SSpcGroup(SGData,ssopt['ssSymb'])
     292    A = G2lat.cell2A(controls[6:12])
     293    Z = controls[1]
     294    Vref = [True if x in ssopt['ssSymb'] else False for x in ['a','b','g']]
     295    values = []
     296    ranges = []   
     297    for v,r in zip(ssopt['ModVec'],Vref):
     298        if r:
     299            ranges += [slice(.02,.98,.05),]
     300            values += [v,]
     301    dmin = getDmin(peaks)-0.005
     302    Peaks = np.copy(np.array(peaks).T)   
     303    result = so.brute(ZSSfunc,ranges,finish=so.fmin_cg,
     304        args=(peaks,dmin,Inst,SGData,SSGData,ssopt['ModVec'],Vref,ssopt['maxH'],A,wave,Z))
     305    return Val2Vec(ssopt['ModVec'],Vref,result)
    263306               
    264307def IndexPeaks(peaks,HKL):
    … …  
    311354    import bisect
    312355    N = len(HKL)
    313     if N == 0: return False,peaks
     356    Peaks = np.copy(peaks)
     357    if N == 0: return False,Peaks
    314358    if len(peaks[0]) == 9:      #add m column if missing
    315         for peak in peaks:
     359        for peak in Peaks:
    316360            peak.insert(7,0)
    317361    hklds = list(np.array(HKL).T[4])+[1000.0,0.0,]
    318362    hklds.sort()                                        # ascending sort - upper bound at end
    319363    hklmax = [0,0,0,0]
    320     for ipk,peak in enumerate(peaks):
     364    for ipk,peak in enumerate(Peaks):
    321365        if peak[2]: #Use
    322366            peak[4:8] = [0,0,0,0]                           #clear old indexing
    … …  
    332376            hkl = HKL[pos]                                 # put in hkl
    333377            if hkl[-1] >= 0:                                 # peak already assigned - test if this one better
    334                 opeak = peaks[hkl[-1]]
     378                opeak = Peaks[hkl[-1]]
    335379                dold = abs(opeak[-2]-hkl[4])
    336380                dnew = min(dm,dp)
    … …  
    343387            peak[4:8] = hkl[:4]
    344388            peak[9] = hkl[4]                                # fill in d-calc
    345     for peak in peaks:
     389    for peak in Peaks:
    346390        peak[3] = False
    347391        if peak[2]:
    … …  
    351395                peak[3] = True
    352396    if hklmax[0]*hklmax[1]*hklmax[2]*hklmax[3] > 0:
    353         return True,peaks
    354     else:
    355         return False,peaks  #nothing indexed!
     397        return True,Peaks
     398    else:
     399        return False,Peaks  #nothing indexed!
    356400       
    357401def Values2A(ibrav,values):
    … …  
    654698    return len(HKL),M20,X20,Aref,Z
    655699   
    656 def refinePeaks(peaks,ibrav,A):
     700def refinePeaks(peaks,ibrav,A,ifX20=True):
    657701    'needs a doc string'
    658702    dmin = getDmin(peaks)
    … …  
    696740            A = oldA
    697741       
    698     M20,X20 = calc_M20(peaks,HKL)
     742    M20,X20 = calc_M20(peaks,HKL,ifX20)
    699743    return len(HKL),M20,X20,A
    700744       
    701 def findBestCell(dlg,ncMax,A,Ntries,ibrav,peaks,V1):
     745def findBestCell(dlg,ncMax,A,Ntries,ibrav,peaks,V1,ifX20=True):
    702746    'needs a doc string'
    703747# dlg & ncMax are used for wx progress bar
    … …  
    715759        if len(HKL) > mHKL[ibrav]:
    716760            peaks = IndexPeaks(peaks,HKL)[1]
    717             Asave.append([calc_M20(peaks,HKL),A[:]])
     761            Asave.append([calc_M20(peaks,HKL,ifX20),A[:]])
    718762    tries = 0
    719763    while tries < Ntries:
    … …  
    727771       
    728772        if IndexPeaks(peaks,HKL)[0] and len(HKL) > mHKL[ibrav]:
    729             Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg)
    730             Asave.append([calc_M20(peaks,HKL),Aref[:]])
     773            Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg,ifX20)
     774            Asave.append([calc_M20(peaks,HKL,ifX20),Aref[:]])
    731775            if ibrav == 9:                          #C-centered orthorhombic
    732776                for i in range(2):
    733777                    Abeg = rotOrthoA(Abeg[:])
    734                     Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg)
     778                    Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg,ifX20)
    735779                    HKL = G2lat.GenHBravais(dmin,ibrav,Aref)
    736780                    peaks = IndexPeaks(peaks,HKL)[1]
    737                     Asave.append([calc_M20(peaks,HKL),Aref[:]])
     781                    Asave.append([calc_M20(peaks,HKL,ifX20),Aref[:]])
    738782            elif ibrav == 11:                      #C-centered monoclinic
    739783                Abeg = swapMonoA(Abeg[:])
    740                 Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg)
     784                Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg,ifX20)
    741785                HKL = G2lat.GenHBravais(dmin,ibrav,Aref)
    742786                peaks = IndexPeaks(peaks,HKL)[1]
    743                 Asave.append([calc_M20(peaks,HKL),Aref[:]])
     787                Asave.append([calc_M20(peaks,HKL,ifX20),Aref[:]])
    744788        else:
    745789            break
    … …  
    754798    X = sortM20(Asave)
    755799    if X:
    756         Lhkl,M20,X20,A = refinePeaks(peaks,ibrav,X[0][1])
     800        Lhkl,M20,X20,A = refinePeaks(peaks,ibrav,X[0][1],ifX20)
    757801        return GoOn,Lhkl,M20,X20,A
    758802       
    … …  
    782826    return A
    783827
    784 def DoIndexPeaks(peaks,controls,bravais):
     828def DoIndexPeaks(peaks,controls,bravais,ifX20=True):
    785829    'needs a doc string'
    786830   
    … …  
    825869                                if not N2:
    826870                                    A = []
    827                                     GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
     871                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1,ifX20)
    828872                                if A:
    829                                     GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A[:],N1s[ibrav],ibrav,peaks,0)
     873                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A[:],N1s[ibrav],ibrav,peaks,0,ifX20)
    830874                            else:
    831                                 GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
     875                                GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1,ifX20)
    832876                            if Nc >= ncMax:
    833877                                GoOn = False
    … …  
    848892                                    a,b,c,alp,bet,gam = G2lat.A2cell(A)
    849893                                    V = G2lat.calc_V(A)
    850                                     print "%10.3f %3d %3d %10.5f %10.5f %10.5f %10.3f %10.3f %10.3f %10.2f %10.2f" % (M20,X20,Nc,a,b,c,alp,bet,gam,V,V1)
    851894                                    if M20 >= 2.0:
     895                                        print "%10.3f %3d %3d %10.5f %10.5f %10.5f %10.3f %10.3f %10.3f %10.2f %10.2f" % (M20,X20,Nc,a,b,c,alp,bet,gam,V,V1)
    852896                                        cells.append([M20,X20,ibrav,a,b,c,alp,bet,gam,V,False,False])
    853897                            if not GoOn:

  • trunk/GSASIIplot.py

    r1586 r1592  
    10711071                    pickIdText = G2frame.PatternTree.GetItemText(G2frame.PickId)
    10721072                    if pickIdText in ['Index Peak List','Unit Cells List','Reflection Lists'] or \
    1073                         'PWDR' in G2frame.PatternTree.GetItemText(PickId):
     1073                        'PWDR' in pickIdText:
    10741074                        indx = -1
    10751075                        if pickIdText in ['Index Peak List','Unit Cells List',]:

  • trunk/GSASIIpwdGUI.py

    r1587 r1592  
    21752175    UnitCellsId = G2gd.GetPatternTreeItemId(G2frame,G2frame.PatternId, 'Unit Cells List')
    21762176    SPGlist = G2spc.spglist
     2177    ifX20 = True
    21772178    bravaisSymb = ['Fm3m','Im3m','Pm3m','R3-H','P6/mmm','I4/mmm',
    21782179        'P4/mmm','Fmmm','Immm','Cmmm','Pmmm','C2/m','P2/m','P1']
    … …  
    22062207        controls[2] = NcNo.GetValue()
    22072208       
     2209    def OnIfX20(event):
     2210        ifX20 = x20.GetValue()
     2211       
    22082212    def OnStartVol(event):
    22092213        try:
    … …  
    22752279       
    22762280    def OnFindMV(event):
    2277         ssopt['ModVec'] = G2indx.findMV(peaks,G2frame.HKL,ssopt)
     2281        Peaks = np.copy(peaks[0])
     2282        ssopt['ModVec'] = G2indx.findMV(Peaks,controls,ssopt,Inst)
    22782283        OnHklShow(event)
    22792284        wx.CallAfter(UpdateUnitCellsGrid,G2frame,data)
    … …  
    23882393        if colLabels[c] == 'M20':
    23892394            cells = G2indx.sortM20(cells)
    2390         elif colLabels[c] in ['Bravais','a','b','c','alpha','beta','gamma','Volume']:
     2395        elif colLabels[c] in ['X20','Bravais','a','b','c','alpha','beta','gamma','Volume']:
     2396            if c == 1:
     2397                c += 1  #X20 before Use
    23912398            cells = G2indx.sortCells(cells,c-1)     #an extra column (Use) not in cells
    23922399        else:
    … …  
    24342441        PickId = G2frame.PickId   
    24352442        peaks = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId, 'Index Peak List'))
    2436         if not peaks[0]:
     2443        if not len(peaks[0]):
    24372444            G2frame.ErrorDialog('No peaks!', 'Nothing to refine!')
    24382445            return       
    … …  
    24972504            for cell in cells:
    24982505                if cell[11]:
     2506                    cell[10] = False    #clear selection flag on keepers
    24992507                    keepcells.append(cell)
    25002508        except IndexError:
    … …  
    25022510        except ValueError:
    25032511            G2frame.ErrorDialog('Error','Need to set controls in Unit Cell List first')
     2512            return
     2513        if ssopt.get('Use',False):
     2514            G2frame.ErrorDialog('Super lattice error','Indexing not available for super lattices')
    25042515            return
    25052516        if True not in bravais:
    … …  
    25092520            G2frame.ErrorDialog('Error','Index Peak List is empty')
    25102521            return
     2522        if len(peaks[0][0]) > 9:
     2523            G2frame.ErrorDialog('Error','You need to reload Index Peaks List first')
     2524            return
    25112525        G2frame.dataFrame.CopyCell.Enable(False)
    25122526        G2frame.dataFrame.RefineCell.Enable(False)
    2513         OK,dmin,newcells = G2indx.DoIndexPeaks(peaks[0],controls,bravais)
     2527        OK,dmin,newcells = G2indx.DoIndexPeaks(peaks[0],controls,bravais,ifX20)
    25142528        cells = keepcells+newcells
    25152529        cells = G2indx.sortM20(cells)
    2516         cells[0][10] = True
     2530        cells[0][10] = True         #select best M20
    25172531        if OK:
    25182532            data = [controls,bravais,cells,dmin,ssopt]
    … …  
    25332547               
    25342548    def RefreshUnitCellsGrid(event):
    2535         data =G2frame.PatternTree.GetItemPyData(UnitCellsId)
     2549        data = G2frame.PatternTree.GetItemPyData(UnitCellsId)
    25362550        cells,dmin = data[2:4]
    25372551        r,c =  event.GetRow(),event.GetCol()
    … …  
    25422556                    UnitCellsTable.SetValue(i,c,False)
    25432557                UnitCellsTable.SetValue(r,c,True)
    2544                 gridDisplay.ForceRefresh()
     2558                gridDisplay.Refresh()
    25452559                cells[r][-2] = True
    25462560                ibrav = cells[r][2]
    … …  
    25612575                    UnitCellsTable.SetValue(r,c,True)
    25622576                gridDisplay.ForceRefresh()
    2563             G2frame.PatternTree.SetItemPyData(UnitCellsId,data)               
     2577            G2frame.PatternTree.SetItemPyData(UnitCellsId,data)
    25642578       
    25652579    def MakeNewPhase(event):
    … …  
    25862600        finally:
    25872601            dlg.Destroy()
    2588            
    25892602    if G2frame.dataDisplay:
     2603        print G2frame.dataDisplay.GetScrollPos(wx.VERTICAL)
    25902604        G2frame.dataFrame.DestroyChildren()
    25912605    G2frame.dataDisplay = wxscroll.ScrolledPanel(G2frame.dataFrame)
    … …  
    26452659    startVol.Bind(wx.EVT_KILL_FOCUS,OnStartVol)
    26462660    littleSizer.Add(startVol,0,WACV)
     2661    x20 = wx.CheckBox(G2frame.dataDisplay,label='Use M20/(X20+1)?')
     2662    x20.SetValue(ifX20)
     2663    x20.Bind(wx.EVT_CHECKBOX,OnIfX20)
     2664    littleSizer.Add(x20,0,WACV)
    26472665    mainSizer.Add(littleSizer,0)
    26482666    mainSizer.Add((5,5),0)
    … …  
    27942812        UnitCellsTable = G2gd.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
    27952813        gridDisplay = G2gd.GSGrid(G2frame.dataDisplay)
    2796         gridDisplay.SetPosition(wx.Point(0,20))               
    27972814        gridDisplay.SetTable(UnitCellsTable, True)
    27982815        G2frame.dataFrame.CopyCell.Enable(True)
Note: See TracChangeset for help on using the changeset viewer.