1
$\begingroup$

I've been trying to code a kinetic monte carlo simulation for a dimerization reaction, 2M ---> D. However I have a kinetic Monte Carlo code for a first order decomposition reaction of AIBN (that's the name of the chemical compound)

The code for AIBN decomposition works perfectly, but when I try to copy it to use it for the M compound, the code fails and Mathematica goes in a infinite loop saying: "Power: Infinte expression 1/0 encountered " and then it says "Infinity:Indeterminite expression o Complexinfinity encountered".

I don't know why the AIBN code works and the dimerization code for M doesn't if is exactly the same structure. Mathematica only graphs the analytical solution for M but doesn't graph the outputstream resulting from the while Loop for kinetic Monte Carlo.

Here is the code that works, for the AIBN decomposition

AIBN0 = 0.01;(*mol/L*)
R = 1.98;TK = 60 + 273.15;
kdi = 2.89*10^(15)*Exp[-31130/(R*TK)];(* s^-1*)
tf = 3600;(*%s*)
ta = {t, 0, tf};
AIBN = (AIBN0)*Exp[-kdi*t];
p1 = Plot[AIBN, ta , PlotStyle -> Red];
NAV = 2*10^6;
XAIBN = AIBN0*NAV;
kdiMC = kdi;
tb = 0.0;
AIBN = AIBN0;
Niter = 0;
file = OpenWrite["Result1.dat", FormatType -> OutputForm]
toc = Timing[
   While[tb < tf,
    Write[file, tb, " ", AIBN];
    Niter = Niter + 1;
    a0 = kdiMC*XAIBN;
    XAIBN = XAIBN - 1;
    AIBN = XAIBN/NAV;
    tb = tb - Log[RandomReal[]]/a0]];
Print[Niter, toc]
Close[file];
Result = Import["Result1.dat"];
p2 = ListPlot[Result, DataRange -> {0, 600}];
Show[p1, p2]

Here the code that doesn't work, for the M compound dimerization

Mo = 5; (*mol/L*) (*Initial concentration of M*)
R = 1.98 (*J/(mol K)*); (*Ideal Gas constant*)
Tk = 383; (*K, Temperature*)
kdim = 2.52(10^4)Exp[(-22347/(R*Tk))] (*L/(mol s)*) (*Kinetic constant*)
tf = 30000(60)(60); (*s, final reaction time, 30'000 hours*)
M = Mo/(1 + 2 kdim *t *Mo); (*Analytical solution*)
p6 = Plot[
       M, {t, 0, tf}, 
       PlotLabel -> "Concentration vs Time", 
       AxesLabel -> {"Time (s)", "Concentration (mol/L)"}, 
       PlotStyle -> {Red}, PlotLegends -> {"[M]"}]
NAV = 2*10^3; (*Avogadro number times Volume, (molecules L)/mol*)
XM = Mo*NAV; (* number of M molecules*)
kdimMC = 2*kdim/NAV; (*Monte Carlo kinetic constant*)
tb = 0.0; (*FIrst Iteration time*)
M = Mo; (*First Concentration Value of M*)
Niter = 0; (*First Iteration*)
file = OpenWrite["Result3.dat", FormatType -> OutputForm]
toc = Timing[
        While[
         tb <= tf, 
         Write[file, tb, " ", M]; 
         Niter = Niter + 1; 
         a0 = kdimMC*XM*(XM - 1)/2; 
         XM = XM - 2; 
         M = XM/NAV; 
         tb = tb - Log[RandomReal[]]/a0
        ]
      ];
Print[Niter, toc]
Close[file];
Result = Import["Result3.dat"];
p9 = ListPlot[Result, DataRange -> {0, 108000000}];
Show[p6, p9]
Improve this question
$\endgroup$
5
  • $\begingroup$ Don't do this: Write[file, tb, " ", M]; , do this instead Write[file, {tb, M}]; and remove , FormatType -> OutputForm when you open the file. $\endgroup$
    – flinty
    Aug 4, 2020 at 12:00
  • $\begingroup$ It's a bit strange writing the results to the file during the loop, and this incurs some overhead. You would be better off Sow-ing them, and then Reap-ing them to a list, then write that list out to a file afterwards. $\endgroup$
    – flinty
    Aug 4, 2020 at 12:06
  • $\begingroup$ I have already implemented the suggestions in the broken code, doing ````Write[file,{tb,M}];``` and removing ,FormatType->OutputForm when I open the file. However, now what happen is that Mathematica only shows p6 when Show[p6,p9] command is executed. If I only put Show[p9] (in order to graph only the numerical solution of resulting from the Kinetic Monte Carlo algorithm), Mathematica returns only the axes. $\endgroup$
    – Nebula
    Aug 4, 2020 at 19:31
  • $\begingroup$ I don't know how to implement Sow and Reap in my code, I'm sorry, could you explain me? Thank you $\endgroup$
    – Nebula
    Aug 4, 2020 at 19:32
  • $\begingroup$ @flinty mathematica should show a dotted graph for M, like the first code, for the AIBN, which represent a random walk $\endgroup$
    – Nebula
    Aug 4, 2020 at 20:02

1 Answer 1

Reset to default
0
$\begingroup$

Use Reap and Sow to build up a list of data, Export it all in one go, and there's no need to re-import it just to plot the results:

Mo = 5;
(*mol/L*)
(*Initial concentration of M*)

R = 1.98
(*J/(mol \K)*);
(*Ideal Gas constant*)

Tk = 383;
(*K,Temperature*)

kdim = 2.52 (10^4) Exp[(-22347/(R*Tk))]
(*L/(mol s)*)
(*Kinetic constant*)

tf = 30000 (60) (60);
(*s,final reaction time,30'000 hours*)

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

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

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

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

kdimMC = 2*kdim/NAV;
(*Monte Carlo kinetic constant*)

tb = 0.0;
(*FIrst Iteration time*)

M = Mo;
(*First Concentration Value of M*)

Niter = 0;
(*First Iteration*)

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:" , Niter, "\ntiming:", First[toc]]

p9 = ListPlot[result, DataRange -> {0, 108000000}];
Show[p6, p9]

plots

Improve this answer
$\endgroup$
1
  • $\begingroup$ it works perfectly, thank you so much!! $\endgroup$
    – Nebula
    Aug 5, 2020 at 18:18

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.

Not the answer you're looking for? Browse other questions tagged or ask your own question.