|
初投稿です。
分子動力学法による箱の中のArの運動のシミュレーションをしています。
以下に添付したdecimal basicのプログラムを動かしてみると、
windows (ver 7.7.8)では比較的滑らかに粒子が動くのですが、
mac (ver 0.6.5.2)では、ぎくしゃくとした動きをします。
対策はあるのでしょうか?
ご教示いただければ幸いです。
---------------
!
! ========= molecular dynamics 2D ==========
!
! ArGasMD2D.bas
! Mitsuru Ikeuchi
!
! ver 0.0.0 2017.01.11
!
DECLARE EXTERNAL NUMERIC md2d.sysTime ! system time
DECLARE EXTERNAL SUB md2d.setInitialCondition
DECLARE EXTERNAL SUB md2d.moveParticles
DECLARE EXTERNAL SUB md2d.drawParticles
DECLARE EXTERNAL FUNCTION md2d.sysTemp
!
CALL setInitialCondition
CALL drawParticles(100)
!
FOR it=1 TO 1000
FOR i=1 TO 20
CALL moveParticles
NEXT i
CALL drawParticles(100)
!IF (MOD(it,10)=0) THEN
! PRINT USING "time =####.## (ps) temp =####.## (K)":sysTime*1.0E12,sysTemp
!END IF
NEXT it
!
END
!
! ---------- md2d module ----------
!
MODULE md2d
PUBLIC NUMERIC sysTime
PUBLIC SUB setInitialCondition,moveParticles,drawParticles
PUBLIC FUNCTION sysTemp
SHARE NUMERIC dt,Nmt, xMax, yMax, sigm, mass, epsi, mag
SHARE NUMERIC xx(100),yy(100),vx(100),vy(100),ffx(100),ffy(100)
! xx(i),yy(i) : position of i-th paryicle
! vx(i),vy(i) : velocity of i-th paryicle
! ffx(i),ffy(i): force of i-th paryicle
LET sysTime = 0.0 ! system time (s)
LET dt = 20.0*1.0e-15 ! time step (s)
LET Nmt=36 ! number of particles
LET xMax =6.0E-9 ! x-Box size (m)
LET yMax =6.0E-9 ! y-Box size (m)
LET mass = 39.95*1.67e-27 ! mass of Ar (kg)
LET sigm = 3.40e-10 ! L-J potential sigma for Ar (m)
LET epsi = 1.67e-21 ! L-J potential epsilon FOR Ar (J)
LET mag = 1.0E9 ! (view scale(m)) = mag x (particle scale(nm))
!
EXTERNAL SUB setInitialCondition
DECLARE EXTERNAL SUB ajustV
! set particles
RANDOMIZE
LET s = xMax/7.0
FOR i=1 TO Nmt
LET xx(i) = s*MOD((i-1),6)+1.0*s
LET yy(i) = s*int((i-1)/6)+1.0*s
LET vx(i) = 500.0*(RND-0.5)
LET vy(i) = 500.0*(RND-0.5)
LET ffx(i) = 0.0
LET ffy(i) = 0.0
NEXT i
CALL ajustV(150)
! set window
LET b=2
SET WINDOW -b,xMax*mag+b,-b,yMax*mag+b
!
END SUB
!
EXTERNAL SUB moveParticles
DECLARE EXTERNAL SUB calcForce
LET a = 0.5*dt/mass
FOR i=1 TO Nmt
!LET a = 0.5*dt/mass
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 Nmt
!LET a = 0.5*dt/mass
LET vx(i) = vx(i)+a*ffx(i)
LET vy(i) = vy(i)+a*ffy(i)
NEXT i
! set boundary
FOR i=1 TO Nmt
if (xx(i) < 0.0) then
LET xx(i) = 0.0 ! -xx(i)
LET vx(i) = -vx(i)
END IF
IF (xx(i) > xMax) THEN
LET xx(i) = xMax ! 2*xMax-xx(i)
LET vx(i) = -vx(i)
END IF
if (yy(i) < 0.0) then
LET yy(i) = 0.0 ! -yy(i)
LET vy(i) = -vy(i)
END IF
if (yy(i) > yMax) then
LET yy(i) = yMax ! 2*yMax-yy(i)
LET vy(i) = -vy(i)
END IF
NEXT i
LET sysTime=sysTime+dt
END sub
!
EXTERNAL SUB calcForce
DECLARE EXTERNAL FUNCTION force
FOR i=1 TO Nmt
LET ffx(i) = 0.0
LET ffy(i) = 0.0
NEXT i
FOR i=1 TO Nmt-1
FOR j=i+1 TO Nmt
LET xij = xx(i)-xx(j)
LET yij = yy(i)-yy(j)
LET rij = SQR(xij*xij+yij*yij)
LET f = force(rij)
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
NEXT j
NEXT i
END sub
!
EXTERNAL function force(r) ! force(r) = -dV(r)/dr
LET ri = sigm/r
LET r6 = ri^6
LET force = (24.0*epsi*r6*(2.0*r6-1.0)/r)
END function
!
! utility
!
EXTERNAL FUNCTION sysTemp
LET k = 1.38e-23 ! Boltzman's constant (J/K)
LET totalEnergy = 0.0
FOR i=1 TO Nmt
LET totalEnergy = totalEnergy + 0.5*mass*(vx(i)^2+vy(i)^2)
NEXT i
LET sysTemp = totalEnergy/(Nmt*k)
END FUNCTION
!
EXTERNAL SUB ajustV(temp)
DECLARE EXTERNAL FUNCTION sysTemp
LET r = sqr(temp/sysTemp)
FOR i=1 TO Nmt
LET vx(i) = r*vx(i)
LET vy(i) = r*vy(i)
NEXT i
END sub
!
! display
!
EXTERNAL SUB drawParticles(vmag)
DECLARE EXTERNAL FUNCTION sysTemp
LET t1 = vmag*dt
LET b = (sigm/2)*mag
!
SET DRAW MODE HIDDEN
CLEAR
SET LINE COLOR 1 ! black : wall
PLOT LINES: -b,-b; xMax*mag+b,-b; xMax*mag+b,yMax*mag+b; -b,yMax*mag+b; -b,-b
!
FOR i=1 TO Nmt
SET COLOR 2 ! blue : position
DRAW circle WITH SCALE(sigm/2.0*1E9)*SHIFT(xx(i)*mag,yy(i)*mag)
SET LINE COLOR 4 ! red : velocity
PLOT LINES: xx(i)*mag,yy(i)*mag;
PLOT LINES: (xx(i)+vx(i)*t1)*mag,(yy(i)+vy(i)*t1)*mag
SET LINE COLOR 1 ! black : force
PLOT LINES: xx(i)*mag,yy(i)*mag;
PLOT LINES: (xx(i)+0.1*ffx(i)*t1*t1/mass)*mag,(yy(i)+0.1*ffy(i)*t1*t1/mass)*mag
NEXT i
!
SET TEXT COLOR 1 ! black
PLOT TEXT, AT 0,-0.7,USING "box size =##.# x ##.# (nm)":xMax*mag,yMax*mag
PLOT TEXT, AT 0,-1 ,USING "time =####.## (ps) temp =####.## (K)":sysTime*1.0E12,sysTemp
PLOT TEXT, AT 0,-1.3 :"Ar in the box (2-dimensional molecular dynamics)"
!
SET LINE COLOR 4 ! red : velocity
LET xPos = xMax*mag*0.5+1
LET yPos = yMax*mag+1
PLOT LINES: xPos,yPos;xPos+0.5,yPos
SET TEXT COLOR 4 ! red
PLOT TEXT, AT xPos+0.7,yPos-0.05: "velosity"
SET LINE COLOR 1 ! black : velocity
PLOT LINES: xPos,yPos-0.3;xPos+0.5,yPos-0.3
SET TEXT COLOR 1 ! black
PLOT TEXT, AT xPos+0.7,yPos-0.35: "force"
!
SET DRAW MODE EXPLICIT
END SUB
!
END MODULE
!
!
! end of file
|
|