subroutine init_umat_rs_hub_transcorr()
! for now do it really stupidly without any concern for the symmetry
! of the integrals
integer :: n_sites, i, j, k, l
#ifdef DEBUG_
character(*), parameter :: this_routine = "init_umat_rs_hub_transcorr"
#endif
real(dp) :: elem
integer :: iunit
ASSERT(associated(lat))
if (allocated(umat_rs_hub_trancorr_hop) .and. .not. t_recalc_umat) then
! already initialized
return
else
if (allocated(umat_rs_hub_trancorr_hop)) deallocate(umat_rs_hub_trancorr_hop)
end if
! try to fetch the stored matrix elements also for each orbital after.
call init_hop_trancorr_fac_cached_vec(trans_corr_param, lat)
n_sites = lat%get_nsites()
! create an fcidump file:
if (t_print_umat) then
iunit = get_free_unit()
open(iunit, file='UMAT')
end if
root_print "initializing UMAT:"
! with the correct header
! and also try to allocate a umat_cache
allocate(umat_rs_hub_trancorr_hop(n_sites, n_sites, n_sites, n_sites))
umat_rs_hub_trancorr_hop = 0.0_dp
do i = 1, n_sites
do j = 1, n_sites
do k = 1, n_sites
do l = 1, n_sites
elem = uhub * sum_hop_transcorr_factor(i, j, k, l)
! write to the dumpfile
if (abs(elem) > matele_cutoff) then
if (t_print_umat) then
write(iunit, *) i, j, k, l, elem
end if
! and also store in the umat
umat_rs_hub_trancorr_hop(i, j, k, l) = elem
end if
end do
end do
end do
end do
if (t_print_umat) then
close(iunit)
root_print "Done"
end if
end subroutine init_umat_rs_hub_transcorr