energy_calc.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module energy_calc_oct_m
  use batch_oct_m
  use comm_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use energy_oct_m
  use electron_space_oct_m
  use exchange_operator_oct_m
  use ext_partner_list_oct_m
  use gauge_field_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use hamiltonian_elec_base_oct_m
  use ions_oct_m
  use ks_potential_oct_m
  use lda_u_oct_m
  use magnetic_constrain_oct_m
  use mesh_oct_m
  use mesh_batch_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use namespace_oct_m
  use interaction_partner_oct_m
  use pcm_oct_m
  use profiling_oct_m
  use smear_oct_m
  use space_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use wfs_elec_oct_m
  use xc_oct_m
 
  implicit none
 
  private
  public ::                               &
    energy_calc_total,                    &
    denergy_calc_electronic,              &
    zenergy_calc_electronic,              &
    energy_calc_eigenvalues,              &
    energy_calc_virial_ex
 
contains
 
  ! ---------------------------------------------------------
  subroutine energy_calc_total(namespace, space, hm, gr, st, ext_partners, iunit, full)
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(grid_t),             intent(in)    :: gr
    type(states_elec_t),      intent(inout) :: st
    type(partner_list_t),     intent(in)    :: ext_partners
    integer, optional,        intent(in)    :: iunit
    logical, optional,        intent(in)    :: full
 
    logical :: full_
    type(states_elec_t) :: xst
    type(gauge_field_t), pointer :: gfield
 
    push_sub(energy_calc_total)
 
    full_ = .false.
    if (present(full)) full_ = full
 
    hm%energy%eigenvalues = states_elec_eigenvalues_sum(st)
 
    if (full_ .or. hm%theory_level == hartree .or. hm%theory_level == hartree_fock &
      .or. hm%theory_level == generalized_kohn_sham_dft) then
      if (states_are_real(st)) then
        hm%energy%kinetic  = denergy_calc_electronic(namespace, hm, gr%der, st, terms = term_kinetic)
        hm%energy%extern_local = denergy_calc_electronic(namespace, hm, gr%der, st, terms = term_local_external)
        hm%energy%extern_non_local   = denergy_calc_electronic(namespace, hm, gr%der, st, terms = term_non_local_potential)
        hm%energy%extern = hm%energy%extern_local + hm%energy%extern_non_local
      else
        hm%energy%kinetic = zenergy_calc_electronic(namespace, hm, gr%der, st, terms = term_kinetic)
 
        hm%energy%extern_local = zenergy_calc_electronic(namespace, hm, gr%der, st, terms = term_local_external)
 
        hm%energy%extern_non_local = zenergy_calc_electronic(namespace, hm, gr%der, st, terms = term_non_local_potential)
        hm%energy%extern   = hm%energy%extern_local   + hm%energy%extern_non_local
      end if
    end if
 
    ! Computes the HF energy when not doing the ACE
    if((hm%theory_level == hartree_fock .or. &
      (hm%theory_level == generalized_kohn_sham_dft .and. family_is_hybrid(hm%xc))) ) then
      if (.not. hm%exxop%useACE) then
        call xst%nullify()
        if (states_are_real(st)) then
          call dexchange_operator_compute_potentials(hm%exxop, namespace, space, gr, st, &
            xst, hm%kpoints)
          hm%energy%exchange_hf = dexchange_operator_compute_ex(gr, st, xst)
        else
          call zexchange_operator_compute_potentials(hm%exxop, namespace, space, gr, st, &
            xst, hm%kpoints)
          hm%energy%exchange_hf = zexchange_operator_compute_ex(gr, st, xst)
        end if
        call states_elec_end(xst)
        hm%energy%exchange_hf = hm%energy%exchange_hf + hm%exxop%singul%energy
      end if
      ! If (useACE) is handled elsewhere
    else
      ! If exact exchange is not required, zero the term regardless
      hm%energy%exchange_hf = m_zero
    end if
 
    if (hm%pcm%run_pcm) then
      hm%pcm%counter = hm%pcm%counter + 1
      if (hm%pcm%localf) then
        call pcm_elect_energy(hm%ions, hm%pcm, hm%energy%int_ee_pcm, hm%energy%int_en_pcm, &
          hm%energy%int_ne_pcm, hm%energy%int_nn_pcm, &
          e_int_e_ext = hm%energy%int_e_ext_pcm,      &
          e_int_n_ext = hm%energy%int_n_ext_pcm       )
      else
        call pcm_elect_energy(hm%ions, hm%pcm, hm%energy%int_ee_pcm, hm%energy%int_en_pcm, &
          hm%energy%int_ne_pcm, hm%energy%int_nn_pcm  )
      end if
    end if
 
    select case (hm%theory_level)
    case (independent_particles)
      hm%energy%total = hm%ep%eii + hm%energy%eigenvalues
 
    case (hartree)
      hm%energy%total = hm%ep%eii + m_half * (hm%energy%eigenvalues + hm%energy%kinetic + hm%energy%extern)
 
    case (hartree_fock, generalized_kohn_sham_dft)
      hm%energy%total = hm%ep%eii + &
        m_half*(hm%energy%eigenvalues + hm%energy%kinetic + hm%energy%extern - hm%energy%intnvxc &
        - hm%energy%int_dft_u) &
        + hm%energy%exchange + hm%energy%exchange_hf + hm%energy%correlation &
        + hm%energy%vdw - hm%energy%intnvstatic + hm%energy%dft_u
 
      ! FIXME: pcm terms are only added to total energy in DFT case
 
    case (kohn_sham_dft)
      hm%energy%total = hm%ep%eii + hm%energy%eigenvalues &
        - hm%energy%hartree + hm%energy%exchange + hm%energy%correlation + hm%energy%vdw - hm%energy%intnvxc &
        - hm%energy%pcm_corr + hm%energy%int_ee_pcm + hm%energy%int_en_pcm &
        + hm%energy%int_nn_pcm + hm%energy%int_ne_pcm &
        + hm%energy%int_e_ext_pcm + hm%energy%int_n_ext_pcm &
        + hm%energy%dft_u -  hm%energy%int_dft_u - hm%energy%intnvstatic &
        + hm%energy%photon_exchange
 
    case default
      hm%energy%total = m_zero
    end select
 
    hm%energy%entropy = smear_calc_entropy(st%smear, st%eigenval, st%nik, st%nst, st%kweights, st%occ)
    if (st%smear%method == smear_fixed_occ) then ! no temperature available
      hm%energy%TS = m_zero
    else
      hm%energy%TS = st%smear%dsmear * hm%energy%entropy
    end if
 
    gfield => list_get_gauge_field(ext_partners)
    if(associated(gfield)) hm%energy%total = hm%energy%total + gauge_field_get_energy(gfield)
 
    if (allocated(hm%vberry)) then
      hm%energy%total = hm%energy%total + hm%energy%berry
    else
      hm%energy%berry = m_zero
    end if
 
    ! Add the magnetic constrain energy
    if (hm%magnetic_constrain%level /= constrain_none) then
      hm%energy%total = hm%energy%total + hm%magnetic_constrain%energy
    end if
 
    if (optional_default(iunit, -1) /= -1) then
      write(message(1), '(6x,a, f18.8)')'Total       = ', units_from_atomic(units_out%energy, hm%energy%total)
      write(message(2), '(6x,a, f18.8)')'Free        = ', units_from_atomic(units_out%energy, hm%energy%total - hm%energy%TS)
      write(message(3), '(6x,a)') '-----------'
      call messages_info(3, iunit)
 
      write(message(1), '(6x,a, f18.8)')'Ion-ion     = ', units_from_atomic(units_out%energy, hm%ep%eii)
      write(message(2), '(6x,a, f18.8)')'Eigenvalues = ', units_from_atomic(units_out%energy, hm%energy%eigenvalues)
      write(message(3), '(6x,a, f18.8)')'Hartree     = ', units_from_atomic(units_out%energy, hm%energy%hartree)
      write(message(4), '(6x,a, f18.8)')'Int[n*v_xc] = ', units_from_atomic(units_out%energy, hm%energy%intnvxc)
      write(message(5), '(6x,a, f18.8)')'Exchange    = ', &
        units_from_atomic(units_out%energy, hm%energy%exchange+hm%energy%exchange_hf)
      write(message(6), '(6x,a, f18.8)')'Correlation = ', units_from_atomic(units_out%energy, hm%energy%correlation)
      write(message(7), '(6x,a, f18.8)')'vanderWaals = ', units_from_atomic(units_out%energy, hm%energy%vdw)
      write(message(8), '(6x,a, f18.8)')'Delta XC    = ', units_from_atomic(units_out%energy, hm%energy%delta_xc)
      write(message(9), '(6x,a, f18.8)')'Entropy     = ', hm%energy%entropy ! the dimensionless sigma of Kittel&Kroemer
      write(message(10), '(6x,a, f18.8)')'-TS         = ', -units_from_atomic(units_out%energy, hm%energy%TS)
      write(message(11), '(6x,a, f18.8)')'Photon ex.  = ', units_from_atomic(units_out%energy, hm%energy%photon_exchange)
      call messages_info(11, iunit)
 
      if (hm%pcm%run_pcm) then
        write(message(1),'(6x,a, f18.8)')'E_e-solvent = ',  units_from_atomic(units_out%energy, hm%energy%int_ee_pcm + &
          hm%energy%int_en_pcm + &
          hm%energy%int_e_ext_pcm)
        write(message(2),'(6x,a, f18.8)')'E_n-solvent = ',  units_from_atomic(units_out%energy, hm%energy%int_nn_pcm + &
          hm%energy%int_ne_pcm + &
          hm%energy%int_n_ext_pcm)
        write(message(3),'(6x,a, f18.8)')'E_M-solvent = ',  units_from_atomic(units_out%energy, &
          hm%energy%int_ee_pcm + hm%energy%int_en_pcm + &
          hm%energy%int_nn_pcm + hm%energy%int_ne_pcm + &
          hm%energy%int_e_ext_pcm + hm%energy%int_n_ext_pcm)
        call messages_info(3, iunit)
      end if
 
      if (full_) then
        write(message(1), '(6x,a, f18.8)')'Kinetic     = ', units_from_atomic(units_out%energy, hm%energy%kinetic)
        write(message(2), '(6x,a, f18.8)')'External    = ', units_from_atomic(units_out%energy, hm%energy%extern)
        write(message(3), '(6x,a, f18.8)')'Non-local   = ', units_from_atomic(units_out%energy, hm%energy%extern_non_local)
        write(message(4), '(6x,a, f18.8)')'Int[n*v_E]  = ', units_from_atomic(units_out%energy, hm%energy%intnvstatic)
        call messages_info(4, iunit)
      end if
      if (allocated(hm%vberry) .and. space%is_periodic()) then
        write(message(1), '(6x,a, f18.8)')'Berry       = ', units_from_atomic(units_out%energy, hm%energy%berry)
        call messages_info(1, iunit)
      end if
      if (hm%lda_u_level /= dft_u_none) then
        write(message(1), '(6x,a, f18.8)')'Hubbard     = ', units_from_atomic(units_out%energy, hm%energy%dft_u)
        write(message(2), '(6x,a, f18.8)')'Int[n*v_U]  = ', units_from_atomic(units_out%energy, hm%energy%int_dft_u)
        call messages_info(2, iunit)
      end if
    end if
 
    pop_sub(energy_calc_total)
  end subroutine energy_calc_total
 
  ! --------------------------------------------------------------------
 
  subroutine energy_calc_eigenvalues(namespace, hm, der, st)
    type(namespace_t),        intent(in)    :: namespace
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(derivatives_t),      intent(in)    :: der
    type(states_elec_t),      intent(inout) :: st
 
    push_sub(energy_calc_eigenvalues)
 
    if (states_are_real(st)) then
      call dcalculate_eigenvalues(namespace, hm, der, st)
    else
      call zcalculate_eigenvalues(namespace, hm, der, st)
    end if
 
    pop_sub(energy_calc_eigenvalues)
  end subroutine energy_calc_eigenvalues
 
  ! ---------------------------------------------------------
  subroutine energy_calc_virial_ex(der, vxc, st, ex)
    type(derivatives_t),      intent(in) :: der
    real(real64),             intent(in) :: vxc(:,:)
    type(states_elec_t),      intent(in) :: st
    real(real64),            intent(out) :: ex
 
    integer :: idir, ip, isp
    real(real64), allocatable :: gradvx(:,:), nrgradvx(:)
    real(real64) :: rr, xx(der%dim)
 
    push_sub(energy_calc_virial_ex)
 
    assert(st%d%ispin /= spinors)
 
    safe_allocate(nrgradvx(1:der%mesh%np_part))
    safe_allocate(gradvx(1:der%mesh%np, 1:der%dim))
 
    ex = m_zero
    do isp = 1, st%d%nspin
      !We output the energy from the virial relation
      nrgradvx(1:der%mesh%np) = vxc(1:der%mesh%np,isp)
      call dderivatives_grad(der, nrgradvx, gradvx)
      nrgradvx = m_zero
      do idir = 1, der%dim
        do ip = 1, der%mesh%np
          call mesh_r(der%mesh, ip, rr, coords=xx)
          nrgradvx(ip) = nrgradvx(ip) - gradvx(ip, idir) * st%rho(ip, isp) * xx(idir)
        end do
      end do
      ex = ex + dmf_integrate(der%mesh, nrgradvx)
    end do
 
    safe_deallocate_a(gradvx)
    safe_deallocate_a(nrgradvx)
 
    pop_sub(energy_calc_virial_ex)
  end subroutine energy_calc_virial_ex
 
 
#include "undef.F90"
#include "real.F90"
#include "energy_calc_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "energy_calc_inc.F90"
 
end module energy_calc_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: