Our electrochemical Monte Carlo simulations require many evaluations to limit the noise in the simulated results and explorations over many different experimental conditions. While our initial implementation was aimed at conceptual simplicity, we will explore here how to make further modifications to achieve an additional 16x speedup in computational efficiency to unlock more ambitious simulations…

Notional Experiment

For simplicity, we will perform a simple cyclic voltammetry experiment to use as a benchmarking example, using the code from Part 2:

Clear[cv]
cv[v0_, time_Integer] := With[
    {v = N@Rest@Subdivide[v0, -v0, time/2] }, 
    Join[v, Rest@Reverse[v]]] 
 
expt = cv[10., 400]; (* demo *)

Monte Carlo Simulation, Revisited

Part 1 introduced the notion of simulating the one-dimensional diffusion problem using a discrete-time-step Monte Carlo simulation. The analyte begins in the oxidized state (M+) in solution at some random position in solution, and then takes a random walk during its trajectory. If it reaches the electrode, then it has the possibility of reacting as governed by the Boltzmann-Nernst probability. The reaction is assumed to be electrochemically reversible and occurs immediately upon impingement upon the electrode. We saw that compiling the code lead to a 10x speed improvement over the interpretted version. Part 2 expanded this slightly to specify the distance range over which to sample the initial particle position and some numerical simplifications. We shall use the Part 2 code as our starting point.

For the sake of conceptual simplicity, the function was defined to perform the time evolution of a single candidate particle. This is often a helpful way to build a Monte Carlo simulation–just focus on one trial. The final simulation is then the sum of these many events, which we performed using a ParallelSum operation. Here is the code from Part 2 to remind you, followed by a timing benchmark (output is in seconds):

sweepF = FunctionCompile@Function[ 
     {Typed[vv, TypeSpecifier["PackedArray"]["Real64", 1]], 
      Typed[maxX, TypeSpecifier["Integer16"]]}, 
     Module[
      {x =  RandomInteger[{1, maxX}],  
       state = +1, 
       obs = ConstantArray[0, Length[vv]], 
       dx = 2 RandomInteger[{0, +1}, Length[vv]] - 1}, 
      Do[
       x += dx[[t]]; 
        If[x == 0, 
         If[RandomReal[] < 1./(1. + Exp[-vv[[t]]]), 
          If[state == -1, obs[[t]] = +1; state = +1;], 
          If[state == +1, obs[[t]] = -1; state = -1;] 
         ]; 
         x = 1; 
        ]; 
       , {t, Length[vv]}]; 
      obs 
     ]] 
 
ParallelSum[sweepF[expt, 86], {i, 10^6}]; // AbsoluteTiming

05a4l8keilu0x

(*{7.60349, Null}*)

As a comment, this ParalleSum evaluation used 8 kernels on the local machine to perform the evaluation.

Make it faster by compiling the loop

While the sweepF function is compiled, ParallelSum is an interpreted function. Furthermore, parallel calculation involves overhead in initializing the parallel kernels and collecting the results, and intermediate memory usage in retaining the intermediate output vectors before summing them. A path towards a faster evaluation is to perform the entire loop inside the compiled function—essentially just to nest the prior trial inside an additional Do loop, resetting some variables (like the initial position, x, the random moves, dx, and the charge state) at the beginning of each trial, and retaining a persistent record of the observed charging processes, obs, which gets incremented or decremented during each trial to built up the final sum over the events:

sweepF2 = FunctionCompile@Function[ 
     {Typed[vv, TypeSpecifier["PackedArray"]["Real64", 1]], 
      Typed[maxX, TypeSpecifier["Integer16"]], 
      Typed[nMC, TypeSpecifier["UnsignedInteger64"]]},  (*!! specify the number of MC trials *)
     Module[
      {x, dx, state, obs, 
       maxTime = Length[vv]}, 
      obs = ConstantArray[0, maxTime];                  (*!! initialize observations *)
      Do[                                               (*!! loop over MC steps *)
       x = RandomInteger[{1, maxX}];                        (*!! reinitialize the trial *)
        dx = 2 RandomInteger[{0, +1}, maxTime] - 1; 
        state = +1; 
        Do[                                                 (*!! loop over simulation time *)
         x += dx[[t]]; 
          If[x == 0, 
           If[RandomReal[] < 1./(1. + Exp[-vv[[t]]]), 
             If[state == -1, obs[[t]]++; state = +1;],      (*!! increment, not assign *)
             If[state == +1, obs[[t]]--; state = -1;]       (*!! decrement, not assign *) 
            ]; 
           x = 1;                                           (* even if no reaction occurs, return to solution *) 
          ]; 
         , {t, maxTime}]; 
       , nMC]; 
      obs 
     ]] 
 
sweepF2[expt, 86, 10^6]; // AbsoluteTiming

03u9ri8hje7sm

(*{3.65595, Null}*)

These modifications speed up the wall time (i.e., elapsed real time) by 2x, but use only one (not 8!) core to achieve that goal, so it is 16x more efficient. We can now use ParallelMap operations over different experimental conditions to advance our electrochemistry studies.

Conclusion

The first goal in scientific computing is to write a simulation that is correct and clear. Only after that is achieved should you try to make it fast. A general pattern in writing Monte Carlo simulations is to focus on writing a single trial first, and then using parallelization to speed up the independent evaluations. But computational efficiency often beats brute force parallelism—one rabbit can pull you faster than 8 turtles.

ImageSynthesize["A race scene, drawn in the style of Beatrix Potter. One capybara (dressed in a red roman soldier outfit) is riding on a chariot pulled by **EIGHT** turtles. The reins must go to a yoke on each of the eight turtles. IMPORTANT that the reins only go to the turtles and nothing else. Another Capybara (dressed in a blue roman soldier outfit) is  riding on a saddle on one rabbit. Both are riding towards the right."]

0iarx5qpkpgmw

ToJekyll["Simulated Electrochemistry, part 5", 
  "chemistry mathematica montecarlo science teaching electrochem chiguiro"]

Monte Carlo Simulations of Electrochemistry