## Search in the document preview

%Demo_IsingPlots3D - Generate plots for various

% thermodynamic quanitities in simulated 3D Ising

% Model

% Warning: This can take hours to run.

clear all; help Demo_IsingPlots3D;

tic

% Lattice size

lsize = [20 20 20];

% Boltzmann constant; easiest to set to 1

kb = 1;

%Particle spin

spin = 1/2;

% Strength of interaction

J = 1;

fprintf('Simulating 20x20x20 spin 1/2 particles with J = 1.\n')

Bz = input('Input external magnetic field strength: ');

% Range of temperatures

T = linspace(2,0,25);

T_plots = linspace(2,0);

% Random Lattice

nspins = 2*spin + 1;

lattice = -spin + floor(rand(lsize)*nspins);

E = totalEnergy(lattice, J, Bz);

beta=1./(kb*T);

beta_plots = 1./(kb*T_plots);

for tstep = 1:length(T)

% Status update

if (mod(tstep,2) == 1)

fprintf('[%g s] Now simulating T = %g (%g%%

done)\n',toc,T(tstep),100*(tstep-1)/length(T));

end

% Number of steps to average over

maxsteps = 200;

% Give the system some time to reach equilibrium

for istep=1:10

[lattice, dE] = MCStep(lattice, beta(tstep), J, Bz, spin);

E = E + dE;

end

E_ct = 0;

M_ct = 0; % Cumulative total for averaging

E_sq = 0;

M_sq = 0; % Cumulative total of squares for std dev

for istep=1:maxsteps

% Take a Monte Carlo step

[lattice, dE] = MCStep(lattice, beta(tstep), J, Bz, spin);

% Update the energy, magnetization, and running totals

E = E + dE;

E_ct = E_ct + E;

E_sq = E_sq + E^2;

M = avgMagnetization(lattice);

M_ct = M_ct + M;

M_sq = M_sq + M^2;

end

% Record the energy, magnetization, and fluctuations in magnetization

E_record(tstep) = E_ct/(istep*prod(lsize));

M_record(tstep) = M_ct/istep;

dM_record(tstep) = ((M_sq/istep) - (M_ct/istep)^2)^(1/2);

end

% Record the specific heat

c_v = diff(E_record)/(T(2)-T(1)); % c_v = d(total E)/d(T)

fprintf('[%g s] Done!',toc);

% Plot Energy

figure(1); clf;

plot(T, E_record, '.')

xlabel('T');

ylabel('E');

title('Average site energy v. temperature')

legend('Simulation',0);

% Plot c_v

figure(2); clf;

plot(T(1:end-1), c_v, '.')

xlabel('T');

ylabel('c_v');

title('Specific heat v. temperature')

legend('Simulation',0);

% Plot Magnetization

figure(3); clf;

plot(T, M_record, '.')

xlabel('T');

ylabel('M');

title('Average site magnetization v. temperature')

legend('Simulation',0);

% Plot fluctuations in magnetization

figure(4); clf;

plot(T, dM_record, '.')

xlabel('T');

ylabel('\Delta{}M');

title('Fluctuations of the average site magnetization v. temperature')

legend('Simulation',0);