File: energy.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (279 lines) | stat: -rw-r--r-- 8,885 bytes parent folder | download
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
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279

! Copyright (C) 2002-2006 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

!BOP
! !ROUTINE: energy
! !INTERFACE:
subroutine energy
! !USES:
use modmain
use moddftu
use modtest
! !DESCRIPTION:
!   Computes the total energy and its individual contributions. The kinetic
!   energy is given by
!   $$ T_s=\sum_i n_i\epsilon_i-\int\rho({\bf r})[v_{\rm C}({\bf r})
!    +v_{\rm xc}({\bf r})]d{\bf r}-\int {\bf m}({\bf r})\cdot
!    ({\bf B}_{\rm xc}({\bf r})+{\bf B}_{\rm ext}({\bf r}))d{\bf r}, $$
!   where $n_i$ are the occupancies and $\epsilon_i$ are the eigenvalues of both
!   the core and valence states; $\rho$ is the density; ${\bf m}$ is the
!   magnetisation density; $v_{\rm C}$ is the Coulomb potential; $v_{\rm xc}$
!   and ${\bf B}_{\rm xc}$ are the exchange-correlation potential and magnetic
!   field, respectively; and ${\bf B}_{\rm ext}$ is the external magnetic field.
!   The Hartree, electron-nuclear and nuclear-nuclear electrostatic energies are
!   combined into the Coulomb energy:
!   \begin{align*}
!    E_{\rm C}&=E_{\rm H}+E_{\rm en}+E_{\rm nn} \\
!             &=\frac{1}{2}V_{\rm C}+E_{\rm Mad},
!   \end{align*}
!   where
!   $$ V_{\rm C}=\int\rho({\bf r})v_{\rm C}({\bf r})d{\bf r} $$
!   is the Coulomb potential energy. The Madelung energy is given by
!   $$ E_{\rm Mad}=\frac{1}{2}\sum_{\alpha}z_{\alpha}R_{\alpha}, $$
!   where
!   $$ R_{\alpha}=\lim_{r\rightarrow 0}\left(v^{\rm C}_{\alpha;00}(r)Y_{00}
!    +\frac{z_{\alpha}}{r}\right) $$
!   for atom $\alpha$, with $v^{\rm C}_{\alpha;00}$ being the $l=0$ component of
!   the spherical harmonic expansion of $v_{\rm C}$ in the muffin-tin, and
!   $z_{\alpha}$ is the nuclear charge. Using the nuclear-nuclear energy
!   determined at the start of the calculation, the electron-nuclear and Hartree
!   energies can be isolated with
!   $$ E_{\rm en}=2\left(E_{\rm Mad}-E_{\rm nn}\right) $$
!   and
!   $$ E_{\rm H}=\frac{1}{2}(E_{\rm C}-E_{\rm en}). $$
!   Finally, the total energy is
!   $$ E=T_s+E_{\rm C}+E_{\rm xc}, $$
!   where $E_{\rm xc}$ is obtained either by integrating the
!   exchange-correlation energy density, or in the case of exact exchange, the
!   explicit calculation of the Fock exchange integral. The energy from the
!   external magnetic fields in the muffin-tins, {\tt bfcmt}, is always removed
!   from the total since these fields are non-physical: their field lines do not
!   close. The energy of the physical external field, {\tt bfieldc}, is also not
!   included in the total because this field, like those in the muffin-tins, is
!   used for breaking spin symmetry and taken to be infintesimal. If this field
!   is intended to be finite, then the associated energy, {\tt engybext}, should
!   be added to the total by hand. See {\tt potxc}, {\tt exxengy} and related
!   subroutines.
!
! !REVISION HISTORY:
!   Created May 2003 (JKD)
!EOP
!BOC
implicit none
! local variables
integer ik,ist,ispn,idm,jdm
integer is,ia,ias,np,n2,i
real(8) cb,sum,f
complex(8) z1
! allocatable arrays
real(8), allocatable :: rfmt(:,:)
complex(8), allocatable :: evecsv(:,:),kmat(:,:),c(:,:)
! external functions
real(8) rfinp
complex(8) zdotc
external rfinp,zdotc
! coupling constant of the external field (g_e/4c)
cb=gfacte/(4.d0*solsc)
!-----------------------------------------------!
!     exchange-correlation potential energy     !
!-----------------------------------------------!
engyvxc=rfinp(rhomt,rhoir,vxcmt,vxcir)
!-----------------------------------------------------!
!     exchange-correlation effective field energy     !
!-----------------------------------------------------!
engybxc=0.d0
do idm=1,ndmag
  engybxc=engybxc+rfinp(magmt(:,:,idm),magir(:,idm),bxcmt(:,:,idm),bxcir(:,idm))
end do
!------------------------------------------!
!     external magnetic field energies     !
!------------------------------------------!
engybext=0.d0
do idm=1,ndmag
  if (ncmag) then
    jdm=idm
  else
    jdm=3
  end if
! energy of physical global field
  engybext=engybext+cb*momtot(idm)*bfieldc(jdm)
end do
!----------------------------------!
!     Coulomb potential energy     !
!----------------------------------!
engyvcl=rfinp(rhomt,rhoir,vclmt,vclir)
!-----------------------!
!     Madelung term     !
!-----------------------!
engymad=0.d0
do ias=1,natmtot
  is=idxis(ias)
  engymad=engymad+0.5d0*spzn(is)*(vclmt(1,ias)-vcln(1,is))*y00
