1! 2! Copyright (C) 2004-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 11 !--------------------------------------------------------------- 12 ! 13 ! atomic total energy in the local-spin-density scheme 14 ! atomic pseudopotentials with nonlinear core correction are allowed 15 ! gradient correction allowed (A. Dal Corso fecit AD 1993) 16 ! 17 use kinds, only: DP 18 use constants, only: fpi 19 use radial_grids, only : ndmx 20 use ld1_parameters, only : nwfsx 21 use ld1inc, only : nlcc, grid, nspin, rhoc, rhos, lsd, vpsloc, vxt, vh, & 22 encl, ehrt, ecxc, evxt, ekin, ecc, epseu, vnl, & 23 etots, pseudotype, phits, ikk, nbeta, betas, bmat, & 24 nwfts, rel, jjts, llts, octs, enlts, jjs, lls, & 25 vxc, exc, excgga 26 use funct, only : dft_is_gradient 27 implicit none 28 real(DP) :: & 29 excc, vxcc(2), & ! exch-corr energy from core charge 30 int_0_inf_dr, & ! the integral function 31 rh0(2), & ! the charge in a given point 32 rhc, & ! core charge in a given point 33 rho_tot, & ! the total charge in one point 34 work(nwfsx) ! auxiliary space (similar to becp) 35 36 real(DP),allocatable :: & 37 f1(:), & ! auxiliary 38 f2(:), & ! auxiliary 39 f3(:), & ! auxiliary 40 f4(:), & ! auxiliary 41 f5(:), & ! auxiliary 42 vgc(:,:), & ! the gga potential 43 egc(:), & ! the gga energy 44 rho_aux(:,:), & ! auxiliary space 45 exccc(:) ! the exchange and correlation energy of the core 46 REAL(dp) :: & ! compatibility with metaGGA - not yet used 47 tau(ndmx) = 0.0_dp, vtau(ndmx) = 0.0_dp 48 49 integer :: & 50 n,i,ns,nst,lam,n1,n2,ikl,ierr,ind 51 52 allocate(f1(grid%mesh), stat=ierr) 53 allocate(f2(grid%mesh), stat=ierr) 54 allocate(f3(grid%mesh), stat=ierr) 55 allocate(f4(grid%mesh), stat=ierr) 56 allocate(f5(grid%mesh), stat=ierr) 57 allocate(exccc(ndmx), stat=ierr) 58 ! 59 ! If there is NLCC we calculate here also the exchange and correlation 60 ! energy of the pseudo core charge. 61 ! This quantity is printed but not added to the total energy 62 ! 63 exccc=0.0_DP 64 ecc=0.0_DP 65 if (nlcc) then 66 rh0(1)=0.0_DP 67 rh0(2)=0.0_DP 68 do i=1,grid%mesh 69 rhc= rhoc(i)/grid%r2(i)/fpi 70 call vxc_t(lsd,rh0,rhc,excc,vxcc) 71 exccc(i) = excc*rhoc(i) 72 enddo 73 if (dft_is_gradient()) then 74 allocate(rho_aux(ndmx,2), stat=ierr) 75 allocate(vgc(ndmx,2),stat=ierr) 76 allocate(egc(ndmx),stat=ierr) 77 vgc=0.0_DP 78 egc=0.0_DP 79 rho_aux=0.0_DP 80 call vxcgc ( ndmx, grid%mesh, nspin, grid%r, grid%r2, rho_aux, & 81 rhoc, vgc, egc, tau, vtau, 1) 82 do i=1,grid%mesh 83 exccc(i) = exccc(i) + egc(i)*fpi*grid%r2(i) 84 enddo 85 deallocate(egc) 86 deallocate(vgc) 87 deallocate(rho_aux) 88 endif 89 ecc= int_0_inf_dr(exccc,grid,grid%mesh,2) 90 endif 91 ! 92 ! Now prepare the integrals 93 ! 94 do i=1,grid%mesh 95 rho_tot=rhos(i,1) 96 if (lsd.eq.1) rho_tot=rho_tot+rhos(i,2) 97 ! 98 ! The integral for the interaction with the local potential 99 ! 100 f1(i)= vpsloc(i) * rho_tot 101 ! 102 ! The integral for the Hartree energy 103 ! 104 f2(i)= vh(i) * rho_tot 105 ! 106 ! The integral for the exchange and correlation energy 107 ! 108 f3(i)= exc(i) * (rho_tot+rhoc(i)) + excgga(i) 109 ! 110 ! The integral for the interaction with the external potential 111 ! 112 f4(i)= vxt(i)*rho_tot 113 ! 114 ! The integral to add to the sum of the eigenvalues to have the 115 ! kinetic energy. 116 ! 117 f5(i) =-vxc(i,1)*rhos(i,1)-f1(i)-f2(i)-f4(i) 118 if (nspin==2) f5(i)=f5(i)-vxc(i,2)*rhos(i,2) 119 enddo 120 ! 121 ! And now compute the integrals 122 ! 123 encl= int_0_inf_dr(f1,grid,grid%mesh,1) 124 ehrt=0.5_DP*int_0_inf_dr(f2,grid,grid%mesh,2) 125 ecxc= int_0_inf_dr(f3,grid,grid%mesh,2) 126 evxt= int_0_inf_dr(f4,grid,grid%mesh,2) 127 ! 128 ! Now compute the nonlocal pseudopotential energy. There are two cases: 129 ! The potential in semilocal form or in fully separable form 130 ! 131 epseu=0.0_DP 132 if (pseudotype == 1) then 133 ! 134 ! Semilocal form 135 ! 136 do ns=1,nwfts 137 if (octs(ns)>0.0_DP) then 138 if ( rel < 2 .or. llts(ns) == 0 .or. & 139 abs(jjts(ns)-llts(ns)+0.5_dp) < 0.001_dp) then 140 ind=1 141 else if ( rel == 2 .and. llts(ns) > 0 .and. & 142 abs(jjts(ns)-llts(ns)-0.5_dp) < 0.001_dp) then 143 ind=2 144 endif 145 f1=0.0_DP 146 lam=llts(ns) 147 do n=1, grid%mesh 148 f1(n) = f1(n) + phits(n,ns)**2 * octs(ns) 149 enddo 150 do n=1,grid%mesh 151 f1(n) = f1(n) * vnl(n,lam,ind) 152 end do 153 if (ikk(ns) > 0) & 154 epseu = epseu + int_0_inf_dr(f1,grid,ikk(ns),2*(lam+1)) 155 endif 156 enddo 157 else 158 ! 159 ! Fully separable form 160 ! 161 do ns=1,nwfts 162 if (octs(ns).gt.0.0_DP) then 163 do n1=1,nbeta 164 if ( llts(ns).eq.lls(n1).and. & 165 abs(jjts(ns)-jjs(n1)).lt.1.e-7_DP) then 166 nst=(llts(ns)+1)*2 167 ikl=ikk(n1) 168 do n=1,ikl 169 f1(n)=betas(n,n1)*phits(n,ns) 170 enddo 171 work(n1)=int_0_inf_dr(f1,grid,ikl,nst) 172 else 173 work(n1)=0.0_DP 174 endif 175 enddo 176 do n1=1,nbeta 177 do n2=1,nbeta 178 epseu=epseu & 179 + bmat(n1,n2)*work(n1)*work(n2)*octs(ns) 180 enddo 181 enddo 182 endif 183 enddo 184 endif 185 ! 186 ! Now compute the kinetic energy 187 ! 188 ekin = int_0_inf_dr(f5,grid,grid%mesh,2) - epseu 189 do ns=1,nwfts 190 if (octs(ns).gt.0.0_DP) then 191 ekin=ekin+octs(ns)*enlts(ns) 192 endif 193 end do 194 ! 195 ! And the total energy 196 ! 197 etots= ekin + encl + epseu + ehrt + ecxc + evxt 198 199 deallocate(f5) 200 deallocate(f4) 201 deallocate(f3) 202 deallocate(f2) 203 deallocate(f1) 204 deallocate(exccc) 205 206 return 207end subroutine elsdps 208