On Fri, Jul 3, 2009 at 7:14 PM, Carletto Rossi
<address@hidden> wrote:
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
Thanks in advance
_______________________________________________
Help-octave mailing list
address@hidden
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave