Semester Project
Molecular Dynamics Simulation
Solid-Liquid Phase Diagram of Argon

Name: 			Gan Sik Hong & Hwang Hsien Shiung
Matrics Number:	    105513         &           105762

Abstract

This project’s aim is to produce the solid-liquid phase diagram of Argon.

The method employed is three dimensional molecular dynamics simulation.

The simulation begins with setting up the solid phase for Argon with fcc structure for particular density of box.
Then, the system is heated manually and repeatedly until liquid phase is achieved.
Temperature and pressure at phase transition are recorded.
The process is repeated for different density of the box, i.e., different pressure.

Few points (pressure versus temperature at transition) are obtained and phase diagram is plotted.
Phase diagram obtained from simulation is compared with experimental data.

Discussion:

Classical approach is chosen, where Newton second law will be used as the core of simulation.
F=ma

However, there is still quantum element within this simulation, where Lennard-Jones potential is employed to determine the interaction potential and force between two particles.
V(r)=4ϵ[((((σ/r)))^12)-((((σ/r)))^6)]

Verlet method is used for time evolving the particle position because it can conserve energy better than Euler type methods.
(Subscript[x, i])(n+1)=2(Subscript[x, i])(n)-(Subscript[x, i])(n-1)+(Subscript[a, i,x])(n)(((Δt))^2)   (same goes for y and z)

Reduced units are used. So, in the end when plotting the phase diagram, the real unit have to be restored.
For reduced unit, we set σ, ϵ, m to 1.
To restore the unit, for pressure it will be multiplied by (ϵ/(σ^3))= 420 bar (the original unit is Pascal, but we convert it to bar)
for temperature it will be multiplied by (ϵ/((Subscript[k, B])))= 120 K

Periodic boundary condition is used to eliminate the effect of wall. Besides that, PBC allow us to obtain the effect of large amount of particles simulation based on only small amount of particles.
In order to speed up calculation, the pair interaction is considered zero for separation larger than 3.2σ. This is called cutting off potential.

Equipartition theorem is used to calculate the temperature.
KE=(3/2)NkT

Virial equation is used to calculate the pressure.
 D is the dimension, in our case, D=3



Figures above show the fcc setup for initial position of the particles of solid Argon. 


Without any heating, we can see from above figure that the total energy goes to equilibrium very fast, at about 100 time step.
So, to ensure the system has really went into equilibrium, we set the time step needed for equilibrium to be 200.

Heating is done by adding kinetic energy to the system manually, through rescaling the velocity.
Mean square displacement, MSD is the parameter employed to measure the solidness and liquidness of the system.

If the system is solid, MSD saturates to a finite value, as shown below.
The particles will have very small vibration, which reflect how atoms behave in solid.


if the system is liquid, MSD grows linearly with time.

The figure above shows one of the MSD versus time curve for heating processes that successfully changes solid into liquid phase.
We can see that the transition to liquid begin at around 500 time step, where MSD starts to grow linearly with time.
So we take the range from 0.10 till 0.15 as an indicator of phase transition.
As the system is heated even more, we can see that the diffusion occurs more rigorously, where increment of MSD becoming larger. (the visualisation is more clear with manipulate plot)

The phase diagram obtained is shown in figure below:

Red line is the experiment data, blue line is the result of simulation.
We can see that this program manage to produce the phase diagram, although it have some deviation from the experimental data.
This shows that the program coding successfully employed the theory and equations.

However, this program is still not stable and may not produce the same result even if same parameters are used for repeated run.

There are few suggestions to improve this program:
1) Use smaller time step and let the program run longer. This will reduce the errors but require longer run time.
2) Use smaller R to rescale the velocity. This will produce a more accurate result but require longer run time.
3) Employed a better technique to capture the transition phase. One of the possible technique is observe the speed distribution of the particles.
4) Simulate with more number of particle. The program will be more accurate with increasing number of particles. However, this will certainly and massively extend the run      
     time. The next number of particle suitable for simulation is 4((3^3)) = 108, which is about three times more than the particle used now (32 particles).

When the program is more stable, it is better to combine part 1 and part 2 together, so that it is not necessary to export and import the data.
Furthermore, another loop for manipulating the pressure can be added, so that user doesn’t have to manually manipulate the pressure.
However, this should only be done if the program is really stable. Otherwise, the program may have took more time but still produce an undesired result.

