Ionization, chemisorption and surface migration: energy calculations

In Maths Lab 4, you studied the phenomenon of chemisorption, focusing on forces and equilibrium. This week, we extend and deepen these ideas, using the mathematical concepts of integration and scalar product to bring in considerations of work and energy.

The preliminary code for this work is virtually identical to that for Lab 4; all we've done is add the function lattice2D, which allows you to generate 2D plots of planar lattices of a given size. Because there's nothing new here, we've put all the code in a closed cell group. If you want to examine it, by all means open this group (double-click on the "hooked" cell bracket on the right). Otherwise, you should simply select this "hooked" bracket and "shift-enter".

DON'T OMIT THIS STEP, however, or nothing will work!

Preliminary Code (select and "shift-enter")

<< Graphics`Arrow3D`

(Subscript[ϵ, 0]) =  8.85419*^-12;

e  =  1.60219*^-19;

rbond =  2.5*^-10;

rlat = 2.815*^-10;

color[chloride] = RGBColor[0, 0.8, 0];

color[sodium] = RGBColor[0.8, 0, 0];

radius[chloride] = 0.05;

radius[sodium] = 0.04;

particlePoint[{particle_, coords_}] := {color[particle], 
        PointSize[radius[particle]], Point[coords]}

sodiumLattice[size_] := 
      Flatten[Table[{sodium, {i rlat, j rlat, 0}}, {i, -size, size}, {
            j, Mod[i + size, 2] - size, size - Mod[i + size, 2], 2}], 1]

chlorideLattice[size_] := 
      Flatten[Table[{chloride, {i rlat, j rlat, 0}}, {i, -size, size}, {
            j, Mod[i + 1 + size, 2] - size, size - Mod[i + 1 + size, 2], 2}], 
        1]

lattice[size_] := Join[sodiumLattice[size], chlorideLattice[size]]

sodiumLattice2D[size_] := 
      Flatten[Table[{sodium, {i rlat, j rlat}}, {i, -size, size}, {j, 
            Mod[i + size, 2] - size, size - Mod[i + size, 2], 2}], 1]

chlorideLattice2D[size_] := 
      Flatten[Table[{chloride, {i rlat, j rlat}}, {i, -size, size}, {j, 
            Mod[i + 1 + size, 2] - size, size - Mod[i + 1 + size, 2], 2}], 1]

lattice2D[size_] := 
      Join[sodiumLattice2D[size], chlorideLattice2D[size]]

mag[vector_] := Sqrt[vector . vector]

(Subscript[F, like])[r_List]  := 
      -(((e^2) r)/(4 π (Subscript[ϵ, 0 ]) 
              mag[r]^3))

(Subscript[F, opp])[r_List]  := 
      (((e^2) r)/(4 π (Subscript[ϵ, 0]))) 
        ((1/mag)[r]^3 - ((((rbond))^8)/mag)[r]^11)

forceBetween[{particle1_, coords1_}, {particle2_, coords2_}] := 
      forceFunction[particle1, particle2][coords2 - coords1]

forceFunction[chloride, chloride] = (Subscript[F, like]);

forceFunction[sodium, sodium] = (Subscript[F, like]);

forceFunction[chloride, sodium] = (Subscript[F, opp]);

forceFunction[sodium, chloride] = (Subscript[F, opp]);

The task

Dissociation energy

In one dimension, calculations of work are relatively straightforward. If we move a particle in the direction of (say) the z-axis, and the particle is subject to a force F(z), also in the direction of the z-axis, then the work done in moving the particle from z = (Subscript[z, 1]) to z = (Subscript[z, 2]) is given by

(Subscript[(Subscript[(∫_z), 1]^z), 2])F(z) \[DifferentialD]z .

We can use this to calculate the work done in moving a chloride ion away from a sodium ion at the origin: from the equilibrium bond separation to twice that distance, say. First, let's calculate the force between the sodium ion and the chloride ion when it's at (0, 0, z), as a vector in terms of z:

vectorforce = 
      forceBetween[{chloride, {0, 0, z}}, {sodium, {0, 0, 0}}]

Now, we pick out the third (that is, the z-) component:

zforce = vectorforce[[3]]

A plot of this force might be helpful:

Plot[Evaluate[zforce], {z, 0.8  rbond, 3  rbond}]

And finally, the integration.

Integrate[zforce, {z, rbond, 2  rbond}]

(i) Using these ideas, calcuate the dissociation energy of the sodium chloride dipole: the work done in removing one of the ions to infinity. Take a note of its magnitude. Is it positive or negative? What is the significance of its sign?

Integrate[zforce, {z, rbond, ∞}]

(ii) For later work we're going to need the numerical NIntegrate command, because our integrals will become too complicated to do quickly using symbolic methods. Try

NIntegrate[zforce, {z, rbond, 10  rbond}]

NIntegrate[zforce, {z, rbond, 100  rbond}]

NIntegrate[zforce, {z, rbond, 1000  rbond}]

etc. How far out do you need to go before you obtain the dissociation energy to three significant figures?

Adsorption energy

We can do similar calculations to obtain estimates for the adsorption energy: the work done in removing an adsorbed ion from its equilibrium separation to infinity. Here's how to do it for a small lattice (remember, lattice[2] is 2 by 2, measured in unit cells). First, we define a function that specifies to contribution to the resultant force from each ion in the lattice:

