Changeset 1669 for trunk


Ignore:
Timestamp:
Feb 24, 2015 11:13:07 AM (10 years ago)
Author:
vondreele
Message:

trap IndexErrors? in GetTypeName? & sortArray
add FindMolecule? to G2math (taken out of G2phsGUI)
SearchMap? now adds dcent (peak-cent distance) to table of results
tighten tolerance to 0.5A in noDuplicate (peaks)
Add dcent column to map peaks GUI - can be sorted

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • TabularUnified trunk/GSASIIgrid.py

    r1667 r1669  
    24332433            if self.data[row][col] is None: return None
    24342434            return self.dataTypes[col]
    2435         except TypeError:
     2435        except (TypeError,IndexError):
    24362436            return None
    24372437

  • TabularUnified trunk/GSASIImath.py

    r1659 r1669  
    205205################################################################################
    206206
     207def FindMolecule(ind,generalData,atomData):                    #uses numpy & masks - very fast even for proteins!
     208
     209    def getNeighbors(atom,radius):
     210        neighList = [] 
     211        Dx = IARS[1]-np.array(atom[cx:cx+3])
     212        dist = ma.masked_less(np.sqrt(np.sum(np.inner(Amat,Dx)**2,axis=0)),0.5) #gets rid of disorder "bonds" < 0.5A
     213        sumR = IARS[2]+radius
     214        return set(ma.nonzero(ma.masked_greater(dist-factor*sumR,0.))[0])                #get indices of bonded atoms
     215
     216    import numpy.ma as ma
     217    indices = (-1,0,1)
     218    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices],dtype='f')
     219    cx,ct,cs,ci = generalData['AtomPtrs']
     220    DisAglCtls = generalData['DisAglCtls']
     221    SGData = generalData['SGData']
     222    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
     223    radii = DisAglCtls['BondRadii']
     224    atomTypes = DisAglCtls['AtomTypes']
     225    factor = DisAglCtls['Factors'][0]
     226    unit = np.zeros(3)
     227    try:
     228        indH = atomTypes.index('H')
     229        radii[indH] = 0.5
     230    except:
     231        pass           
     232    nAtom = len(atomData)
     233    Indx = range(nAtom)
     234    UAtoms = []
     235    SymOp = []
     236    Radii = []
     237    for atom in atomData:
     238        UAtoms.append(np.array(atom[cx:cx+3]))
     239        Radii.append(radii[atomTypes.index(atom[ct])])
     240        SymOp += [[1,0,unit],]
     241    UAtoms = np.array(UAtoms)
     242    Radii = np.array(Radii)
     243    for nOp,Op in enumerate(SGData['SGOps'][1:]):
     244        UAtoms = np.concatenate((UAtoms,(np.inner(Op[0],UAtoms[:nAtom]).T+Op[1])))
     245        Radii = np.concatenate((Radii,Radii[:nAtom]))
     246        SymOp += [[nOp,0,unit] for symop in SymOp]
     247        Indx += Indx[:nAtom]
     248    for icen,cen in enumerate(SGData['SGCen'][1:]):
     249        UAtoms = np.concatenate((UAtoms,(UAtoms+cen)))
     250        Radii = np.concatenate((Radii,Radii))
     251        SymOp += [[symop[0],icen+1,unit] for symop in SymOp]
     252        Indx += Indx[:nAtom]
     253    if SGData['SGInv']:
     254        UAtoms = np.concatenate((UAtoms,-UAtoms))
     255        Radii = np.concatenate((Radii,Radii))
     256        SymOp += [[-symop[0],symop[1],unit] for symop in SymOp]
     257        Indx += Indx
     258    UAtoms %= 1.
     259    mAtoms = len(UAtoms)
     260    for unit in Units:
     261        if np.any(unit):    #skip origin cell
     262            UAtoms = np.concatenate((UAtoms,UAtoms[:mAtoms]+unit))
     263            Radii = np.concatenate((Radii,Radii[:mAtoms]))
     264            SymOp += [[symop[0],symop[1],unit] for symop in SymOp[:mAtoms]]                       
     265            Indx += Indx[:mAtoms]
     266    UAtoms = np.array(UAtoms)
     267    Radii = np.array(Radii)
     268    IARS = [Indx,UAtoms,Radii,SymOp]
     269    newAtoms = [atomData[ind],]
     270    atomData[ind] = None
     271    radius = Radii[ind]
     272    IndB = getNeighbors(newAtoms[-1],radius)
     273    while True:
     274        if not len(IndB):
     275            break
     276        indb = IndB.pop()
     277        if atomData[Indx[indb]] == None:
     278            continue
     279        while True:
     280            try:
     281                jndb = Indb.index(indb)
     282                Indb.remove(jndb)
     283            except:
     284                break
     285        newAtom = copy.copy(atomData[Indx[indb]])
     286        newAtom[cx:cx+3] = UAtoms[indb]     #NB: thermal Uij, etc. not transformed!
     287        newAtoms.append(newAtom)
     288        atomData[Indx[indb]] = None
     289        IndB = set(list(IndB)+list(getNeighbors(newAtoms[-1],radius)))
     290    for atom in atomData:
     291        if atom != None:
     292            newAtoms.append(atom)
     293    return newAtoms
     294       
    207295def FindAtomIndexByIDs(atomData,IDs,Draw=True):
    208296    '''finds the set of atom array indices for a list of atom IDs. Will search
    … …  
    21032191        * dzeros : ndarray
    21042192          the distance of the peaks from  the unit cell origin
     2193        * dcent : ndarray
     2194          the distance of the peaks from  the unit cell center
    21052195
    21062196    '''       
    … …  
    21632253    mags = []
    21642254    dzeros = []
     2255    dcent = []
    21652256    try:
    21662257        mapData = generalData['Map']
    … …  
    22022293            peak = (np.array(x1[1:4])-ind)/incre
    22032294        peak = fixSpecialPos(peak,SGData,Amat)
    2204         rho = rollMap(rho,-ind)       
     2295        rho = rollMap(rho,-ind)
     2296    cent = np.ones(3)*.5     
    22052297    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
     2298    dcent = np.sqrt(np.sum(np.inner(Amat,peaks-cent)**2,axis=0))
    22062299    if Neg:     #want negative magnitudes for negative peaks
    2207         return np.array(peaks),-np.array([mags,]).T,np.array([dzeros,]).T
     2300        return np.array(peaks),-np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
    22082301    else:
    2209         return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T
     2302        return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
    22102303   
    22112304def sortArray(data,pos,reverse=False):
    … …  
    22152308    T = []
    22162309    for i,M in enumerate(data):
    2217         T.append((M[pos],i))
     2310        try:
     2311            T.append((M[pos],i))
     2312        except IndexError:
     2313            return data
    22182314    D = dict(zip(T,data))
    22192315    T.sort()
    … …  
    22712367
    22722368    def noDuplicate(xyz,peaks,Amat):
    2273         if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=1.0) for peak in peaks]:
     2369        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
    22742370            return False
    22752371        return True

  • TabularUnified trunk/GSASIIphsGUI.py

    r1668 r1669  
    17991799               
    18001800    def MakeMolecule(event):     
    1801        
    1802         def FindMolecule(ind):                    #uses numpy & masks - very fast even for proteins!
    1803 
    1804             def getNeighbors(atom,radius):
    1805                 neighList = [] 
    1806                 Dx = IARS[1]-np.array(atom[cx:cx+3])
    1807                 dist = ma.masked_less(np.sqrt(np.sum(np.inner(Amat,Dx)**2,axis=0)),0.5) #gets rid of disorder "bonds" < 0.5A
    1808                 sumR = IARS[2]+radius
    1809                 return set(ma.nonzero(ma.masked_greater(dist-factor*sumR,0.))[0])                #get indices of bonded atoms
    1810        
    1811             import numpy.ma as ma
    1812             indices = (-1,0,1)
    1813             Units = np.array([[h,k,l] for h in indices for k in indices for l in indices],dtype='f')
    1814             cx,ct,cs,ci = generalData['AtomPtrs']
    1815             SGData = generalData['SGData']
    1816             Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
    1817             radii = DisAglCtls['BondRadii']
    1818             atomTypes = DisAglCtls['AtomTypes']
    1819             factor = DisAglCtls['Factors'][0]
    1820             unit = np.zeros(3)
    1821             try:
    1822                 indH = atomTypes.index('H')
    1823                 radii[indH] = 0.5
    1824             except:
    1825                 pass           
    1826             nAtom = len(atomData)
    1827             Indx = range(nAtom)
    1828             UAtoms = []
    1829             SymOp = []
    1830             Radii = []
    1831             for atom in atomData:
    1832                 UAtoms.append(np.array(atom[cx:cx+3]))
    1833                 Radii.append(radii[atomTypes.index(atom[ct])])
    1834                 SymOp += [[1,0,unit],]
    1835             UAtoms = np.array(UAtoms)
    1836             Radii = np.array(Radii)
    1837             for nOp,Op in enumerate(SGData['SGOps'][1:]):
    1838                 UAtoms = np.concatenate((UAtoms,(np.inner(Op[0],UAtoms[:nAtom]).T+Op[1])))
    1839                 Radii = np.concatenate((Radii,Radii[:nAtom]))
    1840                 SymOp += [[nOp,0,unit] for symop in SymOp]
    1841                 Indx += Indx[:nAtom]
    1842             for icen,cen in enumerate(SGData['SGCen'][1:]):
    1843                 UAtoms = np.concatenate((UAtoms,(UAtoms+cen)))
    1844                 Radii = np.concatenate((Radii,Radii))
    1845                 SymOp += [[symop[0],icen+1,unit] for symop in SymOp]
    1846                 Indx += Indx[:nAtom]
    1847             if SGData['SGInv']:
    1848                 UAtoms = np.concatenate((UAtoms,-UAtoms))
    1849                 Radii = np.concatenate((Radii,Radii))
    1850                 SymOp += [[-symop[0],symop[1],unit] for symop in SymOp]
    1851                 Indx += Indx
    1852             UAtoms %= 1.
    1853             mAtoms = len(UAtoms)
    1854             for unit in Units:
    1855                 if np.any(unit):    #skip origin cell
    1856                     UAtoms = np.concatenate((UAtoms,UAtoms[:mAtoms]+unit))
    1857                     Radii = np.concatenate((Radii,Radii[:mAtoms]))
    1858                     SymOp += [[symop[0],symop[1],unit] for symop in SymOp[:mAtoms]]                       
    1859                     Indx += Indx[:mAtoms]
    1860             UAtoms = np.array(UAtoms)
    1861             Radii = np.array(Radii)
    1862             IARS = [Indx,UAtoms,Radii,SymOp]
    1863             newAtoms = [atomData[ind],]
    1864             atomData[ind] = None
    1865             radius = Radii[ind]
    1866             IndB = getNeighbors(newAtoms[-1],radius)
    1867             while True:
    1868                 if not len(IndB):
    1869                     break
    1870                 indb = IndB.pop()
    1871                 if atomData[Indx[indb]] == None:
    1872                     continue
    1873                 while True:
    1874                     try:
    1875                         jndb = Indb.index(indb)
    1876                         Indb.remove(jndb)
    1877                     except:
    1878                         break
    1879                 newAtom = copy.copy(atomData[Indx[indb]])
    1880 #                ops = SymOp[indb]
    1881 #                sop = SGData['SGOps'][abs(ops[0])-1]
    1882 #                xyz = np.array(newAtom[cx:cx+3])
    1883 #                xyz = np.inner(sop[0],xyz)+sop[1]
    1884 #                xyz += SGData['SGCen'][ops[1]]
    1885 #                if ops[0] < 0:
    1886 #                    xyz *= -1
    1887 #                xyz %= 1.
    1888 #                xyz += ops[2]
    1889                 newAtom[cx:cx+3] = UAtoms[indb]
    1890                 newAtoms.append(newAtom)
    1891                 atomData[Indx[indb]] = None
    1892                 IndB = set(list(IndB)+list(getNeighbors(newAtoms[-1],radius)))
    1893             for atom in atomData:
    1894                 if atom != None:
    1895                     newAtoms.append(atom)
    1896             return newAtoms
    1897        
    18981801        indx = Atoms.GetSelectedRows()
    18991802        Oxyz = []
    … …  
    19131816            generalData['DisAglCtls'] = DisAglCtls
    19141817            atomData = copy.deepcopy(data['Atoms'])
    1915             data['Atoms'] = FindMolecule(indx[0])
     1818            data['Atoms'] = G2mth.FindMolecule(indx[0],generalData,atomData)
     1819            OnReloadDrawAtoms(event)           
    19161820            FillAtomsGrid(Atoms)
    19171821           
    … …  
    55575461                mapPeaks = data['Map Peaks']
    55585462                c =  event.GetCol()
    5559                 if colLabels[c] == 'mag':
     5463                if colLabels[c] == 'mag':   #big to small order
    55605464                    mapPeaks = G2mth.sortArray(mapPeaks,c,reverse=True)
    5561                 elif colLabels[c] in ['x','y','z','dzero']:
     5465                elif colLabels[c] in ['x','y','z','dzero','dcent']:     #small to big
    55625466                    mapPeaks = G2mth.sortArray(mapPeaks,c)
    55635467                else:
    … …  
    55675471            G2plt.PlotStructure(G2frame,data)                   
    55685472           
    5569         G2frame.dataFrame.setSizePosLeft([450,300])
     5473        G2frame.dataFrame.setSizePosLeft([500,300])
    55705474        G2frame.dataFrame.SetStatusText('')
    55715475        if 'Map Peaks' in data:
    5572             G2frame.dataFrame.SetStatusText('Select mag or dzero columns to sort')
     5476            G2frame.dataFrame.SetStatusText('Double click any column heading to sort')
    55735477            mapPeaks = data['Map Peaks']                       
    55745478            rowLabels = []
    55755479            for i in range(len(mapPeaks)): rowLabels.append(str(i))
    5576             colLabels = ['mag','x','y','z','dzero']
    5577             Types = 5*[wg.GRID_VALUE_FLOAT+':10,4',]
     5480            colLabels = ['mag','x','y','z','dzero','dcent']
     5481            Types = 6*[wg.GRID_VALUE_FLOAT+':10,4',]
    55785482            G2frame.MapPeaksTable = G2gd.Table(mapPeaks,rowLabels=rowLabels,colLabels=colLabels,types=Types)
    55795483            MapPeaks.SetTable(G2frame.MapPeaksTable, True)
    … …  
    55915495            Ind = MapPeaks.GetSelectedRows()
    55925496            for ind in Ind:
    5593                 mag,x,y,z,d = mapPeaks[ind]
     5497                mag,x,y,z = mapPeaks[ind][:4]
    55945498                AtomAdd(x,y,z,'H',Name='M '+'%d'%(int(100*mag/peakMax)))
    55955499            G2plt.PlotStructure(G2frame,data)
    … …  
    57705674    def OnSearchMaps(event):
    57715675                                   
    5772         peaks = []
    5773         mags = []
    57745676        print ' Begin fourier map search - can take some time'
    57755677        time0 = time.time()
    … …  
    57805682            wx.BeginBusyCursor()
    57815683            try:
    5782                 peaks,mags,dzeros = G2mth.SearchMap(generalData,drawingData)
     5684                peaks,mags,dzeros,dcents = G2mth.SearchMap(generalData,drawingData)
    57835685                if 'N' in mapData['Type']:      #look for negatives in neutron maps
    5784                     npeaks,nmags,ndzeros = G2mth.SearchMap(generalData,drawingData,Neg=True)
     5686                    npeaks,nmags,ndzeros,ndcents = G2mth.SearchMap(generalData,drawingData,Neg=True)
    57855687                    peaks = np.concatenate((peaks,npeaks))
    57865688                    mags = np.concatenate((mags,nmags))
    57875689                    dzeros = np.concatenate((dzeros,ndzeros))
     5690                    dcents = np.concatenate((dcents,ndcents))
    57885691            finally:
    57895692                wx.EndBusyCursor()
    57905693            if len(peaks):
    5791                 mapPeaks = np.concatenate((mags,peaks,dzeros),axis=1)
     5694                mapPeaks = np.concatenate((mags,peaks,dzeros,dcents),axis=1)
    57925695                data['Map Peaks'] = G2mth.sortArray(mapPeaks,0,reverse=True)           
    57935696            print ' Map search finished, time = %.2fs'%(time.time()-time0)

  • TabularUnified trunk/GSASIIplot.py

    r1660 r1669  
    49524952            XYZ = mapPeaks.T[1:4].T
    49534953            mapBonds = FindPeaksBonds(XYZ)
    4954             for ind,[mag,x,y,z,d] in enumerate(mapPeaks):
     4954            for ind,[mag,x,y,z] in enumerate(mapPeaks[:,:4]):
    49554955                if ind in Ind and pageName == 'Map peaks':
    49564956                    RenderMapPeak(x,y,z,Gr,1.0)
Note: See TracChangeset for help on using the changeset viewer.