For more generic code, i.e. solving for other elements, we need to know few essential details of that particular element, such as the interaction potential, σ, m, ϵ.

First Part: Data collecting

SetDirectory["C:\\Users\\WestWind\\Desktop\\data"]; (* set the directory to store the data *)
Needs["PlotLegends`"]; (* required to assign legend for plots *)

(* value to be changed manually for each data collection *)
np=10; (* the np-th data collection *) (* suggestion: start with 1 and increase in step of 1 until 5-6 data obtained *)
ρ=1.9; (* density of box *) (* suggestion: choose ρ between 1.8 to 2.2 *)
R=3.0; (* scaling factor to rescale velocity, i.e. heat the system *) (* might need to change accordingly, different ρ might need different R *)
(* value to be changed manually for each data collection *)


n=32; (* number of particles, must follow the rule n = 4(M^3), where M = 2,3,... *)
L=((n/ρ)^(1/3)); (* length of each edge of box *)

σ=1; (* constant that set the distance scale *)
ϵ=1;(* constant that set the energy scale *)
m=1;(* constant that set the mass scale *)

δr=0.05; (* maximum shifting allowed for initial configuration setting *)
Δt=0.005; (* time step *)

rcut=3.2; (* distance to cut off potential *)
ttMSD={}; (* table to record mean square displacement at each instant *)

nrise=1;(* the first time step to rise temperature *)
nrisetotal=6; (* total number of heating + 1 *)

ntrans=0; (* to record the nrise for transition stage *)

ne=300;(* time step required for equilibrium *)
nm=100;(* time step for measurement *)

tlast=(ne+nm)*nrisetotal; (* total time step *)

(* Preparing initial configuration *)
pos[0]={}; (* record the initial position *)

M=1;
While[4*(M^3)<n,
M++;
];

a=(L/M); (* lattice constant *)

(* configuration for the basis of fcc *)
xCell[1]=0.25;xCell[2]=0.75;xCell[3]=0.75;xCell[4]=0.25;
yCell[1]=0.25;yCell[2]=0.75;yCell[3]=0.25;yCell[4]=0.75;
zCell[1]=0.25;zCell[2]=0.25;zCell[3]=0.75;zCell[4]=0.75;
(* fcc configuration of basis *)

(* assigning all the n particles at respective basis locations *)
i=1;
For[xx=0,xx<M,xx++,
For[yy=0,yy<M,yy++,
For[zz=0,zz<M,zz++,
For[k=1,k≤4,k++,
xini[i]=(xx+xCell[k])*a;
yini[i]=(yy+yCell[k])*a;
zini[i]=(zz+zCell[k])*a;
i++;
];
];
];
];
(* end of assigning all the n particles at respective basis locations *)

For[i=1,i≤n,i++,

(* random shifting of initial configuration *)
x[i,0]= xini[i]+2*(Random[]-0.5)*δr;
y[i,0]=yini[i]+2*(Random[]-0.5)*δr;
z[i,0]=zini[i]+2*(Random[]-0.5)*δr;
(* end of random shifting of initial configuration *)

(* set the pre-initial velocity *)
vx[i,-1]=0;
vy[i,-1]=0;
vz[i,-1]=0;
(* end of set the pre-initial velocity *)

(* calculating pre-initial position, required for Verlet algorithm *)
x[i,-1]=x[i,0]-vx[i,-1]*Δt;
y[i,-1]=y[i,0]-vy[i,-1]*Δt;
z[i,-1]=z[i,0]-vz[i,-1]*Δt;
(* end of calculating pre-initial position *)

pos[0]=Append[pos[0],{x[i,0],y[i,0],z[i,0]}];

(* initialising the correction factor for calculating MSD *)
xf[i]=0; 
yf[i]=0;
zf[i]=0;
(* end of initialising the correction factor for calculating MSD *)

]; (* end for i *)
(* end of preparing initial configuration *)


Print["Time stepping:",ProgressIndicator[Dynamic[t],{0,tlast}],"(",Dynamic[t],"/",tlast,")"];
Print["i iterating:",ProgressIndicator[Dynamic[i],{1,n+1}],"(",Dynamic[i],"/",n+1,")"];
Print["j iterating:",ProgressIndicator[Dynamic[j],{1,n+1}],"(",Dynamic[j],"/",n+1,")"];

