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

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

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

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


-----------

!
! ========= periodic steepest descent method 3D ==========
!
! 020periodicPSD3D.bas
!   Copyright(C) 2017 Mitsuru Ikeuchi
!
! ver 0.0.1   2017.07.10  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB psd3d.setInitialCondition,psd3d.SDiteration,psd3d.drawState
DECLARE EXTERNAL FUNCTION INKEY$

LET Ax = PI/12     !rotate angle around x-axis
LET Ay = -PI/6     !rotate angle around y-axis
LET ddAy = PI/180   !Ay=Ay+ddAy
LET vIndex = 3     !vIndex=0: V=1/r, 1: V=(1/2)r^2, 2:V=0(r<=4) or V=16(r>4) 3:V=0
LET stateMax = 6   !state 0,1,...,stateMax-1
LET drawMode = 0   !drawMode=0:density3D,  1:grid |i> in (x,y,0)
LET ist = 0        !display state |i>
CALL setInitialCondition(stateMax,vIndex)

DO
   LET Ay = Ay + ddAy
   CALL SDiteration(stateMax, 2) !2 - iteration in SDiteration()
   CALL drawState(ist,Ax,Ay,drawMode)
   LET S$=INKEY$
   IF S$="0" OR S$="1" OR S$="2" OR S$="3"  OR S$="4" OR S$="5" THEN
      LET ist = VAL(S$)
   ELSEIF S$="W" OR S$="w" THEN
      if ddAy=0.0 then LET ddAy=-pi/180 else LET ddAy = 0.0
   ELSEIF S$="E" OR S$="e" THEN
      IF ddAy=0.0 THEN LET ddAy= PI/180 ELSE LET ddAy = 0.0
   ELSEIF S$="D" OR S$="d" THEN
      LET drawMode = MOD(drawMode+1,2)
   END IF
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 3D ----------
!
! 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 psd3d
MODULE OPTION ARITHMETIC NATIVE
OPTION BASE 0
PUBLIC SUB setInitialCondition, SDiteration, drawState
SHARE NUMERIC NNx, NNy, NNz, dx, dy, dz, iterCount
SHARE NUMERIC cosAx,sinAx,cosAy,sinAy,cx0,cy0,cz0 !--- 3D graphics
SHARE NUMERIC sdEnergy(10)           ! electron state energy
SHARE NUMERIC sdState(10,65,65,65)   ! electron states 0...20 0:ground state
SHARE NUMERIC wrk(65,65,65)          ! state work space in steepestDescent
SHARE NUMERIC vv(65,65,65)           ! external potential
SHARE NUMERIC xApex(0 TO 7),yApex(0 TO 7),zApex(0 TO 7) ! boxApex x- y- z-coordinate
SHARE NUMERIC pxApex(0 TO 7),pyApex(0 TO 7),pzApex(0 TO 7) ! rotated boxApex x- y- z-coordinate
SHARE NUMERIC boxEdge(0 TO 11,0 TO 2)! boxEdge(i,j) i-th edge j=0,1:j-th apex, j=2:for marking
LET NNx = 16                         ! x-max number of sdState(,NNx,NNy,NNz)
LET NNy = 16                         ! y-max number of sdState(,NNx,NNy,NNz)
LET NNz = 16                         ! z-max number of sdState(,NNx,NNy,NNz)
LET dx = 0.5                         ! (au) x division
LET dy = 0.5                         ! (au) y division
LET dz = 0.5                         ! (au) y division
LET iterCount = 0                    ! sd iteration count

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

EXTERNAL SUB setInitialCondition(stateMax,vIndex)   !public
   DECLARE EXTERNAL SUB setInitialState,setPotential,setBox
   RANDOMIZE
   CALL setInitialState(stateMax)
   CALL setPotential(vIndex)
   CALL setBox
   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 NNy-1
            FOR k=0 TO NNz-1
               LET sdState(ist,i,j,k) = RND-0.5
            NEXT k
         NEXT j
      NEXT i
      CALL normalizeState(ist)
   NEXT ist
