1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
|
!
! Copyright (C) 2004-2007 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!---------------------------------------------------------------
subroutine elsdps
!---------------------------------------------------------------
!
! atomic total energy in the local-spin-density scheme
! atomic pseudopotentials with nonlinear core correction are allowed
! gradient correction allowed (A. Dal Corso fecit AD 1993)
!
use kinds, only: DP
use constants, only: fpi
use radial_grids, only : ndmx
use ld1_parameters, only : nwfsx
use ld1inc, only : nlcc, grid, nspin, rhoc, rhos, lsd, vpsloc, vxt, vh, &
encl, ehrt, ecxc, evxt, ekin, ecc, epseu, vnl, &
etots, pseudotype, phits, ikk, nbeta, betas, bmat, &
nwfts, rel, jjts, llts, octs, enlts, jjs, lls, &
vxc, exc, excgga
use funct, only : dft_is_gradient, exc_t
implicit none
real(DP) :: &
int_0_inf_dr, & ! the integral function
rh0(2), & ! the charge in a given point
rhc, & ! core charge in a given point
rho_tot, & ! the total charge in one point
work(nwfsx) ! auxiliary space (similar to becp)
real(DP),allocatable :: &
f1(:), & ! auxiliary
f2(:), & ! auxiliary
f3(:), & ! auxiliary
f4(:), & ! auxiliary
f5(:), & ! auxiliary
vgc(:,:), & ! the gga potential
egc(:), & ! the gga energy
rho_aux(:,:), & ! auxiliary space
exccc(:) ! the exchange and correlation energy of the core
REAL(dp) :: & ! compatibility with metaGGA - not yet used
tau(ndmx) = 0.0_dp, vtau(ndmx) = 0.0_dp
integer :: &
n,i,ns,nst,lam,n1,n2,ikl,ierr,ind
allocate(f1(grid%mesh), stat=ierr)
allocate(f2(grid%mesh), stat=ierr)
allocate(f3(grid%mesh), stat=ierr)
allocate(f4(grid%mesh), stat=ierr)
allocate(f5(grid%mesh), stat=ierr)
allocate(exccc(ndmx), stat=ierr)
!
! If there is NLCC we calculate here also the exchange and correlation
! energy of the pseudo core charge.
! This quantity is printed but not added to the total energy
!
exccc=0.0_DP
ecc=0.0_DP
if (nlcc) then
rh0(1)=0.0_DP
rh0(2)=0.0_DP
do i=1,grid%mesh
rhc= rhoc(i)/grid%r2(i)/fpi
exccc(i) = exc_t(rh0,rhc,lsd)*rhoc(i)
enddo
if (dft_is_gradient()) then
allocate(rho_aux(ndmx,2), stat=ierr)
allocate(vgc(ndmx,2),stat=ierr)
allocate(egc(ndmx),stat=ierr)
vgc=0.0_DP
egc=0.0_DP
rho_aux=0.0_DP
call vxcgc ( ndmx, grid%mesh, nspin, grid%r, grid%r2, rho_aux, &
rhoc, vgc, egc, tau, vtau, 1)
do i=1,grid%mesh
exccc(i) = exccc(i) + egc(i)*fpi*grid%r2(i)
enddo
deallocate(egc)
deallocate(vgc)
deallocate(rho_aux)
endif
ecc= int_0_inf_dr(exccc,grid,grid%mesh,2)
endif
!
! Now prepare the integrals
!
do i=1,grid%mesh
rho_tot=rhos(i,1)
if (lsd.eq.1) rho_tot=rho_tot+rhos(i,2)
!
! The integral for the interaction with the local potential
!
f1(i)= vpsloc(i) * rho_tot
!
! The integral for the Hartree energy
!
f2(i)= vh(i) * rho_tot
!
! The integral for the exchange and correlation energy
!
f3(i)= exc(i) * (rho_tot+rhoc(i)) + excgga(i)
!
! The integral for the interaction with the external potential
!
f4(i)= vxt(i)*rho_tot
!
! The integral to add to the sum of the eigenvalues to have the
! kinetic energy.
!
f5(i) =-vxc(i,1)*rhos(i,1)-f1(i)-f2(i)-f4(i)
if (nspin==2) f5(i)=f5(i)-vxc(i,2)*rhos(i,2)
enddo
!
! And now compute the integrals
!
encl= int_0_inf_dr(f1,grid,grid%mesh,1)
ehrt=0.5_DP*int_0_inf_dr(f2,grid,grid%mesh,2)
ecxc= int_0_inf_dr(f3,grid,grid%mesh,2)
evxt= int_0_inf_dr(f4,grid,grid%mesh,2)
!
! Now compute the nonlocal pseudopotential energy. There are two cases:
! The potential in semilocal form or in fully separable form
!
epseu=0.0_DP
if (pseudotype == 1) then
!
! Semilocal form
!
do ns=1,nwfts
if (octs(ns)>0.0_DP) then
if ( rel < 2 .or. llts(ns) == 0 .or. &
abs(jjts(ns)-llts(ns)+0.5_dp) < 0.001_dp) then
ind=1
else if ( rel == 2 .and. llts(ns) > 0 .and. &
abs(jjts(ns)-llts(ns)-0.5_dp) < 0.001_dp) then
ind=2
endif
f1=0.0_DP
lam=llts(ns)
do n=1, grid%mesh
f1(n) = f1(n) + phits(n,ns)**2 * octs(ns)
enddo
do n=1,grid%mesh
f1(n) = f1(n) * vnl(n,lam,ind)
end do
if (ikk(ns) > 0) &
epseu = epseu + int_0_inf_dr(f1,grid,ikk(ns),2*(lam+1))
endif
enddo
else
!
! Fully separable form
!
do ns=1,nwfts
if (octs(ns).gt.0.0_DP) then
do n1=1,nbeta
if ( llts(ns).eq.lls(n1).and. &
abs(jjts(ns)-jjs(n1)).lt.1.e-7_DP) then
nst=(llts(ns)+1)*2
ikl=ikk(n1)
do n=1,ikl
f1(n)=betas(n,n1)*phits(n,ns)
enddo
work(n1)=int_0_inf_dr(f1,grid,ikl,nst)
else
work(n1)=0.0_DP
endif
enddo
do n1=1,nbeta
do n2=1,nbeta
epseu=epseu &
+ bmat(n1,n2)*work(n1)*work(n2)*octs(ns)
enddo
enddo
endif
enddo
endif
!
! Now compute the kinetic energy
!
ekin = int_0_inf_dr(f5,grid,grid%mesh,2) - epseu
do ns=1,nwfts
if (octs(ns).gt.0.0_DP) then
ekin=ekin+octs(ns)*enlts(ns)
endif
end do
!
! And the total energy
!
etots= ekin + encl + epseu + ehrt + ecxc + evxt
deallocate(f5)
deallocate(f4)
deallocate(f3)
deallocate(f2)
deallocate(f1)
deallocate(exccc)
return
end subroutine elsdps
|
