tetrahedron_integral Function

public function tetrahedron_integral(F_in, D_in, t_in, hw, type, tet_cutoff, avoid_deg)

Uses

  • proc~~tetrahedron_integral~~UsesGraph proc~tetrahedron_integral tetrahedron_integral module~w90_constants w90_constants proc~tetrahedron_integral->module~w90_constants module~w90_utility w90_utility proc~tetrahedron_integral->module~w90_utility module~w90_utility->module~w90_constants module~w90_comms w90_comms module~w90_utility->module~w90_comms module~w90_comms->module~w90_constants module~w90_error_base w90_error_base module~w90_comms->module~w90_error_base

triangle 1

triangle 2

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), dimension(4) :: F_in
real(kind=dp), intent(in), dimension(4) :: D_in
real(kind=dp), intent(in), dimension(3, 3) :: t_in
real(kind=dp), intent(in) :: hw
integer, intent(in) :: type
real(kind=dp), intent(in) :: tet_cutoff
real(kind=dp), intent(in) :: avoid_deg

Return Value real(kind=dp)


Calls

proc~~tetrahedron_integral~~CallsGraph proc~tetrahedron_integral tetrahedron_integral proc~tetrahedron_jacobian tetrahedron_jacobian proc~tetrahedron_integral->proc~tetrahedron_jacobian proc~tetrahedron_log1p tetrahedron_log1p proc~tetrahedron_integral->proc~tetrahedron_log1p proc~tetrahedron_sort tetrahedron_sort proc~tetrahedron_integral->proc~tetrahedron_sort proc~utility_inv3 utility_inv3 proc~tetrahedron_integral->proc~utility_inv3

Called by

proc~~tetrahedron_integral~~CalledByGraph proc~tetrahedron_integral tetrahedron_integral proc~tetrahedron_fermidirac tetrahedron_fermidirac proc~tetrahedron_fermidirac->proc~tetrahedron_integral proc~tetrahedron_spinhall tetrahedron_spinhall proc~tetrahedron_spinhall->proc~tetrahedron_fermidirac proc~berry_main berry_main proc~berry_main->proc~tetrahedron_spinhall program~postw90 postw90 program~postw90->proc~berry_main

