CalcHashTableStats Subroutine

public subroutine CalcHashTableStats(TotWalkersNew, iter_data)


Type IntentOptional Attributes Name
integer, intent(inout) :: TotWalkersNew
type(fcimc_iter_data), intent(inout) :: iter_data


Source Code

Source Code

    subroutine CalcHashTableStats(TotWalkersNew, iter_data)

        use FciMCData, only: CurrentDets, n_prone_dets
        use LoggingData, only: FCIMCDebug
        use bit_rep_data, only: IlutBits

        integer, intent(inout) :: TotWalkersNew
        type(fcimc_iter_data), intent(inout) :: iter_data

        integer :: i, j, AnnihilatedDet, lbnd, ubnd, part_type
        real(dp) :: CurrentSign(lenof_sign)
        real(dp) :: pRemove, r
        integer :: nI(nel), run
        logical :: tIsStateDeterm
        real(dp) :: scaledOccupiedThresh
        character(*), parameter :: t_r = 'CalcHashTableStats'

        if (.not. bNodeRoot) return

        TotParts = 0.0_dp
        norm_psi_squared = 0.0_dp
        norm_semistoch_squared = 0.0_dp
        iHighestPop = 0
        AnnihilatedDet = 0
        tIsStateDeterm = .false.
        InstNoAtHf = 0.0_dp
        n_prone_dets = 0

        if (TotWalkersNew > 0) then
            do i = 1, TotWalkersNew

                call extract_sign(CurrentDets(:, i), CurrentSign)
                if (tSemiStochastic) tIsStateDeterm = test_flag_multi(CurrentDets(:, i), flag_deterministic)

                if (IsUnoccDet(CurrentSign) .and. (.not. tIsStateDeterm)) then
                    if (.not. tAccumEmptyDet(CurrentDets(:, i))) AnnihilatedDet = AnnihilatedDet + 1

                    ! count the number of walkers that are single-spawns at the threshold
                    if (t_prone_walkers) then
                        if (test_flag(CurrentDets(:, i), flag_prone)) n_prone_dets = n_prone_dets + 1
                    end if

                    if (tEScaleWalkers) then
                        scaledOccupiedThresh = OccupiedThresh * scaleFunction(det_diagH(i))
                        scaledOccupiedThresh = OccupiedThresh
                    end if
                    do j = 1, lenof_sign
                        run = part_type_to_run(j)
                        if (.not. tIsStateDeterm) then
                            if ((abs(CurrentSign(j)) > 1.e-12_dp) .and. (abs(CurrentSign(j)) < scaledOccupiedThresh)) then
                                !We remove this walker with probability 1-RealSignTemp
                                pRemove = (scaledOccupiedThresh - abs(CurrentSign(j))) / scaledOccupiedThresh
                                r = genrand_real2_dSFMT()
                                if (pRemove > r) then
                                    !Remove this walker
                                    NoRemoved(run) = NoRemoved(run) + abs(CurrentSign(j))
                                    iter_data%nremoved(j) = iter_data%nremoved(j) &
                                                            + abs(CurrentSign(j))
                                    CurrentSign(j) = 0.0_dp
                                    call nullify_ilut_part(CurrentDets(:, i), j)
                                    call decode_bit_det(nI, CurrentDets(:, i))
                                    if (IsUnoccDet(CurrentSign) .and. .not. tAccumEmptyDet(CurrentDets(:, i))) then
                                        call RemoveHashDet(HashIndex, nI, i)
                                        ! also update both the number of annihilated dets
                                        AnnihilatedDet = AnnihilatedDet + 1
                                        ! and the number of holes
                                        HolesInList = HolesInList + 1
                                    end if
                                    NoBorn(run) = NoBorn(run) + scaledOccupiedThresh - abs(CurrentSign(j))
                                    iter_data%nborn(j) = iter_data%nborn(j) &
                                                         + scaledOccupiedThresh - abs(CurrentSign(j))
                                    CurrentSign(j) = sign(scaledOccupiedThresh, CurrentSign(j))
                                    call encode_part_sign(CurrentDets(:, i), CurrentSign(j), j)
                                end if
                            end if
                        end if
                    end do

                    TotParts = TotParts + abs(CurrentSign)

                    call addNormContribution(CurrentSign, tIsStateDeterm)

                    if (tCheckHighestPop) then
                        ! If this option is on, then we want to compare the
                        ! weight on each determinant to the weight at the HF
                        ! determinant.
                        ! Record the highest weighted determinant on each
                        ! processor. If double run, only consider set 1 to keep things simple.
                        do run = 1, inum_runs
                            lbnd = min_part_type(run)
                            ubnd = max_part_type(run)
                            if (abs_sign(CurrentSign(lbnd:ubnd)) > iHighestPop(run)) then
                                iHighestPop(run) = int(abs_sign(CurrentSign(lbnd:ubnd)))
                                HighestPopDet(:, run) = CurrentDets(:, i)
                            end if
                        end do
                    end if
                end if

                if (tFillingStochRDMonFly) then
                    if (IsUnoccDet(CurrentSign) .and. (.not. tIsStateDeterm)) then
                        if (DetBitEQ(CurrentDets(:, i), iLutHF_True, nifd)) then
                            ! We have to do this such that AvNoAtHF matches up with AvSign.
                            ! AvSign is extracted from CurrentH, and if the HFDet is unoccupied
                            ! at this moment during annihilation, it's CurrentH entry is removed
                            ! and the averaging information in it is lost.
                            ! In some cases (a successful spawning event) a CurrentH entry will
                            ! be recreated, but with AvSign 0, so we must match this here.
                            AvNoAtHF = 0.0_dp
                            IterRDM_HF = Iter + 1
                        end if
                    end if
                end if

                ! This InstNoAtHF call must be placed at the END of the routine
                ! as the value of CurrentSign can change during it!
                if (DetBitEQ(CurrentDets(:, i), iLutHF_True, nifd)) then
                    InstNoAtHF = CurrentSign
                end if

            end do
        end if

        IFDEBUGTHEN(FCIMCDebug, 6)
            write(stdout, *) "After annihilation: "
            write(stdout, *) "TotWalkersNew: ", TotWalkersNew
            write(stdout, *) "AnnihilatedDet: ", AnnihilatedDet
            write(stdout, *) "HolesInList: ", HolesInList
            write(stdout, "(A,I12)") "Walker list length: ", TotWalkersNew
            write(stdout, "(A)") "TW: Walker  Det"
            do j = 1, int(TotWalkersNew)
                CurrentSign = &
                    transfer(CurrentDets(IlutBits%ind_pop:IlutBits%ind_pop + lenof_sign - 1, j), &
                write(stdout, "(A,I10,a)", advance='no') 'TW:', j, '['
                do part_type = 1, lenof_sign
                    write(stdout, "(f16.3)", advance='no') CurrentSign(part_type)
                end do
                call WriteBitDet(stdout, CurrentDets(:, j), .true.)
                call neci_flush(stdout)
            end do

        ! RT_M_Merge: i have to ask werner why this check makes sense
        ! AnnihilatedDet is only affected by empty dets and emptying a det increses HolesInList
        ! But adding a new det decreases HolesInList and does not affect AnnihilatedDet ->?
        if (AnnihilatedDet /= HolesInList) then
            write(stdout, *) "TotWalkersNew: ", TotWalkersNew
            write(stdout, *) "AnnihilatedDet: ", AnnihilatedDet
            write(stdout, *) "HolesInList: ", HolesInList
            write(stdout, *) "iStartFreeSlot, iEndFreeSlot:", iStartFreeSlot, iEndFreeSlot
            write(stdout, *) "TotParts: ", TotParts
            call neci_flush(stdout)
            call stop_all(t_r, "Error in determining annihilated determinants")
        end if
    end subroutine CalcHashTableStats