3.029 Spring 2022
Lecture 17 - 04/04/2022

Surface Tension & Energy

This week, we will develop the thermodynamics of interfaces in crystalline materials

In particular, we’ll investigate the surface tension terms arising from finite size particles using the Lennard-Jones potential

Next lecture, we’ll investigate the implications of this orientational dependence on the equilibrium shapes of crystals using the Wulff construction

Surface Tension in Lennard-Jones Particles

Finite-Size Particles

Recall the normalized Lennard-Jones potential for a pair of atoms a distance ρ=(r/(Subscript[r, min])) apart is given by:

lennardJonesPotentialNonDimensional[ρ_]=(1/(ρ^12))-(2/(ρ^6));

The sum of pairwise energies each atom in a particle experiences is similarly given by:

Note, we are dividing by two to avoid double-counting, and the 1/2 term is to account for the additional self-energy terms we introduced by setting the diagonal to ρ=1

lennardJonesEnergies[pts_?MatrixQ]:= Block[
{
distanceMatSquared =DistanceMatrix[pts, DistanceFunction->SquaredEuclideanDistance] + IdentityMatrix[Length[pts]],
distanceMatSixth
},
distanceMatSixth= (distanceMatSquared^3);
(1/2)+(1/2)Total[(1/(distanceMatSixth^2)) - (2/distanceMatSixth)]
]

As we saw in L13, non-hexagonal lattices undergo a lattice instability when interacting with pairwise central potentials in 2D

As such, we’ll focus on hexagonal particles

We first make a list of hexagonal lattice sites from which we can “punch-out” round particles of a specified radius

We need to ensure it’s large enough for the largest particle we’d like to punch out, say 50 normalized units

