Demo Ising Plots3D-Computational Physics-Codes, Exercises for Computational Physics

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