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