```julia
using DifferentialEquations
using Plots
```

Consider the simple reaction:

\begin{align}
    A + * &\longleftrightarrow A_*\\
    A_*   &\longleftrightarrow B + *\\
\end{align}

Both are elementary steps that occur in the liquid phase, and we will consider it in a few different solvent environments.


```julia
function batch(du, u, p, t)
    MAR = p["MAR"]
    PAR = p["PAR"]

    k10, k20, K10, K20, Ns, NS = PAR

    NA     = u[:,1]
    NAs    = u[:,2]
    NB     = u[:,3]

    NT     = NA + NB .+ NS
    
    XA   = NA./NT
    XB   = NB./NT
    XAs  = XA
    XTS1 = XA
    XTS2 = XB
    XS   = NS./NT

    thetaA = NAs./Ns
    thetastar = 1 .- thetaA
    
    #For A in solvent
    A12A = MAR[1]
    A21A = MAR[2]

    #For As in solvent
    A12As = MAR[3]
    A21As = MAR[4]
    
    #For B in solvent
    A12B = MAR[5]
    A21B = MAR[6]
    
    #For Transition State 1 in solvent
    A12TS1 = MAR[7]
    A21TS1 = MAR[8]
    
    #For Transition State 2 in solvent
    A12TS2 = MAR[9]
    A21TS2 = MAR[10]

    gammaA  = exp.(XS.^2 .*(A12A  .+ 2*(A21A  - A12A)*XA))
    gammaAs = exp.(XS.^2 .*(A12As .+ 2*(A21As - A12As)*XAs))
    gammaB  = exp.(XS.^2 .*(A12B  .+ 2*(A21B  - A12B)*XB))

    gammaTS1 = exp.(XS.^2 .*(A12TS1 .+ 2*(A21TS1 - A12TS1)*XTS1))
    gammaTS2 = exp.(XS.^2 .*(A12TS2 .+ 2*(A21TS2 - A12TS2)*XTS2))

    z1 = 1/K10*gammaAs./gammaA.*thetaA./XA./thetastar
    z2 = 1/K20*gammaB./gammaAs.*XB.*thetastar./thetaA

    z2[isnan.(z2)] .= 0
    
    r1 = k10*gammaA./gammaTS1.*XA.*thetastar.*(1 .- z1)
    r2 = k20*gammaAs./gammaTS2.*thetaA.*(1 .- z2)
   
    RA  = -r1[1]
    RB  =       r2[1]
    RAs =  r1[1] .- r2[1]
    
    du[1] = RA*Ns
    du[2] = RAs*Ns
    du[3] = RB*Ns
    return du, r1, r2, z1, z2
end
```




    batch (generic function with 1 method)




```julia
K    = 1
Ns   = 0.001
k10  = 1000000
k20  = 0.01
K10  = 1
K20  = K/K10
V    = 1
NTOT = 100

NA0  = 0.1
NAs0 = 0.0
NB0  = 0.0

NS   = NTOT - NA0 - NB0

var0     = [NA0 NAs0 NB0]
span     = (0.0, 1000000/Ns);
```


```julia
#Solvate transition state relative to reactants
MARSET1 = zeros(10,3)
MARSET1[:,1] = [0.0, 0.0, 0.0, 0.0,  0.0, 0.0,  0.0, 0.0,  0.0,  0.0] #no solvation
MARSET1[:,2] = [0.0, 0.0, 0.0, 0.0,  0.0, 0.0,  0.0, 0.0,  0.0,  0.0] #destabilize TS1
MARSET1[:,3] = [0.0, 0.0, 0.0, 0.0,  0.0, 0.0,  0.0, 0.0, 0.0, 0.0] #stabilize TS1

tfine  = 10 .^range(0.0, stop = log10(maximum(span)), length = 1000)
Xout   = zeros(length(tfine), size(MARSET1, 2))
r1out  = zeros(length(tfine), size(MARSET1, 2))
r2out  = zeros(length(tfine), size(MARSET1, 2))
z1out  = zeros(length(tfine), size(MARSET1, 2))
z2out  = zeros(length(tfine), size(MARSET1, 2))
cov    = zeros(length(tfine), size(MARSET1, 2))

for i = 1:size(MARSET1, 2)
    p0   = Dict("MAR" => MARSET1[:,i], "PAR" => [k10, k20, K10, K20, Ns, NS])
    prob = ODEProblem(batch, var0, span, p0)
    sol  = solve(prob, Rodas5(), abstol = 1e-16, reltol = 1e-16)
    solf = sol(tfine)
    NA   = solf[1,:]
    NAs  = solf[2,:]
    NB   = solf[3,:]
    NT   = NA + NB .+ NS
    
    dut, rt1, rt2, zt1, zt2 = batch([0., 0., 0.], [NA NAs NB], p0, tfine)

    ext  = NB/NA0
    cov[:, i]  = NAs/Ns
    Xout[:,i]  = ext
    r1out[:,i] = rt1
    r2out[:,i] = rt2
    z1out[:,i] = zt1
    z2out[:,i] = zt2   
end

plt1 = plot(tfine/maximum(tfine), Xout, xlabel  = "time", ylabel = "Conversion", labels = ["e1" nothing nothing], legend = :bottomright)
plt2 = plot(tfine[2:end]/maximum(tfine), r2out[2:end,:], xlabel  = "time", ylabel = "rate", labels = ["r2" nothing nothing], legend = :topright)
plt3 = plot(Xout[2:end,:], r2out[2:end,:], xlabel  = "extent", ylabel = "rate", labels = ["r1" nothing nothing], legend = :topright)
plt4 = plot(Xout, z1out, xlabel  = "extent", ylabel = "z", labels = ["z1" nothing nothing], legend = :topright)
plt4 = plot!(Xout, z2out, labels = ["z2" nothing nothing])
plt5 = plot(tfine/maximum(tfine), z1out, xlabel  = "time", ylabel = "z", labels = ["z1" nothing nothing], legend = :topright)
plt5 = plot!(tfine/maximum(tfine), z2out, labels = ["z2" nothing nothing])
plt6 = plot(tfine/maximum(tfine), cov, xlabel = "time", ylabel = "Coverage of A")


display(plt1)
display(plt2)
display(plt3)
display(plt4)
display(plt5)
display(plt6)

```


