
%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]


#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)')


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>

Data for Nash reaction, tracking C1¶

Toy model of the formose reaction¶

Toy model with Cannizzaro reaction¶

Data NMR¶