Changeset 1668
 Timestamp:
 Feb 23, 2015 4:52:38 PM (8 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIphsGUI.py
r1663 r1668 1444 1444 # exclList will be 'x' or 'xu' if TLS used in RB 1445 1445 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] 1447 1448 if atomData: 1448 1449 for item in Items: … … 1796 1797 else: 1797 1798 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(distfactor*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") 1798 1923 1799 1924 def OnDistAnglePrt(event): … … 5864 5989 G2frame.dataFrame.Bind(wx.EVT_MENU, AtomModify, id=G2gd.wxID_ATOMSMODIFY) 5865 5990 G2frame.dataFrame.Bind(wx.EVT_MENU, AtomTransform, id=G2gd.wxID_ATOMSTRANSFORM) 5991 G2frame.dataFrame.Bind(wx.EVT_MENU, MakeMolecule, id=G2gd.wxID_MAKEMOLECULE) 5866 5992 G2frame.dataFrame.Bind(wx.EVT_MENU, OnReloadDrawAtoms, id=G2gd.wxID_RELOADDRAWATOMS) 5867 5993 G2frame.dataFrame.Bind(wx.EVT_MENU, OnDistAngle, id=G2gd.wxID_ATOMSDISAGL)
Note: See TracChangeset
for help on using the changeset viewer.