all_wfcs.f90 Source File


Source Code

!
! Copyright (C) 2024 INTW group
!
! This file is part of INTW.
!
! INTW 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 3 of the License, or
! (at your option) any later version.
!
! INTW 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, see <https://www.gnu.org/licenses/>.
!
module intw_allwfcs

  !! display: none
  !!
  !! This module contains variables and subroutines for obtaining wave functions.
  !!
  !! ### Details
  !!
  !! Wave functions for irreducible k points are read and stored, and wave functions
  !! for general k points are obtained using symmetry.
  !!

  use kinds, only: dp

  implicit none

  ! variables
  public :: eig_all_irr, wfc_k_all_irr, list_iG_all_irr, ngk_all_irr

  ! subroutines
  public :: allocate_and_get_all_irreducible_wfc, get_psi_general_k_all_wfc

  private

  real(dp), allocatable :: eig_all_irr(:,:)
  ! Eigenvalues for each band and irreducible k point
  complex(dp), allocatable :: wfc_k_all_irr(:,:,:,:)
  ! Wave function components (also spin component) for each band and irreducible k point
  integer, allocatable :: list_iG_all_irr(:,:)
  ! G vector indices of the wave function components for each irreducible k point
  integer, allocatable :: ngk_all_irr(:)
  ! Number of components of the wave functions for each irreducible k point

