[011]定常状態の電子状態を求めるプログラム（２次元）

 投稿者：mike  投稿日：2017年 3月24日(金)09時28分56秒
  　定常状態の電子のエネルギーと電子状態を求める２次元のプログラム(011harmonicsSD2D.bas)を公開します。
　定常状態の電子のエネルギーと電子状態を求める問題は、系のハミルトニアンをHとして　H |i> = e_i |i>
の最小固有値とその近くの低いエネルギー固有値と固有関数を求める問題になります。
　最急降下法（steepest descent method）によって、この最小固有値付近の固有値と固有関数を求めます。
適当に与えた試行関数（ここでは乱数NNx*NNyの配列）の組{|i>}から出発し、エネルギーei=<i|H|i>を求め、
エネルギーの減少する方向　(H-ei)|i>　から　次世代の関数|i(next)> = |i> - damp(H-ei)|i>を求め、
関数の組を直交化することを繰り返す（逐次近似する）ことで小さい方からN個の固有値と固有関数を求めます。
本プログラムでは、 放物型と井戸型のポテンシャルを選択できます。

表示の説明：
波動関数は黄色が正、水色が負の値になることを表します。緑色はポテンシャルで色が濃いほど値が高くなります。
ポテンシャルが放物型のとき、理論的には固有エネルギーは2次元の場合、E0=1, E1=E2=2, E3=E4=E5=3,.. となります。
|1>,|2>はほぼ同じ固有値になる（縮退している）ため、向きが不定です。また、|3>,|4>,|5>も縮退しているため、
一次従属な関数も固有関数になるので関数の形が試行ごとに異なった形になります。


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


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


!
! ========= steepest descent method 2D ==========
!
! 011harmonicsSD2D.bas
! Mitsuru Ikeuchi (C) Copyleft
!
! ver 0.0.1   2017.03.24  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB sd2d.setInitialCondition,sd2d.SDiteration,sd2d.drawState
DECLARE EXTERNAL FUNCTION INKEY$

LET stateMax = 10  !state 0,1,...,stateMax-1
LET vIndex = 0     !0:harmonic,  1:well
LET iterMax = 5    !5 -iteration in SDiteration()
LET drawMode = 1   !0:draw3DPsi, 1:drawPsi
CALL setInitialCondition(stateMax,vIndex)

DO
   CALL SDiteration(stateMax, iterMax)
   CALL drawState(drawMode)
   LET S$=INKEY$
   IF S$="6" THEN CALL setInitialCondition(stateMax,0) !Hermonics
   IF S$="7" THEN CALL setInitialCondition(stateMax,1) !Well
   IF S$="8" THEN LET drawMode = 0 !0:draw3DPsi
   IF S$="9" THEN LET drawMode = 1 !1:drawPsi
   IF S$="." THEN EXIT DO
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 2D ----------
!
! system Hamiltonian: H = -delta/2 + V(r) , delta r = div grad r
! eigen energy 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 sd2d
MODULE OPTION ARITHMETIC NATIVE
OPTION BASE 0
PUBLIC SUB setInitialCondition, SDiteration, drawState
SHARE NUMERIC NNx, NNy, dx, dy, iterCount
SHARE NUMERIC cosAx,sinAx,cosAy,sinAy,cx0,cy0,cz0 !--- 3D graphics
SHARE NUMERIC sdEnergy(20)        ! electron state energy
SHARE NUMERIC sdState(20,200,200) ! electron states 0...20 0:ground state
SHARE NUMERIC wrk(200,200)        ! state work space in steepestDescent
SHARE NUMERIC vv(200,200)         ! external potential
SHARE NUMERIC psicol(64,64)       ! MAT PLOT CELL matrix for psi(x,y)
SHARE NUMERIC vcol(64,64)         ! MAT PLOT CELL matrix  for potential V(x,y)
LET NNx = 64                      ! max number of sdState(,NNx,NNy)
LET NNy = 64                      ! max number of sdState(,NNx,NNy)
LET dx = 0.25                     ! (au) x division
LET dy = 0.25                     ! (au) y division
LET iterCount = 0                 ! sd iteration count

! ---------- set initial condition

EXTERNAL SUB setInitialCondition(stateMax,vIndex)   !public
   DECLARE EXTERNAL SUB setInitialState,setPotential,initDraw
   LET iterCount = 0
   CALL setInitialState(stateMax)
   CALL setPotential(vIndex)
   CALL initDraw
   ! set window
   LET xMargin = 60
   LET yMargin = 120
   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
         FOR j=1 TO NNx-2
            LET sdState(ist,i,j) = RND-0.5
         NEXT j
      NEXT i
      CALL normalizeState(ist)
   NEXT ist
END SUB

