1!
2! Copyright (C) 2004-2010 Quantum ESPRESSO 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 elsd ( zed, grid, rho, vxt, vh, vxc, exc, excgga, nwf,&
11                  nspin, enl, oc, etot, ekin, encl, ehrt, ecxc, evxt )
12  !---------------------------------------------------------------
13  !
14  !   atomic total energy in the local-spin-density scheme
15  !   atomic pseudopotentials with nonlinear core correction are allowed
16  !   gradient correction allowed (A. Dal Corso fecit AD 1993)
17  !
18  use kinds, only : DP
19  use constants, only: fpi
20  use radial_grids, only: ndmx, radial_grid_type
21  use funct, only: get_iexch, dft_is_meta
22  use ld1inc, only: vx, noscf, tau, vtau
23  implicit none
24  integer, intent(in) :: nwf, nspin
25  type(radial_grid_type),intent(in)::grid
26  real(DP), intent(in)  :: zed
27  real(DP), intent(in)  :: enl(nwf), oc(nwf), rho(ndmx,2), &
28                           vxt(ndmx), vh(ndmx), vxc(ndmx,2), exc(ndmx), &
29                           excgga(ndmx)
30  real(DP), intent(out) :: etot,encl,ekin,ehrt,ecxc,evxt
31  real(DP),allocatable :: f1(:), f2(:), f3(:), f4(:), f5(:)
32  real(DP) :: int_0_inf_dr, rhotot
33  integer:: i,n,is,ierr
34  logical:: oep, meta, kli
35
36  if (noscf) return
37  oep = get_iexch().eq.4
38  kli = get_iexch().eq.10
39
40  meta = dft_is_meta()
41
42  allocate(f1(grid%mesh),stat=ierr)
43  allocate(f2(grid%mesh),stat=ierr)
44  allocate(f3(grid%mesh),stat=ierr)
45  allocate(f4(grid%mesh),stat=ierr)
46  allocate(f5(grid%mesh),stat=ierr)
47
48  do i=1,grid%mesh
49     rhotot = rho(i,1)
50     if (nspin==2) rhotot=rhotot+rho(i,2)
51!
52!   The integral for the energy due to the interaction with nuclei
53!
54     f1(i)=-2.0_DP*zed/grid%r(i) * rhotot
55!
56!   The integral for the Hartree energy
57!
58     f2(i)= vh (i) * rhotot
59!
60!   The integral for the exchange and correlation energy
61!
62     f3(i) = exc(i) * rhotot + excgga(i)
63!
64!   The integral for the interaction with an external potential
65!
66     f4(i)= vxt(i) * rhotot
67!
68!   The integral to be subtracted to the sum of the eigenvalues to
69!   get the kinetic energy
70!
71     f5(i) =-vxc(i,1)*rho(i,1)-f1(i)-f2(i)-f4(i)
72     if (nspin==2) f5(i) =f5(i)-vxc(i,2)*rho(i,2)
73
74     if (oep .or. kli ) then
75        do is = 1, nspin
76           f5(i) = f5(i) - vx(i,is)*rho(i,is)
77        end do
78     end if
79
80     if (meta) THEN
81        do is = 1, nspin
82           f5(i) = f5(i) - vtau(i)*tau(i,is)*fpi*grid%r2(i)
83        end do
84     end if
85  enddo
86!
87!  Now compute the integrals
88!
89  encl=       int_0_inf_dr(f1,grid,grid%mesh,1)
90  ehrt=0.5_DP*int_0_inf_dr(f2,grid,grid%mesh,2)
91  ecxc=       int_0_inf_dr(f3,grid,grid%mesh,2)
92  evxt=       int_0_inf_dr(f4,grid,grid%mesh,2)
93  !
94!
95!  The kinetic energy is the sum of the eigenvalues plus the f5 integral
96!
97  ekin = int_0_inf_dr(f5,grid,grid%mesh,1)
98  do n=1,nwf
99     if (oc(n)>0.0_DP) ekin=ekin+oc(n)*enl(n)
100  enddo
101
102  if (oep .or. kli) call add_exchange (ecxc)
103
104  etot= ekin + encl + ehrt + ecxc + evxt
105
106  deallocate(f5)
107  deallocate(f4)
108  deallocate(f3)
109  deallocate(f2)
110  deallocate(f1)
111
112end subroutine elsd
113