|
|
> No.4268[元記事へ]
MD粒子登録法による高速化プログラム（004registerLJMD2D.bas）の変更
MDの粒子登録法では近くの分子からの力のみを使うことで計算量をO(N^2)からO(N)に減らすことができます。
その際に用いられる力のカットオフは、ver0.0.1ではcutoff(r)=1(r<rCut), 0(r>rCut)としていました。
しかし、力が不連続に変わると、粒子の感じるポテンシャルは浅く評価されてしまい、エネルギーの保存性が
悪くなっていました。今回、cutoff(r)関数を滑らかにすることによって、改善しています。
下記に004registerLJMD2D.basのver0.0.2を掲載します。
--------------------
!
! ========= molecular dynamics 2D ==========
!
! 004registerLJMD2D.bas
! Mitsuru Ikeuchi (C) Copyleft
!
! ver 0.0.1 2017.02.12 created
! ver 0.0.2 2017.02.24 change cutoff -> smooth cutoff
!
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 = 1 !drawMode: 0:ball 1:ball+v+f 2:velocitySpace
DIM menu$(8)
LET menu$(1) = "Ar 36 6x6"
LET menu$(2) = "Ar100 8x8"
LET menu$(3) = "Ar300 12x12"
LET menu$(4) = "Kr100 8x8"
let menu$(5) = "Xe100 8x8"
LET menu$(6) = "Hg100 8x8"
let menu$(7) = "Ar500 15x15"
LET menu$(8) = "continue"
!setInitialCondition(molecKind,nMolecule,xMaximum,yMaximum,contTemp)
CALL setInitialCondition(2,100,8,8,contTemp) !molecKind: 1:Ne 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 10 TO 300 , AT 150 : temp
LOCATE VALUE (1),RANGE 0 TO 300 : 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 !--- Ar N=36 box=6x6nm
CALL setInitialCondition(2,36,6,6,contTemp)
ELSEIF nmenu=2 THEN !--- Ar N=100 box=8x8nm
CALL setInitialCondition(2,100,8,8,contTemp)
ELSEIF nmenu=3 THEN !--- Ar N=300 box=12x12nm
CALL setInitialCondition(2,300,12,12,contTemp)
ELSEIF nmenu=4 THEN !--- Kr N=100 box=8x8nm
CALL setInitialCondition(3,100,8,8,contTemp)
ELSEIF nmenu=5 THEN !--- Xe N=100 box=8x8nm
CALL setInitialCondition(4,100,8,8,contTemp)
ELSEIF nmenu=6 THEN !--- Hg N=100 box=8x8nm
CALL setInitialCondition(5,100,8,8,contTemp)
ELSEIF nmenu=7 THEN !--- Ar N=500 box=15x15nm
CALL setInitialCondition(2,500,15,15,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 !(molecKind,nMolecule,xMaximum,yMaximum,contTemp)
PUBLIC SUB moveParticles !(tempMode,contTemp)
PUBLIC SUB drawParticles !(drawMode:0:realSpace 1:velocitySpace)
SHARE NUMERIC sysTime, dt, nMolec, xMax, yMax, sigma, mass, epsilon, 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 reg(500,0 TO 100) ! 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 forceTable(0 TO 1200) ! force table
LET sysTime = 0.0 ! system time (s) in the module
LET dt = 20.0*1.0e-15 ! time step (s)
LET nMolec = 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.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 set molecData(,)
! ---------- set initial condition
EXTERNAL SUB setInitialCondition(molecKind,nMolecule,xMaximum,yMaximum,contTemp)
DECLARE EXTERNAL SUB setMoreculesData,setForceTable,setMorecules
RANDOMIZE
! set particles
CALL setMoreculesData
LET sysTime = 0.0
LET nMolec = nMolecule
LET xMax = xMaximum*1e-9
LET yMax = yMaximum*1e-9
LET mass = molecData(molecKind,0)
LET epsilon = molecData(molecKind,1)
LET sigma = molecData(molecKind,2)
LET rCutoff = MIN(1e-9, 3*sigma)
CALL setForceTable
CALL setMorecules(nMolecule,contTemp)
! 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
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
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
FOR i=1 TO j-1
IF (xx(i)-xx(j))^2+(yy(i)-yy(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 ffx(i) = 0.0
LET ffy(i) = 0.0
NEXT i
CALL ajustVelocity(contTemp)
END SUB
! ---------- set force table
EXTERNAL SUB setForceTable
DECLARE EXTERNAL FUNCTION cutoff
FOR ir=10 TO 1200
LET r = ir*hh
LET r6 = (sigma/r)^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 - linear interporation
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 20
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(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/mass(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*sigma
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)
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(xMax+s-xx(i))
LET ffy(i) = ffy(i)+boundaryForce(yy(i)+s)-boundaryForce(yMax+s-yy(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 adsorp = 0.5*1.67e-21 ! epsilon(Ar)=1.67e-21
LET ri = sigma/r
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 systemTime
! LET systemTime = sysTime
!END FUNCTION
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)
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 "
PLOT TEXT, AT 10,445 :"------------------------------------------------------"
!--- draw caption
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 50, 10 :"2-dimensional molecular dynamics"
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/mass
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 2 ! blue : molecule
DRAW circle WITH SCALE(sigma/2.0*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 2 ! blue : molecule
DRAW disk WITH SCALE(sigma/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 2 ! blue
DRAW circle WITH SCALE(5)*SHIFT(vx(i)*sc+xp,vy(i)*sc+yp)
NEXT i
END PICTURE
END MODULE
|
|