[024]多電子系の電子密度とエネルギー（２次元）

 投稿者：mike  投稿日：2017年 9月 6日(水)12時20分34秒
  　定常状態の多電子系の電子密度とエネルギーを求める２次元のプログラム(024sampleLDA2D.bas)を公開します。
１次元定常状態の多電子系の電子密度とエネルギーを求める２次元のプログラム(022sampleLDA1D.bas)の２次元版です。

表示の説明：
オリーブ色は電子密度を表します。[D]を押すことで表示を変更できます。
[1] ～　[6] を押すと、系の中の電子の個数が変更できます。
表示のEnは電子軌道のエネルギー、Occは軌道の占有率です。

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


プログラムが500行を超えたため、２つに分割して投稿します。


便宜のため、下記のURLでgoogle driveの一部を公開設定しました。
公開期限はありませんので、いつでもダウンロードできると思います。
本プログラムは024sampleLDA2Dv1.basです。
今までに投稿したプログラム（[001]～[023]）もあります。


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

!
! ========= RSDFT - local density approximation 2D ==========
!
! 024sampleLDA2D.bas
!   Copyright(C) 2017 Mitsuru Ikeuchi
!
! ver 0.0.1   2017.09.06  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB lda2d.setInitialCondition,lda2d.iterateLDA,lda2d.drawState,lda2d.setNumberOfElectron
DECLARE EXTERNAL FUNCTION INKEY$

LET stateMax = 10  !state 0,1,...,stateMax-1
LET vIndex = 0     !0:harmonic potential,  1:quantum well
LET nElectron = 4
LET iterMax = 5    !5 = iteration in iterateLDA
LET dispMode = 0   !0:Vext+rho, 1:Vext+rho+orbit, 2:rho+Vext+Veff+Vh+Vx+Vc
CALL setInitialCondition(stateMax,vIndex,nElectron)

DO
   CALL iterateLDA(stateMax, iterMax)
   CALL drawState(dispMode)
   LET S$=INKEY$
   IF S$="D" OR S$="d" THEN
      LET dispMode = mod(dispMode+1,3)
   ELSEIF S$="1" OR S$="2" OR S$="3" OR S$="4" OR S$="5" OR S$="6" THEN
      LET nElectron = VAL(s$)
      CALL setNumberOfElectron(nElectron)
   ELSEIF S$="7" THEN
      LET vIndex = 0
      CALL setInitialCondition(stateMax,vIndex,nElectron)
   ELSEIF S$="8" THEN
      LET vIndex = 1
      CALL setInitialCondition(stateMax,vIndex,nElectron)
   END IF
LOOP UNTIL S$=chr$(27)
END

EXTERNAL FUNCTION INKEY$
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION

