linear_response.F90

Go to the documentation of this file.

!! Copyright (C) 2004 E.S. Kadantsev, M. Marques
!!
!! 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 linear_response_oct_m
  use batch_oct_m
  use comm_oct_m
  use debug_oct_m
  use electron_space_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use profiling_oct_m
  use restart_oct_m
  use smear_oct_m
  use space_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_calc_oct_m
 
  implicit none
 
  private
 
  public ::            &
    lr_t,              &
    lr_init,           &
    lr_allocate,       &
    lr_zero,           &
    lr_copy,           &
    dlr_orth_vector,   &
    zlr_orth_vector,   &
    dlr_build_dl_rho,  &
    zlr_build_dl_rho,  &
    dlr_orth_response, &
    zlr_orth_response, &
    lr_alpha_j,        &
    lr_dealloc,        &
    lr_is_allocated,   &
    dlr_swap_sigma,    &
    zlr_swap_sigma,    &
    dlr_dump_rho,      &
    zlr_dump_rho,      &
    dlr_load_rho,      &
    zlr_load_rho
 
 
 
  type lr_t
    ! Components are public by default
    logical, private :: is_allocated, is_allocated_rho
 
    real(real64), allocatable :: ddl_rho(:,:)
    real(real64), allocatable :: ddl_psi(:,:,:,:)
 
    complex(real64), allocatable :: zdl_rho(:,:)
    complex(real64), allocatable :: zdl_psi(:,:,:,:)
 
    complex(real64), allocatable :: dl_j(:,:,:)
  end type lr_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine lr_init(lr)
    type(lr_t), intent(out) :: lr
 
    push_sub(lr_init)
 
    lr%is_allocated = .false.
    lr%is_allocated_rho = .false.
 
    pop_sub(lr_init)
 
  end subroutine lr_init
 
 
 
  ! ---------------------------------------------------------
  subroutine lr_allocate(lr, st, mesh, allocate_rho)
    type(lr_t),          intent(inout) :: lr
    type(states_elec_t), intent(in)    :: st
    class(mesh_t),       intent(in)    :: mesh
    logical, optional,   intent(in)    :: allocate_rho
 
    push_sub(lr_allocate)
 
    lr%is_allocated_rho = optional_default(allocate_rho, .true.)
 
    if (states_are_complex(st)) then
      safe_allocate(lr%zdl_psi(1:mesh%np_part, 1:st%d%dim, st%st_start:st%st_end, st%d%kpt%start:st%d%kpt%end))
      if (lr%is_allocated_rho) then
        safe_allocate(lr%zdl_rho(1:mesh%np, 1:st%d%nspin))
      end if
    else
      safe_allocate(lr%ddl_psi(1:mesh%np_part, 1:st%d%dim, st%st_start:st%st_end, st%d%kpt%start:st%d%kpt%end))
      if (lr%is_allocated_rho) then
        safe_allocate(lr%ddl_rho(1:mesh%np, 1:st%d%nspin))
      end if
    end if
 
    lr%is_allocated = .true.
 
    call lr_zero(lr, st)
 
    pop_sub(lr_allocate)
 
  end subroutine lr_allocate
 
 
 
  ! ---------------------------------------------------------
  subroutine lr_zero(lr, st)
    type(lr_t),          intent(inout) :: lr
    type(states_elec_t), intent(in)    :: st
 
    push_sub(lr_zero)
 
    assert(lr%is_allocated)
 
    if (states_are_complex(st)) then
      lr%zdl_psi = m_zero
      if (lr%is_allocated_rho) lr%zdl_rho = m_zero
    else
      lr%ddl_psi = m_zero
      if (lr%is_allocated_rho) lr%ddl_rho = m_zero
    end if
 
    pop_sub(lr_zero)
 
  end subroutine lr_zero
 
 
  ! ---------------------------------------------------------
  subroutine lr_dealloc(lr)
    type(lr_t), intent(inout) :: lr
 
    push_sub(lr_dealloc)
 
    safe_deallocate_a(lr%ddl_psi)
    safe_deallocate_a(lr%zdl_psi)
 
    safe_deallocate_a(lr%ddl_rho)
    safe_deallocate_a(lr%zdl_rho)
 
    safe_deallocate_a(lr%dl_j)
    lr%is_allocated = .false.
    lr%is_allocated_rho = .false.
 
    pop_sub(lr_dealloc)
  end subroutine lr_dealloc
 
 
  ! ---------------------------------------------------------
  subroutine lr_copy(st, mesh, src, dest)
    type(states_elec_t), intent(in) :: st
    class(mesh_t),       intent(in) :: mesh
    type(lr_t),          intent(in) :: src
    type(lr_t),       intent(inout) :: dest
 
    integer :: ik, ist
 
    push_sub(lr_copy)
 
    if (src%is_allocated_rho .and. dest%is_allocated_rho) then
      if (states_are_complex(st)) then
        call lalg_copy(mesh%np, st%d%nspin, src%zdl_rho, dest%zdl_rho)
      else
        call lalg_copy(mesh%np, st%d%nspin, src%ddl_rho, dest%ddl_rho)
      end if
    else
      if (dest%is_allocated_rho) then
        if (states_are_complex(st)) then
          dest%zdl_rho(:, :) = m_zero
        else
          dest%ddl_rho(:, :) = m_zero
        end if
      end if
    end if
 
    do ik = st%d%kpt%start, st%d%kpt%end
      do ist = st%st_start, st%st_end
        if (states_are_complex(st)) then
          call lalg_copy(mesh%np_part, st%d%dim, src%zdl_psi(:, :, ist, ik), dest%zdl_psi(:, :, ist, ik))
        else
          call lalg_copy(mesh%np_part, st%d%dim, src%ddl_psi(:, :, ist, ik), dest%ddl_psi(:, :, ist, ik))
        end if
      end do
    end do
 
    pop_sub(lr_copy)
 
  end subroutine lr_copy
 
 
 
  ! ---------------------------------------------------------
  logical function lr_is_allocated(this)
    type(lr_t), intent(in) :: this
 
    push_sub(lr_is_allocated)
    lr_is_allocated = this%is_allocated
 
    pop_sub(lr_is_allocated)
  end function lr_is_allocated
 
 
 
  ! ---------------------------------------------------------
  real(real64) function lr_alpha_j(st, jst, ik)
    type(states_elec_t), intent(in) :: st
    integer,             intent(in) :: jst
    integer,             intent(in) :: ik
 
    real(real64) :: dsmear
 
    push_sub(lr_alpha_j)
 
    if (st%smear%method == smear_fixed_occ) then
      lr_alpha_j = st%occ(jst, ik) / st%smear%el_per_state
    else
      dsmear = max(1e-14_real64, st%smear%dsmear)
      lr_alpha_j = max(st%smear%e_fermi + m_three*dsmear - st%eigenval(jst, ik), m_zero)
    end if
 
    pop_sub(lr_alpha_j)
  end function lr_alpha_j
 
#include "undef.F90"
#include "real.F90"
#include "linear_response_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "linear_response_inc.F90"
 
end module linear_response_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: