Changeset 1926

Timestamp:

Jul 13, 2015 4:10:42 PM (7 years ago)

Author:

vondreele

Message:

Work on add Hydrogens

Location:

trunk

Files:

• GSASIIgrid.py (modified) (4 diffs)
• GSASIImath.py (modified) (1 diff)
• GSASIIphsGUI.py (modified) (2 diffs)

trunk/GSASIIgrid.py

@@ -351 +351 @@
         mainSizer.Add(wx.StaticText(self.panel,-1,'H atom add controls for phase %s:'%(phase['General']['Name'])),
             0,wx.LEFT|wx.TOP,10)
+        mainSizer.Add(wx.StaticText(self.panel,-1,'NB: Check selections as they may not be correct'),0,WACV|wx.LEFT,10)
         mainSizer.Add(wx.StaticText(self.panel,-1," Atom:  Add # H's          Neighbors, dist"),0,wx.TOP|wx.LEFT,5)
         nHatms = ['0','1','2','3']
@@ -359 +360 @@
             nH = 1      #for O atom
             if 'C' in neigh[0] or 'N' in neigh[0]:
-                nH = 4-len(neigh[1])
-            neigh[2] = nH
+                nH = 4-len(neigh[1][0])
             checks = wx.BoxSizer(wx.HORIZONTAL)
             Ids = []
@@ -374 +374 @@
             dataSizer.Add(checks,0,WACV)
             lineSizer = wx.BoxSizer(wx.HORIZONTAL)
-            for bond in neigh[1]:
+            for bond in neigh[1][0]:
                 lineSizer.Add(wx.StaticText(self.panel,-1,' %s, %.3f'%(bond[0],bond[1])),0,WACV)
-            dataSizer.Add(lineSizer,0,WACV)
+            dataSizer.Add(lineSizer,0,WACV|wx.RIGHT,10)
         mainSizer.Add(dataSizer,0,wx.LEFT,5)
 
@@ -389 +389 @@
         btnSizer.Add(CancelBtn)
         btnSizer.Add((20,20),1)
-        mainSizer.Add(btnSizer,0,wx.EXPAND|wx.BOTTOM|wx.TOP, 10)
+        mainSizer.Add(btnSizer,0,wx.BOTTOM|wx.TOP, 10)
+        size = np.array(self.GetSize())
+        self.panel.SetupScrolling()
         self.panel.SetSizer(mainSizer)
-        self.panel.SetupScrolling()
+        self.panel.SetAutoLayout(1)
+        size = [size[0]-5,size[1]-20]       #this fiddling is needed for older wx!
+        self.panel.SetSize(size)
         
     def GetData(self):

trunk/GSASIImath.py

@@ -385 +385 @@
     atTypes = General['AtomTypes']
     Radii = np.array(General['BondRadii'])
+    DisAglCtls = General['DisAglCtls']    
+    radiusFactor = DisAglCtls['Factors'][0]
     AtInfo = dict(zip(atTypes,Radii)) #or General['BondRadii']
     Orig = atNames.index(FrstName)
+    OId = Atoms[Orig][cia+8]
     OType = Atoms[Orig][ct]
     XYZ = getAtomXYZ(Atoms,cx)        
     Neigh = []
+    Ids = []
     Dx = np.inner(Amat,XYZ-XYZ[Orig]).T
     dist = np.sqrt(np.sum(Dx**2,axis=1))
     sumR = np.array([AtInfo[OType]+AtInfo[atom[ct]] for atom in Atoms])
-    IndB = ma.nonzero(ma.masked_greater(dist-0.85*sumR,0.))
+    IndB = ma.nonzero(ma.masked_greater(dist-radiusFactor*sumR,0.))
     for j in IndB[0]:
         if j != Orig:
             Neigh.append([AtNames[j],dist[j]])
