oct_exchange.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 oct_exchange_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use electron_space_oct_m
  use global_oct_m
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use xc_oct_m
 
  implicit none
 
  private
  public ::                          &
    oct_exchange_t,                  &
    oct_exchange_enabled,            &
    oct_exchange_prepare,            &
    oct_exchange_set,                &
    oct_exchange_remove,             &
    oct_exchange_operator
 
  type oct_exchange_t
    private
    logical :: oct_exchange = .false.
    type(states_elec_t), pointer :: oct_st => null()
    real(real64), allocatable :: oct_fxc(:, :, :)
    real(real64), allocatable :: oct_pot(:, :)
    real(real64), allocatable :: oct_rho(:, :)
  end type oct_exchange_t
 
contains
 
  ! ---------------------------------------------------------
  logical function oct_exchange_enabled(this) result(oct_exchange)
    type(oct_exchange_t), intent(in) :: this
 
    push_sub(oct_exchange_enabled)
 
    oct_exchange = this%oct_exchange
 
    pop_sub(oct_exchange_enabled)
  end function oct_exchange_enabled
 
 
  ! ---------------------------------------------------------
  subroutine oct_exchange_set(this, st, mesh)
    type(oct_exchange_t),     intent(inout) :: this
    type(states_elec_t), target, intent(in) :: st
    class(mesh_t),            intent(in)    :: mesh
 
    integer :: np, nspin
 
    push_sub(oct_exchange_set)
 
    ! In this release, no non-local part for the QOCT Hamiltonian.
    nullify(this%oct_st)
 
    this%oct_st => st
    this%oct_exchange = .true.
    np = mesh%np
    nspin = this%oct_st%d%nspin
 
    safe_allocate(this%oct_fxc(1:np, 1:nspin, 1:nspin))
    safe_allocate(this%oct_pot(1:np, 1:nspin))
    safe_allocate(this%oct_rho(1:np, 1:nspin))
 
    this%oct_fxc = m_zero
    this%oct_pot = m_zero
    this%oct_rho = m_zero
 
    pop_sub(oct_exchange_set)
  end subroutine oct_exchange_set
 
 
  ! ---------------------------------------------------------
  subroutine oct_exchange_prepare(this, mesh, psi, xc, psolver, namespace)
    type(oct_exchange_t), intent(inout) :: this
    class(mesh_t),        intent(in)    :: mesh
    complex(real64),      intent(in)    :: psi(:, :, :, :)
    type(xc_t),           intent(in)    :: xc
    type(poisson_t),      intent(in)    :: psolver
    type(namespace_t),    intent(in)    :: namespace
 
    integer :: jst, ip, ik
    complex(real64), allocatable :: psi2(:, :)
 
    push_sub(oct_exchange_prepare)
 
    safe_allocate(psi2(1:mesh%np, 1:this%oct_st%d%dim))
 
    select case (this%oct_st%d%ispin)
    case (unpolarized)
      assert(this%oct_st%nik == 1)
 
      this%oct_pot = m_zero
      this%oct_rho = m_zero
      do jst = 1, this%oct_st%nst
        call states_elec_get_state(this%oct_st, mesh, jst, 1, psi2)
        do ip = 1, mesh%np
          this%oct_rho(ip, 1) = this%oct_rho(ip, 1) + this%oct_st%occ(jst, 1)*aimag(conjg(psi2(ip, 1))*psi(ip, 1, jst, 1))
        end do
      end do
      call dpoisson_solve(psolver, namespace, this%oct_pot(:, 1), this%oct_rho(:, 1), all_nodes = .false.)
 
    case (spin_polarized)
      assert(this%oct_st%nik == 2)
 
      this%oct_pot = m_zero
      this%oct_rho = m_zero
      do ik = 1, 2
        do jst = 1, this%oct_st%nst
          call states_elec_get_state(this%oct_st, mesh, jst, ik, psi2)
          do ip = 1, mesh%np
            this%oct_rho(ip, ik) = this%oct_rho(ip, ik) + this%oct_st%occ(jst, ik) * aimag(conjg(psi2(ip, 1))*psi(ip, 1, jst, ik))
          end do
        end do
      end do
 
      do ik = 1, 2
        call dpoisson_solve(psolver, namespace, this%oct_pot(:, ik), this%oct_rho(:, ik), all_nodes = .false.)
      end do
 
    end select
 
    this%oct_fxc = m_zero
    call xc_get_fxc(xc, mesh, namespace, this%oct_st%rho, this%oct_st%d%ispin, this%oct_fxc)
 
    safe_deallocate_a(psi2)
    pop_sub(oct_exchange_prepare)
  end subroutine oct_exchange_prepare
 
 
  ! ---------------------------------------------------------
  subroutine oct_exchange_remove(this)
    type(oct_exchange_t), intent(inout) :: this
 
    push_sub(oct_exchange_remove)
 
    nullify(this%oct_st)
    this%oct_exchange = .false.
    safe_deallocate_a(this%oct_fxc)
    safe_deallocate_a(this%oct_pot)
    safe_deallocate_a(this%oct_rho)
 
    pop_sub(oct_exchange_remove)
  end subroutine oct_exchange_remove
 
  ! ---------------------------------------------------------
  subroutine oct_exchange_operator(this, namespace, mesh, hpsi, ist, ik)
    type(oct_exchange_t), intent(in)    :: this
    type(namespace_t),    intent(in)    :: namespace
    class(mesh_t),        intent(in)    :: mesh
    complex(real64),      intent(inout) :: hpsi(:, :)
    integer,              intent(in)    :: ist
    integer,              intent(in)    :: ik
 
    integer :: ik2
    complex(real64), allocatable :: psi(:, :), psi2(:, :)
    integer :: ip
 
    push_sub(oct_exchange_operator)
 
    safe_allocate(psi(1:mesh%np, 1:this%oct_st%d%dim))
    safe_allocate(psi2(1:mesh%np, 1:this%oct_st%d%dim))
 
    select case (this%oct_st%d%ispin)
    case (unpolarized)
      assert(this%oct_st%nik == 1)
      call states_elec_get_state(this%oct_st, mesh, ist, 1, psi2)
      do ip = 1, mesh%np
        hpsi(ip, 1) = hpsi(ip, 1) + m_two*m_zi*psi2(ip, 1)*(this%oct_pot(ip, 1) + this%oct_fxc(ip, 1, 1)*this%oct_rho(ip, 1))
      end do
 
    case (spin_polarized)
      assert(this%oct_st%nik == 2)
 
      call states_elec_get_state(this%oct_st, mesh, ist, ik, psi2)
 
      do ik2 = 1, 2
        do ip = 1, mesh%np
          hpsi(ip, 1) = hpsi(ip, 1) + m_two * m_zi * this%oct_st%occ(ist, ik) * &
            psi2(ip, 1) * (this%oct_pot(ip, ik2) + this%oct_fxc(ip, ik, ik2)*this%oct_rho(ip, ik2))
        end do
      end do
 
    case (spinors)
      call messages_not_implemented("Function oct_exchange_operator for spinors", namespace=namespace)
    end select
 
    safe_deallocate_a(psi)
    safe_deallocate_a(psi2)
    pop_sub(oct_exchange_operator)
  end subroutine oct_exchange_operator
 
end module oct_exchange_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: