Hendrickson-Bray motional narrowing analysis processing pipeline¶
This script processes data by: 1. Reading Bruker format NMR data files 2. Fitting peaks and analyzing lineshapes 3. Calculating widths and errors normalized to Larmor frequency 4. Plotting calculated widths vs temperature 5. Fitting Hendrickson-Bray model to data 6. Exporting the fitted parameters to a CSV file 7. Plotting the data and fitted model
Functions:
get_filepaths(directory): Lists and sorts files in directory by integer value
get_nmr_paths(directory, filepaths): Constructs full paths to NMR data files
get_larmor_freq(path): Extracts Larmor frequency from Bruker data dictionary
process_nmr_data(nmr_paths, processor): Main processing pipeline that:
- Loads NMR data
- Selects region of interest (-500 to 500 ppm)
- Normalizes data
- Fits peaks with given initial parameters -- this is an important step,
as the initial parameters are crucial for the fitting process. If you have multiple peaks, that is another story.
Anyway, ensure that the first initial parameters are that of the peak that correlates to the atomic site of interest.
- Plots and saves results
- Calculates widths normalized to Larmor frequency
def get_filepaths(directory):
"""Sort directory contents numerically"""
filepaths = os.listdir(directory)
filepaths.sort(key=lambda x: int(x))
return filepaths
def get_nmr_paths(directory, filepaths):
nmr_paths = []
for filename in filepaths:
f = os.path.join(directory, filename)
nmr_paths.append(os.path.join(f, 'pdata', '1'))
return nmr_paths
def get_larmor_freq(path):
dic, data = ng.bruker.read(path)
udic = ng.bruker.guess_udic(dic, data)
larmor_freq = udic[0]['obs']
return larmor_freq
def process_nmr_data(nmr_paths, processor):
widths = []
width_errors = []
for i, nmr_path in enumerate(nmr_paths):
processor.load_data(nmr_path)
x_region, y_region = processor.select_region(-500, 500)
x_data, y_data, y_amp, y_ground = processor.normalize_data(x_region, y_region)
initial_params = [-10, 1, 110, 0.5, 0.0]
fixed_x0 = [False]
fixed_amp = [False]
fixed_width = [False]
# fixed_eta = [False] fits eta freely (default); set to [True] to fix it
# at its initial_params value instead
fixed_eta = [False]
popt, metrics, fitted = processor.fit_peaks(x_data, y_data, initial_params, fixed_x0, fixed_amp, fixed_width, fixed_eta)
fig, axes, components = processor.plot_results(x_data, y_data, fitted, popt)
legend = fig.gca().get_legend()
legend.set_frame_on(False)
fig.set_size_inches(5, 4)
path_components = nmr_path.split(os.sep)
exp_folder = path_components[path_components.index('pdata') - 1]
fig.suptitle(f'nmr folder: {exp_folder}', x=0.75, y=0.95)
processor.save_results(nmr_path, x_data, y_data, fitted, metrics, popt, components)
larmor_freq = get_larmor_freq(nmr_path)
width = metrics[0]['width'][0] * larmor_freq / 1000
width_error = metrics[0]['width'][1] * larmor_freq / 1000
widths.append(width)
width_errors.append(width_error)
return widths, width_errors
if __name__ == "__main__":
directory = "../../data/single_peak_HB_analysis/lineshape"
processor = NMRProcessor()
filepaths = get_filepaths(directory)
nmr_paths = get_nmr_paths(directory, filepaths)
widths, width_errors = process_nmr_data(nmr_paths, processor)
x data, which is the temperature in Kelvin, can be created using linspace function from numpy. Where the first argument is the starting temperature in Celsius, the second argument is the ending temperature in Celsius, and the third argument is the number of points to be generated between the starting and ending temperatures (which is the length of the widths/filepaths).
Also, the Temperature in Kelvin can be pre-stored in a csv file such that each filepath in filepaths corresponds to an NMR experimental number and the corresponding temperature in Kelvin. So you have to create a csv file (named: temperature_points.csv) with two columns: column 1 is the filepath and column 2 is the temperature in Kelvin. The csv file should be stored in the same directory as the lineshape folder
def HB_equation (x_data, B, Ea, A, D):
k = 8.62*10**-5
return (A / (1 + ((A / B) - 1) * np.exp(-Ea / (k * x_data)))) + D
def HB_analysis (x_data, y_data, y_error, initial_guesses, file_name):
#initial guesses should be written as [B_guess, Ea_guess, A_guess, D_guess]
#if u guess wrong, the fit would misbehave. try another guess values if such happens
#x_data should be in Kelvin #y_data should be in kHz : that is the imported data
k = 8.61733034e-5 #to be expressed in eV
#now, we are defining the HB equation
#now we have to run the fitting
popt, pcov = curve_fit(HB_equation, x_data, y_data, p0=initial_guesses, maxfev=5000)
#storing the optimised parameters in a tag
B_opt, Ea_opt, A_opt, D_opt = popt
#the diagonal of pcov contains the error of each parameter
perr = np.sqrt(np.diag(pcov))
#store the respective errors in a tag
B_err, Ea_err, A_err, D_err = perr
# determining the quality of the fit
##predict the model, the HB model, using the optimized parameters
#x_fit = np.linspace(np.min(x_data), np.max(x_data), 1000) #would affect determination of r2
y_fit = HB_equation(x_data, B_opt, Ea_opt, A_opt, D_opt)
##predicting the r square
squaredDiffs = np.square(y_data - y_fit)
squaredDiffsFromMean = np.square(y_data - np.mean(y_data))
rsquared = 1 - np.sum(squaredDiffs)/np.sum(squaredDiffsFromMean)
print(f"R² = {rsquared}")
#Plotting the results
#_____________________________________________________________________________#
# inspect the parameters
print(f"w = [{A_opt} / 1 + ({A_opt}/{B_opt} - 1) * e^(-{Ea_opt}/k*T)] + {D_opt}")
#print the errors
perr = np.sqrt(np.diag(pcov))
#store the respective errors in a tag
B_err, Ea_err, A_err, D_err = perr
#print the optimised variables and their errors
print('B:', f'{B_opt} +/- {B_err}')
print('Ea:', f'{Ea_opt} +/- {Ea_err}')
print('A:', f'{A_opt} +/- {A_err}')
print('D:', f'{D_opt} +/- {D_err}')
params_data = pd.DataFrame({'Parameters': ['B', 'Ea', 'A', 'D'],
'Values': [B_opt, Ea_opt, A_opt, D_opt], 'Errors': [B_err, Ea_err, A_err, D_err]})
params_data.to_csv(file_name+'HB_parameters.csv')
return y_fit, Ea_opt, Ea_err, B_opt, A_opt, D_opt
If the visualisation of the fit and optimsed parameters are satisfactory, you can save the lineshape_data dataframe to a csv file just in case you would like to visualise the data in your preferred software