# Demo Ising Plots1D-Computational Physics-Codes, Exercises for Computational Physics. Aligarh Muslim University

PDF (67 KB)
3 pages
1000+Number of visits
Description
Main topics for this course are Brownian dynamics, chaos, fluctuation, genetic algorithm, modelling and simulations, moments and variance, Monte Carlo modelling of neutron motion. This is code for assigned task. Its mai...
20 points
this document
Preview3 pages / 3

%Demo_IsingPlots1D - Generate plots for various

% thermodynamic quanitities in simulated 1D Ising

% Model and compare with exact solution

clear all; help Demo_IsingPlots1D;

tic

% Lattice size

lsize = [40 1 1];

% Boltzmann constant; easiest to set to 1

kb = 1;

%Particle spin

spin = 1/2;

% Strength of interaction

J = 1;

fprintf('Simulating 40 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 = 500;

% 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);

% Exact results for comparison

exactE = -spin^2*tanh(J*spin^2*beta_plots);

exactc_v = kb*(J*spin^2*beta_plots).^2./(cosh(J/4*beta_plots)).^2;

exactM = zeros(size(T_plots));

% 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, exactc_v, '--')

xlabel('T');

ylabel('c_v');

title('Specific heat v. temperature')

legend('Simulation','Exact Solution',0);

% Plot Magnetization

figure(3); clf;

plot(T, M_record, '.', T_plots, exactM, '--')

xlabel('T');

ylabel('M');

title('Average site magnetization v. temperature')

legend('Simulation','Exact Solution',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);