[010]混合ガスの分子運動のプログラム

 投稿者：mike  投稿日：2017年 3月18日(土)13時17分7秒
  混合ガスの分子動力学法(MD2D)のプログラム(010mixLJMD2D.bas)を公開します。
　分子動力学法は、それぞれの分子に掛かる力を現在の分子の位置と分子間のポテンシャルから計算し、
現在の分子の位置と速度と力からNewtonの運動方程式を使って短い時間dt後の位置や速度を求め、
これを繰り返すことで運動を追跡する手法です。
　異種分子に働く力はLennard-Jonesポテンシャルの場合、epsilon（ポテンシャルの深さ）は
２つのポテンシャルの幾何平均、sigma（分子の大きさ）は代数平均にすると現実に近いことが知られています。
平衡状態にある時、分子の質量が小さい方が平均速度は速くなります。（速度空間を図示すると拡がりが大きくなります）

表示の説明：
初期画面では、青色の円がアルゴン分子、緑色の円がキセノン分子を表します。
キーを押すと反応します。[.]で終了、[0]は温度を変えます。出てきた調節バーで温度を変え、OKをクリックします。
[1]は画面表示を変えます。[9]はいくつかのテーマ（メニューに混合ガスの種類と個数が書かれています）の中から選択します。

試験環境：
本プログラムは十進BASIC 6.6.2 / macOS 10.7.5, 十進BASIC Ver 7.8.0 / windows 10でテストしました。


-------------------


!
! ========= molecular dynamics 2D ==========
!
! 010mixLJMD2D.bas
! Mitsuru Ikeuchi (C) Copyleft
!
! ver 0.0.1   2017.03.18  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB md2d.setInitialCondition, md2d.moveParticles, md2d.drawParticles
DECLARE EXTERNAL FUNCTION INKEY$
LET tempMode = 0    !tempMode: 0:adiabatic  1:constant temperature
LET contTemp = 150  !contTemp: controled constant temperature(K)
LET drawMode = 0    !drawMode: 0:ball 1:ball+v+f 2:velocitySpace
DIM menu$(8)
LET menu$(1) = "Ar50 Xe50"
LET menu$(2) = "Ar200 Xe100"
LET menu$(3) = "Ar50 Hg50"
LET menu$(4) = "Ar50 Kr50"
LET menu$(5) = "Kr50 Xe50"
LET menu$(6) = "Kr50 Hg50"
let menu$(7) = "Xe50 Hg50"
LET menu$(8) = "continue"

!setInitialCondition(moleculeKind,nMolecule,xMaximum,yMaximum,contTemp)
CALL setInitialCondition(2,50,4,50,8,8,contTemp) !moleculeKind: 2:Ar 3:Kr 4:Xe 5:Hg

DO
   CALL moveParticles(tempMode,contTemp)
   CALL drawParticles(tempMode,contTemp,drawMode,menu)
   LET S$=INKEY$
   IF S$="." THEN
      EXIT DO
   ELSEIF S$="0" THEN
      LOCATE VALUE (1),RANGE 0 TO 400 : temp
      IF temp>1 THEN
         LET tempMode = 1
         LET contTemp = temp
      ELSE
         LET tempMode = 0
      END IF
   ELSEIF S$="1" THEN
      LET drawMode = MOD(drawMode+1,3)
   ELSEIF S$="9" THEN
      LOCATE CHOICE (menu$) :nmenu
      IF nmenu=1 THEN
         CALL setInitialCondition(2,50,4,50,8,8,contTemp)
      ELSEIF nmenu=2 THEN
         CALL setInitialCondition(2,200,4,100,12,12,contTemp)
      ELSEIF nmenu=3 THEN
         LET tempMode = 1
         LET contTemp = 300
         CALL setInitialCondition(2,50,5,50,8,8,contTemp)
      ELSEIF nmenu=4 THEN
         CALL setInitialCondition(2,50,3,50,8,8,contTemp)
      ELSEIF nmenu=5 THEN
         CALL setInitialCondition(3,50,4,50,8,8,contTemp)
      ELSEIF nmenu=6 THEN
         CALL setInitialCondition(3,50,5,50,8,8,contTemp)
      ELSEIF nmenu=7 THEN
         CALL setInitialCondition(4,50,5,50,8,8,contTemp)
      ELSEIF nmenu=8 THEN  !contine
      !
      END IF

   END IF
LOOP
END

EXTERNAL FUNCTION INKEY$
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION

