Line data    Source code

       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Antisymmetric DC charge response with a scalar relaxation rate.
      10              : !> \author Thomas D. Kuehne, tkuehne@cp2k.org
      11              : ! **************************************************************************************************
      12              : MODULE kubo_hall
      13              :    USE ieee_arithmetic,                 ONLY: ieee_is_finite
      14              :    USE kinds,                           ONLY: dp
      15              :    USE kubo_bloch,                      ONLY: fermi_slope
      16              : 
      17              :    IMPLICIT NONE
      18              :    PRIVATE
      19              :    PUBLIC :: hall_response
      20              : 
      21              :    REAL(dp), PARAMETER :: antihermiticity_tolerance = 1.e-9_dp
      22              : CONTAINS
      23              : ! **************************************************************************************************
      24              : !> \brief Hall tensor, with first index the current and second index the electric field.
      25              : !> \param energies Complete eigenvalue list for one spin channel and k point
      26              : !> \param current Anti-Hermitian -i*hbar*v matrix triplet, in exactly the same eigenframe
      27              : !> \param mu Chemical potential
      28              : !> \param kt Positive temperature in energy units
      29              : !> \param gamma Positive relaxation rate in energy units
      30              : !> \param hall Antisymmetric tensor before e^2/hbar, volume, spin and k weights
      31              : !> \param status Zero on success; negative for invalid input or nonfinite arithmetic
      32              : !> \note The trace uses (gamma+i[H,.])^-1. Pairing n,m gives
      33              : !>       2*f[en,em]*(en-em)*Im[J_a(m,n)*J_b(n,m)]/(gamma^2+(en-em)^2).
      34              : !>       Equal-energy blocks contribute zero without a band-gap division.
      35              : ! **************************************************************************************************
      36        27750 :    SUBROUTINE hall_response(energies, current, mu, kt, gamma, hall, status)
      37              :       REAL(dp), INTENT(IN)                               :: energies(:)
      38              :       COMPLEX(dp), INTENT(IN)                            :: current(:, :, :)
      39              :       REAL(dp), INTENT(IN)                               :: mu, kt, gamma
      40              :       REAL(dp), INTENT(OUT)                              :: hall(3, 3)
      41              :       INTEGER, INTENT(OUT)                               :: status
      42              : 
      43              :       INTEGER                                            :: a, b, m, n, nstate
      44              :       REAL(dp)                                           :: delta, factor, scale, total(3, 3)
      45              : 
      46        27750 :       hall = 0.0_dp
      47        27750 :       status = -1
      48        27750 :       nstate = SIZE(energies)
      49       111000 :       IF (nstate < 1 .OR. ANY(SHAPE(current) /= [nstate, nstate, 3])) RETURN
      50       110994 :       IF (.NOT. ALL(ieee_is_finite([mu, kt, gamma]))) RETURN
      51        27748 :       IF (kt <= 0.0_dp .OR. gamma <= 0.0_dp) RETURN
      52        85628 :       IF (.NOT. ALL(ieee_is_finite(energies))) RETURN
      53      1883544 :       IF (.NOT. ALL(ieee_is_finite(REAL(current, dp))) .OR. &
      54              :           .NOT. ALL(ieee_is_finite(AIMAG(current)))) RETURN
      55       110970 :       DO a = 1, 3
      56       913944 :          scale = MAX(1.0_dp, MAXVAL(ABS(current(:, :, a))))
      57       941686 :          IF (MAXVAL(ABS(current(:, :, a) + CONJG(TRANSPOSE(current(:, :, a))))) > antihermiticity_tolerance*scale) RETURN
      58              :       END DO
      59        27742 :       total = 0.0_dp
      60        85618 :       DO n = 1, nstate
      61       166188 :          DO m = n + 1, nstate
      62        80570 :             delta = energies(n) - energies(m)
      63        80570 :             IF (.NOT. ieee_is_finite(delta)) RETURN
      64        80570 :             IF (delta == 0.0_dp) CYCLE
      65        80492 :             scale = MAX(gamma, ABS(delta))
      66              :             factor = 2.0_dp*fermi_slope(energies(n), energies(m), mu, kt)* &
      67        80492 :                      (delta/scale)/scale/((gamma/scale)**2 + (delta/scale)**2)
      68        80492 :             IF (.NOT. ieee_is_finite(factor)) RETURN
      69       299352 :             DO b = 2, 3
      70       483030 :                DO a = 1, b - 1
      71       402460 :                   total(a, b) = total(a, b) + factor*AIMAG(current(m, n, a)*current(n, m, b))
      72              :                END DO
      73              :             END DO
      74              :          END DO
      75              :       END DO
      76       360646 :       IF (.NOT. ALL(ieee_is_finite(total))) RETURN
      77       360646 :       hall = total - TRANSPOSE(total)
      78        27742 :       status = 0
      79              :    END SUBROUTINE hall_response
      80              : END MODULE kubo_hall