hexagonalLatticeVectors = {{0,1},{-((Sqrt[3])/2),(1/2)}};
hexagonalLatticePoints=With[{range=Ceiling[50/((Sqrt[3])/2)]},N[Tuples[Range[-range,range],2].hexagonalLatticeVectors]];
Graphics[{Point[hexagonalLatticePoints],White,Thick,Circle[{0,0},#]&/@Range[5,50,5]},PlotRange->50{{-1,1},{-1,1}},PlotRangePadding->1]

And write a simple “punch” function

hexagonalParticle[radius_,latticePoints_:hexagonalLatticePoints]:=hexagonalParticle[radius,latticePoints]= Select[Norm[#]<radius&][latticePoints]

To see the effect of surface tension, let’s color each particle by its energy

energyColor[{eMin_,eMax_},power_:0.25,cf_:ColorData["DarkRainbow"]][energy_]:= cf[
(((Rescale[Clip[energy,{eMin,eMax}],{eMin,eMax},{0,1}]))^power)]

particleGraphic[{eMin_,eMax_},power_:0.25,cf_:ColorData["DarkRainbow"]][positions_]:= MapThread[{energyColor[{eMin,eMax},power,cf][#1],Disk[#2,1/2]}&,{lennardJonesEnergies[positions],positions}]

{ljEnergyMinimum,ljEnergyMaximum}=MinMax[lennardJonesEnergies[hexagonalParticle[50]]];
Multicolumn[Graphics@*particleGraphic[{ljEnergyMinimum,ljEnergyMaximum},0.25]@*hexagonalParticle/@{5,10,15,20},4]

Notice we used a sub-linear power law in our color function (n=0.25)

This is to give more contrast to the lower parts of the energy distribution and increase the “diffuse” width of the surface

Plot[{x,(x^0.25),(x^4)},{x,0,1},Frame->True]

As we’ve seen before, atoms on the surface have a higher average energy

Excess Average Surface Energy

We’ll use these particles to estimate the excess energy due to the presence of a surface

The idea is to try and make increasingly larger and larger particles and compute the average energy per particle

As the particle becomes larger, the fraction of particles on the surface becomes smaller, and thus contributing less to the average energy

particles = <|# -> hexagonalParticle[#]&/@Range[5,50] |>;

aggregatedLennardJonesEnergies[positions_]:=aggregatedLennardJonesEnergies[positions]=With[{es=lennardJonesEnergies[positions]},<|"Energies"->es,"Total Energy"->Total[es],"Average Energy"->Mean[es]|>]

The magnitude of the total energy naturally increases with particle radius

Remember, the potential is normalized to have a value of -1 at a normalized distance of ρ=1

With[{totalEnergies=(aggregatedLennardJonesEnergies/@particles)[[All,"Total Energy"]]},
ListPlot[totalEnergies,Frame->True,FrameLabel->{"Particle Radius","Total Energy"},LabelStyle->Directive[Black,16],ImageSize->500,PlotStyle->Directive[PointSize[0.025],Opacity[0.5],Blue]]
]

Without a surface, we would expect the total energy to decrease linearly with the number of particles

Let’s plot this as an average energy per particle to see this better

With[{averageEnergies=(aggregatedLennardJonesEnergies/@particles)[[All,"Average Energy"]]},
ListPlot[averageEnergies,Frame->True,FrameLabel->{"Particle Radius","Average Energy"},LabelStyle->Directive[Black,16],ImageSize->500,PlotStyle->Directive[PointSize[0.025],Opacity[0.5],Red]]
]

We’ll use this data to fit a power-law model, to extrapolate the average energy as the particle size goes to infinity

averageEnergyModel=constantValue + sizeDependentValue /(radius^power);
averageEnergyModelFit=With[{averageEnergies=(aggregatedLennardJonesEnergies/@particles)[[All,"Average Energy"]]},
FindFit[KeyValueMap[List,averageEnergies],averageEnergyModel,{constantValue,sizeDependentValue,power},radius]
]

With[{averageEnergies=(aggregatedLennardJonesEnergies/@particles)[[All,"Average Energy"]],min=constantValue/.averageEnergyModelFit},
Show[Plot[averageEnergyModel/.averageEnergyModelFit,{radius,5,150},Frame->True,FrameLabel->{"Particle Radius","Average Energy"},LabelStyle->Directive[Black,16],PlotRangePadding->{Automatic,Scaled[.05]},PlotRange->{min,Automatic},ImageSize->500,PlotStyle->Blue,GridLinesStyle->Directive[Gray, Dashed],GridLines->{None,{min}}],ListPlot[averageEnergies,PlotStyle->Directive[PointSize[0.025],Opacity[0.5],Red]]
]]

The excess energy per atom goes like

Excess Total Surface Energy

This is related to the Gibbs-Thomson effect

And we can relate this to the curvature of our particle by looking at the rate of surface energy increase with respect to volume

Note: since we’re in 2D, “volume” here refers to the particle’s area, and “area” refers to the circumference

With[{
volume=π (r^2),
area = 2π r
},
((D[area ,r])/(D[volume,r]))
]

In-fact, let’s use this form to provide another estimate for the surface tension based on the total energy

totalEnergyModel=bulkEnergy (Pi (radius^2)) + γ (2 Pi radius);
totalEnergyModelFit=With[{totalEnergies=(aggregatedLennardJonesEnergies/@particles)[[All,"Total Energy"]]},
FindFit[KeyValueMap[List,totalEnergies],totalEnergyModel,{bulkEnergy,γ},radius]
]

With[{totalEnergies=(aggregatedLennardJonesEnergies/@particles)[[All,"Total Energy"]]},
Show[Plot[totalEnergyModel/.totalEnergyModelFit,{radius,5,50},Frame->True,FrameLabel->{"Particle Radius","Total Energy"},LabelStyle->Directive[Black,16],ImageSize->500,PlotStyle->Red],ListPlot[totalEnergies,PlotStyle->Directive[PointSize[0.025],Opacity[0.5],Blue]]
]]

Looking at our two fits we notice a curious discrepancy in the constant (surface-independent) terms

totalEnergyModelFit

averageEnergyModelFit

((constantValue/.averageEnergyModelFit)/(bulkEnergy/.totalEnergyModelFit))

The discrepancy is resolved by noting that the bulk energy fit is per unit-area, while the average energy fit is per atom

Each atom in our hexagonal lattice takes up (Sqrt[3])/2 normalized units of area

We can see this directly from the Determinant of the hexagonal lattice vectors

Det[hexagonalLatticeVectors]//N

As well as more visually/numerically using a Voronoi Mesh

Block[{voronoi=VoronoiMesh[hexagonalParticle[6]],interior,pts=Select[hexagonalParticle[6],Norm[#]<5&]},
interior=MeshPrimitives[voronoi,{2,"Interior"}];
interior=Select[interior,Norm[RegionCentroid[#]]<5&];
Echo[Graphics[{
{EdgeForm[Black],FaceForm[None],interior},
{Red,PointSize[Large],Point[pts]}
}]];
Mean[Area/@interior]
]

Thermodynamics of Surfaces

Notice we used the symbol γ above for the surface-energy dependent term we multiplied by the surface’s “area”

This γ term is called the surface-energy and we quickly summarize its canonical thermodynamic treatment, following Gibbs

To do so, we need to introduce the concept of Gibbs dividing surface, which is most conveniently demonstrated by considering an interface between a bicrystal

Plot[LogisticSigmoid[x],{x,-25,25},Frame->True,Axes->False,GridLines->{{0},{0,1}},GridLinesStyle->Directive[Gray, Dashed],PlotStyle->{Red,Thick},ImageSize->500,LabelStyle->Directive[Black,16],FrameTicks->{{{{0,"Phase α"}},{{1,"Phase β"}}},{{{0,"Gibbs Dividing Surface"}},None}},FrameLabel->{None,"Extensive Quantity"}]

In general, interfaces between various thermodynamic phases are “diffuse”

i.e. there is a transition region separating the two phases (defined using their extensive quantities above)

To avoid the complications arising from this diffuseness, Gibbs defined a mathematical plane with area A, known as the dividing surface

This allowed him to define interfacial extensive quantities by subtracting from the extensive quantity of the system containing the interface a hypothetical reference system consisting of two bulk phases, extending uniformly up to the dividing surface

Recall that the variation of the internal energy of bulk phase α with respect to its extensive variables is given by

and similarly for phase β:

For the system containing the interface we can similarly write

where we note that the interface has no volume, but instead introduces a surface energy term multiplied by the change in area dA

Since the two phases are at equilibrium, we note that , and accumulate the other extensive variables are

Finally, since all the differentials above are extensive, we can integrate the expression to give

Surface Reconstructions

Uniform Strain Relaxation

All of the calculations above assumed a particle lattice constant at ρ=1

This was the normalized equilibrium lattice constant for a pair of atoms

However, the long-range attractive force of the potential suggests that the particle will experience non-negligible second- and third-nearest neighbor interactions too

This will cause the lattice constant to be slightly different

As a first pass, we can try and find this new lattice constant for uniform strain

With[{positions = hexagonalParticle[10]},
Plot[
aggregatedLennardJonesEnergies[newLatticeConstant positions]["Total Energy"],
{newLatticeConstant,0.9,1.1},
Frame->True,FrameLabel->{"Lattice Constant","Total Energy"},ImageSize->500,LabelStyle->Directive[Black,16],PlotStyle->Directive[Thick,Blue]
]
]

We can see visually that the minimum is slightly less than 1

Let’s investigate the total energy for small particles symbolically

Graphics[particleGraphic[{ljEnergyMinimum,ljEnergyMaximum}][hexagonalParticle[2]]]

With[{positions=hexagonalParticle[2]},
Apart[FullSimplify[
aggregatedLennardJonesEnergies[newLatticeConstant positions]["Total Energy"],newLatticeConstant>0]]
]

Graphics[particleGraphic[{ljEnergyMinimum,ljEnergyMaximum}][hexagonalParticle[2.5]]]

With[{positions=hexagonalParticle[2.5]},
Apart[FullSimplify[
aggregatedLennardJonesEnergies[newLatticeConstant positions]["Total Energy"],newLatticeConstant>0]]
]

Graphics[particleGraphic[{ljEnergyMinimum,ljEnergyMaximum}][hexagonalParticle[3]]]

With[{positions=hexagonalParticle[3]},
Apart[FullSimplify[
aggregatedLennardJonesEnergies[newLatticeConstant positions]["Total Energy"],newLatticeConstant>0]]
]

This suggests we can minimize the total energy symbolically as:

symbolicTotalEnergy[newLatticeConstant_] = (twelfthPowerCoefficient/(newLatticeConstant^12))  -2(sixthPowerCoefficient/(newLatticeConstant^6));
First[SolveValues[symbolicTotalEnergy'[newLatticeConstant]==0,newLatticeConstant]]

We can extract the coefficients using a very similar function as our total energy function

minimizingLatticeConstant[pts_]:=
Block[{
distanceMatSquared =DistanceMatrix[pts, DistanceFunction->SquaredEuclideanDistance] + IdentityMatrix[Length[pts]],
twelfthPowerCoefficient,sixthPowerCoefficient
},
twelfthPowerCoefficient=  Total[(1/(distanceMatSquared^6)),2]-Length[pts];
sixthPowerCoefficient=Total[(1/(distanceMatSquared^3)),2]-Length[pts] ;
((( twelfthPowerCoefficient)/sixthPowerCoefficient))^(1/6)

]

latticeConstantData=minimizingLatticeConstant/@particles

ListPlot[latticeConstantData,Frame->True,FrameLabel->{"Particle Radius","Relaxed Lattice Constant"},LabelStyle->Directive[Black,16],ImageSize->500,PlotStyle->Directive[PointSize[0.025],Opacity[0.5],Red]]

The plot above motivates us to use a statistical approach to fit a functional form to our new lattice constant

newLatticeConstantModel=limitingLatticeConstant + powerCoefficient/radius^power;
newLatticeConstantModelFit=With[{newLatticeConstants=KeyValueMap[List,latticeConstantData]},
FindFit[newLatticeConstants,newLatticeConstantModel,{limitingLatticeConstant,powerCoefficient,power},radius]
]

With[{min=limitingLatticeConstant/.newLatticeConstantModelFit},
Show[Plot[newLatticeConstantModel/.newLatticeConstantModelFit,{radius,5,150},Frame->True,FrameLabel->{"Particle Radius","Relaxed Lattice Constant"},LabelStyle->Directive[Black,16],PlotRangePadding->{Automatic,Scaled[.05]},PlotRange->{min,Automatic},ImageSize->500,PlotStyle->Blue,GridLinesStyle->Directive[Gray, Dashed],GridLines->{None,{min}}],ListPlot[latticeConstantData,PlotStyle->Directive[PointSize[0.025],Opacity[0.5],Red]]
]]

The limiting constant appears to be ~1% smaller than unity

The differences are small, but we might as well define a rule to return this as a function of radius

equilibriumLatticeConstant[radius_]=newLatticeConstantModel/.newLatticeConstantModelFit

{ljRelaxedEnergyMinimum,ljRelaxedEnergyMaximum}=MinMax[lennardJonesEnergies[hexagonalParticle[50,hexagonalLatticePoints equilibriumLatticeConstant[50]]]];
Multicolumn[Graphics[particleGraphic[{ljRelaxedEnergyMinimum,ljRelaxedEnergyMaximum}][hexagonalParticle[#,hexagonalLatticePoints equilibriumLatticeConstant[#]]]]&/@{5,10,15,20},4]

Turns out, this has a fairly significant effect on the surface energy!

strainRelaxedParticles = <|# -> hexagonalParticle[#,equilibriumLatticeConstant[#]hexagonalLatticePoints]&/@Range[5,50] |>;

totalEnergyModelFitStrainRelaxed=With[{totalEnergies=(aggregatedLennardJonesEnergies/@strainRelaxedParticles)[[All,"Total Energy"]]},
FindFit[KeyValueMap[List,totalEnergies],totalEnergyModel,{bulkEnergy,γ},radius]
]

Full Relaxation

Since the uniform strain had such an effect on surface energy, it stands to reason whether our approximation for uniform strain is valid

In theory, there’s nothing stopping the particles in relaxing non-uniformly

and in-fact there's a very rich literature of so-called “surface-reconstruction” phenomena

One approach would be to minimize the total energy with respect to all the particles’ degrees of freedom

Let’s start small, by considering a particle with radius 3

Clear[totalLennardJonesEnergyVector]
totalLennardJonesEnergyVector[listOfPositions_?(VectorQ[#,NumericQ]&)]:=Total[lennardJonesEnergies[Partition[listOfPositions,2]]]

equilibriumPositions[3]=hexagonalParticle[3];
variables[3] = Flatten[Table[{Subscript[x,i],Subscript[y,i]},{i,Length[equilibriumPositions[3]]}]];
initialGuess[3] = Thread[{variables[3], Flatten[equilibriumPositions[3]]}];
AbsoluteTiming[{minval[3],minpos[3]}=FindMinimum[totalLennardJonesEnergyVector[variables[3]],initialGuess[3]]]

With[{us=equilibriumPositions[3]-Partition[variables[3]/.minpos[3],2]},
ListVectorPlot[Thread[{equilibriumPositions[3],us}],VectorMarkers->Placed["Arrow","Start"],VectorPoints->equilibriumPositions[3],Epilog->{Circle[#,1/8]&/@equilibriumPositions[3]},VectorScaling->Automatic,PlotRangePadding->Scaled[.1],Frame->False,ImageSize->500]
]

The result seems correct, but notice we get a numerical accuracy warning

We can trying helping Mathematica a bit, by providing the Gradient used in the Quasi-Newton algorithm symbolically

symbolicGradient[variables__]:=symbolicGradient[variables]=Module[{energy,listVars={variables}},
energy[variables]=Total[lennardJonesEnergies[Partition[listVars,2]]]/.((Abs[((Abs[a_])^2)+((Abs[b_])^2)])^(c_)):>((((a^2)+(b^2)))^c);
Grad[energy[variables],listVars]
]

AbsoluteTiming[{minval[3],minpos[3]}=FindMinimum[totalLennardJonesEnergyVector[variables[3]],initialGuess[3],Gradient->symbolicGradient@@variables[3]];]

Not only does this not complain anymore, it is also much faster

However, when the search space is even larger, it still fails

equilibriumPositions[4.5]=hexagonalParticle[4.5];
variables[4.5] = Flatten[Table[{Subscript[x,i],Subscript[y,i]},{i,Length[equilibriumPositions[4.5]]}]];
initialGuess[4.5] = Thread[{variables[4.5], Flatten[equilibriumPositions[4.5]]}];
AbsoluteTiming[{minval[4.5],minpos[4.5]}=FindMinimum[totalLennardJonesEnergyVector[variables[4.5]],initialGuess[4.5],Gradient->symbolicGradient@@variables[4.5]]];

With[{us=equilibriumPositions[4.5]-Partition[variables[4.5]/.minpos[4.5],2]},
ListVectorPlot[Thread[{equilibriumPositions[4.5],us}],VectorMarkers->Placed["Arrow","Start"],VectorPoints->equilibriumPositions[4.5],Epilog->{Circle[#,1/8]&/@equilibriumPositions[4.5]},VectorScaling->Automatic,PlotRangePadding->Scaled[.1],Frame->False,ImageSize->500]
]

equilibriumPositions[5]=hexagonalParticle[5];
variables[5] = Flatten[Table[{Subscript[x,i],Subscript[y,i]},{i,Length[equilibriumPositions[5]]}]];
initialGuess[5] = Thread[{variables[5], Flatten[equilibriumPositions[5]]}];
{minval[5],minpos[5]}=FindMinimum[totalLennardJonesEnergyVector[variables[5]],initialGuess[5],Gradient->symbolicGradient@@variables[5]];

With[{us=equilibriumPositions[5]-Partition[variables[5]/.minpos[5],2]},
ListVectorPlot[Thread[{equilibriumPositions[5],us}],VectorMarkers->Placed["Arrow","Start"],VectorPoints->equilibriumPositions[5],Epilog->{Circle[#,1/8]&/@equilibriumPositions[5]},VectorScaling->Automatic,PlotRangePadding->Scaled[.1],Frame->False,ImageSize->500]
]

We could use a Gradient-free method like Interior Point

However, this is much too slow to be practical

{minval[5],minpos[5]}=FindMinimum[totalLennardJonesEnergyVector[variables[5]],initialGuess[5],Method->"InteriorPoint",MaxIterations->100];

With[{us=equilibriumPositions[5]-Partition[variables[5]/.minpos[5],2]},
ListVectorPlot[Thread[{equilibriumPositions[5],us}],VectorMarkers->Placed["Arrow","Start"],VectorPoints->equilibriumPositions[5],Epilog->{Circle[#,1/8]&/@equilibriumPositions[5]},VectorScaling->Automatic,PlotRangePadding->Scaled[.1],Frame->False,ImageSize->500]
]

Shell Relaxation

Looking at the arrows above, it’s clear that the displacement is mostly on the outermost atoms

Perhaps we can combine the two methods

I.e. uniform strain in the core of the particle, with full relaxation on the outer shell

carveInnerShell[radiusOuter_,dr_:2]:=With[{
uniformStrained=hexagonalParticle[radiusOuter,equilibriumLatticeConstant[radiusOuter] hexagonalLatticePoints],
radiusInner=radiusOuter-dr},
Values[KeySort[GroupBy[uniformStrained,Norm[#]>radiusInner&]]]
]

Graphics@*particleGraphic[{ljRelaxedEnergyMinimum,ljRelaxedEnergyMaximum}]/@carveInnerShell[5]

semiSymbolicGradient[innerPositions_][variables__]:=semiSymbolicGradient[innerPositions][variables]=Module[{energy,listVars={variables}},
energy[variables]=Total[lennardJonesEnergies[Join[innerPositions,Partition[listVars,2]]]]/.((Abs[((Abs[a_])^2)+((Abs[b_])^2)])^(c_)):>((((a^2)+(b^2)))^c);
Grad[energy[variables],listVars]
]

Monitor[
Do[
With[{dr=2},
{innerPositions[radius],outerPositions[radius]}=carveInnerShell[radius,dr];
outerVariables[radius] = Flatten[Table[{Subscript[x,i],Subscript[y,i]},{i,Length[outerPositions[radius]]}]];
initialGuessOuter[radius] = Thread[{outerVariables[radius], Flatten[outerPositions[radius]]}];
{minvalOuter[radius],minposOuter[radius]}=FindMinimum[totalLennardJonesEnergyVector[Join[Flatten[innerPositions[radius]],outerVariables[radius]]],initialGuessOuter[radius],Gradient->semiSymbolicGradient[innerPositions[radius]]@@outerVariables[radius]];
minposShells[radius]=Join[innerPositions[radius],Partition[outerVariables[radius]/.minposOuter[radius],2]];
],{radius,2,9}],radius]

shellRelaxedParticles = <|# -> minposShells[#]&/@Range[2,9] |>;

Let’s visualize the shell-relaxed displacements

Multicolumn[With[{ps=Join[innerPositions[#],outerPositions[#]]},
ListVectorPlot[Thread[{ps,ps-minposShells[#]}],VectorMarkers->Placed["Arrow","Start"],VectorPoints->ps,Epilog->{Circle[#,1/8]&/@ps},VectorScaling->Automatic,PlotRangePadding->Scaled[.1],Frame->False,ImageSize->200]
]&/@Range[2,9],4,Appearance->"Horizontal"]

In-fact the displacements are even less than the 1% strain relaxation

MinMax[Norm/@(Join[innerPositions[9],outerPositions[9]]-shellRelaxedParticles[9])]

Multicolumn[Graphics@*particleGraphic[{ljRelaxedEnergyMinimum,ljRelaxedEnergyMaximum}]/@Values[shellRelaxedParticles],4,Appearance->"Horizontal"]

totalEnergyModelFitStrainRelaxed=With[{totalEnergies=(aggregatedLennardJonesEnergies/@shellRelaxedParticles)[[All,"Total Energy"]]},
FindFit[KeyValueMap[List,totalEnergies],totalEnergyModel,{bulkEnergy,γ},radius]
]

As expected, the surface reduces even further

Note however we didn’t use too large particles for this