|
|
定常状態の電子のエネルギーと電子状態を求めるプログラム(008harmonicsSD1D.bas)を公開します。
定常状態の電子のエネルギーと電子状態を求める問題は、系のハミルトニアンをHとして H |i> = e_i |i>
の最小固有値とその近くの低いエネルギー固有値と固有関数を求める問題になります。
最急降下法（steepest descent method）によって、この最小固有値付近の固有値と固有関数を求めます。
適当に与えた関数（ここでは乱数NNx個の列）の組{|i>}から出発し、エネルギーei=<i|H|i>を求め、
エネルギーの減少する方向 (H-ei)|i> から 次世代の関数|i(next)> = |i> - damp(H-ei)|i>を求め、
関数の組を直交化することを繰り返す（逐次近似する）ことで小さい方からN個の固有値と固有関数を求めます。
本プログラムでは、方物型と井戸型のポテンシャルを選択できます。
表示の説明：
はじめに関数が上から降ってくるのは、近似を繰り返すごとにエネルギーが低下していく様子をみるため、
試行波動関数にエネルギーの下駄を履かせているためです。
逐次近似が進んでくると、固有値、固有関数（波動関数）はほとんど変わらなくなります。
方物型のとき、理論的には固有エネルギーは1次元の場合、(1/2)+n, n=0,1,2,.. となります。
右側には基底状態|0>と各励起状態|1>,|2>,..の内積を表示しています。16桁程度の精度で直交していることが解ります。
試験環境：
本プログラムは十進BASIC 6.6.0 / macOS 10.7.5, 十進BASIC Ver 7.7.8 / windows 10 でテストしました。
----------------
!
! ========= steepest descent method 1D ==========
!
! 008harmonicsSD1D.bas
! Mitsuru Ikeuchi (C) Copyleft
!
! ver 0.0.1 2017.03.06 created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB sd1d.setInitialCondition,sd1d.SDiteration,sd1d.drawState
DECLARE EXTERNAL FUNCTION INKEY$
LET stateMax = 10 !state 0,1,...,stateMax-1
LET vIndex = 0 !0:harmonic potential, 1:quantum well
LET iterMax = 10 ! 10 = iteration in SDiteration()
CALL setInitialCondition(stateMax,vIndex)
DO
CALL SDiteration(stateMax,iterMax)
CALL drawState
LET S$=INKEY$
IF S$="0" THEN
LET vIndex = 0
CALL setInitialCondition(stateMax,vIndex)
ELSEIF S$="1" THEN
LET vIndex = 1
CALL setInitialCondition(stateMax,vIndex)
ELSEIF S$="." THEN
EXIT DO
END IF
LOOP
END
EXTERNAL FUNCTION INKEY$
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION
! ---------- steepest descent method 1D ----------
!
! system Hamiltonian: H = -delta/2 + V(r) , delta r = div grad r
! eigen energy set { Ei }, eigen function set { |i> }
!
! procedure : successive approximation
! (i) trial function set { |0>,|1>,..,|i>,.. }
! (2) energy of |i> : ei = <i|H|i>/<i|i>
! (3) steepest gradient direction (H-ei)|i>
! (4) next generation : |i(next)> = |i> - dampingFactor*(H-ei)|i>
! (5) orthogonalization { |0>,|1>,..,|i>,.. } (Gram-Schmidt)
! (6) goto (2)
!
MODULE sd1d
MODULE OPTION ARITHMETIC NATIVE
OPTION BASE 0
PUBLIC SUB setInitialCondition, SDiteration, drawState
SHARE NUMERIC NNx, dx, iterCount
SHARE NUMERIC sdEnergy(20) ! electron state energy
SHARE NUMERIC sdState(20,400) ! electron states 0...20 0:ground state
SHARE NUMERIC wrk(400) ! state work space in steepestDescent
SHARE NUMERIC vv(400) ! external potential
LET NNx = 256 ! max number of sdState(,NNx,NNy)
LET dx = 1/16 ! (au) x-division
LET iterCount = 0 ! sd iteration count
! ---------- set initial condition
EXTERNAL SUB setInitialCondition(stateMax,vIndex) !public
DECLARE EXTERNAL SUB setInitialState,setPotential
LET iterCount = 0
CALL setInitialState(stateMax)
CALL setPotential(vIndex)
! set window
SET WINDOW 0,500, 0,500
END SUB
EXTERNAL SUB setInitialState(stateMax)
DECLARE EXTERNAL SUB normalizeState
RANDOMIZE
FOR ist=0 TO stateMax-1
FOR i=1 TO NNx-2
LET sdState(ist,i) = RND-0.5
NEXT i
LET sdState(ist,0) = 0
LET sdState(ist,NNx-1) = 0
CALL normalizeState(ist)
NEXT ist
END SUB
EXTERNAL SUB setPotential(vIndex)
LET x0 = 0.5*NNx*dx
IF vIndex=0 THEN !--- hermonic
FOR i=0 TO NNx-1
LET x = i*dx
LET vv(i) = MIN(0.5*(x-x0)^2,24.5)
NEXT i
ELSEIF vIndex=1 THEN !--- well
FOR i=0 TO NNx-1
LET x = i*dx
IF ABS(x-x0)<4 THEN LET vv(i) = 0 ELSE LET vv(i) = 18
NEXT i
END IF
END SUB
! ---------- steepest descent iteration
EXTERNAL SUB SDiteration(stateMax, iterMax) !public
DECLARE EXTERNAL FUNCTION steepestDescent
DECLARE EXTERNAL SUB GramSchmidt
LET damp = 0.003 !damping factor in steepest descent
FOR i=0 TO iterMax-1
FOR ist=0 TO stateMax-1
LET sdEnergy(ist) = steepestDescent(ist, damp)
NEXT ist
CALL GramSchmidt(stateMax)
LET iterCount = iterCount + 1
NEXT i
END SUB
EXTERNAL FUNCTION steepestDescent(ist,damp) !--- steepest descent method
DECLARE EXTERNAL FUNCTION energyOfState
DECLARE EXTERNAL SUB normalizeState
LET h2 = 2*dx*dx
LET ei = energyOfState(ist) !--- E_ist = <ist|H|ist>
FOR i=1 TO NNx-2 !--- |wrk> = (H-E_ist)|ist>
LET wrk(i) = (2*sdState(ist,i)-sdState(ist,i+1)-sdState(ist,i-1))/h2+(vv(i)-ei)*sdState(ist,i)
NEXT i
FOR i=1 TO NNx-2 !--- |ist(next)> = |ist> - damp*|wrk> ( norm(|ist(next)>) <>1 )
LET sdState(ist,i) = sdState(ist,i)-damp*wrk(i)
NEXT i
CALL normalizeState(ist)
LET steepestDescent = ei
END FUNCTION
EXTERNAL FUNCTION energyOfState(ist) !--- E_ist = <ist|H|ist>/<ist|ist>
LET h2 = 2*dx*dx
LET s = 0
LET sn=0
FOR i=1 TO NNx-2
LET s = s+sdState(ist,i)*((2*sdState(ist,i)-sdState(ist,i+1)-sdState(ist,i-1))/h2+vv(i)*sdState(ist,i))
LET sn = sn + sdState(ist,i)*sdState(ist,i)
next i
LET energyOfState = s/sn
END FUNCTION
EXTERNAL SUB GramSchmidt(stateMax) !--- Gram-Schmidt orthonormalization
DECLARE EXTERNAL FUNCTION innerProduct
DECLARE EXTERNAL SUB normalizeState
CALL normalizeState(0)
FOR istate=1 TO stateMax-1
FOR ist=0 TO istate-1
LET s = innerProduct(ist,istate)
FOR i=1 TO NNx-2
LET sdState(istate,i) = sdState(istate,i) - s*sdState(ist,i)
NEXT i
NEXT ist
CALL normalizeState(istate)
NEXT iState
END SUB
EXTERNAL FUNCTION innerProduct(ist,jst) !--- <ist|jst>
LET s = 0
FOR i=1 TO NNx-2
LET s = s + sdState(ist,i)*sdState(jst,i)
NEXT i
LET innerProduct = s*dx
END FUNCTION
EXTERNAL SUB normalizeState(ist)
LET s = 0
FOR i=1 TO NNx-2
LET s = s + sdState(ist,i)*sdState(ist,i)*dx
NEXT i
LET a = SQR(1/s)
FOR i=1 TO NNx-2
LET sdState(ist,i) = a*sdState(ist,i)
NEXT i
END SUB
! ---------- drawState
EXTERNAL SUB drawState !public
DECLARE EXTERNAL SUB dispInnerProduct
LET sc = 20
LET xp = 50
LET yp = 180
LET vMag = 10
LET stMag = 100
SET DRAW MODE HIDDEN
CLEAR
SET LINE COLOR 1 ! black : PLOT x-axis
PLOT LINES: xp,yp;dx*NNx*sc+xp,yp
!---plot V(x)
SET LINE COLOR 10 ! dark green : PLOT potential V(x);
FOR i=0 TO NNx-1
PLOT LINES: dx*i*sc+xp,vv(i)*vMag+yp;
NEXT i
PLOT LINES
SET TEXT HEIGHT 5
FOR ist=0 TO 9
IF sdEnergy(ist)<12 THEN
SET LINE COLOR 1+ist
SET TEXT COLOR 1+ist
FOR i=0 TO NNx-1 !plot wave function |psi(x,t)>
PLOT LINES: i*dx*sc+xp, sdState(ist,i)*stMag+sdEnergy(ist)*20+yp;
NEXT i
PLOT LINES
PLOT TEXT, AT xp-20,sdEnergy(ist)*20+yp :"|"&STR$(ist)&">"
END IF
NEXT ist
CALL dispInnerProduct(0,dx*NNx*sc+xp+20,yp)
!--- caption
SET TEXT HEIGHT 10
SET TEXT COLOR 1 ! black
PLOT TEXT, AT 50,100 ,USING "count =###### ":iterCount
PLOT TEXT, AT 50, 85 ,USING "E0 =###.########## E5 =###.##########(au)":sdEnergy(0),sdEnergy(5)
PLOT TEXT, AT 50, 70 ,USING "E1 =###.########## E6 =###.##########(au)":sdEnergy(1),sdEnergy(6)
PLOT TEXT, AT 50, 55 ,USING "E2 =###.########## E7 =###.##########(au)":sdEnergy(2),sdEnergy(7)
PLOT TEXT, AT 50, 40 ,USING "E3 =###.########## E8 =###.##########(au)":sdEnergy(3),sdEnergy(8)
PLOT TEXT, AT 50, 25 ,USING "E4 =###.########## E9 =###.##########(au)":sdEnergy(4),sdEnergy(9)
PLOT TEXT, AT 50, 10 :"steepest descent method 1D"
PLOT TEXT, AT 50,470 :"'.' exit "
PLOT TEXT, AT 50,455 :"'0' hermonics k*x^2 '1' quantum well"
SET DRAW MODE EXPLICIT
END SUB
EXTERNAL SUB dispInnerProduct(ist,xp,yp)
DECLARE EXTERNAL FUNCTION innerProduct
SET TEXT HEIGHT 5
SET TEXT COLOR 1 ! black
FOR jst=0 TO 9
PLOT TEXT, AT xp,yp+15*jst ,USING "("&STR$(ist)&"|"&STR$(jst)&") = -%.###^^^^":innerProduct(ist,jst)
NEXT jst
PLOT TEXT, AT xp,yp+15*10 :"(i|j) inner product"
END SUB
END MODULE
|
|