(* commence time stepping *)
For[t=0,t<tlast,t++,
pos[t+1]={};(* to record the position of all particle at the instant *)
PEsumi[t]=0;(* sum of PE of all i, will avoid double counting later *)
pcomp[t]=0;(* component for calculating pressure *)
sumΔrsq=0;(* sum of square displacement *)

(* iterate through n particles *)
For[i=1,i≤n,i++,
ax=0; (* resetting x component of total force acting upon i at time t *)
ay=0;(* resetting y component of total force acting upon i at time t *)
az=0;(* resetting z component of total force acting upon i at time t *)

sumΔrsq=sumΔrsq+(((x[i,t]+xf[i]-x[i,0]))^2)+(((y[i,t]+yf[i]-y[i,0]))^2)+(((z[i,t]+zf[i]-z[i,0]))^2);

For[j=1,j≤n,j++,
If[j≠ i, (* exclude calculation of force that particle exert upon itself *)

(* calculating separation between particle i and j *)
Δxtest[1]=x[i,t]-x[j,t];
Δytest[1]=y[i,t]-y[j,t];
Δztest[1]=z[i,t]-z[j,t];

(* checking for periodic boundary condition, minimum separation rule *)
Δxtest[2]= x[i,t]-x[j,t]+L; 
Δytest[2]= y[i,t]-y[j,t]+L; 
Δztest[2]= z[i,t]-z[j,t]+L; 
Δxtest[3]=x[i,t]-x[j,t]-L;
Δytest[3]=y[i,t]-y[j,t]-L;
Δztest[3]=z[i,t]-z[j,t]-L;
(* end of checking for periodic boundary condition, minimum separation rule *)

(* take the separation with smallest magnitude *)
Δx=Min[Abs[Δxtest[1]],Abs[Δxtest[2]],Abs[Δxtest[3]]];
Δy=Min[Abs[Δytest[1]],Abs[Δytest[2]],Abs[Δytest[3]]];
Δz=Min[Abs[Δztest[1]],Abs[Δztest[2]],Abs[Δztest[3]]];

(* restoring the negative sign if it had it *)
For[k=1,k≤3,k++,
If[Δx==-Δxtest[k],Δx=-Δx];
If[Δy==-Δytest[k],Δy=-Δy];
If[Δz==-Δztest[k],Δz=-Δz];
];(* end for k *)
(* end of restoring the negative sign *)

(* end of calculating separation between particle i and j *)

(* calculating force and potential *)
rji=Sqrt[((Δx)^2)+((Δy)^2)+((Δz)^2)];

If[rji≤rcut, (* cutting off potential for rji larger than rcut *)
ax=ax+24*((2/(rji^14))-(1/(rji^8)))*Δx;
ay=ay+24*((2/(rji^14))-(1/(rji^8)))*Δy;
az=az+24*((2/(rji^14))-(1/(rji^8)))*Δz;

If[j>i, (* avoid double counting *)
PEsumi[t]=PEsumi[t]+4*((1/(rji^12))-(1/(rji^6)));
pcomp[t]=pcomp[t]+rji*24*((2/(rji^13))-(1/(rji^7)));
];(* end if for avoiding double counting *)

]; (* end if for cutting off potential *)
(* end of calculating force and potential *)

]; (* end if for excluding calculation of force that particle exert upon itself *)

]; (* end for j *)
(* end of iterate through other particles *)

(* Verlet method *)
x[i,t+1]=2*x[i,t]-x[i,t-1]+ax*((Δt)^2);
y[i,t+1]=2*y[i,t]-y[i,t-1]+ay*((Δt)^2);
z[i,t+1]=2*z[i,t]-z[i,t-1]+az*((Δt)^2);
(* end of Verlet method *)

(* rescaling the velocity to heat the system *)
If[t==IntegerPart[(tlast/nrisetotal)]*nrise,
x[i,t]=x[i,t+1]-R*(x[i,t+1]-x[i,t]);
y[i,t]=y[i,t+1]-R*(y[i,t+1]-y[i,t]);
nrise++;
];

vx[i,t]=((x[i,t+1]-x[i,t])/(2*Δt));
vy[i,t]=((y[i,t+1]-y[i,t])/(2*Δt));
vz[i,t]=((z[i,t+1]-z[i,t])/(2*Δt));

(* end of rescaling the velocity to heat the system *)

(* checking PBC to see if new position from Verlet method is out of the box, shift it accordingly if yes, and shift its position in previous time step too for calculation at next time step *)
If[x[i,t+1]>L,x[i,t+1]=x[i,t+1]-L; x[i,t]=x[i,t]-L;xf[i]=xf[i]+L];
If[x[i,t+1]<0,x[i,t+1]=x[i,t+1]+L; x[i,t]=x[i,t]+L;xf[i]=xf[i]-L];
If[y[i,t+1]>L,y[i,t+1]=y[i,t+1]-L; y[i,t]=y[i,t]-L;yf[i]=yf[i]+L];
If[y[i,t+1]<0,y[i,t+1]=y[i,t+1]+L; y[i,t]=y[i,t]+L;yf[i]=yf[i]-L];
If[z[i,t+1]>L,z[i,t+1]=z[i,t+1]-L; z[i,t]=z[i,t]-L;zf[i]=zf[i]+L];
If[z[i,t+1]<0,z[i,t+1]=z[i,t+1]+L; z[i,t]=z[i,t]+L;zf[i]=zf[i]-L];

pos[t+1]=Append[pos[t+1],{x[i,t+1],y[i,t+1],z[i,t+1]}]; (* record all the new positions at time t+1 *)

(* give warning and quit kernal if the particle fly too far out of the box *)
If[x[i,t+1]>2*L||y[i,t+1]>2*L||z[i,t+1]>2*L,Print["Warning, got particle fly too far"];Quit[];];

];(* end for i *) 
(* end of iterate through n particles *)

If[ ((sumΔrsq)/n)≥0.1&&((sumΔrsq)/n)≤0.15,  (* capture transition stage *)
ntrans=nrise-1;
];
ttMSD=Append[ttMSD,{t,((sumΔrsq)/n)}];

]; (* end of time stepping *)

