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