经典分子动力学
Classical molecular dynmaics
谁能解释一下如何解释系综平均概念，因为我在下面的 Scilab 中有工作代码，它演示了粒子在 e= 彼此相遇时经历吸引和排斥（如 LJ 势），并且可以看到在我的代码中获得的速度分布与麦克斯韦玻尔兹曼曲线相匹配。
拿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 在 while i < n .... end 指令中没有递增这解释了为什么你有一个无限循环。
我不知道你的算法打算做什么。可能您必须在 while 循环结束之前添加类似 i=i+1 的内容（参见下面的代码："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
谁能解释一下如何解释系综平均概念，因为我在下面的 Scilab 中有工作代码，它演示了粒子在 e= 彼此相遇时经历吸引和排斥（如 LJ 势），并且可以看到在我的代码中获得的速度分布与麦克斯韦玻尔兹曼曲线相匹配。
拿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 在 while i < n .... end 指令中没有递增这解释了为什么你有一个无限循环。 我不知道你的算法打算做什么。可能您必须在 while 循环结束之前添加类似 i=i+1 的内容（参见下面的代码："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