forceFrom[ion_] := forceBetween[{chloride, {0, 0, z}}, ion]

Then, we calculate these forces:

forces = Map[forceFrom, lattice[2]]

Now, we add them up...

vectorforce = Apply[Plus,  forces]

... and take the z-component:

zforce = Apply[Plus, forces][[3]]

We use this to calculate the equilibrium separation (as in Maths Lab 4):

zadsorp = 
      z /. FindRoot[zforce == 0, {z, rbond}, AccuracyGoal -> 15]

Plot...

Plot[Evaluate[zforce], {z, 0.8  zadsorp, 3  zadsorp}]

... and calculate:

NIntegrate[Evaluate[zforce], {z, zadsorp, 10000  zadsorp}]

(i) Repeat the above calculation, but use a larger lattice: say lattice[5]. Important: for larger lattices, the expressions get very big, and it will be necessary to put semicolons after your assignments to forces, vectorforce and zforce in order to keep the size of your document, and your execution time, manageable. How does your equilibrium distance compare? What about your value for the adsorption energy?

(ii) Make your lattice larger still: look at lattice[10], say, and lattice[20]. How rapidly does the equilibrium distance appear to converge to a limit as the lattice size increases? What about the adsorption energy? Can you explain any discrepancy in these rates of convergence?

Surface migration

In more than one dimension, calculations of work become slightly more involved, and vector methods come into the picture. If we move a particle along a (possibly curved) path in xyz-space, and it is subject to a variable vector force F, then the work done in moving a small distance δs along that path is

F.n δs,

where n is the unit vector (that is, the vector of length 1) in the direction of travel, and the "dot" symbol represents the scalar product.

Now, by Pythagoras, ((δs)^2)=((δx)^2)+((δy)^2)+((δz)^2), where δx, δy and δz are the small changes in the coordinates along the small section of path. It follows that the work done is

F.n (Sqrt[δ(x^2)+δ(y^2)+δ(z^2)])

But think about n for a moment. It's the unit vector in the direction of travel, which means, for short paths, that

n= (δx,δy,δz)/(Sqrt[δ(x^2)+δ(y^2)+δ(z^2)]).

So the work done is actually

F . (δx,δy,δz),

or, alternatively

F . (1,δy/δx,δz/δx)δx.

The total work done along the entire path is therefore

∫ F . (1,dy/dx,dz/dx)dx.

We'll restrict ourselves to motion in a horizontal plane, so for us the formula will be

∫ F . (1,dy/dx,0)dx.

Let's apply this analysis to the problem of moving a chloride ion, at a distance of (Subscript[z, 0]) = 2.95725×(10^(-10)) m (the adsorption equilibrium distance) from the lattice plane, along a straight line from the point A  (-3 (Subscript[r, lat]), -3 (Subscript[r, lat]), (Subscript[z, 0])) to B  (3  (Subscript[r, lat]), 0, (Subscript[z, 0])).  (Recall that (Subscript[r, lat]) is the interionic separation in the sodium chloride lattice). This path has equation

y = (1/2) (x - 3  (Subscript[r, lat])), z = (Subscript[z, 0]).

Let's use lattice[5]:

latticeFigure = 
      Show[Graphics[Map[particlePoint, lattice2D[5]]], 
        AspectRatio -> Automatic, Axes -> True]

First, we set up our constant z-coordinate...

z0 = 2.95725×(10^(-10))

... and specify the contribution to the force from each ion in the lattice when the chloride ion is at the point (x, y, (Subscript[z, 0])):

forceFrom[ion_] := forceBetween[{chloride, {x, y, z0}}, ion]

Calculate these forces...

forces = Map[forceFrom, lattice[5]];

... and add them up.

vectorforce = Apply[Plus,  forces];

Now let's define our path equation...

pathEquation = ((x - 3  rlat)/2)

... plot it...

pathFigure = 
      Plot[pathEquation, {x, -3 rlat, 3  rlat}, 
        PlotStyle -> {Thickness[0.02], RGBColor[0, 0, 1]}, 
        AspectRatio -> Automatic]

... and superimpose this plot on our lattice graphic.

Show[{latticeFigure, pathFigure}]

This is what we have to integrate, namely F . (1,dy/dx,0).

pathIntegrand = 
      vectorforce . {1, D[pathEquation, x], 0} /. y -> pathEquation;

and this is the work done along the path:

pathWork = 
      NIntegrate[Evaluate[pathIntegrand], {x, -3 rlat, 3  rlat}]

(i) Show (using Mathematica if you like) that the parabolic path

\(TraditionalForm
    \`y = \(1\/\(12  r\_lat\)\) \((x + 3\ r\_lat)\)\^2 - 3  r\_lat, 
    z = z\_0\)

(for \(TraditionalForm
      \`\(-3\) r\_lat ≤ x ≤ 3  r\_lat\)) also joins the points A and B. Generate a plot of this new path, superimposed upon the lattice. Calculate the work done along the path. How does its value compare with the work done along the straight line path?

(ii) Devise at least one more path of your own that joins A and B. Display it, and calculate the work done in moving a chloride ion along it. What is your tentative conclusion from your findings? Can you explain what's going on?