Компьютерное моделирование термодинамических состояний

Определение параметров моделирования.
Число частиц:
>N:=60;
Длина грани кубической МД-ой ячейки:
> L:=50;
Шаг по времени:
> dt:=0.01;
Число шагов по времени межу усреднениями:
> Nt:=100;
Число наборов усреднения:
> Nn:=150;
 
ЗНАЧЕНИЕ НАЧАЛЬНЫХ СКОРОСТЕЙ И КООРДИНАТ
> r1:=rand(1000);
 
  r1 := proc()
    (proc() option builtin = RandNumberInterface;  end proc)(
    6, 1000, 10)
end proc
 
> for i to N do
> x[i]:=L*0.001*r1():y[i]:=0.1*L*0.001*r1():z[i]:=L*0.001*r1():
> end do:
> r2:=rand(-1000..1000);
 
  r2 := proc()
    (proc() option builtin = RandNumberInterface;  end proc)(
    6, 2001, 11) - 1000
end proc
 
> for i to N do
> Vx[i]:=3*0.001*r2():Vy[i]:=3*0.001*r2():Vz[i]:=3*0.001*r2():
> end do:
> Vxsum:=sum(Vx[k],k=1..N)/N;
 
                        Vxsum := 0.3553000000
 
> Vysum:=sum(Vy[k],k=1..N)/N;
 
                        Vysum := 0.2970000000
 
> Vzsum:=sum(Vz[k],k=1..N)/N;
 
                        Vzsum := -0.1954500000
 
> for i to N do Vx[i]:=Vx[i]-Vxsum:end do:
> for i to N do Vy[i]:=Vy[i]-Vysum:end do:
> for i to N do Vz[i]:=Vz[i]-Vzsum:end do:
> r:=sqrt(dx^2+dy^2+dz^2);
 
                              2     2     2 1/2
                      r := (dx  + dy  + dz )
 
> u:=4*(1/r^12-1/r^6);
 
                          4                    4
             u := ------------------ - ------------------
                     2     2     2 6      2     2     2 3
                  (dx  + dy  + dz )    (dx  + dy  + dz )
 
> fx:=-diff(u,dx):fy:=-diff(u,dy):fz:=-diff(u,dz);
 
                        48 dz                24 dz
            fz := ------------------ - ------------------
                     2     2     2 7      2     2     2 4
                  (dx  + dy  + dz )    (dx  + dy  + dz )
 
> Vmin:=0:Vmax:=10:nstep:=20:
> 
> step:=(Vmax-Vmin)/nstep:
> for j to nstep do particle[j]:=0:end do:
> a:=(0.136*10^7)/(165.324*3.405^3*6.02^2);
 
                           a := 5.749878567
 
> b:=(3.22*100)/(3.405^3*6.02);
 
                           b := 1.354902462
Нахождение ускорения частиц
> for i to N do ax[i]:=0:ay[i]:=0:az[i]:=0:end do:
> for i to N-1 do
> for j from (i+1) to N do
> dx:=x[i]-x[j]:dy:=y[i]-y[j]:dz:=z[i]-z[j]:
> if abs(dx)>0.5*L then dx:=dx-sign(dx)*L end if:
> if abs(dy)>0.5*L then dy:=dy-sign(dy)*L end if:
> if abs(dz)>0.5*L then dz:=dz-sign(dz)*L end if:
> ax[i]:=ax[i]+fx:ay[i]:=ay[i]+fy:az[i]:=az[i]+fz:
> ax[j]:=ax[j]-fx:ay[j]:=ay[j]-fy:az[j]:=az[j]-fz:
> end do;
> end do;
 
