function IsTruncated(nI) result(truncated)
use double_occ_mod, only: get_doublons_holons, count_double_orbs
use lattice_mod, only: lat
use Determinants, only: DefDet
use CalcData, only: tTruncateHubbard, pairing_criterion, repel_criterion, &
af_criterion, af_criterion_bracket, max_doublons, max_holons
integer, intent(in) :: nI(nel)
integer(n_int) :: ilut(0:niftot), ilutRef(0:niftot)
logical :: truncated
integer :: doublons(1:nbasis/2), holons(1:nbasis/2), doublons_count, holons_count, neel_state(nel), exlevel, &
unpaired_count
integer, allocatable :: neighbors(:)
integer :: i, j, af_param, max_diff
logical :: paired
procedure(check_close), pointer :: are_close
if (lat%get_name() == 'square' .or. lat%get_name() == 'rectangle' .or. &
lat%get_name() == 'hexagonal' .or. lat%get_name() == 'honeycomb') then
are_close => are_close_direct
else
are_close => are_close_via_neighbors
end if
truncated = .false.
call EncodeBitDet (nI, iLut)
if(tTruncateHubbard)then
call get_doublons_holons(ilut(:NIfD), doublons, holons, doublons_count, holons_count)
if((max_doublons>=0 .and. doublons_count>max_doublons) .or. (max_holons>=0 .and. holons_count>max_holons))then
truncated = .true.
return
endif
!-----------------------------------------------------------------!
if(pairing_criterion>=0)then
unpaired_count = 0
do i=1, doublons_count
paired = .false.
if(holons_count == 0) paired = .true.
do j=1, holons_count
if(are_close(doublons(i), holons(j), pairing_criterion))then
paired = .true.
exit
endif
enddo
if(.not. paired) then
unpaired_count = unpaired_count + 1
if(tTruncHubbAllClose)then
max_diff = 0
else
max_diff = max(doublons_count-holons_count,0)
endif
if(unpaired_count > max_diff)then
!More unpaired than excess doublons, truncate
truncated = .true.
return
endif
endif
enddo
unpaired_count = 0
do j=1, holons_count
paired = .false.
if(doublons_count == 0) paired = .true.
do i=1, doublons_count
if(are_close(doublons(i), holons(j), pairing_criterion))then
paired = .true.
exit
endif
enddo
if(.not. paired) then
if(tTruncHubbAllClose)then
max_diff = 0
else
max_diff = max(holons_count-doublons_count,0)
endif
unpaired_count = unpaired_count + 1
if(unpaired_count > max_diff)then
!More unpaired than excess holons, truncate
truncated = .true.
return
endif
endif
enddo
endif
!-----------------------------------------------------------------!
if(repel_criterion>=0)then
do i=1, doublons_count
!Check if there is another doublon nearby
do j=i+1, doublons_count
if(are_close(doublons(i), doublons(j), repel_criterion))then
truncated = .true.
return
endif
enddo
enddo
do i=1, holons_count
!Check if there is another holon nearby
do j=i+1, holons_count
if(are_close(holons(i), holons(j), repel_criterion))then
truncated = .true.
return
endif
enddo
enddo
endif
!-----------------------------------------------------------------!
if(af_criterion>=0)then
if (tExLevelAFCriterion) then
call EncodeBitDet(DefDet, ilutRef)
exlevel = FindBitExcitLevel(ilutRef, ilut)
if(min(exlevel, nel-exlevel)-doublons_count>af_criterion)then
truncated = .true.
endif
else
af_param = af_func(ilut)
if (tAFTruncBracket .and. af_criterion_bracket>=0) then
if (af_param > af_criterion .and. &
af_param < af_criterion_bracket) &
truncated = .true.
else
if (af_param > af_criterion) truncated = .true.
end if
end if
end if
endif
contains
function af_func(ilut) result(val)
integer(n_int), intent(in) :: ilut(0:niftot)
integer :: i, j, idx, val
integer, allocatable :: neighbors(:)
val = lat%get_total_links()
do i = 1, lat%get_nsites()
neighbors = lat%get_neighbors(i)
do idx = 1, size(neighbors)
j = neighbors(idx)
if (j > i) cycle
if (IsOcc(ilut, 2*i-1) .and. IsOcc(ilut, 2*j)) val = val - 1
if (IsOcc(ilut, 2*i-1) .and. IsOcc(ilut, 2*j-1)) val = val + 1
if (IsOcc(ilut, 2*i) .and. IsOcc(ilut, 2*j)) val = val + 1
if (IsOcc(ilut, 2*i) .and. IsOcc(ilut, 2*j-1)) val = val - 1
end do
end do
end function af_func
function are_close_direct(site1, site2, criterion) result(nearby)
integer, intent(in) :: site1, site2
integer, intent(in) :: criterion
integer :: r1(3), r2(3), x1, x2, y1, y2, dx, dy, lx, ly
logical :: nearby, t_hc
character(*), parameter :: this_routine = 'are_close_direct'
t_hc = ((lat%get_name() == 'hexagonal') .or. (lat%get_name() == 'honeycomb'))
if(t_hc)then
lx = 4*length_x
ly = 4*length_y
else
lx = length_x
ly = length_y
endif
r1 = lat%get_r_vec(site1)
r2 = lat%get_r_vec(site2)
x1 = r1(1)
x2 = r2(1)
y1 = r1(2)
y2 = r2(2)
if(t_open_bc_x)then
dx = abs(x1-x2)
else
dx = min(modulo(x1-x2, lx), modulo(x2-x1, lx))
endif
if(t_open_bc_y)then
dy = abs(y1-y2)
else
dy = min(modulo(y1-y2, ly), modulo(y2-y1, ly))
endif
if(t_hc)then
select case(criterion)
case (0)
nearby = ((dx<=1) .and. (dy<=1))
case (1)
call stop_all(this_routine, "PC1 not working for honeycomb lattice!")
case (2)
nearby = (((dx+dy)<=2) .or. ((dx<=1) .and. (dy<=2)))
case default
nearby = .false.
endselect
else
select case(criterion)
case (0)
nearby = (dx+dy)<=1
case (1)
nearby = (dx<=1 .and. dy<=1)
case (2)
nearby = (dx+dy)<=2
case default
nearby = .false.
endselect
endif
end function are_close_direct
pure function are_close_via_neighbors(site1, site2, criterion) result(nearby)
integer, intent(in) :: site1, site2
integer, intent(in) :: criterion
integer, allocatable :: diag_neighs(:)
logical :: nearby
if (criterion == 0) then
nearby = any(lat%get_neighbors(site1) == site2)
else if (criterion == 1) then
diag_neighs = find_diag_neighs(&
lat%get_nn_neighbors(site1), &
lat%get_num_neighbors(site1)**2)
nearby = any(lat%get_neighbors(site1) == site2) &
.or. any(diag_neighs == site2)
else if (criterion == 2) then
nearby = any(lat%get_neighbors(site1) == site2) &
.or. any(lat%get_nn_neighbors(site1) == site2)
end if
end function are_close_via_neighbors
pure function find_diag_neighs(nn_neighs, num_nn_neighs) result(diag_neighs)
integer, intent(in) :: num_nn_neighs, nn_neighs(num_nn_neighs)
integer :: res(num_nn_neighs)
integer, allocatable :: diag_neighs(:)
integer :: i, j, ind, k
k = 0
do i = 1, size(nn_neighs)
ind = nn_neighs(i)
do j = i + 1, size(nn_neighs)
if (nn_neighs(j) == ind) then
k = k + 1
res(k) = ind
end if
end do
end do
diag_neighs = res(1:k)
end function
end function IsTruncated