subroutine ab_initio_find_tau_from_refdet_conn()
type(excit_gen_store_type) :: store, store2
logical :: tAllExcitFound, tParity, tSameFunc, tSwapped, tSign
character(len=*), parameter :: this_routine = "find_tau_from_refdet_conn"
integer :: ex(2, maxExcit), ex2(2, maxExcit), exflag, iMaxExcit, nStore(6), nExcitMemLen(1)
integer, allocatable :: Excitgen(:)
real(dp) :: nAddFac, MagHel, pGen, pGenFac
HElement_t(dp) :: hel
integer :: ic, nJ(nel), nJ2(nel), iExcit, ex_saved(2, maxExcit)
integer(kind=n_int) :: iLutnJ(0:niftot), iLutnJ2(0:niftot)
real(dp) :: new_tau
type(ExcitGenSessionType) :: session
ASSERT(.not. tGUGA)
ASSERT(.not. t_k_space_hubbard)
new_tau = huge(new_tau)
nAddFac = MaxWalkerBloom
tAllExcitFound = .false.
Ex_saved(:, :) = 0
exflag = 3
tSameFunc = .false.
call init_excit_gen_store(store)
call init_excit_gen_store(store2)
store%tFilled = .false.
store2%tFilled = .false.
CALL construct_class_counts(ProjEDet(:, 1), store%ClassCountOcc, &
store%ClassCountUnocc)
store%tFilled = .true.
if (tKPntSym) then
!TODO: It REALLY needs to be fixed so that we don't need to do this!!
!Setting up excitation generators that will work with kpoint sampling
iMaxExcit = 0
nStore(:) = 0
CALL gensymexcitit2par_worker(ProjEDet(:, 1), NEl, G1, nBasis, .TRUE., nExcitMemLen, nJ, iMaxExcit, nStore, exFlag, 1, nEl)
allocate(EXCITGEN(nExcitMemLen(1)))
EXCITGEN(:) = 0
CALL gensymexcitit2par_worker(ProjEDet(:, 1), NEl, G1, nBasis, .TRUE., EXCITGEN, nJ, iMaxExcit, nStore, exFlag, 1, nEl)
end if
do while (.not. tAllExcitFound)
if (tKPntSym) then
call gensymexcitit2par_worker(ProjEDet(:, 1), nel, G1, nBasis, .false., EXCITGEN, nJ, iExcit, nStore, exFlag, 1, nEl)
if (nJ(1) == 0) exit
!Calculate ic, tParity and Ex
call EncodeBitDet(nJ, iLutnJ)
Ex(:, :) = 0
ic = FindBitExcitlevel(iLutnJ, iLutRef(:, 1), 2)
ex(1, 1) = ic
call GetExcitation(ProjEDet(:, 1), nJ, Nel, ex, tParity)
else
if (tReltvy) then
call GenExcitations4(session, ProjEDet(:, 1), nJ, exflag, ex_saved, tParity, tAllExcitFound, .false.)
else
CALL GenExcitations3(ProjEDet(:, 1), iLutRef(:, 1), nJ, exflag, Ex_saved, tParity, tAllExcitFound, .false.)
end if
IF (tAllExcitFound) EXIT
Ex(:, :) = Ex_saved(:, :)
if (Ex(2, 2) == 0) then
ic = 1
else
ic = 2
end if
call EncodeBitDet(nJ, iLutnJ)
end if
! Exclude an excitation if it isn't symmetry allowed.
! Note that GenExcitations3 is not perfect, especially if there
! additional restrictions, such as LzSymmetry.
if (.not. SymAllowedExcit(ProjEDet(:, 1), nJ, ic, ex)) &
cycle
if (tHPHF) then
if (.not. TestClosedShellDet(iLutnJ)) then
CALL ReturnAlphaOpenDet(nJ, nJ2, iLutnJ, iLutnJ2, .true., .true., tSwapped)
if (tSwapped) then
!Have to recalculate the excitation matrix.
ic = FindBitExcitLevel(iLutnJ, iLutRef(:, 1), 2)
ex(:, :) = 0
if (ic <= max_ex_level) then
ex(1, 1) = ic
call GetBitExcitation(iLutRef(:, 1), iLutnJ, Ex, tParity)
end if
end if
end if
hel = hphf_off_diag_helement_norm(ProjEDet(:, 1), nJ, iLutRef(:, 1), iLutnJ)
else
hel = get_helement(ProjEDet(:, 1), nJ, ic, ex, tParity)
end if
MagHel = abs(hel)
!Find pGen (nI -> nJ)
if (tHPHF) then
call CalcPGenHPHF(ProjEDet(:, 1), iLutRef(:, 1), nJ, iLutnJ, ex, store%ClassCountOcc, &
store%ClassCountUnocc, pDoubles, pGen, tSameFunc)
else
call CalcNonUnipGen(ProjEDet(:, 1), ilutRef(:, 1), ex, ic, store%ClassCountOcc, store%ClassCountUnocc, pDoubles, pGen)
end if
if (tSameFunc) cycle
if (MagHel > 0.0_dp) then
pGenFac = pGen * nAddFac / MagHel
new_tau = clamp(&
merge(new_tau, min(pGenFac, new_tau), near_zero(pGenFac)), &
min_tau, max_tau)
end if
!Find pGen(nJ -> nI)
CALL construct_class_counts(nJ, store2%ClassCountOcc, &
store2%ClassCountUnocc)
store2%tFilled = .true.
if (tHPHF) then
ic = FindBitExcitLevel(iLutnJ, iLutRef(:, 1), 2)
ex2(:, :) = 0
if (ic <= max_ex_level) then
ex2(1, 1) = ic
call GetBitExcitation(iLutnJ, iLutRef(:, 1), Ex2, tSign)
end if
call CalcPGenHPHF(nJ, iLutnJ, ProjEDet(:, 1), iLutRef(:, 1), ex2, store2%ClassCountOcc, &
store2%ClassCountUnocc, pDoubles, pGen, tSameFunc)
else
ex2(1, :) = ex(2, :)
ex2(2, :) = ex(1, :)
call CalcNonUnipGen(nJ, ilutnJ, ex2, ic, store2%ClassCountOcc, store2%ClassCountUnocc, pDoubles, pGen)
end if
if (tSameFunc) cycle
if (MagHel > 0.0_dp) then
pGenFac = pGen * nAddFac / MagHel
new_tau = clamp(&
merge(new_tau, min(pGenFac, new_tau), near_zero(pGenFac)), &
min_tau, max_tau)
end if
end do
call clean_excit_gen_store(store)
call clean_excit_gen_store(store2)
if (tKPntSym) deallocate(EXCITGEN)
if (new_tau > 0.075_dp) then
new_tau = clamp(0.075_dp, min_tau, max_tau)
write(stdout, "(A,F8.5,A)") "Small system. Setting initial timestep to be ", Tau, " although this &
&may be inappropriate. Care needed"
end if
call assign_value_to_tau(new_tau, this_routine)
end subroutine ab_initio_find_tau_from_refdet_conn