contains

  subroutine allocate_and_get_all_irreducible_wfc()

    use intw_reading, only: nGk_max, nkpoints_QE, get_K_folder_data, nspin, num_bands_intw
    use intw_useful_constants, only: cmplx_0, ZERO

    implicit none

    integer :: list_iG(nGk_max)
    real(dp) :: QE_eig(num_bands_intw)
    complex(dp) :: wfc_g(nGk_max,num_bands_intw,nspin)
    integer :: i_folder, ispin, iG, ibnd


    !
    ! Allocate eigenvalue and wfc related arrays
    !
    if (allocated(eig_all_irr)) deallocate (eig_all_irr)
    allocate(eig_all_irr(nkpoints_QE,num_bands_intw))
    eig_all_irr = ZERO
    !
    if (allocated(wfc_k_all_irr)) deallocate (wfc_k_all_irr)
    allocate(wfc_k_all_irr(nkpoints_QE,nGk_max,num_bands_intw,nspin))
    wfc_k_all_irr = cmplx_0
    !
    if (allocated(list_iG_all_irr)) deallocate(list_iG_all_irr)
    allocate(list_iG_all_irr(nkpoints_QE,nGk_max))
    list_iG_all_irr = 0
    !
    if (allocated(ngk_all_irr)) deallocate(ngk_all_irr)
    allocate(ngk_all_irr(nkpoints_QE))
    ngk_all_irr = 0

    !
    ! Read all irreducible k points
    !
    do i_folder = 1, nkpoints_QE
      !
      call get_K_folder_data(i_folder, list_iG, wfc_g, QE_eig, ngk_all_irr(i_folder))
      !
      eig_all_irr(i_folder,:) = QE_eig(:)
      !
      list_iG_all_irr(i_folder, 1:ngk_all_irr(i_folder)) = list_iG(1:ngk_all_irr(i_folder))
      !
      do iG = 1, ngk_all_irr(i_folder)
        do ispin = 1, nspin
          do ibnd = 1, num_bands_intw
            !
            wfc_k_all_irr(i_folder,iG,ibnd,ispin) = wfc_g(iG,ibnd,ispin)
            !
          enddo ! ibnd
        enddo ! ispin
      enddo ! iG
      !
    enddo ! i_folder

  end subroutine allocate_and_get_all_irreducible_wfc


  subroutine get_psi_general_k_all_wfc(kpoint, ngk, list_iG_k, wfc_k, eig_k)

    use intw_reading, only: s, ftau, nGk_max, nspin, kpoints_QE, num_bands_intw
    use intw_input_parameters, only: nk1, nk2, nk3
    use intw_symmetries, only: full_mesh, symlink, QE_folder_sym, QE_folder_nosym, &
                               apply_TR_to_wfc, rotate_wfc
    use intw_utility, only: find_k_1BZ_and_G, triple_to_joint_index_g
    use intw_fft, only: wfc_by_expigr
    use intw_useful_constants, only : eps_5, ZERO

    implicit none

    !I/O variables
    real(dp), intent(in) :: kpoint(3) ! General k point in crystal coordinates for which wave functions are desired
    integer, intent(out) :: ngk ! Number of components of the wave functions
    integer, intent(out) :: list_iG_k(nGk_max) ! G vector indices of the wave function components
    complex(dp), intent(out) :: wfc_k(nGk_max,num_bands_intw,nspin) ! Wave function components (also spin component) for all bands
    real(dp), intent(out), optional :: eig_k(num_bands_intw) ! Eigenvalues for all band

    !local variables
    integer :: ikpt, i_folder
    integer :: i_sym, TR
    integer :: i_1bz, j_1bz, k_1bz
    integer :: sym(3,3), G_sym(3), G_plus(3)
    real(dp) :: ftau_sym(3)
    real(dp) :: kpoint_1BZ(3), k_irr(3), ktest(3)
    integer :: list_iG_irr(nGk_max)
    complex(dp) :: wfc_k_irr(nGk_max,num_bands_intw,nspin)


    call find_k_1BZ_and_G(kpoint, nk1, nk2, nk3, i_1bz, j_1bz, k_1bz, kpoint_1bz, G_plus)
    !
    call triple_to_joint_index_g(nk1, nk2, nk3, ikpt, i_1bz, j_1bz, k_1bz)
    !
    if (full_mesh) then
      !
      ! Use the full BZ, no symmetry!
      !
      i_folder = QE_folder_nosym(ikpt)
      k_irr = kpoints_QE(:,i_folder)
      !
      ngk = ngk_all_irr(ikpt)
      list_iG_k(:) = list_iG_all_irr(ikpt,:)
      wfc_k(:,:,:) = wfc_k_all_irr(ikpt,:,:,:)
      if (present(eig_k)) eig_k(:) = eig_all_irr(ikpt,:)
      !
      G_sym = nint(kpoint - k_irr)
      !
      call wfc_by_expigr(num_bands_intw, nspin, G_sym, list_iG_k, wfc_k)
      !
    else
      !
      ! Use the IBZ and symmetry!
      !
      ! This is the irreducible point
      ! connected to aimed kpoint
      i_folder = QE_folder_sym(ikpt)
      k_irr = kpoints_QE(:,i_folder)
      !
      ! The symmetry which takes kpoints_QE(:,i_folder) into aimed kpoint.
      ! sym * kpoints_QE = kpoint
      i_sym = symlink(ikpt,1)
      TR = symlink(ikpt,2)
      ftau_sym = ftau(:,(i_sym))
      sym = s(:,:,(i_sym))
      !
      ! Load the corresponding irreducible wfcs in kpoints_QE
      ngk = ngk_all_irr(i_folder)
      list_iG_irr(:) = list_iG_all_irr(i_folder,:)
      wfc_k_irr(:,:,:) = wfc_k_all_irr(i_folder,:,:,:)
      if (present(eig_k)) eig_k(:) = eig_all_irr(i_folder,:)
      !
      ktest = matmul(sym ,k_irr)
      !
      if (TR==1) then
        ! If TR needed
        ! G_sym = aimed_kpoint -TR[S*QE_kpoint] = aimed_kpoint + S*QE_kpoint, such that
        ! aimed_kpoint = -S*QE_kpoint + G_ sym
        G_sym = nint(kpoint + ktest)
      else
        ! G_sym = aimed_kpoint - S*QE_kpoint, such that
        ! aimed_kpoint = S*QE_kpoint + G_ sym
        G_sym = nint(kpoint - ktest)
      end if
      !
      call rotate_wfc(wfc_k_irr, list_iG_irr, wfc_k, list_iG_k, i_sym, sym, ftau_sym, (/0, 0, 0/))
      !
      ! If time-reversal is present, the wavefunction currently stored
      ! in wfc_k is actually for (-k). Complex conjugation must now
      ! be applied to recover wfc_k.
      !
      if (TR==1) then
        ktest = -ktest
        !
        call apply_TR_to_wfc(wfc_k, list_iG_k)
        !
      endif
      !
      if ( sum( abs( ktest + dble(G_sym) - kpoint ) ) > eps_5 ) then
          write(*,*)            "ERROR in get_psi_general_k_all_wfc:"
          write(*,"(a,3f12.6)") "Aimed kpoint is     :", kpoint
          write(*,"(a,3f12.6)") "Symmetry induced is :", ktest + G_sym
      end if
      !
      call wfc_by_expigr(num_bands_intw, nspin, G_sym, list_iG_k, wfc_k)
      !
    endif

  end subroutine get_psi_general_k_all_wfc

end module intw_allwfcs