totale.f90 11.6 KB
Newer Older
1
2
3
4
5
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
6
7
MODULE m_totale
CONTAINS
Matthias Redies's avatar
Matthias Redies committed
8
9
  SUBROUTINE totale(atoms,sphhar,stars,vacuum,dimension,sym,input,noco, &
       cell,oneD, xcpot,hybrid,vTot,vCoul,it,den,results)
10
11
    !
    !     ***************************************************
12
    !     subroutine calculates the total energy 
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
    !     ***************************************************
    !     single particle energies
    !     SEIGC  sum of the eigenvalues of the core states
    !            calculated in cdngen.f
    !     SEIGSCV  sum of the eigenvalues of the semicore and valence states
    !              calculated in fermie.f 
    !     TS         : entropy contribution to the free energy
    !     SEIGC,SEIGSCV, TS are calculated in fermie.f
    !     ***************************************************
    !     TE_VCOUL  :   charge density-coulomb potential integral
    !     TE_VEFF:   charge density-effective potential integral
    !     TE_EXC :   charge density-ex-corr.energy density integral
    !                 exchange-correlation energy
    !     TE_VCOUL,TE_VEFF,TE_EXC are calculated in vgen.f
    !     VMD :   Madelung term
    !     ***************************************************
    !     TOTE    :   total energy due to all electrons
    !     TOTE = SEIGC + SEIGSCV + TE_VCOUL/2 -TE_VEFF + TE_EXC + VMD
    !
    !     if HF calculation/hybrid-functional calculation :
    !     TOTE = SEIGC + SEIGSCV + TE_VCOUL/2 -TE_VEFF + TE_EXC_loc + VMD - 1/2 E_FOCK
    !
    !     E_FOCK: sum of diagonal elements of fock matrix
    !
    !     ***************************************************
    !     FREE ENRGY: F = TOTE - TS
    !     total electron energy extrapolated for T->0
    !     E0 = TOTE - TS/2
    !     ***************************************************
    !
    USE m_intgr    , ONLY : intgr3 
44
    USE m_constants
45
46
47
    USE m_force_a4
    USE m_force_a3
    USE m_forcew
48
    USE m_cdn_io
49
    USE m_types
Gregor Michalicek's avatar
Gregor Michalicek committed
50
    USE m_xmlOutput
51
52
    IMPLICIT NONE

53
    CLASS(t_xcpot),INTENT(IN)       :: xcpot
54
55
56
57
58
59
60
61
62
63
    TYPE(t_oneD),INTENT(IN)         :: oneD
    TYPE(t_hybrid),INTENT(IN)       :: hybrid
    TYPE(t_input),INTENT(IN)        :: input
    TYPE(t_vacuum),INTENT(IN)       :: vacuum
    TYPE(t_noco),INTENT(IN)         :: noco
    TYPE(t_sym),INTENT(IN)          :: sym
    TYPE(t_stars),INTENT(IN)        :: stars
    TYPE(t_cell),INTENT(IN)         :: cell
    TYPE(t_sphhar),INTENT(IN)       :: sphhar
    TYPE(t_atoms),INTENT(IN)        :: atoms
64
    TYPE(t_dimension),INTENT(IN)    :: dimension
65
    TYPE(t_potden),INTENT(IN)       :: vTot,vCoul
66
    TYPE(t_potden),INTENT(IN)       :: den
Matthias Redies's avatar
Matthias Redies committed
67
    TYPE(t_results),INTENT(INOUT)   :: results
68
69
70
    !     ..
    !     .. Scalar Arguments ..
    INTEGER,INTENT (IN) :: it      
71
72
73

    ! Local type instances

74
    !     .. Local Scalars ..
75
    REAL rhs,totz, eigSum, fermiEnergyTemp
76
    INTEGER n,j,nt,i, archiveType
77
    LOGICAL l_qfix
78
79

    !     .. Local Arrays ..
80
    REAL vmd(atoms%ntype),zintn_r(atoms%ntype)
