Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 : MODULE m_totale
7 : CONTAINS
8 668 : SUBROUTINE totale(fmpi,atoms,sphhar,stars,vacuum, &
9 : sym,input,noco,cell , xcpot,hybdat,vTot,vCoul,it,den,results)
10 : !
11 : ! ***************************************************
12 : ! subroutine calculates the total energy
13 : ! ***************************************************
14 : ! single particle energies
15 : ! SEIGC sum of the eigenvalues of the core states
16 : ! calculated in cdngen.f
17 : ! SEIGV sum of the eigenvalues of the semicore and valence states
18 : ! calculated in fermie.f
19 : ! TS : entropy contribution to the free energy
20 : ! SEIGC,SEIGV, TS are calculated in fermie.f
21 : ! ***************************************************
22 : ! TE_VCOUL : charge density-coulomb potential integral
23 : ! TE_VEFF: charge density-effective potential integral
24 : ! TE_EXC : charge density-ex-corr.energy density integral
25 : ! exchange-correlation energy
26 : ! TE_VCOUL,TE_VEFF,TE_EXC are calculated in vgen.f
27 : ! VMD : Madelung term
28 : ! ***************************************************
29 : ! TOTE : total energy due to all electrons
30 : ! TOTE = SEIGC + SEIGV + TE_VCOUL/2 -TE_VEFF + TE_EXC + VMD
31 : !
32 : ! if HF calculation/hybinp-functional calculation :
33 : ! TOTE = SEIGC + SEIGV + TE_VCOUL/2 -TE_VEFF + TE_EXC_loc + VMD - 1/2 E_FOCK
34 : !
35 : ! E_FOCK: sum of diagonal elements of fock matrix
36 : !
37 : ! ***************************************************
38 : ! FREE ENRGY: F = TOTE - TS
39 : ! total electron energy extrapolated for T->0
40 : ! E0 = TOTE - TS/2
41 : ! ***************************************************
42 : !
43 : USE m_intgr , ONLY : intgr3
44 : USE m_constants
45 : USE m_force_a4
46 : USE m_force_a3
47 : USE m_force_a4_add ! Klueppelberg (force level 1)
48 : USE m_force_sf ! Klueppelberg (force level 3)
49 : USE m_forcew
50 : USE m_cdn_io
51 : USE m_types
52 : USE m_xmlOutput
53 : use m_judft
54 : USE m_vdWfleur_grimme
55 :
56 : IMPLICIT NONE
57 : TYPE(t_mpi),INTENT(IN) :: fmpi
58 : TYPE(t_results),INTENT(INOUT) :: results
59 : CLASS(t_xcpot),INTENT(IN) :: xcpot
60 :
61 : TYPE(t_hybdat),INTENT(IN) :: hybdat
62 : TYPE(t_input),INTENT(IN) :: input
63 : TYPE(t_vacuum),INTENT(IN) :: vacuum
64 : TYPE(t_noco),INTENT(IN) :: noco
65 : TYPE(t_sym),INTENT(IN) :: sym
66 : TYPE(t_stars),INTENT(IN) :: stars
67 : TYPE(t_cell),INTENT(IN) :: cell
68 : TYPE(t_sphhar),INTENT(IN) :: sphhar
69 : TYPE(t_atoms),INTENT(IN) :: atoms
70 :
71 : TYPE(t_potden),INTENT(IN) :: vTot,vCoul
72 : TYPE(t_potden),INTENT(IN) :: den
73 : ! ..
74 : ! .. Scalar Arguments ..
75 : INTEGER,INTENT (IN) :: it
76 :
77 : ! Local type instances
78 :
79 : ! .. Local Scalars ..
80 : REAL rhs,totz, eigSum, fermiEnergyTemp
81 : INTEGER n,j,nt,i, archiveType,jsp
82 : LOGICAL l_qfix
83 :
84 : ! .. Local Arrays ..
85 668 : REAL vmd(atoms%ntype),zintn_r(atoms%ntype)
86 668 : REAL dpj(atoms%jmtd),mt(atoms%jmtd,atoms%ntype)
87 : CHARACTER(LEN=20) :: attributes(3)
88 :
89 : !CALL den%init(stars,atoms,sphhar,vacuum,noco ,input%jspins,.FALSE.,POTDEN_TYPE_DEN)
90 668 : IF (fmpi%irank==0) THEN
91 334 : WRITE (oUnit,FMT=8000)
92 : 8000 FORMAT (/,/,/,5x,'t o t a l e n e r g y')
93 : !
94 : ! ---> sum of eigenvalues (core, semicore and valence states)
95 : !
96 334 : eigSum = results%seigv + results%seigc
97 334 : results%tote = eigSum
98 334 : WRITE (oUnit,FMT=8010) results%tote
99 : 8010 FORMAT (/,10x,'sum of eigenvalues =',t40,f20.10)
100 : !
101 : ! ---> add contribution of coulomb potential
102 : !
103 334 : results%tote = results%tote + 0.5e0*results%te_vcoul
104 334 : WRITE (oUnit,FMT=8020) results%te_vcoul
105 : 8020 FORMAT (/,10x,'density-coulomb potential integral =',t40,f20.10)
106 : !
107 : ! ---> add contribution of effective potential
108 : !
109 334 : results%tote = results%tote - results%te_veff
110 334 : WRITE (oUnit,FMT=8030) results%te_veff
111 : 8030 FORMAT (/,10x,'density-effective potential integral =',t40,f20.10)
112 : !
113 : ! ---> add contribution of exchange-correlation energy
114 : !
115 334 : results%tote = results%tote + results%te_exc
116 334 : WRITE (oUnit,FMT=8040) results%te_exc
117 : 8040 FORMAT (/,10x,'charge density-ex.-corr.energy density integral=', t40,f20.10)
118 : !
119 : ! ---> Fock exchange contribution
120 : !
121 334 : IF (xcpot%is_hybrid()) THEN
122 : !IF (xcpot%is_name("exx")) THEN
123 : ! results%tote = results%tote + 0.5e0*results%te_hfex%valence
124 : !ELSE
125 93 : results%tote = results%tote - 0.5e0*results%te_hfex%valence + 0.5e0*results%te_hfex%core
126 : !END IF
127 93 : write (oUnit,*) 'Fock-exchange energy (valence)= ' // float2str(0.5e0*results%te_hfex%valence)
128 93 : write (oUnit,*) 'Fock-exchange energy (core)= ' // float2str(0.5e0*results%te_hfex%core)
129 : ENDIF
130 :
131 :
132 : ! ----> VM terms
133 : ! ---> reload the density
134 : !
135 : !archiveType = CDN_ARCHIVE_TYPE_CDN1_const
136 : !IF (noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_CDN_const
137 :
138 : !CALL readDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,archiveType,&
139 : ! CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,den)
140 :
141 :
142 : ! CLASSICAL HELLMAN-FEYNMAN FORCE
143 334 : CALL force_a3(atoms,sym,sphhar, input, den%mt,vCoul%mt, results%force)
144 :
145 334 : IF (input%l_f) THEN
146 : ! core contribution to force: needs TOTAL POTENTIAL and core charge
147 :
148 29 : IF (input%ctail.AND.(input%f_level.GE.1)) THEN
149 : ! Add core correction to forces from tails of core states
150 : ! Klueppelberg, Sep'12 (force level 1)
151 1 : CALL force_a4_add(atoms,input,results)
152 : END IF
153 :
154 29 : CALL force_a4(atoms,sym,sphhar,input, vTot%mt, results%force)
155 :
156 : END IF
157 :
158 : !-for
159 : ! ---> add spin-up and spin-down charge density for lh=0
160 : !
161 418588 : mt=0.0
162 919 : DO n = 1,atoms%ntype
163 401036 : DO i = 1,atoms%jri(n)
164 400702 : mt(i,n) = den%mt(i,0,n,1) + den%mt(i,0,n,input%jspins)
165 : ENDDO
166 : ENDDO
167 191407 : IF (input%jspins.EQ.1) mt=mt/2 !we just added the same value twice
168 :
169 : !
170 : ! ----> coulomb interaction between electrons and nuclei of different m.t.s
171 : !
172 919 : DO n = 1,atoms%ntype
173 400702 : DO j = 1,atoms%jri(n)
174 400702 : dpj(j) = mt(j,n)/atoms%rmsh(j,n)
175 : ENDDO
176 585 : CALL intgr3(dpj,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),rhs)
177 : !
178 585 : zintn_r(n) = atoms%neq(n)*atoms%zatom(n)*sfp_const*rhs/2.
179 585 : results%tote = results%tote - zintn_r(n)
180 : !
181 585 : WRITE (oUnit,FMT=8045) zintn_r(n)
182 585 : CALL intgr3(mt(1,n),atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),totz)
183 585 : vmd(n) = atoms%rmt(n)*vCoul%mt(atoms%jri(n),0,n,1)/sfp_const + atoms%zatom(n) - totz*sfp_const
184 585 : vmd(n) = -atoms%neq(n)*atoms%zatom(n)*vmd(n)/ (2.*atoms%rmt(n))
185 585 : WRITE (oUnit,FMT=8050) n,vmd(n)
186 2089 : results%tote = results%tote + vmd(n)
187 : ENDDO
188 334 : IF (atoms%n_u+atoms%n_hia.GT.0) THEN
189 32 : WRITE (oUnit,FMT=8090) results%e_ldau
190 32 : results%tote = results%tote - results%e_ldau ! gu test
191 : ENDIF
192 : ! print 'HF' before total energy to make it grepable
193 334 : IF ( .NOT. hybdat%l_calhf ) THEN
194 328 : WRITE (oUnit,FMT=8060) results%tote
195 : ELSE
196 6 : WRITE (oUnit,FMT=8061) results%tote
197 : END IF
198 : !
199 : ! ---> calculate the free energy and the ground state energy,
200 : ! extrapolated for T->0
201 : !
202 : ! print 'HF' before all energies to make them grepable
203 334 : IF ( .NOT. hybdat%l_calhf ) THEN
204 328 : WRITE (oUnit,FMT=8065) results%ts
205 328 : WRITE (oUnit,FMT=8070) results%tote-results%ts
206 328 : WRITE (oUnit,FMT=8080) results%tote-0.5e0*results%ts
207 : ELSE
208 6 : WRITE (oUnit,FMT=8066) results%ts
209 6 : WRITE (oUnit,FMT=8071) results%tote-results%ts
210 6 : WRITE (oUnit,FMT=8081) results%tote-0.5e0*results%ts
211 : END IF
212 :
213 : ! vdW D3 Grimme contribution
214 334 : IF (btest(input%vdW,0)) THEN
215 0 : if (.not.allocated(results%force_vdw)) ALLOCATE(results%force_vdW(3,atoms%ntype))
216 0 : call vdW_fleur_grimme(input,atoms,sym,cell,results%e_vdW,results%force_vdw)
217 : ENDIF
218 : !ADD vdW contribution to energy (is zero if no vdW was evaluated)
219 334 : results%tote=results%tote+results%e_vdW
220 334 : WRITE(attributes(1),'(f20.10)') results%tote
221 334 : WRITE(attributes(2),'(a)') 'Htr'
222 334 : WRITE(attributes(3),'(a)') 'HF'
223 334 : IF (hybdat%l_calhf) THEN
224 24 : CALL openXMLElementForm('totalEnergy',(/'value ','units ','comment'/),attributes,reshape((/40,20/),(/1,2/)))
225 : ELSE
226 984 : CALL openXMLElementForm('totalEnergy',(/'value','units'/),attributes(1:2),reshape((/40,20/),(/1,2/)))
227 : END IF
228 1002 : CALL openXMLElementFormPoly('sumOfEigenvalues',(/'value'/),(/eigSum/),reshape((/32,20/),(/1,2/)))
229 1002 : CALL writeXMLElementFormPoly('coreElectrons',(/'value'/),(/results%seigc/),reshape((/32,20/),(/1,2/)))
230 1002 : CALL writeXMLElementFormPoly('valenceElectrons',(/'value'/),(/results%seigv/),reshape((/29,20/),(/1,2/)))
231 334 : CALL closeXMLElement('sumOfEigenvalues')
232 1002 : CALL writeXMLElementFormPoly('densityCoulombPotentialIntegral',(/'value'/),(/results%te_vcoul/),reshape((/17,20/),(/1,2/)))
233 1002 : CALL writeXMLElementFormPoly('densityEffectivePotentialIntegral',(/'value'/),(/results%te_veff/),reshape((/15,20/),(/1,2/)))
234 1002 : CALL writeXMLElementFormPoly('chargeDenXCDenIntegral',(/'value'/),(/results%te_exc/),reshape((/26,20/),(/1,2/)))
235 1002 : CALL writeXMLElementFormPoly('FockExchangeEnergyValence',(/'value'/),(/0.5e0*results%te_hfex%valence/),reshape((/23,20/),(/1,2/)))
236 1002 : CALL writeXMLElementFormPoly('FockExchangeEnergyCore',(/'value'/),(/0.5e0*results%te_hfex%core/),reshape((/26,20/),(/1,2/)))
237 334 : if (btest(input%vdw,0).or.btest(input%vdW,1)) call writeXMLElementFormPoly('vdWEnergy',(/'value'/),(/results%e_vdW/),reshape((/17,20/),(/1,2/)))
238 919 : DO n = 1,atoms%ntype
239 1755 : CALL openXMLElementPoly('atomTypeDependentContributions',(/'atomType'/),(/n/))
240 1755 : CALL writeXMLElementFormPoly('electronNucleiInteractionDifferentMTs',(/'value'/),(/zintn_r(n)/),reshape((/8,20/),(/1,2/)))
241 1755 : CALL writeXMLElementFormPoly('MadelungTerm',(/'value'/),(/vmd(n)/),reshape((/33,20/),(/1,2/)))
242 919 : CALL closeXMLElement('atomTypeDependentContributions')
243 : END DO
244 334 : IF (atoms%n_u+atoms%n_hia.GT.0) THEN
245 96 : CALL writeXMLElementFormPoly('dftUCorrection',(/'value'/),(/results%e_ldau/),reshape((/34,20/),(/1,2/)))
246 : END IF
247 1002 : CALL writeXMLElementFormPoly('tkbTimesEntropy',(/'value'/),(/results%ts/),reshape((/33,20/),(/1,2/)))
248 1002 : CALL writeXMLElementFormPoly('freeEnergy',(/'value'/),(/results%tote-results%ts/),reshape((/38,20/),(/1,2/)))
249 1002 : CALL writeXMLElementFormPoly('extrapolationTo0K',(/'value'/),(/results%tote-0.5e0*results%ts/),reshape((/31,20/),(/1,2/)))
250 334 : CALL closeXMLElement('totalEnergy')
251 : 8060 FORMAT (/,/,' ----> total energy=',t40,f20.10,' htr')
252 : 8061 FORMAT (/,/,' ----> HF total energy=',t40,f20.10,' htr')
253 : 8050 FORMAT (/,10x,'Madelung term for atom type:',i3,t40,f20.10)
254 : 8045 FORMAT (/,10x,'el.-nucl. inter. diff. m.t.',t40,f20.10)
255 : 8065 FORMAT (/,/,' ----> (input%tkb*entropy) TS=',t40,f20.10,' htr')
256 : 8066 FORMAT (/,/,' ----> HF (input%tkb*entropy) TS=',t40,f20.10,' htr')
257 : 8070 FORMAT (/,/,' ----> free energy=',t40,f20.10,' htr')
258 : 8071 FORMAT (/,/,' ----> HF free energy=',t40,f20.10,' htr')
259 : 8080 FORMAT (/,/,' extrapolation for T->0',&
260 : /,' ----> total electron energy=',t40,f20.10,' htr')
261 : 8081 FORMAT (/,/,' extrapolation for T->0',&
262 : /,' ----> HF total electron energy=',t40,f20.10,' htr')
263 : 8090 FORMAT (/,/,' ----> correction for lda+U =',t40,f20.10,' htr')
264 : ENDIF
265 :
266 : ! Klueppelberg (force level 3)
267 668 : IF (input%l_f.AND.(input%f_level.GE.3)) THEN
268 4 : DO jsp=1,input%jspins
269 4 : CALL exit_sf(jsp,atoms,results%force)
270 : END DO
271 : END IF
272 668 : CALL force_w(fmpi,input,atoms,sym,results,cell ,vacuum)
273 :
274 664 : END SUBROUTINE totale
275 : END MODULE m_totale
|
