[023]イオン結晶（２次元）投稿者：mike 投稿日：2017年 8月 4日(金)09時33分26秒 |
|
イオン結晶の分子動力学法プログラム023NaClIonMD2D.basを公開します。 NaClのような陽イオンと陰イオンからなる系において、電子殻の相互作用による引力や斥力に加え、静電的な 引力や斥力が働きます。静電気力は1/r^2に比例するため、比較的遠くまでその力が及びます。このため、分子動力学の 高速化（O(N)）に必要な近くの粒子だけで力の計算をすることが難しくなります。 しかしながら、多くの場合、陽イオンと陰イオンは近くに存在し、中性に近くなるので、互いに静電気力を打ち消す方向に 働くため、静電気力は急激に減少します。これをDebye遮蔽(1/r^2)*exp(-ar)といいます。 本プログラムではDebye遮蔽を前提として力の計算をしています。（粒子登録法によるO(N)の高速化を行っています） 表示の説明： 赤色は陽イオン、青色は陰イオンをあらわします。 [K]を押し、温度を上げていくと結晶は不安定になりやがて融解します。(融点は2次元のため現実とは異なります) [J]を押し、温度を下げていくと再び結晶に戻ります。 [1]-[7]はイオン結晶の種類を選択します。 表示は[D]のキーを押すことで変更できます。 試験環境： 本プログラムは十進BASIC 6.6.3.3 / macOS 10.7.5 でテストしました。 ---------- ! ! ========= molecular dynamics 2D - ion ========== ! ! 023NaClIonMD2D.bas ! Copyright(C) 2017 Mitsuru Ikeuchi ! ! ver 0.0.1 2017.08.03 created ! OPTION ARITHMETIC NATIVE DECLARE EXTERNAL SUB imd2d.setInitialCondition, imd2d.moveParticles, imd2d.drawParticles DECLARE EXTERNAL FUNCTION INKEY$ LET tempMode = 0 !tempMode: 0:adiabatic 1:constant temperature LET contTemp = 300 !contTemp: controled constant temperature(K) LET ddTemp = 0 !contTemp = contTemp + ddTemp LET drawMode = 1 !drawMode: 0:bond 1:circle+bond 2:velocitySpace LET pauseFlag = 0 !if pauseFlag=0, CALL moveParticles(tempMode,contTemp) !setInitialCondition(material,boxSizeInNM,contTemp) CALL setInitialCondition(1,6.0,contTemp) !material = 1:NaCl 2:MgO 3:CaO 4:BaO 5:NaF 6:KF 7:KCl DO LET contTemp = contTemp + ddTemp IF contTemp>3000 THEN LET contTemp = 3000 IF contTemp<10 THEN LET contTemP = 10 IF pauseFlag=0 THEN CALL moveParticles(tempMode,contTemp) ELSE WAIT DELAY 0.05 CALL drawParticles(tempMode,contTemp,drawMode) LET S$=INKEY$ IF S$="1" OR S$="2" OR S$="3" OR S$="4" OR S$="5" OR S$="6" OR S$="7" THEN LET material = VAL(S$) LET tempMode = 0 LET contTemp = 300 CALL setInitialCondition(material,6.0,contTemp) ELSEIF S$="T" OR S$="t" THEN LET tempMode = MOD(tempMode+1,2) IF tempMode=0 THEN LET ddTemp = 0 ELSEIF S$="K" OR S$="k" THEN LET tempMode = 1 IF ddtemp=0 THEN LET ddTemp = 1 ELSE LET ddTemp = 0 ELSEIF S$="J" OR S$="j" THEN LET tempMode = 1 IF ddTemp=0 THEN LET ddTemp = -1 ELSE LET ddTemp = 0 ELSEIF S$="D" OR S$="d" THEN LET drawMode = MOD(drawMode+1,3) ELSEIF S$=" " THEN LET pauseFlag = MOD(pauseFlag+1,2) END IF LOOP UNTIL S$=CHR$(27) END EXTERNAL FUNCTION INKEY$ !--- from decimal BASIC library inkey$.bas OPTION ARITHMETIC NATIVE SET ECHO "OFF" LET S$="" CHARACTER INPUT NOWAIT: S$ LET INKEY$=S$ END FUNCTION ! ---------- molecular dynamics 2D - ion ---------- ! ! 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) ! (5) goto (1) ! ! force: ion f(r) = fc + fr + fa ! fc = eForceConst*zi*zj*(EXP(-r/6.5e-10)/r)*(1.0/r+1.0/6.5e-10) !Debye-screened Coulomb force ! fr = 6.9742e-11*EXP((a-r)/b) !repulsive force ! fa = -6.9742e-21*(c/r^6) !attractive force ! MODULE imd2d 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 ionKind1,ionKind2, sysTime, dt, nMolec, xMax, yMax, rCutoff, hh SHARE NUMERIC xx(2000),yy(2000) ! (xx(i),yy(i)) : position of i-th particle SHARE NUMERIC vx(2000),vy(2000) ! (vx(i),vy(i)) : velocity of i-th particle SHARE NUMERIC ffx(2000),ffy(2000) ! (ffx(i),ffy(i)): total force of i-th particle SHARE NUMERIC kind(2000),mas(2000) ! kind(i),mas(i) : molec kind, mass of i-th particle SHARE NUMERIC reg(2000,0 TO 100) ! register near molec reg(i,0):number of near i-th molec SHARE NUMERIC ion_m(0 TO 17),ion_z(0 TO 17) ! ion mass,ion charge SHARE NUMERIC ion_a(0 TO 17),ion_b(0 TO 17) ! ion force param a,ion force param b SHARE NUMERIC ion_c(0 TO 17),ion_r(0 TO 17) ! ion force param c, ion radius SHARE NUMERIC ion_col(0 TO 17) ! ion draw color SHARE STRING ion_str$(0 TO 17) ! ion string SHARE NUMERIC forceTable(0 TO 17, 0 TO 17,0 TO 1001) ! force table LET ionKind1 = 3 ! ion kind1 - 3:Na+ LET ionKind2 = 7 ! ion kind2 - 7:Cl- LET sysTime = 0.0 ! system time (s) in the module LET dt = 2.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 ion_z(0)=0.0 ! if ion_z(0)=0.0 then set ion data ! ---------- set initial condition EXTERNAL SUB setInitialCondition(material,boxSizeInNM,contTemp) DECLARE EXTERNAL SUB setIonData,ajustVelocity,setForceTable DECLARE EXTERNAL FUNCTION setNaClTypeBlock RANDOMIZE LET sysTime = 0.0 CALL setIonData CALL setForceTable LET xMax = boxSizeInNM*1.0e-9 LET yMax = boxSizeInNM*1.0e-9 IF material=1 THEN !NaCl LET ionKind1 = 3 !Na+ LET ionKind2 = 7 !Cl- LET lattice = 5.6407e-10 ELSEIF material=2 THEN !MgO LET ionKind1 = 4 !Mg++ LET ionKind2 = 1 !O-- LET lattice = 4.212e-10*0.94 !0.94: correction factor ELSEIF material=3 THEN !CaO LET ionKind1 = 9 !Ca++ LET ionKind2 = 1 !O-- LET lattice = 4.80e-10*0.94 !0.94: correction factor ELSEIF material=4 THEN !BaO LET ionKind1 = 12 !Ba++ LET ionKind2 = 1 !O-- LET lattice = 5.536e-10 ELSEIF material=5 THEN !NaF LET ionKind1 = 3 !Na+ LET ionKind2 = 2 !F- LET lattice = 4.62e-10 ELSEIF material=6 THEN !KF LET ionKind1 = 8 !K+ LET ionKind2 = 2 !F- LET lattice = 5.34e-10 ELSEIF material=7 THEN !KCl LET ionKind1 = 8 !K+ LET ionKind2 = 7 !Cl- LET lattice = 6.29e-10 END IF LET s = 0.5*(xMax - lattice*6) ! setNaClTypeBlock(ii, knd1, knd2, nx, ny, lattice, xPos, yPos) LET nMolec = setNaClTypeBlock(1, ionKind1, ionKind2, 6, 6, lattice, s, s) CALL ajustVelocity(contTemp) ! set window SET WINDOW 0,500,0,500 END SUB EXTERNAL SUB setIonData ! ion potential data ! 0 mass,1 charge, 2 a , 3 b , 4 c , 5 r-ion, 6 color 7 str$ DATA 10.81, 3.0, 0.720e-10, 0.080e-10, 0.0, 0.23e-10, 13 , "B+++" !0 DATA 16.00, -2.0, 1.626e-10, 0.085e-10, 20.0, 1.40e-10, 2 , "O--" !1 DATA 19.00, -1.0, 1.565e-10, 0.085e-10, 20.0, 1.33e-10, 2 , "F-" !2 DATA 22.99, 1.0, 1.260e-10, 0.080e-10, 20.0, 1.02e-10, 4 , "Na+" !3 DATA 24.31, 2.0, 1.161e-10, 0.080e-10, 10.0, 0.72e-10, 4 , "Mg++" !4 DATA 26.98, 3.0, 1.064e-10, 0.080e-10, 2.0, 0.53e-10, 13 , "Al+++" !5 DATA 28.09, 4.0, 1.012e-10, 0.080e-10, 0.0, 0.40e-10, 7 , "Si++++" !6 DATA 35.45, -1.0, 1.950e-10, 0.090e-10, 30.0, 1.81e-10, 2 , "Cl-" !7 DATA 39.10, 1.0, 1.595e-10, 0.080e-10, 15.0, 1.38e-10, 4 , "K+" !8 DATA 40.08, 2.0, 1.414e-10, 0.080e-10, 10.0, 1.00e-10, 4 , "Ca++" !9 DATA 47.88, 4.0, 1.235e-10, 0.080e-10, 0.0, 0.61e-10, 7 , "Ti++++" !10 DATA 87.62, 2.0, 1.632e-10, 0.080e-10, 15.0, 1.16e-10, 4 , "Sr++" !11 DATA 137.3, 2.0, 1.820e-10, 0.080e-10, 20.0, 1.36e-10, 4 , "Ba++" !12 DATA 4.003, 0.0, 1.200e-10, 0.110e-10, 4.76, 1.28e-10, 3 , "He" !13 DATA 20.18, 0.0, 1.415e-10, 0.112e-10,11.03, 1.37e-10, 3 , "Ne" !14 DATA 39.95, 0.0, 1.878e-10, 0.117e-10,38.53, 1.70e-10, 3 , "Ar" !15 DATA 83.80, 0.0, 2.041e-10, 0.130e-10,55.33, 1.83e-10, 3 , "Kr" !16 DATA 131.3, 0.0, 2.258e-10, 0.145e-10,85.55, 1.99e-10, 3 , "Xe" !17 IF ion_z(0)=0.0 THEN FOR i=0 TO 17 READ ion_m(i),ion_z(i),ion_a(i),ion_b(i),ion_c(i),ion_r(i),ion_col(i),ion_str$(i) LET ion_m(i) = ion_m(i)*1.661e-27 LET ion_z(i) = ion_z(i)*1.602e-19 NEXT i END IF END SUB EXTERNAL FUNCTION setNaClTypeBlock(ii, knd1, knd2, nx, ny, lattice, xPos, yPos) DECLARE EXTERNAL SUB setParticle LET ipp = ii LET a = lattice/2.0 FOR i=0 TO 2*nx-1 FOR j=0 TO 2*ny-1 LET x = xPos + a*i LET y = yPos + a*j IF MOD(i+j,2)=0 THEN LET knd = knd1 ELSE LET knd = knd2 END IF CALL setParticle(ipp, knd, x, y) LET ipp = ipp + 1 NEXT j NEXT i LET setNaClTypeBlock = ipp - 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) = ion_m(knd) LET kind(i) = knd END SUB EXTERNAL SUB setForceTable DECLARE EXTERNAL FUNCTION cutoff LET eForceConst = 1.0/(4.0*PI*8.8542e-12) !epsilon0=8.8542e-12 FOR ki=0 TO 17 FOR kj=0 TO 17 LET a = ion_a(ki)+ion_a(kj) LET b = ion_b(ki)+ion_b(kj) LET c = ion_c(ki)*ion_c(kj)*1.0e-60 LET zi = ion_z(ki) LET zj = ion_z(kj) FOR ir=10 TO 1001 LET r = ir*hh LET fc = eForceConst*zi*zj*(EXP(-r/6.5e-10)/r)*(1.0/r+1.0/6.5e-10) !Debye-screened Coulomb force LET fr = 6.9742e-11*EXP((a-r)/b) !repulsive force LET fa = -6.9742e-21*(c/(r*r*r*r*r*r)) !attractive force LET forceTable(ki,kj,ir) = cutoff(r)*(fc + fr + fa) 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 :ion_str$(ionKind1)&","&ion_str$(ionKind2) PLOT TEXT, AT 150, 55 ,USING "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 ion - molecular dynamics" PLOT TEXT, AT 50,480 :"[esc] exit [space]pause/go [D]change draw mode" PLOT TEXT, AT 50,460 :"[1]NaCl [2]MgO [3]CaO [4]BaO [5]NaF [6]KF [7]KCl" PLOT TEXT, AT 50,440 :"[T] toggle 0:adiabatic/1:constant temperature" PLOT TEXT, AT 50,420 :"temp control -> [J]coolDown/stop [k]heatUp/stop" 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 SET LINE COLOR ion_col(kind(i)) DRAW circle WITH SCALE(ion_r(kind(i))*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 = (ion_r(kind(i))+ion_r(kind(j))) IF r<1.1*r0 AND kind(i)<>kind(j) THEN SET LINE COLOR 8 !gray 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)=ionKind1 THEN SET LINE COLOR 4 ELSE SET LINE COLOR 2 !4:red, 2:blue DRAW circle WITH SCALE(5)*SHIFT(vx(i)*mag+boxSize/2+xp,vy(i)*mag+boxSize/2+yp) NEXT i END sub END MODULE 
|