3.029 Spring 2022
Lecture 20 - 04/13/2022

Statistical Mechanics II

Ensembles Reminder

Last lecture, we introduced the concept of statistical ensembles, and how we can use them to extract meaningful thermodynamic quantities

We outlined the following general steps

Decide the relevant thermodynamic potential

This is often set by “experimental” conditions.
“Experimental” here might refer to conditions you can control in your numerical simulation (e.g. the number and total energy of the particles)

The natural variables of your thermodynamic potential define your statistical ensemble

Entropy  			->	NVE ensemble, aka Micro-canonical ensemble

Helmholtz free energy 	->	NVT ensemble, aka Canonical ensemble

Grand potential 		->	μVT ensemble, aka Grand-canonical ensemble

Compute the partition function for the statistical ensemble

NVE ensemble		->	“Degeneracy”,

NVE ensemble		->	“Partition Function”,

μVT ensemble		->	“Grand Partition Function”,

Calculate the relevant thermodynamic potential using the partition function

“Degeneracy”			->	Entropy,

“Partition Function”		->	Helmholtz free energy,

“Grand Partition Function”	->	Grand potential,

Use the differential form of the thermodynamic potential to obtain all other thermodynamic quantities

Entropy,			 		->	Inverse temperature, internal energy, etc..

Helmholtz free energy,				->	Entropy, internal energy, etc..

Grand potential,				->	Chemical potential, internal energy, etc..

Lennard Jones Dimer

In the staircase example we saw last time, the different microstates all had discrete energy levels

This is actually rather reminiscent of quantum-mechanical systems, which you will be study next semester

Let’s investigate a system with continuous energy levels, using our favourite interatomic potential

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

Plot[{lennardJonesPotentialNonDimensional[ρ],lennardJonesForce[ρ]},{ρ,3/4,3/2},PlotStyle->{Directive[Blue,Thick],Directive[Red,Thick]},Frame->True,FrameStyle->Directive[Black,Thick],FrameLabel->{"Separation distance, ρ","Energy, (Subscript[U, ljp])(ρ)"},BaseStyle->16,ImageSize->375]

We want to analyze the possible configurations a Lennard-Jones dimer can take in a Microcanonical ensemble

I.e. we need to evolve the system at constant energy

To do this, we will use the symplectic integration scheme we introduced in L13

parametricNDSolveValue=ParametricNDSolveValue[{
p'[t]==lennardJonesForce[q[t]], 
q'[t] == p[t],
 p[0] == 0,
q[0]==q0,
count[0]==0,
WhenEvent[p[t]==0,{count[t]->count[t]+1,If[count[t]>= 2,"StopIntegration"]}]
},{q,p},{t,0,100},q0,DiscreteVariables->{count},Method->{"SymplecticPartitionedRungeKutta","DifferenceOrder"->4,"PositionVariables"->{q[t]}},
StartingStepSize->1/1000,
AccuracyGoal->24,
PrecisionGoal->24,
WorkingPrecision->24]

Notes:

We used an auxiliary “count” variable to count the zero-momentum crossings and stop the integration scheme

We used ParametricNDSolveValue to vary the initial position of the particle

We can know plot the dimer’s position, q[t] vs its momentum, p[t]

I.e. the phase - space

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

Clear[phaseVariables]
phaseVariables[initialPosition_]:=phaseVariables[initialPosition]=Block[{psol,qsol,time},
{qsol,psol}=parametricNDSolveValue[initialPosition];
time=InterpolatingFunctionDomain[psol][[1,-1]];
{{psol,qsol},time}
]

If we only start at initial displacements (Subscript[ρ, 0])>(2^(-1/6)), we’re guaranteed to have a “bound state” which will oscillate around the minimum

Plot[lennardJonesPotentialNonDimensional[ρ],{ρ,3/4,3/2},PlotStyle->Directive[Blue,Thick],Frame->True,FrameStyle->Directive[Black,Thick],GridLines->{{(2^(-1/6))},None},FrameLabel->{"Separation distance, ρ","Energy, (Subscript[U, ljp])(ρ)"},BaseStyle->16,ImageSize->375]

Indeed, if we plot the momentum vs the position for a range of initial displacements

we obtain closed orbits

