1!
2! Copyright (C) 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 elsd_highv (nc)
11  !---------------------------------------------------------------
12  !
13  !   additional information on the atomic total energy. The valence
14  !   and core contributions to the different terms are calculated
15  !
16  use kinds, only : DP
17  use constants, only : e2
18  use radial_grids, only: ndmx, hartree
19  use ld1inc, only: grid, aeccharge, aevcharge, nwf, nspin, enl, oc, v0, &
20                    vxcts, excts, excggats, nlcc, enclc, enclv, ehrtvv, &
21                    ehrtcv, ehrtcc, ekinc, ekinv, ecxc, ae_fc_energy, &
22                    core_state, ekinc0, etot, frozen_core, iswitch
23  implicit none
24  integer, intent(in) :: nc
25  real(DP),allocatable :: f2vv(:), f2cv(:), f2vc(:), f2cc(:), f1c(:), f1v(:),  &
26                          f5c(:), f5v(:), vhval(:), vhcore(:), vnew(:,:)
27
28  real(DP) :: int_0_inf_dr, rhotot, fact
29  integer:: i,n,ierr
30!
31!  when iswitch=1 this routine cannot work because the code does not
32!  know which are the core and valence states.
33!
34  if (iswitch==1) return
35
36  allocate(f1c(grid%mesh),stat=ierr)
37  allocate(f1v(grid%mesh),stat=ierr)
38
39  allocate(f2vv(grid%mesh),stat=ierr)
40  allocate(f2cv(grid%mesh),stat=ierr)
41  allocate(f2vc(grid%mesh),stat=ierr)
42  allocate(f2cc(grid%mesh),stat=ierr)
43
44  allocate(vnew(ndmx,2),stat=ierr)
45  allocate(vhval(ndmx),stat=ierr)
46  allocate(vhcore(ndmx),stat=ierr)
47
48  allocate(f5c(grid%mesh),stat=ierr)
49  allocate(f5v(grid%mesh),stat=ierr)
50
51  call set_rc_rv()
52  call v_of_rho_at (aevcharge,aeccharge,vhval,vxcts,excts,excggats, &
53                   vnew,.true.,1)
54  call hartree(0,2,grid%mesh,grid,aeccharge,vhcore)
55  vhcore=e2*vhcore
56
57  fact=1.0_DP/nspin
58  do i=1,grid%mesh
59     rhotot=aevcharge(i,1)
60     if (nspin==2) rhotot=rhotot+aevcharge(i,2)
61!
62!   The integral for the energy due to the interaction with the nuclei
63!
64     f1c(i)= v0(i) * aeccharge(i)
65     f1v(i)= v0(i) * rhotot
66!
67!   The integrals for the Hartree energy
68!
69     f2vv(i)= vhval (i) * rhotot
70     f2cv(i)= vhcore(i) * rhotot
71     f2vc(i)= vhval(i) *  aeccharge(i)
72     f2cc(i)= vhcore(i) * aeccharge(i)
73!
74!   The integral to be subtracted to the sum of the eigenvalues to
75!   get the kinetic energy
76!
77     f5c(i) =-vxcts(i,1)*aeccharge(i)*fact-f1c(i)-f2vc(i)-f2cc(i)
78     if (nspin==2) f5c(i) =f5c(i)-vxcts(i,2)*aeccharge(i)*fact
79     f5v(i) =-vxcts(i,1)*aevcharge(i,1)-f1v(i)-f2cv(i)-f2vv(i)
80     if (nspin==2) f5v(i) =f5v(i)-vxcts(i,2)*aevcharge(i,2)
81  enddo
82!
83!  Now compute the integrals
84!
85  enclc=      int_0_inf_dr(f1c,grid,grid%mesh,1)
86  enclv=      int_0_inf_dr(f1v,grid,grid%mesh,1)
87  ehrtvv=0.5_DP*int_0_inf_dr(f2vv,grid,grid%mesh,2)
88  ehrtcc=0.5_DP*int_0_inf_dr(f2cc,grid,grid%mesh,2)
89  ehrtcv=     int_0_inf_dr(f2cv,grid,grid%mesh,2)
90
91  ekinc=      int_0_inf_dr(f5c,grid,grid%mesh,1)
92  ekinv=      int_0_inf_dr(f5v,grid,grid%mesh,1)
93  do n=1,nwf
94     if (oc(n)>0.0_DP.and.core_state(n)) ekinc=ekinc+oc(n)*enl(n)
95     if (oc(n)>0.0_DP.and..not.core_state(n)) ekinv=ekinv+oc(n)*enl(n)
96  enddo
97  if (nc==1) ekinc0=ekinc
98  if (frozen_core.and.nc>1) then
99     etot=etot-ekinc+ekinc0
100     ekinc=ekinc0
101  endif
102
103  ae_fc_energy=ekinv+ehrtvv+ehrtcv+ecxc+enclv
104
105  deallocate(f5c)
106  deallocate(f5v)
107  deallocate(f2vc)
108  deallocate(f2cv)
109  deallocate(f2cc)
110  deallocate(f2vv)
111  deallocate(f1c)
112  deallocate(f1v)
113  deallocate(vnew)
114  deallocate(vhval)
115  deallocate(vhcore)
116
117  return
118end subroutine elsd_highv
119