cont_time_fcimc.F90 Source File


Source Code

Source Code

#include "macros.h"
module cont_time

    use CalcData, only: tContTimeFull, tOrthogonaliseReplicas, &
                        orthogonalise_iter, use_spawn_hash_table, &
                        tTruncInitiator, DiagSft, tPairedReplicas
    use fcimc_helper, only: rezero_iter_stats_each_iter, CalcParentFlag, &
                            create_particle, create_particle_with_hash_table, &
                            SumEContrib, end_iter_stats
    use cont_time_rates, only: spawn_rate_full, cont_time_gen_excit_full, &
                               oversample_factors, cont_time_gen_excit, ostag, &
    use hash, only: remove_hash_table_entry, clear_hash_table
    use DetBitOps, only: FindBitExcitLevel, count_open_orbs
    use global_det_data, only: det_diagH, det_offdiagH, get_spawn_rate
    use GenRandSymExcitNUMod, only: init_excit_gen_store
    use Determinants, only: get_helement
    use DeterminantData, only: write_det
    use orthogonalise, only: orthogonalise_replicas
    use dSFMT_interface, only: genrand_real2_dSFMT
    use AnnihilationMod, only: DirectAnnihilation, communicate_and_merge_spawns
    use bit_reps, only: nullify_ilut, encode_sign
    use fcimc_iter_utils, only: update_iter_data
    use hphf_integrals, only: hphf_diag_helement
    use SystemData, only: nel, tHPHF, LMS
    use bit_reps, only: extract_bit_rep, writebitdet
    use LoggingData, only: FCIMCDebug
    use bit_rep_data, only: NIfTot
    use rdm_data, only: rdm_definitions
    use FciMCData
    use constants
    use util_mod
    use tau_main, only: tau
    implicit none


    ! TODO: Set the initial spawning rate of the first particle (In set global
    !       data).

    subroutine iterate_cont_time(iter_data)

        type(fcimc_iter_data), intent(inout) :: iter_data
        character(*), parameter :: this_routine = 'iterate_cont_time'

        real(dp) :: sgn(lenof_sign), rate, hdiag
        HElement_t(dp) :: hoffdiag
        integer :: sgn_abs, iunused, flags, det(nel), j, p, TotWalkersNew

        integer :: part_type, ic_hf, nopen, err, MaxIndex
        logical :: survives

        if (lenof_sign /= 1) then
            call stop_all(this_routine, "Complex and multi-run calculations &
                         &not supported using CONT-TIME")
        end if

        ! Reset free position list
        ValidSpawnedList = InitialSpawnedSlots
        FreeSlot(1:iEndFreeSlot) = 0
        iStartFreeSlot = 1
        iEndFreeSlot = 0
        if (use_spawn_hash_table) call clear_hash_table(spawn_ht)

        call rezero_iter_stats_each_iter(iter_data, rdm_definitions)
        call set_timer(walker_time)

        ! Loop over the main list
        do j = 1, int(TotWalkers)

            fcimc_excit_gen_store%tFilled = .false.
            call extract_bit_rep(CurrentDets(:, j), det, sgn, flags, j, &

            ! As the main list is not contiguous, skip (but store) empty sites.
            if (IsUnoccDet(sgn)) then
                iEndFreeSlot = iEndFreeSlot + 1
                FreeSlot(iEndFreeSlot) = j
            end if

            IFDEBUG(FCIMCDebug, 3) then
                write(stdout, "(A,I10,a)", advance='no') 'TW:', j, '['
                do part_type = 1, lenof_sign
                    write(stdout, "(f10.5)", advance='no') sgn(part_type)
                end do
                write(stdout, '(a)', advance='no') '] '
                call WriteBitDet(stdout, CurrentDets(:, j), .true.)
                call neci_flush(stdout)
            end if

            ! Global stored data to make things efficient
            hdiag = det_diagH(j)
            hoffdiag = det_offdiagH(j)
            if (tContTimeFull) then
                rate = get_spawn_rate(j)
                ASSERT(rate.isclose.spawn_rate_full(det, CurrentDets(:, j)))
            end if

            ! Calculate the flags that ought to be carried through
            if (tTruncInitiator) then
                call CalcParentFlag(j, iunused)
            end if

            ! Sum in the energy terms, yeah!
            ic_hf = FindBitExcitLevel(ilutRef(:, 1), CurrentDets(:, j))
            call SumEContrib(det, ic_hf, sgn, CurrentDets(:, j), hdiag, &
                             hoffdiag, 1.0_dp, tPairedReplicas, j)

            ! Needed for calculating oversample factors
            nopen = count_open_orbs(CurrentDets(:, j))

            ! Loop over determinants, and the particles on the determinant
            do part_type = 1, lenof_sign

                ! n.b. only using integer particles atm.
                sgn_abs = int(abs(sgn(part_type)))
                do p = 1, sgn_abs

                    survives = process_part_cont_time( &
                               CurrentDets(:, j), det, sgn(part_type), &
                               part_type, rate, totimagtime, &
                               totimagtime + tau, hdiag, iter_data, &
                               nopen, fcimc_excit_gen_store)

                    if (.not. survives) then
                        iter_data%ndied = iter_data%ndied + 1
                        NoDied = NoDied + 1
                        sgn(part_type) = sgn(part_type) &
                                         - sign(1.0_dp, sgn(part_type))
                    end if
                end do
            end do

            ! If this particle has been completely destroyed, then remove it
            ! from the simulation
            if (all(near_zero(sgn))) then
                call remove_hash_table_entry(HashIndex, det, j)

                ! And add to the "freeslot" list
                iEndFreeSlot = iEndFreeSlot + 1
                FreeSlot(iEndFreeSlot) = j
                call nullify_ilut(CurrentDets(:, j))
                call encode_sign(CurrentDets(:, j), sgn)
            end if

        end do
        IFDEBUG(FCIMCDebug, 2) write(stdout, '("Finnished loop over sites")')
        call halt_timer(walker_time)

        ! Update statistics. This is done before annihilation as the output
        ! statistics normally apply to the iteration _before_ the one that has
        ! just run (as they are accumulated during the normal loop)
        TotWalkersNew = int(TotWalkers)
        call end_iter_stats(TotWalkersNew)

        ! Send walkers to the correct nodes, and annihilate
        call set_timer(annihil_time)
        call communicate_and_merge_spawns(MaxIndex, iter_data, .false.)
        call DirectAnnihilation(TotWalkersNew, MaxIndex, iter_data, err)

        TotWalkers = TotWalkersNew
        call halt_timer(annihil_time)
        IFDEBUG(FCIMCDebug, 2) write(stdout, '("Finished annihilation")')

        ! If we are orthogonalising the replica wavefunctions, to generate
        ! excited states, then do that here.
        if (tOrthogonaliseReplicas .and. iter > orthogonalise_iter) then
            call orthogonalise_replicas(iter_data)
        end if

        ! Update iteration counters
        call update_iter_data(iter_data)

    end subroutine

    recursive function process_part_cont_time(ilut, det, sgn, part_type, &
                                              rate, curr_time, annihil_time, &
                                              hdiag, iter_data, nopen, &
                                              store) &

        integer(n_int), intent(in) :: ilut(0:NifTot)
        integer, intent(in) :: det(nel), part_type, nopen
        real(dp), intent(in) :: sgn, curr_time, annihil_time, rate, hdiag
        type(fcimc_iter_data), intent(inout) :: iter_data
        logical :: survives
        type(excit_gen_store_type), intent(inout) :: store
        character(*), parameter :: this_routine = 'process_part_cont_time'

        real(dp) :: time, child(lenof_sign), dt, rate_adj, rate_spwn
        real(dp) :: hdiag_spwn, spwn_sgn
        integer :: nspawn, spawn_sgn, det_spwn(nel), ic, i, y, nopen_spwn, err
        integer(n_int) :: ilut_spwn(0:NIfTot)
        logical :: child_survives
        HElement_t(dp) :: hoffdiag, htmp, hdiag_spawn

        ! A quick sanity check that we have calculated the spawning rate
        ! reasonably
        ASSERT(rate > 0 .or. .not. tContTimeFull)

        ! The starting time for the spawning
        time = curr_time

        survives = .true.
        do while (.true.)

            ! Calculate the _live_ spawning rate inside the particle loop so
            ! that the live rate gets updated as soon as the oversampling
            ! factors are tweaked
            if (tContTimeFull) then
                rate_adj = rate + abs(hdiag - DiagSft(part_type))
                rate_adj = sum(oversample_factors(:, nopen)) &
                           + abs(hdiag - DiagSft(part_type))
            end if

            ! Generate the next spawning event
            dt = -1.0_dp * log(genrand_real2_dSFMT()) / rate_adj
            time = time + dt
            if (time > annihil_time) exit

            if (tContTimeFull) then
                nspawn = 1
                call cont_time_gen_excit_full(det, ilut, rate_adj, hdiag, &
                                              det_spwn, ilut_spwn, &
                                              hoffdiag, ic, part_type)
                call cont_time_gen_excit(det, ilut, rate_adj, hdiag, &
                                         det_spwn, ilut_spwn, hoffdiag, &
                                         ic, part_type, nopen, nspawn, store)
            end if

            if (.not. IsNullDet(det_spwn)) then

                ! If we generate the starting det, this is the same as particle
                ! death in normal FCIQMC. Simplify everything by just killing
                ! the current particle, rather than generating antiparticles.
                if (ic == 0) then
                    survives = .false.
                    IFDEBUG(FCIMCDebug, 3) write(stdout, '("Particle died")')
                end if

                ! We need the diagonal matrix element for calculating further
                ! spawning rates
                if (tHPHF) then
                    htmp = hphf_diag_helement(det_spwn, ilut_spwn)
                    htmp = get_helement(det_spwn, det_spwn, 0)
                end if
                hdiag_spwn = real(htmp, dp) - Hii

                ! Calculate full rate if needed
                if (tContTimeFull) then
                    rate_spwn = spawn_rate_full(det_spwn, ilut_spwn)
                end if

                ! Particle accountancy
                iter_data%nborn = iter_data%nborn + nspawn
                NoBorn = NoBorn + nspawn

                ! This is the spawning coefficient. n.b. it is possible to
                ! spawn more than one particle whilst adjusting the overspawn
                ! factors
#ifndef CMPLX_
                spwn_sgn = -sign(1.0_dp, sgn) * sign(real(nspawn, dp), hoffdiag)
                call stop_all(this_routine, "Not implemented")

                ! Need this for calculating further oversampling factors
                nopen_spwn = count_open_orbs(ilut_spwn)

                ! We can't store data in this structure recursively...
                secondary_gen_store%tFilled = .false.

                do i = 1, nspawn

                    child_survives = process_part_cont_time( &
                                     ilut_spwn, det_spwn, spwn_sgn, &
                                     part_type, rate_spwn, time, &
                                     annihil_time, hdiag_spwn, iter_data, &
                                     nopen_spwn, secondary_gen_store)

                    if (.not. child_survives) then
                        iter_data%ndied = iter_data%ndied + 1
                        NoDied = NoDied + 1
                        spwn_sgn = spwn_sgn - sign(1.0_dp, spwn_sgn)
                    end if
                end do

                ! And now create the generated particle (if it has survived
                ! this far)
                ! TODO: Ensure that the child iluts passed into this routine
                !       don't have the initiator flag set
                if (.not. near_zero(spwn_sgn)) then
                    child = 0
                    child(part_type) = spwn_sgn

                    IFDEBUG(FCIMCDebug, 3) then
                        write(stdout, '(a)', advance='no') 'SP: ['
                        do y = 1, lenof_sign
                            write(stdout, '(f12.5)', advance='no') &
                        end do
                        write(stdout, '("] ")', advance='no')
                        call write_det(stdout, det_spwn, .true.)
                        call neci_flush(stdout)
                    end if

                    if (use_spawn_hash_table) then
                        call create_particle_with_hash_table( &
                            det_spwn, ilut_spwn, child, &
                            part_type, ilut, iter_data, err)
                        call create_particle(det_spwn, ilut_spwn, child, &
                                             part_type, hdiag_spawn, err, ilut)
                    end if
                end if

            end if ! Spawn not aborted

        end do ! Loop until annihil time

    end function

end module