init_tmat_rs_hub_spin_transcorr Subroutine

    subroutine init_tmat_rs_hub_spin_transcorr()
#ifdef DEBUG_
        character(*), parameter :: this_routine = "init_tmat_rs_hub_spin_transcorr"
#endif
        integer :: n_sites, i, j
        real(dp) :: elem
        integer :: iunit

        ASSERT(associated(lat))

        if (allocated(tmat_rs_hub_spin_transcorr) .and. .not. t_recalc_tmat) then
            return
        else
            if (allocated(tmat_rs_hub_spin_transcorr)) deallocate(tmat_rs_hub_spin_transcorr)
        end if

        call init_hop_trancorr_fac_cached_vec(trans_corr_param, lat)

        n_sites = lat%get_nsites()
        if (t_print_tmat) then
            iunit = get_free_unit()
            open(iunit, file="TMAT")
        end if
        root_print "initializing spin-dependent TMAT"

        ! tmat is stored with spin-orbitals!
        allocate(tmat_rs_hub_spin_transcorr(2 * n_sites, 2 * n_sites))
        tmat_rs_hub_spin_transcorr = 0.0_dp

        do i = 1, n_sites
            do j = 1, n_sites
                ! the alpha spin are the transcorrelated ones! even numbers!
                ! but the sum_ function gets accessed with spatial orbitals!
                elem = nOccBeta * uhub * sum_spin_transcorr_factor(i, j)
                if (abs(elem) > matele_cutoff) then
                    if (t_print_tmat) then
                        write(iunit, *) 2 * i, 2 * j, elem
                    end if
                    tmat_rs_hub_spin_transcorr(2 * i, 2 * j) = elem
                end if
            end do
        end do

        if (t_print_tmat) then
            close(iunit)
        end if
        root_print "Done!"

    end subroutine init_tmat_rs_hub_spin_transcorr