r11302 - /1.3/generic_fns/pcs.py -- July 15, 2010 - 15:26

Author: bugman
Date: Thu Jul 15 15:26:35 2010
New Revision: 11302

URL: http://svn.gna.org/viewcvs/relax?rev=11302&view=rev
Log:
Wrote the pcs.back_calc() user function back end.


Modified:
    1.3/generic_fns/pcs.py

Modified: 1.3/generic_fns/pcs.py
URL: 
http://svn.gna.org/viewcvs/relax/1.3/generic_fns/pcs.py?rev=11302&r1=11301&r2=11302&view=diff
==============================================================================
--- 1.3/generic_fns/pcs.py (original)
+++ 1.3/generic_fns/pcs.py Thu Jul 15 15:26:35 2010
@@ -24,17 +24,82 @@
 """Module for the manipulation of pseudocontact shift data."""
 
 # Python module imports.
-from math import sqrt
-from numpy import array, float64, zeros
+from math import pi, sqrt
+from numpy import array, float64, ones, zeros
+from numpy.linalg import norm
 import sys
 from warnings import warn
 
 # relax module imports.
 from generic_fns import grace, pipes
+from generic_fns.align_tensor import get_tensor_index
 from generic_fns.mol_res_spin import exists_mol_res_spin_data, return_spin, 
spin_loop
+from maths_fns.pcs import ave_pcs_tensor
+from physical_constants import g1H, pcs_constant
 from relax_errors import RelaxError, RelaxNoPdbError, RelaxNoSequenceError, 
RelaxNoSpinError
 from relax_io import open_write_file, read_spin_data, write_spin_data
 from relax_warnings import RelaxWarning
+
+
+def back_calc(align_id=None):
+    """Back calculate the PCS from the alignment tensor.
+
+    @keyword align_id:      The alignment tensor ID string.
+    @type align_id:         str
+    """
+
+    # Arg check.
+    if align_id and align_id not in cdp.align_ids:
+        raise RelaxError, "The alignment ID '%s' is not in the alignment ID 
list %s." % (align_id, cdp.align_ids)
+
+    # Convert the align IDs to an array, or take all IDs.
+    if align_id:
+        align_ids = [align_id]
+    else:
+        align_ids = cdp.align_ids
+
+    # The weights.
+    weights = ones(cdp.N, float64) / cdp.N
+
+    # Unit vector data structure init.
+    unit_vect = zeros((cdp.N, 3), float64)
+
+    # Loop over the spins.
+    for spin in spin_loop():
+        # Skip spins with no position.
+        if not hasattr(spin, 'pos'):
+            continue
+
+        # Atom positions.
+        pos = spin.pos
+        if type(pos[0]) in [float, float64]:
+            pos = [pos]
+
+        # Loop over the alignments.
+        for id in align_ids:
+            # Vectors.
+            vect = zeros((cdp.N, 3), float64)
+            r = zeros(cdp.N, float64)
+            dj = zeros(cdp.N, float64)
+            for c in range(cdp.N):
+                # The vector.
+                vect[c] = pos - cdp.paramagnetic_centre
+
+                # The length.
+                r[c] = norm(vect[c])
+
+                # Normalise.
+                vect[c] = vect[c] / r[c]
+
+                # Calculate the PCS constant.
+                dj[c] = pcs_constant(cdp.temperature[id], cdp.frq[id] * 2.0 
* pi / g1H, r[c]/1e10)
+
+            # Initialise if necessary.
+            if not hasattr(spin, 'pcs_bc'):
+                spin.pcs_bc = {}
+
+            # Calculate the PCSs (in ppm).
+            spin.pcs_bc[id] = ave_pcs_tensor(dj, vect, cdp.N, 
cdp.align_tensors[get_tensor_index(id)].A, weights=weights) * 1e6
 
 
 def centre(pos=None, atom_id=None, pipe=None, verbosity=1, ave_pos=False, 
force=False):