Our past electrochemical simulations have focused on reporting experimental observables, specifically the amount (and direction) of current transferred during the reaction at the electrode. But a distinct advantage of simulations is that they also allow us to observe quantities that are less experimentally evident. In this post, we explore a simple modification that lets us observe the concentration profile of oxidized and reduced species as a function of time…

Coding strategy

The key idea is that we want to maintain an array of {position, time} dimensions that records the distribution of oxidized and reduced species during the course of the simulation. In the simulation each particle behaves independently, and so this distribution will be obtained by summing over multiple independent trajectories. For brevity we will refer to the oxidized and reduced states as “ions”, with charge states of +1 and -1. When oxidized molecules are dominant, a positive value will be seen, when reduced molecules are dominant a negative value, and when they are equal then we will observe zero.. This need not be the case in experiment, but makes the modifications to the code from part 5 relatively minor:

sweepF3 = FunctionCompile@Function[
    {Typed[vv, TypeSpecifier["PackedArray"]["Real64", 1]], 
     Typed[maxX, TypeSpecifier["Integer64"]],               (*!! increase size to be comparable to implied type of vv *)
     Typed[nMC, TypeSpecifier["UnsignedInteger64"]]}, 
    Module[{x, dx, state, obs, chargeDistbn, 
     maxTime = Length[vv]}, 
     obs = ConstantArray[0, maxTime];                       (*!! initialize observations *)
     chargeDistbn = ConstantArray[0, {maxX - 1, maxTime}];  (*!! initialize charge state *) 
      Do[(*!!loop over MC steps*)
       x = RandomInteger[{1, maxX}]; 
       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;], 
           If[state == +1, obs[[t]]--; state = -1;] ]; 
          x = 1;  (* return to solution regardless of reaction outcome *)]; 
         If[x < maxX, chargeDistbn [[x, t]] += state];     (*!! record ion location *) 
        , {t, maxTime}];, 
       nMC]; 
     {obs, chargeDistbn}] 
   ]

11dw20jzi3elz

A few salient points:

  • FunctionCompile demands that the arguments to ConstantArray have a consistent type (Integer64), leading us to change the type of maxX

  • The only substantive modification is to add the chargeDistbn array (initialized to zero), and then update it at every timestep with the current state of the system.

  • Our simulation allows particles to diffuse to arbitrarily large values of position (x) away from the electrode, but the chargeDistbn array only has a finite size, so we only record information within the finite region.

  • Because we have diffusion away from the farther positions (but not inward from beyond the starting boundary), we will want to restrict our attention to a narrower window of distances when plotting the charge distribution so as to represent the bulk.

Demonstration

Consider a basic cyclic-voltammetry experiment of the type introduced in part 1:

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

Then visualize the results by generating a list of plots of corresponding time points in the observed charge transfer (obs) and charge distribution (charge) arrays and animating the collection:

ListAnimate@ Table[
    GraphicsColumn[
     {ListPlot[obs, 
       Frame -> True, FrameLabel -> {"time", "current"}, 
       Epilog -> {Red, PointSize[0.025], Point[{t, obs[[t]]}]}], 
      
      ListPlot[charge[[All, t]], 
       Frame -> True, FrameLabel -> {"distance from interface", "net charge"}, 
       PlotRange -> { {0, 60}, MinMax[charge]}]}], 
    {t, 1, 399}]

We are unlikely to win an Academy Award for this video, but it is indeed interesting to see how the concentration distribution varies over time when aggregated over our particles. Perhaps this suggests a future study of electrochemical impedance.

ImageSynthesize["Cartoon style image of a capybara as a 1930s movie director.  He is pointing a movie camera at an electrochemical experiment in a laboratory."]

0qfck3rp24h1e

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

Monte Carlo Simulations of Electrochemistry