5. Computing in Chemistry: Diffusion and Brownian Motion

Most chemical processes are too complicated to be modelled simply, but with diffusion and Brownian motion, the computer can help to explain and make vivid how molecules move. The first model here is a sort of cellular automaton – as discussed in the previous chapter – and it makes uses of the “bitwise” operator xor which was introduced in §4.4. But apart from that, it can be understood quite independently, because it works so straightforwardly.

5.1 A Model of Diffusion {#5-1-a-model-of-diffusion}

The “Diffusion” example program (pictured here at an early stage in its run) models diffusion between two liquids in a tapering tube. The taper makes things more interesting, because it affects the numbers of the different molecules that come into contact with each other along the tube. To keep the program relatively straightforward, the canvas dimensions and resolution are set such that even though the tube is visually wider than it is high, in pixels its width and height are both exactly 100 (so the pixels are flattish rectangles rather than squares). Having such big pixels affects the quality of the onscreen display, but the text is at least adequate for labelling the two graphs at the top. The upper graph shows the proportion of red molecules at each horizontal point in the tube. The other shows the number of red molecules at each point. Both graphs are continuously updated as diffusion occurs between the red and blue liquids, while two small markers continue to record the original point of division between the liquids.

Most of the program is concerned with setting up the image and handling the graphs, but only the mechanics of diffusion are of interest to us here. The relevant code is as follows:

x1 := random(width);

y1 := random(width);

if y1<=x1 then

begin

y1 := bottom-y1;

x2 := x1;

y2 := y1;

if random(2)=0 then

inc(x2)

else

dec(y2);

temp := pixcol(x1+leftaxis,y1);

if (pixcol(x2+leftaxis,y2) xor temp)=xorval then

begin

pixset(x1+leftaxis,y1,pixcol(x2+leftaxis,y2));

pixset(x2+leftaxis,y2,temp);

if (x2<>x1) then

showswitch(x1,temp)

end

end;

All this is set within a never-ending “repeat” loop to ensure that diffusion will continue indefinitely (“until 0=1”). The loop starts by assigning random values to the variables x1 and y1, in the range 0 to 99 inclusive (note that the constant width is set to 100, adjustable if desired). But if y1 is greater than x1, the loop is abandoned and starts again. (If instead we first chose x1 and then made y1 equal to a random value between 0 and x1, that would disrupt the probabilities, making the molecules at the thin end of the tube over-active compared with the rest.) We next treat y1 as a distance from the bottom of the canvas, and reassign it accordingly (with “y1 := bottom-_y1”). From now on, _x1 and y1 can together be considered as identifying some random pixel within the tube, with x1 being measured from its thin end (rather than the canvas origin), so the pixel is physically located at (x1+leftaxis, y1) in canvas coordinates.

Having identified this pixel, x2 and y_2 are made equal to _x1 and y1 respectively, then either x2 is incremented or y2 decremented (at random). We have now identified a random pair of pixels, which are adjacent either horizontally (if x2 was incremented) or vertically (if y2 was decremented). The pixel at (x2+leftaxis, y2) is not guaranteed to be within the tube, as it could be in the black border (either above or to the right). But if one of the two pixels is red and the other blue, we now swap them round and record the resulting diffusion. To achieve this, we first note the colour of the first pixel, and call this temp. Then we xor this with the colour of the second pixel and check whether the result equals xorval, which is set to #FF00FF, the value of (red xor blue). This gives a neat way of testing that one of the pixels is red and the other blue (in any order), since if either of them were black, the result would be different. Assuming that the test is passed, we swap the pixel colours. Finally, if x2 is different from x1 – that is, if the swap has been horizontal along the tube rather than vertical within it – then we graph the change through the showswap procedure, which is given both the x1 coordinate and the colour of the first pixel (so it knows where the change has taken place, and what kind of change it was). In the full example program, these molecule swapping and graph adjustment commands take place while the screen is temporarily “frozen” (using noupdate to freeze it, and update to unfreeze it). This improves the smoothness of the visual changes, but doesn’t affect the logic.

This mechanism for modelling diffusion – selecting a random pair of molecules and swapping them – might seem very crude. But in fact it does give a fairly realistic picture of the overall dynamics, and the program can help to give an intuitive understanding of this important process.

5.2 Brownian Motion {#5-2-brownian-motion}

In 1827, the Scottish botanist Robert Brown, using a microscope to examine pollen grains suspended in water, noticed that the grains seemed to move around randomly, as though being hit by some tiny invisible objects. This phenomenon – named Brownian motion in his honour – was fully explained by Albert Einstein in one of his great papers of 1905, thus providing the strongest evidence up to that time of the real existence of molecules and the atomic theory of matter.

The “BrownianMotion” example program, pictured here, does not try to model reality accurately (because to scale, the molecules of water would be invisible and their speeds enormous), but instead, aims to give a general understanding of how Brownian motion occurs. We see one blue “pollen grain”, which starts out in the centre of the canvas, and lots of red “molecules” which are given random directions and speeds (for simplicity, when they shoot off the canvas at one side, they continue back onto it with the same direction and speed at the other side). The motion of the pollen grain is tracked in a similar way to the projectile and rocket of Chapter 3, by recording its x- and y-position, and its x- and y-velocity, using integer variables:

VAR px, py: integer;

pxvel, pyvel: integer;

The same sort of tracking has to be done for the molecules, but since there are lots of them – in fact 400, given by the value of the constant molecules – it is impractical to have separately named variables for each of them. Instead, we define appropriate array variables:

mx, my, ms, md: array[1..molecules] of integer;

The array mx now contains 400 integer variables, named in turn mx[1], mx[2], mx[3], and so on up to mx[400]. So these can be used to store the x-coordinate of each of the molecules, and likewise my[1] to my[400] will store the corresponding y-coordinates. To record the velocities of the molecules, it is simpler to store their speed (in canvas units per cycle) and direction (in degrees) rather than x- and y-components, because this will enable us to use the turtle to calculate their movements. Hence we define the arrays ms and md for this purpose.

The setup procedure of the program puts the pollen grain in the centre of the canvas and initially at rest (so px and py are 500; pxvel and pyvel are 0). Then a “halo” blot is drawn around it, whose radius – called hitradius – is equal to the pollen radius plus the molecule radius:

setxy(px,py);

colour(halocolour);

blot(hitradius);

The colour halocolour has code #FFFFFE (set as a constant at the beginning of the program) and is thus indistinguishable from white, so this blot will not be visible on the canvas. Its role is to mark out prohibited places for molecules to be located, avoiding overlap between the circle marking the pollen grain, and the circle marking the molecule. So when we now run a loop to place the (centres of the) molecules, we make sure that only white pixels are chosen:

for n:=1 to molecules do

begin

repeat

mx[n] := random(1000-2*molradius)+molradius;

my[n] := random(1000-2*molradius)+molradius

until pixcol(mx[n],my[n])=white;

ms[n] := random(highspeed-slowspeed+1)+slowspeed;

md[n] := random(360);

setxy(mx[n],my[n]);

blot(2*molradius)

end;

First, random values are chosen for mx[n] and mx[y], ensuring that both are far enough away from the edges of the canvas to avoid overlap (which isn’t really necessary here – just “:= random(1000)” would be fine – but it illustrates how to do this sort of thing in your own programs). Then the corresponding pixel colour is tested, and if it’s anything other than white, new choices are made. Once a white location has been identified, ms[n] is given a random value between slowspeed and highspeed inclusive, md[n] is given a random value between 0 and 359, and a blot (still coloured halocolour) is drawn in the relevant position, to stop any future molecule overlapping with this one. (Again, avoidance of overlap isn’t really necessary here, so this blot drawing could be omitted, but the technique may be useful in other programs.) By the end of the setup procedure, non-overlapping locations will have been identified for the pollen grain and for all 400 molecules, but so far, the canvas will still appear to be completely white, with the halocolour blots indistinguishable.

Drawing of the pollen grain and molecules is done by the draw routine, which demonstrates another useful technique:

Procedure draw(positive: boolean);

Var n: integer;

Begin

if positive then

colour(molcolour)

else

colour(white);

for n:=1 to molecules do

begin

setxy(mx[n],my[n]);

blot(molradius)

end;

if positive then

colour(polcolour);

setxy(px,py);

blot(polradius)

End;

This procedure takes a “Boolean” parameter (named after the British mathematician and logician George Boole) – that is, a parameter that is either true or false. If true, then molcolour (i.e. red) is selected as the colour for the molecules, and polcolour (i.e. blue) as the colour for the pollen grain. But if the parameter is false, then both pollen grain and molecules will be drawn as white – thus they will be erased. This enables easy expression of the main loop of the program, following the usual sequence of operations as explained in §2.1: freezing the canvas, erasing the objects with draw(false), moving them, redrawing them with draw(true), and then unfreezing.

5.3 Moving and Handling Collisions {#5-3-moving-and-handling-collisions}

To move molecule number m, we take advantage of the turtle (just as we did in §2.3 earlier):

setxy(mx[m],my[m]);

direction(md[m]);

forward(ms[m]);

…

mx[m]:=(turtx+1000) mod 1000;

my[m]:=(turty+1000) mod 1000;

First, the turtle is placed in the current location of the m^th molecule using setxy, then the turtle’s direction is set to the current direction of that molecule (i.e. md[m]), and it moves forward by that molecule’s speed (i.e. ms[m]). After moving, turtx and turty will yield the new position of the turtle, so that the new molecule coordinates – the updated mx[m] and my[m] – can be set accordingly. But to do this we don’t simply make mx[m] equal to turtx and my[m] to turty, because if the movement has carried the turtle off the canvas, then we want to put the molecule back on the other side of the canvas. This is achieved by adding the canvas width or height (i.e. 1000) and taking the remainder when the result is divided by 1000 – so values between 0 and 999 remain unchanged, while -5, say, becomes 995 and 1010 becomes 10. (This technique, using mod, was introduced in §4.5.)

But there is a major complication to be taken into account, which is precisely what Brownian motion is all about: the molecule might, in its movement, collide with the pollen grain. The gap in the code above – shown as “…” – contains code to deal with such a collision, which is conceptually quite tricky, but demonstrates some useful techniques:

if hypot(turtx-px,turty-py,1)<=hitradius then

begin

while hypot(turtx-px,turty-py,1)<hitradius do

back(1);

turnxy(px-turtx,py-turty);

degturn := turtd-md[m];

md[m] := (180+(turtd+degturn)) mod 360;

impact := cos(degturn,1,ms[m]);

pxvel := pxvel+sin(turtd,1,impact);

pyvel := pyvel-cos(turtd,1,impact)

end;

A collision occurs when the distance between the centre of the molecule and the centre of the pollen grain (as calculated by Pythagoras’s theorem using the hypot function) is equal to the sum of their two radii, for which we have the convenient constant hitradius (equal to polradius+molradius). But of course the movement of the molecule, as described above, is quite likely to “overshoot”, so that its calculated position turns out to be overlapping with the pollen grain rather than just touching it. In this case, we need to “backtrack” – pulling the molecule backwards – until the two circles are just touching. Here our use of the turtle to do the molecular movement is really helpful, because we can use the back(1) command to pull it back, one step at a time, until the calculated distance between the centres is just enough to remove the overlap.

We have now identified the position where the molecule collides with the pollen grain, and our task is to work out how this will affect the direction of the molecule after it bounces off the grain, and the effect of the collision on the pollen grain itself. For simplicity, we here assume that the impact is completely “elastic” and that the speed of the pollen grain is negligible compared with that of the molecule – hence the molecule’s speed will be unaffected by the bounce, and its lines of approach and departure will be symmetrical around the line of centres:

pollen grain

line of centres

molecule

The x-distance between the centre of the molecule (at turtx) and the centre of the pollen grain (at px) is (px-_turtx_), while the y-distance is (py-_turty_). Hence the turtle can be made to point along the direction of this line using the command turnxy(px-turtx,py-turty). This adjusts turtd accordingly, so if we now set variable degturn equal to (turtd-_md[m]), this will tell us how far the molecule’s initial direction of motion _md[m] has to turn to be pointing directly along the line of centres – this turn is shown by the curved arrow in the diagram above. If we then add degturn to turtd, and reverse the direction (by adding 180° and then, if necessary, subtracting 360°), that gives us the direction in which the molecule should bounce back from the collision.

As for the impulse of the molecule on the pollen grain – that is, the resulting change of momentum – this will act down the line of centres, and its magnitude will be proportional to the molecule’s velocity down that line, which is equal to the molecule’s speed ms[m] multiplied by the cosine of the angle degturn (as shown by the curved arrow in the diagram). For the sake of simplicity, we record the velocity of the pollen grain in these same proportional units, so no conversion is required at this point (it takes place later, when the velocity is divided by the constant speedratio before being applied to the position of the pollen grain). Since the impulse acts along the direction turtd, its x- and y-components are found by multiplying respectively by sin(turtd) and cos(turtd). (The latter is negative because when turtd is zero – pointing north – this is in the negative y direction.) Then pxvel and pyvel are adjusted accordingly, reaching the end of the code above. This whole process is repeated for every molecule, so pxvel and pyvel will be changed repeatedly if there are several collisions in a single cycle. Finally, at the end of this cycle, the pollen grain is moved to reflect the overall resulting velocity, by px:=px+pxvel/speedratio and py:=py+pyvel/speedratio, before being redrawn in its new position. The constant speedratio, set at 10, can be adjusted to change the jerkiness or smoothness of the pollen grain’s movement.

5.3 Ideas for Independent Exploration: Laws of Diffusion, and the Diffusion of Disease

One interesting challenge for ambitious students would be to try developing the diffusion model of §5.1 (or something like it) to derive the laws of diffusion that were first stated by Adolf Fick in the nineteenth century.

Another promising investigation, as suggested in the previous chapter, would be introduce diffusion into the SIR model of infection discussed at §4.2. This might be expected to have a dramatic impact on the spread of disease, since infected individuals will then be free to move beyond the neighbourhood of those previously infected, thus coming into contact with a new range of susceptible individuals. It would be interesting to see how “quarantine” – putting up barriers to movement in various places, depending on the severity of an outbreak – can influence the spread of infection, and to explore the interplay of this with inoculation as a strategy for controlling the spread and duration of epidemics.