This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Numerics Volume, Chapter 1

In its simplest form, the Cauchy–Born hypothesis states that the positions of the lattice points (atoms) of a linearly deformed elastic body is given by the linear interpolations corresponding to the bodies faces [✶], [✶], [✶], [✶], [✶]. For the simplest example of a rectangular piece of a square lattice, carry out numerical tests of the Cauchy–Born hypothesis. To model the 2D elastic body, assume linear springs between all nearest and second nearest neighbors.

shear[{x_, y_}, {{αx_, βx_}, {αy_, βy_}}] := {αx x + βx y, αy x + βy y}

Given a discretized elastic 2D body in the form of a matrix (the elements being the positions of the atoms), the function energy calculates the potential energy of the underlying spring system with lattice constants ax, ay and spring constants k1, k2. We assume that the potential energy of each spring is of the form k/2 ((Δl)^2) where Δl is the length change of the spring from its equilibrium length.

(* potential energy of a lattice with spring constants 𝓀1, 𝓀2
   and lattice constant 𝒶x_, 𝒶y *)
energy[ℒ_, {ax_, ay_}, {k1_, k2_}] :=
With[{nx = Length[ℒ[[1]]], ny = Length[ℒ]}, 
 (* horizontal and vertical springs *)
 Sum[k1/2 (𝒹[ℯ[ℒ, {ix, jy}] - ℯ[ℒ, {ix + 1, jy}]] - ax)^2, 
     {jy, 0, ny - 1}, {ix, 0, nx - 2}] +
 Sum[k1/2 (𝒹[ℯ[ℒ, {ix, jy}] - ℯ[ℒ, {ix, jy + 1}]] - ay)^2, 
     {ix, 0, nx - 1}, {jy, 0, ny - 2}] +
(* diagonal springs connecting second nearest neighbors *)
 Sum[k2/2 (𝒹[ℯ[ℒ, {ix, jy}] - ℯ[ℒ, {ix + 1, jy + 1}]] - 𝒹[{ax, ay}])^2, 
     {ix, 0, nx - 2}, {jy, 0, ny - 2}] +
 Sum[k2/2 (𝒹[ℯ[ℒ, {ix, jy}] - ℯ[ℒ, {ix - 1, jy + 1}]] - 𝒹[{ax, ay}])^2, 
     {ix, 1, nx - 1}, {jy, 0, ny - 2}]];

(* some supporting functions *) 

(* distance between two lattice points *)
𝒹[{x_, y_}] := Sqrt[x^2 + y^2];

(* extract element from the lattice *)
ℯ[ℒ_, {ix_, jy_}] := ℒ[[jy + 1, ix + 1]];

(* visualize a spring lattice *) 
latticeLines[points_] := {Line /@ points, Line /@ Transpose[points],
                          Map[Point, points, {2}]}

To test the Cauchy–Born hypothesis, we will do the following: For a given deformed spring lattice, we calculate the positions of the nonboundary atoms by linear interpolation. Then, we disturb these positions slightly and use these perturbed positions as the starting points for an energy minimization carried out by FindRoot. Because the total potential energy is just the sum of the potential energy of all springs, the potential energy is a sum of squares and we will use the especially-suited-to-this-case "LevenbergMarquardt" method option setting. The function relaxedStretchedLatticeGraphics carries out these steps. It returns a list containing the potential energy of the linearly interpolated atom positions, the potential energy found by explicit minimization, and a graphics showing the two lattice configurations. (Because we will reuse some of the variables used inside relaxedStretched2DLatticeGraphics later, we do not declare all of the variables as local to this function).

relaxedStretched2DLatticeGraphics[
 {{αx_, βx_}, {αy_, βy_}}, {{nx_, ny_}, {ax_, ay_}, {k1_, k2_}}] := 
