|
|
3D分子動力学法(MD)の雛形プログラム(006ArGasMD3D.bas)の高速化プログラム(007fastLJMD3D.bas)
を公開します。下記に示すように計算時間はおよそO(N)になりました。
高速化は、遠い分子からの力を無視する力のカットオフと、それぞれの分子について近くの分子を登録し、
その登録表をもとに近くの分子についてのみ力の計算をすることによって行っています。
また、滑らかなカットオフ関数の計算が煩雑なため、あらかじめ距離rと力fの表を作り、距離rから
線型補間により力fを得るようにしました。
表示の説明：
[006]のプログラムと同様に、緑色の円がアルゴン分子で、深い色ほど奥にあることを示します。
箱とその中の分子は回転に伴って回っているように見えますが、実際に箱が回転しているわけではなく、
視点が回転しているとになります。箱の中の分子たちは回転の影響を受けません。
試験環境：
本プログラムは十進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 24.77s 15.31s 83.49s
Ar 300個 Box 6x6x6nm 102.83s 50.44s 149.98s
Ar 500個 Box 8x8x8nm 186.37s 89.14s 232.08s
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 ==========
!
! 007fastLJMD3D.bas
! Mitsuru Ikeuchi (C) Copyleft
!
! ver 0.0.1 2017.03.02 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(molecKind,nMolecule,xMaximum,yMaximum,zMaximum,contTemp)
!CALL setInitialCondition(2,100,5,5,5,200) !molecKind=2 :Ar
CALL setInitialCondition(2,300,6,6,6,200) !molecKind=2 :Ar
!CALL setInitialCondition(2,500,8,8,8,200) !molecKind=2 :Ar
LET t0 = TIME
FOR it=1 TO 1000
LET Ay = Ay+dAy
CALL moveParticles(0,200) !(tempMode,contTemp) : tempMode 0:adiabatic 1:constant-temp
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, molecKind, mass, sigma, epsilon, 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 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 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
SHARE NUMERIC forceTable(0 TO 1200) ! force table
SHARE STRING molecStr$(0 TO 12)
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 molecKind = 2 ! molecule kind ( 2:Ar )
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)
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(moleculeKind,nMolecule,xMaximum,yMaximum,zMaximum,contTemp)
DECLARE EXTERNAL SUB setMoreculesData,setForceTable,setMorecules,setBox
RANDOMIZE
CALL setMoreculesData
LET sysTime = 0.0
LET molecKind = moleculeKind
LET nMolec = nMolecule
LET xMax = xMaximum*1e-9
LET yMax = yMaximum*1e-9
LET zMax = zMaximum*1e-9
LET mass = molecData(molecKind,0)
LET sigma = molecData(molecKind,2)
LET epsilon = molecData(molecKind,1)
LET rCutoff = MIN(1e-9, 3*sigma)
CALL setForceTable
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 setMoreculesData
! 0:mass(in AU) 1:eps(in kB) 2:sigma(m) 3:dt(s)
DATA 4.003 , 10.2 , 2.576e-10 , 5.0e-15 ! 0 He
DATA 20.183 , 36.2 , 2.976e-10 , 10.0e-15 ! 1 Ne
DATA 39.948 , 124.0 , 3.418e-10 , 20.0e-15 ! 2 Ar
DATA 83.50 , 190.0 , 3.610e-10 , 20.0e-15 ! 3 Kr
DATA 131.30 , 229.0 , 4.055e-10 , 20.0e-15 ! 4 Xe
DATA 200.59 , 851.0 , 2.898e-10 , 20.0e-15 ! 5 Hg
DATA 2.016 , 33.3 , 2.968e-10 , 5.0e-15 ! 6 H2
DATA 28.013 , 91.5 , 3.681e-10 , 10.0e-15 ! 7 N2
DATA 31.999 , 113.0 , 3.433e-10 , 10.0e-15 ! 8 O2
DATA 18.015 , 809.1 , 2.641e-10 , 10.0e-15 ! 9 H2O
DATA 16.043 , 137.0 , 3.822e-10 , 10.0e-15 ! 10 CH4
DATA 44.010 , 190.0 , 3.996e-10 , 20.0e-15 ! 11 CO2
DATA 28.011 , 110.0 , 3.590e-10 , 10.0e-15 ! 12 CO
!
! 0 1 2 3 4 5 6 7 8 9 10 11 12
DATA "He", "Ne", "Ar", "Kr", "Xe", "Hg", "H2", "N2", "O2", "H2O", "CH4", "CO2", "CO"
IF molecData(0,0)=0 THEN
MAT READ molecData
FOR i=0 TO 12
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 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
! ---------- set force table
EXTERNAL SUB setForceTable
DECLARE EXTERNAL FUNCTION cutoff
FOR ir=10 TO 1200
LET r = ir*hh
LET ri = sigma/r
LET r6 = ri^6
LET forceTable(ir) = cutoff(r,rCutoff)*(24*epsilon*r6*(2*r6-1)/r)
NEXT ir
FOR ir=0 TO 9
LET forceTable(ir) = forceTable(10)
NEXT ir
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) !force(r) <-- forceTable
LET ir = INT(r/hh)
LET a = r - ir*hh
LET force = ((hh-a)*forceTable(ir) + a*forceTable(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 10
CALL moveParticlesDT
NEXT i
END SUB
EXTERNAL SUB moveParticlesDT ! 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
LET rCut2 = rCutoff^2 !(3*sigma)^2 ! force cutoff rij>3*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 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)
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(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
EXTERNAL SUB registerNearMolec
LET rCut = rCutoff+10*2000*dt !3*sigma+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*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,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
LET Ay=Ay+dAy
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,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$(molecKind)&" 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
|
|