Author: bugman
Date: Fri May 31 15:02:40 2013
New Revision: 19827
URL: http://svn.gna.org/viewcvs/relax?rev=19827&view=rev
Log:
Improvements for the phi_ex and dw relaxation dispersion model parameters.
These are now stored with the units of ppm^2 and ppm respectively. The
conversion to (rad/s)^2 and
rad/s units respectively now is spin specific, allowing mixed spin types (1H,
13C, 15N, etc.) to be
analysed simultaneously.
Modified:
branches/relax_disp/specific_analyses/relax_disp/__init__.py
branches/relax_disp/specific_analyses/relax_disp/disp_data.py
branches/relax_disp/specific_analyses/relax_disp/parameters.py
branches/relax_disp/target_functions/relax_disp.py
Modified: branches/relax_disp/specific_analyses/relax_disp/__init__.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/specific_analyses/relax_disp/__init__.py?rev=19827&r1=19826&r2=19827&view=diff
==============================================================================
--- branches/relax_disp/specific_analyses/relax_disp/__init__.py (original)
+++ branches/relax_disp/specific_analyses/relax_disp/__init__.py Fri May 31
15:02:40 2013
@@ -95,7 +95,7 @@
self.PARAMS.add('r2', scope='spin', default=15.0, desc='The
transversal relaxation rate', set='params', py_type=list,
grace_string='\\qR\\s2\\N\\Q (rad.s\\S-1\\N)', err=True, sim=True)
self.PARAMS.add('pA', scope='spin', default=0.5, desc='The
population for state A', set='params', py_type=float,
grace_string='\qp\sA\N\Q', err=True, sim=True)
self.PARAMS.add('pB', scope='spin', default=0.5, desc='The
population for state B', set='params', py_type=float,
grace_string='\qp\sB\N\Q', err=True, sim=True)
- self.PARAMS.add('phi_ex', scope='spin', default=5.0, desc='The
pA.pB.dw**2 value scaled by wH (phi_ex = pA * pB * Delta_omega**2 /
omega_H**2)', set='params', py_type=float, grace_string='\\xF\\B\\sex\\N\\q
(p\\sA\\N.p\\sB\\N.\\xDw\\B\\S2\\N / \\xw\\B\\sH\\N\\S2\\N)', err=True,
sim=True)
+ self.PARAMS.add('phi_ex', scope='spin', default=5.0, desc='The
phi_ex = pA.pB.dw**2 value (ppm^2)', set='params', py_type=float,
grace_string='\\xF\\B\\sex\\N = \\q p\\sA\\N.p\\sB\\N.\\xDw\\B\\S2\\N\\Q
(ppm\\S2\\N)', err=True, sim=True)
self.PARAMS.add('dw', scope='spin', default=0.0, desc='The chemical
shift difference between states A and B (in ppm)', set='params',
py_type=float, grace_string='\q\\xDw\f{}\Q (ppm)', err=True, sim=True)
self.PARAMS.add('kex', scope='spin', default=10000.0, desc='The
exchange rate', set='params', py_type=float, grace_string='\\qk\\sex\\N\\Q
(rad.s\\S-1\\N)', err=True, sim=True)
self.PARAMS.add('r2a', scope='spin', default=15.0, desc='The
transversal relaxation rate for state A', set='params', py_type=float,
grace_string='\\qR\\s2,A\\N\\Q (rad.s\\S-1\\N)', err=True, sim=True)
@@ -124,10 +124,10 @@
scaling_matrix = assemble_scaling_matrix(spins=[spin], scaling=False)
# Initialise the data structures for the target function.
- values, errors, missing = return_r2eff_arrays(spins=[spin],
spin_ids=[spin_id], fields=cdp.spectrometer_frq_list,
field_count=cdp.spectrometer_frq_count)
+ values, errors, missing, frqs = return_r2eff_arrays(spins=[spin],
spin_ids=[spin_id], fields=cdp.spectrometer_frq_list,
field_count=cdp.spectrometer_frq_count)
# Initialise the relaxation dispersion fit functions.
- model = Dispersion(model=cdp.model,
num_params=param_num(spins=[spin]), num_spins=1,
num_frq=cdp.spectrometer_frq_count, num_disp_points=cdp.dispersion_points,
values=values, errors=errors, missing=missing,
frqs=cdp.spectrometer_frq_list, cpmg_frqs=return_cpmg_frqs(ref_flag=False),
spin_lock_nu1=return_spin_lock_nu1(ref_flag=False),
scaling_matrix=scaling_matrix)
+ model = Dispersion(model=cdp.model,
num_params=param_num(spins=[spin]), num_spins=1,
num_frq=cdp.spectrometer_frq_count, num_disp_points=cdp.dispersion_points,
values=values, errors=errors, missing=missing, frqs=frqs,
cpmg_frqs=return_cpmg_frqs(ref_flag=False),
spin_lock_nu1=return_spin_lock_nu1(ref_flag=False),
scaling_matrix=scaling_matrix)
# Make a single function call. This will cause back calculation and
the data will be stored in the class instance.
model.func(param_vector)
@@ -379,15 +379,15 @@
lower.append(0.0)
upper.append(1.0)
- # The pA.pB.dw**2/wH**2 parameter.
+ # The pA.pB.dw**2 parameter.
elif spin.params[i] == 'phi_ex':
- lower.append(1e-20)
- upper.append(1e-17)
+ lower.append(0.0)
+ upper.append(10.0)
# Chemical shift difference between states A and B.
elif spin.params[i] == 'dw':
- lower.append(1.0)
- upper.append(10000.0)
+ lower.append(0.0)
+ upper.append(10.0)
# Exchange rates.
elif spin.params[i] in ['kex', 'ka']:
@@ -1105,7 +1105,7 @@
# Loop over the spin blocks.
for spins, spin_ids in self.model_loop():
# The R2eff/R1rho data.
- values, errors, missing = return_r2eff_arrays(spins=spins,
spin_ids=spin_ids, fields=cdp.spectrometer_frq_list,
field_count=cdp.spectrometer_frq_count)
+ values, errors, missing, frqs = return_r2eff_arrays(spins=spins,
spin_ids=spin_ids, fields=cdp.spectrometer_frq_list,
field_count=cdp.spectrometer_frq_count)
# Create the initial parameter vector.
param_vector = assemble_param_vector(spins=spins)
@@ -1138,7 +1138,7 @@
print("Unconstrained grid search size: %s (constraints
may decrease this size).\n" % grid_size)
# Initialise the function to minimise.
- model = Dispersion(model=cdp.model,
num_params=param_num(spins=spins), num_spins=len(spins),
num_frq=cdp.spectrometer_frq_count, num_disp_points=cdp.dispersion_points,
values=values, errors=errors, missing=missing,
frqs=cdp.spectrometer_frq_list, cpmg_frqs=return_cpmg_frqs(ref_flag=False),
spin_lock_nu1=return_spin_lock_nu1(ref_flag=False),
scaling_matrix=scaling_matrix)
+ model = Dispersion(model=cdp.model,
num_params=param_num(spins=spins), num_spins=len(spins),
num_frq=cdp.spectrometer_frq_count, num_disp_points=cdp.dispersion_points,
values=values, errors=errors, missing=missing, frqs=frqs,
cpmg_frqs=return_cpmg_frqs(ref_flag=False),
spin_lock_nu1=return_spin_lock_nu1(ref_flag=False),
scaling_matrix=scaling_matrix)
# Grid search.
if search('^[Gg]rid', min_algor):
Modified: branches/relax_disp/specific_analyses/relax_disp/disp_data.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/specific_analyses/relax_disp/disp_data.py?rev=19827&r1=19826&r2=19827&view=diff
==============================================================================
--- branches/relax_disp/specific_analyses/relax_disp/disp_data.py (original)
+++ branches/relax_disp/specific_analyses/relax_disp/disp_data.py Fri May 31
15:02:40 2013
@@ -32,12 +32,13 @@
"""
# Python module imports.
-from math import sqrt
+from math import pi, sqrt
from numpy import float64, int32, ones, zeros
# relax module imports.
-from lib.errors import RelaxError, RelaxNoSpectraError
+from lib.errors import RelaxError, RelaxNoSpectraError, RelaxSpinTypeError
from lib.list import count_unique_elements, unique_elements
+from lib.physical_constants import g1H, return_gyromagnetic_ratio
from pipe_control.mol_res_spin import return_spin, spin_loop
from specific_analyses.relax_disp.variables import CPMG_EXP, FIXED_TIME_EXP,
R1RHO_EXP
@@ -539,8 +540,8 @@
@type field_count: int
@keyword sim_index: The index of the simulation to return the data
of. This should be None if the normal data is required.
@type sim_index: None or int
- @return: The numpy array structures of the R2eff/R1rho
values, errors and missing data. For each structure, the first dimension
corresponds to the spins of a spin block, the second to the spectrometer
field strength, and the third is the dispersion points.
- @rtype: numpy rank-3 float array, numpy rank-3 float
array, numpy rank-3 int array
+ @return: The numpy array structures of the R2eff/R1rho
values, errors, missing data, and corresponding Larmor frequencies. For each
structure, the first dimension corresponds to the spins of a spin block, the
second to the spectrometer field strength, and the third is the dispersion
points. For the Larmor frequency structure, the third dimension is omitted.
+ @rtype: numpy rank-3 float array, numpy rank-3 float
array, numpy rank-3 int array, numpy rank-2 int array
"""
# The spin count.
@@ -550,6 +551,7 @@
values = zeros((spin_num, field_count, cdp.dispersion_points), float64)
errors = ones((spin_num, field_count, cdp.dispersion_points), float64)
missing = ones((spin_num, field_count, cdp.dispersion_points), int32)
+ frqs = zeros((spin_num, field_count), float64)
# Pack the R2eff/R1rho data.
data_flag = False
@@ -562,6 +564,10 @@
continue
data_flag = True
+ # No isotope information.
+ if not hasattr(spin, 'isotope'):
+ raise RelaxSpinTypeError(spin_id=spin_ids[spin_index])
+
# The keys.
keys = list(spin.r2eff.keys())
@@ -574,6 +580,9 @@
# The key.
key = return_param_key_from_data(frq=frq, point=point)
+ # The Larmor frequency for this spin and field strength (in
MHz*2pi to speed up the ppm to rad/s conversion).
+ frqs[spin_index, frq_index] = 2.0 * pi * frq / g1H *
return_gyromagnetic_ratio(spin.isotope) * 1e-6
+
# Missing data.
if key not in spin.r2eff:
continue
@@ -595,7 +604,7 @@
raise RelaxError("No R2eff/R1rho data could be found for the spin
cluster %s." % spin_ids)
# Return the structures.
- return values, errors, missing
+ return values, errors, missing, frqs
def return_spin_lock_nu1(ref_flag=True):
Modified: branches/relax_disp/specific_analyses/relax_disp/parameters.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/specific_analyses/relax_disp/parameters.py?rev=19827&r1=19826&r2=19827&view=diff
==============================================================================
--- branches/relax_disp/specific_analyses/relax_disp/parameters.py (original)
+++ branches/relax_disp/specific_analyses/relax_disp/parameters.py Fri May 31
15:02:40 2013
@@ -134,7 +134,7 @@
else:
param_vector.append(spin.pA)
- # The pA.pB.dw**2/wH**2 parameter.
+ # The pA.pB.dw**2 parameter.
elif spin.params[i] == 'phi_ex':
if sim_index != None:
param_vector.append(spin.phi_ex_sim[sim_index])
@@ -236,9 +236,9 @@
elif spin.params[i] == 'pA':
scaling_matrix[param_index, param_index] = 1
- # The pA.pB.dw**2/wH**2 parameter.
+ # The pA.pB.dw**2 parameter.
elif spin.params[i] == 'phi_ex':
- scaling_matrix[param_index, param_index] = 1e-18
+ scaling_matrix[param_index, param_index] = 1
# Chemical shift difference between states A and B scaling.
elif spin.params[i] == 'dw':
@@ -365,7 +365,7 @@
else:
spin.pA = param_vector[param_index]
- # The pA.pB.dw**2/wH**2 parameter.
+ # The pA.pB.dw**2 parameter.
if spin.params[i] == 'phi_ex':
if sim_index != None:
spin.phi_ex_sim[sim_index] =
param_vector[param_index]
@@ -500,7 +500,7 @@
b.append(0.5 / scaling_matrix[i, i])
j += 2
- # The pA.pB.dw**2/wH**2 parameter (phi_ex >= 0).
+ # The pA.pB.dw**2 parameter (phi_ex >= 0).
elif spin.params[i] == 'phi_ex':
A.append(zero_array * 0.0)
A[j][i] = 1.0
Modified: branches/relax_disp/target_functions/relax_disp.py
URL:
http://svn.gna.org/viewcvs/relax/branches/relax_disp/target_functions/relax_disp.py?rev=19827&r1=19826&r2=19827&view=diff
==============================================================================
--- branches/relax_disp/target_functions/relax_disp.py (original)
+++ branches/relax_disp/target_functions/relax_disp.py Fri May 31 15:02:40
2013
@@ -63,8 +63,8 @@
@type errors: numpy rank-3 float array
@keyword missing: The data structure indicating missing
R2eff/R1rho data. The three dimensions must correspond to those of the
values argument.
@type missing: numpy rank-3 int array
- @keyword frqs: The spectrometer frequencies in Hz.
- @type frqs: numpy rank-1 float array
+ @keyword frqs: The spin Larmor frequencies (in MHz*2pi
to speed up the ppm to rad/s conversion). The dimensions correspond to the
first two of the value, error and missing structures.
+ @type frqs: numpy rank-2 float array
@keyword cpmg_frqs: The CPMG frequencies in Hertz for each
separate dispersion point. This will be ignored for R1rho experiments.
@type cpmg_frqs: numpy rank-1 float array
@keyword spin_lock_nu1: The spin-lock field strengths in Hertz
for each separate dispersion point. This will be ignored for CPMG
experiments.
@@ -142,7 +142,7 @@
# Loop over the spectrometer frequencies.
for frq_index in range(self.num_frq):
# Convert dw from ppm to rad/s.
- dw_frq = dw * self.frqs[frq_index]
+ dw_frq = dw * self.frqs[spin_index, frq_index]
# Back calculate the R2eff values.
r2eff_CR72(r20=R20[frq_index], pA=pA, dw=dw_frq, kex=kex,
cpmg_frqs=self.cpmg_frqs, back_calc=self.back_calc[spin_index, frq_index],
num_points=self.num_disp_points)
@@ -185,8 +185,8 @@
for spin_index in range(self.num_spins):
# Loop over the spectrometer frequencies.
for frq_index in range(self.num_frq):
- # Spectrometer field strength scaling.
- phi_ex_scaled = phi_ex * self.frqs[frq_index]**2
+ # Convert phi_ex from ppm^2 to (rad/s)^2.
+ phi_ex_scaled = phi_ex * self.frqs[spin_index, frq_index]**2
# Back calculate the R2eff values.
r2eff_LM63(r20=R20[frq_index], phi_ex=phi_ex_scaled,
kex=kex, cpmg_frqs=self.cpmg_frqs, back_calc=self.back_calc[spin_index,
frq_index], num_points=self.num_disp_points)
r19827 - in /branches/relax_disp: specific_analyses/relax_disp/ target_functions/