[022]多電子系の電子密度とエネルギー

 投稿者：mike  投稿日：2017年 7月28日(金)09時45分33秒
  定常状態の多電子系の電子密度とエネルギーを求めるプログラム(022sampleLDA1D.bas)を公開します。
　１電子の定常状態の電子のエネルギーと電子状態を求める問題は、[008]最急降下法で解くことができました。
多電子系の場合、厳密に電子状態を求めるのは、電子数が多くなるとほとんど不可能になります。
ここでは、１つの電子に注目し他の電子は平均場として扱う近似により多電子系の問題を解く、密度汎関数法により
電子密度とエネルギーを求めます。
　密度汎関数法（Density Functional Theory）では、電子密度はそれぞれの電子軌道から決まる電子の存在確率の
総計から求められ、１電子はSchroedingerの方程式に似たKohn-Sham方程式（ポテンシャルの中に多電子の寄与がある）
に従うと仮定します。実効ポテンシャルVeffへの多電子の寄与は電子密度の関数になります。
        Veff = Vext + VH + Vx + Vc
Vextは外部ポテンシャル、VHは静電ポテンシャル、Vxは交換ポテンシャル、Vcは相関ポテンシャルです。
　本プログラムでは電子密度を仮定し、Kohn-Sham方程式を最急降下法で解き、各電子軌道の総計から電子密度を求め、
これを繰り返すことで電子密度と電子軌道、エネルギーを求めます。

表示の説明：
黒の線は外部ポテンシャル、灰色は電子密度を表します。[D]を押すことで表示を変更できます。

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

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

!
! ========= RSDFT - local density approximation 1D ==========
!
! 022sampleLDA1D.bas
!   Copyright(C) 2017 Mitsuru Ikeuchi
!
! ver 0.0.1   2017.07.28  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB lda1d.setInitialCondition,lda1d.setNumberOfElectron,lda1d.iterateLDA,lda1d.drawState
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 = 10   ! 10 = iteration in iterateLDA
LET dispMode = 1   !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$ !from decimal BASIC library
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION

! ---------- RSDFT - local density approximation 1D ----------
!
! - 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 lda1d
MODULE OPTION ARITHMETIC NATIVE
OPTION BASE 0
PUBLIC SUB setInitialCondition, setNumberOfElectron, iterateLDA, drawState
SHARE NUMERIC NNx, dx, lylz, iterCount
SHARE NUMERIC numberOfElectron, numberOfElectronBounds, numberOfOrbit, errorDecisionOrbit
SHARE NUMERIC energyMem, iterationError, convergenceErrorMax, dampingFactor
SHARE NUMERIC mixing,broadening
SHARE NUMERIC sdEnergy(20)    ! electron state energy
SHARE NUMERIC sdState(20,400) ! electron states 0...20 0:ground state
SHARE NUMERIC occupation(20)  ! occupation of orbit
SHARE NUMERIC wrk(400)        ! state work space in steepestDescent
SHARE NUMERIC vv(400)         ! effective potential
SHARE NUMERIC vvext(400)      ! external potential
SHARE NUMERIC vvh(400)        ! Hartree potential
SHARE NUMERIC vvx(400)        ! exchange potential
SHARE NUMERIC vvc(400)        ! correlation potential
SHARE NUMERIC rho(400)        ! charge density
LET NNx = 64                  ! max number of sdState(,NNx,NNy)
LET dx = 0.25                 ! (au) x-division
LET lylz = 16.0*16.0          ! (au) v=dx*ly*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.03     ! steepest descent method
LET mixing = 0.5             ! charge mixing in setRho()
LET broadening = 0.001       ! (au) level broadening IN setOccupation

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

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

EXTERNAL SUB setNumberOfElectron(nElectron)   !public
   LET iterCount = 0
   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
         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 setExternalPotential(vIndex)
   LET x0 = 0.5*NNx*dx
   IF vIndex=0 THEN !--- hermonic
      FOR i=0 TO NNx-1
         LET x = i*dx
         LET vvext(i) = MIN(0.5*(x-x0)^2,24.5)
      NEXT i
   ELSEIF vIndex=1 THEN !--- well
      FOR i=0 TO NNx-1
         LET x = i*dx
         IF ABS(x-x0)<4 THEN LET vvext(i) = 0 ELSE LET vvext(i) = 18
      NEXT i
   END IF
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

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