end do
!---------------------------------------------!
!     electron-nuclear interaction energy     !
!---------------------------------------------!
engyen=2.d0*(engymad-engynn)
!------------------------!
!     Hartree energy     !
!------------------------!
engyhar=0.5d0*(engyvcl-engyen)
!------------------------!
!     Coulomb energy     !
!------------------------!
engycl=engynn+engyen+engyhar
!-------------------------!
!     exchange energy     !
!-------------------------!
if ((xctype(1).lt.0).or.(task.eq.5)) then
! exact exchange for OEP-EXX or Hartree-Fock on last self-consistent loop
  if (tlast) then
    call exxengy
! mix exact and DFT exchange energies for hybrid functionals
    if (hybrid) then
      engyx=engyx*hybridc
      engyx=engyx+rfinp(rhomt,rhoir,exmt,exir)
    end if
  else
    engyx=0.d0
  end if
else
! exchange energy from the density
  engyx=rfinp(rhomt,rhoir,exmt,exir)
end if
!----------------------------!
!     correlation energy     !
!----------------------------!
if (task.eq.5) then
  if (hybrid) then
! fraction of DFT correlation energy for hybrid functionals
    engyc=rfinp(rhomt,rhoir,ecmt,ecir)
  else
! zero correlation energy for pure Hartree-Fock
    engyc=0.d0
  end if
else
! correlation energy from the density
  engyc=rfinp(rhomt,rhoir,ecmt,ecir)
end if
!----------------------!
!     DFT+U energy     !
!----------------------!
engydu=0.d0
if (dftu.ne.0) then
  do i=1,ndftu
    is=idftu(1,i)
    do ia=1,natoms(is)
      engydu=engydu+engyadu(ia,i)
    end do
  end do
end if
!----------------------------!
!     sum of eigenvalues     !
!----------------------------!
! core eigenvalues
evalsum=0.d0
do ias=1,natmtot
  is=idxis(ias)
  do ist=1,nstsp(is)
    if (spcore(ist,is)) evalsum=evalsum+occcr(ist,ias)*evalcr(ist,ias)
  end do
end do
! valence eigenvalues
do ik=1,nkpt
  do ist=1,nstsv
    evalsum=evalsum+wkpt(ik)*occsv(ist,ik)*evalsv(ist,ik)
  end do
end do
!------------------------!
!     kinetic energy     !
!------------------------!
! core electron kinetic energy
call energykncr
! total electron kinetic energy
if (task.eq.5) then
! Hartree-Fock case
  engykn=engykncr
! kinetic energy from valence states
  allocate(evecsv(nstsv,nstsv),kmat(nstsv,nstsv),c(nstsv,nstsv))
  do ik=1,nkpt
    call getevecsv(filext,ik,vkl(:,ik),evecsv)
    call getkmat(ik,kmat)
    call zgemm('N','N',nstsv,nstsv,nstsv,zone,kmat,nstsv,evecsv,nstsv,zzero,c, &
     nstsv)
    do ist=1,nstsv
      z1=zdotc(nstsv,evecsv(:,ist),1,c(:,ist),1)
      engykn=engykn+wkpt(ik)*occsv(ist,ik)*dble(z1)
    end do
  end do
  deallocate(evecsv,kmat,c)
else
! Kohn-Sham case
  allocate(rfmt(npmtmax,natmtot))
! remove magnetic field contribution
  sum=0.d0
  do idm=1,ndmag
    do ias=1,natmtot
      is=idxis(ias)
      call rfsht(nrcmt(is),nrcmti(is),bsmt(:,ias,idm),rfmt(:,ias))
    end do
    call rfmtctof(rfmt)
    sum=sum+rfinp(magmt(:,:,idm),magir(:,idm),rfmt,bsir(:,idm))
  end do
! remove integral of w_xc times tau for generalised Kohn-Sham meta-GGA
  if (xcgrad.eq.4) then
    do ispn=1,nspinor
      do ias=1,natmtot
        is=idxis(ias)
        np=npmt(is)
        rfmt(1:np,ias)=taumt(1:np,ias,ispn)-taucr(1:np,ias,ispn)
      end do
      sum=sum+rfinp(rfmt,tauir,wxcmt,wxcir)
    end do
  end if
! remove fixed tensor moment potential matrix contribution
  if (ftmtype.ne.0) then
    n2=(lmmaxdm*nspinor)**2
    do ias=1,natmtot
      z1=zdotc(n2,dmatmt(:,:,:,:,ias),1,vmftm(:,:,:,:,ias),1)
      sum=sum+dble(z1)
    end do
  end if
  engykn=evalsum-engyvcl-engyvxc-sum
  deallocate(rfmt)
end if
!-------------------------------!
!     entropic contribution     !
!-------------------------------!
entrpy=0.d0
engyts=0.d0
! non-zero only for the Fermi-Dirac smearing function
if (stype.eq.3) then
  sum=0.d0
  do ik=1,nkpt
    do ist=1,nstsv
      f=occsv(ist,ik)/occmax
      if ((f.gt.0.d0).and.(f.lt.1.d0)) then
        sum=sum+wkpt(ik)*(f*log(f)+(1.d0-f)*log(1.d0-f))
      end if
    end do
  end do
! entropy
  entrpy=-occmax*kboltz*sum
! contribution to free energy
  engyts=-swidth*entrpy/kboltz
end if
!----------------------!
!     total energy     !
!----------------------!
engytot=engykn+0.5d0*engyvcl+engymad+engyx+engyc+engyts
! add the DFT+U correction if required
if (dftu.ne.0) engytot=engytot+engydu
! write total energy to test file
call writetest(0,'total energy',tol=1.d-4,rv=engytot)
return
end subroutine
!EOC