Changeset 1926


Ignore:
Timestamp:
Jul 13, 2015 4:10:42 PM (7 years ago)
Author:
vondreele
Message:

Work on add Hydrogens

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIgrid.py

    r1925 r1926  
    351351        mainSizer.Add(wx.StaticText(self.panel,-1,'H atom add controls for phase %s:'%(phase['General']['Name'])),
    352352            0,wx.LEFT|wx.TOP,10)
     353        mainSizer.Add(wx.StaticText(self.panel,-1,'NB: Check selections as they may not be correct'),0,WACV|wx.LEFT,10)
    353354        mainSizer.Add(wx.StaticText(self.panel,-1," Atom:  Add # H's          Neighbors, dist"),0,wx.TOP|wx.LEFT,5)
    354355        nHatms = ['0','1','2','3']
     
    359360            nH = 1      #for O atom
    360361            if 'C' in neigh[0] or 'N' in neigh[0]:
    361                 nH = 4-len(neigh[1])
    362             neigh[2] = nH
     362                nH = 4-len(neigh[1][0])
    363363            checks = wx.BoxSizer(wx.HORIZONTAL)
    364364            Ids = []
     
    374374            dataSizer.Add(checks,0,WACV)
    375375            lineSizer = wx.BoxSizer(wx.HORIZONTAL)
    376             for bond in neigh[1]:
     376            for bond in neigh[1][0]:
    377377                lineSizer.Add(wx.StaticText(self.panel,-1,' %s, %.3f'%(bond[0],bond[1])),0,WACV)
    378             dataSizer.Add(lineSizer,0,WACV)
     378            dataSizer.Add(lineSizer,0,WACV|wx.RIGHT,10)
    379379        mainSizer.Add(dataSizer,0,wx.LEFT,5)
    380380
     
    389389        btnSizer.Add(CancelBtn)
    390390        btnSizer.Add((20,20),1)
    391         mainSizer.Add(btnSizer,0,wx.EXPAND|wx.BOTTOM|wx.TOP, 10)
     391        mainSizer.Add(btnSizer,0,wx.BOTTOM|wx.TOP, 10)
     392        size = np.array(self.GetSize())
     393        self.panel.SetupScrolling()
    392394        self.panel.SetSizer(mainSizer)
    393         self.panel.SetupScrolling()
     395        self.panel.SetAutoLayout(1)
     396        size = [size[0]-5,size[1]-20]       #this fiddling is needed for older wx!
     397        self.panel.SetSize(size)
    394398       
    395399    def GetData(self):
  • trunk/GSASIImath.py

    r1925 r1926  
    385385    atTypes = General['AtomTypes']
    386386    Radii = np.array(General['BondRadii'])
     387    DisAglCtls = General['DisAglCtls']   
     388    radiusFactor = DisAglCtls['Factors'][0]
    387389    AtInfo = dict(zip(atTypes,Radii)) #or General['BondRadii']
    388390    Orig = atNames.index(FrstName)
     391    OId = Atoms[Orig][cia+8]
    389392    OType = Atoms[Orig][ct]
    390393    XYZ = getAtomXYZ(Atoms,cx)       
    391394    Neigh = []
     395    Ids = []
    392396    Dx = np.inner(Amat,XYZ-XYZ[Orig]).T
    393397    dist = np.sqrt(np.sum(Dx**2,axis=1))
    394398    sumR = np.array([AtInfo[OType]+AtInfo[atom[ct]] for atom in Atoms])
    395     IndB = ma.nonzero(ma.masked_greater(dist-0.85*sumR,0.))
     399    IndB = ma.nonzero(ma.masked_greater(dist-radiusFactor*sumR,0.))
    396400    for j in IndB[0]:
    397401        if j != Orig:
    398402            Neigh.append([AtNames[j],dist[j]])
    399     return Neigh
     403            Ids.append(Atoms[j][cia+8])
     404    return Neigh,[OId,Ids]
     405   
     406def AddHydrogens(AtLookUp,General,Atoms,AddHydId):
     407    cx,ct,cs,cia = General['AtomPtrs']
     408    Cell = General['Cell'][1:7]
     409    Amat,Bmat = G2lat.cell2AB(Cell)
     410    nBonds = AddHydId[2]+len(AddHydId[1])
     411    Oatom = GetAtomsById(Atoms,AtLookUp,[AddHydId[0],])[0]
     412    OXYZ = np.array(Oatom[cx:cx+3])
     413    Tatoms = GetAtomsById(Atoms,AtLookUp,AddHydId[1])
     414    TXYZ = np.array([tatom[cx:cx+3] for tatom in Tatoms]) #3 x xyz
     415    DX = np.inner(Amat,TXYZ-OXYZ).T
     416    if nBonds == 4:
     417        if AddHydId[2] == 1:
     418            Vec = TXYZ-OXYZ
     419            Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2,axis=0))
     420            Vec = np.sum(Vec/Len,axis=0)
     421            Len = np.sqrt(np.sum(Vec**2))
     422            Hpos = OXYZ-0.98*np.inner(Bmat,Vec).T/Len
     423            return [Hpos,]
     424        elif AddHydId[2] == 2:
     425            print 'add 2 H'
     426            return []
     427        else:
     428            print 'add 3 H'
     429            return []           
     430    elif nBonds == 3:
     431        if AddHydId[2] == 1:
     432            Vec = np.sum(TXYZ-OXYZ,axis=0)               
     433            Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2))
     434            Vec = -0.93*Vec/Len
     435            Hpos = OXYZ+Vec
     436            return [Hpos,]
     437        elif AddHydId[2] == 2:
     438            print 'add 2 H'
     439            return []
     440    else:   #2 bonds
     441        if 'C' in Oatom[ct]:
     442            print 'sp atom',Oatom[ct-1]
     443            return []
     444        elif 'O' in Oatom[ct]:
     445            #idea - construct ring at 0.82 from O atom & find high spot?
     446            print 'sp3 atom ',Oatom[ct-1]
     447            print 'add 1 H'
     448            return []
     449    return []
     450   
     451   
     452       
     453   
    400454       
    401455def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj):   #unfinished & not used
  • trunk/GSASIIphsGUI.py

    r1925 r1926  
    15581558            atomData = data['Atoms']
    15591559            AtNames = [atom[ct-1] for atom in atomData]
    1560             colLabels = [Atoms.GetColLabelValue(c) for c in range(Atoms.GetNumberCols())]
     1560            AtLookUp = G2mth.FillAtomLookUp(atomData,cia+8)
    15611561            Neigh = []
     1562            AddHydIds = []
    15621563            for ind in indx:
    15631564                atom = atomData[ind]
     
    15651566                    continue
    15661567                neigh = [atom[ct-1],G2mth.FindNeighbors(data,atom[ct-1],AtNames),0]
    1567                 if len(neigh[1]) > 3 or (atom[ct] == 'O' and len(neigh[1]) > 1):
     1568                if len(neigh[1][0]) > 3 or (atom[ct] == 'O' and len(neigh[1][0]) > 1):
    15681569                    continue
     1570                nH = 1      #for O atom
     1571                if 'C' in neigh[0] or 'N' in neigh[0]:
     1572                    nH = 4-len(neigh[1][0])
     1573                bonds = dict(neigh[1][0])
     1574                for bond in bonds:
     1575                    if 'C' in neigh[0]:
     1576                        if 'C' in bond and bonds[bond] < 1.42:
     1577                            nH -= 1
     1578                            break
     1579                        elif 'O' in bond and bonds[bond] < 1.3:
     1580                            nH -= 1
     1581                            break
     1582                    elif 'O' in neigh[0] and 'C' in bonds and bonds[bond] < 1.3:
     1583                        nH -= 1
     1584                        break
     1585                neigh[2] = max(0,nH)  #set expected no. H's needed
     1586                AddHydIds.append(neigh[1][1])
    15691587                Neigh.append(neigh)
    15701588            if Neigh:
    15711589                dlg = G2gd.AddHatomDialog(G2frame,Neigh,data)
    15721590                if dlg.ShowModal() == wx.ID_OK:
     1591                    Nat = len(atomData)
    15731592                    Neigh = dlg.GetData()
    1574                     for neigh in Neigh:
    1575                         print neigh
    1576                    
     1593                    for ineigh,neigh in enumerate(Neigh):
     1594                        AddHydIds[ineigh].append(neigh[2])
     1595                        Hxyz = G2mth.AddHydrogens(AtLookUp,generalData,atomData,AddHydIds[ineigh])
     1596                        for X in Hxyz:
     1597                            Nat += 1
     1598                            AtomAdd(X[0],X[1],X[2],'H','H(%d)'%(Nat))
     1599
     1600                SetupGeneral()
     1601                FillAtomsGrid(Atoms)
    15771602                dlg.Destroy()
     1603                G2plt.PlotStructure(G2frame,data)
    15781604            else:
    15791605                wx.MessageBox('No candidates found',caption='Add H atom Error',style=wx.ICON_EXCLAMATION)
Note: See TracChangeset for help on using the changeset viewer.