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