END SUB

EXTERNAL SUB setPotential(vIndex)
   FOR i=0 TO NNx-1
      FOR j=0 TO NNy-1
         FOR k=0 TO NNz-1
            LET x = i*dx
            LET y = j*dy
            LET z = k*dz
            IF vIndex=0 THEN
               LET v0 = 10.0
               IF SQR((x-8)*(x-8)+(z-8)*(z-8))<2 THEN LET V0 = 0.0
               LET vv(i,j,k) = v0
            ELSEIF vIndex=1 THEN
               LET vv(i,j,k) = MIN(0.5*r*r,18)
            ELSEIF vIndex=2 THEN
               IF r<=4 THEN LET vv(i,j,k) = 0 ELSE LET vv(i,j,k) = 16
            ELSE
               LET vv(i,j,k) = 0
            END IF
         NEXT k
      NEXT j
   NEXT i
END SUB

EXTERNAL SUB setBox
   IF boxEdge(11,1)<>7 THEN
      DATA 0,0,0, 1,0,0, 0,1,0, 1,1,0, 0,0,1, 1,0,1, 0,1,1, 1,1,1
      FOR i=0 TO 7
         READ x,y,z
         LET xApex(i) = x*NNx*dx
         LET yApex(i) = y*NNy*dy
         LET zApex(i) = z*NNz*dz
      NEXT i
      DATA 0,1,9, 0,2,9, 0,4,9, 1,3,9, 1,5,9, 2,3,9, 2,6,9, 3,7,9, 4,5,9, 4,6,9, 5,7,9, 6,7,9
      MAT READ boxEdge !(0 TO 11,0 TO 2)
   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.05    !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)
         FOR k=0 TO NNz-1
            LET kp1 = MOD(k+1,NNz)
            LET km1 = MOD(k-1+NNz,NNz)
            LET kesdState = (6*sdState(ist,i,j,k)-sdState(ist,ip1,j,k)-sdState(ist,im1,j,k)&
&-sdState(ist,i,jp1,k)-sdState(ist,i,jm1,k)-sdState(ist,i,j,kp1)-sdState(ist,i,j,km1))/h2
            LET wrk(i,j,k) = kesdState+(vv(i,j,k)-ei)*sdState(ist,i,j,k)
         NEXT k
      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
         FOR k=0 TO NNz-1
            LET sdState(ist,i,j,k) = sdState(ist,i,j,k)-damp*wrk(i,j,k)
         NEXT k
      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                          !---  |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)
         FOR k=0 TO NNz-1
            LET kp1 = MOD(k+1,NNz)
            LET km1 = MOD(k-1+NNz,NNz)
            LET kesdState = (6*sdState(ist,i,j,k)-sdState(ist,ip1,j,k)-sdState(ist,im1,j,k)&
&-sdState(ist,i,jp1,k)-sdState(ist,i,jm1,k)-sdState(ist,i,j,kp1)-sdState(ist,i,j,km1))/h2
            LET s = s + sdState(ist,i,j,k)*(kesdState+vv(i,j,k)*sdState(ist,i,j,k))
            LET sn = sn+sdState(ist,i,j,k)*sdState(ist,i,j,k)
         NEXT k
      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
               FOR k=0 TO NNz-1
                  LET sdState(istate,i,j,k) = sdState(istate,i,j,k) - s*sdState(ist,i,j,k)
               NEXT k
            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
               FOR k=0 TO NNz-1
                  LET w = sdState(ist,i,j,k)
                  LET sdState(ist,i,j,k) = sdState(ist+1,i,j,k)
                  LET sdState(ist+1,i,j,k) = w
               NEXT k
            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
         FOR k=0 TO NNz-1
            LET s = s + sdState(ist,i,j,k)*sdState(jst,i,j,k)
         NEXT k
      NEXT j
   NEXT i
   LET innerProduct = s*dx*dy*dz
END FUNCTION

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

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