Module[{pRule = P[xy_] :> {x[xy], y[xy]}, 
        startLattice, distortedLattice, selectVertices,
        innerVerticesRules, fm, relaxedLattice},
(* general square lattice *)
topologicalLattice = Table[P[{x, y}], {y, 0, ny - 1}, {x, 0, nx - 1}];
(* square lattice with concrete lattice constants *)
startLattice = Map[#[[1]] {ax, ay}&,  topologicalLattice, {2}];
(* stretched square lattice *)
distortedLattice = Map[shear[#, {{αx, βx}, {αy, βy}}]&, 
                       startLattice, {2}];
(* select outer and inner vertices *)
selectVertices[cdc_] := Rule @@@ cdc[Flatten[
      MapThread[List, {topologicalLattice, distortedLattice}, 2], 1],
         {P[{0 | nx - 1, _}], _} | {P[{_, 0 | ny - 1}], _}];
(* rules for fixed boundary vertices *)
boundaryVerticesRules = selectVertices[Cases];
(* rules for movable inner vertices *)
innerVerticesRules = selectVertices[DeleteCases];        
(* potential energy *)
ℰ = energy[topologicalLattice /. boundaryVerticesRules /. pRule,
           {ax, ay}, {k1, k2}];
(* initial Cauchy-Born positions *)
inits = Flatten[Transpose /@ Flatten[Take[#, {2, -2}]& /@
 Take[MapThread[List, {topologicalLattice, distortedLattice}, 2] /. 
                                           pRule, {2, -2}], 1], 1];
(* minimize potential energy *)
fm = FindMinimum[Evaluate[ℰ], Evaluate[Sequence @@ 
          ({#1, (1 + Random[Real, 10^-0 {-1, 1}]) #2}& @@@ inits)],
          Method -> "LevenbergMarquardt", MaxIterations -> 1000];
(* the relaxed lattice *)
relaxedLattice = topologicalLattice /. pRule /. fm[[2]] /.
                 (boundaryVerticesRules /. pRule);   
(* return Cauchy–Born energy, minimized energy and graphic *)          
{{ℰ /. (Rule @@@ inits), fm[[1]]},
Graphics[(* Cauchy–Born lattice in gray *)
   {{PointSize[0.022], Thickness[0.015], GrayLevel[0.8],
     latticeLines[distortedLattice]},
         (* energy-minimized lattice in blue *)
    {PointSize[0.001], Thickness[0.002], RGBColor[0, 0, 1],
     latticeLines[relaxedLattice]},
         (* lattice point distortions in red *)
    {Thickness[0.002], Hue[0], MapThread[Line[{#1, #2}]&, 
     {distortedLattice, relaxedLattice}, 2]}}, 
    AspectRatio -> Automatic, Frame -> True]}]

Using the result returned from relaxedStretched2DLatticeGraphics, the function showResult prints the two energy values and displays a graphics showing the two lattices.

showResult[args__] := 
(Print["Cauchy–Born: " <> ToString[InputForm[#1[[1]]]] <> "\n" <>
       "FindMinimum: ", ToString[InputForm[#1[[2]]]]];    
   Show[#2])& @@ relaxedStretched2DLatticeGraphics[args]

Here is a first example. We use a 6×6 lattice. The deformation is relatively small. The linearly interpolated lattice (gray lines) and the lattice found by explicit energy minimization (in blue) are identical.

showResult[{{1.1, 0.1}, {-0.1, 0.9}}, 
           {{6, 6}, {1, 1}, {1, 1}}]

Here is a larger deformation. The two lattices again agree.

showResult[{{1.6, 0.3}, {0.5, 1.5}}, {{6, 6}, {1.5, 1}, {2.5, 2}}]

But it is not difficult to find a counter example to the Cauchy–Born hypothesis. In the next example, we use again a 6×6 lattice. This time the linearly interpolated lattice (gray lines) has a considerably higher potential energy. The red lines indicate the deviations of the energy-minimized lattices positions from the linearly interpolated ones.

SeedRandom[3];
showResult[{{0.62, 0.44}, {-0.34, 1.06}}, 
           {{6, 6}, {1.98, 0.56}, {2.24, 2.34}}]

Using NMinimize we can further minimize the energy of the lattice. (FindMinimum found one local potential energy minimum, lower minima might exist elsewhere.)

NMinimize[Evaluate[ℰ], First /@ inits][[1]]

We can watch the approach to the minimum by adding damping to the spring system and letting it relax into its equilibrium state. The function solveSpringSystemODEs creates the corresponding coupled differential equations for the lattice point coordinates and solves them.

solveSpringSystemODEs[T_, γ_, δ_] :=
Module[{vars, indices, odes}, 
 (* ODE variables *)
 vars = Cases[ℰ, _x | _y, Infinity] // Union;
 indices = Union[First /@ vars];
 (* equations of motion; with damping terms *)
 odes = (#''[t] == - γ #'[t] + (-D[ℰ, #] /. 
                         xy:(_x | _y) :> xy[t]))& /@ vars;
 NDSolve[Join[odes, (* randomized initial conditions *)
              (#1[0] == (1 + δ Random[Real, {-1, 1}]) #2)& @@@ inits, 
              (#'[0] == 0)& /@ vars], vars, {t, 0, T}]];

(* the time\[Hyphen]dependent relaxing lattices *)
relaxingLattice[τ_, ndsol_] :=  topologicalLattice /. 
   boundaryVerticesRules /. P[xy_] :> {x[xy], y[xy]} /.
   xy:(_x | _y) :> xy[τ] /. ndsol[[1]]

For four different random initial conditions, we let the last spring system now relax into the equilibrium shape. Starting from different initial conditions, we see different final configurations with different final energies. The second graphic shows the time-dependence of the potential energy.

Module[{T = 20, frames = 4, o = 4, γ = 1, δ = 0.2, ndsol,
        ax = 1.98, ay = 0.56, k1 = 2.24, k2 = 2.34,
        ℰCB = 27.485, minℰ = 18.386},
  Show[GraphicsArray[#]]& /@ Transpose[Table[
  (* relax four times *)
  Table[ndsol[j] = solveSpringSystemODEs[T, γ, δ];
       (* snapshots of the relaxation process *)
        Graphics[latticeLines[relaxingLattice[τ, ndsol[j]]],
                 Frame -> False], {τ, 0, T, T/frames}], {j, o}]];
  (* time\[Hyphen]dependent potential energy *)
  Show[Graphics[
  {{Thickness[0.006], GrayLevel[0.5], Line[{{0, minℰ}, {T, minℰ}}]},
   {Thickness[0.006], GrayLevel[0.0], Line[{{0, ℰCB}, {T, ℰCB}}]},
   Table[{(* color each run differently *) Hue[0.78 (j - 1)/(o - 1)], 
         Line[Table[{τ, energy[relaxingLattice[τ, ndsol[j]],  
                               {ax, ay}, {k1, k2}]},
                    {τ, 0, T, T/25}]]}, {j, o}]},
                    Frame -> True,
        PlotRange -> {1/2 minℰ, 2 minℰ}]]]

We end with some graphics showing violations of the Cauchy–Born hypothesis for some larger lattices. (For a detailed analysis of the local shape of the deformations, see [✶].)

Show[GraphicsArray[
Block[{$DisplayFunction = Identity}, showResult @@@
 {{{{0.63383, -0.10383}, {-0.422062, 1.36884}},
   {{16, 16}, {0.85177, 0.71945}, {1.24497, 1.48163}}},
  {{{1.11238, -0.25194}, {-0.023546, 0.519724}},
   {{16, 16}, {0.84189, 1.81386}, {0.43108, 0.99715}}},
  {{{0.89331, -0.03642}, { 0.0533719, 0.751591}},
   {{16, 16}, {0.69912, 0.68504}, {2.03383, 2.53343}}}}]]]