apply_perturbation Subroutine

    subroutine apply_perturbation(perturb, ndets, dets_in, dets_out, phase)

        ! Take in a list of determinants (dets_in) and apply a pertubation
        ! to each determinant. As we go we shuffle down determinants to fill in
        ! gaps opened up by removed determinants.

        ! WARNING: This routine uses SpanwedParts and the annihilation routine
        ! SendProcNewParts to perform communication of the perturbed
        ! determinants. It will overwrite whatever is in SpawnedParts.

        use AnnihilationMod, only: SendProcNewParts
        use bit_rep_data, only: NIfTot, nifd, extract_sign
        use bit_reps, only: encode_sign, decode_bit_det
        use DetBitOps, only: ilut_lt, ilut_gt
        use load_balance_calcnodes, only: DetermineDetNode
        use FciMCData, only: SpawnedParts, SpawnedParts2
        use FciMCData, only: ValidSpawnedList, InitialSpawnedSlots
        use sort_mod, only: sort
        use SystemData, only: nel

        type(perturbation), intent(in) :: perturb
        integer, intent(inout) :: ndets
        integer(n_int), intent(in) :: dets_in(0:, :) ! First dimension must be 0:NIfTot
        integer(n_int), intent(out) :: dets_out(0:, :) ! First dimension must be 0:NIfTot
        real(dp), intent(in), optional :: phase ! add some phase to the perturbation

        integer(n_int) :: ilut(0:NIfTot)
        integer(n_int), pointer :: PointTemp(:, :)
        integer :: i, nremoved, proc, run
        integer :: nI(nel)
        real(dp) :: tmp_sign(lenof_sign), tmp_real
        character(*), parameter :: this_routine = 'apply_perturbation'

        ! If the perturbation is the identity operator then just return.
        ! rneci_consitency: Possible optimization: Define behaviour in this case as copying
        if (perturb%nannihilate == 0 .and. perturb%ncreate == 0) then
            dets_out = dets_in
            return
        end if

        nremoved = 0
        ! Reset the spawning slot positions in SpawnedParts.
        ValidSpawnedList = InitialSpawnedSlots
        do i = 1, ndets
            ! Copy the ilut so that we don't alter the input list.
            ilut = dets_in(:, i)
            call perturb_det(ilut, perturb)

            if (all(ilut(0:nifd) == 0_n_int)) then
                nremoved = nremoved + 1
            else
                call decode_bit_det(nI, ilut)
                proc = DetermineDetNode(nel, nI, 0)
                ! If a phase factor is to be added, do it now
                if (present(phase)) then
                    call extract_sign(ilut, tmp_sign)
                    do run = 1, inum_runs
                        ! multiply by exp(i*phase)
                        tmp_real = tmp_sign(min_part_type(run))
                        tmp_sign(min_part_type(run)) = cos(phase) * tmp_sign(min_part_type(run)) &
                                                       - sin(phase) * tmp_sign(max_part_type(run))
                        tmp_sign(max_part_type(run)) = sin(phase) * tmp_real + cos(phase) * &
                                                       max_part_type(run)
                    end do
                    call encode_sign(ilut, tmp_sign)
                end if
                SpawnedParts(0:NIfTot, ValidSpawnedList(proc)) = ilut
                ValidSpawnedList(proc) = ValidSpawnedList(proc) + 1
                if (checkValidSpawnedList(proc)) then
                    call stop_all(this_routine, "Not enough memory to apply perturbation")
                end if
            end if
        end do

        if (allocated(perturb%ann_orbs) .and. allocated(perturb%crtn_orbs)) &
            write(stdout, *) "Transfering from orbital ", perturb%ann_orbs(1), &
            " to ", perturb%crtn_orbs(1)
        ndets = ndets - nremoved

        write(stdout, *) "Communicating perturbed dets"
        ! Send perturbed determinants to their new processors.
        call SendProcNewParts(ndets, tSingleProc=.false.)

        ! The result of SendProcNewParts is now stored in an array pointed to
        ! by SpawnedParts2. We want it in an array pointed to by SpawnedParts,
        ! so swap the pointers around.
        ! Why not just use SpawnedParts2 below?
        PointTemp => SpawnedParts2
        SpawnedParts2 => SpawnedParts
        SpawnedParts => PointTemp
        nullify (PointTemp)

        ! Now move the contents of SpawnedParts to ilut_list.
        do i = 1, ndets
            dets_out(:, i) = SpawnedParts(0:NIfTot, i)
        end do
        call sort(dets_out(:, 1:ndets), ilut_lt, ilut_gt)

    end subroutine apply_perturbation