[016]１次元結晶の電子状態を求める

 投稿者：mike  投稿日：2017年 5月29日(月)09時31分3秒
  定常状態の１次元結晶の電子状態とエネルギーを求めるプログラム(016periodicPSD1D.bas)を公開します。
　定常状態の電子のエネルギーと電子状態を求める問題は、系のハミルトニアンをHとして　H |i> = e_i |i>
の最小固有値とその近くの低いエネルギー固有値と固有関数を求める問題になります。
　008harmonicsSD1D.basを改造して周期的境界条件を入れ、矩形の繰り返しポテンシャルのもとで
最小固有値とその近くの低いエネルギー固有値と固有関数を求めます。

表示の説明：
はじめに関数が上から降ってくるのは、近似を繰り返すごとにエネルギーが低下していく様子をみるため、
試行波動関数にエネルギーの下駄を履かせているためです。
周期構造があると、ほとんど同じエネルギーをもったいくつかの状態ができることがわかります。
これが結晶のバンド構造の起源です。

試験環境：
本プログラムは十進BASIC 6.6.3.0 / macOS 10.7.5, 十進BASIC Ver 7.8.0 / windows 10　でテストしました。


------------------


!
! ========= steepest descent method 1D ==========
!
! 016periodicPSD1D.bas
!   Copyright(C) 2017 Mitsuru Ikeuchi
!
! ver 0.0.1   2017.05.29  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB psd1d.setInitialCondition,psd1d.SDiteration,psd1d.drawState
DECLARE EXTERNAL FUNCTION INKEY$

LET stateMax = 10  !state 0,1,...,stateMax-1
LET vIndex = 1     !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$ !from decimal BASIC library
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION

! ---------- periodic 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 psd1d
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=0 TO NNx-1
         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 !--- free space
      FOR i=0 TO NNx-1
         LET vv(i) = 0
      NEXT i
   ELSEIF vIndex=1 THEN !--- well
      LET el = NNx*dx
      FOR i=0 TO NNx-1
         LET x = i*dx
         IF sin(2*PI*4*x/el)>0 THEN LET vv(i) = 10 ELSE LET vv(i) = 0
      NEXT i
   END IF
END SUB

! ---------- steepest descent iteration

EXTERNAL SUB SDiteration(stateMax, iterMax) !public
   DECLARE EXTERNAL FUNCTION steepestDescent
   DECLARE EXTERNAL SUB GramSchmidt,sortState
   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)
      CALL sortState(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=0 TO NNx-1                          !---  |wrk> = (H-E_ist)|ist>
      LET ip1 = MOD(i+1,NNx)
      LET im1 = mod(i-1+NNx,NNx)
      LET wrk(i) = (2*sdState(ist,i)-sdState(ist,ip1)-sdState(ist,im1))/h2+(vv(i)-ei)*sdState(ist,i)
   NEXT i
   FOR i=0 TO NNx-1                          !---  |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=0 TO NNx-1
      LET ip1 = MOD(i+1,NNx)
      LET im1 = mod(i-1+NNx,NNx)
      LET s = s+sdState(ist,i)*((2*sdState(ist,i)-sdState(ist,ip1)-sdState(ist,im1))/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=0 TO NNx-1
            LET sdState(istate,i) = sdState(istate,i) - s*sdState(ist,i)
         NEXT i
      NEXT ist
      CALL normalizeState(istate)
   NEXT iState
END SUB

EXTERNAL SUB sortState(stateMax)
   FOR ist=stateMax-2 TO 0 STEP -1
      IF sdEnergy(ist)>sdEnergy(ist+1)+0.00001 THEN
         FOR i=0 TO NNx-1
            LET w = sdState(ist,i)
            LET sdState(ist,i) = sdState(ist+1,i)
            LET sdState(ist+1,i) = w
         NEXT i
         LET w = sdEnergy(ist)
         LET sdEnergy(ist) = sdEnergy(ist+1)
         LET sdEnergy(ist+1) = w
      END IF
   NEXT ist
END SUB

EXTERNAL FUNCTION innerProduct(ist,jst)  !--- <ist|jst>
   LET s = 0
   FOR i=0 TO NNx-1
      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=0 TO NNx-1
      LET s = s + sdState(ist,i)*sdState(ist,i)*dx
   NEXT i
   LET a = SQR(1/s)
   FOR i=0 TO NNx-1
      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 :"periodic steepest descent method 1D"
   PLOT TEXT, AT 50,470 :"[.] exit "
   PLOT TEXT, AT 50,450 :"[0] free space   [1] 1d-crystal "
   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
 

戻る