communicate_kp_matrices Subroutine

    subroutine communicate_kp_matrices(overlap_matrix, hamil_matrix, spin_matrix)

        ! Add all the overlap and projected Hamiltonian matrices together, with
        ! the result being held only on the root node.

        use MPI_wrapper, only: root

        real(dp), intent(inout) :: overlap_matrix(:, :)
        real(dp), intent(inout) :: hamil_matrix(:, :)
        real(dp), optional, intent(inout) :: spin_matrix(:, :)

        real(dp), allocatable :: mpi_mat_in(:, :), mpi_mat_out(:, :)
        integer :: nrow, ncol

        nrow = size(hamil_matrix, 1)

        if (present(spin_matrix)) then
            ncol = 3 * nrow
        else
            ncol = 2 * nrow
        end if

        allocate(mpi_mat_in(nrow, ncol))
        allocate(mpi_mat_out(nrow, ncol))

        mpi_mat_in(1:nrow, 1:nrow) = overlap_matrix
        mpi_mat_in(1:nrow, nrow + 1:2 * nrow) = hamil_matrix
        if (present(spin_matrix)) mpi_mat_in(1:nrow, 2 * nrow + 1:3 * nrow) = spin_matrix

        call MPISum(mpi_mat_in, mpi_mat_out)

        if (iProcIndex == root) then
            overlap_matrix = mpi_mat_out(1:nrow, 1:nrow)
            hamil_matrix = mpi_mat_out(1:nrow, nrow + 1:2 * nrow)
            if (present(spin_matrix)) spin_matrix = mpi_mat_out(1:nrow, 2 * nrow + 1:3 * nrow)
        end if

        deallocate(mpi_mat_in)
        deallocate(mpi_mat_out)

    end subroutine communicate_kp_matrices