subroutine tetrahedron_P_matrix_init(P_matrix)
use w90_constants, only: dp
implicit none
real(kind=dp), dimension(4, 20), intent(inout) :: P_matrix
! correction, Kawamura PRB 89 094515
P_matrix(1, 1:4) = real((/1440, 0, 30, 0/), dp)
P_matrix(2, 1:4) = real((/0, 1440, 0, 30/), dp)
P_matrix(3, 1:4) = real((/30, 0, 1440, 0/), dp)
P_matrix(4, 1:4) = real((/0, 30, 0, 1440/), dp)
!
P_matrix(1, 5:8) = real((/-38, 7, 17, -28/), dp)
P_matrix(2, 5:8) = real((/-28, -38, 7, 17/), dp)
P_matrix(3, 5:8) = real((/17, -28, -38, 7/), dp)
P_matrix(4, 5:8) = real((/7, 17, -28, -38/), dp)
!
P_matrix(1, 9:12) = real((/-56, 9, -46, 9/), dp)
P_matrix(2, 9:12) = real((/9, -56, 9, -46/), dp)
P_matrix(3, 9:12) = real((/-46, 9, -56, 9/), dp)
P_matrix(4, 9:12) = real((/9, -46, 9, -56/), dp)
!
P_matrix(1, 13:16) = real((/-38, -28, 17, 7/), dp)
P_matrix(2, 13:16) = real((/7, -38, -28, 17/), dp)
P_matrix(3, 13:16) = real((/17, 7, -38, -28/), dp)
P_matrix(4, 13:16) = real((/-28, 17, 7, -38/), dp)
!
P_matrix(1, 17:20) = real((/-18, -18, 12, -18/), dp)
P_matrix(2, 17:20) = real((/-18, -18, -18, 12/), dp)
P_matrix(3, 17:20) = real((/12, -18, -18, -18/), dp)
P_matrix(4, 17:20) = real((/-18, 12, -18, -18/), dp)
!
P_matrix(1:4, 1:20) = P_matrix(1:4, 1:20)/1260.0_dp
end subroutine tetrahedron_P_matrix_init
