Changeset 1668


Ignore:
Timestamp:
Feb 23, 2015 4:52:38 PM (7 years ago)
Author:
vondreele
Message:

Attempt at a molecule builder in Atoms - seems to work in some cases

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIphsGUI.py

    r1663 r1668  
    14441444        # exclList will be 'x' or 'xu' if TLS used in RB
    14451445        Items = [G2gd.wxID_ATOMSEDITINSERT,G2gd.wxID_ATOMSEDITDELETE,G2gd.wxID_ATOMSREFINE,
    1446             G2gd.wxID_ATOMSMODIFY,G2gd.wxID_ATOMSTRANSFORM,G2gd.wxID_ATOMVIEWINSERT,G2gd.wxID_ATOMMOVE]
     1446            G2gd.wxID_ATOMSMODIFY,G2gd.wxID_ATOMSTRANSFORM,G2gd.wxID_MAKEMOLECULE,
     1447            G2gd.wxID_ATOMVIEWINSERT,G2gd.wxID_ATOMMOVE]
    14471448        if atomData:
    14481449            for item in Items:   
     
    17961797            else:
    17971798                Atoms.ForceRefresh()
     1799               
     1800    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       
     1898        indx = Atoms.GetSelectedRows()
     1899        Oxyz = []
     1900        xyz = []
     1901        DisAglCtls = {}
     1902        if len(indx) == 1:
     1903            generalData = data['General']
     1904            if 'DisAglCtls' in generalData:
     1905                DisAglCtls = generalData['DisAglCtls']
     1906            dlg = G2gd.DisAglDialog(G2frame,DisAglCtls,generalData)
     1907            if dlg.ShowModal() == wx.ID_OK:
     1908                DisAglCtls = dlg.GetData()
     1909            else:
     1910                dlg.Destroy()
     1911                return
     1912            dlg.Destroy()
     1913            generalData['DisAglCtls'] = DisAglCtls
     1914            atomData = copy.deepcopy(data['Atoms'])
     1915            data['Atoms'] = FindMolecule(indx[0])
     1916            FillAtomsGrid(Atoms)
     1917           
     1918           
     1919#                G2frame.ErrorDialog('Distance/Angle calculation','try again but do "Reset" to fill in missing atom types')
     1920        else:
     1921            print "select one atom"
     1922            G2frame.ErrorDialog('Select one atom',"select one atom to begin molecule build then redo")
    17981923
    17991924    def OnDistAnglePrt(event):
     
    58645989        G2frame.dataFrame.Bind(wx.EVT_MENU, AtomModify, id=G2gd.wxID_ATOMSMODIFY)
    58655990        G2frame.dataFrame.Bind(wx.EVT_MENU, AtomTransform, id=G2gd.wxID_ATOMSTRANSFORM)
     5991        G2frame.dataFrame.Bind(wx.EVT_MENU, MakeMolecule, id=G2gd.wxID_MAKEMOLECULE)
    58665992        G2frame.dataFrame.Bind(wx.EVT_MENU, OnReloadDrawAtoms, id=G2gd.wxID_RELOADDRAWATOMS)
    58675993        G2frame.dataFrame.Bind(wx.EVT_MENU, OnDistAngle, id=G2gd.wxID_ATOMSDISAGL)
Note: See TracChangeset for help on using the changeset viewer.