1! 2! Copyright (C) 2004 PWSCF group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8! 9!--------------------------------------------------------------- 10subroutine sic_correction(n,vhn1,vhn2,egc) 11 !--------------------------------------------------------------- 12 ! set up the orbital-dependent selfconsistent potential generated 13 ! by the n-th wavefunction - for self-interaction correction 14 ! 15 use kinds, only : dp 16 use radial_grids, only : ndmx 17 use constants, only: e2, fpi 18 use ld1inc, only : nspin, lsd, rel, nlcc, rhoc, grid, psi 19 use funct, only: dft_is_gradient 20 use radial_grids, only: hartree 21 implicit none 22 integer :: n 23 real(DP):: vhn1(ndmx),vhn2(ndmx), egc(ndmx) 24 REAL(dp) :: & ! compatibility with metaGGA - not yet used 25 tau(ndmx) = 0.0_dp, vtau(ndmx) = 0.0_dp 26 ! 27 integer :: i 28 real(DP):: rh(2), rhc, vxcp(2), excp 29 real(DP):: vgc(ndmx,2), egc0(ndmx), rhotot(ndmx,2) 30 logical :: gga 31 32 vhn1=0.0_dp 33 vhn2=0.0_dp 34 gga=dft_is_gradient() 35 nspin=1 36 if (lsd.eq.1) nspin=2 37 ! 38 ! compute hartree potential with the charge of orbital n 39 ! 40 rhotot=0.0_dp 41 if (rel.eq.2) then 42 do i=1,grid%mesh 43 rhotot(i,1)=psi(i,1,n)**2+psi(i,2,n)**2 44 enddo 45 else 46 do i=1,grid%mesh 47 rhotot(i,1)=psi(i,1,n)**2 48 enddo 49 endif 50 !call hartree(0,2*(ll(n)+1),grid%mesh,grid,rhotot,vhn1) 51 call hartree(0,2,grid%mesh,grid,rhotot,vhn1) 52 ! 53 ! add exchange and correlation potential: LDA or LSDA terms 54 ! 55 rhc=0.0_dp 56 rh=0.0_dp 57 do i=1,grid%mesh 58 vhn1(i) = e2*vhn1(i) 59 rh(1) = rhotot(i,1)/grid%r2(i)/fpi 60 if (nlcc) rhc = rhoc(i)/grid%r2(i)/fpi 61 call vxc_t(lsd,rh,rhc,excp,vxcp) 62 vhn2(i)= vhn1(i)+vxcp(1) 63 egc(i)= excp*rhotot(i,1) 64 end do 65 66 if (.not.gga) return 67 ! 68 ! add gradient-correction terms to exchange-correlation potential 69 ! 70 egc0=egc 71 call vxcgc ( ndmx, grid%mesh, nspin, grid%r, grid%r2, rhotot, rhoc, & 72 vgc, egc, tau, vtau, 1) 73 do i=1,grid%mesh 74 vhn2(i)=vhn2(i)+vgc(i,1) 75 egc(i)=egc(i)*grid%r2(i)*fpi+egc0(i) 76 enddo 77 return 78end subroutine sic_correction 79