Author: bugman
Date: Thu Feb 18 11:41:19 2010
New Revision: 10710
URL: http://svn.gna.org/viewcvs/relax?rev=10710&view=rev
Log:
Shifted the RDC and PCS Q-factor calculation code.
The methods previously in the N-state model specific code have been converted
to functions of the
generic_fns.rdc and generic_fns.pcs modules.
Modified:
1.3/generic_fns/pcs.py
1.3/generic_fns/rdc.py
1.3/specific_fns/n_state_model.py
Modified: 1.3/generic_fns/pcs.py
URL:
http://svn.gna.org/viewcvs/relax/1.3/generic_fns/pcs.py?rev=10710&r1=10709&r2=10710&view=diff
==============================================================================
--- 1.3/generic_fns/pcs.py (original)
+++ 1.3/generic_fns/pcs.py Thu Feb 18 11:41:19 2010
@@ -1,6 +1,6 @@
###############################################################################
#
#
-# Copyright (C) 2003-2009 Edward d'Auvergne
#
+# Copyright (C) 2003-2010 Edward d'Auvergne
#
#
#
# This file is part of the program relax.
#
#
#
@@ -25,6 +25,7 @@
# Python module imports.
from copy import deepcopy
+from math import sqrt
from numpy import array, float64, zeros
# relax module imports.
@@ -450,6 +451,46 @@
# Return the index.
return index
+
+
+def q_factors():
+ """Calculate the Q-factors for the PCS data."""
+
+ # Q-factor list.
+ cdp.q_factors_pcs = []
+
+ # Loop over the alignments.
+ for i in xrange(len(cdp.align_tensors)):
+ # Init.
+ pcs2_sum = 0.0
+ sse = 0.0
+
+ # Spin loop.
+ for spin in spin_loop():
+ # Skip deselected spins.
+ if not spin.select:
+ continue
+
+ # Skip spins without PCS data.
+ if not hasattr(spin, 'pcs') or not hasattr(spin, 'pcs_bc') or
spin.pcs[i] == None:
+ continue
+
+ # Sum of squares.
+ sse = sse + (spin.pcs[i] - spin.pcs_bc[i])**2
+
+ # Sum the PCSs squared (for normalisation).
+ pcs2_sum = pcs2_sum + spin.pcs[i]**2
+
+ # The Q-factor for the alignment.
+ Q = sqrt(sse / pcs2_sum)
+ cdp.q_factors_pcs.append(Q)
+
+ # The total Q-factor.
+ cdp.q_pcs = 0.0
+ for Q in cdp.q_factors_pcs:
+ cdp.q_pcs = cdp.q_pcs + Q**2
+ cdp.q_pcs = cdp.q_pcs / len(cdp.q_factors_pcs)
+ cdp.q_pcs = sqrt(cdp.q_pcs)
def read(align_id=None, file=None, dir=None, file_data=None,
spin_id_col=None, mol_name_col=None, res_num_col=None, res_name_col=None,
spin_num_col=None, spin_name_col=None, data_col=None, error_col=None,
sep=None, spin_id=None):
Modified: 1.3/generic_fns/rdc.py
URL:
http://svn.gna.org/viewcvs/relax/1.3/generic_fns/rdc.py?rev=10710&r1=10709&r2=10710&view=diff
==============================================================================
--- 1.3/generic_fns/rdc.py (original)
+++ 1.3/generic_fns/rdc.py Thu Feb 18 11:41:19 2010
@@ -1,6 +1,6 @@
###############################################################################
#
#
-# Copyright (C) 2003-2009 Edward d'Auvergne
#
+# Copyright (C) 2003-2010 Edward d'Auvergne
#
#
#
# This file is part of the program relax.
#
#
#
@@ -25,7 +25,7 @@
# Python module imports.
from copy import deepcopy
-from math import pi
+from math import pi, sqrt
from numpy import float64, ones, zeros
from numpy.linalg import norm
import sys
@@ -408,6 +408,76 @@
# Return the index.
return index
+
+
+def q_factors():
+ """Calculate the Q-factors for the RDC data."""
+
+ # Q-factor list.
+ cdp.q_factors_rdc = []
+ cdp.q_factors_rdc_norm2 = []
+
+ # Loop over the alignments.
+ for i in xrange(len(cdp.align_tensors)):
+ # Init.
+ D2_sum = 0.0
+ sse = 0.0
+
+ # Spin loop.
+ dj = None
+ N = 0
+ for spin in spin_loop():
+ # Skip deselected spins.
+ if not spin.select:
+ continue
+
+ # Skip spins without RDC data.
+ if not hasattr(spin, 'rdc') or not hasattr(spin, 'rdc_bc') or
spin.rdc[i] == None:
+ continue
+
+ # Sum of squares.
+ sse = sse + (spin.rdc[i] - spin.rdc_bc[i])**2
+
+ # Sum the RDCs squared (for one type of normalisation).
+ D2_sum = D2_sum + spin.rdc[i]**2
+
+ # Gyromagnetic ratios.
+ gx = return_gyromagnetic_ratio(spin.heteronuc_type)
+ gh = return_gyromagnetic_ratio(spin.proton_type)
+
+ # Calculate the RDC dipolar constant (in Hertz, and the 3 comes
from the alignment tensor), and append it to the list.
+ dj_new = 3.0/(2.0*pi) * dipolar_constant(gx, gh, spin.r)
+ if dj and dj_new != dj:
+ raise RelaxError("All the RDCs must come from the same
nucleus type.")
+ else:
+ dj = dj_new
+
+ # Increment the number of data sets.
+ N = N + 1
+
+ # Normalisation factor of 2Da^2(4 + 3R)/5.
+ D = dj * cdp.align_tensors[i].A_diag
+ Da = 1.0/3.0 * (D[2, 2] - (D[0, 0]+D[1, 1])/2.0)
+ Dr = 1.0/3.0 * (D[0, 0] - D[1, 1])
+ R = Dr / Da
+ norm = 2.0 * (Da)**2 * (4.0 + 3.0*R**2)/5.0
+ if Da == 0.0:
+ norm = 1e-15
+
+ # The Q-factor for the alignment.
+ Q = sqrt(sse / N / norm)
+ cdp.q_factors_rdc.append(Q)
+ cdp.q_factors_rdc_norm2.append(sqrt(sse / D2_sum))
+
+ # The total Q-factor.
+ cdp.q_rdc = 0.0
+ cdp.q_rdc_norm2 = 0.0
+ for Q in cdp.q_factors_rdc:
+ cdp.q_rdc = cdp.q_rdc + Q**2
+ for Q in cdp.q_factors_rdc_norm2:
+ cdp.q_rdc_norm2 = cdp.q_rdc_norm2 + Q**2
+ cdp.q_rdc = sqrt(cdp.q_rdc / len(cdp.q_factors_rdc))
+ cdp.q_rdc_norm2 = sqrt(cdp.q_rdc_norm2 / len(cdp.q_factors_rdc_norm2))
def read(align_id=None, file=None, dir=None, file_data=None,
spin_id_col=None, mol_name_col=None, res_num_col=None, res_name_col=None,
spin_num_col=None, spin_name_col=None, data_col=None, error_col=None,
sep=None, spin_id=None):
Modified: 1.3/specific_fns/n_state_model.py
URL:
http://svn.gna.org/viewcvs/relax/1.3/specific_fns/n_state_model.py?rev=10710&r1=10709&r2=10710&view=diff
==============================================================================
--- 1.3/specific_fns/n_state_model.py (original)
+++ 1.3/specific_fns/n_state_model.py Thu Feb 18 11:41:19 2010
@@ -39,7 +39,7 @@
from float import isNaN, isInf
import generic_fns
from generic_fns.mol_res_spin import return_spin, spin_loop
-from generic_fns import pipes
+from generic_fns import pcs, pipes, rdc
import generic_fns.structure.geometric
from generic_fns.structure.internal import Internal
import generic_fns.structure.mass
@@ -1099,116 +1099,6 @@
# Return the param number.
return num
-
-
- def _q_factors_rdc(self):
- """Calculate the Q-factors for the RDC data."""
-
- # Q-factor list.
- cdp.q_factors_rdc = []
- cdp.q_factors_rdc_norm2 = []
-
- # Loop over the alignments.
- for i in xrange(len(cdp.align_tensors)):
- # Init.
- D2_sum = 0.0
- sse = 0.0
-
- # Spin loop.
- dj = None
- N = 0
- for spin in spin_loop():
- # Skip deselected spins.
- if not spin.select:
- continue
-
- # Skip spins without RDC data.
- if not hasattr(spin, 'rdc') or not hasattr(spin, 'rdc_bc')
or spin.rdc[i] == None:
- continue
-
- # Sum of squares.
- sse = sse + (spin.rdc[i] - spin.rdc_bc[i])**2
-
- # Sum the RDCs squared (for one type of normalisation).
- D2_sum = D2_sum + spin.rdc[i]**2
-
- # Gyromagnetic ratios.
- gx = return_gyromagnetic_ratio(spin.heteronuc_type)
- gh = return_gyromagnetic_ratio(spin.proton_type)
-
- # Calculate the RDC dipolar constant (in Hertz, and the 3
comes from the alignment tensor), and append it to the list.
- dj_new = 3.0/(2.0*pi) * dipolar_constant(gx, gh, spin.r)
- if dj and dj_new != dj:
- raise RelaxError("All the RDCs must come from the same
nucleus type.")
- else:
- dj = dj_new
-
- # Increment the number of data sets.
- N = N + 1
-
- # Normalisation factor of 2Da^2(4 + 3R)/5.
- D = dj * cdp.align_tensors[i].A_diag
- Da = 1.0/3.0 * (D[2, 2] - (D[0, 0]+D[1, 1])/2.0)
- Dr = 1.0/3.0 * (D[0, 0] - D[1, 1])
- R = Dr / Da
- norm = 2.0 * (Da)**2 * (4.0 + 3.0*R**2)/5.0
- if Da == 0.0:
- norm = 1e-15
-
- # The Q-factor for the alignment.
- Q = sqrt(sse / N / norm)
- cdp.q_factors_rdc.append(Q)
- cdp.q_factors_rdc_norm2.append(sqrt(sse / D2_sum))
-
- # The total Q-factor.
- cdp.q_rdc = 0.0
- cdp.q_rdc_norm2 = 0.0
- for Q in cdp.q_factors_rdc:
- cdp.q_rdc = cdp.q_rdc + Q**2
- for Q in cdp.q_factors_rdc_norm2:
- cdp.q_rdc_norm2 = cdp.q_rdc_norm2 + Q**2
- cdp.q_rdc = sqrt(cdp.q_rdc / len(cdp.q_factors_rdc))
- cdp.q_rdc_norm2 = sqrt(cdp.q_rdc_norm2 /
len(cdp.q_factors_rdc_norm2))
-
-
- def _q_factors_pcs(self):
- """Calculate the Q-factors for the PCS data."""
-
- # Q-factor list.
- cdp.q_factors_pcs = []
-
- # Loop over the alignments.
- for i in xrange(len(cdp.align_tensors)):
- # Init.
- pcs2_sum = 0.0
- sse = 0.0
-
- # Spin loop.
- for spin in spin_loop():
- # Skip deselected spins.
- if not spin.select:
- continue
-
- # Skip spins without PCS data.
- if not hasattr(spin, 'pcs') or not hasattr(spin, 'pcs_bc')
or spin.pcs[i] == None:
- continue
-
- # Sum of squares.
- sse = sse + (spin.pcs[i] - spin.pcs_bc[i])**2
-
- # Sum the PCSs squared (for normalisation).
- pcs2_sum = pcs2_sum + spin.pcs[i]**2
-
- # The Q-factor for the alignment.
- Q = sqrt(sse / pcs2_sum)
- cdp.q_factors_pcs.append(Q)
-
- # The total Q-factor.
- cdp.q_pcs = 0.0
- for Q in cdp.q_factors_pcs:
- cdp.q_pcs = cdp.q_pcs + Q**2
- cdp.q_pcs = cdp.q_pcs / len(cdp.q_factors_pcs)
- cdp.q_pcs = sqrt(cdp.q_pcs)
def _ref_domain(self, ref=None):
@@ -1755,11 +1645,11 @@
# Calculate the RDC Q-factors.
if 'rdc' in data_types:
- self._q_factors_rdc()
+ rdc.q_factors()
# Calculate the PCS Q-factors.
if 'pcs' in data_types:
- self._q_factors_pcs()
+ pcs.q_factors()
def model_statistics(self, model_info=None, spin_id=None,
global_stats=None):
r10710 - in /1.3: generic_fns/pcs.py generic_fns/rdc.py specific_fns/n_state_model.py