1! 2! Copyright (C) 2004-2009 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 vxc_t(lsd,rho,rhoc,exc,vxc) 11 !--------------------------------------------------------------- 12 ! 13 ! this function returns the XC potential and energy in LDA or 14 ! LSDA approximation 15 ! 16 use kinds, only : DP 17 use xc_lda_lsda, only: xc 18 implicit none 19 integer, intent(in) :: lsd ! 1 in the LSDA case, 0 otherwise 20 real(DP), intent(in) :: rho(2), rhoc ! the system density 21 real(DP), intent(out):: exc(1), vxc(2) 22 ! 23 integer, parameter :: length=1 24 real(DP), dimension(length) :: ex, ec , arho 25 REAL(DP), dimension(length,2) :: rhoaux, vx, vc 26 ! 27 real(DP), parameter :: e2=2.0_dp, eps=1.e-30_dp 28 ! 29 vxc(1) = 0.0_dp 30 exc = 0.0_dp 31 ! 32 if (lsd == 0) then 33 ! 34 ! LDA case 35 ! 36 rhoaux(1,1) = abs(rho(1) + rhoc) 37 if (rhoaux(1,1) > eps) then 38 ! 39 CALL xc( length, 1, 1, rhoaux, ex, ec, vx(:,1), vc(:,1) ) 40 ! 41 vxc(1) = e2 * ( vx(1,1) + vc(1,1) ) 42 exc = e2 * ( ex(1) + ec(1) ) 43 endif 44 else 45 ! 46 ! LSDA case 47 ! 48 vxc(2)=0.0_dp 49 ! 50 rhoaux(1,1) = rho(1) + rho(2) + rhoc 51 rhoaux(1,2) = rho(1) - rho(2) 52 ! 53 CALL xc( length, 2, 2, rhoaux, ex, ec, vx, vc ) 54 ! 55 vxc(1) = e2 * ( vx(1,1) + vc(1,1) ) 56 vxc(2) = e2 * ( vx(1,2) + vc(1,2) ) 57 exc = e2 * ( ex(1) + ec(1) ) 58 ! 59 endif 60 ! 61 return 62 ! 63end subroutine vxc_t 64! 65! 66!--------------------------------------------------------------- 67subroutine vxcgc( ndm, mesh, nspin, r, r2, rho, rhoc, vgc, egc, & 68 tau, vtau, iflag ) 69 !--------------------------------------------------------------- 70 ! 71 ! 72 ! This routine computes the exchange and correlation potential and 73 ! energy to be added to the local density, to have the first 74 ! gradient correction. 75 ! In input the density is rho(r) (multiplied by 4*pi*r2). 76 ! 77 ! The units of the potential are Ry. 78 ! 79 use kinds, only : DP 80 use constants, only : fpi, e2 81 use funct, only : dft_is_meta 82 use xc_gga, only : xc_gcx, change_threshold_gga 83 use metagga, only : tpsscxc 84 implicit none 85 integer, intent(in) :: ndm,mesh,nspin,iflag 86 real(DP), intent(in) :: r(mesh), r2(mesh), rho(ndm,2), rhoc(ndm) 87 real(DP), intent(out):: vgc(ndm,2), egc(ndm) 88 real(DP), intent(in) :: tau(ndm,2) 89 real(DP), intent(out):: vtau(mesh) 90 91 integer :: i, is, ierr 92 real(DP) :: sx, sc, v2c, v1x, v2x, v1c 93 ! 94 REAL(DP) :: grho_v(3,mesh,nspin) 95 REAL(DP), ALLOCATABLE, DIMENSION(:) :: sx_v, sc_v, v2c_ud 96 REAL(DP), ALLOCATABLE, DIMENSION(:,:) :: v1x_v, v2x_v, v1c_v, v2c_v 97 ! 98 real(DP) :: v1xup, v1xdw, v2xup, v2xdw, v1cup, v1cdw 99 real(DP) :: v3x, v3c, de_cc, dv1_cc,dv2_cc 100 real(DP) :: segno, arho 101 real(DP) :: rh(1), zeta(1), grh2(1), grho2(2) 102 real(DP),parameter :: eps=1.e-12_dp, small = 1.E-10_DP 103 104 real(DP), allocatable :: grho(:,:), h(:,:), dh(:), rhoaux(:,:) 105 ! 106 ! First compute the charge and the charge gradient, assumed 107 ! to have spherical symmetry. The gradient is the derivative of 108 ! the charge with respect to the modulus of r. 109 ! 110 allocate(rhoaux(mesh,nspin),stat=ierr) 111 allocate(grho(mesh,nspin),stat=ierr) 112 allocate(h(mesh,2),stat=ierr) 113 allocate(dh(mesh),stat=ierr) 114 115 egc=0.0_dp 116 vgc=0.0_dp 117 118 do is=1,nspin 119 do i=1, mesh 120 rhoaux(i,is)=(rho(i,is)+rhoc(i)/nspin)/fpi/r2(i) 121 enddo 122 call radial_gradient(rhoaux(1,is),grho(1,is),r,mesh,iflag) 123 enddo 124 ! 125 do is=1,nspin 126 do i=1, mesh 127 grho_v(:,i,is) = grho(i,is)/SQRT(3.d0) 128 enddo 129 enddo 130 ! 131 allocate( sx_v(mesh) , sc_v(mesh) ) 132 allocate( v1x_v(mesh,nspin), v2x_v(mesh,nspin) ) 133 allocate( v1c_v(mesh,nspin), v2c_v(mesh,nspin) ) 134 IF (nspin==2) allocate( v2c_ud(mesh) ) 135 ! 136 if (nspin.eq.1) then 137 ! 138 IF ( dft_is_meta () ) THEN 139 ! 140 ! meta-GGA case 141 ! 142 ! for core correction - not implemented 143 de_cc = 0.0_dp 144 dv1_cc= 0.0_dp 145 dv2_cc= 0.0_dp 146 ! 147 vtau(:) = 0.0_dp 148 ! 149 do i=1,mesh 150 arho=abs(rhoaux(i,1)) 151 segno=sign(1.0_dp,rhoaux(i,1)) 152 if (arho.gt.eps.and.abs(grho(i,1)).gt.eps) then 153! 154! currently there is a single meta-GGA implemented (tpss) 155! that calculates all needed terms (LDA, GGA, metaGGA) 156! 157 call tpsscxc ( arho, grho(i,1)**2, tau(i,1)+tau(i,2), & 158 sx, sc, v1x, v2x, v3x, v1c, v2c, v3c ) 159 ! 160 egc(i)=sx+sc+de_cc 161 vgc(i,1)= v1x+v1c + dv1_cc 162 h(i,1) =(v2x+v2c)*grho(i,1)*r2(i) 163 vtau(i) = v3x+v3c 164 else 165 vgc(i,1)=0.0_dp 166 egc(i)=0.0_dp 167 h(i,1)=0.0_dp 168 vtau(i)=0.0_dp 169 endif 170 end do 171 172 ELSE 173 ! 174 ! GGA case 175 ! 176 CALL change_threshold_gga( small, eps**2 ) 177 ! 178 CALL xc_gcx( mesh, nspin, rhoaux, grho_v, sx_v, sc_v, v1x_v, v2x_v, v1c_v, v2c_v ) 179 ! 180 egc(1:mesh) = sx_v + sc_v 181 vgc(1:mesh,1) = v1x_v(1:mesh,1) + v1c_v(1:mesh,1) 182 h(1:mesh,1) = ( v2x_v(1:mesh,1) + v2c_v(1:mesh,1) ) * grho(1:mesh,1)*r2(1:mesh) 183 ! 184 END IF 185 ! 186 ELSE 187 ! 188 ! this is the \sigma-GGA case 189 ! 190 CALL change_threshold_gga( small, small ) 191 ! 192 CALL xc_gcx( mesh, 2, rhoaux, grho_v, sx_v, sc_v, v1x_v, v2x_v, v1c_v, v2c_v, v2c_ud ) 193 ! 194 do i = 1, mesh 195 egc(i)=sx+sc 196 vgc(i,1) = v1x_v(i,1)+v1c_v(i,1) 197 vgc(i,2) = v1x_v(i,2)+v1c_v(i,2) 198 h(i,1) =((v2x_v(i,1)+v2c_v(i,1))*grho(i,1)+v2c_v(i,1)*grho(i,2))*r2(i) 199 h(i,2) =((v2x_v(i,2)+v2c_v(i,1))*grho(i,2)+v2c_v(i,1)*grho(i,1))*r2(i) 200 ! if (i.lt.4) write(6,'(f20.12,e20.12,2f20.12)') & 201 ! rho(i,1)*2.0_dp, grho(i,1)**2*4.0_dp, & 202 ! vgc(i,1), h(i,2) 203 enddo 204! 205 endif 206 ! 207 deallocate( sx_v , sc_v ) 208 deallocate( v1x_v , v2x_v ) 209 deallocate( v1c_v , v2c_v ) 210 IF (nspin==2) deallocate( v2c_ud ) 211 ! 212 ! We need the gradient of h to calculate the last part of the exchange 213 ! and correlation potential. 214 ! 215 do is=1,nspin 216 call radial_gradient(h(1,is),dh,r,mesh,iflag) 217 ! 218 ! Finally we compute the total exchange and correlation energy and 219 ! potential. We put the original values on the charge and multiply 220 ! by e^2 = two to have as output Ry units. 221 222 do i=1, mesh 223 vgc(i,is)=vgc(i,is)-dh(i)/r2(i) 224 vgc(i,is)=e2*vgc(i,is) 225 if (is.eq.1) egc(i)=e2*egc(i) 226 ! if (is.eq.1.and.i.lt.4) write(6,'(3f20.12)') & 227 ! vgc(i,1) 228 enddo 229 enddo 230 IF ( dft_is_meta() ) vtau(:) = e2*vtau(:) 231 232 deallocate(dh) 233 deallocate(h) 234 deallocate(grho) 235 deallocate(rhoaux) 236 237 return 238end subroutine vxcgc 239