ПРОГРАММА МОЛЕКУЛЯРНОЙ ДИНАМИКИ
> t:=0:
> for k to Nn do
> virial:=0:n:=0:
> px:=0:py:=0:pz:=0:
> ke:=0:pe:=0:
> for it to Nt do
> for i to N do
> x[i]:=x[i]+Vx[i]*dt+0.5*ax[i]*dt^2:
> y[i]:=y[i]+Vy[i]*dt+0.5*ay[i]*dt^2:
> z[i]:=z[i]+Vz[i]*dt+0.5*az[i]*dt^2:
> Vx[i]:=Vx[i]+0.5*ax[i]*dt:
> Vy[i]:=Vy[i]+0.5*ay[i]*dt:
> Vz[i]:=Vz[i]+0.5*az[i]*dt:
> while x[i]<0 do x[i]:=x[i]+L: px:=px-Vx[i]:end do:
> while x[i]>L do x[i]:=x[i]-L: px:=px+Vx[i]:end do:
> while y[i]<0 do y[i]:=y[i]+L: py:=py-Vy[i]:end do:
> while y[i]>L do y[i]:=y[i]-L: py:=py+Vy[i]:end do:
> while z[i]<0 do z[i]:=z[i]+L: pz:=pz-Vz[i]:end do:
> while z[i]>L do z[i]:=z[i]-L: pz:=pz+Vz[i]:end do:
> end do:
> for i to N do
> if y[i]<=0.1*L then n:=n+1 end if:
> ax[i]:=0:ay[i]:=0:az[i]:=0:
> end do:
> for i to N-1 do
> for j from (i+1) to N do
> dx:=x[i]-x[j]:dy:=y[i]-y[j]:dz:=z[i]-z[j]:
> if abs (dx)>0.5*L then dx:=dx-sign(dx)*L end if:
> if abs (dy)>0.5*L then dy:=dy-sign(dy)*L end if:
> if abs (dz)>0.5*L then dz:=dz-sign(dz)*L end if:
> ax[i]:=ax[i]+fx:ay[i]:=ay[i]+fy:az[i]:=az[i]+fz:
> ax[i]:=ax[i]-fx:ay[i]:=ay[i]-fy:az[i]:=az[i]-fz:
> pe:=pe+u:
> end do:
> end do:
> for i to N do
> Vx[i]:=Vx[i]+0.5*ax[i]*dt:
> Vy[i]:=Vy[i]+0.5*ay[i]*dt:
> Vz[i]:=Vz[i]+0.5*az[i]*dt:
> ke:=ke+0.5*(Vx[i]^2+Vy[i]^2+Vz[i]^2):
> virial:=virial+x[i]*ax[i]+y[i]*ay[i]+z[i]*az[i]:
> if k>=Nn-1 then
> V:=sqrt(Vx[i]^2+Vy[i]^2+Vz[i]^2):
> for j to nstep do
> if (V>=(Vmin+step*(j-1)))and(V<(Vmin+step*j)) then particle[j]:=particle[j]+1 end if:
> end do:
> end if:
> end do:
> end do:
> t:=t+dt*Nt:
> conc[k]:=(n*10^5)/(Nt*0.1*(L*3.405)^3):
> Ke[k]:=(ke*165.324*10^(-4))/Nt:
> Pe[k]:=(pe*165.324*10^(-4))/Nt:
> E[k]:=(Pe[k]+Ke[k])/N:
> T[k]:=(2*ke*119.8)/(3*N*Nt):Pid:=N*1653.24*T[k]/(3.405*L)^3:
> PVdV:=((T[k]/(L^3/N-b))-((N^2*a)/L^6))*(1653.24/3.405^3):
> Pflux[k]:=((px+py+pz)*1653.24)/(3*L^2*3.405^3*dt*Nt):
> Pvirial[k]:=((N*T[k])/L^3+virial/(3*Nt*L^3))*(1653.24/3.405^3):
> dPflux1[k]:=Pflux[k]-Pid:dPvirial1[k]:=Pvirial[k]-Pid:
> dPflux2[k]:=Pflux[k]-PVdV:dPvirial2[k]:=Pvirial[k]-PVdV:
> end do:
Warning, computation interrupted
 
> for j to nstep do particle[j]:=particle[j]/(N*2*Nt):end do:
> EE:=[seq(E[i],i=1..Nn)]:TT:=[seq(T[i],i=1..Nn)]:
> PE:=[seq(Pe[i],i=1..Nn)]:KE:=[seq(Ke[i],i=1..Nn)]:
> PFLUX:=[seq(Pflux[i],i=1..Nn)]:
> PVIRIAL:=[seq(Pvirial[i],i=1..Nn)]:
> Conc:=[seq(conc[i],i=1..Nn)]:
> Particle:=[seq(particle[i],i=1..nstep)]:
> DPflux2:=[seq(dPflux2[i],i=1..Nn)]:
> DPvirial2:=[seq(dPvirial2[i],i=1..Nn)]:
> DPflux1:=[seq(dPflux1[i],i=1..Nn)]:
> DPvirial1:=[seq(dPvirial1[i],i=1..Nn)]:
> with(plots):
> listplot(PE);
> 
 
                             listplot(PE)
 
> listplot(KE);
 
                             listplot(KE)
 
> listplot(EE);
 
                             listplot(EE)
 
> listplot(PE);
 
                             listplot(PE)
 
> listplot(TT);
 
                             listplot(TT)
 
> listplot(PFLUX);
 
                           listplot(PFLUX)
 
> listplot(DPflux1);
 
                          listplot(DPflux1)
 
> listplot(DPflux2);
 
                          listplot(DPflux2)
 
> listplot(PVIRIAL);
 
                          listplot(PVIRIAL)
 
> listplot(DPvirial1);
 
                         listplot(DPvirial1)
 
> listplot(DPvirial2);
 
                         listplot(DPvirial2)
 
> listplot(Conc);
 
                            listplot(Conc)
 
> listplot(Particle);