Energetics (Maths Lab 5)

Preliminary Code (select the cell bracket 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]);

Task Code

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

zforce = vectorforce[[3]]

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

Integrate[zforce, z]

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

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

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

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

forces = Map[forceFrom, lattice[2]]

vectorforce = Apply[Plus,  forces]

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

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

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

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

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

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

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

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

vectorforce = Apply[Plus,  forces];

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

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

Show[{latticeFigure, pathFigure}]

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

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