EXTERNAL SUB setPotential(vIndex)
   LET x0 = 0.5*NNx*dx
   LET y0 = 0.5*NNy*dy
   FOR i=0 TO NNx-1
      FOR j=0 TO NNy-1
         LET x = i*dx
         LET y = j*dy
         IF vIndex=0 THEN
            LET vv(i,j) = 0.5*((x-x0)*(x-x0)+(y-y0)*(y-y0))
         ELSEIF vIndex=1 THEN
            LET r = SQR((x-x0)*(x-x0)+(y-y0)*(y-y0))
            IF r>5 THEN LET vv(i,j)=20 ELSE LET vv(i,j)=0
         ELSE
            LET vv(i,j)=0
         END IF
      NEXT j
   NEXT i
END SUB

EXTERNAL SUB initDraw !--- set set color pallet and MAT PLOT matrix vcol(,)
   FOR i = 0 TO 50 !--- set color pallet
      SET COLOR MIX( 40+i) 0.02*i,0.02*i,0     ! red for |psi> +
      SET COLOR MIX(100+i) 0,0.02*i,0.02*i     ! blue for |psi> -
      SET COLOR MIX(160+i) 0,0.02*i,0     ! green for V(x)
   NEXT i
   FOR i=0 TO NNx-1 !--- set vcol(,)
      FOR j=0 TO NNy-1
         LET col = 0.02*vv(i,j)
         IF col>1 THEN LET col = 1
         IF col<0 THEN LET col = 0
         LET vcol(i,j) = 160+INT(col*50)
      NEXT j
   NEXT i
END SUB

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

EXTERNAL SUB SDiteration(stateMax, iterMax) !public
   DECLARE EXTERNAL FUNCTION steepestDescent
   DECLARE EXTERNAL SUB GramSchmidt,sortState
   LET damp = 0.01 !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=1 TO NNx-2                          !---  |wrk> = (H-E_ist)|ist>
      FOR j=1 TO NNy-2
         LET pij = sdState(ist,i,j)
         LET wrk(i,j) = (4*pij-sdState(ist,i+1,j)-sdState(ist,i-1,j)-sdState(ist,i,j-1)-sdState(ist,i,j+1))/h2+(vv(i,j)-ei)*pij
      NEXT j
   NEXT i
   FOR i=1 TO NNx-2                          !---  |ist(next)> = |ist> - damp*|wrk>  ( norm(|ist(next)>) <>1 )
      FOR j=1 TO NNy-2
         LET sdState(ist,i,j) = sdState(ist,i,j)-damp*wrk(i,j)
      NEXT j
   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
      FOR j=1 TO NNy-2
         LET pij = sdState(ist,i,j)
         LET s = s+pij*((4*pij-sdState(ist,i+1,j)-sdState(ist,i-1,j)-sdState(ist,i,j-1)-sdState(ist,i,j+1))/h2+vv(i,j)*pij)
         LET sn = sn+pij*pij
      NEXT j
   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
            FOR j=1 TO NNy-2
               LET sdState(istate,i,j) = sdState(istate,i,j) - s*sdState(ist,i,j)
            NEXT j
         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
            FOR j=0 TO NNy-1
               LET w = sdState(ist,i,j)
               LET sdState(ist,i,j) = sdState(ist+1,i,j)
               LET sdState(ist+1,i,j) = w
            NEXT j
         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=1 TO NNx-2
      FOR j=1 TO NNy-2
         LET s = s + sdState(ist,i,j)*sdState(jst,i,j)
      NEXT j
   NEXT i
   LET innerProduct = s*dx*dy
END FUNCTION

EXTERNAL SUB normalizeState(ist)
   LET s = 0
   FOR i=1 TO NNx-2
      FOR j=1 TO NNy-2
         LET s = s + sdState(ist,i,j)*sdState(ist,i,j)*dx*dy
      NEXT j
   NEXT i
   LET a = SQR(1/s)
   FOR i=1 TO NNx-2
      FOR j=1 TO NNy-2
         LET sdState(ist,i,j) = a*sdState(ist,i,j)
      NEXT j
   NEXT i
END SUB

! ---------- utility

!EXTERNAL FUNCTION iterationCount
!   LET iterationCount = iterCount
!END FUNCTION

EXTERNAL FUNCTION stateEnergy(ist)
   LET stateEnergy = sdEnergy(ist)
END FUNCTION

! ---------- 3D graphics

EXTERNAL SUB setRotateXYParameters(angleX,angleY,xCenter,yCenter,zCenter)
   LET cosAx = COS(angleX)
   LET sinAx = SIN(angleX)
   LET cosAy = COS(angleY)
   LET sinAy = SIN(angleY)
   LET cx0 = xCenter
   LET cy0 = yCenter
   LET cz0 = zCenter
END SUB

EXTERNAL SUB plotLines3D(x,y,z,xShift,yShift) !shift*xRotateAx*yRotateAy*(shift^-1)
   LET x1 = cosAy*(x-cx0)+sinAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
   LET y1 = cosAx*(y-cy0)-sinAx*(z-cz0)
   LET z1 =-sinAy*(x-cx0)+cosAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
   PLOT LINES: x1+cx0+xShift, y1+cy0+yShift; !z=z1+cz0
