perturb_det Subroutine

private subroutine perturb_det(ilut, perturb)


Type IntentOptional Attributes Name
integer(kind=n_int), intent(inout) :: ilut(0:NIfTot)
type(perturbation), intent(in) :: perturb


Source Code

Source Code

    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
                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
                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))
                        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))
                        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