81
    REAL dpj(atoms%jmtd),mt(atoms%jmtd,atoms%ntype)
82
    CHARACTER(LEN=20) :: attributes(3)
83

84
    !CALL den%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,.FALSE.,POTDEN_TYPE_DEN)
85

86
87
88
89
90
    WRITE (6,FMT=8000)
    WRITE (16,FMT=8000)
8000 FORMAT (/,/,/,5x,'t o t a l  e n e r g y')
    !
    !      ---> sum of eigenvalues (core, semicore and valence states)
91
92
93
    !
    eigSum = results%seigscv + results%seigc
    results%tote = eigSum
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
    WRITE (6,FMT=8010) results%tote
    WRITE (16,FMT=8010) results%tote
8010 FORMAT (/,10x,'sum of eigenvalues =',t40,f20.10)
    !
    !      ---> add contribution of coulomb potential
    !
    results%tote = results%tote + 0.5e0*results%te_vcoul
    WRITE (6,FMT=8020) results%te_vcoul
    WRITE (16,FMT=8020) results%te_vcoul
8020 FORMAT (/,10x,'density-coulomb potential integral =',t40,f20.10)
    !
    !      ---> add contribution of effective potential
    !
    results%tote = results%tote - results%te_veff
    WRITE (6,FMT=8030) results%te_veff
    WRITE (16,FMT=8030) results%te_veff
8030 FORMAT (/,10x,'density-effective potential integral =',t40,f20.10)
    !
    !      ---> add contribution of exchange-correlation energy
    !
    results%tote = results%tote + results%te_exc
    WRITE (6,FMT=8040) results%te_exc
    WRITE (16,FMT=8040) results%te_exc
117
8040 FORMAT (/,10x,'charge density-ex.-corr.energy density integral=', t40,ES20.10)
118
119
120
    !
    !      ---> Fock exchange contribution 
    !
121
    IF (xcpot%is_hybrid()) THEN
122
123
124
       !IF (xcpot%is_name("exx")) THEN
       !   results%tote = results%tote + 0.5e0*results%te_hfex%valence
       !ELSE
125
          results%tote = results%tote - 0.5e0*results%te_hfex%valence + 0.5e0*results%te_hfex%core
126
       !END IF
127
    ENDIF
128
129
130
131
132
133
134
135
136
137
138
    WRITE (6,FMT=8100)  0.5e0*results%te_hfex%valence
    WRITE (16,FMT=8100) 0.5e0*results%te_hfex%valence
    WRITE (6,FMT=8101)  0.5e0*results%te_hfex%core
    WRITE (16,FMT=8101) 0.5e0*results%te_hfex%core
8100 FORMAT (/,10x,'Fock-exchange energy (valence)=',t40,f20.10)
8101 FORMAT (10x,'Fock-exchange energy (core)=',t40,f20.10)


    !     ----> VM terms
    !     ---> reload the density
    !
139
140
    !archiveType = CDN_ARCHIVE_TYPE_CDN1_const
    !IF (noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_CDN_const
141

142
143
    !CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,&
    !                 CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,den)
144

145
146

    ! CLASSICAL HELLMAN-FEYNMAN FORCE
147
    CALL force_a3(atoms,sphhar, input, den%mt,vCoul%mt, results%force)
148

149
    IF (input%l_f) THEN
150
151
       ! core contribution to force: needs TOTAL POTENTIAL and core charge
       CALL force_a4(atoms,sphhar,input,dimension, vTot%mt, results%force)
152

153
154
155
156
157
    ENDIF

    !-for
    !     ---> add spin-up and spin-down charge density for lh=0
    !
158
159
160
    DO  n = 1,atoms%ntype
       DO  i = 1,atoms%jri(n)
          mt(i,n) = den%mt(i,0,n,1) + den%mt(i,0,n,input%jspins)
161
       ENDDO