! ---------- Lennard-Jones md2d module ----------
!
!  method: velocity Verlet ( F=m*d^2r/dt^2 -> r(t+dt)=r(t)+v*dt,v=v+(F/m)*dt )
!    (1) vi = vi + (Fi/mi)*(0.5dt)
!    (2) ri = ri + vi*dt
!    (3) calculation Fi <- {r1,r2,...,rn} Fi=sum(Fij,j=1 to n),Fij=F(ri-rj)
!    (4) vi = vi + (Fi/mi)*(0.5dt)
!    (6) goto (1)
!  potential: Lennard-Jones V(r) = 4*epsilon*((sigma/r)^12-(sigma/r)^6)
!                     force F(r) = -dV(r)/dr
!                                = 24*epsilon*r6*(2*r6-1)/r, r6 = (sigma/r)^6
MODULE md2d
MODULE OPTION ARITHMETIC NATIVE
PUBLIC SUB setInitialCondition !(moleculeKind,nMolecule,xMaximum,yMaximum,contTemp)
PUBLIC SUB moveParticles !(tempMode,contTemp)
PUBLIC SUB drawParticles !(tempMode,contTemp,drawMode,menu)
SHARE NUMERIC sysTime,dt,nMolec,xMax,yMax, molecKind1,molecKind2, rCutoff,hh
SHARE NUMERIC xx(500),yy(500)   ! (xx(i),yy(i))  : position of i-th particle
SHARE NUMERIC vx(500),vy(500)   ! (vx(i),vy(i))  : velocity of i-th particle
SHARE NUMERIC ffx(500),ffy(500) ! (ffx(i),ffy(i)): total force of i-th particle
SHARE NUMERIC kind(500),mas(500)! kind(i),mas(i) : molec kind, mass of i-th particle
SHARE NUMERIC reg(500,0 TO 100) ! register near molec reg(i,0):number of near i-th molec
SHARE NUMERIC molecData(0 TO 5,0 TO 4) ! molecule 0:mass, 1:epsilon, 2:sigma, 3:dt  4:color
SHARE NUMERIC forceTable(0 TO 5, 0 TO 5,0 TO 1001)     ! force table
SHARE STRING molecStr$(0 TO 5)
LET sysTime = 0.0          ! system time (s) in the module
LET dt = 20.0*1.0e-15      ! time step (s)
LET nMolec = 36            ! total number of particles
LET xMax = 6.0E-9          ! x-Box size (m)
LET yMax = 6.0E-9          ! y-Box size (m)
LET molecKind1 = 2         ! molecule kind1 ( 2:Ar )
LET molecKind2 = 4         ! molecule kind2 ( 4:Xe )
LET rCutoff = 1e-9         ! force cutoff radius (m)
LET hh = 1e-12             ! force table step (m)
LET molecData(0,0)=0       ! if molecData(0,0)=0 then set molecData(,)

! ---------- set initial condition

EXTERNAL SUB setInitialCondition(kind1,nMol1,kind2,nMol2,xMaximum,yMaximum,contTemp)
   DECLARE EXTERNAL SUB setMoleculesData,setForceTable,setMolecules
   RANDOMIZE
   ! set particles
   CALL setMoleculesData
   LET sysTime = 0.0
   LET nMolec = nMol1+nMol2
   LET molecKind1 = kind1
   LET molecKind2 = kind2
   LET xMax = xMaximum*1e-9
   LET yMax = yMaximum*1e-9
   LET rCutoff = 1e-9
   CALL setForceTable
   CALL setMolecules(kind1,nMol1,kind2,nMol2,contTemp)
   ! set window
   SET WINDOW 0,500, 0,500
END SUB

EXTERNAL SUB setMoleculesData
!        0:mass(in AU) 1:eps(in kB) 2:sigma(m)      3:dt(s)   4:color
   DATA     4.003  ,     10.2 ,     2.576e-10 ,      5.0e-15,   1 !  0 He  black
   DATA     20.183 ,     36.2 ,     2.976e-10 ,     10.0e-15,  13 !  1 Ne  olieve
   DATA     39.948 ,    124.0 ,     3.418e-10 ,     20.0e-15,   2 !  2 Ar  blue
   DATA     83.50  ,    190.0 ,     3.610e-10 ,     20.0e-15,   3 !  3 Kr  green
   DATA    131.30  ,    229.0 ,     4.055e-10 ,     20.0e-15,  10 !  4 Xe  dark green
   DATA    200.59  ,    851.0 ,     2.898e-10 ,     20.0e-15,  12 !  5 Hg  brown

   !      0     1     2     3     4     5
   DATA  "He", "Ne", "Ar", "Kr", "Xe", "Hg"

   IF molecData(0,0)=0 THEN
      MAT READ molecData
      FOR i=0 TO 5
         LET molecData(i,0) = molecData(i,0)*1.661e-27 !mass(AU)--> (kg)
         LET molecData(i,1) = molecData(i,1)*1.38e-23  !eps(kB) --> (J)
      NEXT i
      MAT READ molecStr$
   END IF
