r9651 - /1.3/generic_fns/structure/internal.py -- October 08, 2009 - 09:24

Author: bugman
Date: Thu Oct  8 09:24:50 2009
New Revision: 9651

URL: http://svn.gna.org/viewcvs/relax?rev=9651&view=rev
Log:
Significant improvements to the __find_bonded_atoms() method.

The maximum number of bonds an element can have is now taken into account so 
that protons are not
thought to be attached to 2 carbons within a 2 Angstrom radius!  The atoms 
closest to the atom of
interest are now connected.  Geometry constraints or other advanced 
techniques are not yet
implemented.


Modified:
    1.3/generic_fns/structure/internal.py

Modified: 1.3/generic_fns/structure/internal.py
URL: 
http://svn.gna.org/viewcvs/relax/1.3/generic_fns/structure/internal.py?rev=9651&r1=9650&r2=9651&view=diff
==============================================================================
--- 1.3/generic_fns/structure/internal.py (original)
+++ 1.3/generic_fns/structure/internal.py Thu Oct  8 09:24:50 2009
@@ -137,6 +137,9 @@
         centre = array([mol.x[index], mol.y[index], mol.z[index]], float64)
 
         # Atom loop.
+        dist_list = []
+        connect_list = {}
+        element_list = {}
         for i in xrange(len(mol.atom_num)):
             # Skip proton to proton bonds!
             if mol.element[index] == 'H' and mol.element[i] == 'H':
@@ -148,9 +151,34 @@
             # The distance from the centre.
             dist = linalg.norm(centre-pos)
 
-            # Connect the atoms if within the radius value.
+            # The atom is within the radius.
             if dist < radius:
-                mol.atom_connect(index, i)
+                # Store the distance.
+                dist_list.append(dist)
+
+                # Store the atom index.
+                connect_list[dist] = i
+
+                # Store the element type.
+                element_list[dist] = mol.element[i]
+
+        # The maximum number of allowed covalent bonds.
+        max_conn = 1000   # Ridiculous default!
+        if mol.element[index] == 'H':
+            max_conn = 1
+        elif mol.element[index] == 'O':
+            max_conn = 2
+        elif mol.element[index] == 'N':
+            max_conn = 3
+        elif mol.element[index] == 'C':
+            max_conn = 4
+
+        # Sort.
+        dist_list.sort()
+
+        # Loop over the max number of connections.
+        for i in range(max_conn):
+            mol.atom_connect(index, connect_list[dist_list[i]])
 
 
     def __get_chemical_name(self, hetID):