1!
2! Copyright (C) 2007 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 elsdps_paw( )
11  !---------------------------------------------------------------
12  !
13  !   total paw energy in the local-spin-density scheme
14  !
15  use kinds, only: DP
16  use constants, only: fpi
17  use radial_grids, only : ndmx
18  use ld1_parameters, only : nwfsx
19  use ld1inc, only : nlcc, grid, nspin, rhoc, lsd, &
20       encl, ehrt, ecxc, evxt, ekin, ecc, epseu,  &
21       nwfts, enlts, octs, paw_energy
22  use funct, only : dft_is_gradient
23  implicit none
24  real(DP) :: &
25       excc, vxcc(2), &   ! exch-corr energy from core charge
26       int_0_inf_dr,  &   ! the integral function
27       rh0(2),        &   ! the charge in a given point
28       rhc,           &   ! core charge in a given point
29       edcts              ! auxiliary energy
30
31  real(DP),allocatable :: &
32       vgc(:,:),   &   ! the gga potential
33       egc(:),     &   ! the gga energy
34       rho_aux(:,:), & ! auxiliary space
35       exccc(:)        ! the exchange and correlation energy of the core
36  REAL(dp) :: & ! compatibility with metaGGA - not yet used
37       tau(ndmx) = 0.0_dp, vtau(ndmx) = 0.0_dp
38
39  integer :: &
40       i,ns,ierr
41
42  !
43  !  If there is NLCC we calculate here also the exchange and correlation
44  !  energy of the pseudo core charge.
45  !  This quantity is printed but not added to the total energy
46  !
47  ecc=0.0_DP
48  if (nlcc) then
49     allocate(exccc(ndmx), stat=ierr)
50     exccc=0.0_DP
51     rh0(1)=0.0_DP
52     rh0(2)=0.0_DP
53     do i=1,grid%mesh
54        rhc= rhoc(i)/grid%r2(i)/fpi
55        call vxc_t(lsd,rh0,rhc,excc,vxcc)
56        exccc(i) = excc*rhoc(i)
57     enddo
58     if (dft_is_gradient()) then
59        allocate(rho_aux(ndmx,2), stat=ierr)
60        allocate(vgc(ndmx,2),stat=ierr)
61        allocate(egc(ndmx),stat=ierr)
62        vgc=0.0_DP
63        egc=0.0_DP
64        rho_aux=0.0_DP
65        call vxcgc ( ndmx, grid%mesh, nspin, grid%r, grid%r2, rho_aux, &
66             rhoc, vgc, egc, tau, vtau, 1)
67        do i=1,grid%mesh
68           exccc(i) = exccc(i) + egc(i)*fpi*grid%r2(i)
69        enddo
70        deallocate(egc)
71        deallocate(vgc)
72        deallocate(rho_aux)
73     endif
74     ecc=  int_0_inf_dr(exccc,grid,grid%mesh,2)
75     deallocate(exccc)
76  endif
77  !
78  !  Add the three contributions for each energy
79  !
80  encl= paw_energy(5,1)+paw_energy(5,2)-paw_energy(5,3)
81  ehrt= paw_energy(2,1)+paw_energy(2,2)-paw_energy(2,3)
82  ecxc= paw_energy(3,1)+paw_energy(3,2)-paw_energy(3,3)
83  edcts=  paw_energy(4,1)+paw_energy(4,2)-paw_energy(4,3)
84  !
85  !  The nonlocal pseudopotential energy is not computed.
86  !
87  epseu=0.0_DP
88  !
89  !  Now compute the kinetic energy
90  !
91  ekin = -encl-edcts
92  do ns=1,nwfts
93     if (octs(ns) > 0.0_DP) then
94        ekin=ekin+octs(ns)*enlts(ns)
95     endif
96  end do
97  return
98end subroutine elsdps_paw
99