Initial commit

ercid_doublecol.py

+#!/usr/bin/env python3
+from rawprasslib import load_raw
+from rawautoparams import load_params
+from matplotlib import pyplot as plt
+import errno
+import numpy as np
+import os
+import time
+import scipy.optimize as spopt
+import matplotlib.pyplot as plt
+
+# SCRIPT DOES NOT CHECK FOR SANITY OF THE OUTCOME !!!
+# THAT'S YOUR JOB !!!
+
+# AE == appearance energy from LCQ obtained by sigmoid fitting
+# E_tcid == threshold collision induced dissociation energy from Armentrout
+# red dots .. residual intensity (TIC -fragment - parent)
+# black dots .. relative intensity of the fragment
+
+# Ab-initio conversion ratios derivation
+ev = 1.602176634e-19
+mol = 6.02214076e23
+farad = ev*mol
+ev2kj = farad/1000
+kj2ev = 1/ev2kj
+
+# Energy-resolved CID reference dissociation energies are taken from:
+#
+# BH_HH, BH_HOMe, BH_MeMe:
+# ========================
+# https://pubs.acs.org/doi/pdf/10.1021/acs.analchem.9b02257
+# Benzhydrylpyridinium Ions: A New Class of Thermometer Ions for the
+# Characterization of Electrospray-Ionization Mass Spectrometers
+#
+# R. Rahrt, T. Auth, M. Demireva, P. B. Armentrout, K. Koszinowski, Anal. Chem.
+# 2019, 91, 11703–11711.
+# Took from the Abstract - the values are the same as the one in the Table2
+#
+# p-H, p-Me, p-OMe, p-NO2, p-NO2_alt:
+# ===================================
+# https://link.springer.com/content/pdf/10.1007/s13361-017-1693-0.pdf
+# How Hot are Your Ions Really? A Threshold Collision-Induced Dissociation
+# Study of Substituted Benzylpyridinium “Thermometer” Ions
+#
+# J. E. Carpenter, C. P. McNary, A. Furin, A. F. Sweeney, P. B. Armentrout, J.
+# Am. Soc. Mass Spectrom. 2017, 28, 1876–1888.
+# Took from the Table2, p-NO2_alt means the dissociation of the NO2
+#
+#
+# Due to the argumentation used by Armenrout in the abovementioned paper (2017)
+# we ommited the values presented in:
+# https://www.sciencedirect.com/science/article/pii/S1387380616303608
+# Experimental bond dissociation energies of benzylpyridinium thermometer ions
+# determined by threshold-CID and RRKM modeling
+#
+# D. Gatineau, A. Memboeuf, A.  Milet, R. B. Cole, H. Dossmann, Y. Gimbert, D.
+# Lesage, International Journal of Mass Spectrometry 2017, 417, 69–75.
+
+
+# energies, uncertainty (eV) - from Armentrout experimental values
+energs = {"p-H": (2.58, 0.15),
+          "p-Me": (2.26, 0.13),
+          "p-OMe": (1.93, 0.08),
+          "p-NO2": (3.04, 0.12),
+          "p-NO2_alt": (2.68, 0.13),
+          "BH_HH": (1.79, 0.11),
+          "BH_MeMe": (1.55, 0.13),
+          "BH_HOMe": (1.37, 0.14)}
+
+# masses to select - parent, fragment
+ions = {"p-H": (170, 91),
+        "p-Me": (184, 105),
+        "p-OMe": (200, 121),
+        "p-NO2": (215, 136),
+        "p-NO2_alt": (215, 169),
+        "BH_HH": (246, 167),
+        "BH_MeMe": (274, 195),
+        "BH_HOMe": (276, 197)}
+
+# a bit generous peakwidth
+peakwidth = 1.5
+
+def boltzmann(x, bottom, top, xhalf, slope):
+    """boltzmann sigmoid curve function"""
+    y = bottom + (top - bottom) / (1 + np.exp((x-xhalf)/slope))
+    return y
+
+def fitline(x, a):
+    return x*a
+
+
+def get_tangentparams(parameters):
+    """generate a, b for y = ax + b from sigomid fit parameters"""
+    bottom, top, xhalf, slope = parameters[0]
+    a = -(top - bottom) / (4*slope)
+    yhalf = top + (bottom - top) / 2
+    b = yhalf - a * xhalf
+    return a, b
+
+
+def get_r_sqrd(function, xdata, ydata, popt):
+    """routine for getting fit probability"""
+    residuals = ydata - function(xdata, *popt)
+    ss_res = np.sum(residuals**2)
+    ss_tot = np.sum((ydata-np.mean(ydata))**2)
+    r_sqrd = 1 - (ss_res / ss_tot)
+    return r_sqrd
+
+
+figure = plt.figure(figsize=(8, 8), dpi=100, constrained_layout=True)
+grid = figure.add_gridspec(int(np.ceil(len(energs)/2))+1, 2)
+
+
+xAE, yAE, dates, yerrs = [], [], [], []
+for n,i in enumerate(energs.keys()):
+    print("processing %s..." % i)
+
+    # data readout and processing
+    filename = i+".RAW" if not i == "p-NO2_alt" else "p-NO2.RAW"
+    if not os.path.isfile(filename):
+        raise FileNotFoundError(
+            errno.ENOENT, os.strerror(errno.ENOENT), filename)
+    dates.append(os.stat(filename).st_mtime)
+    matrix = load_raw(filename)[0]
+    legends = load_params(filename)
+    normenergs = np.asarray(legends[1]).T[4]
+    times = matrix[0][0]
+    assert len(times) == len(normenergs)
+    masses = matrix[1]
+    parentargs = np.where((masses > (ions[i][0]-peakwidth/2)) &
+                          (masses < (ions[i][0]+peakwidth/2)))
+    daughterargs = np.where((masses > (ions[i][1]-peakwidth/2)) &
+                            (masses < (ions[i][1]+peakwidth/2)))
+    #tic == zero -> bad scan
+    tic = np.sum(matrix[2].T, axis=0)
+    assert len(times) == len(tic)
+    good = np.where(tic > 0)
+    parentints = np.divide(np.sum(matrix[2].T[parentargs], axis=0), tic)[good]
+    daughterints = np.divide(np.sum(matrix[2].T[daughterargs], axis=0), tic)[good]
+    rest = 1 - parentints - daughterints
+    
+    # fitting
+    x = normenergs[good]
+    y = daughterints
+    try:
+       guess = [min(y), max(y), x[np.argmax(y>(max(y)/2))], 0.1]
+       parameters = spopt.curve_fit(boltzmann, x, y, guess)
+       print(parameters[0])
+    except:
+        print("optimization failed, giving up")
+        raise Exception
+ 
+    a, b = get_tangentparams(parameters)
+    zerocross = -b/a
+    yline = [0, max(y)]
+    xline = [zerocross, (yline[1]-b)/a]
+    rsqrd = get_r_sqrd(boltzmann, x, y, parameters[0])
+    xfit = np.linspace(x[0],x[-1],500)
+    yfit = boltzmann(xfit, *parameters[0])
+
+    plot = figure.add_subplot(grid[int(n/2),n % 2], ylabel="$I_{frag}\ /\ \Sigma I$",
+                              xlim=(min(x), max(x)), ylim=(-0.1,1.1))
+    plot.spines['top'].set_visible(False)
+    plot.spines['right'].set_visible(False)
+    plot.plot(x, rest, '.', color = [1, 0.8, 0.8, 0.1])
+    plot.plot(x, daughterints, '.', color = [0,0,0,0.1])
+    plot.plot(xfit, yfit, color=[0,0,0,0.4])
+    plot.plot(xline, yline, color=[0,0,0,0.7])
+    plot.text(0.01, 1, 
+              "{}\n$AE={:.2f}$\n$E_{{tcid}}={}\ eV$\n$r^2={:.4f}$".format(
+        i, zerocross, energs[i][0], rsqrd), 
+        transform=plot.axes.transAxes, va='top')
+
+    xAE.append(zerocross)
+    yAE.append(energs[i][0])
+    yerrs.append(energs[i][1])
+    print(zerocross)
+    print("DONE! (%s)" % i)
+
+finalfit = spopt.curve_fit(fitline, xAE, yAE)
+k = finalfit[0][0]
+rsq = get_r_sqrd(fitline, np.asarray(xAE), np.asarray(yAE), finalfit[0])
+
+figure.suptitle("Calibration plots\n{} - {}".format(
+    time.ctime(min(dates)), time.ctime(max(dates))))
+
+print("Plotting the fitted values")
+plot = figure.add_subplot(grid[-1,:], xlabel="$AE$",
+                          ylabel="$E_{{tcid}}$")
+plot.errorbar(xAE, yAE, yerr=yerrs, color = [0,0,0,0.4], linestyle='None' )
+plot.plot(xAE, yAE, '.', color = 'black')
+plot.plot((min(xAE), max(xAE)), (min(xAE)*k, max(xAE)*k), color = 'red')
+for n,i in enumerate(energs.keys()):
+    plot.annotate(i, (xAE[n]+0.1, yAE[n]), rotation=0, va='top')
+plot.text(0.01, 1, "$k={:.4f}\ eV\ ({:.4f}\ kJ/mol),\ R^2 ={:.2f}$".format(
+    k, k * farad/1000, rsq), transform=plot.axes.transAxes, va='top')
+plt.savefig("doublecol.png")