! ---------- RSDFT - local density approximation 2D ----------
!
! - real space density functional theory - local density approximation
! - solve Kohn-Sham equation - successive approximation
! - Vxc : LDA(local density approximation)
!         J. P. Perdew and A. Zunger; Phys. Rev., B23, 5048 (1981)
!
! procedure
! (1) given: trial |i>, occupation(i)
!
! (2) set electron density rho
!     rho <-- |i>, occupation[(i), mixing rho(iter-1)
!
! (3) set effective potential
!      Veff = Vext + Vh + Vx + Vc
!      Vh <-- rho (Poisson eq. ,SOR iteration)
!      Vx,Vc <-- rho (LDA:Perdew-Zunger)
!
! (4) solve Kohn-Sham equation (successive approximation)
!      |i> steepest descent method: |i(next)> = |i> - dump{H-E}|i>
!      E(i) <-- <i|H|i>
!      {|0>,..,|i>,..,|N>} orthogonallization : Gram-Schmidt
!
! (5) sort state
!      sort orbit by E(i)
!
! (6) set occupation
!      occupation(i) <-- E(i)
!
! (7) goto (2)
!
MODULE lda2d
MODULE OPTION ARITHMETIC NATIVE
OPTION BASE 0
PUBLIC SUB setInitialCondition, setNumberOfElectron, iterateLDA, drawState
SHARE NUMERIC NNx,NNy, dx,dy, lz, iterCount
SHARE NUMERIC numberOfElectron, numberOfElectronBounds, numberOfOrbit, errorDecisionOrbit
SHARE NUMERIC energyMem, iterationError, convergenceErrorMax, dampingFactor
SHARE NUMERIC mixing,broadening
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 occupation(20)      ! occupation of orbit
SHARE NUMERIC wrk(200,200)        ! state work space in steepestDescent
SHARE NUMERIC vv(200,200)         ! effective potential
SHARE NUMERIC vvext(200,200)      ! external potential
SHARE NUMERIC vvh(200,200)        ! Hartree potential
SHARE NUMERIC vvx(200,200)        ! exchange potential
SHARE NUMERIC vvc(200,200)        ! correlation potential
SHARE NUMERIC rho(200,200)        ! charge density
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                  !
LET dx = 0.25                 ! (au) x-division
LET dy = 0.25
LET lz = 12.0                 ! (au) v=dx*dy*lz
LET iterCount = 0             ! sd iteration count
LET numberOfElectron = 4      !
LET numberOfElectronBounds = 6! selection bounds OF numberOfElectron
LET numberOfOrbit = 10        !
LET energyMem = 0.0
LET iterationError = 1.0
LET convergenceErrorMax = 1.0e-5
LET dampingFactor = 0.01     ! steepest descent method
LET mixing = 0.5             ! charge mixing in setRho()
LET broadening = 0.01        ! (au) level broadening IN setOccupation

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

EXTERNAL SUB setInitialCondition(stateMax,vIndex,nElectron)   !public
   DECLARE EXTERNAL SUB setInitialState,setExternalPotential,initDraw
   LET iterCount = 0
   LET numberOfElectron = nElectron
   CALL setInitialState(stateMax)
   CALL setExternalPotential(vIndex)
   CALL initDraw
   ! set window
   SET WINDOW 0,500, 0,500
END SUB

EXTERNAL SUB setNumberOfElectron(nElectron)   !public
   LET numberOfElectron = nElectron
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 setExternalPotential(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 vvext(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 vvext(i,j)=20 ELSE LET vvext(i,j)=0
         ELSE
            LET vvext(i,j)=0
         END IF
      NEXT j
   NEXT i
END SUB

EXTERNAL SUB setInitialOccupation(nOrbit, nElectron)
   LET occ = 1.0*nElectron/nOrbit
   FOR iState=0 TO nOrbit
      LET occupation(iState) = occ
   NEXT iState
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*vvext(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

! ---------- iterate LDA

EXTERNAL SUB iterateLDA(stateMax, iterMax) !public
   DECLARE EXTERNAL FUNCTION steepestDescent
   DECLARE EXTERNAL SUB setElectronDensity,setEffectivePotential,solveKohnSham,sortState,setOccupation
   LET errorDecisionOrbit = int((numberOfElectron-1)/2)
   CALL setElectronDensity
   CALL setEffectivePotential
   CALL solveKohnSham(numberOfOrbit,iterMax)
   CALL sortState(numberOfOrbit)
   CALL setOccupation(numberOfOrbit,numberOfElectron)
   LET iterationError = sdEnergy(errorDecisionOrbit) - energyMem
   LET energyMem = sdEnergy(errorDecisionOrbit)
END SUB

!--- (2) set electron density rho <-- sdState(), occupation()
EXTERNAL SUB setElectronDensity
   FOR i=0 TO NNx-1
      FOR j=0 TO NNy-1
         LET rho(i,j) = rho(i,j)*(1.0-mixing)
         FOR ie=0 TO numberOfOrbit-1
            IF occupation(ie)>0.0 THEN
               LET rho(i,j) = rho(i,j) + mixing*occupation(ie)*(sdState(ie,i,j)*sdState(ie,i,j))/lz
            END IF
         NEXT ie
      NEXT j
   NEXT i
END SUB

!--- (3) set effective potential <-- electron density
EXTERNAL SUB setEffectivePotential
   DECLARE EXTERNAL SUB poisson,setVxc

   CALL poisson(20) !setVh
   CALL setVxc
   FOR i=0 TO NNx-1
      FOR j=0 TO NNy-1
         LET vv(i,j) = vvext(i,j)+vvh(i,j)+vvx(i,j)+vvc(i,j)
      NEXT j
   NEXT i
END SUB

EXTERNAL SUB poisson(iterMax)
   LET h2 = 4.0*PI*dx*dx
   LET omegav4 = 1.8/4.0
   FOR iter=0 TO iterMax-1
      FOR i=1 TO NNx-2
         FOR j=1 TO NNy-2
            LET vvh(i,j) = vvh(i,j)+omegav4*(vvh(i+1,j)+vvh(i-1,j)+vvh(i,j+1)+vvh(i,j-1)-4.0*vvh(i,j) +h2*rho(i,j))
         NEXT j
      NEXT i
   NEXT iter
END SUB

EXTERNAL SUB setVxc  !set exchage and correlation potential (Perdew and Zunger)
   LET c1 = -0.984745022
   FOR i=0 TO NNx-1
      FOR j=1 TO NNy-2
         LET rh = rho(i,j)
         LET rh3 = rh^0.33333333
         LET vvx(i,j) = c1*rh3
         LET rs = 0.6204/(rh3+1.0e-20)
         IF rs>=1.0 THEN
            LET sqrtrs = SQR(rs)
            LET ec = -0.1423/(1.0+1.0529*sqrtrs+0.3334*rs)
            LET vvc(i,j) = ec*(1.0+1.22838*sqrtrs+0.4445*rs)/(1.0+1.0529*sqrtrs+0.3334*rs)
         ELSE
            LET vvc(i,j) = -0.05837-0.0084*rs +(0.0311+0.00133*rs)*LOG(rs)
         END IF
      NEXT j
   NEXT i
END SUB

EXTERNAL FUNCTION eeCorrelation(rh)
   LET r = 0.6204/(rh^0.33333333+1.0e-20)
   IF r>=1.0 THEN
      LET ec = -0.1423/(1.0+1.0529*SQR(r)+0.3334*r)
   ELSE
      LET ec = -0.0480-0.0116*r+(0.0311+0.0020*r)*LOG(r)
   END IF
   LET eeCorrelation = ec
END FUNCTION

!--- (4) solve Kohn-Sham equation
EXTERNAL SUB solveKohnSham(stateMax, iterMax)
   DECLARE EXTERNAL FUNCTION steepestDescent
   DECLARE EXTERNAL SUB GramSchmidt,sortState
   FOR i=0 TO iterMax-1
      FOR ist=0 TO stateMax-1
         LET sdEnergy(ist) = steepestDescent(ist,dampingFactor)
      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>
      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.0*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

!--- (5) sort state
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

!--- (6) set occupation
EXTERNAL SUB setOccupation(stateMax, nElectron)
   DECLARE EXTERNAL FUNCTION trialOcc,FermiDirac
   LET eUpper = sdEnergy(stateMax-1)+1.0
   LET eLower = sdEnergy(0)-1.0
   FOR i=0 TO stateMax-1
      IF sdEnergy(i)>eUpper THEN LET eUpper = sdEnergy(i)
      IF sdEnergy(i)<eLower THEN LET eLower = sdEnergy(i)
   NEXT i

   DO WHILE (eUpper-eLower>1.0e-12)
      LET eFermi = (eUpper+eLower)/2.0
      LET ntrial = trialOcc(stateMax, eFermi)
      IF ntrial<nElectron THEN
         LET eLower = eFermi
      ELSE
         LET eUpper = eFermi
      END IF
   LOOP
   LET eFermi = (eUpper+eLower)/2.0
   FOR i=0 TO stateMax-1
      LET occupation(i) = 2.0*FermiDirac(sdEnergy(i), eFermi)
      IF (occupation(i)<0.0001) THEN LET occupation(i) = 0.0
      IF (2.0-occupation(i)<0.0001) THEN LET occupation(i) = 2.0
   NEXT i
END SUB

EXTERNAL FUNCTION trialOcc(stateMax, eFermi)
   DECLARE EXTERNAL FUNCTION FermiDirac
   LET s = 0.0
   FOR i=0 TO stateMax-1
      LET s = s + 2.0*FermiDirac(sdEnergy(i), eFermi)
   NEXT i
   LET trialOcc = s
END FUNCTION

EXTERNAL FUNCTION FermiDirac(ee, ef)
   LET et = broadening !level width
   LET a = (ee-ef)/et
   IF a>100 THEN LET ret = 0.0 ELSE LET ret = 1.0/(EXP(a)+1.0)
   LET FermiDirac = ret
END FUNCTION

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

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

EXTERNAL FUNCTION totalEnergy
   DECLARE EXTERNAL FUNCTION eeCorrelation
   LET sei = 0.0
   FOR i=0 TO numberOfOrbit-1
      LET sei = sei + occupation(i)*sdEnergy(i)
   NEXT i
   LET s = 0.0
   FOR i=1 TO NNx-2
      FOR j=1 TO NNy-2
         LET s = s + (-0.5*vvh(i,j)-0.25*vvx(i,j)+eeCorrelation(rho(i,j))-vvc(i,j))*rho(i,j)
      NEXT j
   NEXT i
   LET s = s*dx*dy
   LET totalEnergy = sei + s
END FUNCTION

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

戻る