[012]混合ガスの分子運動(3次元)のプログラム

 投稿者：mike  投稿日：2017年 3月31日(金)11時58分51秒
  　3次元の混合ガスの分子動力学法(MD3D)のプログラム(012mixLJMD3D.bas)を公開します。
　分子動力学法は、それぞれの分子に掛かる力を現在の分子の位置と分子間のポテンシャルから計算し、
現在の分子の位置と速度と力からNewtonの運動方程式を使って短い時間dt後の位置や速度を求め、
これを繰り返すことで運動を追跡する手法です。
　異種分子に働く力はLennard-Jonesポテンシャルの場合、epsilon（ポテンシャルの深さ）は
２つのポテンシャルの幾何平均、sigma（分子の大きさ）は代数平均にすると現実に近いことが知られています。

表示の説明：
緑色の円がアルゴン分子、黄色の円がキセノン分子を表します。また、色の明度が低いほど奥にあることを表します。
箱の回転は立体感をだすためで、実際に箱が回転しているわけではありません。

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


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


!
! ========= molecular dynamics 3D ==========
!
! 0XXmixLJMD3D.bas
! Mitsuru Ikeuchi (C) Copyleft
!
! ver 0.0.1   2017.03.31  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB md3d.setInitialCondition, md3d.moveParticles, md3d.drawParticles
LET Ax = PI/12         ! rotate angle around x-axis
LET Ay = -PI/6         ! rotate angle around y-axis
LET ddAy = 0.5*(PI/180) ! Ay <- Ay+ddAy,  Ay:rotate angle around y-axis
LET contTemp = 200     ! controled temperature at initial and tempMode=1
LET tempMode = 0       ! tempMode 0:adiabatic  1:constant-temp

!setInitialCondition(kind1,nMol1,kind2,nMol2,xMaximum,yMaximum,zMaximum,contTemp)
!CALL setInitialCondition(2,150,4,150,6,6,6,contTemp) !molecKind=2 :Ar
CALL setInitialCondition(2,150,4,150,6,6,6,contTemp) !molecKind=4 :Xe
!CALL setInitialCondition(2,250,5,250,6,6,6,contTemp) !molecKind=5 :Hg

LET t0 = TIME
FOR it=1 TO 1000
   LET Ay = Ay+ddAy
   IF Ay>PI THEN LET Ay = Ay - 2*PI
   IF Ay<-PI THEN LET Ay = Ay + 2*PI
   CALL moveParticles(tempMode,contTemp)
   CALL drawParticles(Ax,Ay)
NEXT it
!PRINT time-t0
END

