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