Source Code

  function tetrahedron_integral(F_in, D_in, t_in, hw, type, tet_cutoff, avoid_deg)
    !=============================================================================
    ! Calculates contribution from a single tetrahedron, for the Kubo formula.   !
    ! To do: extend to complex numbers                                           !
    ! type 1: nondissipative, linear/linear                                      !
    ! type 2: dissipative                                                        !
    ! type 3: nondissipative, linear/(linear)**2                                 !
    !                                                                            !
    ! Ref: Minsu Ghim and Cheol-Hwan Park, PRB 106 075126                        !
    !=============================================================================

    use w90_constants, only: dp
    use w90_utility, only: utility_inv3

    implicit none

    real(kind=dp) :: tetrahedron_integral
    real(kind=dp), dimension(4), intent(in) :: F_in, D_in
    real(kind=dp), dimension(3, 3), intent(in) :: t_in
    real(kind=dp), intent(in) :: hw
    integer, intent(in) :: type
    real(kind=dp), intent(in) :: tet_cutoff
    !small parameter to avoid a problem of degenearcy
    real(kind=dp), intent(in) :: avoid_deg
    !result
    real(kind=dp) :: Ans, Det_t

    !intermediate variables
    integer :: i, j, k, l
    real(kind=dp) :: D(4), F(4), t(3, 3), DAV, D_small_prev, D_large_prev
    real(kind=dp) :: dd(3), ll(3), ff, bb(4), cc(4, 3)
    integer :: a, b, c
    REAL(kind=dp) :: t_inverse(3, 3), dummy(4)
    REAL(kind=dp) :: GradD, Jac, y
    REAL(kind=dp) :: x(3)
    REAL(kind=dp), DIMENSION(0:2) :: F_uv

    D = D_in; F = F_in; t = t_in
    call tetrahedron_sort(D, F, dummy, t)
    Ans = 0.0_dp

    !case 1 and 3: nondissipative part, case 2: dissipative part
    if (type == 1 .or. type == 3) then
      if (type == 3) then
        !treatment for accidental small band spliting(but degenerate actually) cases
        do j = 1, 4
          if (abs(D(j)) < avoid_deg) then
            D(j) = avoid_deg*(abs(D(j))/D(j))
            F = 0
          end if
        end do
      end if

      !cutoff treatment, hw == 0.0 for case 3
      DAV = (D(2) + D(3))/2.0_dp
      if (abs((D(2) - D(3))/(DAV + hw)) < tet_cutoff) then
        D_small_prev = D(2); D_large_prev = D(3)
        D(3) = DAV + 0.5_dp*abs(DAV + hw)*tet_cutoff
        D(2) = DAV - 0.5_dp*abs(DAV + hw)*tet_cutoff
        if (D(1) > D(2)) D(1) = D(1) + (D(2) - D_small_prev)
        if (D(3) > D(4)) D(4) = D(4) + (D(3) - D_large_prev)
      end if
      DAV = (D(1) + D(2))/2.0_dp
      if (abs((D(1) - D(2))/(DAV + hw)) < tet_cutoff) then
        if (D(2) > 0) then
          D(1) = D(2)*(2.0_dp - tet_cutoff)/(2.0_dp + tet_cutoff)
        else
          D(1) = D(2)*(2.0_dp + tet_cutoff)/(2.0_dp - tet_cutoff)
        end if
      end if
      DAV = (D(3) + D(4))/2.0_dp
      if (abs((D(3) - D(4))/(DAV + hw)) < tet_cutoff) then
        if (D(3) > 0) then
          D(4) = D(3)*(2.0_dp + tet_cutoff)/(2.0_dp - tet_cutoff)
        else
          D(4) = D(3)*(2.0_dp - tet_cutoff)/(2.0_dp + tet_cutoff)
        end if
      end if

      !intermediate variables for case 1 and 3
      do i = 1, 3
        dd(i) = (D(4) - D(i))/(D(i) + hw)
        ll(i) = tetrahedron_log1p(dd(i))
      end do
      ff = 1.0_dp

      !determinant factor from parametrisation(tetrahedron volume)
      Det_t = ABS(t(1, 1)*t(2, 2)*t(3, 3) + t(1, 2)*t(2, 3)*t(3, 1) + t(1, 3)*t(2, 1)*t(3, 2) &
                  - t(1, 1)*t(2, 3)*t(3, 2) - t(1, 3)*t(2, 2)*t(3, 1) - t(1, 2)*t(2, 1)*t(3, 3))
    end if

    select case (type)
    case (1)
      do i = 1, 3
        a = i
        b = mod(i, 3) + 1
        c = mod(i + 1, 3) + 1
        cc(a, a) = -(1.0_dp + dd(a))*(3.0_dp*dd(a)**2 - 2.0_dp*(dd(b) + dd(c))*dd(a) + dd(b)*dd(c)) &
                   *((dd(b) - dd(c))*dd(b)*dd(c))**2
        cc(b, a) = -dd(a)*(1.0_dp + dd(b))*(dd(c) - dd(a))*((dd(b) - dd(c))*dd(b)*dd(c))**2
        cc(c, a) = dd(a)*(1.0_dp + dd(c))*(dd(a) - dd(b))*((dd(b) - dd(c))*dd(b)*dd(c))**2
        cc(4, a) = -(dd(a) - dd(b))*(dd(c) - dd(a))*((dd(b) - dd(c))*dd(b)*dd(c))**2
        bb(a) = cc(4, a)*dd(a)
        ff = ff*(1.0_dp + dd(a))/(dd(a)*(dd(a) - dd(b)))**2
      end do
      bb(4) = -dd(1)*dd(2)*dd(3)*((dd(1) - dd(2))*(dd(2) - dd(3))*(dd(3) - dd(1)))**2
      ff = -ff/6.0_dp

      do i = 1, 4
        Ans = Ans + F(i)*(cc(i, 1)*ll(1) + cc(i, 2)*ll(2) + cc(i, 3)*ll(3) + bb(i))
      end do
      Ans = Ans*ff/(D(4) + hw)

      tetrahedron_integral = Ans*Det_t

    case (2)
      ! calculate integ d2k / |grad D| * F(k)

      ! obtaining |grad D|
      call utility_inv3(t, t_inverse, Det_t)
      t_inverse = t_inverse/Det_t

      GradD = 0.0_dp
      DO i = 1, 3
        DO j = 1, 3
          DO k = 1, 3
            GradD = GradD + t_inverse(i, k)*t_inverse(j, k)*(D(i + 1) - D(1))*(D(j + 1) - D(1))
          END DO
        END DO
      END DO
      GradD = SQRT(ABS(GradD))

      IF (hw < D(1)) THEN
        Ans = 0.0_dp
      ELSE IF (hw < D(2)) THEN
        ! parametrization
        x(1) = (hw - D(1))/(D(2) - D(1))
        x(2) = (hw - D(1))/(D(3) - D(1))
        x(3) = (hw - D(1))/(D(4) - D(1))
        ! Jacobian factor
        Jac = tetrahedron_jacobian(t, x, 1)
        ! integration formula
        F_uv(0) = F(1) + (F(2) - F(1))*x(1)
        F_uv(1) = (F(3) - F(1))*x(2) - (F(2) - F(1))*x(1)
        F_uv(2) = (F(4) - F(1))*x(3) - (F(2) - F(1))*x(1)
        Ans = Jac*(F_uv(0)/2.0_dp + (F_uv(1) + F_uv(2))/6.0_dp)/GradD
      ELSE IF (hw < D(3)) THEN
        ! parametrization
        x(1) = (hw - D(4))/(D(2) - D(4))
        x(2) = (hw - D(1))/(D(3) - D(1))
        x(3) = (hw - D(1))/(D(4) - D(1))

        !! triangle 1
        ! Jacobian factor
        Jac = tetrahedron_jacobian(t, x, 2)
        ! integration formula
        F_uv(0) = F(1) + (F(4) - F(1))*x(3)
        F_uv(1) = (F(2) - F(1))*x(1) + (F(4) - F(1))*(1.0_dp - x(1) - x(3))
        F_uv(2) = (F(3) - F(1))*x(2) - (F(4) - F(1))*x(3)
        Ans = Jac*(F_uv(0)/2.0_dp + (F_uv(1) + F_uv(2))/6.0_dp)

        !! triangle 2
        y = (hw - D(3))/(D(2) - D(3))
        ! Jacobian factor
        Jac = tetrahedron_jacobian(t, x, 3)
        ! integration formula
        F_uv(0) = F(1) + (F(2) - F(1))*y + (F(3) - F(1))*(1.0_dp - y)
        F_uv(1) = (F(2) - F(1))*(x(1) - y) + &
                  (F(3) - F(1))*(y - 1.0_dp) + (F(4) - F(1))*(1.0_dp - x(1))
        F_uv(2) = -(F(2) - F(1))*y + (F(3) - F(1))*(y - 1.0_dp + x(2))
        Ans = Ans + Jac*(F_uv(0)/2.0_dp + (F_uv(1) + F_uv(2))/6.0_dp)

        Ans = Ans/GradD
      ELSE IF (hw < D(4)) THEN
        ! parametrization
        x(1) = (hw - D(4))/(D(2) - D(4))
        x(2) = (hw - D(4))/(D(3) - D(4))
        x(3) = (hw - D(1))/(D(4) - D(1))
        ! Jacobian factor
        Jac = tetrahedron_jacobian(t, x, 4)
        ! integration formula
        F_uv(0) = F(1) + (F(4) - F(1))*x(3)
        F_uv(1) = (F(2) - F(1))*x(1) + (F(4) - F(1))*(1.0_dp - x(1) - x(3))
        F_uv(2) = (F(3) - F(1))*x(2) + (F(4) - F(1))*(1.0_dp - x(2) - x(3))
        Ans = Jac*(F_uv(0)/2.0_dp + (F_uv(1) + F_uv(2))/6.0_dp)/GradD
      ELSE
        Ans = 0.0_dp
      END IF
      tetrahedron_integral = Ans

    case (3)
      do i = 1, 3
        a = i
        b = mod(i, 3) + 1
        c = mod(i + 1, 3) + 1
        cc(a, a) = -(1.0_dp + dd(a))*(2.0_dp*dd(a)**3 + (3.0_dp - dd(b) - dd(c))*dd(a)**2 &
                                      - 2.0_dp*(dd(b) + dd(c))*dd(a) + dd(b)*dd(c))*((dd(b) - dd(c))*dd(b)*dd(c))**2
        cc(b, a) = -(1.0_dp + dd(a))*dd(a)*(1.0_dp + dd(b))*(dd(c) - dd(a))*((dd(b) - dd(c))*dd(b)*dd(c))**2
        cc(c, a) = (1.0_dp + dd(a))*dd(a)*(1.0_dp + dd(c))*(dd(a) - dd(b))*((dd(b) - dd(c))*dd(b)*dd(c))**2
        cc(4, a) = -(1.0_dp + dd(a))*(dd(a) - dd(b))*(dd(c) - dd(a))*((dd(b) - dd(c))*dd(b)*dd(c))**2
        bb(a) = cc(4, a)*dd(a)
        ff = ff*(1.0_dp + dd(a))/(dd(a)*(dd(a) - dd(b)))**2
      end do
      bb(4) = -dd(1)*dd(2)*dd(3)*((dd(1) - dd(2))*(dd(2) - dd(3))*(dd(3) - dd(1)))**2
      ff = ff/2.0_dp
      do i = 1, 4
        Ans = Ans + F(i)*(cc(i, 1)*ll(1) + cc(i, 2)*ll(2) + cc(i, 3)*ll(3) + bb(i))
      end do
      Ans = Ans*ff/(D(4) + hw)**2

      tetrahedron_integral = Ans*Det_t
    case default
      tetrahedron_integral = 0.0_dp
    end select

    return
  end function tetrahedron_integral