Sci-lab进行的经典分子动力学模拟

Question

我正在进行经典的分子动力学模拟，我们利用两个分子之间的相互作用这一事实，即在我们的情况下，稀有气体是由LJ电势曲线确定的。要解决这个问题，我们使用了速度velret方法。以下是在scilab中完成的代码和参数，其中N=32,v0=1,tf=15我的代码

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
 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

问题是我的scilab卡在了无限循环中，或者代码中出现了某些问题>任何人都会给它带来一些火花，这将是非常好的。

Answer 1

i在i