162
163
164
    ENDDO
    IF (input%jspins.EQ.1) mt=mt/2 !we just added the same value twice
    
165
166
167
168
169
    !
    ! ----> coulomb interaction between electrons and nuclei of different m.t.s
    !
    DO  n = 1,atoms%ntype
       DO  j = 1,atoms%jri(n)
170
          dpj(j) = mt(j,n)/atoms%rmsh(j,n)
171
172
173
174
175
       ENDDO
       CALL intgr3(dpj,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),rhs)
       !
       results%tote = results%tote - atoms%neq(n)*atoms%zatom(n)*sfp_const*rhs/2.
       !
176
177
178
       zintn_r(n) = atoms%neq(n)*atoms%zatom(n)*sfp_const*rhs/2.
       WRITE (6,FMT=8045) zintn_r(n)
       WRITE (16,FMT=8045) zintn_r(n)
179
       CALL intgr3(mt(1,n),atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),totz)
180
       vmd(n) = atoms%rmt(n)*vCoul%mt(atoms%jri(n),0,n,1)/sfp_const + atoms%zatom(n) - totz*sfp_const
181
182
183
184
       vmd(n) = -atoms%neq(n)*atoms%zatom(n)*vmd(n)/ (2.*atoms%rmt(n))
       WRITE (6,FMT=8050) n,vmd(n)
       WRITE (16,FMT=8050) n,vmd(n)
       results%tote = results%tote + vmd(n)
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
    ENDDO
    IF (atoms%n_u.GT.0) THEN
       WRITE ( 6,FMT=8090) results%e_ldau
       WRITE (16,FMT=8090) results%e_ldau
       results%tote = results%tote - results%e_ldau             ! gu test
    ENDIF
    ! print 'HF' before total energy to make it grepable
    IF ( .NOT. hybrid%l_calhf ) THEN
       WRITE ( 6,FMT=8060) results%tote
       WRITE (16,FMT=8060) results%tote
    ELSE
       WRITE ( 6,FMT=8061) results%tote
       WRITE (16,FMT=8061) results%tote
    END IF

200
    CALL force_w(input,atoms,sym,results,cell,oneD,vacuum)
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
    !
    !     ---> calculate the free energy and the ground state energy,
    !          extrapolated for T->0
    !
    ! print 'HF' before all energies to make them grepable
    IF ( .NOT. hybrid%l_calhf ) THEN
       WRITE ( 6,FMT=8065) results%ts
       WRITE (16,FMT=8065) results%ts
       WRITE ( 6,FMT=8070) results%tote-results%ts
       WRITE (16,FMT=8070) results%tote-results%ts
       WRITE ( 6,FMT=8080) results%tote-0.5e0*results%ts
       WRITE (16,FMT=8080) results%tote-0.5e0*results%ts
    ELSE
       WRITE ( 6,FMT=8066) results%ts
       WRITE (16,FMT=8066) results%ts
       WRITE ( 6,FMT=8071) results%tote-results%ts
       WRITE (16,FMT=8071) results%tote-results%ts
       WRITE ( 6,FMT=8081) results%tote-0.5e0*results%ts
       WRITE (16,FMT=8081) results%tote-0.5e0*results%ts
    END IF
221

222
223
224
    WRITE(attributes(1),'(f20.10)') results%tote
    WRITE(attributes(2),'(a)') 'Htr'
    WRITE(attributes(3),'(a)') 'HF'
225
    IF (hybrid%l_calhf) THEN
226
       CALL openXMLElementForm('totalEnergy',(/'value  ','units  ','comment'/),attributes,reshape((/40,20/),(/1,2/)))
227
    ELSE
