|
|
> No.4237[元記事へ]
mikeさんへのお返事です。
! 変なリストを出した お詫びに、貴方のプログラムを、無断で、
! ざっと、複素数モードの 3D に、移して見ました。(２進モードより速い)
! 何もかも、変更、上書きされ、気分悪いですね、ご容赦・・
!-------------------------------------------
! modify 2017.1.29
!
OPTION ARITHMETIC COMPLEX
DECLARE FUNCTION sysTemp
!
LET side=3
LET Nmt=side^3 ! number of particles
LET vss = 1.0E-9 ! m : view scale step
LET xMax = side * vss ! : x-Box size
LET yMax = side * vss ! : y-Box size
LET zMax = side * vss ! : z-Box size
!
DIM q1(Nmt),q2(Nmt), Dpxy(Nmt),Dvxy(Nmt),Dfxy(Nmt), apex(8,3),Dape(8,2),Dlin(12,2)
DIM Vi(4),Vo(4),Axys(4,4),Mxys(4,4),rotx(4,4),shxyz(4,4),shres(4,4)
DIM pxy(Nmt),pz(Nmt), vxy(Nmt),vz(Nmt), fxy(Nmt),fz(Nmt)
!
! pxy(),pz() : position of i-th paryicle
! vxy(),vz() : velocity of i-th paryicle
! fxy(),fz() : force of i-th paryicle
!
LET temp=150 ! K : initial temperature
LET sysTime = 0 ! s : system time
LET dt =0.02e-12 ! s : time step
!
LET mass = 39.95*1.67e-27 ! kg : mass of Ar
LET sigm = 3.40e-10 ! m : Lennard-Jones potential sigma for Ar
LET epsi = 1.67e-21 ! J : Lennard-Jones potential epsilon for Ar
LET k_ = 1.38e-23 ! J/K : Boltzman's constant
!
LET a_ = 0.5*dt/mass !s/kg : ⊿(m/s) / (kg・m/s^2) trans force to velocity
!
LET vmag = 100 ! view magnitude of velocity & force
LET t1 = vmag*dt ! view size /velocity
LET t2 = 0.1*t1^2 /mass ! view size /force
LET note0 = COMPLEX( xMax*0.9, yMax*1.35 ) ! view Notes Anchor
LET text0 = COMPLEX( -0.9*vss, -1.0*vss ) ! view Text Anchor
!
!--------------
MAT Axys=IDN
MAT rotx=IDN
MAT shxyz=IDN
MAT shres=IDN
LET Vi(4)=1
CALL mat_shxyz(shxyz,-xMax/2,-yMax/2,-zMax/2) !重心を原点へ移動する行列 shxyz 作成
CALL mat_shxyz(shres, xMax/2, yMax/2, zMax/2) !重心を原点から元へ戻す行列 shres 作成
!
LET Ax=-PI/2.5 !開始のｚ軸方向( 画面垂直０度からｘ軸回転成分)
LET Ay=0 ! 〃 〃 ( 〃 〃 ｙ軸回転成分)
LET Az=PI/6 !ｚ軸回転成分
!
!-----------------------
LET b=2*vss
SET WINDOW -b, xMax+b, -b, yMax+b
CALL setInitialCondition
CALL drawParticles
!
PRINT USING "time =####.## (ps) temp =####.## (K)" :sysTime*1e12, sysTemp
LET t0=TIME
DO
FOR i_=1 TO 10
IF mlb=0 THEN CALL moveParticles
NEXT i_
CALL drawParticles
!------------------ print / 2.0 psec.
IF mlb=0 AND MOD( ROUND(sysTime*10e12), 20)=0 THEN PRINT USING "time =####.## (ps) temp =####.## (K)" :sysTime*1e12, sysTemp
mouse poll mx,my,mlb,mrb
!--
WAIT DELAY t22 !t22: 制御出力の休止秒。
LET t11=TIME !t11: 直前の周期の終り。
LET t22=MAX(0,t22+(.08-MOD(t11-t00,86400))/20) !80ms－検出周期(t11-t00)＝偏差 →t22(積分 Gain=1/20)
LET t00=t11 !t00: 次の周期の始め＝ 前の周期の終り
!--
LOOP UNTIL mrb=1 !Right Click to Stop
PRINT TIME-t0
SUB setInitialCondition
! set particles
RANDOMIZE 5
LET s = xMax/(side+1)
LET n=0
FOR z=1 TO side
FOR y=1 TO side
FOR x=1 TO side
LET n=n+1
LET pxy(n) = s*COMPLEX( x, y)
LET pz(n) = s*COMPLEX( 0, z)
LET vxy(n) = 500*COMPLEX( RND-0.5, RND-0.5)
LET vz(n) = 500*COMPLEX( 0, RND-0.5)
LET fxy(n) = 0
LET fz(n) = 0
NEXT x
NEXT y
NEXT z
!-- ajust V to temp --
LET w=SQR( temp / sysTemp)
MAT vxy= w*vxy
MAT vz= w*vz
!-- wall cube apex --
LET i=1
FOR z=0 TO zMax STEP zMax
LET apex(i,1)=0
LET apex(i,2)=0
LET apex(i,3)=z
LET i=i+1
LET apex(i,1)=xMax
LET apex(i,2)=0
LET apex(i,3)=z
LET i=i+1
LET apex(i,1)=xMax
LET apex(i,2)=yMax
LET apex(i,3)=z
LET i=i+1
LET apex(i,1)=0
LET apex(i,2)=yMax
LET apex(i,3)=z
LET i=i+1
NEXT z
END SUB
FUNCTION sysTemp
LET v2 = 0
FOR i=1 TO Nmt
LET v2 = v2 + ABS(vxy(i))^2-vz(i)^2 ! vx^2 + vy^2 + vz^2
!LET v2 = v2 + ABS(ABS(vxy(i))+vz(i))^2 ! vx^2 + vy^2 + vz^2
NEXT i
LET sysTemp = mass*v2/Nmt/3/k_ ! {∑(.5*m*v^2) /N} *2/3 /k
END FUNCTION
SUB moveParticles
FOR i=1 TO Nmt
LET vxy(i) = vxy(i)+a_*fxy(i) !a_ = 0.5*dt/mass
LET vz(i) = vz(i)+ a_*fz(i)
LET pxy(i) = pxy(i)+vxy(i)*dt
LET pz(i) = pz(i)+vz(i)*dt
NEXT i
CALL calcForce
FOR i=1 TO Nmt
LET vxy(i) = vxy(i)+a_*fxy(i) !a_ = 0.5*dt/mass
LET vz(i) = vz(i)+ a_*fz(i)
NEXT i
LET sysTime=sysTime+dt
LET Az=Az+dt*1e11
END SUB
SUB calcForce
MAT fxy=ZER
MAT fz=ZER
FOR i=1 TO Nmt-1
FOR j=i+1 TO Nmt
LET Rij = pxy(i)-pxy(j)
LET Rzij = pz(i)-pz(j)
!--
LET r= ABS( ABS(Rij)+Rzij )
LET r6 = (sigm/r)^6
LET f = 24*epsi*r6*(2*r6-1)/r ! f(r) = (-dφ(r)/dr)
!--
LET fxy(i) = fxy(i) + f*Rij/r
LET fz(i) = fz(i) + f*Rzij/r
LET fxy(j) = fxy(j) - f*Rij/r
LET fz(j) = fz(j) - f*Rzij/r
NEXT j
NEXT i
!--
LET s = sigm *0.5
FOR i=1 TO Nmt ! boundary force
CALL force2( -s, -s, -s)
CALL force2( xMax+s, yMax+s, zMax+s)
NEXT i
END SUB
SUB force2( wx, wy, wz)
LET xy = pxy(i) - COMPLEX( wx, wy)
LET r = re(xy)
LET j = im(xy)
LET z = im( pz(i) - COMPLEX( 0, wz) )
LET r6 = (sigm/r)^6
LET j6 = (sigm/j)^6
LET z6 = (sigm/z)^6
LET fxy(i) = fxy(i) + 24*epsi*COMPLEX( r6*(2*r6-1)/r, j6*(2*j6-1)/j)
LET fz(i) = fz(i) + 24*epsi*COMPLEX( 0, z6*(2*z6-1)/z)
END SUB
SUB edge(z1,z2, p1,p2) !z1,z2：両端のｚ座標。 p1,p2：両端の x+yi 座標
IF z1+z2< zMAx THEN
PLOT LINES: p1;p2 !draw Far edge ：最初に描く。
ELSE
LET j9=j9+1
LET Dlin(j9,1)=p1 !save Near edge ：最後に描く。
LET Dlin(j9,2)=p2
END IF
END SUB
! display
SUB drawParticles
SET DRAW MODE HIDDEN
CLEAR
CALL control_
!--
FOR i=1 TO 8 !rotate box 8 apex
LET Vi(1)=apex(i,1)
LET Vi(2)=apex(i,2)
LET Vi(3)=apex(i,3)
MAT Vo=Vi*Mxys
LET Dape(i,1)=COMPLEX(Vo(1),Vo(2))
LET Dape(i,2)=Vo(3)
NEXT i
!--
SET LINE width 2
LET j9=0 !priority box 12 lines
FOR i=1 TO 3
CALL edge(Dape(i ,2),Dape(i+1,2), Dape(i ,1),Dape(i+1,1))
CALL edge(Dape(i ,2),Dape(i+4,2), Dape(i ,1),Dape(i+4,1))
CALL edge(Dape(i+4,2),Dape(i+5,2), Dape(i+4,1),Dape(i+5,1))
NEXT i
CALL edge(Dape(i ,2),Dape(i-3,2), Dape(i ,1),Dape(i-3,1))
CALL edge(Dape(i ,2),Dape(i+4,2), Dape(i ,1),Dape(i+4,1))
CALL edge(Dape(i+4,2),Dape(i+1,2), Dape(i+4,1),Dape(i+1,1))
SET LINE width 1
!--
FOR i=1 TO Nmt !rotate particles Nmt points
LET w=pxy(i)+fxy(i)*t2
LET Vi(1)=re(w)
LET Vi(2)=im(w)
LET Vi(3)=im( pz(i)+fz(i)*t2 )
MAT Vo=Vi*Mxys
LET Dfxy(i)=COMPLEX(Vo(1),Vo(2))
!--
LET w=pxy(i)+vxy(i)*t1
LET Vi(1)=re(w)
LET Vi(2)=im(w)
LET Vi(3)=im( pz(i)+vz(i)*t1 )
MAT Vo=Vi*Mxys
LET Dvxy(i)=COMPLEX(Vo(1),Vo(2))
!--
LET Vi(1)=re(pxy(i))
LET Vi(2)=im(pxy(i))
LET Vi(3)=im(pz(i))
MAT Vo=Vi*Mxys
LET Dpxy(i)=COMPLEX(Vo(1),Vo(2))
!--
LET q1(i)=Vo(3)
LET q2(i)=i
NEXT i
CALL Qsort00(1,Nmt) !priority particles Nmt points
!--
FOR n=1 TO Nmt !draw particles Nmt disk,verocity,force
LET i=q2(n) !i= ｚ最小(奥) からの particle 番号。
SET AREA COLOR i
DRAW disk WITH SCALE(.5*sigm)*SHIFT(Dpxy(i))
SET LINE COLOR 4 ! red : velocity
PLOT LINES: Dpxy(i); Dvxy(i)
SET LINE COLOR 1 ! black : force
PLOT LINES: Dpxy(i); Dfxy(i)
NEXT n
!--
SET LINE width 2
FOR i=1 TO j9 !draw Near edge
PLOT LINES: Dlin(i,1); Dlin(i,2)
NEXT i
SET LINE width 1
!--
SET COLOR 4 ! red : velocity line & text
PLOT LINES: note0; note0+vss*COMPLEX(0.5, 0)
PLOT label,AT note0+vss*COMPLEX(0.7 ,-0.1) : "velosity"
SET COLOR 1 ! black : force line & text
PLOT LINES: note0+vss*COMPLEX(0,-0.3); note0+vss*COMPLEX(0.5 ,-0.3)
PLOT label,AT note0+vss*COMPLEX(0.7,-0.4): "force"
!--
PLOT label, AT text0 ,USING "box size =##.# x##.# x##.# (nm)" :xMax*1e9,yMax*1e9,zMax*1e9
PLOT label, AT text0-vss*COMPLEX(0,0.30),USING "time =####.## (ps) temp =####.## (K)" :sysTime*1e12, sysTemp
PLOT label, AT text0-vss*COMPLEX(0,0.6 ) :"Ar in the box (3-dimensional molecular dynamics)"
!--
SET DRAW MODE EXPLICIT
END SUB
SUB control_
PLOT label,AT -vss, yMax+1.5*vss:"左 click 一時停止、drag 手動回転。右 click 終了。"
!-----click_drag-----
IF mlb=1 THEN
LET Ax= -(my-mybak)*PI/2E-9 !ドラッグ方向から、軸方向と回転量
LET Ay= +(mx-mxbak)*PI/2E-9
END IF
LET mxbak=mx
LET mybak=my
!-----
LET ar0=SQR(Ax^2+Ay^2) !回転の角度(∝マウス・ドラッグの長さ)
IF ar0<>0 THEN
LET DIRar0=ANGLE(Ax,Ay) !軸の角度
CALL mat_rotx(rotx, ar0)
MAT Axys=Axys*ROTATE(-DIRar0)*rotx*ROTATE(DIRar0) !ドラッグ累積 (方向,回転)
LET Ax=0
LET Ay=0
END IF
MAT Mxys=shxyz*ROTATE(Az)*Axys*shres
END SUB
!---------------------------------
! ｘ軸で 回転する行列 → 配列引数
!(x,y,z,1)| 1, 0, 0, 0 |
! | 0, cos(a), sin(a), 0 |
! | 0,-sin(a), cos(a), 0 |
! | 0, 0, 0, 1 |
!---------------------------------
SUB mat_rotx(m(,), a)
LET m(2,2)=COS(a)
LET m(3,2)=-SIN(a)
LET m(2,3)=SIN(a)
LET m(3,3)=COS(a) !他の要素は、呼出し側で管理
END SUB
!-----------------------------
! 平行移動。(sx,sy,sz)
!(x,y,z,1)| 1, 0, 0, 0 |
! | 0, 1, 0, 0 |
! | 0, 0, 1, 0 |
! | sx, sy, sz, 1 |
!-----------------------------
SUB mat_shxyz(m(,), sx,sy,sz)
LET m(4,1)=sx
LET m(4,2)=sy
LET m(4,3)=sz !他の要素は、呼出し側で管理
END SUB
!---------------------------
! Quick Sort q2() by q1()
!---------------------------
SUB Qsort00(L,R) !昇順にセット。
local i,j
LET i=L
LET j=R
LET w=q1((L+R)/2)
DO
DO WHILE q1(i)< w ![< ]昇順 [>]降順
LET i=i+1
LOOP
DO WHILE w< q1(j) ![< ]昇順 [>]降順
LET j=j-1
LOOP
IF j< i THEN EXIT DO !等号付 j<=i は、暴走。
SWAP q1(i),q1(j)
SWAP q2(i),q2(j)
LET i=i+1
LET j=j-1
LOOP UNTIL j< i !等号付 j<=i は、低速。
IF L< j THEN CALL Qsort00(L,j)
IF i< R THEN CALL Qsort00(i,R)
END SUB
END
|
|