-    return Neigh
+            Ids.append(Atoms[j][cia+8])
+    return Neigh,[OId,Ids]
+    
+def AddHydrogens(AtLookUp,General,Atoms,AddHydId):
+    cx,ct,cs,cia = General['AtomPtrs']
+    Cell = General['Cell'][1:7]
+    Amat,Bmat = G2lat.cell2AB(Cell)
+    nBonds = AddHydId[2]+len(AddHydId[1])
+    Oatom = GetAtomsById(Atoms,AtLookUp,[AddHydId[0],])[0]
+    OXYZ = np.array(Oatom[cx:cx+3])
+    Tatoms = GetAtomsById(Atoms,AtLookUp,AddHydId[1])
+    TXYZ = np.array([tatom[cx:cx+3] for tatom in Tatoms]) #3 x xyz
+    DX = np.inner(Amat,TXYZ-OXYZ).T
+    if nBonds == 4:
+        if AddHydId[2] == 1:
+            Vec = TXYZ-OXYZ
+            Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2,axis=0))
+            Vec = np.sum(Vec/Len,axis=0)
+            Len = np.sqrt(np.sum(Vec**2))
+            Hpos = OXYZ-0.98*np.inner(Bmat,Vec).T/Len
+            return [Hpos,]
+        elif AddHydId[2] == 2:
+            print 'add 2 H'
+            return []
+        else:
+            print 'add 3 H'
+            return []            
+    elif nBonds == 3:
+        if AddHydId[2] == 1:
+            Vec = np.sum(TXYZ-OXYZ,axis=0)                
+            Len = np.sqrt(np.sum(np.inner(Amat,Vec).T**2))
+            Vec = -0.93*Vec/Len
+            Hpos = OXYZ+Vec
+            return [Hpos,]
+        elif AddHydId[2] == 2:
+            print 'add 2 H'
+            return []
+    else:   #2 bonds
+        if 'C' in Oatom[ct]:
+            print 'sp atom',Oatom[ct-1]
+            return []
+        elif 'O' in Oatom[ct]:
+            #idea - construct ring at 0.82 from O atom & find high spot?
+            print 'sp3 atom ',Oatom[ct-1]
+            print 'add 1 H'
+            return []
+    return []
+    
+    
+        
+    
         
 def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj):   #unfinished & not used

trunk/GSASIIphsGUI.py

@@ -1558 +1558 @@
             atomData = data['Atoms']
             AtNames = [atom[ct-1] for atom in atomData]
-            colLabels = [Atoms.GetColLabelValue(c) for c in range(Atoms.GetNumberCols())]
+            AtLookUp = G2mth.FillAtomLookUp(atomData,cia+8)
             Neigh = []
+            AddHydIds = []
             for ind in indx:
                 atom = atomData[ind]
@@ -1565 +1566 @@
                     continue
                 neigh = [atom[ct-1],G2mth.FindNeighbors(data,atom[ct-1],AtNames),0]
-                if len(neigh[1]) > 3 or (atom[ct] == 'O' and len(neigh[1]) > 1):
+                if len(neigh[1][0]) > 3 or (atom[ct] == 'O' and len(neigh[1][0]) > 1):
                     continue
+                nH = 1      #for O atom
+                if 'C' in neigh[0] or 'N' in neigh[0]:
+                    nH = 4-len(neigh[1][0])
+                bonds = dict(neigh[1][0])
+                for bond in bonds:
+                    if 'C' in neigh[0]:
+                        if 'C' in bond and bonds[bond] < 1.42:
+                            nH -= 1
+                            break
+                        elif 'O' in bond and bonds[bond] < 1.3:
+                            nH -= 1
+                            break
+                    elif 'O' in neigh[0] and 'C' in bonds and bonds[bond] < 1.3:
+                        nH -= 1
+                        break
+                neigh[2] = max(0,nH)  #set expected no. H's needed
+                AddHydIds.append(neigh[1][1])
                 Neigh.append(neigh)
             if Neigh:
                 dlg = G2gd.AddHatomDialog(G2frame,Neigh,data)
                 if dlg.ShowModal() == wx.ID_OK:
+                    Nat = len(atomData)
                     Neigh = dlg.GetData()
-                    for neigh in Neigh:
-                        print neigh
-                    
+                    for ineigh,neigh in enumerate(Neigh):
+                        AddHydIds[ineigh].append(neigh[2])
+                        Hxyz = G2mth.AddHydrogens(AtLookUp,generalData,atomData,AddHydIds[ineigh])
+                        for X in Hxyz:
+                            Nat += 1
+                            AtomAdd(X[0],X[1],X[2],'H','H(%d)'%(Nat))
+
+                SetupGeneral()
+                FillAtomsGrid(Atoms)
                 dlg.Destroy()
+                G2plt.PlotStructure(G2frame,data)
             else:
                 wx.MessageBox('No candidates found',caption='Add H atom Error',style=wx.ICON_EXCLAMATION)