Download 2D Ising Model Simulation: Comparing with Exact Solution - Energy, Heat, Magnetization and more Exercises Computational Physics in PDF only on Docsity!
%Demo_IsingPlots2D - Generate plots for various % thermodynamic quanitities in simulated 2D Ising % Model and compare with exact solution
clear all; help Demo_IsingPlots2D;
tic
% Lattice size lsize = [20 20 1];
% Boltzmann constant; easiest to set to 1 kb = 1;
%Particle spin spin = 1/2;
% Strength of interaction J = 1;
fprintf('Simulating 20x20 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 = 2spin + 1; lattice = -spin + floor(rand(lsize)nspins); E = totalEnergy(lattice, J, Bz);
beta=1./(kbT); beta_plots = 1./(kbT_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 = 500;
% Give the system some time to reach equilibrium for istep=1: [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);
% Exact results for comparison Tc = 2.269185Jspin^2/kb; for istep = 1:length(T_plots) z = beta_plots(istep)Jspin^2; K = 2/(cosh(2z)coth(2*z));
EIntFun = inline('sin(x).^2./(sqrt(1-K.^2.sin(x).^2).(1+sqrt(1- K.^2.*sin(x).^2)))','x','K');
exactE(istep) = -2Jspin^2tanh(2z) +... K/(2pi)-8exp(2z)(exp(8z)- 6exp(4z)+1)Jspin^2/(exp(4z)+1)^3... quadl(EIntFun,0,pi,[],[],K); end exactc_v = diff(exactE)/(T_plots(2) - T_plots(1)); exactM = spin(0 . (T_plots >= Tc) + ... (1 - (sinh(2beta_plotsJspin^2)).^-4).^(1/8) . (T_plots < Tc));
% Plot Energy figure(1); clf; plot(T, E_record, '.', T_plots, exactE, '--') xlabel('T'); ylabel('E'); title('Average site energy v. temperature') legend('Simulation','Exact Solution',0);
% Plot c_v figure(2); clf; plot(T(1:end-1), c_v, '.', T_plots(1:end-1), exactc_v, '--')