(***** calculation *****)
ttT={}; (* table to record temperature at each instant *)
ttKE={};(* table to record kinetic energy at each instant *)
ttPE={};(* table to record potential energy at each instant *)
ttE={};(* table to record total energy at each instant *)
ttp={}; (* table to record pressure at each instant *)


For[t=0,t<tlast,t++,
KEsumi[t]=0;(* resetting sum of KE for all i *)

For[i=1,i≤n,i++,
KEsumi[t]=KEsumi[t]+(m/2)*(((vx[i,t])^2)+((vy[i,t])^2)+((vz[i,t])^2));(* sum of KE for all i at particular instant *)

];(* end for i *)


instantT[t]=(2/3)*((KEsumi[t])/n);(* instantaneous temperature *)

Etotalaveri[t]=((KEsumi[t]+PEsumi[t])/n) ;(* total energy averaged over all i at time step t *)

pressure[t]=((n*instantT[t]+(1/(3*n))*pcomp[t])/(L^3));(* pressure at time step t *)



ttT=Append[ttT,{t,instantT[t]}];
ttKE=Append[ttKE,{t,((KEsumi[t])/n)}];
ttPE=Append[ttPE,{t,((PEsumi[t])/n)}];
ttE=Append[ttE,{t,Etotalaveri[t]}];
ttp=Append[ttp,{t,pressure[t]}];

];(* end for t *)

(* calculating the transition temperature and pressure *)
If[ntrans≠ 0,
ncal=0;
sumTcal=0;
sumTcalsq=0;
sumpcal=0;
sumpcalsq=0;

For[t=(ne+nm)*ntrans+ne,t<(ne+nm)*(ntrans+1),t++, (* start calculating after it reach equilibrium and before the next stage *)
sumTcal=sumTcal+instantT[t]; (* sum of all instantaneous temperature at the stage *)
sumTcalsq=sumTcalsq+((instantT[t])^2);(* sum of all square of instantaneous temperature at the stage *)
sumpcal=sumpcal+pressure[t];(* sum of all instantaneous pressure at the stage *)
sumpcalsq=sumpcalsq+((pressure[t])^2);(* sum of all square of instantaneous pressure at the stage *)
ncal++ (* to calculate the number of time step in the stage *)
];(* end for t, for summing data at transition stage *)

Tcal=(sumTcal/ncal);(* mean temperature at transition *)
Tcalsq=(sumTcalsq/ncal);
σTcal=Sqrt[Tcalsq-(Tcal^2)];(* standard error of temperature at transition *)
pcal=(sumpcal/ncal);(* mean pressure at transition *)
pcalsq=(sumpcalsq/ncal);
σpcal=Sqrt[pcalsq-(pcal^2)];(* standard error of pressure at transition *)

];(* end if *)
(* end of calculating the transition temperature and pressure *)

