Changeset 157


Ignore:
Timestamp:
Sep 27, 2010 3:05:00 PM (12 years ago)
Author:
vondreele
Message:

much faster bond finding

Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIphsGUI.py

    r155 r157  
    11481148            G2plt.PlotStructure(self,data)
    11491149           
    1150     def FindBonds6():
     1150    def FindBondsToo():                         #works but slow for large structures - keep as reference
     1151        cx,ct,cs = data['Drawing']['atomPtrs']
     1152        atomData = data['Drawing']['Atoms']
     1153        generalData = data['General']
     1154        Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
     1155        radii = generalData['BondRadii']
     1156        atomTypes = generalData['AtomTypes']
     1157        try:
     1158            indH = atomTypes.index('H')
     1159            radii[indH] = 0.5
     1160        except:
     1161            pass           
     1162        for atom in atomData:
     1163            atom[-1] = []
     1164        Atoms = []
     1165        for i,atom in enumerate(atomData):
     1166            Atoms.append([i,np.array(atom[cx:cx+3]),atom[cs],radii[atomTypes.index(atom[ct])]])
     1167        for atomA in Atoms:
     1168            if atomA[2] in ['lines','sticks','ellipsoids','balls & sticks','polyhedra']:
     1169                for atomB in Atoms:                   
     1170                    Dx = atomB[1]-atomA[1]
     1171                    DX = np.inner(Amat,Dx)
     1172                    dist = np.sqrt(np.sum(DX**2))
     1173                    sumR = atomA[3]+atomB[3]
     1174                    if 0.5 < dist <= 0.85*sumR:
     1175                        i = atomA[0]
     1176                        if atomA[2] == 'polyhedra':
     1177                            atomData[i][-1].append(DX)
     1178                        elif atomB[1] != 'polyhedra':
     1179                            j = atomB[0]
     1180                            atomData[i][-1].append(Dx*atomA[3]/sumR)
     1181                            atomData[j][-1].append(-Dx*atomB[3]/sumR)
     1182                   
     1183    def FindBonds():                    #uses numpy & masks - very fast even for proteins!
     1184        import numpy.ma as ma
    11511185        cx,ct,cs = data['Drawing']['atomPtrs']
    11521186        atomData = data['Drawing']['Atoms']
     
    11621196        for atom in atomData:
    11631197            atom[-1] = []               #clear out old bonds
    1164        
    1165            
    1166     def FindBonds():
    1167         cx,ct,cs = data['Drawing']['atomPtrs']
    1168         atomData = data['Drawing']['Atoms']
    1169         generalData = data['General']
    1170         Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
    1171         radii = generalData['BondRadii']
    1172         atomTypes = generalData['AtomTypes']
    1173         try:
    1174             indH = atomTypes.index('H')
    1175             radii[indH] = 0.5
    1176         except:
    1177             pass           
     1198        Indx = range(len(atomData))
     1199        Atoms = []
     1200        Types = []
     1201        Radii = []
    11781202        for atom in atomData:
    1179             atom[-1] = []               #clear out old bonds
    1180         for i,atomA in enumerate(atomData):
    1181             if atomA[cs] in ['lines','sticks','ellipsoids','balls & sticks','polyhedra']:
    1182                 xyzA = np.array(atomA[cx:cx+3])
    1183                 indA = atomTypes.index(atomA[ct])
    1184                 for j,atomB in enumerate(atomData):
    1185                     xyzB = np.array(atomB[cx:cx+3])
    1186                     indB = atomTypes.index(atomB[ct])
    1187                     Dx = xyzB-xyzA
    1188                     DX = np.inner(Amat,Dx)
    1189                     dist = np.sqrt(np.sum(DX**2))
    1190                     sumR = radii[indA]+radii[indB]
    1191                     if 0 < dist <= 0.85*sumR:
    1192                         inc = i+1
    1193                         if atomA[cs] == 'polyhedra':
    1194                             atomA[-1].append(np.inner(Amat,Dx))
    1195                         elif atomB[cs] != 'polyhedra':
    1196                             atomA[-1].append(Dx*radii[indA]/sumR)
    1197                             atomB[-1].append(-Dx*radii[indB]/sumR)
     1203            Atoms.append(np.array(atom[cx:cx+3]))
     1204            Types.append(atom[cs])
     1205            Radii.append(radii[atomTypes.index(atom[ct])])
     1206        Atoms = np.array(Atoms)
     1207        Radii = np.array(Radii)
     1208        IATR = zip(Indx,Atoms,Types,Radii)
     1209        for atomA in IATR:
     1210            if atomA[2] in ['lines','sticks','ellipsoids','balls & sticks','polyhedra']:
     1211                Dx = Atoms-atomA[1]
     1212                dist = ma.masked_less(np.sqrt(np.sum(np.inner(Amat,Dx)**2,axis=0)),0.5) #gets rid of self & disorder "bonds" < 0.5A
     1213                sumR = atomA[3]+Radii
     1214                IndB = ma.nonzero(ma.masked_greater(dist-0.85*sumR,0.))                 #get indices of bonded atoms
     1215                i = atomA[0]
     1216                for j in IndB[0]:
     1217                    if Types[i] == 'polyhedra':
     1218                        atomData[i][-1].append(np.inner(Amat,Dx[j]))
     1219                    elif Types[j] != 'polyhedra':
     1220                        atomData[i][-1].append(Dx[j]*Radii[i]/sumR[j])
     1221                        atomData[j][-1].append(-Dx[j]*Radii[j]/sumR[j])
    11981222
    11991223    def DrawAtomsDelete(event):   
  • trunk/GSASIIplot.py

    r155 r157  
    14221422            RenderUnitVectors()
    14231423            glPopMatrix()
     1424        time0 = time.time()
    14241425        for iat,atom in enumerate(drawingData['Atoms']):
    14251426            x,y,z = atom[cx:cx+3]
     
    14631464                    RenderSphere(x,y,z,radius,color)
    14641465            elif 'lines' in atom[cs]:
    1465                 radius = 0.5
     1466                radius = 0.1
    14661467                RenderLines(x,y,z,Bonds,color)
    14671468            elif atom[cs] == 'sticks':
    1468                 radius = 0.5
     1469                radius = 0.1
    14691470                RenderBonds(x,y,z,Bonds,bondR,color)
    14701471            elif atom[cs] == 'polyhedra':
     
    14931494            elif atom[cs+1] == 'chain' and atom[ct-1] == 'CA':
    14941495                RenderLabel(x,y,z,atom[ct-2],radius)
    1495            
     1496        print time.time()-time0   
    14961497        Page.canvas.SwapBuffers()
    14971498       
Note: See TracChangeset for help on using the changeset viewer.