228
       CALL openXMLElementForm('totalEnergy',(/'value','units'/),attributes(1:2),reshape((/40,20/),(/1,2/)))
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
    END IF
    CALL openXMLElementFormPoly('sumOfEigenvalues',(/'value'/),(/eigSum/),reshape((/32,20/),(/1,2/)))
    CALL writeXMLElementFormPoly('coreElectrons',(/'value'/),(/results%seigc/),reshape((/32,20/),(/1,2/)))
    CALL writeXMLElementFormPoly('valenceElectrons',(/'value'/),(/results%seigscv/),reshape((/29,20/),(/1,2/)))
    CALL closeXMLElement('sumOfEigenvalues')
    CALL writeXMLElementFormPoly('densityCoulombPotentialIntegral',(/'value'/),(/results%te_vcoul/),reshape((/17,20/),(/1,2/)))
    CALL writeXMLElementFormPoly('densityEffectivePotentialIntegral',(/'value'/),(/results%te_veff/),reshape((/15,20/),(/1,2/)))
    CALL writeXMLElementFormPoly('chargeDenXCDenIntegral',(/'value'/),(/results%te_exc/),reshape((/26,20/),(/1,2/)))
    CALL writeXMLElementFormPoly('FockExchangeEnergyValence',(/'value'/),(/0.5e0*results%te_hfex%valence/),reshape((/23,20/),(/1,2/)))
    CALL writeXMLElementFormPoly('FockExchangeEnergyCore',(/'value'/),(/0.5e0*results%te_hfex%core/),reshape((/26,20/),(/1,2/)))
    DO  n = 1,atoms%ntype
       CALL openXMLElementPoly('atomTypeDependentContributions',(/'atomType'/),(/n/))
       CALL writeXMLElementFormPoly('electronNucleiInteractionDifferentMTs',(/'value'/),(/zintn_r(n)/),reshape((/8,20/),(/1,2/)))
       CALL writeXMLElementFormPoly('MadelungTerm',(/'value'/),(/vmd(n)/),reshape((/33,20/),(/1,2/)))
       CALL closeXMLElement('atomTypeDependentContributions')
    END DO
    IF (atoms%n_u.GT.0) THEN
246
       CALL writeXMLElementFormPoly('dftUCorrection',(/'value'/),(/results%e_ldau/),reshape((/34,20/),(/1,2/)))
247
248
249
250
    END IF
    CALL writeXMLElementFormPoly('tkbTimesEntropy',(/'value'/),(/results%ts/),reshape((/33,20/),(/1,2/)))
    CALL writeXMLElementFormPoly('freeEnergy',(/'value'/),(/results%tote-results%ts/),reshape((/38,20/),(/1,2/)))
    CALL writeXMLElementFormPoly('extrapolationTo0K',(/'value'/),(/results%tote-0.5e0*results%ts/),reshape((/31,20/),(/1,2/)))
Gregor Michalicek's avatar
Gregor Michalicek committed
251
    CALL closeXMLElement('totalEnergy')
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
8060 FORMAT (/,/,' ---->    input%total energy=',t40,f20.10,' htr')
8061 FORMAT (/,/,' ----> HF input%total energy=',t40,f20.10,' htr')
8050 FORMAT (/,10x,'Madelung term for atom type:',i3,t40,f20.10)
8045 FORMAT (/,10x,'el.-nucl. inter. diff. m.t.',t40,f20.10)
8065 FORMAT (/,/,' ---->    (input%tkb*entropy) TS=',t40,f20.10,' htr')
8066 FORMAT (/,/,' ----> HF (input%tkb*entropy) TS=',t40,f20.10,' htr')
8070 FORMAT (/,/,' ---->    free energy=',t40,f20.10,' htr')
8071 FORMAT (/,/,' ----> HF free energy=',t40,f20.10,' htr')
8080 FORMAT (/,/,'      extrapolation for T->0',&
               /,' ---->    input%total electron energy=',t40,f20.10,' htr')
8081 FORMAT (/,/,'      extrapolation for T->0',&
               /,' ----> HF input%total electron energy=',t40,f20.10,' htr')
8090 FORMAT (/,/,' ---->    correction for lda+U =',t40,f20.10,' htr')

  END SUBROUTINE totale
END MODULE m_totale