|
投稿者：mike
投稿日：2017年 1月25日(水)10時35分29秒
|
|
|
|
実変数xと時間tの複素関数psi(x,t)をグラフィック画面上plotするため、
3次元の直交座標系(x,y,z)に
x -> x
Re(psi(x,t0)) -> y
Im(psi(x,t0)) -> z
を対応させ、これをグラフィック画面上に(x',y')の所にplotし、
時間tを少しずつ動かして、動画として表現したいのです。
見る方位の単位ベクトル（ex,ey,ez)あるいは(1,theta,fai)が与えられた時、
(x,y,z)から(x',y')への写像は一意に決まると思うのですが、
具体的にどのように考えたらよいか教えてください。
x' = f(x,y,z,ex,ey,ez)
y' = g(x,y,z,ex,ey,ez)
よろしくお願いします。
|
|
|
|
投稿者：白石 和夫
投稿日：2017年 1月25日(水)15時51分26秒
|
|
|
|
投稿者：mike
投稿日：2017年 1月26日(木)07時13分5秒
|
|
|
|
> No.4238[元記事へ]
白石 和夫さんへのお返事です。
白石先生、プログラムをありがとうございます。
MATの計算で MAT POINT=POINT*m のように左からベクトルPOINTをかけるということは
POINTは行ベクトルなのでしょうか？
また、 MAT m=TRANSFORM ステートメントは the current transformation matrix can be obtained
とありますので、MAT m=TRANSFORMに代入されのは時点の変換行列ということになりますが
どこで定義されたものなのでしょうか？
電子波束psi(x,t)が時間発展しながらポテンシャル中を動くプログラムを書いています。
以下の電子波束の動くプログラムを動かしますと、青色の曲線が複素関数psi(x,t)を表すのですが、
ご覧のように一見3D風ですが、かなり違和感があります。これを改善したいのです。
プログラム中で「!****************」を付けたところがインチキ3Dの部分です。
! ========= Quantum Electron Dynamics 1D ==========
!
! wavePacketQED1D.bas
! Mitsuru Ikeuchi
! ver 0.0.0 2017.01.25
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB QED1D.setInitialCondition,QED1D.evolveWave,QED1D.drawWave
DECLARE EXTERNAL FUNCTION QED1D.systemTime,QED1D.norm
!
CALL setInitialCondition(1) !(menu)
!
FOR it=1 TO 1000
FOR i=1 TO 10
CALL evolveWave
NEXT i
CALL drawWave
!IF (MOD(it,10)=0) THEN
! PRINT USING "time =#####.## (au) norm =#.###############":systemTime,norm
!END IF
NEXT it
!
END
! ---------- QED1D module (QED: Quantum Electron Dynamics) ----------
!
! - time dependent Schrodinger equation: i(d/dt)psi(r,t) = H psi(r,t)
! - time evolution
! psi(r,t+dt) = exp(-i dt H) psi(r,t), (H:Hamiltonian of the system)
! H = -delta/2 + V(r), delta = d^2/dx^2 + d^2/dy^2 + d^2/dz^2
! psi(r,t+dt) = exp(-i dt H) psi(r,t) nearly=
! {exp(-i(dt/2)V} {exp(i dt(delta/2)} {exp(-i(dt/2)V} psi(r,t)
! - algorism: {exp(i dt(delta/2)}
! QED: Watanabe's algorism (semi-implicit method)
! Naoki Watanabe, Masaru Tsukada; arXiv:physics/0011068v1
! (Published from Physical Review E. 62, 2914, (2000).)
!
MODULE QED1D
MODULE OPTION ARITHMETIC NATIVE
OPTION BASE 0
PUBLIC SUB setInitialCondition, evolveWave, drawWave
PUBLIC FUNCTION systemTime, norm
SHARE NUMERIC sysTime, dt, dx, auLength, auTime, auEnergy, NNx
SHARE NUMERIC phRe(500),phIm(500) ! wave function
SHARE NUMERIC vv(500) ! external potential
SHARE NUMERIC bRe(500),bIm(500) ! b vector in kxStep
SHARE NUMERIC uRe(500),uIm(500) ! u vector in kxStep
LET sysTime = 0.0 ! (au) au : atomic unit (hbar=1, me=1, e=1)
LET dx = 0.5 ! (au) x division
LET dt = 1.0*dx*dx ! (au) time division dt/(dx*dx)>2 ~ unstable
LET auLength = 5.29177e-11 ! (m) 1(au) = auLength (m)
LET auTime = 2.41888e-17 ! (s) 1(au) = auTime (s)
LET auEnergy = 4.38975e-18 ! (J) 1(au) = auEnergy (J) (= 27.2114 eV)
LET NNx = 400 ! x division
!
! ----- set initial condition
!
EXTERNAL SUB setInitialCondition(menu) !public
DECLARE EXTERNAL SUB setGaussianWave,setparabolicPotential
LET sysTime = 0.0
IF menu=1 THEN
CALL setGaussianWave(50.0,10.0,0.0) !(wavePos,waveWidth,momentum)
CALL setparabolicPotential(0.5*(NNx-1)*dx,0.0001) ! v(x)=k*(x-x0)^2
ELSEIF menu=2 THEN
CALL setGaussianWave(50.0,10.0,1.0) !(wavePos,waveWidth,momentum)
CALL setparabolicPotential(0.5*(NNx-1)*dx,0.0) ! v(x)=k*(x-x0)^2
END IF
! set window
LET w = NNx*dx
LET b = 20
SET WINDOW -b,w+b,-b,w+b
END SUB
!
EXTERNAL SUB setGaussianWave(wavePos,waveWidth,momentum)
DECLARE EXTERNAL SUB normalize
FOR i=1 TO NNx-2
LET x = i*dx
LET phAb = EXP(-((x-wavePos)/(2.0*waveWidth))^2)
LET phPh = momentum*x
LET phRe(i) = phAb*cos(phPh)
LET phIm(i) = phAb*sin(phPh)
NEXT i
LET phRe(0) = 0.0
LET phIm(0) = 0.0
LET phRe(NNx-1) = 0.0
LET phIm(NNx-1) = 0.0
CALL normalize
END SUB
!
EXTERNAL SUB setparabolicPotential(x0,k) ! v(x)=k*(x-x0)^2
FOR i=0 TO NNx-1
LET x = i*dx
LET vv(i) = k*(x-x0)^2
NEXT i
END SUB
!
! ----- evolve Wave
!
EXTERNAL SUB evolveWave !piblic
DECLARE EXTERNAL SUB phaseStep,kxStep
CALL phaseStep
CALL kxStep
CALL phaseStep
LET sysTime = sysTime + dt
END SUB
!
EXTERNAL SUB kxStep !{exp(i dt(delta/2)} semi-implicit method
LET ai = 4*dx*dx/dt
LET aaAb = 4+ai*ai
FOR i=1 TO NNx-2 !set b-vector
LET bRe(i) = 2*phRe(i)-ai*phIm(i) - phRe(i+1) - phRe(i-1)
LET bIm(i) = 2*phIm(i)+ai*phRe(i) - phIm(i+1) - phIm(i-1)
NEXT i
! forward elimination
LET uRe(1) = -2/aaAb
LET uIm(1) = -ai/aaAb
LET phRe(1) = bRe(1)*uRe(1) - bIm(1)*uIm(1)
LET phIm(1) = bIm(1)*uRe(1) + bRe(1)*uIm(1)
FOR i=2 TO NNx-2
LET auAb = (-2-uRe(i-1))*(-2-uRe(i-1))+(ai-uIm(i-1))*(ai-uIm(i-1))
LET uRe(i) = (-2-uRe(i-1))/auAb
LET uIm(i) = -(ai-uIm(i-1))/auAb
LET phRe(i) = (bRe(i)-phRe(i-1))*uRe(i) - (bIm(i)-phIm(i-1))*uIm(i)
LET phIm(i) = (bRe(i)-phRe(i-1))*uIm(i) + (bIm(i)-phIm(i-1))*uRe(i)
NEXT i
! backward substitution
FOR i=NNx-3 TO 1 STEP -1
LET phRe(i) = phRe(i) - (phRe(i+1)*uRe(i) - phIm(i+1)*uIm(i))
LET phIm(i) = phIm(i) - (phRe(i+1)*uIm(i) + phIm(i+1)*uRe(i))
NEXT i
END SUB
!
EXTERNAL SUB phaseStep !{exp(-i(dt/2)V} evolve 0.5dt
FOR i=1 TO NNx-2
LET th = 0.5*dt*vv(i)
LET cs = COS(th)
LET sn = SIN(th)
LET phr = phRe(i)
LET phi = phIm(i)
LET phRe(i) = cs*phr+sn*phi
LET phIm(i) = cs*phi-sn*phr
NEXT i
END SUB
!
! ----- utility
!
EXTERNAL FUNCTION systemTime !public
LET systemTime = sysTime
END FUNCTION
!
EXTERNAL SUB normalize
LET a = norm !public function
FOR i=1 TO NNx-2
LET phRe(i) = phRe(i)/a
LET phIm(i) = phIm(i)/a
NEXT i
END SUB
!
EXTERNAL FUNCTION norm !public <psi|psi>
!LOCAL i,a2
LET a2 = 0.0
FOR i=1 TO NNx-2
LET a2 = a2 + (phRe(i)*phRe(i)+phIm(i)*phIm(i))*dx
NEXT i
LET norm = SQR(a2)
END FUNCTION
!
EXTERNAL FUNCTION kineticEnergy !<psi| d^2/dx^2 |psi>
LET s = 0.0
FOR i=1 TO NNx-2
LET hphRe = (2.0*phRe(i)-phRe(i+1)-phRe(i-1))/(2.0*dx*dx)
LET hphIm = (2.0*phIm(i)-phIm(i+1)-phIm(i-1))/(2.0*dx*dx)
LET s = s + (phRe(i)*hphRe + phIm(i)*hphIm)*dx
NEXT i
LET kineticEnergy = s
END FUNCTION
!
EXTERNAL FUNCTION potentialEnergy !<psi| V(x) |psi>
LET s = 0.0
FOR i=1 TO NNx-2
LET s = s + vv(i)*(phRe(i)*phRe(i)+phIm(i)*phIm(i))*dx
NEXT i
LET potentialEnergy = s
END FUNCTION
!
! ----- draw wave
!
EXTERNAL SUB drawWave !LOCAL i,yh,sth,cth,ph2,pe,ke
DECLARE EXTERNAL FUNCTION kineticEnergy,potentialEnergy
LET yh = (NNx*dx+40)*0.4
LET sth = SIN(PI/6) !****************
LET cth = COS(PI/6) !****************
SET DRAW MODE HIDDEN
CLEAR
!plot V(x),{psi>,<psi|psi>
SET LINE COLOR 1 ! black : PLOT potential V(x);
FOR i=0 TO NNx-1
PLOT LINES: dx*i, vv(i)*100+yh;
NEXT i
PLOT LINES
SET LINE COLOR 4 ! red : plot probability <psi|psi>
FOR i=0 TO NNx-1
LET ph2 = (phRe(i)*phRe(i)+phIm(i)*phIm(i))*dx
PLOT LINES: dx*i, ph2*1000+yh;
NEXT i
PLOT LINES
SET LINE COLOR 2 ! blue : PLOT wave function |psi> !****************
FOR i=0 TO NNx-1
PLOT LINES: dx*i-sth*phIm(i)*50, (phRe(i)-cth*phIm(i))*50+yh; !****************
NEXT i
PLOT LINES
!caption
LET ke = kineticEnergy
LET pe = potentialEnergy
SET TEXT COLOR 1 ! black
PLOT TEXT, AT 0,48 ,USING "time =#####.## (au) =####.#### (femto s)":sysTime,sysTime*auTime*1e15
PLOT TEXT, AT 0,40 ,USING "norm of wave function =#.###############":norm
PLOT TEXT, AT 0,32 ,USING "kinetic energy =#.###### (au) =##.### (eV)":ke,ke*27.2114
PLOT TEXT, AT 0,24 ,USING "potential energy =#.###### (au) =##.### (eV)":pe,pe*27.2114
PLOT TEXT, AT 0,16 ,USING "total energy =#.###### (au) =##.### (eV)":ke+pe,(ke+pe)*27.2114
PLOT TEXT, AT 0, 8 ,USING "xLength =####.# (au) =###.### (nm)":NNx*dx,NNx*dx*auLength*1e9
PLOT TEXT, AT 0, 0 :"wave packet - quantum electron dynamics 1D"
!
SET TEXT COLOR 1 ! black
PLOT TEXT, AT 0, yh+22 :"potential V(x)"
SET TEXT COLOR 4 ! red
PLOT TEXT, AT 0, yh+6 :"<psi|psi>"
SET TEXT COLOR 2 ! blue
PLOT TEXT, AT 0, yh-10 :"|psi(x,t)>"
SET DRAW MODE EXPLICIT
END SUB
!
END MODULE
|
|
|
|
投稿者：白石和夫
投稿日：2017年 1月26日(木)08時52分51秒
|
|
|
|
> No.4239[元記事へ]
> MATの計算で MAT POINT=POINT*m のように左からベクトルPOINTをかけるということは
> POINTは行ベクトルなのでしょうか？
Full BASICの座標変換は行ベクトルです。左から右に演算を行う流儀です。
> また、 MAT m=TRANSFORM ステートメントは the current transformation matrix can be obtained
> とありますので、MAT m=TRANSFORMに代入されのは時点の変換行列ということになりますが
> どこで定義されたものなのでしょうか？
変換行列は、DRAW文でPICTUREを呼び出したときにBASICの内部で定義されます。
PICTUREを再帰呼び出しすると、変換行列の積が計算されます。
行列による変形
変形指示MAT文
|
|
|
|
投稿者：mike
投稿日：2017年 1月26日(木)12時22分7秒
|
|
|
|
白石和夫さんへのお返事です。
白石先生、ご回答、ありがとうございます。
> Full BASICの座標変換は行ベクトルです。左から右に演算を行う流儀です。
BASICの式評価で、演算の優先順位が同じ場合、左から行われる原則に合わせているのでしょうか。
今の懸案について、
(x’,y’,z’,1) = (x,y,z,1)*(x軸th回転行列)*(y軸fi回転行列) の計算を成分ごとに行い、書き下すと
x’ = COS(fi)*x+SIN(fi)*(SIN(th)*y+COS(th)*z)
y’ = COS(th)*y-SIN(th)*z
z’ = -SIN(fi)*x+COS(fi)*(SIN(th)*y+COS(th)*z)
となりましたので、下記のように関数定義しました。
DEF xRot(x,y,z,th,fi) = COS(fi)*x+SIN(fi)*(SIN(th)*y+COS(th)*z)
DEF yRot(x,y,z,th,fi) = COS(th)*y-SIN(th)*z
!DEF zRot(x,y,z,th,fi) = -SIN(fi)*x+COS(fi)*(SIN(th)*y+COS(th)*z)
使うときは、下記のようにしました。
LET rx = xRot(dx*i-x0,phRe(i)*50,phIm(i)*50,0,fi)
LET ry = yRot(dx*i-x0,phRe(i)*50,phIm(i)*50,0,fi)
!LET rz = zRot(dx*i-x0,phRe(i)*50,phIm(i)*50,0,fi)
PLOT LINES: rx+x0,ry+yh;
結果は、まともに見えますが、残念ながら修正前の見栄えと同じように思えます。
下記に電子波束psi(x,t)が時間発展しながらポテンシャル中を動くプログラムの修正部分を載せます。
!
! ----- draw wave
!
EXTERNAL SUB drawWave
DECLARE EXTERNAL FUNCTION kineticEnergy,potentialEnergy
!rotate th around x-axis and rotate fi around y-axis !*************
DEF xRot(x,y,z,th,fi) = COS(fi)*x+SIN(fi)*(SIN(th)*y+COS(th)*z) !*************
DEF yRot(x,y,z,th,fi) = COS(th)*y-SIN(th)*z !*************
!DEF zRot(x,y,z,th,fi) = -SIN(fi)*x+COS(fi)*(SIN(th)*y+COS(th)*z) !*************
LET x0 = dx*NNx/2 !LET y0 = 0
LET yh = (NNx*dx+40)*0.4
SET DRAW MODE HIDDEN
CLEAR
!plot V(x),{psi>,<psi|psi>
SET LINE COLOR 1 ! black : PLOT potential V(x);
FOR i=0 TO NNx-1
PLOT LINES: dx*i, vv(i)*100+yh;
NEXT i
PLOT LINES
SET LINE COLOR 4 ! red : plot probability <psi|psi>
FOR i=0 TO NNx-1
LET ph2 = (phRe(i)*phRe(i)+phIm(i)*phIm(i))*dx
PLOT LINES: dx*i, ph2*1000+yh;
NEXT i
PLOT LINES
SET LINE COLOR 2 ! blue : PLOT wave function |psi>
LET fi = PI/12 !*************
FOR i=0 TO NNx-1
LET rx = xRot(dx*i-x0,phRe(i)*50,phIm(i)*50,0,fi) !*************
LET ry = yRot(dx*i-x0,phRe(i)*50,phIm(i)*50,0,fi) !*************
!LET rz = zRot(dx*i-x0,phRe(i)*50,phIm(i)*50,0,fi) !*************
PLOT LINES: rx+x0,ry+yh; !*************
NEXT i
PLOT LINES
!caption
LET ke = kineticEnergy
LET pe = potentialEnergy
SET TEXT COLOR 1 ! black
PLOT TEXT, AT 0,48 ,USING "time =#####.## (au) =####.#### (femto s)":sysTime,sysTime*auTime*1e15
PLOT TEXT, AT 0,40 ,USING "norm of wave function =#.###############":norm
PLOT TEXT, AT 0,32 ,USING "kinetic energy =#.###### (au) =##.### (eV)":ke,ke*27.2114
PLOT TEXT, AT 0,24 ,USING "potential energy =#.###### (au) =##.### (eV)":pe,pe*27.2114
PLOT TEXT, AT 0,16 ,USING "total energy =#.###### (au) =##.### (eV)":ke+pe,(ke+pe)*27.2114
PLOT TEXT, AT 0, 8 ,USING "xLength =####.# (au) =###.### (nm)":NNx*dx,NNx*dx*auLength*1e9
PLOT TEXT, AT 0, 0 :"wave packet - quantum electron dynamics 1D"
!
SET TEXT COLOR 1 ! black
PLOT TEXT, AT 0, yh+22 :"potential V(x)"
SET TEXT COLOR 4 ! red
PLOT TEXT, AT 0, yh+6 :"<psi|psi>"
SET TEXT COLOR 2 ! blue
PLOT TEXT, AT 0, yh-10 :"|psi(x,t)>"
SET DRAW MODE EXPLICIT
END SUB
|
|
|
|
投稿者：SECOND
投稿日：2017年 1月30日(月)19時12分47秒
|
|
|
|
> 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
|
|
|
|
投稿者：mike
投稿日：2017年 1月31日(火)05時57分22秒
|
|
|
|
SECONDさんへのお返事です。
SECONDさん、こんにちは。
>! 変なリストを出した お詫びに、貴方のプログラムを、無断で、
>! ざっと、複素数モードの 3D に、移して見ました。
複素数を使ったエレガントなコードですね。３Dの表現に感動しました。
まだ私は複素数を使いこなせていないので、熟読し、勉強させていただきます。
コードを掲載していただき、ありがとうございました。
>(２進モードより速い)
この原因はどこにあるとお考えですか？ ご教示いただければ幸いです。
>! 何もかも、変更、上書きされ、気分悪いですね、ご容赦・・
とんでもありません。
雛形プログラムの公開により、プログラムが改良されていくのが望ましいと思います。
この後も、いくつか雛形プログラム（おもに物理シミュレーション）を公開していきますので、
改良していただけれま幸いです。
よろしくお願いします。
|
|
|
|
投稿者：mike
投稿日：2017年 1月31日(火)09時14分14秒
|
|
|
|
> No.4244[元記事へ]
SECONDさんへのお返事です。
SECONDさんこんにちは。
下記はSECONDさんのプログラムについてのコメントです。
（１）複素数で(x,y,z)を表すことの利点？
! pxy(),pz() : position of i-th paryicle
! vxy(),vz() : velocity of i-th paryicle
! fxy(),fz() : force of i-th paryicle
x、y、zの扱いが対称でないので、コードの見やすさを阻害しているように思えます。
複素数のBASIC組み込みの計算が使え、計算速度が上がるのでしょうか。
（２）表示間隔のコントロール
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: 次の周期の始め＝ 前の周期の終り
これはすばらしいアイデアですね。
(.08-MOD(t11-t00,86400))/20は負帰還に成っていそうですが、安定性については理解が及びません。
（３）ちょっと問題かも
SUB moveParticles
FOR i=1 TO Nmt
LET vxy(i) = vxy(i)+a_*fxy(i)
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
LET sysTime=sysTime+dt
LET Az=Az+dt*1e11
END SUB
これはオイラー法であり、誤差の累積で、真の解からだんだん離れていきます。
EXTERNAL SUB moveParticles ! velocity Verlet method
DECLARE EXTERNAL SUB calcForce
LET a = 0.5*dt/mass
FOR i=1 TO nMolec
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 vx(i) = vx(i)+a*ffx(i)
LET vy(i) = vy(i)+a*ffy(i)
NEXT i
LET sysTime=sysTime+dt
END SUB
は速度Verlet法であり、本来は
EXTERNAL SUB moveParticles ! velocity Verlet method
DECLARE EXTERNAL SUB calcForce
LET a = 0.5*dt/mass
FOR i=1 TO nMolec
LET vx(i) = vx(i)+a*ffx(i)
LET vy(i) = vy(i)+a*ffy(i)
LET xx(i) = xx(i)+vx(i)*dt !xx(t+0.5dt) <-- xx(t-0.5dt)
LET yy(i) = yy(i)+vy(i)*dt
NEXT i
LET sysTime=sysTime+0.5*dt
CALL calcForce ！ ffx,ffy <--xx(t+0.5dt),yy(t+0.5dt)
FOR i=1 TO nMolec
LET vx(i) = vx(i)+a*ffx(i)
LET vy(i) = vy(i)+a*ffy(i)
NEXT i
LET sysTime=sysTime+0.5*dt
END SUB
とすべきところですが、sysTimeは計算と直接関わらないので1つにまとめたため、誤解されたかもしれません。
速度Verlet法はシンプレックス性があり誤差の累積による真の解との乖離が小さいという利点があります。
たとえば下記をご覧ください。
（４）箱や中の粒子の3D表示については、まだ理解できていません
のでこれから解読させていただきます。
|
|
|
戻る