subroutine csf_to_sds_ilut(csf, sds, weights, ms, coeff)
integer(n_int), intent(in) :: csf(0:GugaBits%len_tot)
real(dp), intent(in), optional :: ms
real(dp), intent(in), optional :: coeff
integer(n_int), intent(out), allocatable :: sds(:,:)
real(dp), intent(out), allocatable :: weights(:)
character(*), parameter :: this_routine = "csf_to_sds_ilut"
real(dp) :: ms_
integer :: spin, n_alpha, n_beta
integer :: i, j, step(nSpatorbs), delta_k
integer :: nI(nel)
real(dp) :: bVec(nSpatorbs), aVec(nSpatorbs), lambda_k, x, coeff_
integer(n_int), allocatable :: all_sds(:,:)
real(dp), allocatable :: all_weights(:)
def_default(coeff_, coeff, 1.0_dp)
! only works for heisenberg model for now..
if (any(calcOcc_vector_int(csf) /= 1)) then
call stop_all(this_routine, "only implemented for heisenberg for now")
end if
if (near_zero(coeff_)) then
allocate(sds(0:GugaBits%len_tot, 0))
allocate(weights(0))
return
end if
spin = abs(return_ms(csf))
def_default(ms_, ms, spin/2.)
! construct SDs by attaching (N/2 + ms) up spins and (N/2 - ms) down spins
n_alpha = nint(nSpatorbs / 2. + ms_)
n_beta = nint(nSpatorbs / 2. - ms_)
all_sds = create_all_open_shell_dets(nSpatorbs, n_beta, n_alpha)
allocate(all_weights(size(all_sds,2)), source = 0.0_dp)
step = calcStepvector(csf(0:GugaBits%len_orb))
bVec = calcB_vector_ilut(csf(0:GugaBits%len_orb))
aVec = real(([(i, i = 1,nSpatorbs)] - bVec),dp) / 2.0
do i = 1, size(all_sds,2)
x = 1.0_dp
call decode_bit_det(nI, all_sds(:,i))
do j = 1, nSpatorbs
if (.not. near_zero(x)) then
delta_k = mod(get_spin(nI(j)),2)
lambda_k = get_preceeding_opposites(nI, j)
select case (step(j))
case (1)
x = x * sqrt((aVec(j) + bVec(j) - lambda_k)/bVec(j))
case (2)
x = x * (-1.0_dp)**(bVec(j)+delta_k) * &
sqrt((lambda_k - aVec(j) + 1.0_dp)/(bVec(j)+2.0_dp))
case (3)
x = x * (-1.0_dp)**bVec(j)
end select
end if
end do
all_weights(i) = x
end do
j = 1
do i = 1, size(all_sds,2)
if (.not. near_zero(all_weights(i))) then
all_sds(:,j) = all_sds(:,i)
all_weights(j) = all_weights(i)
j = j + 1
end if
end do
! i might have to normalize
all_weights(1:j-1) = coeff_ * all_weights(1:j-1) / sqrt(sum(all_weights(1:j-1)**2))
allocate(sds(0:GugaBits%len_tot, j - 1), source = 0_n_int)
allocate(weights(j-1), source = all_weights(1:j-1))
do i = 1, j - 1
sds(0:0, i) = all_sds(:, i)
call encode_matrix_element(sds(:,i), all_weights(i), 1)
end do
end subroutine csf_to_sds_ilut