```julia
function rate(XA, p)
    MAR = p["MAR"]
    PAR = p["PAR"]

    k10, k20, K10, K20, Ns, NS = PAR
    
    XS   = 0.99
    XB   = 1.0 - XS .- XA
    XAs  = XA
    XTS1 = XA
    XTS2 = XB

    #For A in solvent
    A12A = MAR[1]
    A21A = MAR[2]

    #For As in solvent
    A12As = MAR[3]
    A21As = MAR[4]
    
    #For B in solvent
    A12B = MAR[5]
    A21B = MAR[6]
    
    #For Transition State 1 in solvent
    A12TS1 = MAR[7]
    A21TS1 = MAR[8]
    
    #For Transition State 2 in solvent
    A12TS2 = MAR[9]
    A21TS2 = MAR[10]

    gammaA  = exp.(XS.^2 .*(A12A  .+ 2*(A21A  - A12A)*XA))
    gammaAs = exp.(XS.^2 .*(A12As .+ 2*(A21As - A12As)*XAs))
    gammaB  = exp.(XS.^2 .*(A12B  .+ 2*(A21B  - A12B)*XB))

    gammaTS1 = exp.(XS.^2 .*(A12TS1 .+ 2*(A21TS1 - A12TS1)*XTS1))
    gammaTS2 = exp.(XS.^2 .*(A12TS2 .+ 2*(A21TS2 - A12TS2)*XTS2))

    thetastar = 1 ./(1 .+ K10*gammaA./gammaAs.*XA)
    thetaA    = K10*gammaA./gammaAs.*XA.*thetastar
    
    z2 = 1/K20*gammaB./gammaAs.*XB.*thetastar./thetaA

    z2[isnan.(z2)] .= 0
    z2[isinf.(z2)] .= 0
    r2 = k20*gammaAs./gammaTS2.*thetaA.*(1 .- z2)
    return r2
end

MARSET2 = zeros(10,3)
MARSET2[:,1] = [0.0, 0.0, 0.0, 0.0,  0.0, 0.0,  0.0, 0.0,  0.0,  0.0] #no solvation
MARSET2[:,2] = [0.0, 0.0, 1.0, 1.0,  0.0, 0.0,  0.0, 0.0,  0.0,  0.0] #destabilize TS1
MARSET2[:,3] = [0.0, 0.0, -1.0, -1.0,  0.0, 0.0,  0.0, 0.0,  0.0, 0.0] #stabilize TS1
K    = 1
Ns   = 0.001
k10  = 1000000
k20  = 0.01
K10  = 10000000
K20  = K/K10


XA = 0.0:0.0001:0.01
r2out = zeros(length(XA), size(MARSET2, 2))
for i = 1:size(MARSET2, 2)
    p0   = Dict("MAR" => MARSET2[:,i], "PAR" => [k10, k20, K10, K20, Ns, NS])
    r2out[:,i] = rate(XA, p0)
end

plot(XA, r2out)
```

c98271502725cfd264200ec6c7b8c50bdbceb42f
ipynb
Jupyter Notebook
2021_JCAT_DeDonder_Solvents/Case Study 3.ipynb
jqbond/Research_Public
a6eb581e4e3e72f40fd6c7e900b6f4b30311076f
MIT
2021_JCAT_DeDonder_Solvents/Case Study 3.ipynb
jqbond/Research_Public
a6eb581e4e3e72f40fd6c7e900b6f4b30311076f
MIT
2021_JCAT_DeDonder_Solvents/Case Study 3.ipynb
jqbond/Research_Public
a6eb581e4e3e72f40fd6c7e900b6f4b30311076f
MIT
Qwen/Qwen-72B
1. YES
2. YES