[018]周期的境界条件下の電子状態（２次元）

 投稿者：mike  投稿日：2017年 6月26日(月)09時33分32秒
  　２次元-周期的境界条件下の下で定常状態の電子のエネルギーと電子状態を求めるプログラム
(018harmonicsPSD2D.bas)を公開します。これは、最急降下法による２次元の定常状態の電子状態を求めるプログラム
[011]に周期的境界条件を付けたものです。

表示の説明：
　２次元-周期的境界条件下の自由空間（V(x,y)=0）では基底状態は|0＞=constant, E(0)=0 になります。
x、y方向の周期が同じ場合、E1=E2=E3=E4で縮退しているため、|1＞..|4＞はこの固有空間の中の
任意の直交する関数の組でよいため、形が一意には決まりません。
E5=E6=E7=E8も縮退しています。

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


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

!
! ========= periodic steepest descent method 2D ==========
!
! 018periodicPSD2D.bas
!   Copyright(C) 2017 Mitsuru Ikeuchi
!
! ver 0.0.1   2017.06.26  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB psd2d.setInitialCondition,psd2d.SDiteration,psd2d.drawState
DECLARE EXTERNAL FUNCTION INKEY$

LET stateMax = 10  !state 0,1,...,stateMax-1
LET vIndex = 0     !0:free space,  1:2d-crystal
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,vIndex)
   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
LOOP UNTIL S$=CHR$(27)
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 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 psd2d
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=0 TO NNx-1
         FOR j=0 TO NNx-1
            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
         ELSEIF vIndex=1 THEN
            IF (x-x0)*(y-y0)>0 THEN LET vv(i,j)=10 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=0 TO NNx-1                          !---  |wrk> = (H-E_ist)|ist>
      LET ip1 = MOD(i+1,NNx)
      LET im1 = MOD(i-1+NNx,NNx)
      FOR j=0 TO NNy-1
         LET jp1 = MOD(j+1,NNy)
         LET jm1 = MOD(j-1+NNy,NNy)
         LET pij = sdState(ist,i,j)
         LET wrk(i,j) = (4*pij-sdState(ist,ip1,j)-sdState(ist,im1,j)-sdState(ist,i,jm1)-sdState(ist,i,jp1))/h2+(vv(i,j)-ei)*pij
      NEXT j
   NEXT i
   FOR i=0 TO NNx-1                          !---  |ist(next)> = |ist> - damp*|wrk>  ( norm(|ist(next)>) <>1 )
      FOR j=0 TO NNy-1
         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=0 TO NNx-1
      LET ip1 = MOD(i+1,NNx)
      LET im1 = MOD(i-1+NNx,NNx)
      FOR j=0 TO NNy-1
         LET jp1 = MOD(j+1,NNy)
         LET jm1 = MOD(j-1+NNy,NNy)
         LET pij = sdState(ist,i,j)
         LET pij = sdState(ist,i,j)
         LET s = s+pij*((4*pij-sdState(ist,ip1,j)-sdState(ist,im1,j)-sdState(ist,i,jm1)-sdState(ist,i,jp1))/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=0 TO NNx-1
            FOR j=0 TO NNy-1
               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=0 TO NNx-1
      FOR j=0 TO NNy-1
         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=0 TO NNx-1
      FOR j=0 TO NNy-1
         LET s = s + sdState(ist,i,j)*sdState(ist,i,j)*dx*dy
      NEXT j
   NEXT i
   LET a = SQR(1/s)
   FOR i=0 TO NNx-1
      FOR j=0 TO NNy-1
         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,vIndex)
   DECLARE EXTERNAL SUB setRotateXYParameters,drawState3D,draw3DPsix6,drawPsix6
   IF vIndex=0 THEN LET v$="free space" ELSE LET v$="2D-crystal"
   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 :"periodic steepest descent method 2D"
   PLOT TEXT, AT 50,480 :"potential -> [6] free space [7] 2d-crystal "
   PLOT TEXT, AT 50,465 :"display   -> [8] psi3D      [9] psi  "
   PLOT TEXT, AT 50,450 :"[esc] exit "
   PLOT TEXT, AT 50,420 :"potential : "&v$
   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
 

戻る