END SUB

! ---------- drawState

EXTERNAL SUB drawState(drawMode)
   DECLARE EXTERNAL SUB setRotateXYParameters,drawState3D,draw3DPsix6,drawPsix6
   SET DRAW MODE HIDDEN
   CLEAR
   IF drawMode=0 THEN
      CALL setRotateXYParameters(-PI/6,-PI/12,NNx/2,NNy/2,0)
      CALL draw3DPsix6
   ELSEIF drawMode=1 THEN
      CALL drawPsix6
   END IF
   !--- caption
   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT 50,100 ,USING "iteration 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 2D"
   PLOT TEXT, AT 50,480 :"potential -> [6] Hermonics  [7] well "
   PLOT TEXT, AT 50,465 :"display   -> [8] psi3D      [9] psi  "
   PLOT TEXT, AT 50,450 :"[.] exit "
   SET DRAW MODE EXPLICIT
END SUB

EXTERNAL SUB draw3DPsix6
   DECLARE EXTERNAL SUB drawState3D
   CALL drawState3D(0,2,0.5, 30,120) !(ist,sc,zMag,xShift,yShift)
   CALL drawState3D(1,2,0.5,180,120)
   CALL drawState3D(2,2,0.5,330,120)
   CALL drawState3D(3,2,0.5, 30,270)
   CALL drawState3D(4,2,0.5,180,270)
   CALL drawState3D(5,2,0.5,330,270)
   SET TEXT HEIGHT 10
   PLOT TEXT, AT  30,130:"|0>"
   PLOT TEXT, AT 180,130:"|1>"
   PLOT TEXT, AT 330,130:"|2>"
   PLOT TEXT, AT  30,280:"|3>"
   PLOT TEXT, AT 180,280:"|4>"
   PLOT TEXT, AT 330,280:"|5>"
END SUB

EXTERNAL SUB drawState3D(ist,sc,zMag,xShift,yShift)
   DECLARE EXTERNAL SUB plotLines3D
   FOR j=0 TO NNy-1 STEP 1
      FOR i=0 TO NNx-1
         LET psi = sdState(ist,i,j)*200
         IF psi>1 THEN
            SET LINE COLOR 4 ! red
         ELSEIF psi<-1 THEN
            SET LINE COLOR 2 ! blue
         ELSE
            SET LINE COLOR 3 ! potential:green
         END IF
         CALL plotLines3D(i*sc,j*sc,zMag*(psi + 2*vv(i,j)),xShift,yShift) !(x,y,z)
      NEXT i
      PLOT LINES
   NEXT j
   FOR i=0 TO NNx-1 STEP 1
      FOR j=0 TO NNy-1
         LET psi = sdState(ist,i,j)*200
         IF psi>1 THEN
            SET LINE COLOR 4 ! red
         ELSEIF psi<-1 THEN
            SET LINE COLOR 2 ! blue
         ELSE
            SET LINE COLOR 3 ! potential:green
         END IF
         CALL plotLines3D(i*sc,j*sc,zMag*(psi + 2*vv(i,j)),xShift,yShift) !(x,y,z)
      NEXT j
      PLOT LINES
   NEXT i
END SUB

EXTERNAL SUB drawPsix6
   DECLARE EXTERNAL SUB drawPsi
   CALL drawPsi(0, 30,135) !(ist,xPos,yPos)
   CALL drawPsi(1,180,135)
   CALL drawPsi(2,330,135)
   CALL drawPsi(3, 30,285)
   CALL drawPsi(4,180,285)
   CALL drawPsi(5,330,285)
   SET TEXT HEIGHT 10
   PLOT TEXT, AT  30,120:"|0>"
   PLOT TEXT, AT 180,120:"|1>"
   PLOT TEXT, AT 330,120:"|2>"
   PLOT TEXT, AT  30,270:"|3>"
   PLOT TEXT, AT 180,270:"|4>"
   PLOT TEXT, AT 330,270:"|5>"
END SUB

EXTERNAL SUB drawPsi(ist,xPos,yPos)  !drawMode=0
   MAT psicol = vcol !--- psicol(,) <- Vcol(,) setted SUB initDraw
   FOR i=0 TO NNx-1 !--- set psicol(,) for MAT PLOT CELLS
      FOR j=0 TO NNy-1
         LET psi = sdState(ist,i,j)
         IF abs(psi)>0.002 THEN
            IF psi>=0 THEN
               LET col = psi*10
               IF col>1 THEN LET col = 1
               IF col<0 THEN LET col = 0
               LET psicol(i,j) = 40+INT(col*50)
            ELSE
               LET col = -psi*10
               IF col>1 THEN LET col = 1
               IF col<0 THEN LET col = 0
               LET psicol(i,j) = 100+INT(col*50)
            END IF
         END IF
      NEXT j
   NEXT i
   MAT PLOT CELLS,IN xPos,yPos; xPos+128,yPos+128 :psicol
END SUB

END MODULE
 

戻る