Trusted, Bureaucrats
4,223
edits
Changes
→References: Switched to a labelled section transclusion for the citations.
__TOC__
== About ==
The production to these figures relates to the Suppport Request:<br>
[https://gna.org/support/?3124 sr #3124: Grace graphs production for R1rho analysis with R2_eff as function of Omega_eff]
== References ==
[http://www.nmr-relax.com/manual/Dispersion_model_summary.html Refer to the manual for parameter explanation]
* {{#lst:Citations|Evenäs01}}
* {{#lst:Citations|KempfLoria04}}
* {{#lst:Citations|Massi05}}
* {{#lst:Citations|Palmer01}}
* {{#lst:Citations|PalmerMassi06}}
=== Figures ===
Ref [1], Figure 1.b.
'''The bell-curves''' As function of angle calculation.
Ref [1], Figure 1.c.
The wanted graph. No clear "name" for the calculated parameter.
Ref [2], Equation 27.
Here the calculated value is noted as: {{:Reff}} = {{:R1rho}} / sin<sup>2</sup>(θ) - {{:R1}} / tan<sup>2</sup>(θ) = {{:R2zero}} + {{:Rex}}, where {{:R2zero}} refers to {{:R1rhoprime}} as seen at [[DPL94]]
Ref [3], Equation 20.
Here the calculated value is noted as: {{:R2}} = {{:R1rho}} / sin<sup>2</sup>(θ) - {{:R1}} / tan<sup>2</sup>(θ). Figure 11+16, would be the reference.
Ref [4], Equation 43. {{:Reff}} = {{:R1rho}} / sin<sup>2</sup>(θ) - {{:R1}} / tan<sup>2</sup>(θ).
Ref [5], Material and Methods, page 740. Here the calculated value is noted as: {{:R2}}: {{:R2}} = {{:R2zero}} + {{:Rex}}. Figure 4 would be the wished graphs.
A little table of conversion then gives
<source lang="text">
Relax equation | Relax store | Articles
---------------------------------------------------------------
R1rho' spin.r2 R^{0}_2 or Bar{R}_2
Fitted pars Not stored R_ex
R1rho spin.r2eff R1rho
R_1 spin.ri_data['R1'] R_1 or Bar{R}_1
</source>
The parameter is called R_2 or R_eff in the articles.
Since reff is not used in relax, this could be used?
A description could be:
* The effective rate
* The effective transverse relaxation rate constant
* The effective relaxation rate constant.
== Make graphs ==
=== The outcome ===
[[File:Matplotlib 52 N R1 rho theta sep.png|center|upright=2|Figure 1]]
[[File:Matplotlib 52 N R1 rho R2eff w eff.png|center|upright=2|Figure 2]]
[[File:Matplotlib 52 N R1 rho R2eff disp.png|center|upright=2|Figure 3]]
== To run ==
<source lang="bash">
relax -p r1rhor2eff.py
</source>
=== Code ===
{{Collapsible script| title == Code ==File: '''r1rhor2eff.py'''Python script<source | lang ="python">| script =
### python imports
import sys
import os
from math import cos, sin, sqrt, pifrom numpy import array, float64
### plotting facility.
import matplotlib.pyplot as plt
# Ordered dictionary
import collections
### relax modules
# Import some tools to loop over the spins.
from pipe_control.mol_res_spin import return_spin, spin_loop
# Import method to calculate the R1_rho offset data
from specific_analyses.relax_disp.disp_data import calc_rotating_frame_params, generate_r20_key, loop_exp_frq, loop_exp_frq_offset, loop_point, return_param_key_from_data, return_spin_lock_nu1from specific_analyses.relax_disp import optimisationfrom lib.nmr import frequency_to_Hz, frequency_to_ppm, frequency_to_rad_per_s
###############
# You have to provide a DPL94 results state file
res_folder = "resultsR1"
#res_folder = "results_clustering"
res_state = os.path.join(res_folder, "DPL94", "results")
spin_inte = ":5244@N"# Make a fake spin, from the spin of interestfake_spin_inte = spin_inte.replace("N","X") # Interpolate graph settings#num_points=1000, extend=500.0num_points=100extend=5000.0 ################
spin_inte_rep = spin_inte.replace('#', '_').replace(':', '_').replace('@', '_')
# Load the state
state.load(res_state, force=True)
# Get the dictionary key
for exp_type, frq in loop_exp_frq():
r20_key = generate_r20_key(exp_type=exp_type, frq=frq)
# Show pipes
pipe.display()
pipe.current()
# Get the spin of interest and save it in cdp, to access it after execution of script.
cdp.myspin = return_spin(spin_inte)
# Copy the parameters from spin of interest to a fake spin to be modified.
spin.copy(spin_from=spin_inte, spin_to=fake_spin_inte)
# Returnspin
cdp.fakespin = return_spin(fake_spin_inte)
# Modify data
if spin_inte == ":52@N":
# Set reference data
cdp.fakespin.r2[r20_key] = 6.51945
cdp.fakespin.kex = 13193.82986
cdp.fakespin.kex_err = 2307.09152
phi_ex_rad2_s2 = 93499.92172
phi_ex_err_rad2_s2 = 33233.23039
scaling_rad2_s2 = frequency_to_ppm(frq=1/(2*pi), B0=cdp.spectrometer_frq_list[0], isotope='15N')**2
print scaling_rad2_s2
cdp.fakespin.phi_ex = phi_ex_rad2_s2*scaling_rad2_s2
cdp.fakespin.phi_ex_err = phi_ex_err_rad2_s2*scaling_rad2_s2
print cdp.myspin.ri_data['R1'], cdp.myspin.ri_data_err['R1'], cdp.myspin.r2[r20_key], cdp.myspin.kex, cdp.myspin.phi_ex
print cdp.fakespin.ri_data['R1'], cdp.fakespin.ri_data_err['R1'], cdp.fakespin.r2[r20_key], cdp.fakespin.kex, cdp.fakespin.phi_ex
# Calculate the offset data
theta_spin_dic, Domega_spin_dic, w_eff_spin_dic, dic_key_list = calc_rotating_frame_params(spin=cdp.myspin, spin_id=spin_inte, verbosity=10)
# Save the data in cdp to access it after execution of script.
cdp.myspin.theta_spin_dic = theta_spin_dic
cdp.myspin.w_eff_spin_dic = w_eff_spin_dic
cdp.myspin.dic_key_list = dic_key_list
############################
# First creacte back calculated R2eff data for interpolated plots.
############################
# Return the original structure for frq, offset
spin_lock_nu1 = return_spin_lock_nu1(ref_flag=False)
# Create spinBack calculate R2eff data for the set parameters.r2 keysr2keys cdp.fakespin.back_calc = []for exp_type, frq in loop_exp_frq(): r2keysoptimisation.appendback_calc_r2eff(generate_r20_key(exp_type=exp_type, frqspin=frq) )# Save the keyscdp.r2keys fakespin, spin_id= r2keys x_w_eff fake_spin_inte, spin_lock_nu1= []x_theta = []x_disp_point = []y = []y_R1_rho = []spin_lock_nu1)
# Prepare list to hold new data
spin_lock_nu1_new = []
# Loop over the structures to generate data
for ei in range(len(spin_lock_nu1)):
# Add a new dimension.
spin_lock_nu1_new.append([])
# Then loop over the spectrometer frequencies.
for mi in range(len(spin_lock_nu1[ei])):
# Add a new dimension.
spin_lock_nu1_new[ei].append([])
# Finally the offsets.
for oi in range(len(spin_lock_nu1[ei][mi])):
# Add a new dimension.
spin_lock_nu1_new[ei][mi].append([])
# No data.
if not len(spin_lock_nu1[ei][mi][oi]):
continue
# Interpolate (adding the extended amount to the end).
for di in range(num_points):
point = (di + 1) * (max(spin_lock_nu1[ei][mi][oi])+extend) / num_points
spin_lock_nu1_new[ei][mi][oi].append(point)
# Intersert field 0
#spin_lock_nu1_new[ei][mi][oi][0] = 0.0
# Convert to a numpy array.
spin_lock_nu1_new[ei][mi][oi] = array(spin_lock_nu1_new[ei][mi][oi], float64)
# Then back calculate R2eff data for the interpolated points.
cdp.myspin.back_calc = optimisation.back_calc_r2eff(spin=cdp.myspin, spin_id=spin_inte, spin_lock_nu1=spin_lock_nu1_new)
# Calculate the offset data, interpolated
theta_spin_dic_inter, Domega_spin_dic_inter, w_eff_spin_dic_inter, dic_key_list_inter = calc_rotating_frame_params(spin=cdp.myspin, spin_id=spin_inte, fields = spin_lock_nu1_new, verbosity=0)
###### Store the data before plotting
# Create a dictionary to hold data
cdp.mydic = collections.OrderedDict()
# Loop over the data structures and save to dictionaryfor dic_key exp_type, frq, offset, ei, mi, oi in cdploop_exp_frq_offset(return_indices=True): # This is not used, but could be used to get Rex.myspin.dic_key_list: R1_rho_prime = cdp.myspin.r2[cdp.r2keys[0]r20_key] #print R1_rho_prime # Get R1
R1 = cdp.myspin.ri_data['R1']
# Add to dic
if exp_type not in cdp.mydic:
cdp.mydic[exp_type] = collections.OrderedDict()
if frq_1e6 frq not in cdp.mydic[exp_type]: cdp.mydic[exp_type][frq_1e6frq] = collections.OrderedDict() if offset not in cdp.mydic[exp_type][frq_1e6frq]: cdp.mydic[exp_type][frq_1e6frq][offset] = collections.OrderedDict() # X val cdp.mydic[exp_type][frq_1e6frq][offset]['point'] = [] cdp.mydic[exp_type][frq_1e6frq][offset]['thetapoint_inter'] = [] cdp.mydic[exp_type][frq_1e6frq][offset]['w_efftheta'] = [] cdp.mydic[exp_type][frq_1e6frq][offset]['R1_rhotheta_inter'] = [] cdp.mydic[exp_type][frq_1e6frq][offset]['R1_rho_R2effw_eff'] = [] cdp.mydic[exp_type][frq_1e6frq][offset]['pointw_eff_inter'].append(float(point))= [] # Y val cdp.mydic[exp_type][frq_1e6frq][offset]['thetaR1_rho'].append(float(theta))= [] cdp.mydic[exp_type][frq_1e6frq][offset]['w_effR1_rho_err'].append(float(w_eff))= [] cdp.mydic[exp_type][frq_1e6frq][offset]['R1_rhoR1_rho_bc'].append(float(R1_rho))= [] cdp.mydic[exp_type][frq_1e6frq][offset]['R1_rho_R2effR1_rho_inter'].append(float(R1_rho_R2eff))= []
# Y val fake cdp.mydic[exp_type][frq][offset]['fake_R1_rho'] = [] # Print values in dicY2 val cdp.mydic[exp_type][frq][offset]['R1_rho_R2eff'] = [] cdp.mydic[exp_type][frq][offset]['R1_rho_R2eff_err'] = [] cdp.mydic[exp_type][frq][offset]['R1_rho_R2eff_bc'] = [] cdp.mydic[exp_type][frq][offset]['R1_rho_R2eff_inter'] = [] # Loop over the original dispersion points. for exptypepoint, frq_dic di in loop_point(exp_type=exp_type, frq=frq, offset=offset, return_indices=True): param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) # X val cdp.mydic[exp_type][frq][offset]['point'].itemsappend(point): for theta = theta_spin_dic[param_key] cdp.mydic[exp_type][frq, offset_dic in frq_dic][offset]['theta'].itemsappend(theta): for w_eff = w_eff_spin_dic[param_key] cdp.mydic[exp_type][frq][offset, val_dics in offset_dic]['w_eff'].itemsappend(w_eff): # Average resonance spin_lock_offset #print exptype, Domega_spin_dic[param_key] # Y val R1_rho = cdp.myspin.r2eff[param_key] cdp.mydic[exp_type][frq, ][offset, val_dics]['R1_rho'], val_dics.append(R1_rho) R1_rho_err = cdp.myspin.r2eff_err[param_key] cdp.mydic[exp_type][frq][offset]['R1_rho_err'].append(R1_rho_err) R1_rho_bc = cdp.myspin.r2eff_bc[param_key] cdp.mydic[exp_type][frq][offset]['thetaR1_rho_bc'].append(R1_rho_bc)
# Modify dataY val, fakew_eff_div fake_R1_rho = 10**4rem_points = 2cdp.fakespin.back_calc[ei][0][mi][oi][di]x_w_eff_mod = cdp.mydic[x/w_eff_div for x in x_w_effexp_type][:-rem_pointsfrq][offset]y_mod = y[:-rem_points'fake_R1_rho'].append(fake_R1_rho)
# Y2 val # Calc R1_rho_R2eff R1_rho_R2eff = (R1_rho - R1*cos(theta)*cos(theta)) / (sin(theta) * sin(theta)) cdp.mydic[exp_type][frq][offset]['R1_rho_R2eff'].append(R1_rho_R2eff) R1_rho_R2eff_err = (R1_rho_err - R1_err*cos(theta)*cos(theta)) / (sin(theta) * sin(theta)) cdp.mydic[exp_type][frq][offset]['R1_rho_R2eff_err'].append(R1_rho_R2eff_err) R1_rho_R2eff_bc = (R1_rho_bc - R1*cos(theta)*cos(theta)) / (sin(theta) * sin(theta)) cdp.mydic[exp_type][frq][offset]['R1_rho_R2eff_bc'].append(R1_rho_R2eff_bc) ## Loop over the new dispersion points. for di in range(len(cdp.myspin.back_calc[ei][0][mi][oi])): point = spin_lock_nu1_new[ei][mi][oi][di] param_key = return_param_key_from_data(exp_type=exp_type, frq=frq, offset=offset, point=point) # X val cdp.mydic[exp_type][frq][offset]['point_inter'].append(point) theta = theta_spin_dic_inter[param_key] cdp.mydic[exp_type][frq][offset]['theta_inter'].append(theta) w_eff = w_eff_spin_dic_inter[param_key] cdp.mydic[exp_type][frq][offset]['w_eff_inter'].append(w_eff) # Y val R1_rho = cdp.myspin.back_calc[ei][0][mi][oi][di] cdp.mydic[exp_type][frq][offset]['R1_rho_inter'].append(R1_rho) # Y2 val # Calc R1_rho_R2eff R1_rho_R2eff = (R1_rho - R1*cos(theta)*cos(theta)) / (sin(theta) * sin(theta)) cdp.mydic[exp_type][frq][offset]['R1_rho_R2eff_inter'].append(R1_rho_R2eff) #if oi == 0: #print exp_type, frq, offset, point, theta, w_eff ####### PLOT #### ## Define labels for plottingplotlabel_R1_rho_R2eff filesave_R1_rho_R2eff = 'R1_rho_R2eff'plotlabel_R1_rho filesave_R1_rho = 'R1_rho'xlabel_w_eff # For writing math in matplotlib, see# http://matplotlib.org/1.3.1/users/mathtext.html ylabel_R1_rho = r'Effective field in rotating frame R$_{1\rho}$ [%rad s $^{-1}$]'ylabel_R1_rho_R2eff = r'R$_{1\rho}$, R$_{2,eff}$ [rad.s$^{-1}$]'%(str(w_eff_div))
xlabel_theta = 'Rotating frame tilt angle [rad]'
xlabel_w_eff = r'Effective field in rotating frame [rad s$^{-1}$]'
xlabel_lock = 'Spin-lock field strength [Hz]'
plt.xlabel(xlabel_w_eff)
plt.ylabel(ylabel_R1_rho_R2eff)
plt.legend(loc='best', prop={'size':legend_size})
plt.grid(True)
#plt.ylim([0,16])#plt.xlim([0,20000*20000])plt.xlim([0,20000])
plt.title("%s \n %s as function of %s"%(spin_inte, ylabel_R1_rho_R2eff, xlabel_w_eff))
plt.xlabel(xlabel_lock)
plt.ylabel(ylabel_R1_rho)
plt.legend(loc='best', prop={'size':legend_size})
plt.grid(True)
#plt.ylim([0,16])plt.title("%s \n %s as function of %s"%(spin_inte, ylabel_R1_rho, xlabel_theta))plt.savefig("matplotlib_%s_%s_theta.png"%(spin_inte_rep, plotlabel_R1_rho) ) # Plot R1_rho as function of theta.plt.figure()for exptype, frq_dic in cdp.mydic.items(): for frq, offset_dic in frq_dic.items(): for offset, val_dics in offset_dic.items(): graphlabel = "%s_%s"%(frq, offset) plt.plot(val_dics['theta'], val_dics['R1_rho'], '-o', label=graphlabel)plt.ylabel(ylabel_R1_rho)plt.xlabel(xlabel_theta)plt.legend(loc='best')plt.grid(Truexlabel_lock)plt.ylim([0,16])pltfig3.titleset_size_inches("%s \n %s as function of %s"%(spin_inteimg_inch_x, ylabel_R1_rho, xlabel_theta)img_inch_y)plt.savefig("matplotlib_%s_%s_theta_seps_disp.png"%(spin_inte_rep, plotlabel_R1_rhofilesave_R1_rho_R2eff) )
plt.show()
== To run Bugs ? ==<source lang="bash">relax -p r1rhor2eff.py</source>Do you get an error with matplotlib about dateutil? Then see [[Matplotlib_dateutil_bug]]
== See also ==
[[Category:Matplotlib]]
[[Category:Relaxation_dispersion analysis]]