EXTERNAL SUB setRotateXY(angleX,angleY)
   LET cosAx = COS(angleX)
   LET sinAx = SIN(angleX)
   LET cosAy = COS(angleY)
   LET sinAy = SIN(angleY)
   LET cx0 = 0.5*NNx*dx
   LET cy0 = 0.5*NNy*dy
   LET cz0 = 0.5*NNz*dz
END SUB

EXTERNAL SUB rotateXY !particles and box apex
   FOR i=0 TO 7
      LET pxApex(i) = cosAy*(xApex(i)-cx0)+sinAy*(sinAx*(yApex(i)-cy0)+cosAx*(zApex(i)-cz0)) + cx0
      LET pyApex(i) = cosAx*(yApex(i)-cy0)-sinAx*(zApex(i)-cz0) + cy0
      LET pzApex(i) =-sinAy*(xApex(i)-cx0)+cosAy*(sinAx*(yApex(i)-cy0)+cosAx*(zApex(i)-cz0)) + cz0
   NEXT i
END SUB

EXTERNAL SUB drawDisk3D(x,y,z,r,mag,xShift,yShift)
   LET x1 = cosAy*(x-cx0)+sinAy*(sinAx*(y-cy0)+cosAx*(z-cz0)) + cx0
   LET y1 = cosAx*(y-cy0)-sinAx*(z-cz0) + cy0
   LET z1 =-sinAy*(x-cx0)+cosAy*(sinAx*(y-cy0)+cosAx*(z-cz0)) + cz0
   DRAW disk WITH SCALE(r)*SHIFT(x1*mag+xShift,y1*mag+yShift)
END SUB

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

EXTERNAL SUB markFarEdge
!  seek far apex --> iMin
   LET zMin = pzApex(0)
   LET iMin = 0
   FOR i=1 TO 7
      IF zMin>pzApex(i) THEN
         LET zMin = pzApex(i)
         LET iMin = i
      END IF
   NEXT i
   !mark far edge
   FOR iEdge = 0 TO 11
      LET boxEdge(iEdge,2) = 0
      IF (boxEdge(iEdge,0)=iMin OR boxEdge(iEdge,1)=iMin) THEN LET boxEdge(iEdge,2) = 1
   NEXT iEdge
END SUB

EXTERNAL SUB plotNearEdge(mag,xp,yp,lineColor)
   DECLARE EXTERNAL SUB plotEdge
   FOR iEdge=0 TO 11
      IF boxEdge(iEdge,2)=0 THEN !far edge mark = 0
         CALL plotEdge(iEdge,mag,xp,yp,lineColor)
      END IF
   NEXT iEdge
END SUB

EXTERNAL SUB plotFarEdge(mag,xp,yp,lineColor)
   DECLARE EXTERNAL SUB plotEdge
   FOR iEdge=0 TO 11
      IF boxEdge(iEdge,2)=1 THEN !far edge mark = 1
         CALL plotEdge(iEdge,mag,xp,yp,lineColor)
      END IF
   NEXT iEdge
END SUB

EXTERNAL SUB plotEdge(iEdge,mag,xp,yp,lineColor)
   SET LINE COLOR lineColor
   FOR i=0 TO 1
      LET iApex = boxEdge(iEdge,i)
      PLOT LINES: pxApex(iApex)*mag+xp, pyApex(iApex)*mag+yp;
   NEXT i
   PLOT LINES
END SUB

! ---------- draw state

