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