0
$\begingroup$

I've been trying to code a kinetic monte carlo simulation for a system of consecutive chemical reactions using a fixed (repaired) code that the user flinty gives me in the last questions I post: Strange error not found in a monte carlo code, Infinite expression 1/0 encountered. In that post, the original code for one chemical reaction fail because of the use of Import and FormatType function, and it was fixed using Reap and Sow. Now I have a code that works perfectly using again Import and FormatType but it fails when I try to implement Reap and Sow on it.

The aim is to obtain the concentration of every chemical compound involved in a system of 2 chemical reactions.

The Reason why I don't want to use Import and FormatType is because when I try to add mor chemical reactions it fails, I mean, If you change something of that code, it gets broken also.

I don't undesrtand whats wrong if the Reap and Sow structure is similar to the code for one chemical reaction

Here the Code that works using Reap and Sow for One chemical Reaction, Dimerization, 2M-->D

Mo = 5; (*mol/L*) (*Initial concentration of M, Monomer*)
DDo = 0; (*mol/L*) (*Initial concentration of D, Dimer*)
R = 1.98 ;(*J/(mol K), Ideal Gas constant*)
Tk = 383; (*K, Tmperature of the reaction in kelvin*)
kdim = 2.52*10^4*E^(-22347/(R*Tk)); (*L/(mol s), kinetic constant*)
tf = 30000*60*60; (*s, final time in seconds, 30'000 hours of reaction*)

M = Mo/(1 + 2*kdim*t*Mo); (*Analytical solution for the ODE system, for M*)

DD = (2*kdim*t* Mo^2 + DDo + 2*kdim*t* Mo*DDo)/(1 + 2*kdim*t*Mo); (*Analytical solution for the ODE system, for D*)

p8 = Plot[{M, DD}, {t, 0, tf}, PlotLabel -> "Concentration vs Time", 
  AxesLabel -> {"Time (s)", "Concentration mol/L"}, PlotStyle -> {Red, Green}, PlotLegends -> {"[M]", "[D]"}];

(*Kinetic Monte Carlo loop starts here*)

Mo = 5;
R = 1.98;
Tk = 383;
kdim = 2.52*10^4*Exp[-22347/(R*Tk)];

tf = 30000*60*60

M = Mo/(1 + 2*kdim*t*Mo);

p6 = Plot[M, {t, 0, tf}, PlotStyle -> Red, 
   AxesLabel -> {"time (s)", "Concentration (mol/L)"}, 
   PlotLegends -> {"[M]"}, PlotLabel -> "Concentration of M"];

NAV = 2*10^3; (*Avogadro number times system volume*)

XM = Mo*NAV; (*Number of M molecules*)

kdimMC = 2*(kdim/NAV); (*Kinetic Monte carlo constant*)

tb = 0.0;

M = Mo;

Niter = 0;

toc = Timing[
   result = 
    Reap[While[tb <= tf, Sow[{tb, M}]; Niter = Niter + 1; 
       a0 = kdimMC*XM*(XM - 1)/2; XM = XM - 2; M = XM/NAV; 
       tb = tb - Log[RandomReal[]]/a0]][[2, 1]]];

Export["Result3.dat", result];

Print["iterations for M:", Niter, "\ntiming for M", First[toc]]

p9 = ListPlot[result, DataRange -> {0, tf}, PlotStyle -> Blue, 
   PlotLegends -> {"kMC [M]"}];

Mo = 5;
DDo = 0;
R = 1.98;
Tk = 383;

kdim = 2.52*10^4*Exp[-22347/(R*Tk)];

tf = 30000*60*60

M = Mo/(1 + 2*kdim*t*Mo);

DD = (2*kdim*t*Mo^2 + DDo + 2*kdim*t*Mo*DDo)/(1 + 2*kdim*t*Mo);

p7 = Plot[DD, {t, 0, tf}, PlotStyle -> Green, 
   AxesLabel -> {"time (s)", "Concentration (mol/L)"}, 
   PlotLegends -> {"[D]"}, PlotLabel -> "Concentration of D"];

NAV = 2*10^3; (*Avogadro number times system volume*)

XM = Mo*NAV; (*Initial number of molecules of M*)

XDD = DDo*NAV/2; (*It's divide by 2 because there must be half of the number of molecules of D regarding the number of M molecules, due to D is a Dimer of M*)

kdimMC = 2*(kdim/NAV); (*kinetic monte carlo constant*)

tb = 0.0;

M = Mo;

DD = DDo;

Niter = 0;

toc = Timing[result = Reap[While[tb <= tf, Sow[{tb, DD}];
       Niter = Niter + 1;
       a0 = kdimMC*XM*(XM - 1)/2;
       XM = XM - 2;
       XDD = XDD + 1;
       M = XM/NAV;
       DD = Mo - M;(*by stoichiometry, we can rewrite D concentration like this, just a mass balance*)
       tb = tb - Log[RandomReal[]]/a0]][[2, 1]]];

Export["Result3.dat", result];

Print["iterations for D:", Niter, "\ntiming for D", First[toc]]

p10 = ListPlot[result, DataRange -> {0, tf}, PlotStyle -> Orange, 
   PlotLegends -> {"kMC [D]"}];

Show[p8, p9, p10, PlotRange -> {{0, tf}, {0, 5}}]

Here the code that works using Import and FormatType for the consecutive reactions A-->B; B-->C

A0 = 1.00;(*mol/L, initial concentration of A*)
B0 = 0.00;(*mol/L, initial concentration of B*)
C0 = 0.00;(*mol/L, Initial concentration of C*)

k1 = 2*10^(-4);(* s^-1, kinetic constant for the reaction A-->B*)
k2 = 1*10^(-4);(* s^-1, kinetic constant for the reaction B-->C*)

tf = 3600*5;(*s, final time*)

ta = {t, 0, tf}; (*Interval of time*)

A = A0*Exp[-k1*t]; (*Analytical solution for A*)

B = B0*Exp[-k2*t] + (A0*k1/(k2 - k1))*(Exp[-k1*t] - Exp[-k2*t]); (*Analytical solution for B*)

Cs = A0 + B0 + C0 - A - B; (*Analytical solution for C*)

p1 = Plot[A, ta ,  PlotStyle -> {Dashed, Red}];
p2 = Plot[B, ta , PlotStyle -> {Dashed, Blue}];
p3 = Plot[Cs, ta , PlotStyle -> {Dashed, Gray}];

NAV = 2*10^3; (*Avogadro number times volume system*)

XA = A0*NAV; (*Number of molecules of A*)
XB = B0*NAV; (*Number of molecules of B*)
XC = C0*NAV; (*Number of molecules of C*)

k1MC = k1; (*Kinetic monte carlo constant*)
k2MC = k2; (*Kinetic Monte Carlo constant*)

Maxr = 2;  (*Maximum number of reactions*)

kmax = Ceiling[Log2[Maxr]]; (*Number of Steps, or movements, needed to occur one of the two possible reactions*)

VR = Table[0.0, {i, 1, 1 kmax + 1}, {j, 1, 2^kmax}]; (*Vector of Reactions*)

tb = 0.0;
Niter = 0;

file = OpenWrite["Result.dat", FormatType -> OutputForm]
toc = Timing[
   While[tb < tf,
    Niter = Niter + 1;
    
    
    VR[[kmax + 1, 1]] = k1MC*XA; (*Reaction 1*)
    
    VR[[kmax + 1, 2]] = k2MC*XB; (*Reaction 2*)

    For[i = 1, i <= kmax + 1, i++,
     For[j = 1, j <= 2^(kmax - i), j++,
       VR[[kmax + 1 - i, j]] = 
         VR[[kmax + 2 - i, 2*j]] + VR[[kmax + 2 - i, 2*j - 1]];
       ];];
    
    a0 = VR[[1, 1]]; (*Total reaction rate*)
    ra0 = RandomReal[]*a0;
    mv = 1;
    
    For[i = 1, i <= kmax, i++,
     If[ra0 <= VR[[i + 1, 2*mv - 1]], mv = 2*mv - 1, 
      ra0 = ra0 - VR[[i + 1, 2*mv - 1]]; mv = 2*mv]];
    
    If [mv == 1, XA = XA - 1; XB = XB + 1;];
    If [mv == 2, XC = XC + 1; XB = XB - 1;];

    A = XA/NAV;
    B = XB/NAV;
    Cs = XC/NAV;

    Write[file, tb, " ", A, " ", B, " ", Cs];
    tb = tb - Log[RandomReal[]]/a0]];

Close[file];
Print[Niter];
Result = Import["Result.dat"];

tb = Result[[All, 1]];
A = Result[[All, 2]];
B = Result[[All, 3]];
Cs = Result[[All, 4]];

p4 = ListPlot[Transpose@{tb, A}, PlotStyle -> {Dotted, Red}];
p5 = ListPlot[Transpose@{tb, B}, PlotStyle -> {Dotted, Blue}];
p6 = ListPlot[Transpose@{tb, Cs}, PlotStyle -> {Dotted, Gray}];
Show[p1, p2, p3, p4, p5, p6, PlotRange -> {0, 1}]

**Here the Broken Code for the system A-->B; B-->C using Reap and Sow; Updated version: I removed Write [tb," ", A, " ", B, " ", Cs]; and Close[file]; Now doesn't graph the random walk Updated Version 2: I have implemented the suggestions of flenty, Now the random walk doesn't fit to the analytical curve Updated 3: Sorry I put B=A-Ao-CS at the line before tb=tb-Log, and the correct line should be B=XB/NAV, that gives the graph at the end **

A0 = 1.00; (*mol/L*)
B0 = 0.00;(*mol/L*)
C0 = 0.00;(*mol/L*)

k1 = 2*10^-4; (*s^-1*)
k2 = 1*10^-4; (*s^-1*)

tf = 3600*5; (*s*)
ta = {t, 0, tf};
A = A0*Exp[-k1*t];
B = B0*E^(-k2*t) + (A0*k1)/(k2 - k1)*(E^(-k1*t) - E^(-k2*t));
Cs = A0 + B0 + C0 - A - B;

p1 = Plot[A, ta, PlotStyle -> {Dashed, Red}];
p2 = Plot[B, ta, PlotStyle -> {Dashed, Blue}];
p3 = Plot[Cs, ta, PlotStyle -> {Dashed, Gray}];

NAV = 2*10^3; (* (molecules L / mol) Avogadro number times Volume*)

XA = A0*NAV;  (* Number of A molecules *)
XB = B0*NAV;   (* Number of B molecules *)
XC = C0*NAV;   (* Number of C momlecules *)

k1MC = k1; (*Kinetic monte carlo reaction constants*)
k2MC = k2;

Maxr = 2;  (* Number of reactions *)

kmax = Ceiling[
   Log2[Maxr]];   (*Number of steps (or movements) to determine which \
chemical reactions proceeds, according to the binary tree*)

VR = Table[0.0, {i, 1, kmax + 1}, {j, 1, 2^kmax}];

tb = 0.0;

Niter = 0;

ClearAll[A, B, Cs]; A = A0; B = B0; Cs = A0 + B0 + C0 - A - B;

toc = Timing[result = Reap[
      While[tb < tf, Sow[{tb, A, B, Cs}];
       Niter = Niter + 1;
       
       VR[[kmax + 1, 1]] = k1MC*XA; (*R1, reaction 1, 
       VR is Vector of reactions*)
       
       VR[[kmax + 1, 2]] = k2MC*XB; (*R2, reaction 2*)
       
       For[i = 1, i <= kmax + 1, i++,
        For[j = 1, j <= 2^(kmax - i), j++,
          
          VR[[kmax + 1 - i, j]] = 
            VR[[kmax + 2 - i, 2*j]] + VR[[kmax + 2 - i, 2*j - 1]];
          ];];
       a0 = VR[[1, 1]]; (*Total vector of reactionr, VR*)
       ra0 = RandomReal[]*a0;
       mv = 1;
       
       For[i = 1, i <= kmax, i++,
        If[ra0 <= VR[[i + 1, 2*mv - 1]], mv = 2*mv - 1,
         ra0 = ra0 - VR[[i + 1, 2*mv - 1]]; mv = 2*mv]];
       
       If[mv == 1, XA = XA - 1; XB = XB + 1;];
       If[mv == 2, XC = XC + 1; XB = XB - 1;];
       
       A = XA/NAV; (*Concentration of A*)
       B = XB/NAV;
       Cs = XC/NAV;
       
       tb = tb - Log[RandomReal[]]/a0]][[2, 1]]];

Export["Result3.dat", result];
Print["iterations: ", Niter, "\ntiming ", First[toc]]

p4 = ListPlot[Transpose[result], DataRange -> {0, tf}];

Show[p1, p2, p3, p4, PlotRange -> {0, 1}]

enter image description here

Improve this question
$\endgroup$
4
  • $\begingroup$ Sow[{tb, , A, B, Cs}] you have a Null in your list i.e two commas {... , , ...} Also if you take a look at result you have unevaluated t. $\endgroup$
    – flinty
    Aug 16, 2020 at 22:42
  • $\begingroup$ Quick fix: add this ClearAll[A, B, Cs]; A = A0; B = B0; Cs = A0 + B0 + C0 - A - B; on the line just before toc = and when you do the list plot, instead do ListPlot[Transpose[result], DataRange -> {0, tf}]; $\endgroup$
    – flinty
    Aug 16, 2020 at 22:51
  • $\begingroup$ I have already removed the Null in Sow, and added the Quick fix on the line just before toc = and theListPlot[Transpose[result],DataRange -> {0,tf}]; However, now mathematica shows the random walk, but this random walk doesn't fit to the analytical solution curve just like the first two codes does, @flinty, I don't understand what are you saying about t because the variable of time is tb and the final time is tf $\endgroup$
    – Nebula
    Aug 17, 2020 at 0:28
  • $\begingroup$ The reason you had to Clear A,B,Cs is because your old definitions of A,B,Cs had a t symbol in them. This propagated through the code and you were trying to ListPlot them with unevaluated t. ListPlot can only take numerical values - that's why I recommended clearing the old definitions before the loop body. I'm afraid I don't know enough chemistry / physics so I can't help with the other problems. $\endgroup$
    – flinty
    Aug 17, 2020 at 13:04

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.