Can anyone explain how can one explain ensemble average concept as i have the working code in Scilab below which demonstrates particles experiencing attraction and repulsion when they comes e=to each other (like the LJ potential) and one sees that the speed distribution obtained in my code matches with Maxwell Boltzmann curve.
Take N=25 tf=25 dt=0.02 and V0=1
function r=randint(a, b) //Generating random number for any interval for given value of a and b.
r = a + (b-a)*rand();
endfunction
function [x0,y0,L]= initialise(N,V0,dt,tf) //Initialising variables
//N represents number of Argon atoms
//V0 represents initial velocity of Argon atoms
//tf represents final time of simulation
// dt step size to increment in time
L=round(sqrt(N)*3); // Length of square box
x=1:3:L; // creating x-axis
y=1:3:L; // creating y-axis
[gx,gy] = meshgrid(x,y); // places the atoms in the grid
i=1;
while i <= N,
x0(i) = gx(i) + randint(-1,1); // displacing atoms randomly in x-direction
y0(i) = gy(i) + randint(-1,1); // displacing atoms randomly in y-direction
i=i+1;
end
plot(gx(1:N),gy(1:N),'.r') //plot of atoms when placed initally at unifrom distance in red colour
xclick() // pauses the screen for next initiation for plotting another graph
plot(x0,y0,'.b') // plot of atoms when placed initally at unifrom distance in blue colour
i=1;
while i <= N,
theta = randint(0,2*%pi); // creates random value of angle for the ith atom
vx1(i)= V0*cos(theta); // generates random value of velocity in x direction
vy1(i)= V0*sin(theta);// generates random value of velocity in y direction
V(i)=sqrt(vx1(i)^2+vy1(i)^2);
i=i+1;
end
//Implementation of Velocity Verlet Algorithm
x1=x0; // initial value of x-axis
y1=y0; // initial value of y-axis
t=0; // Starting time of simulation
i=1;
while i<=N
ax(i)=0; // initial acceleration in x direction
ay(i)=0; //initial acceleration in x direction
j=1 // keeping record of interaction of other atoms
while j<=N
if i<>j // neglects interaction among the atoms themselves
dx=x1(i)-x1(j); // displacement between atoms i and j in y-direction
dy=y1(i)-y1(j); // displacement between atoms i and j in y-direction
r=sqrt(dx^2+dy^2); // distance between atoms i and j
if r<1 // condtition of repulsive forces
r=1;
end
if r<=3 // cut-off region for attarctive forces
ax(i)=ax(i)+24*(2/r^13-1/r^7)*(-dx/r); // calculation of acceleration in x-direction of ith atom
ay(i)=ay(i)+24*(2/r^13-1/r^7)*(-dy/r);// calculation of acceleration in x-direction of ith atom
end
end
j=j+1;
end
i=i+1;
end
// calculation of velocities and position further
f=1; // indexing ensemble
g=1; // indexing modulo
while t<=tf;
vx2=vx1+ax*dt/2; // calculating velocity at half time step in x direction
x1=x1+vx2*dt; // calculating positions further time in x-direction
vy2=vy1+ay*dt/2; // calculating velocity at half time step in y direction
y1=y1+vy2*dt; // calculating positions further time in y-direction
// Apllying the Periodic Boundary condition
k=1;
while k<=N //
if x1(k)<0 //chechking condition whether atom crosses the boundary along x-axis
x1(k)=x1(k)+L;// If condition satisfies changes the direction of atoms in x-axis
x0(k)=x0(k)+L;
end;
if x1(k)>L// chechking condition whether atom crosses the boundary along x-axis
x1(k)=x1(k)-L;// If condition satisfies changes the direction of atoms in x-axis
x0(k)=x0(k)-L;
end;
if y1(k)<0// chechking condition whether atom crosses the boundary along y-axis
y1(k)=y1(k)+L;// If condition satisfies changes the direction of atoms in y-axis
y0(k)=y0(k)+L;
end;
if y1(k)>L// chechking condition whether atom crosses the boundary along y-axis
y1(k)=y1(k)-L;// If condition satisfies changes the direction of atoms in y-axis
y0(k)=y0(k)-L;
end;
k=k+1
end
plot(x1,y1,'.g')
// calculation of new acceleration atfter dt time step
i=1;// index iteration used to calculate at later time
while i<=N
ax(i)=0;
ay(i)=0;
j=1
while j<=N
if i<>j
dx=x1(i)-x1(j);
dy=y1(i)-y1(j);
r=sqrt(dx^2+dy^2);
if r<1
r=1;
end
if r<=3
ax(i)=ax(i)+24*(2/r^13-1/r^7)*(-dx/r); // calculation of acceleration in x-direction of ith atom
ay(i)=ay(i)+24*(2/r^13-1/r^7)*(-dy/r);// calculation of acceleration in y-direction of ith atom
end
end
j=j+1;
end
i=i+1;
end
vx1=vx2+ax*dt/2; // velocities calculated at later time of step size dt in x-direction
vy1=vy2+ax*dt/2;// velocities calculated at later time of step size dt in y-direction
if(modulo(g,10)==0) // Generates ensemble copies after every 10 iteration
n=1;
while n<=N
vxe(f,n)=vx1(n); //redefining variable of velocities in x-direction
vye(f,n)=vy1(n);//redefining variable of velocities in y-direction
ve(f,n)=sqrt((vx1(n))^2+(vy1(n))^2);
n=n+1;
end
f=f+1;
end
g=g+1;
t=t+dt; //incrementation of time
end
// Calculation of Ensemble copies
vbin=0:.4:5;
va=length(vbin);
vt1=zeros(1,va-1);
k=1;
while k<=(f-1);
vt1=vt1+histc(ve(k,:),vbin,"countsNorm");
k=k+1;
end
scf(1); //Sets the current graphic figure
clf(); // clear and resets the figure
vbin=vbin+0.2;
plot(vbin(1:va-1),vt1/(f-1),'*',vbin(1:va-1),vt1/(f-1),'r');
vbin1=-3:.3:3;
vc1=length(vbin1);
vt2=zeros(1,vc1-1);
vt3=zeros(1,vc1-1);
k=1;
while k<=(f-1)
vt2=vt2+histc(vxe(k,:),vbin1,"countsNorm");
vt3=vt3+histc(vye(k,:),vbin1,"countsNorm");
k=k+1;
end
vbin1=vbin1+0.15;
scf(2);
clf();
subplot(1,2,1)
plot(vbin1(1:vc1-1),vt2/(f-1),'*',vbin1(1:vc1-1),vt2/(f-1),'r');
subplot(1,2,2);
plot(vbin1(1:vc1-1),vt3/(f-1),'*',vbin1(1:vc1-1),vt3/(f-1),'g');
endfunction
i is not incremented in the while i < n .... end instruction this explains why you got an infinite loop. I do not know what your algorithm is intended to do. May be you have to add something like i=i+1 just before the end of the while loop (see in the code below: "may be here")
function [x,y]=artrial(N,v0,tf)
dt=0.02;
t=0;
m=1;
L=sqrt(N)*4;
x=1:3:L
y=1:3:L
gx=zeros(1,N);
gy=zeros(1,N);
[gx,gy]=meshgrid(x,y)
i=1;
while i<=N
x0(i)=gx(i) + randint(-1,1);
y0(i)=gy(i) + randint(-1,1);
theta = randint(0,2*%pi)
vx0(i)=v0*cos(theta);
vy0(i)=v0*sin(theta);
v(i)=sqrt(vx0(i)^2+vy0(i)^2)
i=i+1;
end
plot(gx(1:N),gy(1:N),'*b')
xclick()
plot(x0,y0,'.r')
//force calculation
while t<tf
i=1
while i<N
j=1
fx(i)=0
fy(i)=0
if i<>j
while j<=N
dx=x0(j)-x0(i)
dy=y0(j)-y0(i)
if abs(dx)>L/2
dx=L-abs(dx)
end
if abs(dy)>=L/2
dy=L-abs(dy)
end
r=sqrt(dx^2+dy^2)
if r<1
r=1
end
if r<=3
f=24*((2./r^.13)-(1./r.^7))
fx(i)=fx(i)+f*dx./r
fy(i)=fy(i)+f*dy./r
end
j=j+1
end
i=i+1;// May be here,
end
ax(i)=fx(i)./m
ay(i)=fy(i)./m
end
for i=1:N
vhx(i)=vx0(i)+ax(i)*dt/2// velocity at half time step
vhy(i)=vy0(i)+ay(i)*dt/2
x1(i)=x0(i)+vhx(i)*dt
y1(i)=y0(i)+vhy(i)*dt
end
//boundary conditions
if x1(i)<0
x1(i)=x1(i)+L
x0(i)=x0(i)+L
end
if x1(i)>L
x1(i)=x1(i)-L
x0(i)=x0(i)-L
end
if y1(i)<0
y1(i)=y1(i)+L
y0(i)=y0(i)+L
end
if y1(i)>L
y1(i)=y1(i)-L
y0(i)=y0(i)-L
end
vx(i)=vhx(i)+ax(i)*dt/2
vy(i)=vhy(i)+ay(i)*dt/2
v(i)=sqrt(vx(i)^2+vy(i)^2)
x(i)=x0(i)+vx(i)*dt
y(i)=y0(i)+vy(i)*dt
i=i+1
plot(x,y,'.g')
//plot(Vx,Vy,'.y')
t=t+1
end
endfunction