Hi guys. I would to convert this Matlab script in an Octave file:
N = 4; %number of particles n = 1000; %size of cell iteration = 500; %times of measurement step = 20; %number of steps per-iteration T = 1.0; %temperature
K = 1.0; %Bolztmann constant x = zeros(N,1); y = zeros(N,1); p = 1; too_near = 0; % ###################### initiation ############################ while(p <= N) %generate trial position x_trial = rand*n;
y_trial = rand*n; for q = 1:p d = sqrt( (x_trial-x(q))^2 + (y_trial-y(q))^2 ); if(d <= 1.0) %the new position is too "near" too_near = 1; break; %there's no need to continue the loop
else too_near = 0; end
end if(too_near == 0) %the new position is "far", accept it x(p) = x_trial; y(p) = y_trial; p = p+1; %step increment
end end % ######################### metropolis ############################# for k=1:iteration W_initial = calculate_potential(N,x,y); for j=1:step %how many steps per-iteration for i=1:N %for all particles
%step range between [-0.5,0.5] trial_x = rand-0.5; trial_y = rand-0.5; %save the old values lest the step is not accepted x_temp = x(i); y_temp = y(i);
%trial step, bounded by the periodic boundary condition x(i) = x(i) + trial_x; if (x(i) > n) x(i) = x(i) - n; elseif (x(i) < 0.0) x(i) = n + x(i);
end y(i) = y(i) + trial_y; if (y(i) > n) y(i) = y(i) - n; elseif (y(i) < 0.0) y(i) = n + y(i); end
W_final = calculate_potential(N,x,y);
if( (W_final-W_initial) < 0.0 ) %move is accepted W_initial = W_final; else %apply metropolis criterion if (rand < exp(W_initial-W_final)/(K*T)) W_initial = W_final;
else %move is not accepted, return to initial position x(i) = x_temp; y(i) = y_temp; end end end end scatter(x,y,3,[0 0 1], 'filled') %plot 2D scatter-graph
F(k) = getframe; %capture every plot end save simulation.mat F; movie2avi(F, 'simulation.avi'); % ################# calculate potential ################### function E = calculate_potential(N,x,y)
E = 0;
for i=1:N %for all particles for j=(i+1):N %compute the potential of j with respect to i d = sqrt( (x(j)-x(i))^2 + (y(j)-y(i))^2); U = d^(-12) - d^(-6); %Lennard-Jones potential
E = E + U; end end