EXTERNAL SUB drawState(ist,xRotateAngle,yRotateAngle,drawMode)
   DECLARE EXTERNAL SUB drawDensity3D,setRotateXY,drawDensity3D,drawStateGrid
   SET DRAW MODE HIDDEN
   CLEAR
   IF drawMode=0 THEN
      CALL drawDensity3D(ist,xRotateAngle,yRotateAngle,200/(NNx*dx),140,150) !!(ist,Ax,Ay,sc,xp,yp)
   ELSEIF drawMode=1 THEN
      CALL drawStateGrid(ist,200/(NNx*dx),0.1,140,200) !(ist,sc,zMag,xp,yp)
   END IF
   !--- caption
   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT 50, 70 ,USING "iterarion count =##### ":iterCount
   PLOT TEXT, AT 50, 55 ,USING "E0 =###.######### E3 =###.#########(au)":sdEnergy(0),sdEnergy(3)
   PLOT TEXT, AT 50, 40 ,USING "E1 =###.######### E4 =###.#########(au)":sdEnergy(1),sdEnergy(4)
   PLOT TEXT, AT 50, 25 ,USING "E2 =###.######### E5 =###.#########(au)":sdEnergy(2),sdEnergy(5)
   PLOT TEXT, AT 50, 10 :"periodic steepest descent method 3D"
   PLOT TEXT, AT 50,470 :"[0],[1],...,[5] view state   [esc] exit "
   PLOT TEXT, AT 50,455 :"[D] change drawMode"
   PLOT TEXT, AT 50,440 :"[W] rotateLeft/stop [E] rotateRight/stop"
   SET DRAW MODE EXPLICIT
END SUB

EXTERNAL SUB drawDensity3D(ist,xRotateAngle,yRotateAngle,sc,xp,yp) !----- drawMode=0
   DECLARE EXTERNAL SUB setRotateXY,rotateXY,markFarEdge,plotFarEdge,drawDisk3D,plotNearEdge
   CALL setRotateXY(xRotateAngle,yRotateAngle)
   CALL rotateXY    !rotateXY BOX
   CALL markFarEdge ! boxEdge(iEdge,2)=1:far side edge or 0:near side edge
   CALL plotFarEdge(sc,xp,yp,8) !8:gray
   FOR ii=0 TO NNx-1
      LET i=ii
      IF (pzApex(1)-pzApex(0)<0) THEN LET i=NNx-ii-1
      FOR jj=0 TO NNy-1
         LET j=jj
         IF (pzApex(2)-pzApex(0)<0) THEN LET j=NNy-jj-1
         FOR kk=0 TO NNz-1
            LET k=kk
            IF (pzApex(4)-pzApex(0)<0) THEN LET k=NNz-kk-1
            LET psi = sdState(ist,i,j,k)
            LET psi2 = psi*psi
            IF psi2>0.001 THEN
               IF psi>0 THEN SET AREA COLOR 4 ELSE SET AREA COLOR 2
               CALL drawDisk3D(i*dx,j*dy,k*dz,ABS(psi)*10,sc,xp,yp)
            END IF
         NEXT kk
      NEXT jj
   NEXT ii
   !SET AREA COLOR 1
   !CALL drawDisk3D(0,0,0,5,sc,xp,yp)
   CALL plotNearEdge(sc,xp,yp,1) !1:black
   SET TEXT COLOR 1
   SET TEXT HEIGHT 10
   PLOT TEXT, AT 50,100:"|"&STR$(ist)&">"
   PLOT TEXT, AT 50, 85 ,USING "Ax =###.#(deg) Ay =###.#(deg)":MOD(xRotateAngle*180/PI,360),MOD(yRotateAngle*180/PI,360)
END SUB

EXTERNAL SUB drawStateGrid(ist,sc,zMag,xp,yp)  !--- drawMode=1
   DECLARE EXTERNAL SUB setRotateXY,drawGridXY
   CALL setRotateXY(-PI/3,-PI/6)
   CALL drawGridXY(ist,sc,zMag,xp,yp)
   SET TEXT COLOR 1
   SET TEXT HEIGHT 10
   PLOT TEXT, AT 50,100:"|"&STR$(ist)&"> in (x,y,0)"
END SUB

EXTERNAL SUB drawGridXY(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,NNz/2)*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*dx,j*dx,zMag*(psi + 2*vv(i,j,NNz/2)),sc,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,NNz/2)*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*dx,j*dx,zMag*(psi + 2*vv(i,j,NNz/2)),sc,xShift,yShift) !(x,y,z)
      NEXT j
      PLOT LINES
   NEXT i
END SUB
END MODULE
 

戻る