Atomic Clusters

APAM 1601
Columbia University

Introduction

In this notebook, we continue to explore the use of random walks and the Rosenbluth-Teller-Metropolis algorithm to understand complicated physical systems in "thermal equilibrium".  For this analysis, we will explore the three-dimensional structure of clusters of "atoms" interacting through the Lennard-Jones potential. Sir John Edward Lennard-Jones (1894-1954) was a leading theoretical chemist of his time. His simple (yet highly nonlinear) model for the potential energy between two atoms or molecules has become a much studied system of molecular configuration. The potential energy between two atoms is a combination of an attractive van der Waals force (due to the electric polarization of atoms) and a repulsive term resulting from the overlap of electron orbits. Today's research trying to understand the configurations found in inorganic nanoscience and in protein folding continue to make these sorts of models an active area of research.

lj[r_] := 1/r^12 - 2 /r^6

lj2[r2_] := 1/r2^6 - 2 /r2^3 (* where r2 = r . r *)

Plot[lj[r],{r,0.001,2.},PlotRange->{-1,1},AxesLabel->{"r","ϕ(r)"},PlotLabel->"Lennard-Jones Potential"]

Atoms are attracted to each other, and a minimum energy is found when separated at r~ 1. When the atoms move closer together, a strong potential gradient separates them.

Monte Carlo Simulation Rules

The "rules" for our numerical simulation are similar to those used for many problems in statistical physics (like an ideal non-interacting atmosphere)–with one difference. For particles interacting by way of the Lennard-Jones potential, all particles are strongly interacting with each other. Many large problems in statistical physics assume that particles only interact with nearest neighbors. When a particle moves, we need to compute the change in system energy by summing over all interactions.

We will follow the following set of rules:
Step 1) Initialize the positions of the particles within a sphere of given radius..
Step 2) Select a random particle. Make a potential random step  in any direction. However, if the step causes the particle to escapee 
Step 3) Calculate the energy change to the entire system if the particle makes the step.
Step 4a) If the energy change is negative, accept the step.
Step 4b) If the energy change is positive, pick a random number between 0 and 1. If the random number is less than Exp[ - (energy change)/kT], then accept the step.
Step 4c) Otherwise, reject the step and leave the position (and energy) unchanged.
Step 5) Go to Step 2 and repeat (for a very long time....)

These rules comprise the "Metropolis" method for computational statistical physics. It is problably the most often used computational method in physics! See online biography for Nick Metropolis. The real insight behind the algorythm was developed by late Marshall Rosenbluth, the father of modern plasma physics.

The Potential Energy of a Lennard-Jones Cluster

The total potential energy for N particles will involve the sum between all pairs of particles:
	(Subscript[U, tot]) = (∑_(i=2)^N)(∑_(j=1)^(i-1))ϕ(|(Subscript[x, i]) - (Subscript[x, j])|)
where ϕ(r) is the Lennard-Jones potential. For N particles, this sum has N (N – 1)/2 pairings. (For this reason, we will not be able to simulate many particles!)

When only one particle moves, the change in potential energy, ΔU, requires only those pairings involving the particle that moved. The change in energy is
	ΔU = (∑_(i≠k))[ϕ((Subscript[x, i]) - ((Subscript[x, k]) + Δx)) - ϕ((Subscript[x, i]) - (Subscript[x, k]))]

If ΔU > 0, the  step is accepted only if a random number between 0 and 1 is less than Exp[-ΔU / T]. If ΔU < 0, the step is always accepted.

Initializing the Cluster

maxR = 4.0; (* all particles must remain within maxR *)

nWalkers=12;
walkers=Table[RandomReal[{-(maxR/2),(maxR/2)},3],{nWalkers}];

walkers // Dimensions

ListPointPlot3D[walkers,PlotStyle->{PointSize[0.02]},PlotRange->{{-maxR,maxR},{-maxR,maxR},{-maxR,maxR}}]

Evaluating the total potential

ϕlj[rvec_] := lj2[rvec.rvec] (* ϕlj[r] takes a vector as argument *)

Plot[ϕlj[{x,0,0.5}],{x,-2.,2.},PlotRange->{-1,1}]

ϕlj[walkers[[1]] - walkers[[3]]]

potE[w_List] := Sum[Sum[ϕlj[w[[i]] - w[[j]]],{j,1,i-1}], {i, 2, nWalkers}]

potE[walkers]

Evaluating the Cluster Size

The cluster size can be "measured" by the average radius from the center-of-mass (COM).