EXTERNAL SUB iterateLDA(stateMax, iterMax) !public
   DECLARE EXTERNAL FUNCTION steepestDescent
   DECLARE EXTERNAL SUB setElectronDensity,setEffectivePotential,solveKohnSham,sortState,setOccupation
   LET errorDecisionOrbit = (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
      LET rho(i) = rho(i)*(1.0-mixing)
      FOR ie=0 TO numberOfOrbit-1
         IF occupation(ie)>0.0 THEN
            LET rho(i) = rho(i) + mixing*occupation(ie)*(sdState(ie,i)*sdState(ie,i))/lylz
         END IF
      NEXT ie
   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
      LET vv(i) = vvext(i)+vvh(i)+vvx(i)+vvc(i)
   NEXT i
END SUB

EXTERNAL SUB poisson(iterMax)
   LET h2 = 4.0*PI*dx*dx
   LET omegav2 = 0.5*1.8
   FOR iter=0 TO iterMax-1
      FOR i=1 TO NNx-2
         LET vvh(i) = vvh(i)+omegav2*(vvh(i+1)+vvh(i-1)-2.0*vvh(i) +h2*rho(i))
      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
      LET rh = rho(i)
      LET rh3 = rh^0.33333333
      LET vvx(i) = 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) = ec*(1.0+1.22838*sqrtrs+0.4445*rs)/(1.0+1.0529*sqrtrs+0.3334*rs)
      ELSE
         LET vvc(i) = -0.05837-0.0084*rs +(0.0311+0.00133*rs)*LOG(rs)
      END IF
   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>
      LET wrk(i) = (2*sdState(ist,i)-sdState(ist,i+1)-sdState(ist,i-1))/h2+(vv(i)-ei)*sdState(ist,i)
   NEXT i
   FOR i=1 TO NNx-2                          !---  |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=1 TO NNx-2
      LET s = s+sdState(ist,i)*((2*sdState(ist,i)-sdState(ist,i+1)-sdState(ist,i-1))/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=1 TO NNx-2
            LET sdState(istate,i) = sdState(istate,i) - s*sdState(ist,i)
         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
            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

!--- (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
      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=1 TO NNx-2
      LET s = s + sdState(ist,i)*sdState(ist,i)*dx
   NEXT i
   LET a = SQR(1/s)
   FOR i=1 TO NNx-2
      LET sdState(ist,i) = a*sdState(ist,i)
   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-1
      LET s = s + (-0.5*vvh(i)-0.25*vvx(i)+eeCorrelation(rho(i))-vvc(i))*rho(i)
   NEXT i
   LET s = s*dx
   LET totalEnergy = sei + s
END FUNCTION

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

EXTERNAL SUB drawState(dispMode)  !public
   DECLARE EXTERNAL SUB dispInnerProduct,plotfn
   DECLARE EXTERNAL FUNCTION totalEnergy
   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 rho
   SET AREA COLOR 8 ! gray : PLOT rho;
   FOR i=0 TO NNx-2
      PLOT AREA: dx*i*sc+xp,yp;dx*i*sc+xp,rho(i)*20000+yp;dx*(i+1)*sc+xp,rho(i+1)*20000+yp;dx*(i+1)*sc+xp,yp
   NEXT i
   CALL plotFN(vvext,sc,vMag,xp,yp,0,1) !black, plot external potential vvext(x)
   !
   IF dispmode=2 THEN !---plot Vext,Veff(),vvh(),vvxc()x10
      SET TEXT HEIGHT 6
      CALL plotFN(vv,sc,vMag,xp,yp,0,10)    !dark green, plot effective potential vv(x)
      CALL plotFN(vvh,sc,vMag,xp,yp,0,2)    !blue, plot Hartree potential vvh(x)
      CALL plotFN(vvx,sc,vMag*10,xp,yp,0,4) !red, plot exchange potential vvx(x)
      CALL plotFN(vvc,sc,vMag*10,xp,yp,0,7) !magenta, plot correlation potential vvc(x)
      SET TEXT COLOR 1
      PLOT TEXT, AT 50, yp+65 :"Vext()"
      SET TEXT COLOR 10
      PLOT TEXT, AT 50, yp+50 :"Veff()"
      SET TEXT COLOR 2
      PLOT TEXT, AT 50, yp+35 :"Vh()"
      SET TEXT COLOR 4
      PLOT TEXT, AT 50, yp+20 :"Vx() x 10"
      SET TEXT COLOR 7
      PLOT TEXT, AT 50, yp+5 :"Vc() x 10"
   END IF
   IF dispMode=1 THEN !--- plot Vext(),|i>
      SET TEXT HEIGHT 5
      FOR ist=0 TO 4
         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
   END IF
   CALL dispInnerProduct(0,dx*NNx*sc+xp+20,yp)
   !--- caption
   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT 50,115 ,USING "iteration count =###### error =#.##########":iterCount,iterationError
   PLOT TEXT, AT 50,100 ,USING "total energy =###.##########":totalEnergy
   PLOT TEXT, AT 50, 85 ,USING "E0 =###.########## Occ =#.#####":sdEnergy(0),occupation(0)
   PLOT TEXT, AT 50, 70 ,USING "E1 =###.########## Occ =#.#####":sdEnergy(1),occupation(1)
   PLOT TEXT, AT 50, 55 ,USING "E2 =###.########## Occ =#.#####":sdEnergy(2),occupation(2)
   PLOT TEXT, AT 50, 40 ,USING "E3 =###.########## Occ =#.#####":sdEnergy(3),occupation(3)
   PLOT TEXT, AT 50, 25 ,USING "E4 =###.########## Occ =#.#####":sdEnergy(4),occupation(4)
   PLOT TEXT, AT 50, 10 :"RS-DFT - Local Density Approximation 1D"
   PLOT TEXT, AT 50,470 :"[esc] exit            [D] chage dispMode"
   PLOT TEXT, AT 50,455 :"[1] ... [6] number of electron"
   PLOT TEXT, AT 50,440 :"[7] hermonics k*x^2   [8] quantum well"
   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 numberOfOrbit-1
      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

EXTERNAL SUB plotFN(a(),sc,mag,xp,yp,offset,col)
   SET LINE COLOR col
   FOR i=0 TO NNx-1
      PLOT LINES: dx*i*sc+xp,a(i)*mag+offset+yp;
   NEXT i
   PLOT LINES
END SUB

END MODULE
 

戻る