%pylab inline
import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize
from scipy.optimize import curve_fit
Populating the interactive namespace from numpy and matplotlib
# Data for bulk Nash experiment, 1M C1 0.2M C2 at t=0
# Remark: separate attempts to convert C2 and other species to Nash products did not yield
# measurable absorption. We thus attribute the full absorption to C1 (formaldehyde)
# Measured intensity in A.U. after subtraction of background, every minute, 3 repeats (placed one after another)
C1MM=[1.09233333333333 ,0.728733333333333 ,0.668233333333333 ,0.476533333333333 ,0.539433333333333 ,
0.508733333333333 ,0.535533333333333 ,0.359333333333333 ,0.349933333333333 ,2.88083333333333 ,0.290233333333333 ,
0.294533333333333 ,0.251933333333333 ,0.196033333333333 ,0.185833333333333 ,0.184333333333333 ,0.180233333333333 ,
0.146333333333333 ,0.120133333333333 ,0.0991333333333333 ,0.0741333333333333 ,0.933633333333333 ,0.729933333333333 ,
0.558133333333333 ,0.595233333333333 ,0.531133333333333 ,0.454833333333333 ,0.438733333333333 ,0.393433333333333 ,
0.322933333333333 ,0.315433333333333 ,0.318033333333333 ,0.267433333333333 ,0.285433333333333 ,0.242233333333333 ,
0.230833333333333 ,0.187633333333333 ,0.184333333333333 ,0.156133333333333 ,0.114533333333333 ,0.100033333333333 ,
0.0772333333333333 ,0.945933333333333 ,0.765833333333333 ,0.679433333333333 ,0.548433333333333 ,0.524333333333333 ,
0.455633333333333 ,0.440933333333333 ,0.391533333333333 ,0.313433333333333 ,0.217033333333333 ,0.312933333333333 ,
0.312933333333333 ,0.261033333333333 ,0.192933333333333 ,0.249333333333333 ,0.198033333333333 ,0.173433333333333 ,
0.137833333333333 ,0.126433333333333 ,0.102533333333333 ,0.0731333333333333]
BB=0.895705255 #Calibration factor for formaldehyde, M / A.U.
C1M,TTEXPM=[],[] # Concentration, Time in Minutes
nreps=int(len(C1MM)/3)
for k in range(3):
for i in range(nreps):
TTEXPM.append(0.5+i)
for i in range(len(C1MM)):
C1M.append(C1MM[i]*BB) #Concentration from intensity
#remove obvious outlier at index 9
C1M = [x for i,x in enumerate(C1M) if i!=9]
TTEXPM = [x for i,x in enumerate(TTEXPM) if i!=9]
#reorder time points and add initial condition (simplifies integration and fitting)
I = numpy.argsort(TTEXPM)
TTEXPM = [0]+[TTEXPM[i] for i in I]
C1M = [1]+[C1M[i] for i in I]
If the data has loaded correctly, the following yields plots for formaldehyde at 3M and 1M, (and 0.2M C2)
#matplotlib.pyplot.subplots_adjust(left=None, bottom=None, right=2.0, top=2.0, wspace=0.3, hspace=0.5)
plot(TTEXPM,C1M,'o')#
ylabel('Concentration (M)',size=20)
xlabel('Time (min)',size=20)
Text(0.5, 0, 'Time (min)')
bla
def TOYFORM(CC,t,K): #equations of toy formose
[C1A,C2A,C3A,C4A]=CC #concentrations in mol/L
[k12,k13,k4,k22]=K
R12A=k12*C1A*C2A
R13A=k13*C1A*C3A
R42A=k4*C4A - k22*C2A*C2A #reverse reaction made explicit for comparison, note that k22 is set to 0
f=[]
f.append((- R12A - R13A)) #dN1A/dt
f.append((- R12A + 2.0*R42A)) #dN2A/dt
f.append((-R13A + R12A)) #dN3A/dt
f.append((+R13A - R42A)) #dN4A/dt
return f
def INTToy(time,kadd,kcleave,C1,C2): #time integration of toy formose
abserr = 1.0e-11
relerr = 1.0e-11
k=[kadd,kadd,kcleave,kadd]
NL0 = [C1,C2,0,0]
NLF = odeint(TOYFORM, NL0, time, args=(k,),atol=abserr, rtol=relerr) # Call the ODE solver.
N1MM=[]
for i in range(len(NL0)):
N1MM.append([])
for k in range(len(time)):
N1MM[i].append(NLF[k,i])
return N1MM[0],N1MM[1],N1MM[2],N1MM[3] # time courses of C1 C2 C3 C4
def TOYtoFIT(time,kadd,kcleave): #formats for curve_fit
return INTToy(time,kadd,kcleave,1,0.2)[0] #C1=1M C2=0.2M outputs time course of C1
Kfit=curve_fit(TOYtoFIT,TTEXPM,C1M,absolute_sigma=True)#,method='trf')
C:\Users\Philippe\anaconda3\lib\site-packages\scipy\integrate\odepack.py:247: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information. warnings.warn(warning_msg, ODEintWarning)
kadd = Kfit[0][0]#rate constant of C1 addition (M^-2.min^-1)
kcleave = Kfit[0][1]#rate constant of C4 cleavage into 2 C2 (M^-1.min^-1)
display(kadd)
display(kcleave)
display(kadd/kcleave)#ratio found 47 in the most reliable and similar (higher initial C2) estimate of Socha 1981
2.832822500149278
0.05966631190084524
47.47775436257773
N1=INTToy(TTEXPM,kadd,kcleave,1.0,0.2)[0]
N2=INTToy(TTEXPM,kadd,kcleave*2,1.0,0.2)[0]
N3=INTToy(TTEXPM,kadd,kcleave/2,1.0,0.2)[0]
N4=INTToy(TTEXPM,kadd*2,kcleave,1.0,0.2)[0]
N5=INTToy(TTEXPM,kadd/2,kcleave,1.0,0.2)[0]
N6=INTToy(TTEXPM,kcleave*187,kcleave,1.0,0.2)[0]
fig = plt.figure()
ax = plt.subplot(111)
plot(TTEXPM,C1M, 'o', label='Nash Data')
plot(TTEXPM,N1,'k',label='best fit')
plot(TTEXPM,N2,'g--',label='$k_{cleave}$ x2')
plot(TTEXPM,N3,'g--',label='$k_{cleave}$ /2')
plot(TTEXPM,N4,'m--',label='$k_{add}$ x2')
plot(TTEXPM,N5,'m--',label='$k_{add}$ /2')
plot(TTEXPM,N6,'k--',label='$k_{add}$ x4')
xlim(-0.5,21)
ylim(-0.05,1.05)
ylabel('Concentration',size=20)
xlabel('Time (min)',size=20)
xticks(size=16)
yticks(size=16)
legend(fontsize=12)
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
#plt.savefig('NashFitVariation.eps',dpi=1000,bbox_inches='tight')
<matplotlib.legend.Legend at 0x1ec45ddfeb0>
def TOYFORMCAN(CC,t,K): #equations of toy formose with Cannizzaro on C1
[k12,k13,k4,k22,kcan]=K
[C1A,C2A,C3A,C4A,meth]=CC #concentrations in mol/L
R12A=k12*C1A*C2A
R13A=k13*C1A*C3A
R42A=k4*C4A - k22*C2A*C2A #reverse reaction made explicit for comparison, note that k22 is set to 0
Rc=kcan*C1A**2 #cannizzaro reaction of C1 gives formate and methanol for fixed [OH-]
f=[]
f.append((- R12A - R13A -2*Rc)) #dN1A/dt
f.append((- R12A + 2.0*R42A)) #dN2A/dt
f.append((-R13A + R12A)) #dN3A/dt
f.append((+R13A - R42A)) #dN4A/dt
f.append((+Rc)) #d[methanol]/dt
return f
def INTToyCAN(time,kadd,kcleave,kcan,C1,C2): #time integration of toy formose
abserr = 1.0e-11
relerr = 1.0e-11
k=[kadd,kadd,kcleave,kadd,kcan]
N0methanol=0.226 #Methanol
NL0 = [C1,C2,0,0,N0methanol]
NLF = odeint(TOYFORMCAN, NL0, time, args=(k,),atol=abserr, rtol=relerr) # Call the ODE solver.
N1MM=[]
for i in range(len(NL0)):
N1MM.append([])
for k in range(len(time)):
N1MM[i].append(NLF[k,i])
return N1MM[0],N1MM[1],N1MM[2],N1MM[3],N1MM[4] # time courses of C1 C2 C3 C4 formate methanol
oh=0.44
kcan=0.0056*oh # M**-1.min**-1 from [OH-]*k4 value in Socha1981
Ncan=INTToyCAN(TTEXPM,kadd,kcleave,kcan,1.0,0.2)[0]
plot(TTEXPM,C1M, 'o', label='Nash Data')
plot(TTEXPM,N1,label='no Cannizzaro')
plot(TTEXPM,Ncan,label='Cannizzaro')
plot(0,1,'ko')
legend()
xlim(-0.5,21)
ylim(-0.05,1.05)
ylabel('Concentration',size=20)
xlabel('Time (min)',size=20)
xticks(size=16)
yticks(size=16)
legend(fontsize=12)
#plt.savefig('NashFitCannizzaro.eps',dpi=1000,bbox_inches='tight')
<matplotlib.legend.Legend at 0x1ec471c2850>
#Data for NMR experiments, 1: 3M C1, 0.2M C2 2: 1M C1, 0.2M C2 3: 1M C1
ME1 = np.genfromtxt("NMRN\\from_2.28_to_1.97_normalizedf1.txt")
ME2 = np.genfromtxt("NMRN\\from_2.34_to_2.05_normalizedf2.txt")
ME3 = np.genfromtxt("NMRN\\from_2.39_to_2.05_normalizedf3.txt")
TMG1 = np.genfromtxt("NMRN\\from_2.92_to_2.61_normalizedf1.txt")
TMG2 = np.genfromtxt("NMRN\\from_2.93_to_2.68_normalizedf2.txt")
TMG3 = np.genfromtxt("NMRN\\from_3.05_to_2.56_normalizedf3.txt")
MEOH1 = np.genfromtxt("NMRN\\from_3.27_to_3.04_normalizedf1.txt")
MEOH2 = np.genfromtxt("NMRN\\from_3.29_to_3.24_normalizedf2.txt")
MEOH3 = np.genfromtxt("NMRN\\from_3.41_to_3.12_normalizedf3.txt")
FORM1 = np.genfromtxt("NMRN\\from_8.36_to_8.27_normalizedf1.txt")
FORM2 = np.genfromtxt("NMRN\\from_8.4_to_8.35_normalizedf2.txt")
FORM3 = np.genfromtxt("NMRN\\from_8.45_to_8.35_normalizedf3.txt")
# Timepoints for NMR in hours
TT1=linspace(0,len(ME1)*17.72/3600.0,len(ME1))
TT2=linspace(0,len(ME2)*11.28/3600.0,len(ME2))
TT3=linspace(0,len(ME3)*11.28/3600.0,len(ME3))
# Timepoints for NMR in minutes
TT1M=linspace(0,len(ME1)*17.72/60.0,len(ME1))
TT2M=linspace(0,len(ME2)*11.28/60.0,len(ME2))
TT3M=linspace(0,len(ME3)*11.28/60.0,len(ME3))
Nmeth=INTToyCAN(TT2M,kadd,kcleave,kcan,1.0,0.2)[4]
plt.figure()
plot(TT2M,(MEOH2-MEOH2[1])/(0.25*TMG2[0]),label='methanol')
plot(TT2M,Nmeth-Nmeth[0],label='model')
plot(TT2M,12*(FORM2-FORM2[0])/(TMG2[0]),label='formate')#-FORM2[0]
yticks(size=15)
xticks(size=15)
ylabel('Concentration',size= 18)
xlabel('Time (min)',size=18)
title(r'1M C1, 0.2M C2',size=15)
legend(fontsize=15)
#plt.savefig('Methanol.eps',dpi=1000,bbox_inches='tight')
<matplotlib.legend.Legend at 0x1ec473af550>