Octopus
energy_calc.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
22 use batch_oct_m
23 use comm_oct_m
24 use debug_oct_m
26 use energy_oct_m
31 use global_oct_m
32 use grid_oct_m
35 use ions_oct_m
37 use lda_u_oct_m
39 use mesh_oct_m
45 use pcm_oct_m
47 use smear_oct_m
48 use space_oct_m
52 use unit_oct_m
55 use xc_oct_m
56
57 implicit none
58
59 private
60 public :: &
66
67contains
68
69 ! ---------------------------------------------------------
73 subroutine energy_calc_total(namespace, space, hm, gr, st, ext_partners, iunit, full)
74 type(namespace_t), intent(in) :: namespace
75 class(space_t), intent(in) :: space
76 type(hamiltonian_elec_t), intent(inout) :: hm
77 type(grid_t), intent(in) :: gr
78 type(states_elec_t), intent(inout) :: st
79 type(partner_list_t), intent(in) :: ext_partners
80 integer, optional, intent(in) :: iunit
81 logical, optional, intent(in) :: full
82
83 logical :: full_
84 type(states_elec_t) :: xst
86 type(gauge_field_t), pointer :: gfield
87
88 push_sub(energy_calc_total)
89
90 full_ = .false.
91 if (present(full)) full_ = full
92
93 hm%energy%eigenvalues = states_elec_eigenvalues_sum(st)
94
95 if (full_ .or. hm%theory_level == hartree .or. hm%theory_level == hartree_fock &
96 .or. hm%theory_level == generalized_kohn_sham_dft) then
97 if (states_are_real(st)) then
98 hm%energy%kinetic = denergy_calc_electronic(namespace, hm, gr%der, st, terms = term_kinetic)
99 hm%energy%extern_local = denergy_calc_electronic(namespace, hm, gr%der, st, terms = term_local_external)
100 hm%energy%extern_non_local = denergy_calc_electronic(namespace, hm, gr%der, st, terms = term_non_local_potential)
101 hm%energy%extern = hm%energy%extern_local + hm%energy%extern_non_local
102 else
103 hm%energy%kinetic = zenergy_calc_electronic(namespace, hm, gr%der, st, terms = term_kinetic)
104
105 hm%energy%extern_local = zenergy_calc_electronic(namespace, hm, gr%der, st, terms = term_local_external)
106
107 hm%energy%extern_non_local = zenergy_calc_electronic(namespace, hm, gr%der, st, terms = term_non_local_potential)
108 hm%energy%extern = hm%energy%extern_local + hm%energy%extern_non_local
109 end if
110 end if
111
112 ! Computes the HF energy when not doing the ACE
113 if( full_ .and. (hm%theory_level == hartree_fock .or. &
114 (hm%theory_level == generalized_kohn_sham_dft .and. family_is_hybrid(hm%xc))) &
115 .and. .not. hm%exxop%useACE ) then
116 call xst%nullify()
117 if (states_are_real(st)) then
118 call dexchange_operator_compute_potentials(hm%exxop, namespace, space, gr, st, &
119 xst, hm%kpoints, hm%energy%exchange_hf)
120 else
121 call zexchange_operator_compute_potentials(hm%exxop, namespace, space, gr, st, &
122 xst, hm%kpoints, hm%energy%exchange_hf)
123 end if
124 call states_elec_end(xst)
125 else
126 if(.not. hm%exxop%useACE) hm%energy%exchange_hf = m_zero
127 end if
128
129 if (hm%pcm%run_pcm) then
130 hm%pcm%counter = hm%pcm%counter + 1
131 if (hm%pcm%localf) then
132 call pcm_elect_energy(hm%ions, hm%pcm, hm%energy%int_ee_pcm, hm%energy%int_en_pcm, &
133 hm%energy%int_ne_pcm, hm%energy%int_nn_pcm, &
134 e_int_e_ext = hm%energy%int_e_ext_pcm, &
135 e_int_n_ext = hm%energy%int_n_ext_pcm )
136 else
137 call pcm_elect_energy(hm%ions, hm%pcm, hm%energy%int_ee_pcm, hm%energy%int_en_pcm, &
138 hm%energy%int_ne_pcm, hm%energy%int_nn_pcm )
139 end if
140 end if
141
142 select case (hm%theory_level)
144 hm%energy%total = hm%ep%eii + hm%energy%eigenvalues
145
146 case (hartree)
147 hm%energy%total = hm%ep%eii + m_half * (hm%energy%eigenvalues + hm%energy%kinetic + hm%energy%extern)
148
150 hm%energy%total = hm%ep%eii + &
151 m_half*(hm%energy%eigenvalues + hm%energy%kinetic + hm%energy%extern - hm%energy%intnvxc &
152 - hm%energy%int_dft_u) &
153 + hm%energy%exchange + hm%energy%exchange_hf + hm%energy%correlation &
154 + hm%energy%vdw - hm%energy%intnvstatic + hm%energy%dft_u
155
156 ! FIXME: pcm terms are only added to total energy in DFT case
157
158 case (kohn_sham_dft)
159 hm%energy%total = hm%ep%eii + hm%energy%eigenvalues &
160 - hm%energy%hartree + hm%energy%exchange + hm%energy%correlation + hm%energy%vdw - hm%energy%intnvxc &
161 - hm%energy%pcm_corr + hm%energy%int_ee_pcm + hm%energy%int_en_pcm &
162 + hm%energy%int_nn_pcm + hm%energy%int_ne_pcm &
163 + hm%energy%int_e_ext_pcm + hm%energy%int_n_ext_pcm &
164 + hm%energy%dft_u - hm%energy%int_dft_u - hm%energy%intnvstatic &
165 + hm%energy%photon_exchange
167 end select
168
169 hm%energy%entropy = smear_calc_entropy(st%smear, st%eigenval, st%nik, st%nst, st%kweights, st%occ)
170 if (st%smear%method == smear_fixed_occ) then ! no temperature available
171 hm%energy%TS = m_zero
172 else
173 hm%energy%TS = st%smear%dsmear * hm%energy%entropy
174 end if
175
176 gfield => list_get_gauge_field(ext_partners)
177 if(associated(gfield)) hm%energy%total = hm%energy%total + gauge_field_get_energy(gfield)
178
179 if (allocated(hm%vberry)) then
180 hm%energy%total = hm%energy%total + hm%energy%berry
181 else
182 hm%energy%berry = m_zero
183 end if
184
185 ! Add the magnetic constrain energy
186 if (hm%magnetic_constrain%level /= constrain_none) then
187 hm%energy%total = hm%energy%total + hm%magnetic_constrain%energy
188 end if
189
190 if (optional_default(iunit, -1) /= -1) then
191 write(message(1), '(6x,a, f18.8)')'Total = ', units_from_atomic(units_out%energy, hm%energy%total)
192 write(message(2), '(6x,a, f18.8)')'Free = ', units_from_atomic(units_out%energy, hm%energy%total - hm%energy%TS)
193 write(message(3), '(6x,a)') '-----------'
194 call messages_info(3, iunit)
195
196 write(message(1), '(6x,a, f18.8)')'Ion-ion = ', units_from_atomic(units_out%energy, hm%ep%eii)
197 write(message(2), '(6x,a, f18.8)')'Eigenvalues = ', units_from_atomic(units_out%energy, hm%energy%eigenvalues)
198 write(message(3), '(6x,a, f18.8)')'Hartree = ', units_from_atomic(units_out%energy, hm%energy%hartree)
199 write(message(4), '(6x,a, f18.8)')'Int[n*v_xc] = ', units_from_atomic(units_out%energy, hm%energy%intnvxc)
200 write(message(5), '(6x,a, f18.8)')'Exchange = ', &
201 units_from_atomic(units_out%energy, hm%energy%exchange+hm%energy%exchange_hf)
202 write(message(6), '(6x,a, f18.8)')'Correlation = ', units_from_atomic(units_out%energy, hm%energy%correlation)
203 write(message(7), '(6x,a, f18.8)')'vanderWaals = ', units_from_atomic(units_out%energy, hm%energy%vdw)
204 write(message(8), '(6x,a, f18.8)')'Delta XC = ', units_from_atomic(units_out%energy, hm%energy%delta_xc)
205 write(message(9), '(6x,a, f18.8)')'Entropy = ', hm%energy%entropy ! the dimensionless sigma of Kittel&Kroemer
206 write(message(10), '(6x,a, f18.8)')'-TS = ', -units_from_atomic(units_out%energy, hm%energy%TS)
207 write(message(11), '(6x,a, f18.8)')'Photon ex. = ', units_from_atomic(units_out%energy, hm%energy%photon_exchange)
208 call messages_info(11, iunit)
209
210 if (hm%pcm%run_pcm) then
211 write(message(1),'(6x,a, f18.8)')'E_e-solvent = ', units_from_atomic(units_out%energy, hm%energy%int_ee_pcm + &
212 hm%energy%int_en_pcm + &
213 hm%energy%int_e_ext_pcm)
214 write(message(2),'(6x,a, f18.8)')'E_n-solvent = ', units_from_atomic(units_out%energy, hm%energy%int_nn_pcm + &
215 hm%energy%int_ne_pcm + &
216 hm%energy%int_n_ext_pcm)
217 write(message(3),'(6x,a, f18.8)')'E_M-solvent = ', units_from_atomic(units_out%energy, &
218 hm%energy%int_ee_pcm + hm%energy%int_en_pcm + &
219 hm%energy%int_nn_pcm + hm%energy%int_ne_pcm + &
220 hm%energy%int_e_ext_pcm + hm%energy%int_n_ext_pcm)
221 call messages_info(3, iunit)
222 end if
223
224 if (full_) then
225 write(message(1), '(6x,a, f18.8)')'Kinetic = ', units_from_atomic(units_out%energy, hm%energy%kinetic)
226 write(message(2), '(6x,a, f18.8)')'External = ', units_from_atomic(units_out%energy, hm%energy%extern)
227 write(message(3), '(6x,a, f18.8)')'Non-local = ', units_from_atomic(units_out%energy, hm%energy%extern_non_local)
228 write(message(4), '(6x,a, f18.8)')'Int[n*v_E] = ', units_from_atomic(units_out%energy, hm%energy%intnvstatic)
229 call messages_info(4, iunit)
230 end if
231 if (allocated(hm%vberry) .and. space%is_periodic()) then
232 write(message(1), '(6x,a, f18.8)')'Berry = ', units_from_atomic(units_out%energy, hm%energy%berry)
233 call messages_info(1, iunit)
234 end if
235 if (hm%lda_u_level /= dft_u_none) then
236 write(message(1), '(6x,a, f18.8)')'Hubbard = ', units_from_atomic(units_out%energy, hm%energy%dft_u)
237 write(message(2), '(6x,a, f18.8)')'Int[n*v_U] = ', units_from_atomic(units_out%energy, hm%energy%int_dft_u)
238 call messages_info(2, iunit)
239 end if
240 end if
241
242 pop_sub(energy_calc_total)
243 end subroutine energy_calc_total
244
245 ! --------------------------------------------------------------------
246
247 subroutine energy_calc_eigenvalues(namespace, hm, der, st)
248 type(namespace_t), intent(in) :: namespace
249 type(hamiltonian_elec_t), intent(inout) :: hm
250 type(derivatives_t), intent(in) :: der
251 type(states_elec_t), intent(inout) :: st
252
254
255 if (states_are_real(st)) then
256 call dcalculate_eigenvalues(namespace, hm, der, st)
257 else
258 call zcalculate_eigenvalues(namespace, hm, der, st)
259 end if
260
262 end subroutine energy_calc_eigenvalues
263
264 ! ---------------------------------------------------------
265 subroutine energy_calc_virial_ex(der, vxc, st, ex)
266 type(derivatives_t), intent(in) :: der
267 real(real64), intent(in) :: vxc(:,:)
268 type(states_elec_t), intent(in) :: st
269 real(real64), intent(out) :: ex
270
271 integer :: idir, ip, isp
272 real(real64), allocatable :: gradvx(:,:), nrgradvx(:)
273 real(real64) :: rr, xx(der%dim)
274
275 push_sub(energy_calc_virial_ex)
276
277 assert(st%d%ispin /= spinors)
278
279 safe_allocate(nrgradvx(1:der%mesh%np_part))
280 safe_allocate(gradvx(1:der%mesh%np, 1:der%dim))
281
282 ex = m_zero
283 do isp = 1, st%d%nspin
284 !We output the energy from the virial relation
285 nrgradvx(1:der%mesh%np) = vxc(1:der%mesh%np,isp)
286 call dderivatives_grad(der, nrgradvx, gradvx)
287 nrgradvx = m_zero
288 do idir = 1, der%dim
289 do ip = 1, der%mesh%np
290 call mesh_r(der%mesh, ip, rr, coords=xx)
291 nrgradvx(ip) = nrgradvx(ip) - gradvx(ip, idir) * st%rho(ip, isp) * xx(idir)
292 end do
293 end do
294 ex = ex + dmf_integrate(der%mesh, nrgradvx)
295 end do
296
297 safe_deallocate_a(gradvx)
298 safe_deallocate_a(nrgradvx)
299
300 pop_sub(energy_calc_virial_ex)
301 end subroutine energy_calc_virial_ex
302
303
304#include "undef.F90"
305#include "real.F90"
306#include "energy_calc_inc.F90"
307
308#include "undef.F90"
309#include "complex.F90"
310#include "energy_calc_inc.F90"
311
312end module energy_calc_oct_m
313
314!! Local Variables:
315!! mode: f90
316!! coding: utf-8
317!! End:
This module implements batches of mesh functions.
Definition: batch.F90:133
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public dderivatives_grad(der, ff, op_ff, ghost_update, set_bc, to_cartesian)
apply the gradient to a mesh function
integer, parameter, public spinors
subroutine, public energy_calc_total(namespace, space, hm, gr, st, ext_partners, iunit, full)
This subroutine calculates the total energy of the system. Basically, it adds up the KS eigenvalues,...
real(real64) function, public zenergy_calc_electronic(namespace, hm, der, st, terms)
subroutine dcalculate_eigenvalues(namespace, hm, der, st)
calculates the eigenvalues of the orbitals
subroutine, public energy_calc_virial_ex(der, vxc, st, ex)
subroutine zcalculate_eigenvalues(namespace, hm, der, st)
calculates the eigenvalues of the orbitals
real(real64) function, public denergy_calc_electronic(namespace, hm, der, st, terms)
subroutine, public energy_calc_eigenvalues(namespace, hm, der, st)
subroutine, public zexchange_operator_compute_potentials(this, namespace, space, gr, st, xst, kpoints, ex, F_out)
subroutine, public dexchange_operator_compute_potentials(this, namespace, space, gr, st, xst, kpoints, ex, F_out)
type(gauge_field_t) function, pointer, public list_get_gauge_field(partners)
real(real64) function, public gauge_field_get_energy(this)
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public m_half
Definition: global.F90:194
This module implements the underlying real-space grid.
Definition: grid.F90:117
integer, parameter, public term_local_external
integer, parameter, public term_non_local_potential
integer, parameter, public term_kinetic
This module defines classes and functions for interaction partners.
A module to handle KS potential, without the external potential.
integer, parameter, public hartree
integer, parameter, public hartree_fock
integer, parameter, public independent_particles
integer, parameter, public generalized_kohn_sham_dft
integer, parameter, public kohn_sham_dft
integer, parameter, public dft_u_none
Definition: lda_u.F90:201
This modules implements the routines for doing constrain DFT for noncollinear magnetism.
integer, parameter, public constrain_none
This module defines functions over batches of mesh functions.
Definition: mesh_batch.F90:116
This module defines various routines, operating on mesh functions.
real(real64) function, public dmf_integrate(mesh, ff, mask, reduce)
Integrate a function on the mesh.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
pure subroutine, public mesh_r(mesh, ip, rr, origin, coords)
return the distance to the origin for a given grid point
Definition: mesh.F90:336
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
subroutine, public pcm_elect_energy(ions, pcm, E_int_ee, E_int_en, E_int_ne, E_int_nn, E_int_e_ext, E_int_n_ext)
Calculates the solute-solvent electrostatic interaction energy .
Definition: pcm.F90:1544
real(real64) function, public smear_calc_entropy(this, eigenvalues, nik, nst, kweights, occ)
Definition: smear.F90:593
integer, parameter, public smear_fixed_occ
Definition: smear.F90:171
pure logical function, public states_are_real(st)
This module handles spin dimensions of the states and the k-point distribution.
subroutine, public states_elec_end(st)
finalize the states_elec_t object
real(real64) function, public states_elec_eigenvalues_sum(st, alt_eig)
function to calculate the eigenvalues sum using occupations as weights
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:132
This module defines the unit system, used for input and output.
type(unit_system_t), public units_out
class representing derivatives
The states_elec_t class contains all electronic wave functions.