r5838 - /1.3/generic_fns/structure/main.py -- April 20, 2008 - 13:26

Author: bugman
Date: Sun Apr 20 13:26:25 2008
New Revision: 5838

URL: http://svn.gna.org/viewcvs/relax?rev=5838&view=rev
Log:
Shifted the XH vector calculation from the xh_vector() structural API method 
to the vector() fn.

This calculation can be done independently from the type of the structural 
data object.


Modified:
    1.3/generic_fns/structure/main.py

Modified: 1.3/generic_fns/structure/main.py
URL: 
http://svn.gna.org/viewcvs/relax/1.3/generic_fns/structure/main.py?rev=5838&r1=5837&r2=5838&view=diff
==============================================================================
--- 1.3/generic_fns/structure/main.py (original)
+++ 1.3/generic_fns/structure/main.py Sun Apr 20 13:26:25 2008
@@ -21,6 +21,8 @@
 
###############################################################################
 
 # Python module imports.
+from math import sqrt
+from numpy import dot
 from os import F_OK, access
 import sys
 from warnings import warn
@@ -33,7 +35,7 @@
 from generic_fns.structure.scientific import Scientific_data
 from relax_errors import RelaxError, RelaxFileError, RelaxNoPipeError, 
RelaxNoSequenceError, RelaxPdbError
 from relax_io import get_file_path
-from relax_warnings import RelaxNoPDBFileWarning
+from relax_warnings import RelaxNoPDBFileWarning, RelaxZeroVectorWarning
 
 
 
@@ -228,33 +230,79 @@
             print "\nCalculating the unit XH vectors from the structure."
 
     # Loop over all the structural data.
+    first_model = None
+    spin_list = []
     for atom in cdp.structure.atom_loop(atom_id=spin_id, 
model_num_flag=True, mol_name_flag=True, res_num_flag=True, 
res_name_flag=True, atom_num_flag=True, atom_name_flag=True, 
element_flag=True, pos_flag=True):
         # Unpack the data.
         model_num, mol_name, res_num, res_name, atom_num, atom_name, 
element, pos = atom
 
         # The spin identification string.
         spin_id = generate_spin_id(mol_name, res_num, res_name, atom_num, 
atom_name)
-        print spin_id
 
         # Get the corresponding spin.
         spin = return_spin(selection=spin_id)
 
-    return
-
-    # Loop over the sequence.
-    for spin in spin_loop(spin_id):
-        # Skip unselected residues.
+        # Proton name does not match, so skip the atom.
+        if proton != atom_name:
+            continue
+
+        # No matching spin, so skip the atom.
+        if not spin:
+            continue
+
+        # Skip deselected spins.
         if not spin.select:
             continue
 
+        # The first model.
+        if first_model == None:
+            first_model = model_num
+
+        # The XH vector already exists.
+        if first_model == model_num and hasattr(spin, 'xh_vect'):
+            warn(RelaxWarning("The XH vector for the spin " + `spin_id` + " 
already exists"))
+            continue
+
+        # Store the spin_id.
+        if first_model == model_num:
+            spin_list.append(spin_id)
+
         # Set the attached proton name.
-        spin.attached_proton = proton
-
-        # Calculate the vector.
-        vector = cdp.structure.xh_vector(spin)
+        if not hasattr(spin, 'attached_proton'):
+            spin.attached_proton = proton
+        elif spin.attached_proton != proton:
+            raise RelaxError, "The attached proton " + 
`spin.attached_proton` + " does not match the proton argument " + `proton` + 
"."
+
+        # Set the heteronucleus and proton positions.
+        posX = spin.pos
+        posH = pos
+
+        # Calculate the XH bond vector.
+        vector = posH - posX
+
+        # Unit vector.
+        if unit:
+            # Normalisation factor.
+            norm_factor = sqrt(dot(vector, vector))
+
+            # Test for zero length.
+            if norm_factor == 0.0:
+                warn(RelaxZeroVectorWarning(spin_id))
+
+            # Calculate the normalised vector.
+            else:
+                vector = vector / norm_factor
 
         # Set the vector and deselect the spin if the vector doesn't exist.
-        if vector != None:
-            spin.xh_vect = vector    
+        if not hasattr(spin, 'xh_vect'):
+            spin.xh_vect = vector
         else:
-            spin.select = False
+            spin.xh_vect = spin.xh_vect + vector
+
+    # Average the XH vectors for all models.
+    for spin_id in spin_list:
+        # Get the spin.
+        spin = return_spin(selection=spin_id)
+
+        # Average the vector.
+        spin.xh_vect = spin.xh_vect / cdp.structure.num_models()