Length[range=Select[Subdivide[0.89175,2.125,102],lennardJonesPotentialNonDimensional[#]>-0.99&]]

plots =With[{solutions = phaseVariables/@range},
Show[
ParametricPlot[Through[First[#][t]],{t,0,Last[#]},PlotStyle->RandomColor[]]&/@solutions
,Frame->True,FrameStyle->Directive[Black,Thick],AspectRatio->Automatic,FrameLabel->{"Momentum, p[t]","Position, q[t]"},LabelStyle->Directive[14,Black],ImageSize->500,PlotRange->All]
]

Let’s ensure our simulations were indeed conserving energy

Seems like our percent error is at-most (10^(-6)) -- often much smaller

Manipulate[
With[{sol=phaseVariables[q0],ϵ0=lennardJonesPotentialNonDimensional[q0]},
Plot[((ϵ0-(lennardJonesPotentialNonDimensional[#1[[2]][t]]+((#1[[1]][t])^2)/2))/(ϵ0))100,{t,0,#2},GridLinesStyle->Directive[Gray,Dashed],PlotStyle->Red,Frame->True,PlotRange->All]&@@sol
],{{q0,1.4},(2^(-(1/6)))+0.05,1.5},Paneled->False]

Since we’re indeed evolving at constant-energy, each closed loop above is a microcanonical ensemble

I.e. each point on a closed loop, is a particular (continuous) microstate of the microcanonical ensemble with that particular energy

In L19, we mentioned that the degeneracy of the microcanonical ensemble can also be given by the “surface area” of the constant-energy surface

Since our phase-space is only 2D in this case, the degeneracy is given by the arclength of our trajectories

?ArcLength

ArcLength[{p[t],q[t]},{t,0,t0}]

Note however we need to rescale this such that the zero arc-length point, i.e. at ρ=1, has the ground state degeneracy of 1

As such, we’ll add 1.0 to our arc-length values

Clear[energyAndDegeneracy]
energyAndDegeneracy[initialPosition_]:=energyAndDegeneracy[initialPosition]=With[{sol=phaseVariables[initialPosition]},
{lennardJonesPotentialNonDimensional[initialPosition],1.0+Quiet@ArcLength[Through[#1[t]],{t,0,#2}]&@@sol}]

{minEnergy,maxEnergy}=MinMax[lennardJonesPotentialNonDimensional[range]]
lennardJonesDegeneracy=Interpolation[energyAndDegeneracy/@range]

Let’s plot it!

Plot[lennardJonesDegeneracy[ϵ],{ϵ,minEnergy,maxEnergy},Frame->True,LabelStyle->Directive[Black,Thick,16],ImageSize->500,PlotStyle->Directive[Red,Thick],FrameLabel->{"Internal Energy, U","Degeneracy, Ω(U,V,N)"}]

We can know use Boltzmann’s formula to get the entropy

lennardJonesEntropy[ϵ_]=Log[lennardJonesDegeneracy[ϵ]]

Quiet[Plot[lennardJonesEntropy[ϵ],{ϵ,minEnergy,maxEnergy},Frame->True,LabelStyle->Directive[Black,Thick,16],ImageSize->500,PlotStyle->Directive[Red,Thick],FrameLabel->{"Internal Energy, U","Entropy, S(U,V,N)"}]]

Let’s plot internal energy as a function of temperature expression parametrically

lennardJonesTemperature[ϵ_]=1/lennardJonesEntropy'[ϵ];

ParametricPlot[{lennardJonesTemperature[ϵ],ϵ},{ϵ,minEnergy,maxEnergy},Frame->True,LabelStyle->Directive[Black,Thick,16],ImageSize->500,PlotStyle->Directive[Blue,Thick],FrameLabel->{"Normalized Temperature, τ","Internal Energy, U"}]

We see that our internal energy function is multi-valued (and thus non-invertible)

This is due to the change of convexity in our entropy function

Quiet[Plot[lennardJonesEntropy''[ϵ],{ϵ,minEnergy,maxEnergy},Frame->True,LabelStyle->Directive[Black,Thick,16],ImageSize->500,PlotStyle->Directive[Red,Thick],FrameLabel->{"Internal Energy, U","Entropy Second Derivative"}]]

cutoffEnergy=ϵ/.FindRoot[lennardJonesEntropy''[ϵ]==0,{ϵ,-0.2}]

Will investigate the lower-energy part using a canonical-ensemble approach instead

Remember, the ensembles should all be equivalent!

lennardJonesPartitionFunction[τ_]=With[{energiesAndDegeneracies=Select[energyAndDegeneracy/@range,First[#]<cutoffEnergy&]},Total[#2Exp[-#1/τ]&@@@energiesAndDegeneracies]];
lennardJonesFreeEnergy[τ_]=-τ Log[lennardJonesPartitionFunction[τ]];
lennardJonesEntropyCanonical[τ_]=-lennardJonesFreeEnergy'[τ];
lennardJonesInternalEnergy[τ_]=lennardJonesFreeEnergy[τ]+τ lennardJonesEntropyCanonical[τ];
lennardJonesHeatCapacity[τ_]=lennardJonesInternalEnergy'[τ];

{minTemperature,maxTemperature}=lennardJonesTemperature/@{minEnergy,cutoffEnergy};
Quiet@Plot[lennardJonesHeatCapacity[τ],{τ,minTemperature,maxTemperature},PlotStyle->Directive[Blue,Thick],Frame->True,LabelStyle->Directive[Black,Thick,16],ImageSize->500,FrameLabel->{"Normalized Temperature","Heat Capacity"}]

We recover the familiar finite-size effect at high temperatures

But our numerical investigation is hampered at low temperatures

It’s fair to speculate this would approach 0 as temperature -> 0

Frenkel Defects

Consider a perfect square lattice (ground state)

groundState[numberOfAtoms_]:=Graphics[{Disk[#,1/(2numberOfAtoms)]&/@Tuples[Subdivide[-1,1,numberOfAtoms-1],2]},PlotRange->((numberOfAtoms+1)/(numberOfAtoms-1)){{-1,1},{-1,1}},ImageSize->100]

groundState[10]

Create Frenkel defects (vacancy - interstitial pair)

The ground state has energy 0, while each Frenkel defect ‘costs’ an energy of ϵ

frenkelDefect[numberOfAtoms_,numberOfDefects_:1]:=With[{gs=Tuples[Subdivide[-1,1,numberOfAtoms-1],2]},
Graphics[{
Disk[#,1/(2numberOfAtoms)]&/@Union[RandomSample[gs,Length[gs]-numberOfDefects]],
Red,Disk[#+{1,1}/(numberOfAtoms-1),1/(2numberOfAtoms)]&/@Union[RandomSample[gs,numberOfDefects]]},
PlotRange->((numberOfAtoms+1)/(numberOfAtoms-1)){{-1,1},{-1,1}},ImageSize->100]]

Multicolumn[{groundState[10],frenkelDefect[10,10]},2]

Since we can choose the vacancy in one of N positions, and similarly the interstitial in N, the # of different configurations is given by:

Ω(N)=(N^2)

binomial[n_,k_]:=ExpandAll[((n!)/((n-k)!k!))]

((binomial[9,1])^2)
Table[frenkelDefect[3],5000]//Union//Length

Multicolumn[Table[frenkelDefect[3],5000]//Union,9,Frame->All,Appearance->"Horizontal"]

If we allow for m Frenkel defects in a lattice with N atoms, then this would be given by:

Ω(N,m)=(((({{N}, {m}})))^2)

((binomial[9,2])^2)
Table[frenkelDefect[3,2],10000]//Union//Length

FunctionExpand[((binomial[n,m])^2)]

The canonical partition function for a system with m Frenkel defects is given by

partitionFunction["frenkel"][n_,m_,ϵ_][τ_]=Sum[(((Gamma[1+n])^2)/(((Gamma[1+M])^2) ((Gamma[1-M+n])^2)))Exp[-M ϵ/τ],{M,0,m}]

Where we’ve used the Gamma function to expand the degeneracy of multiple Frenkel defects.

We then proceed as usual to recover the free energy, entropy, and internal energy

helmholtzFreeEnergy ["frenkel"][n_,m_,ϵ_][τ_]= -τ Log[partitionFunction["frenkel"][n,m,ϵ][τ]];

entropy["frenkel"][n_,m_,ϵ_][τ_]=-helmholtzFreeEnergy ["frenkel"][n,m,ϵ]'[τ];

internalEnergy["frenkel"][n_,m_,ϵ_][τ_]=helmholtzFreeEnergy ["frenkel"][n,m,ϵ][τ]+ τ entropy["frenkel"][n,m,ϵ][τ];

heatCapacity["frenkel"][n_,m_,ϵ_][τ_]=internalEnergy["frenkel"][n,m,ϵ]'[τ];

Plotting the heat capacity we get the familiar “finite-system” form

Quiet@Plot[heatCapacity["frenkel"][10,1,1][τ],{τ,0,1},Frame->True,FrameStyle->Directive[Black,Thick],PlotRange->All,PlotStyle->Directive[Red,Thick],BaseStyle->18,ImageSize->500,FrameLabel->{"Normalized Temperature","Heat Capacity"}]

Allowing more Frenkel defects broadens the heat capacity, up to a point

Quiet@Plot[heatCapacity["frenkel"][10,3,1][τ],{τ,0,1},Frame->True,FrameStyle->Directive[Black,Thick],PlotRange->All,PlotStyle->Directive[Blue,Thick],BaseStyle->18,ImageSize->500,FrameLabel->{"Normalized Temperature","Heat Capacity"}]

Quiet@Plot[heatCapacity["frenkel"][10,5,1][τ],{τ,0,1},Frame->True,FrameStyle->Directive[Black,Thick],PlotRange->All,PlotStyle->Directive[RGBColor[0,1/3,0],Thick],BaseStyle->18,ImageSize->500,FrameLabel->{"Normalized Temperature","Heat Capacity"}]

Quiet@Plot[heatCapacity["frenkel"][10,10,1][τ],{τ,0,1},Frame->True,FrameStyle->Directive[Black,Thick],PlotRange->All,PlotStyle->Directive[Black,Thick],BaseStyle->18,ImageSize->500,FrameLabel->{"Normalized Temperature","Heat Capacity"}]