|
Morseポテンシャルによる異種金属結晶の接触の分子動力学法プログラム017mixMMD2D.basを公開します。
異なった種類のナノ結晶を接触させます。接触すると表面エネルギーを解放し、温度が高くなります。
一般に表面エネルギー(Morseポテンシャルの場合は解離エネルギーに相当)の大きい方の表面に
小さい方の原子が覆う形が安定になります。
表示の説明:
鉄原子は暗い水色、アルミニウム原子は茶色の円であらわします。
表示は[D]のキーを押すことで変更できます。
試験環境:
本プログラムは十進BASIC 6.6.3.1 / macOS 10.7.5, 十進BASIC Ver 7.8.0 / windows 10 でテストしました。
---------------------
!
! ========= molecular dynamics 2D - Morse potential ==========
!
! 017mixMMD2D.bas
! Copyright(C) 2017 Mitsuru Ikeuchi
!
! ver 0.0.0 2017.04.22 created
! ver 0.0.1 2017.06.11 bug fixed
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB mmd2d.setInitialCondition, mmd2d.moveParticles, mmd2d.drawParticles
DECLARE EXTERNAL FUNCTION INKEY$
LET tempMode = 0 !tempMode: 0:adiabatic 1:constant temperature
LET contTemp = 300 !contTemp: controled constant temperature(K)
LET drawMode = 0 !drawMode: 0:bond 1:circle+bond 2:velocitySpace
DIM menu$(8)
LET menu$(1) = "Fe Al"
LET menu$(2) = "Fe Mo"
LET menu$(3) = "Fe Hg"
LET menu$(4) = "W Pb"
LET menu$(5) = "Ni Cr"
LET menu$(6) = "Ag Cu"
LET menu$(7) = "Al Ar"
LET menu$(8) = "continue"
!setInitialCondition(kind1,kind2,boxSizeInNM,contTemp)
!molKind: 0:W 1:Mo 2:Cr 3:Fe 4:Ni 5:Al 6:Pb 7:Cu 8:Ag 17:Ar 20:Hg
CALL setInitialCondition(3,5,6.0,contTemp)
DO
CALL moveParticles(tempMode,contTemp)
CALL drawParticles(tempMode,contTemp,drawMode)
LET S$=INKEY$
IF S$="0" THEN
LOCATE VALUE (1),RANGE 0 TO 4000 : temp
IF temp>1 THEN
LET tempMode = 1
LET contTemp = temp
ELSE
LET tempMode = 0
END IF
ELSEIF S$="D" OR S$="d" THEN
LET drawMode = MOD(drawMode+1,3)
ELSEIF S$="9" THEN
LOCATE CHOICE (menu$) :nmenu
IF nmenu=1 THEN !--- Fe Al box=6x6nm
CALL setInitialCondition(3,5,6.0,contTemp)
ELSEIF nmenu=2 THEN !--- Fe Mo
CALL setInitialCondition(3,1,6.0,contTemp)
ELSEIF nmenu=3 THEN !--- Fe Hg
CALL setInitialCondition(3,20,6.0,contTemp)
ELSEIF nmenu=4 THEN !--- W Pb
CALL setInitialCondition(0,6,6.0,contTemp)
ELSEIF nmenu=5 THEN !--- Ni Cr
CALL setInitialCondition(4,2,6.0,contTemp)
ELSEIF nmenu=6 THEN !--- Ag Cu
CALL setInitialCondition(8,7,6.0,contTemp)
ELSEIF nmenu=7 THEN !--- Al Ar
CALL setInitialCondition(5,17,6.0,contTemp)
ELSEIF nmenu=8 THEN !contine
END IF
END IF
LOOP UNTIL S$=CHR$(27)
END
EXTERNAL FUNCTION INKEY$ !from decimalBASIC library
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION
! ---------- molecular dynamics 2D - Morse potential ----------
!
! method: velocity Verlet algorithm
! (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: Morse V(r) = D*((1-EXP(-A*(r-r0)))^2-1)
! = D*(EXP(-2*A*(r-r0))-2*EXP(-A*(r-r0)))
! (D:dissociation energy, r0:bond length, A:width parameter { A=SQR(k/(2*D)) }
! force F(r) = -dV(r)/dr
! = 2*D*A*y*(y-1), y=EXP(-A*(r-r0))
MODULE mmd2d
MODULE OPTION ARITHMETIC NATIVE
PUBLIC SUB setInitialCondition !(molKind,boxSizeInNM,xtalSizeInNM,contTemp)
PUBLIC SUB moveParticles !(tempMode,contTemp)
PUBLIC SUB drawParticles !(drawMode:0:realSpace 1:velocitySpace)
SHARE NUMERIC molecKind1,molecKind2, sysTime, dt, nMolec, xMax, yMax, 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 kind(500),mas(500) ! kind(i),mas(i) : molec kind, mass 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 20,0 TO 3) ! molecule 0:mass, 1:epsilon, 2:sigma, 3:dt
SHARE STRING molecSTR$(0 TO 20) ! molec string
SHARE NUMERIC forceTable(0 TO 20, 0 TO 20,0 TO 1001) ! force table
LET molecKind1 = 3 ! molecule kind1 - 3:Fe
LET molecKind2 = 1 ! molecule kind2 - 1:Mo
LET sysTime = 0.0 ! system time (s) in the module
LET dt = 5.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 rCutoff = 1.0e-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(kind1,kind2,boxSizeInNM,contTemp)
DECLARE EXTERNAL SUB setMolecData,setMolecules,ajustVelocity,setForceTable
DECLARE EXTERNAL FUNCTION setCrystalBlock
RANDOMIZE
CALL setMolecData
LET molecKind1 = kind1
LET molecKind2 = kind2
LET sysTime = 0.0
LET xMax = boxSizeInNM*1.0e-9
LET yMax = boxSizeInNM*1.0e-9
! set particles
LET s = 0.1*xMax
LET xtalSize = 0.35*xMax
LET ss = 0.55*xMax
LET ii = setCrystalBlock(1, kind1, s, s, xtalSize, xtalSize, PI/4)
LET ii = setCrystalBlock(ii+1, kind2, s, ss, xtalSize, xtalSize, 0)
LET ii = setCrystalBlock(ii+1, kind2, ss, s, xtalSize, xtalSize, 0)
LET nMolec = setCrystalBlock(ii+1, kind1, ss, ss, xtalSize, xtalSize, 0)
CALL ajustVelocity(contTemp)
LET rCutoff = MIN(1.0e-9, 3.0*MAX(molecData(kind1,3),molecData(kind2,3)))
CALL setForceTable
SET COLOR MIX(240) 1.0,0.2,0.2 !for bond direction color
SET COLOR MIX(241) 0.6,0.6,0.0
SET COLOR MIX(242) 0.0,1.0,0.0
SET COLOR MIX(243) 0.0,0.6,0.6
SET COLOR MIX(244) 0.2,0.2,1.0
SET COLOR MIX(245) 0.8,0.0,0.8
! set window
SET WINDOW 0,500,0,500
END SUB
EXTERNAL SUB setMolecData
! Morse potential data
! 0:mass(AU) 1:D(eV) 2:A(m^-1) 3:r0(m)
DATA 183.85 , 0.9906, 1.4116e10, 3.032e-10 ! 0 W
DATA 95.94 , 0.8032, 1.5079e10, 2.976e-10 ! 1 Mo
DATA 51.996, 0.4414, 1.5721e10, 2.754e-10 ! 2 Cr
DATA 55.847, 0.4174, 1.3885e10, 2.845e-10 ! 3 Fe
DATA 58.71 , 0.4205, 1.4199e10, 2.780e-10 ! 4 Ni
DATA 26.98 , 0.2703, 1.1646e10, 3.253e-10 ! 5 Al
DATA 207.19 , 0.2348, 1.1836e10, 3.733e-10 ! 6 Pb
DATA 63.54 , 0.3429, 1.3588e10, 2.866e-10 ! 7 Cu
DATA 107.87 , 0.3323, 1.3690e10, 3.115e-10 ! 8 Ag
DATA 40.08 , 0.1623, 0.8054e10, 4.569e-10 ! 9 Ca
DATA 87.62 , 0.1513, 0.7878e10, 4.988e-10 ! 10 Sr
DATA 137.34 , 0.1416, 0.6570e10, 5.373e-10 ! 11 Ba
DATA 22.99 , 0.0633, 0.5900e10, 5.336e-10 ! 12 Na
DATA 39.102, 0.0542, 0.4977e10, 6.369e-10 ! 13 K
DATA 85.47 , 0.0464, 0.4298e10, 7.207e-10 ! 14 Rb
DATA 132.905, 0.0449, 0.4157e10, 7.557e-10 ! 15 Cs
DATA 20.183, 0.0031, 1.6500e10, 3.076e-10 ! 16 Ne
DATA 39.948, 0.0104, 1.3400e10, 3.816e-10 ! 17 Ar
DATA 83.80 , 0.0141, 1.2500e10, 4.097e-10 ! 18 Kr
DATA 131.30 , 0.0200, 1.2400e10, 4.467e-10 ! 19 Xe
DATA 200.59 , 0.0734, 1.4900e10, 3.255e-10 ! 20 Hg
DATA "W" ,"Mo","Cr","Fe","Ni","Al","Pb","Cu","Ag","Ca","Sr"
DATA "Ba","Na","K" ,"Rb","Cs","Ne","Ar","Kr","Xe","Hg"
IF molecData(0,0)=0 THEN
MAT READ molecData
FOR i=0 TO 20
LET molecData(i,0) = molecData(i,0)*1.661e-27 !mass(AU)--> (kg)
LET molecData(i,1) = molecData(i,1)*1.602e-19 !eps(eV) --> (J)
NEXT i
MAT READ molecSTR$
END IF
END SUB
EXTERNAL SUB setGas(nMolecule)
DECLARE EXTERNAL SUB ajustVelocity
LET r0 = molecData(knd,3)
FOR j=1 TO nMolecule
LET loopCount = 0
DO
LET xx(j) = (xMax-2*r0)*RND + r0
LET yy(j) = (yMax-2*r0)*RND + r0
FOR i=1 TO j-1
IF (xx(i)-xx(j))^2+(yy(i)-yy(j))^2 < 2*r0^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
END SUB
EXTERNAL FUNCTION setCrystalBlock(ii, knd, x0, y0, xLen, yLen, theta)
DECLARE EXTERNAL SUB setParticle
LET iip = ii
LET a = 0.98*molecData(knd,3) !r0 of knd
LET b = 0.866025*a
LET leng = xLen
IF (leng<yLen) THEN LET leng = yLen
LET leng = 1.5*leng
LET nx = INT(leng/b) + 1
LET ny = INT(leng/a) + 1
LET sth = SIN(theta)
LET cth = COS(theta)
FOR i=0 TO nx-1
LET x = b*i - leng/2.0
FOR j=0 TO ny-1
LET y = a*j - leng/2.0
IF MOD(i,2)=1 THEN LET y = y + 0.5*a
LET xp = x0 + xLen/2.0 + cth*x - sth*y
LET yp = y0 + yLen/2.0 + sth*x + cth*y
IF (xp>=x0 AND xp<=x0+xLen AND yp>=y0 AND yp<=y0+yLen) THEN
CALL setParticle(iip, knd, xp, yp)
LET iip = iip + 1
END IF
NEXT j
NEXT i
LET setCrystalBlock = iip - 1
end function
EXTERNAL SUB setParticle(i, knd, x, y)
LET xx(i) = x
LET yy(i) = y
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
LET mas(i) = molecData(knd,0)
LET kind(i) = knd
END SUB
EXTERNAL SUB setForceTable
DECLARE EXTERNAL FUNCTION cutoff
FOR ki=0 TO 20
FOR kj=0 TO 20
LET dd = SQR(molecData(ki,1)*molecData(kj,1))
LET aa = 0.5*(molecData(ki,2)+molecData(kj,2))
LET r0 = 0.5*(molecData(ki,3)+molecData(kj,3))
FOR ir=10 TO 1001
LET r = ir*hh
LET y = EXP(-aa*(r-r0))
LET forceTable(ki,kj,ir) = cutoff(r)*2.0*dd*aa*y*(y-1)
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 moveParticlesDT,ajustVelocity,registerNearMolec
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
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 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/mas(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*3.418e-10
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 rij = SQR(xij*xij+yij*yij)
IF rij<rCutoff THEN
LET f = force(rij,kind(i),kind(j))
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(xx(i)-xMax-s)
LET ffy(i) = ffy(i)+boundaryForce(yy(i)+s)+boundaryForce(yy(i)-yMax-s)
NEXT i
END SUB
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 = (3.418e-10/r)^6
LET boundaryForce = (24.0*(0.5*1.711e-21)*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 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)
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)
DECLARE EXTERNAL FUNCTION systemTemprature,maxNearMolec
DECLARE EXTERNAL sub realSpace,velocitySpace,plotBond
SET DRAW MODE HIDDEN
CLEAR
IF drawMode=0 OR drawMode=1 THEN !--- 0:circle 1:circle+bond
call plotBond(drawMode)
ELSEIF drawMode=2 THEN
call velocitySpace
END IF
!--- draw caption
SET TEXT HEIGHT 10
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 molecSTR$(molecKind1)&","&molecSTR$(molecKind2)&" N =####":nMolec
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"
PLOT TEXT, AT 50,480 :"[esc] exit [0] temp > value "
PLOT TEXT, AT 50,460 :"[D] dispMode [9] menu > select"
SET DRAW MODE EXPLICIT
END SUB
EXTERNAL sub plotBond(drawMode)
LET boxSize = 300
LET mag = boxSize/xMax
LET xp = 100
LET yp = 100
SET LINE COLOR 1 ! black : !--- box
PLOT LINES: xp,yp; boxSize+xp,yp; boxSize+xp,boxSize+yp; xp,boxSize+yp; xp,yp
SET TEXT HEIGHT 6
SET TEXT COLOR 1 ! black
PLOT TEXT, AT xp,boxSize+2+yp ,USING "box size =##.# x ##.# (nm)":xMax*1e9,yMax*1e9
FOR i=1 TO nMolec
IF kind(i)=molecKind1 THEN SET LINE COLOR 11 ELSE SET LINE COLOR 12
DRAW circle WITH SCALE(molecData(kind(i),3)/2.0*mag)*SHIFT(xx(i)*mag+xp,yy(i)*mag+yp)
NEXT i
IF drawMode=1 THEN
FOR i=1 TO nMolec-1
FOR k=1 TO reg(i,0)-1
LET j = reg(i,k)
LET r = SQR((xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j)))
LET r0 = 0.5*(molecData(kind(i),3)+molecData(kind(j),3))
IF r<1.2*r0 THEN
LET yij = yy(i)-yy(j)
IF yij=0.0 THEN LET yij = 1e-20
LET th = 3.0*(ATN((xx(i)-xx(j))/yij)+PI/2.0)/PI
LET col = 240 + INT((th-INT(th))*6)
SET LINE COLOR col
PLOT LINES: xx(i)*mag+xp,yy(i)*mag+yp;xx(j)*mag+xp,yy(j)*mag+yp
END IF
NEXT k
NEXT i
END IF
END sub
EXTERNAL sub velocitySpace
LET boxSize = 300
LET xp = 100
LET yp = 100
SET LINE COLOR 1 !black : axis
PLOT LINES: xp,boxSize/2+yp; boxSize+xp,boxSize/2+yp !vx-axis
PLOT LINES: boxSize/2+xp,yp; boxSize/2+xp,boxSize+yp !vy-axis
SET TEXT HEIGHT 6
SET TEXT COLOR 1 ! black
PLOT TEXT, AT boxSize+xp,boxSize/2+yp: "vx"
PLOT TEXT, AT boxSize+xp,boxSize/2-12+yp: "2000m/s"
PLOT TEXT, AT boxSize/2-12+xp,boxSize+yp: "vy 2000m/s"
PLOT TEXT, AT boxSize/2-8+xp,boxSize/2-10+yp: "0"
PLOT TEXT, AT xp,boxSize+8+yp: "velocity space (vx,vy)"
LET mag = boxSize/4000
FOR i=1 TO nMolec
IF kind(i)=molecKind1 THEN SET LINE COLOR 11 ELSE SET LINE COLOR 12
DRAW circle WITH SCALE(5)*SHIFT(vx(i)*mag+boxSize/2+xp,vy(i)*mag+boxSize/2+yp)
NEXT i
END sub
END MODULE
|
|