(***** end of calculation *****)

ListPlot[{ttKE,ttPE,ttE},PlotRange->All,Joined->True,AxesLabel->{"Time step","Energy"},PlotLegend->{"KE","PE","E"},LegendPosition->{1.1,-0.4},ImageSize->700]

ListPlot[ttE,PlotRange->All,Joined->True,AxesLabel->{"Time step","Total Energy"}]
ListPlot[ttKE,PlotRange->All,Joined->True,AxesLabel->{"Time step","Kinetic Energy"}]
ListPlot[ttT,PlotRange->All,Joined->True,AxesLabel->{"Time step","Temperature"}]
ListPlot[ttp,PlotRange->All,Joined->True,AxesLabel->{"Time step","pressure"}]
ListPlot[ttMSD,PlotRange->All,Joined->True,AxesLabel->{"Time step","MSD"}]

Manipulate[
ListPointPlot3D[pos[t],PlotRange->{{0,L},{0,L},{0,L}},Axes->False,AspectRatio->1,PlotLabel->Style[Framed["MD in 3D"],Background->Lighter[Yellow]]],
{t,0,tlast,1}
]

If[ntrans==0,
Print["rescaling failed to heat the system into liquid state"]
,
Print["Transition temperature = ",Tcal," ± ", σTcal," at pressure of ",pcal, " ± ",σpcal];
];

Print["Time used = ",((TimeUsed[])/60)," minutes" ];

EmitSound[Sound[SoundNote[]]];

(* exporting data for plotting the phase diagram *)
filename="Tp"<>ToString[np]<>".dat";
Export[filename,{Tcal,pcal},"Data"];

filename2="errorT"<>ToString[np]<>".dat";
Export[filename2,σTcal,"Data"];

filename3="errorp"<>ToString[np]<>".dat";
Export[filename3,σpcal,"Data"];
(* end of exporting data for plotting the phase diagram *)

Second Part: Phase diagram plotting

SetDirectory["C:\\Users\\WestWind\\Desktop\\data"]; (* set the directory to import the data *)
Needs["ErrorBarPlots`"];  (* needed for errorlistplot *)

np=8 ; (* specify the number of points available to be imported *)

Tunit=120 ; (* to restore the original unit of temperature, unit: Kelvin *)
Punit=420; (* to restore the original unit of pressure, unit: bar *)

(* Importing data *)
For[n=1,n≤np,n++,
filename="Tp"<>ToString[n]<>".dat";
point[n]=Flatten[Import[filename,"Data"]];
point[n]={point[n][[1]]*Tunit,point[n][[2]]*Punit};

filename2="errorT"<>ToString[n]<>".dat";
σT[n]=Import[filename2,"Data"];
σT[n]=(σT[n][[1]][[1]])*Tunit;

filename3="errorp"<>ToString[n]<>".dat";
σp[n]=Import[filename3,"Data"];
σp[n]=(σp[n][[1]][[1]])*Punit;
];(* end for *)
(* end of importing data *)

tTp=Table[{point[n],ErrorBar[σT[n],σp[n]]},{n,1,np}];
L1=ErrorListPlot[tTp,Joined->True];

tTptheory={{100,641},{150,2906},{184,5000},{200,5983},{250,9487},{300,13504},{350,17863},{400,22521}};
L2=ListPlot[tTptheory,AxesLabel->{"Temperature(Kelvin)","Pressure(bar)"},Joined->True,PlotStyle->Red,PlotRange->All];

Show[L2,L1]

Conclusion

Phase diagram is successfully obtained (although there are deviations)

This means that the program is able to achieve the following:
1) prepare to setup for particles in Argon solid phase (fcc)
2) utilise the relevant equations (Newton’s second law, LJ potential, Verlet algorithm, MSD etc..)
3) vary the temperature of system (heating)
4) vary the pressure of system
5) capture the stage where phase transition occured

Basically, the program is considered as successful.
The only thing left is refinement of the coding (refinement of parameters used such as time step, scaling factor, etc)