Monte Carlo simulation of Leonnard-Jones fluid

List of variables

Config = list of particle positions

boxL = system linear dimension

nTot = total number of particles

bete = beta epsilon

rMax = truncation range

nCell = number of cells in one direction

wid = width of a cell

cellist = 2d lists of particles in cells

eTotal = total energy

pExcess = excess pressure

iSave = index into eValues and pValues

eValues = list of stored energy values

pValues = list of stored excess pressure values

 

List of functions

1) Leonnard-Jones interaction energy 

rMax = 2.0;

iLJtrunc = 1/rMax^12 - 2/rMax^6;

iLJ[r_] := (1/r^12 - 2/r^6) – iLJtrunc

 

2) List of neighbor separations

nearestImage[pos_] := Mod[pos, boxL, -boxL/2];

neighborSep[k_, pos_] := Module[{i, j, sepVecs}, {i, j} = cell[pos];

neighbors = Table[cellList[[Mod[i1, nCell, 1], Mod[j1, nCell, 1]]], {i1, i - 1, i + 1}, {j1, j - 1, j + 1}] // Flatten // DeleteCases[#, k] &;

sepVecs = nearestImage[# - pos] & /@ config[[neighbors]];

Select[Norm /@ sepVecs, (# < rMax) &] ]

 

3) Total energy between all pairs of particles

singleE[k_, pos_] := Total[iLJ /@ neighborSeps[k, pos]];

totalE[] := (1/2) Sum[singleE[i, config[[i]]], {i, 1, nTot}]

*delta energy upon move

deltaE [iMove_, newPos_] := singleE[iMove, newPos] - singleE[iMove, config[[iMove]]]

 

4) excess pressure

fLJ[r_] := 12 (1/r^13 - 1/r^7)

singleP[k_, pos_] := Total[# fLJ[#] & /@ neighborSeps[k, pos]]/(2 boxL^2)

excessP[] := (1/2) Sum[singleP[i, config[[i]]], {i, 1, nTot}] 

*delta excess pressure upon move

deltaexcessP[iMove_, newPos_] := singleP[iMove, newPos] - singleP[iMove, config[[iMove]]]

 

Simulation tutorial

1) Initialize simulation parameters

initSim[nParticles_, areaF_, be_] := Module[{},

  bete = be;

  nTot = nParticles;

  boxL = Sqrt[nTot/areaF];

  range = boxL;

  config = initConfig[nTot];

  nCell = Floor[boxL/rMax];

  wid = boxL/nCell;

  buildCellList[];

  eTotal = totalE [];

pExcess = excessP [];]

 

2) Create initial configurations

initConfig[n_] := Module[{},

nside = Ceiling[Sqrt[n]];

dx = boxL/(nside + 1);

Take[Flatten[dx Table[{i, j}, {i, nside}, {j, nside}], 1], n]]

 

3) Set up cell lists

cell[pos_] := Ceiling[pos/wid];

buildCellList[] := Module[{i, j, k},

cellList = Table[{}, {i, 1, nCell}, {j, 1, nCell}];

Do[{i, j} = cell[config[[k]]];

   cellList[[i, j]] = Append[cellList[[i, j]], k],

   {k, 1, nTot}]]

removeFromCellList[k_] := Module[{i, j},

{i, j} = cell[config[[k]]];

cellList[[i, j]] = DeleteCases[cellList[[i, j]], k]]

addToCellList[k_] := Module[{i, j},

{i, j} = cell[config[[k]]];

cellList[[i, j]] = Append[cellList[[i, j]], k]]

 

4) Monte Carlo move

* First, we randomly select a particle index iMove to move, find the old position oldLoc of the particle. Then, we randomly select a new position newLoc. We also use deltaEnergy[ ] to compute the energy change dE if the particle were moved. Last, we use the Metropolis criterion to determine if succeed is true (if the move succeeds).*

attemptMove[] := Module[{newLoc, iMove, dE, succeed},

iMove = RandomInteger[{1, nTot}];

newLoc = boxL {RandomReal[], RandomReal[]};

dE = deltaE [iMove, newLoc];

succeed = If[dE < 0, True, RandomReal[] < Exp[-bete dE]];

 

*If the move succeeds, we should update the total energy eTotal and use deltaPExcess[ ] to update the excess pressure pExcess. Then, we use removeFromCellList[ ] to remove the particle from the cell list, move the particle to its new location, and use addToCellList[ ] to add the particle back to the cell list.*

If[succeed,

  eTotal += dE;

  pExcess += deltaexcessP[iMove, newLoc];

  removeFromCellList[iMove];

  config[[iMove]] = newLoc;

addToCellList[iMove];];]

 

runSimulation[nMonteCarlosteps_] := Module[{},

iSave = 0;

eValues = Array[0 &, nTot* nMonteCarlosteps];

pValues = Array[0 &, nTot* nMonteCarlosteps];

Do[ attemptMove[];

   eValues[[++iSave]] = eTotal;

   pValues[[iSave]] = pExcess,

   {istep, 1, nTot nMonteCarlosteps }] ]

 

drawConfig[] := Graphics[Disk[#, 0.5] & /@ config, Frame -> True,

PlotRange -> {{-0.5, boxL + 0.5}, {-0.5, boxL + 0.5}},

AspectRatio -> 1]

 

 

5) Obtain average energy and pressure

aveE = Mean[eValues]

aveP = Mean[pValues]

 

6) Determine Cv and Z in units of kB and εσ2, respectively

* Ideal heat capacity is NK, which in reduced units is 1. The excess heat capacity is correction to the ideal one. Ideal gas pressure is P = N kT/L2. In reduced units, P/P* is just areaF/bete. The excess pressure is correction to the ideal one.*

excessCv=Variance[eVaules] bete^2 / nTot;

Cv=1+excessCv

Z = (areaF/bete + aveP)/(areaF/bete)