Author: tlinnet
Date: Fri Jul 25 17:50:27 2014
New Revision: 24759
URL: http://svn.gna.org/viewcvs/relax?rev=24759&view=rev
Log:
Added an intermediate attempt to show the back_calculated data in the graph
for R1rho R2 as function of the effective field in rotating frame: w_eff.
The graph is aiming for the representation of Figure 2 in Kjaergaard et al
2013. (http://dx.doi.org/10.1021/bi4001062).
The figure can be seen at https://gna.org/support/download.php?file_id=20208.
It becomes clear, that it is not neccessary interpolate through the spin-lock
offset, but it is suffucient to interpolate through the spin-lock field
strengths.
The necessary step was the extraction of the effective effective field in
rotating frame, w_eff.
In earlier attempt is shown at:
http://wiki.nmr-relax.com/File:Matplotlib_52_N_R1_rho_R2eff_w_eff.png
This though show lines for 6 offset values.
The question is how to show the single line of interpolation.
sr #3124(https://gna.org/support/?3124): Grace graphs production for R1rho
analysis with R2_eff as function of Omega_eff.
sr #3138(https://gna.org/support/?3138): Interpolating theta through
spin-lock offset [Omega], rather than spin-lock field strength [w1].
Modified:
branches/r1rho_plotting/specific_analyses/relax_disp/data.py
Modified: branches/r1rho_plotting/specific_analyses/relax_disp/data.py
URL:
http://svn.gna.org/viewcvs/relax/branches/r1rho_plotting/specific_analyses/relax_disp/data.py?rev=24759&r1=24758&r2=24759&view=diff
==============================================================================
--- branches/r1rho_plotting/specific_analyses/relax_disp/data.py
(original)
+++ branches/r1rho_plotting/specific_analyses/relax_disp/data.py Fri
Jul 25 17:50:27 2014
@@ -63,7 +63,7 @@
from lib.errors import RelaxError, RelaxNoSpectraError, RelaxNoSpinError,
RelaxSpinTypeError
from lib.float import isNaN
from lib.io import extract_data, get_file_path, open_write_file, strip,
write_data
-from lib.nmr import frequency_to_ppm, frequency_to_rad_per_s
+from lib.nmr import frequency_to_ppm, frequency_to_ppm_from_rad,
frequency_to_rad_per_s
from lib.physical_constants import g1H, return_gyromagnetic_ratio
from lib.sequence import read_spin_data, write_spin_data
from lib.software.grace import write_xy_data, write_xy_header,
script_grace2images
@@ -1751,11 +1751,10 @@
if not spin.model in [MODEL_R2EFF]:
# Interpolate through spin-lock offset points.
- interpolated_flag, back_calc, spin_lock_offset_new,
chemical_shifts, spin_lock_fields_inter, offsets, tilt_angles, Delta_omega,
w_eff = plot_disp_curves_interpolate_sl_offset(spin=spin, spin_id=spin_id,
si=si, num_points=num_points, extend=extend)
+ interpolated_flag, back_calc, spin_lock_offset_new,
chemical_shifts, spin_lock_fields_inter, offsets_inter, tilt_angles_inter,
Delta_omega_inter, w_eff_inter =
plot_disp_curves_interpolate_sl_offset(spin=spin, spin_id=spin_id, si=si,
num_points=num_points, extend=extend)
else:
back_calc = None
- spin_lock_offset_new = None
# Open the file for writing.
file_path = get_file_path(file_name, dir)
@@ -1785,7 +1784,7 @@
current_spin = proton
# Loop over the spectrometer frequencies and offsets.
- err, data, set_labels, set_colours, x_axis_type_zero, symbols,
symbol_sizes, linetype, linestyle, axis_labels =
plot_disp_curves_loop_frq(exp_type=exp_type, ei=ei,
current_spin=current_spin, spin_id=spin_id, si=si, back_calc=back_calc,
chemical_shifts=chemical_shifts,
spin_lock_fields_inter=spin_lock_fields_inter, offsets=offsets,
tilt_angles=tilt_angles, Delta_omega=Delta_omega, w_eff=w_eff,
interpolated_flag=interpolated_flag, graph_index=graph_index,
colour_order=colour_order, data=data, set_labels=set_labels,
set_colours=set_colours, x_axis_type_zero=x_axis_type_zero, symbols=symbols,
symbol_sizes=symbol_sizes, linetype=linetype, linestyle=linestyle,
axis_labels=axis_labels)
+ err, data, set_labels, set_colours, x_axis_type_zero, symbols,
symbol_sizes, linetype, linestyle, axis_labels =
plot_disp_curves_loop_frq(exp_type=exp_type, ei=ei,
current_spin=current_spin, spin_id=spin_id, si=si, back_calc=back_calc,
chemical_shifts=chemical_shifts,
spin_lock_fields_inter=spin_lock_fields_inter, offsets_inter=offsets_inter,
tilt_angles_inter=tilt_angles_inter, Delta_omega_inter=Delta_omega_inter,
w_eff_inter=w_eff_inter, interpolated_flag=interpolated_flag,
graph_index=graph_index, colour_order=colour_order, data=data,
set_labels=set_labels, set_colours=set_colours,
x_axis_type_zero=x_axis_type_zero, symbols=symbols,
symbol_sizes=symbol_sizes, linetype=linetype, linestyle=linestyle,
axis_labels=axis_labels)
# Increment the graph index.
graph_index += 1
@@ -2026,7 +2025,7 @@
return interpolated_flag, back_calc, spin_lock_offset_new,
chemical_shifts, spin_lock_fields_inter, offsets, tilt_angles, Delta_omega,
w_eff
-def plot_disp_curves_loop_frq(exp_type=None, ei=None, current_spin=None,
spin_id=None, si=None, back_calc=None, chemical_shifts=None,
spin_lock_fields_inter=None, offsets=None, tilt_angles=None,
Delta_omega=None, w_eff=None, interpolated_flag=None, graph_index=None,
colour_order=None, data=None, set_labels=None, set_colours=None,
x_axis_type_zero=None, symbols=None, symbol_sizes=None, linetype=None,
linestyle=None, axis_labels=None):
+def plot_disp_curves_loop_frq(exp_type=None, ei=None, current_spin=None,
spin_id=None, si=None, back_calc=None, chemical_shifts=None,
spin_lock_fields_inter=None, offsets_inter=None, tilt_angles_inter=None,
Delta_omega_inter=None, w_eff_inter=None, interpolated_flag=None,
graph_index=None, colour_order=None, data=None, set_labels=None,
set_colours=None, x_axis_type_zero=None, symbols=None, symbol_sizes=None,
linetype=None, linestyle=None, axis_labels=None):
"""Loop function over the spectrometer frequencies and offsets for 2D
Grace plotting function.
@keyword exp_type: The experiment type.
@@ -2047,14 +2046,14 @@
@type chemical_shifts: rank-3 list of floats
@keyword spin_lock_fields_inter: The interpolated spin-lock field
strengths in Hertz {Ei, Mi, Oi}
@type spin_lock_fields_inter: rank-3 list of floats
- @keyword offsets: The interpolated spin-lock offsets
in rad/s {Ei, Si, Mi, Oi}
- @type offsets: rank-4 list of floats
- @keyword tilt_angles: The rotating frame tilt angles {Ei,
Si, Mi, Oi, Di}
- @type tilt_angles: rank-5 list of floats
- @keyword Delta_omega: The average resonance offset in the
rotating frame in rad/s {Ei, Si, Mi, Oi, Di}
- @type Delta_omega: rank-5 list of floats
- @keyword w_eff: The effective field in rotating
frame in rad/s {Ei, Si, Mi, Oi, Di}.
- @type w_eff: rank-5 list of floats
+ @keyword offsets_inter: The interpolated spin-lock offsets
in rad/s {Ei, Si, Mi, Oi}
+ @type offsets_inter: rank-4 list of floats
+ @keyword tilt_angles_inter: The interpolated rotating frame tilt
angles {Ei, Si, Mi, Oi, Di}
+ @type tilt_angles_inter: rank-5 list of floats
+ @keyword Delta_omega_inter: The interpolated average resonance
offset in the rotating frame in rad/s {Ei, Si, Mi, Oi, Di}
+ @type Delta_omega_inter: rank-5 list of floats
+ @keyword w_eff_inter: The interpolated effective field in
rotating frame in rad/s {Ei, Si, Mi, Oi, Di}.
+ @type w_eff_inter: rank-5 list of floats
@keyword interpolated_flag: Flag telling if the graph should be
interpolated.
@type interpolated_flag: bool
@keyword graph_index: Graph index for xmgrace.
@@ -2227,6 +2226,137 @@
set_index += 1
colour_index += 1
+ # Add the interpolated back calculated data.
+ if interpolated_flag:
+ colour_index = 0
+ for frq, mi in loop_frq(return_indices=True):
+ # Loop over interpolated offset.
+ for oi, offset_rad in enumerate(offsets_inter[ei][si][mi]):
+ # Add a new set for the data at each frequency and offset.
+ #data[graph_index].append([])
+
+ # Convert offset to ppm from rad/s.
+ offset = frequency_to_ppm_from_rad(frq=offset_rad, B0=frq,
isotope=current_spin.isotope)
+
+ # Add a new label.
+ if exp_type in EXP_TYPE_LIST_CPMG:
+ label = "R\\s2eff\\N interpolated curve"
+ else:
+ label = "R\\s1\\xr\\B\\N interpolated curve"
+ if offset != None and frq != None:
+ label += " (%.1f MHz, %.3f ppm)" % (frq / 1e6, offset)
+ elif frq != None:
+ label += " (%.1f MHz)" % (frq / 1e6)
+ elif offset != None:
+ label += " (%.3f ppm)" % (offset)
+ #set_labels[ei].append(label)
+
+ # The other settings.
+ #set_colours[graph_index].append(colour_order[colour_index])
+ #x_axis_type_zero[graph_index].append(True)
+ #if current_spin.model in MODEL_LIST_NUMERIC_CPMG:
+ # symbols[graph_index].append(8)
+ #else:
+ # symbols[graph_index].append(0)
+ #symbol_sizes[graph_index].append(0.20)
+ #linetype[graph_index].append(1)
+ #linestyle[graph_index].append(1)
+
+ # Loop over the dispersion points.
+ for di, point in
enumerate(spin_lock_fields_inter[ei][mi][oi]):
+ # Assign r2eff.
+ r2eff = back_calc[ei][si][mi][oi][di]
+
+ # Skip invalid points (values of 1e100).
+ if r2eff > 1e50:
+ continue
+
+ # Set x_point.
+ x_point = w_eff_inter[ei][si][mi][oi][di]
+
+ theta = tilt_angles_inter[ei][si][mi][oi][di]
+
+ # R_2 = R1rho / sin^2(theta) - R_1 / tan^2(theta) =
(R1rho - R_1 * cos^2(theta) ) / sin^2(theta)
+ y_point = ( r2eff - r1[si][mi]*cos(theta)**2 ) /
sin(theta)**2
+
+ # Add the data.
+ #data[graph_index][set_index].append([x_point, y_point])
+
+ # Handle the errors.
+ #if err:
+ # data[graph_index][set_index][-1].append(None)
+
+ # Increment the graph set index.
+ #set_index += 1
+ #colour_index += 1
+
+ # Add the residuals for statistical comparison.
+ colour_index = 0
+ for frq, offset, mi, oi in loop_frq_offset(exp_type=exp_type,
return_indices=True):
+ # Add a new set for the data at each frequency and offset.
+ data[graph_index].append([])
+
+ # Add a new label.
+ label = "Residuals"
+ if offset != None and frq != None:
+ label += " (%.1f MHz, %.3f ppm)" % (frq / 1e6, offset)
+ elif frq != None:
+ label += " (%.1f MHz)" % (frq / 1e6)
+ elif offset != None:
+ label += " (%.3f ppm)" % (offset)
+ set_labels[ei].append(label)
+
+ # The other settings.
+ set_colours[graph_index].append(colour_order[colour_index])
+ x_axis_type_zero[graph_index].append(True)
+ symbols[graph_index].append(9)
+ symbol_sizes[graph_index].append(0.45)
+ linetype[graph_index].append(1)
+ linestyle[graph_index].append(3)
+
+ # Loop over the dispersion points.
+ for point, di in loop_point(exp_type=exp_type, frq=frq,
offset=offset, return_indices=True):
+ # The data key.
+ key = return_param_key_from_data(exp_type=exp_type, frq=frq,
offset=offset, point=point)
+
+ # No data present.
+ if key not in current_spin.r2eff or not hasattr(current_spin,
'r2eff_bc') or key not in current_spin.r2eff_bc:
+ continue
+
+ # Convert offset to rad/s from ppm.
+ offset_rad = frequency_to_rad_per_s(frq=offset, B0=frq,
isotope=current_spin.isotope)
+ # Calculate the tilt angle.
+ omega1 = point * 2.0 * pi
+ Delta_omega = chemical_shifts[ei][si][mi] - offset_rad
+ if Delta_omega == 0.0:
+ theta = pi / 2.0
+ else:
+ theta = atan2(omega1, Delta_omega)
+
+ # Calculate effective field in rotating frame
+ w_eff = sqrt( Delta_omega*Delta_omega + omega1*omega1 )
+
+ # Set x_point.
+ x_point = w_eff
+
+ # Set y_point. When R_1 is set 0.0.
+ # R_2 = R1rho / sin^2(theta) - R_1 / tan^2(theta) = (R1rho - R_1
* cos^2(theta) ) / sin^2(theta)
+ y_point = ( current_spin.r2eff[key] - r1[si][mi]*cos(theta)**2 )
/ sin(theta)**2
+ y_point_bc = ( current_spin.r2eff_bc[key] -
r1[si][mi]*cos(theta)**2 ) / sin(theta)**2
+ y_point_residual = y_point - y_point_bc
+
+ # Add the data.
+ data[graph_index][set_index].append([x_point, y_point_residual])
+
+ # Handle the errors.
+ if err:
+ err = True
+ y_err_point = ( current_spin.r2eff_err[key] -
r1_err[si][mi]*cos(theta)**2 ) / sin(theta)**2
+ data[graph_index][set_index][-1].append(y_err_point)
+
+ # Increment the graph set index.
+ set_index += 1
+ colour_index += 1
# The axis labels.
axis_labels.append(['\\qEffective field in rotating frame: w_eff \\Q',
'\\qR\\s2\\N\\Q (rad.s\\S-1\\N)'])
r24759 - /branches/r1rho_plotting/specific_analyses/relax_disp/data.py