subroutine perturb_det(ilut, perturb)
! This routine takes a determinant encoded in ilut and applies a
! collection of creation and annihilation operators, which together is
! referred to as the perturbation operator.
! ***IMPORTANT*** All annihilation operators are applied before (ie to
! the right of) all creation operators.
! The a_bits(i)'th bit of ilut(a_elems(i)) should be the bit of the
! orbital for the i'th annihilation operator, and similarly for
! c_bits and c_elems for creation operators.
use bit_rep_data, only: NIfTot, nifd, NIfD, extract_sign
use bit_reps, only: encode_sign
use DetBitOps, only: CountBits
integer(n_int), intent(inout) :: ilut(0:NIfTot)
type(perturbation), intent(in) :: perturb
integer(n_int) :: combined_mask(0:NIfD), new_mask(0:NIfD), ones
integer(n_int) :: original_ilut(0:NIfTot), masked_ilut(0:NIfD)
integer :: i, j, num_minus_signs, nswap
real(dp) :: new_sign_factor, real_sign(lenof_sign)
! If the perturbation is the identity operator then just return.
if (perturb%nannihilate == 0 .and. perturb%ncreate == 0) return
original_ilut = ilut
associate (a_elems => perturb%ann_elems, a_bits => perturb%ann_bits, &
c_elems => perturb%crtn_elems, c_bits => perturb%crtn_bits)
! First, loop over all creation and annihilation operators and see if they
! destroy the determinant encoded in ilut.
do i = 1, perturb%nannihilate
if (.not. btest(ilut(a_elems(i)), a_bits(i))) then
ilut(0:nifd) = 0_n_int
return
end if
ilut(a_elems(i)) = ibclr(ilut(a_elems(i)), a_bits(i))
end do
do i = 1, perturb%ncreate
if (btest(ilut(c_elems(i)), c_bits(i))) then
ilut(0:nifd) = 0_n_int
return
end if
ilut(c_elems(i)) = ibset(ilut(c_elems(i)), c_bits(i))
end do
! If we get to this point then the determinant was not destroyed by the
! creation and annihilation operators, else we would have returned.
! Now, determine whether or not the application of these operators
! introduces a minus sign.
! This is done in two steps. First, we find the number of
! minus signs introduced assuming all of the operators in the
! perturbation operator are already ordered by their NECI orbital
! number. Second, we find the number of permutations necessary to
! reorder the annihilation and creation operators into this order.
! We do the first part here:
! We need to consider moving the creation and annihilation operators
! in the perturbing operator through the creation operators which
! create the original determinant (by acting on the vacuum).
! For a creation or annihilation operator for orbital p in the
! perturbing operator, we want to consider moving it through creation
! operators in the original determinant for orbitals less than p.
! Suppose there are N such operators. Then the ordering introduces a
! factor of (-1)^N. So we just need to count the number of orbitals
! before orbital p in the determinant. We can do this by creating a bit
! mask with 1's for all orbitals *before* (and not including) p in the
! original determinant, and then performing an and operation between
! this mask and the ilut, and then counting the number of bits set in
! the resulting bitstring.
! However, in general we want to consider several annihilation and
! creation operators in the perturbing operator. Instead of doing the
! above for each operator, which is slow, we can do it all in one go.
! If a creation operator in the original ilut is 'passed' an odd number
! of times then we want to count it, otherwise we don't. So we want to
! create a bit mask that has 1's for all orbitals that would be passed
! an odd number of times by the applied creation and annihilation
! operators in the perturbation operator, and 0's for all other
! orbitals. This can be done by performing an ieor on the bit masks
! for induvidual operators. Once an exclusive or has been performed for
! each annihilation and creation operator in the perturbing operator,
! we can finally perform an and with the original ilut,and then count
! the number of set bits. If this is odd, include a minus sign.
ones = not(0_n_int)
combined_mask = 0_n_int
do i = 1, perturb%nannihilate
do j = 0, NIfD
! Create a mask which has 0's in all orbitals which might be
! passed in the ordering, and 1's elsewhere.
if (j < a_elems(i)) then
new_mask(j) = 0_n_int
else if (j == a_elems(i)) then
new_mask(j) = ishft(ones, a_bits(i))
else
new_mask(j) = ones
end if
end do
! Now make it so that it has 1's for all orbitals that might be
! passed and 0's elsewhere.
new_mask = not(new_mask)
! Now perform an ieor to combine this mask with the combined masks
! for all previously considered orbitals.
combined_mask = ieor(combined_mask, new_mask)
end do
! Now do exactly the same for the creation operators.
do i = 1, perturb%ncreate
do j = 0, NIfD
if (j < c_elems(i)) then
new_mask(j) = 0_n_int
else if (j == c_elems(i)) then
new_mask(j) = ishft(ones, c_bits(i))
else
new_mask(j) = ones
end if
end do
new_mask = not(new_mask)
combined_mask = ieor(combined_mask, new_mask)
end do
! Finally apply the mask and count the resulting number of set bits.
masked_ilut = iand(combined_mask, original_ilut(0:NIfD))
num_minus_signs = CountBits(masked_ilut, NIfD)
end associate
! The above finds the number of permutations assuming that the orbitals
! were all ordered according to NECI's orbital ordering convention
! (beta, alpha, beta, alpha...). So, to get the final permutation, we
! need to see if there is an extra minus sign from sorting the
! requested ordering into NECI's ordering.
! For now this is done in the completley naive way which is order N^2
! for N orbitals. Hopefully this code isn't used any performance
! critical, but the speed of this can be improved later by using a
! merge sort (look up algorithms for counting inversions!)
! In this naive approach, we simply loop over all orbitals and count
! the number of orbitals to the right which have a smaller index. This
! is the total number of permutations (inversions) which must be made.
! This is done for both annihilation and creation operators.
associate(a_orbs => perturb%ann_orbs, c_orbs => perturb%crtn_orbs)
nswap = 0
do i = 1, perturb%nannihilate - 1
do j = i + 1, perturb%nannihilate
if (a_orbs(i) > a_orbs(j)) nswap = nswap + 1
end do
end do
do i = 1, perturb%ncreate - 1
do j = i + 1, perturb%ncreate
if (c_orbs(i) > c_orbs(j)) nswap = nswap + 1
end do
end do
end associate
num_minus_signs = num_minus_signs + nswap
! Finally, calculate and apply the new sign.
new_sign_factor = (-1.0_dp)**num_minus_signs
call extract_sign(ilut, real_sign)
real_sign = real_sign * new_sign_factor
call encode_sign(ilut, real_sign)
end subroutine perturb_det