avgR[w_List] := Module[{com, δr},
com = (Plus @@ w)/Length[w];
δr = (# - com)& /@ w;
Plus @@ (Sqrt[# . #]& /@ δr)/Length[w]]

avgR[walkers]

Evaluating the change in potential

Δx = 0.1;

s[]:=RandomReal[{-Δx,Δx},3]

dpotE[w_List,k_Integer,step_] := Module[{rkOld,rkNew},
rkOld = w[[k]];
rkNew = rkOld + step;
Sum[If[i == k, 0,ϕlj[w[[i]] - rkNew] - ϕlj[w[[i]] - rkOld]], {i, 1, nWalkers}]]

Table[dpotE[walkers,RandomInteger[{1,nWalkers}],s[]],{5}]

A Random Walker Biased by Lennard-Jones Potential

We now have everything needed to define a random walker that can thermally agitate particles interacting through the Lennard-Jones potential.

step[w_]:=Module[{k,trialStep,ΔU},k=RandomInteger[{1,nWalkers}];trialStep=s[];If[(#1.#1&)[w[[k]]+trialStep]>(maxR^2),reject+=1;Return[w]];ΔU=dpotE[w,k,trialStep];If[ΔU<0,ReplacePart[w,w[[k]]+trialStep,k],If[RandomReal[]<(E^(-((ΔU)/T))),ReplacePart[w,w[[k]]+trialStep,k],reject+=1;w]]]

T = 0.1; 
Δx = 0.05; (* adjust to keep rejection rate < 50% *)
reject = 0;
w1 = NestList[step, walkers, 10];
w1 // Dimensions

reject

Run the simulation…

reject = 0;
(w1 = NestList[step,walkers,12000];Print["Rejection Rate = ", reject/12000 100.0 PerCent];) // Timing

Evaluate the total potential energy…

(potN = potE /@ w1[[Range[1,12000,nWalkers]]]) // Shallow

ListPlot[potN,PlotLabel->"Total Potential Energy",AxesLabel->{"n/N","∑ϕ"}]

Observe how the cluster size changes…

radN = avgR /@ w1[[Range[1,12000,nWalkers]]];

ListPlot[radN,PlotLabel->"Avg Cluster Radius",AxesLabel->{"n/N","⟨δr⟩"}]

cl1=ListPointPlot3D[w1[[1]],PlotStyle->{RGBColor[1,0,0],PointSize[0.02]},PlotRange->{{-maxR,maxR},{-maxR,maxR},{-maxR,maxR}}]

cl2=ListPointPlot3D[w1[[10000]],PlotStyle->{PointSize[0.02]},PlotRange->{{-(maxR/2),(maxR/2)},{-(maxR/2),(maxR/2)},{-(maxR/2),(maxR/2)}}]

A annimation of the final state…

ListAnimate[Table[ListPointPlot3D[w1[[n]],PlotStyle->PointSize[0.02],PlotRange->{{-(maxR/2),(maxR/2)},{-(maxR/2),(maxR/2)},{-(maxR/2),(maxR/2)}}],{n,12000-50 nWalkers,12000,nWalkers}], AnimationRunning-> False]

Ice

T = 0.05; 
Δx = 0.02; (* adjust to keep rejection rate < 50% *)

reject = 0;
(w2 = NestList[step,w1[[12000]],12000];Print["Rejection Rate = ", reject/12000 100.0 PerCent];) // Timing

potN = potE /@ w2[[Range[1,12000,nWalkers]]];

ListPlot[potN,PlotLabel->"Total Potential Energy",AxesLabel->{"n/N","∑ϕ"}]

radN = avgR /@ w2[[Range[1,12000,nWalkers]]];

ListPlot[radN,PlotLabel->"Avg Cluster Radius",AxesLabel->{"n/N","⟨δr⟩"}]

cl3=ListPointPlot3D[w2[[12000]],PlotStyle->{PointSize[0.02]},PlotRange->{{-(maxR/2),(maxR/2)},{-(maxR/2),(maxR/2)},{-(maxR/2),(maxR/2)}}]

Warming

T = 0.2; 
Δx = 0.05; (* adjust to keep rejection rate < 50% *)

reject = 0;
(w3 = NestList[step,w1[[12000]],6000];Print["Rejection Rate = ", reject/6000 100.0 PerCent];) // Timing

potN = potE /@ w3[[Range[1,6000,nWalkers]]];

ListPlot[potN,PlotLabel->"Total Potential Energy",AxesLabel->{"n/N","∑ϕ"}]

radN = avgR /@ w3[[Range[1,6000,nWalkers]]];

ListPlot[radN,PlotLabel->"Avg Cluster Radius",AxesLabel->{"n/N","⟨δr⟩"}]

Melting

T = 0.5; 
Δx = 0.1; (* adjust to keep rejection rate < 50% *)

reject = 0;
(w4 = NestList[step,w1[[12000]],6000];Print["Rejection Rate = ", reject/6000 100.0 PerCent];) // Timing

potN = potE /@ w4[[Range[1,6000,nWalkers]]];

ListPlot[potN,PlotLabel->"Total Potential Energy",AxesLabel->{"n/N","∑ϕ"}]

radN = avgR /@ w4[[Range[1,6000,nWalkers]]];

ListPlot[radN,PlotLabel->"Avg Cluster Radius",AxesLabel->{"n/N","⟨δr⟩"}]

Liquid

T = 1.0; 
Δx = 0.2; (* adjust to keep rejection rate < 50% *)

reject = 0;
(w5 = NestList[step,w1[[12000]],6000];Print["Rejection Rate = ", reject/6000 100.0 PerCent];) // Timing

potN = potE /@ w5[[Range[1,6000,nWalkers]]];

ListPlot[potN,PlotLabel->"Total Potential Energy",AxesLabel->{"n/N","∑ϕ"}]

radN = avgR /@ w5[[Range[1,6000,nWalkers]]];

ListPlot[radN,PlotLabel->"Avg Cluster Radius",AxesLabel->{"n/N","⟨δr⟩"}]

A annimation of the final state…

ListAnimate[Table[ListPointPlot3D[w5[[n]],PlotStyle->PointSize[0.02],PlotRange->{{-maxR,maxR},{-maxR,maxR},{-maxR,maxR}}],{n,6000-50 nWalkers,6000,nWalkers}], AnimationRunning-> False]

Summary

We've used the Rosenbluth-Teller-Metropolis algorithm to describe the equilibrium potential (or average binding energy) and equilibrium dimensions of a small cluster of atoms interacting by way of the Lennard-Jones potential. When the temperature of the system is low (T ≤ 0.2), then the cluster is tightly bound. As the temperature increases, the atoms within the cluster separate. Finally, for high temperatures (T ≥ 1), the atoms are no longer "clustered" but appear more like un-bound randomly-moving particles or pairs of particles.