! ---------- Lennard-Jones md3d 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 md3d
MODULE OPTION ARITHMETIC NATIVE
PUBLIC SUB setInitialCondition,moveParticles,drawParticles
SHARE NUMERIC sysTime, dt, nMolec, xMax, yMax, zMax, molecKind1,molecKind2, rCutoff, hh
SHARE NUMERIC cosAx,sinAx,cosAy,sinAy, cx0,cy0,cz0 !--- 3D graphics
SHARE NUMERIC xx(500),yy(500),zz(500)    ! (xx(i),yy(i),zz(i))   : position of i-th particle
SHARE NUMERIC vx(500),vy(500),vz(500)    ! (vx(i),vy(i),vz(i))   : velocity of i-th particle
SHARE NUMERIC ffx(500),ffy(500),ffz(500) ! (ffx(i),ffy(i),ffz(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 ppx(500),ppy(500),ppz(500) ! (ppx(i),ppy(i),ppz(i)): rotated particle position
SHARE NUMERIC srtzix(500)                ! sort ppz(i) : sort z --> ppz(srtzix(i))
SHARE NUMERIC reg(500,0 TO 300)          ! 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 xApex(0 TO 7),yApex(0 TO 7),zApex(0 TO 7) ! boxApex x- y- z-coordinate
SHARE NUMERIC pxApex(0 TO 7),pyApex(0 TO 7),pzApex(0 TO 7) ! rotated boxApex x- y- z-coordinate
SHARE NUMERIC boxEdge(0 TO 11,0 TO 2)    ! boxEdge(i,j) i-th edge j=0,1:j-th apex, j=2:for marking
SHARE NUMERIC forceTable(0 TO 5, 0 TO 5,0 TO 1001)  ! force table
SHARE STRING molecStr$(0 TO 5) !molecStr$(2) = "Ar"
LET sysTime = 0.0          ! system time (s)
LET dt = 20.0*1.0e-15      ! time step (s)
LET nMolec = 100           ! number of particles
LET xMax = 6.0E-9          ! x-Box size (m)
LET yMax = 6.0E-9          ! y-Box size (m)
LET zMax = 6.0e-9          ! z-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 mat read molecData

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

EXTERNAL SUB setInitialCondition(kind1,nMol1,kind2,nMol2,xMaximum,yMaximum,zMaximum,contTemp)
   DECLARE EXTERNAL SUB setMoleculesData,setForceTable,setMolecules,setBox
   RANDOMIZE
   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 zMax = zMaximum*1e-9
   CALL setForceTable
   CALL setMolecules(kind1,nMol1,kind2,nMol2,contTemp)
   CALL setBox
   !set color mix
   FOR i = 0 TO 49 !--- set color pallet
      SET COLOR MIX(50+i) 0,0.02*i,0        !green, the darker green, the deeper z-position
      SET COLOR MIX(100+i) 0.02*i,0.02*i,0  !yellow the darker yellow, the deeper z-position
   NEXT i
   ! 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
         LET zz(j) = (zMax-2*sigmax)*RND + sigmax
         FOR i=1 TO j-1
            IF (xx(i)-xx(j))^2+(yy(i)-yy(j))^2+(zz(i)-zz(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 vz(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
      LET ffx(i) = 0.0
      LET ffy(i) = 0.0
      LET ffz(i) = 0.0
      LET mas(i) = molecData(kind(i),0)
   NEXT i
   CALL ajustVelocity(contTemp)
END SUB

EXTERNAL SUB setBox
   IF boxEdge(11,1)<>7 THEN
      DATA 0,0,0, 1,0,0, 0,1,0, 1,1,0, 0,0,1, 1,0,1, 0,1,1, 1,1,1
      FOR i=0 TO 7
         READ x,y,z
         LET xApex(i) = x*xMax
         LET yApex(i) = y*yMax
         LET zApex(i) = z*zMax
      NEXT i
      DATA 0,1,9, 0,2,9, 0,4,9, 1,3,9, 1,5,9, 2,3,9, 2,6,9, 3,7,9, 4,5,9, 4,6,9, 5,7,9, 6,7,9
      MAT READ boxEdge !(0 TO 11,0 TO 2)
   END IF
END SUB

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)*(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)
   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

! ---------- 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 10
      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 vz(i) = vz(i)+a*ffz(i)
      LET xx(i) = xx(i)+vx(i)*dt
      LET yy(i) = yy(i)+vy(i)*dt
      LET zz(i) = zz(i)+vz(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)
      LET vz(i) = vz(i)+a*ffz(i)
   NEXT i
   LET sysTime=sysTime+dt
END SUB

EXTERNAL SUB calcForce
   DECLARE EXTERNAL FUNCTION force,boundaryForce
   LET s = 0.5*molecData(2,2)
   LET rCut2 = rCutoff^2
   FOR i=1 TO nMolec
      LET ffx(i) = 0.0
      LET ffy(i) = 0.0
      LET ffz(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 zij = zz(i)-zz(j)
         LET r2ij = xij*xij+yij*yij+zij*zij
         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 fzij = f*zij/rij
            LET ffx(i) = ffx(i)+fxij
            LET ffy(i) = ffy(i)+fyij
            LET ffz(i) = ffz(i)+fzij
            LET ffx(j) = ffx(j)-fxij
            LET ffy(j) = ffy(j)-fyij
            LET ffz(j) = ffz(j)-fzij
         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)
      LET ffz(i) = ffz(i)+boundaryForce(zz(i)+s)+boundaryForce(zz(i)-zMax-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.0*epsilon*r6*(2.0*r6-1.0)/r)
!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

EXTERNAL FUNCTION boundaryForce(r)
   LET r6 = (molecData(2,2)/r)^6
   LET boundaryForce = (24.0*(0.5*molecData(2,1))*r6*(2.0*r6-1.0)/r)
END FUNCTION

EXTERNAL SUB registerNearMolec
   LET rCut = rCutoff+10*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))+(zz(i)-zz(j))*(zz(i)-zz(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)^2+vy(i)^2+vz(i)^2)
   NEXT i
   LET systemTemprature = 2*ek/(3*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)
      LET vz(i) = r*vz(i)
   NEXT i
END SUB

! ---------- 3D graphics aid

EXTERNAL SUB setRotateXY(angleX,angleY)
   LET cosAx = COS(angleX)
   LET sinAx = SIN(angleX)
   LET cosAy = COS(angleY)
   LET sinAy = SIN(angleY)
   LET cx0 = 0.5*xMax
   LET cy0 = 0.5*yMax
   LET cz0 = 0.5*zMax
END SUB

EXTERNAL SUB rotateXY !particles and box apex
   FOR i=1 TO nMolec
      LET ppx(i) = cosAy*(xx(i)-cx0)+sinAy*(sinAx*(yy(i)-cy0)+cosAx*(zz(i)-cz0)) + cx0
      LET ppy(i) = cosAx*(yy(i)-cy0)-sinAx*(zz(i)-cz0) + cy0
      LET ppz(i) =-sinAy*(xx(i)-cx0)+cosAy*(sinAx*(yy(i)-cy0)+cosAx*(zz(i)-cz0)) + cz0
   NEXT i
   FOR i=0 TO 7
      LET pxApex(i) = cosAy*(xApex(i)-cx0)+sinAy*(sinAx*(yApex(i)-cy0)+cosAx*(zApex(i)-cz0)) + cx0
      LET pyApex(i) = cosAx*(yApex(i)-cy0)-sinAx*(zApex(i)-cz0) + cy0
      LET pzApex(i) =-sinAy*(xApex(i)-cx0)+cosAy*(sinAx*(yApex(i)-cy0)+cosAx*(zApex(i)-cz0)) + cz0
   NEXT i
END SUB

EXTERNAL SUB markFarEdge
!  seek far apex --> iMin
   LET zMin = pzApex(0)
   LET iMin = 0
   FOR i=1 TO 7
      IF zMin>pzApex(i) THEN
         LET zMin = pzApex(i)
         LET iMin = i
      END IF
   NEXT i
   !mark far edge
   FOR iEdge = 0 TO 11
      LET boxEdge(iEdge,2) = 0
      IF (boxEdge(iEdge,0)=iMin OR boxEdge(iEdge,1)=iMin) THEN LET boxEdge(iEdge,2) = 1
   NEXT iEdge
END SUB

EXTERNAL SUB sortz
   DECLARE EXTERNAL SUB qSort
   FOR i=1 TO nMolec
      LET srtzix(i) = i
   NEXT i
   CALL qSort(ppz,srtzix,1,nMolec)
END SUB

EXTERNAL SUB qSort(m(),a(),le,ri) !modified decimal BASIC library
   IF ri>le THEN
      LET i=le-1
      LET j=ri
      LET pv=m(a(ri))
      DO
         DO
            LET i=i+1
         LOOP UNTIL pv<=m(a(i))
         DO
            LET j=j-1
         LOOP UNTIL j<=i OR m(a(j))<=pv
         IF j<=i THEN EXIT DO
         LET t=a(i)
         LET a(i)=a(j)
         LET a(j)=t
      LOOP
      LET t=a(i)
      LET a(i)=a(ri)
      LET a(ri)=t
      CALL qSort(m,a,le,i-1)
      CALL qSort(m,a,i+1,ri)
   END IF
END SUB

EXTERNAL SUB plotNearEdge(mag,xp,yp,lineColor)
   DECLARE EXTERNAL SUB plotEdge
   FOR iEdge=0 TO 11
      IF boxEdge(iEdge,2)=0 THEN !far edge mark = 0
         CALL plotEdge(iEdge,mag,xp,yp,lineColor)
      END IF
   NEXT iEdge
END SUB

EXTERNAL SUB plotFarEdge(mag,xp,yp,lineColor)
   DECLARE EXTERNAL SUB plotEdge
   FOR iEdge=0 TO 11
      IF boxEdge(iEdge,2)=1 THEN !far edge mark = 1
         CALL plotEdge(iEdge,mag,xp,yp,lineColor)
      END IF
   NEXT iEdge
END SUB

EXTERNAL SUB plotEdge(iEdge,mag,xp,yp,lineColor)
   SET LINE COLOR lineColor
   FOR i=0 TO 1
      LET iApex = boxEdge(iEdge,i)
      PLOT LINES: pxApex(iApex)*mag+xp, pyApex(iApex)*mag+yp;
   NEXT i
   PLOT LINES
END SUB

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

EXTERNAL SUB drawParticles(Ax,Ay) !Ax:rotate angle around s-axis, Ay:rotate angle around y-axis
   DECLARE EXTERNAL FUNCTION systemTemprature,maxNearMolec
   DECLARE EXTERNAL SUB setRotateXY,rotateXY,sortz,markFarEdge,plotFarEdge,plotNearEdge
   LET sc = 200/xMax   ! 200=boxsize in graphic window
   LET xp = 150
   LET yp = 120
   CALL setRotateXY(Ax,Ay)
   CALL rotateXY    !xx(i),yy(i),zz(i) rotate--> ppx(i),ppy(i),ppz(i)
   CALL sortz       !sort ppz(i) : ppz(srtzix(1))<ppz(srtzix(2))<...<ppz(srtzix(nMolec))
   CALL markFarEdge ! boxEdge(iEdge,2)=1:far side edge or 0:near side edge

   SET DRAW MODE HIDDEN
   CLEAR
   CALL plotFarEdge(sc,xp,yp,8) !8:gray
   FOR i=1 TO nMolec
      LET j = srtzix(i)
      IF kind(j)=molecKind1 THEN
         SET AREA COLOR 50+INT((ppz(j)/zMax+1)*20) ! the darker green, the deeper z-position
      ELSEIF kind(j)=molecKind2 THEN
         SET AREA COLOR 100+INT((ppz(j)/zMax+1)*20) ! the darker yellow, the deeper z-position
      END IF
      DRAW disk WITH SCALE(0.5*molecData(kind(j),2)*sc)*SHIFT(ppx(j)*sc+xp,ppy(j)*sc+yp)
   NEXT i
   CALL plotNearEdge(sc,xp,yp,1) !1:black
   !draw caption
   SET TEXT HEIGHT 8
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT 50,55 ,USING "time =#####.## (ps) temp =####.## (K)":sysTime*1E12,systemTemprature
   PLOT TEXT, AT 50,40 ,USING  "box size =##.# x ##.# x ##.# (nm)":xMax*1e9,yMax*1e9,zMax*1e9
   PLOT TEXT, AT 50,25 ,USING  "N=###    max number of registerd near molec=###":nMolec,maxNearMolec
   PLOT TEXT, AT 50,10 :molecStr$(molecKind1)&molecStr$(molecKind2)&" in the box (3D molecular dynamics)"
   PLOT TEXT, AT 50,380 ,USING "Ax =####.#(deg) Ay =####.#(deg)":Ax*180/PI,Ay*180/PI
   SET DRAW MODE EXPLICIT
END SUB

END MODULE
 

戻る