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