|
|
3D分子動力学法(MD)の雛形プログラム(006ArGasMD3D.bas)を公開します。
なるべくプログラムが読みやすいように、高速化や汎用性を狙わずに書いたつもりです。
分子動力学法は、それぞれの分子に掛かる力を現在の分子の位置と分子間のポテンシャルから計算し、
現在の分子の位置と速度と力からNewtonの運動方程式を使って短い時間dt後の位置や速度を求め、
これを繰り返すことで運動を追跡する手法です。
３次元の分子位置の２次元の画面上への表示はSECONDさんが掲載されたプログラムを熟読しながら、
理解できた範囲で書いてみました。
速度と力の矢線表示は見切れ(奥方向に向かう矢は見えないなど)の処理方法がわからず、実装できていません。
表示の説明：
緑色の円がアルゴン分子で、深い色ほど奥にあることを示します。
3Dでは分子が交差しても、衝突が起こるとは限りません。衝突は濃さの近い色の分子に限られます。
試験環境：
本プログラムは十進BASIC 6.6.0 / macOS 10.7.5, 十進BASIC Ver 7.7.8 / windows 10,
BASIC Acc 0.9.8.1 / windows 10でテストしました。
計算時間の比較： BASIC778 BASICAcc BASIC660
Ar 100個 Box 5x5x5nm 111.64s 33.05s 149.98s
Ar 300個 Box 6x6x6nm 895.07s 211.82s 850.09s
Ar 500個 Box 8x8x8nm 2428.13s 545.67s 2240.29s
BASIC778: BASIC 7.7.8 / windows 10 64bit / FMV-A561D CPU intel core-i5 2400M (2.5GHz) 4GB
BASICAcc: BASIC Acc 0.9.8.1 / windows 10 64bit / FMV-A561D CPU intel core-i5 2400M (2.5GHz) 4GB
BASIC660: BASIC 6.6.0 / MAC OS 10.7.5 / mac mini Intel Core-i7 (2.7GHz) 4GB
高速化の手法を用いていないため、分子数が増えると急激に計算時間がかかるようになります。
---------------------
!
! ========= molecular dynamics 3D ==========
!
! 006ArGasMD3D.bas
! Mitsuru Ikeuchi (C) Copyleft
!
! ver 0.0.1 2017.02.26 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 dAy = 0.5*(PI/180) ! Ay <- Ay+dAy, Ay:rotate angle around y-axis
!setInitialCondition(nMolecule,xMaximum,yMaximum,zMaximum,contTemp)
CALL setInitialCondition(100,5,5,5,200)
!CALL setInitialCondition(300,6,6,6,200)
!CALL setInitialCondition(500,8,8,8,200)
LET t0 = TIME
FOR it=1 TO 1000
LET Ay = Ay+dAy
FOR i=1 TO 10
CALL moveParticles
NEXT i
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, mass, sigma, epsilon
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 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 molecData(0 TO 12,0 TO 3) ! molecule 0:mass, 1:epsilon, 2:sigma, 3:dt
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
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 = 5.0E-9 ! x-Box size (m)
LET yMax = 5.0E-9 ! y-Box size (m)
LET zMax = 5.0e-9 ! z-Box size (m)
LET mass = 39.948*1.661e-27 ! mass of Ar (kg)
LET sigma = 3.418e-10 ! Lennard-Jones potential sigma for Ar (m)
LET epsilon =1.711e-21 ! Lennard-Jones potential epsilon FOR Ar (J)
! ---------- set initial condition
EXTERNAL SUB setInitialCondition(nMolecule,xMaximum,yMaximum,zMaximum,contTemp)
DECLARE EXTERNAL SUB setMorecules,setBox
RANDOMIZE
LET sysTime = 0.0
LET nMolec = nMolecule
LET xMax = xMaximum*1e-9
LET yMax = yMaximum*1e-9
LET zMax = zMaximum*1e-9
CALL setMorecules(nMolecule,contTemp)
CALL setBox
!set color mix
FOR i = 0 TO 100 !--- set color pallet
SET COLOR MIX(150+i) 0,0.01*i,0 !the darker green, the deeper z-position
NEXT i
! set window
SET WINDOW 0,500,0,500
END SUB
EXTERNAL SUB setMorecules(nMolecule,contTemp)
DECLARE EXTERNAL SUB ajustVelocity
FOR j=1 TO nMolecule
LET loopCount = 0
DO
LET xx(j) = (xMax-2*sigma)*RND + sigma
LET yy(j) = (yMax-2*sigma)*RND + sigma
LET zz(j) = (zMax-2*sigma)*RND + sigma
FOR i=1 TO j-1
IF (xx(i)-xx(j))^2+(yy(i)-yy(j))^2+(zz(i)-zz(j))^2 < 2*sigma^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
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
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
! ---------- move particles
EXTERNAL SUB moveParticles ! velocity Verlet method
DECLARE EXTERNAL SUB calcForce
LET a = 0.5*dt/mass
FOR i=1 TO nMolec
!LET a = 0.5*dt/mass
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/mass
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*sigma
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 j=i+1 TO nMolec
LET xij = xx(i)-xx(j)
LET yij = yy(i)-yy(j)
LET zij = zz(i)-zz(j)
LET rij = SQR(xij*xij+yij*yij+zij*zij)
LET f = force(rij)
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
NEXT j
NEXT i
FOR i=1 TO nMolec ! boundary force
LET ffx(i) = ffx(i)+boundaryForce(xx(i)+s)-boundaryForce(xMax+s-xx(i))
LET ffy(i) = ffy(i)+boundaryForce(yy(i)+s)-boundaryForce(yMax+s-yy(i))
LET ffz(i) = ffz(i)+boundaryForce(zz(i)+s)-boundaryForce(zMax+s-zz(i))
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 boundaryForce(r)
LET ri = sigma/r
LET r6 = ri^6
LET boundaryForce = (24.0*(0.5*epsilon)*r6*(2.0*r6-1.0)/r)
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*mass*(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
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)
SET AREA COLOR 150+int((ppz(j)/zMax+1)*40) ! the darker green, the deeper z-position
DRAW disk WITH SCALE(0.5*sigma*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,40 ,USING "time =#####.## (ps) temp =####.## (K)":sysTime*1E12,systemTemprature
PLOT TEXT, AT 50,25 ,USING "box size =##.# x ##.# x ##.# (nm) N=###":xMax*1e9,yMax*1e9,zMax*1e9,nMolec
PLOT TEXT, AT 50,10 :"Ar in the box (3D molecular dynamics)"
PLOT TEXT, AT 50,380 ,USING "Ax =####.#(deg) Ay =####.#(deg)":MOD(Ax*180/PI,360),MOD(Ay*180/PI,360)
SET DRAW MODE EXPLICIT
END SUB
END MODULE
|
|