|
|
分子動力学法の雛形プログラムを公開します。
なるべくプログラムが読みやすいように、高速化や汎用性を狙わずに書いたつもりです。
分子動力学法は、それぞれの分子に掛かる力を現在の分子の位置と分子間のポテンシャルから計算し、
現在の分子の位置と速度と力からNewtonの運動方程式を使って短い時間dt後の位置や速度を求め、
これを繰り返すことで運動を追跡する手法です。（粒子数N一定、体積V一定、全エネルギE一定の系）
分子の位置や速度、力のデータから、いろいろな分子の統計量を求めますが、ここでは温度を表示しています。
温度が揺らいでいるのは、分子の位置関係が刻々と変わるため、分子間ポテンシャルが変化するためです。
表示の説明：
黒い四角は箱を表します。箱と分子の間にも力が働きます。
青色の丸がアルゴン分子で、その中心から引かれた赤線の長さは速度に比例します。
黒線の長さはその分子にかかる力の大きさに比例します。
試験環境：
本プログラムは十進BASIC 0.6.5.2 / macOS 10.7.5と
十進BASICBASIC Ver 7.7.8 / windows 10でテストしました。
------------------------
!
! ========= molecular dynamics 2D ==========
!
! ArGasMD2D001.bas
! Mitsuru Ikeuchi (C) Copyleft
!
! ver 0.0.1 2017.01.19 created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB md2d.setInitialCondition, md2d.moveParticles, md2d.drawParticles
DECLARE EXTERNAL FUNCTION md2d.systemTime, md2d.systemTemprature
!
CALL setInitialCondition
CALL drawParticles
!
FOR it=1 TO 1000
FOR i=1 TO 20
CALL moveParticles
NEXT i
CALL drawParticles
!IF (MOD(it,10)=0) THEN
! PRINT USING "time =####.## (ps) temp =####.## (K)":systemTime*1.0E12,systemTemprature
!END IF
NEXT it
!
END
!
! ---------- md2d module ----------
!
! method: velocity Verlet
! potential: Lennard-Jones V(r) = 4*epsilon*((sigma/r)^12-(sigma/r)^6)
! force F(r) = -dV(r)/dr
MODULE md2d
MODULE OPTION ARITHMETIC NATIVE
PUBLIC SUB setInitialCondition,moveParticles,drawParticles
PUBLIC FUNCTION systemTime, systemTemprature
SHARE NUMERIC sysTime, dt, nMolec, xMax, yMax, sigma, mass, epsilon, boxSize
SHARE NUMERIC xx(100),yy(100) ! (xx(i),yy(i)) : position of i-th particle
SHARE NUMERIC vx(100),vy(100) ! (vx(i),vy(i)) : velocity of i-th particle
SHARE NUMERIC ffx(100),ffy(100) ! (ffx(i),ffy(i)): total force of i-th particle
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.67e-27 ! mass of Ar (kg)
LET sigma = 3.40e-10 ! Lennard-Jones potential sigma for Ar (m)
LET epsilon = 1.67e-21 ! Lennard-Jones potential epsilon FOR Ar (J)
LET boxSize = 6 ! world box size in the graphic window
!
EXTERNAL SUB setInitialCondition
DECLARE EXTERNAL SUB ajustVelocity
! set particles
RANDOMIZE
LET s = xMax/7.0
FOR i=1 TO nMolec
LET xx(i) = s*MOD((i-1),6)+1.0*s
LET yy(i) = s*int((i-1)/6)+1.0*s
LET vx(i) = 500.0*(RND-0.5)
LET vy(i) = 500.0*(RND-0.5)
LET ffx(i) = 0.0
LET ffy(i) = 0.0
NEXT i
CALL ajustVelocity(150) !(temp)
! set window
LET b=2
SET WINDOW -b,boxSize+b,-b,boxSize+b
END SUB
!
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 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
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
FOR i=1 TO nMolec
LET ffx(i) = 0.0
LET ffy(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 rij = SQR(xij*xij+yij*yij)
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
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))
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 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)
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
!
! display
!
EXTERNAL SUB drawParticles
DECLARE EXTERNAL FUNCTION systemTemprature
LET mag = boxSize/xMax
LET vScate = 100*dt
LET fScale = 1000*dt*dt/mass
!
SET DRAW MODE HIDDEN
CLEAR
SET LINE COLOR 1 ! black : wall
PLOT LINES: 0,0; boxSize,0; boxSize,boxSize; 0,boxSize; 0,0
!draw Balls, velocity and force
FOR i=1 TO nMolec
SET COLOR 2 ! blue : molecule
DRAW circle WITH SCALE(sigma/2.0*mag)*SHIFT(xx(i)*mag,yy(i)*mag)
SET LINE COLOR 4 ! red : velocity
PLOT LINES: xx(i)*mag,yy(i)*mag;
PLOT LINES: (xx(i)+vx(i)*vScate)*mag,(yy(i)+vy(i)*vScate)*mag
SET LINE COLOR 1 ! black : force
PLOT LINES: xx(i)*mag,yy(i)*mag;
PLOT LINES: (xx(i)+ffx(i)*fScale)*mag,(yy(i)+ffy(i)*fScale)*mag
NEXT i
!draw caption
SET TEXT COLOR 1 ! black
PLOT TEXT, AT 0,-0.5 ,USING "time =####.## (ps) temp =####.## (K)":sysTime*1E12,systemTemprature
PLOT TEXT, AT 0,-1.0 ,USING "box size =##.# x ##.# (nm)":xMax*1e9,yMax*1e9
PLOT TEXT, AT 0,-1.5 :"Ar in the box (2-dimensional molecular dynamics)"
!
LET xp = boxSize*0.5+1
LET yp = boxSize+1
SET LINE COLOR 4 ! red : velocity
PLOT LINES: xp,yp;xp+0.5,yp
SET TEXT COLOR 4 ! red
PLOT TEXT, AT xp+0.7,yp-0.05: "velosity"
SET LINE COLOR 1 ! black : force
PLOT LINES: xp,yp-0.3;xp+0.5,yp-0.3
SET TEXT COLOR 1 ! black
PLOT TEXT, AT xp+0.7,yp-0.35: "force"
!
SET DRAW MODE EXPLICIT
END SUB
!
END MODULE
|
|