(* Student Name *)
(* July 30, 2018 *)
(* Physics III - Introduction to Statistical Physics *)
(* Assignment 6 - Problem 1*)

(a) Function Definitions

f[x_,y_]:=Sqrt[x^2+y^2]; (*Defines function*)
f[3,4] (*Computes function at value of x = 4*)

? ConstantArray

?Total

?Exp

?For

?RandomInteger

?RandomReal

?If

?AppendTo

?ListPlot

?Part

? Mean

?N

? Tanh

(b) Low Temperature Spin Simulations

(* Parameter Definitions *)
kBT =0.10; (*Temperature value*)
mu = 1.0 ; (* Magnetic dipole moment value *)
H = 1.0; (*Magnetic field value*)
Nspns = 1000; (*Number of spins in system*)
Nstps = 10000; (*Number of steps in simulation*)

(* Definition of initial Spin Configuration*)
S0 = ConstantArray[-1,Nspns];
Savg = {}; (* Empty vector representing total spin after each simulation *)

(* Energy  of Microststate S*) 
Energy[S_] := - mu*H*Total[S] ;

(* Boltzmann Factor Ratio; Use Mathematical operator "N" to make sure you get numerical value  *)
BoltzRatio[Sfinal_,Sinitial_]:=N[Exp[-(Energy[Sfinal] -Energy[Sinitial])/kBT]]; ;

(* Monte Carlo Simulation *)
For[i = 0, i< Nstps, i++, x = RandomInteger[{1,Nspns}]; u = RandomReal[]; AppendTo[Savg, N[Total[S0]/Nspns]];Snew = S0; Snew[[x]] = Snew[[x]]*(-1);
 If[u≤BoltzRatio[Snew,S0], S0 = Snew,S0 =S0 ]];

(* Plots the average spin of the system as a function of steps*)
ListPlot[Savg]

(* Selects final 1000 elements of Mtot and computes the average *)
N[Mean[ Part[Savg, Nstps-1000;;Nstps]]]

(*Calculates theoretical value of Mtot based on equilibrium statistical physics*)
Savgtheory[mu_, H_, kBT_] := N[Tanh[mu H/kBT]];
Savgtheory[mu, H, kBT]

(c) High Temperature Spin Simulations

(* Parameter Definitions *)
kBT = 10.0; (*Temperature value*)
mu = 1.0 ; (* Magnetic dipole moment value *)
H = 1.0; (*Magnetic field value*)
Nspns = 1000; (*Number of spins in system*)
Nstps = 10000; (*Number of steps in simulation*)

(* Definition of initial Spin Configuration*)
S0 = ConstantArray[-1,Nspns];
Savg2 = {}; (* Empty vector representing total spin after each simulation *)

(* Energy  of Microststate S*) 
Energy[S_] := - mu*H*Total[S] ;

(* Boltzmann Factor Ratio; Use Mathematical operator "N" to make sure you get numerical value  *)
BoltzRatio[Sfinal_,Sinitial_]:=N[Exp[-(Energy[Sfinal] -Energy[Sinitial])/kBT]]; ;

(* Monte Carlo Simulation *)
For[i = 0, i< Nstps, i++, x = RandomInteger[{1,Nspns}]; u = RandomReal[]; AppendTo[Savg2, N[Total[S0]/Nspns]];Snew = S0; Snew[[x]] = Snew[[x]]*(-1);
 If[u≤BoltzRatio[Snew,S0], S0 = Snew,S0 =S0 ]];

(* Plots the average spin of the system as a function of steps*)
ListPlot[Savg2]

(* Selects final 100 elements of Mtot and computes the average *)
N[Mean[ Part[Savg2, Nstps-1000;;Nstps]]]

(*Calculates theoretical value of Mtot based on equilibrium statistical physics*)
Savgtheory[mu_, H_, kBT_] := N[Tanh[mu H/kBT]];
Savgtheory[mu, H, kBT]

(d) We take the last 1000 steps in order to ensure the system has equilibrated (i.e., that it has entered a dynamical regime in which microstate probabilities are no longer time-dependent). We need to ensure that the system has equilibrated because our results  from statistical physics (in particular the average spin) only apply to equilibrium systems.

(e) We notice that the the average results computed from the simulation and the theoretical average are approximately equal to one another. The correspondence between the computational results from the low temperature simulation and the high temperature simulation is better because at low temperatures, the fluctuations around the equilibrium average are smaller than they are at high temperature. We can check this by computing the standard deviation of  the final 1000 elements in both parts. Doing so we find

N[StandardDeviation[ Part[Savg, Nstps-1000;;Nstps]]]
N[StandardDeviation[ Part[Savg2, Nstps-1000;;Nstps]]]

which affirms our claim that the standard deviation is lower for lower temperature. Indeed, the temperature is so low in part (a), that the fluctuations in the total spin are too small to be registered by Mathematica. 

If we wanted to improve the accuracy of our results (at the given temperatures) we could either consider more spins (i.e., increase Nspns) or run our simulation for longer (i.e., increase Nstps).

Optional Bonus

(* Simulation of Mean-Field Ising Model at two temperatures *)

(* Parameter Definitions *)
J= 1.0 ; (* Energy Coupling *)
Nspns = 1000; (*Number of spins in system*)
Nstps = 20000; (*Number of steps in simulation*)

(* Definition of initial Spin Configuration*)
S0 = ConstantArray[-1,Nspns];
Savg15 = {}; (* Empty vector representing total spin after each simulation, T = 1.5 *)
Savg05 = {}; (* Empty vector representing total spin after each simulation, T = 0.5 *)

(* Energy  of Microststate S*) 
Energy[S_] := -(J/(2* Nspns))*(Total[S])^2 ;

(* Boltzmann Factor Ratio; Use Mathematical operator "N" to make sure you get numerical value  *)
BoltzRatio[Sfinal_,Sinitial_, kBT_]:=N[Exp[-(Energy[Sfinal] -Energy[Sinitial])/kBT]]; ;

(* Monte Carlo Simulation at T = 1.5 *)
For[i = 0, i< Nstps, i++, x = RandomInteger[{1,Nspns}]; u = RandomReal[]; AppendTo[Savg15, N[Total[S0]/Nspns]];Snew = S0; Snew[[x]] = Snew[[x]]*(-1);
 If[u≤BoltzRatio[Snew,S0, 1.5], S0 = Snew,S0 =S0 ]];

(* Monte Carlo Simulation at T = 0.5 *)
For[i = 0, i< Nstps, i++, x = RandomInteger[{1,Nspns}]; u = RandomReal[]; AppendTo[Savg05, N[Total[S0]/Nspns]];Snew = S0; Snew[[x]] = Snew[[x]]*(-1);
 If[u≤BoltzRatio[Snew,S0, 0.5], S0 = Snew,S0 =S0 ]];

(* Plots the average spin of the system as a function of steps*)
ListPlot[Savg15]
ListPlot[Savg05]

(* Selects final 100 elements of Mtot and computes the average *)
N[Mean[ Part[Savg15, Nstps-1000;;Nstps]]]
N[Mean[ Part[Savg05, Nstps-1000;;Nstps]]]

We see that the simulation that is at a temperature T = 1.5 produces an average spin that is close to zero, but the simulation that is at the temperature T=0.5 produces an average spin which is close to 1. The temperature T = 1.0, is the critical temperature of this system so we expect (and we see) that above this temperature, the system has roughly zero average spin.