END SUB

EXTERNAL SUB setMolecules(kind1,nMol1,kind2,nMol2,contTemp)
   DECLARE EXTERNAL SUB ajustVelocity
   LET sigmax = MAX(molecData(kind1,2),molecData(kind2,2))
   FOR j=1 TO nMol1+nMol2
      LET loopCount = 0
      DO
         LET xx(j) = (xMax-2*sigmax)*RND + sigmax
         LET yy(j) = (yMax-2*sigmax)*RND + sigmax
         FOR i=1 TO j-1
            IF (xx(i)-xx(j))^2+(yy(i)-yy(j))^2 < 2*sigmax^2 THEN EXIT FOR
         NEXT i
         LET loopCount = loopCount + 1
         IF loopCount>1000 THEN EXIT DO
      LOOP UNTIL i>=j
      IF loopCount>1000 THEN
         LET nMolec = j - 1
         EXIT FOR
      END IF
      IF j<=nMol1 THEN LET kind(j) = kind1 ELSE LET kind(j) = kind2
   NEXT j
   FOR i=1 TO nMolec
      LET vx(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
      LET vy(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
      LET ffx(i) = 0.0
      LET ffy(i) = 0.0
      LET mas(i) = molecData(kind(i),0)
   NEXT i
   CALL ajustVelocity(contTemp)
END SUB

! ---------- set force table

EXTERNAL SUB setForceTable
   DECLARE EXTERNAL  FUNCTION cutoff
   FOR ki=0 TO 5
      FOR kj=0 TO 5
         LET epsi = SQR(molecData(ki,1)*molecData(kj,1))
         LET sigm = 0.5*(molecData(ki,2)+molecData(kj,2))
         FOR ir=10 TO 1001
            LET r = ir*hh
            LET r6 = (sigm/r)^6
            LET forceTable(ki,kj,ir) = cutoff(r,rCutoff)*(24*epsi*r6*(2*r6-1)/r)
         NEXT ir
         FOR ir=0 TO 9
            LET forceTable(ki,kj,ir) = forceTable(ki,kj,10)
         NEXT ir
      NEXT kj
   NEXT ki
END SUB

EXTERNAL FUNCTION cutoff(r,rCutoff)
   IF r>0 AND r<0.8*rCutoff THEN
      LET ret = 1
   ELSEIF r>=0.8*rCutoff AND r<rCutoff THEN
      LET ret = 0.5+0.5*COS(PI*(r-0.8*rCutoff)/(0.2*rCutoff))
   else
      LET ret = 0
   END IF
   LET cutoff = ret
END FUNCTION

EXTERNAL FUNCTION force(r,ki,kj) !force(r) <-- forceTable - linear interporation
   LET ir = INT(r/hh)
   LET a = r - ir*hh
   LET force = ((hh-a)*forceTable(ki,kj,ir) + a*forceTable(ki,kj,ir+1))/hh
END FUNCTION

! ---------- move particles

EXTERNAL SUB moveParticles(tempMode,contTemp) !tempMode 0:adiabatic  1:constant-temp
   DECLARE EXTERNAL SUB registerNearMolec,moveParticlesDT,ajustVelocity
   IF (tempMode=1) THEN CALL ajustVelocity(contTemp)
   CALL registerNearMolec
   FOR i=1 TO 20
      CALL moveParticlesDT
   NEXT i
END SUB

EXTERNAL SUB moveParticlesDT ! velocity Verlet method
   DECLARE EXTERNAL SUB calcForce
   FOR i=1 TO nMolec
      LET a = 0.5*dt/mas(i)
      LET vx(i) = vx(i)+a*ffx(i)
      LET vy(i) = vy(i)+a*ffy(i)
      LET xx(i) = xx(i)+vx(i)*dt
      LET yy(i) = yy(i)+vy(i)*dt
   NEXT i
   CALL calcForce
   FOR i=1 TO nMolec
      LET a = 0.5*dt/mas(i)
      LET vx(i) = vx(i)+a*ffx(i)
      LET vy(i) = vy(i)+a*ffy(i)
   NEXT i
   LET sysTime=sysTime+dt
END SUB

EXTERNAL SUB calcForce
   DECLARE EXTERNAL FUNCTION force,boundaryForce
   LET s = 0.5*molecData(2,2) !molecData(2,2)=sigma for Ar
   LET rCut2 = rCutoff^2 ! force cutoff radius^2
   FOR i=1 TO nMolec
      LET ffx(i) = 0.0
      LET ffy(i) = 0.0
   NEXT i
   FOR i=1 TO nMolec-1
      FOR k=1 TO reg(i,0)-1
         LET j = reg(i,k)
         LET xij = xx(i)-xx(j)
         LET yij = yy(i)-yy(j)
         LET r2ij = xij*xij+yij*yij
         IF (r2ij<rCut2) THEN
            LET rij = SQR(r2ij)
            LET f = force(rij,kind(i),kind(j))
            LET fxij = f*xij/rij
            LET fyij = f*yij/rij
            LET ffx(i) = ffx(i)+fxij
            LET ffy(i) = ffy(i)+fyij
            LET ffx(j) = ffx(j)-fxij
            LET ffy(j) = ffy(j)-fyij
         END IF
      NEXT k
   NEXT i
   FOR i=1 TO nMolec  ! boundary force
      LET ffx(i) = ffx(i)+boundaryForce(xx(i)+s)+boundaryForce(xx(i)-xMax-s)
      LET ffy(i) = ffy(i)+boundaryForce(yy(i)+s)+boundaryForce(yy(i)-yMax-s)
   NEXT i
END SUB

!EXTERNAL FUNCTION force(r) ! force(r) = -dV(r)/dr
!   LET ri = sigma/r
!   LET r6 = ri^6
!   LET force = (24*epsilon*r6*(2*r6-1)/r)
!END FUNCTION

EXTERNAL FUNCTION boundaryForce(r)
   LET adsorp = 0.5*molecData(2,1) !molecData(2,1)=epsilon for Ar
   LET ri = molecData(2,2)/r !molecData(2,2)=sigma for Ar
   LET r6 = ri^6
   LET boundaryForce = (24.0*adsorp*r6*(2.0*r6-1.0)/r)
END FUNCTION

EXTERNAL SUB registerNearMolec
   LET rCut = rCutoff+20*2000*dt
   LET rcut2 = rCut*rCut
   FOR i=1 TO nMolec-1
      LET k = 1
      FOR j=i+1 TO nMolec
         LET r2 = (xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j))
         IF (r2<rcut2) THEN
            LET reg(i,k) = j
            LET k = k + 1
         END IF
      NEXT j
      LET reg(i,0) = k
   NEXT i
END SUB

EXTERNAL FUNCTION maxNearMolec
   LET mx = 0
   FOR i=1 TO nMolec-1
      IF mx<reg(i,0) THEN LET mx = reg(i,0)
   NEXT i
   LET maxNearMolec = mx-1
END FUNCTION

! ---------- utility

EXTERNAL FUNCTION systemTemprature
   LET kB = 1.38e-23 ! Boltzman's constant (J/K)
   LET ek= 0.0       !kinetic energy (J)
   FOR i=1 TO nMolec
      LET ek = ek + 0.5*mas(i)*(vx(i)*vx(i)+vy(i)*vy(i))
   NEXT i
   LET systemTemprature = ek/(nMolec*kB) !2D: E/N=kT, 3D: E/N=(3/2)kT
END FUNCTION

EXTERNAL SUB ajustVelocity(temp)
   DECLARE EXTERNAL FUNCTION systemTemprature
   LET r = sqr(temp/systemTemprature)
   FOR i=1 TO nMolec
      LET vx(i) = r*vx(i)
      LET vy(i) = r*vy(i)
   NEXT i
END SUB

! ---------- draw particles

EXTERNAL SUB drawParticles(tempMode,contTemp,drawMode,menu)
   DECLARE EXTERNAL FUNCTION systemTemprature,maxNearMolec
   DECLARE EXTERNAL PICTURE realSpace,velocitySpace
   SET DRAW MODE HIDDEN
   CLEAR
   IF drawMode=0 OR drawMode=1 THEN !--- 0:disk  1:circle+V+F
      DRAW realSpace(drawMode)
   ELSEIF drawMode=2 THEN
      DRAW velocitySpace
   END IF
   !--- control key guide
   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT 10,480 :"'.':exit  '0':Temp control (Temp<1 adiabatic mode)"
   PLOT TEXT, AT 10,460 :"'1':changeGraph  '9':select theme  "
   !--- draw caption
   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   LET tmp$ = "adiabatic   constantTemp  "
   PLOT TEXT, AT  50, 70 ,USING "time =#####.## (ps)   temp =####.## (K)":sysTime*1E12,systemTemprature
   PLOT TEXT, AT  50, 55 ,USING "N =####  max number of near molec =###":nMolec,maxNearMolec
   PLOT TEXT, AT  50, 40 ,USING "tempMode = ## ":tempMode
   PLOT TEXT, AT 200, 40 :tmp$(tempMode*12+1:tempMode*12+12)
   PLOT TEXT, AT  50, 25 ,USING "controled Temperature =####.# (K)":contTemp
   PLOT TEXT, AT  90, 10 :" in the box(2D molecular dynamics)"
   SET TEXT COLOR molecData(molecKind1,4) ! color of molec kind1
   PLOT TEXT, AT  50, 10 :molecStr$(molecKind1)
   SET TEXT COLOR molecData(molecKind2,4) ! color of molec kind2
   PLOT TEXT, AT  70, 10 :molecStr$(molecKind2)
   SET DRAW MODE EXPLICIT
END SUB

EXTERNAL PICTURE realSpace(drawMode)
   LET boxSize = 300
   LET deltat = 2e-14 !(s)
   LET sc = boxSize/xMax
   LET xp = 100
   LET yp = 100
   LET vScate = 100*deltat  !velocity line length = v*100*deltat
   LET fScale = 1000*deltat*deltat/molecData(2,0) !mass of Ar
   SET LINE COLOR 1 ! black : !--- box
   PLOT LINES: xp,yp; xp+boxSize,yp; xp+boxSize,yp+boxSize; xp,yp+boxSize; xp,yp
   SET TEXT HEIGHT 6
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT xp,yp+boxSize+2 ,USING  "box size =##.# x ##.# (nm)":xMax*1e9,yMax*1e9
   FOR i=1 TO nMolec
      IF drawMode=1 THEN  !--- draw circle, velocity and force
         SET LINE COLOR molecData(kind(i),4) ! molec color
         DRAW circle WITH SCALE(molecData(kind(i),2)/2*sc)*SHIFT(xp+xx(i)*sc,yp+yy(i)*sc)
         SET LINE COLOR 4 ! red : velocity
         PLOT LINES: xp+xx(i)*sc,yp+yy(i)*sc;
         PLOT LINES: xp+(xx(i)+vx(i)*vScate)*sc,yp+(yy(i)+vy(i)*vScate)*sc
         SET LINE COLOR 1 ! black : force
         PLOT LINES: xp+xx(i)*sc,yp+yy(i)*sc;
         PLOT LINES: xp+(xx(i)+ffx(i)*fScale)*sc,yp+(yy(i)+ffy(i)*fScale)*sc
      ELSE  !--- draw disk
         SET AREA COLOR molecData(kind(i),4) ! molec color
         DRAW disk WITH SCALE(molecData(kind(i),2)/2*sc)*SHIFT(xp+xx(i)*sc,yp+yy(i)*sc)
      END IF
   NEXT i
   IF drawMode=1 THEN
      LET xp = 100+boxSize*0.6
      LET yp = 100+boxSize+25
      SET LINE COLOR 4 ! red : velocity
      PLOT LINES: xp,yp;xp+23,yp
      SET TEXT COLOR 4 ! red
      PLOT TEXT, AT xp+30,yp-4: "velosity"
      SET LINE COLOR 1 ! black : force
      PLOT LINES: xp,yp-15;xp+23,yp-15
      SET TEXT COLOR 1 ! black
      PLOT TEXT, AT xp+30,yp-19: "force"
   END IF
END PICTURE

EXTERNAL PICTURE velocitySpace
   LET boxSize = 300
   LET xp = 100+boxSize/2
   LET yp = 100+boxSize/2
   SET LINE COLOR 1 !black : axis
   PLOT LINES: 100,yp; 100+boxSize,yp !vx-axis
   PLOT LINES: xp,100; xp,100+boxSize !vy-axis
   SET TEXT HEIGHT 6
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT 100+boxSize,yp: "vx"
   PLOT TEXT, AT 100+boxSize,yp-12: "1000m/s"
   PLOT TEXT, AT xp-12,100+boxSize: "vy  1000m/s"
   PLOT TEXT, AT xp-8,yp-10: "0"
   PLOT TEXT, AT 100,100+boxSize+8: "velocity space (vx,vy)"
   LET sc = boxSize/2000
   FOR i=1 TO nMolec
      SET LINE COLOR molecData(kind(i),4) ! molec color
      DRAW circle WITH SCALE(5)*SHIFT(vx(i)*sc+xp,vy(i)*sc+yp)
   NEXT i
END PICTURE

END MODULE
 

戻る