perform_kp_fciqmc Subroutine

public subroutine perform_kp_fciqmc(kp)


Type IntentOptional Attributes Name
type(kp_fciqmc_data), intent(inout) :: kp


Source Code

Source Code

    subroutine perform_kp_fciqmc(kp)

        use orthogonalise, only: calc_replica_overlaps

        type(kp_fciqmc_data), intent(inout) :: kp
        integer :: iiter, idet, ireplica, ispawn, ierr, err
        integer :: iconfig, irepeat, ivec, nlowdin
        integer :: nspawn, parent_flags, unused_flags
        integer :: ex_level_to_ref, ex_level_to_hf
        integer :: TotWalkersNew, determ_ind, ic, ex(2, maxExcit), ms_parent
        integer :: nI_parent(nel), nI_child(nel), MaxIndex
        integer(n_int) :: ilut_child(0:NIfTot)
        integer(n_int), pointer :: ilut_parent(:)
        real(dp) :: prob, unused_rdm_real, parent_hdiag
        real(dp) :: child_sign(lenof_sign), parent_sign(lenof_sign)
        real(dp) :: unused_sign(lenof_sign), precond_fac
        real(dp), allocatable :: lowdin_evals(:, :)
        logical :: tChildIsDeterm, tParentIsDeterm, tParentUnoccupied
        logical :: tParity, tSingBiasChange, tWritePopsFound
        HElement_t(dp) :: HElGen, parent_hoffdiag

        ! Stores of the overlap and projected Hamiltonian matrices.
        real(dp), pointer :: overlap_matrices(:, :, :)
        real(dp), pointer :: hamil_matrices(:, :, :)
        ! Pointers to the matrices for a given repeat only.
        real(dp), pointer :: overlap_matrix(:, :)
        real(dp), pointer :: hamil_matrix(:, :)

        ! Variables to hold information output for the test suite.
        real(dp) :: s_sum, h_sum

        integer(n_int) :: int_sign(lenof_all_signs)
        real(dp) :: test_sign(lenof_all_signs)
        type(ll_node), pointer :: temp_node

        ! Unused factor
        precond_fac = 1.0_dp

        call init_kp_fciqmc(kp)
        if (.not. tAllSymSectors) ms_parent = lms

        allocate(overlap_matrices(kp%nvecs, kp%nvecs, kp%nrepeats), stat=ierr)
        allocate(hamil_matrices(kp%nvecs, kp%nvecs, kp%nrepeats), stat=ierr)
        allocate(lowdin_evals(kp%nvecs, kp%nvecs), stat=ierr)
        overlap_matrices = 0.0_dp
        hamil_matrices = 0.0_dp

        outer_loop: do iconfig = 1, kp%nconfigs

            do irepeat = 1, kp%nrepeats

                ! Point to the regions of memory where the projected Hamiltonian
                ! and overlap matrices for this repeat will be accumulated and stored.
                overlap_matrix => overlap_matrices(:, :, irepeat)
                hamil_matrix => hamil_matrices(:, :, irepeat)
                overlap_matrix(:, :) = 0.0_dp
                hamil_matrix(:, :) = 0.0_dp

                call init_kp_fciqmc_repeat(iconfig, irepeat, kp%nrepeats, kp%nvecs, iter_data_fciqmc)

                do ivec = 1, kp%nvecs

                    ! Copy the current state of CurrentDets to krylov_vecs.
                    call store_krylov_vec(ivec)

                    ! Calculate the overlap of the perturbed ground state vector
                    ! with the new Krylov vector, if requested.
                    if (tOverlapPert) call calc_perturbation_overlap(ivec)

                    do iiter = 1, kp%niters(ivec)

                        call set_timer(walker_time)

                        iter = iter + 1
                        call init_kp_fciqmc_iter(iter_data_fciqmc, determ_ind)

                        !if (iter < 10) then
                        !    write(stdout,*) "CurrentDets before:"
                        !    do idet = 1, int(TotWalkers)
                        !        call extract_bit_rep(CurrentDets(:, idet), nI_parent, parent_sign, unused_flags, &
                        !                              fcimc_excit_gen_store)
                        !        if (tUseFlags) then
                        !            write(stdout,'(i7, i12, 4x, f18.7, 4x, f18.7, 4x, l1)') idet, CurrentDets(0,idet), parent_sign, &
                        !                test_flag(CurrentDets(:,idet), flag_deterministic)
                        !        else
                        !            write(stdout,'(i7, i12, 4x, f18.7, 4x, f18.7)') idet, CurrentDets(0,idet), parent_sign
                        !        end if
                        !    end do

                        !    write(stdout,"(A)") "Hash Table: "
                        !    do idet = 1, nWalkerHashes
                        !        temp_node => HashIndex(idet)
                        !        if (temp_node%ind /= 0) then
                        !            write(stdout,'(i9)',advance='no') idet
                        !            do while (associated(temp_node))
                        !                write(stdout,'(i9)',advance='no') temp_node%ind
                        !                temp_node => temp_node%next
                        !            end do
                        !            write(stdout,'()',advance='yes')
                        !        end if
                        !    end do
                        !end if

                        do idet = 1, int(TotWalkers)

                            ! The 'parent' determinant from which spawning is to be attempted.
                            ilut_parent => CurrentDets(:, idet)
                            parent_flags = 0_n_int

                            ! Indicate that the scratch storage used for excitation generation from the
                            ! same walker has not been filled (it is filled when we excite from the first
                            ! particle on a determinant).
                            fcimc_excit_gen_store%tFilled = .false.

                            call extract_bit_rep(ilut_parent, nI_parent, parent_sign, unused_flags, &
                                                 idet, fcimc_excit_gen_store)

                            ex_level_to_ref = FindBitExcitLevel(iLutRef(:, 1), ilut_parent, &
                            if (tRef_Not_HF) then
                                ex_level_to_hf = FindBitExcitLevel(iLutHF_true, ilut_parent, max_calc_ex_level)
                                ex_level_to_hf = ex_level_to_ref
                            end if

                            tParentIsDeterm = check_determ_flag(ilut_parent, core_run)
                            tParentUnoccupied = IsUnoccDet(parent_sign)

                            ! If this determinant is in the deterministic space then store the relevant
                            ! data in arrays for later use.
                            if (tParentIsDeterm) then
                                ! Store the index of this state, for use in annihilation later.
                                associate(rep => cs_replicas(core_run))
                                    rep%indices_of_determ_states(determ_ind) = idet
                                    ! Add the amplitude to the deterministic vector.
                                    rep%partial_determ_vecs(:, determ_ind) = parent_sign
                                end associate
                                determ_ind = determ_ind + 1

                                ! The deterministic states are always kept in CurrentDets, even when
                                ! the amplitude is zero. Hence we must check if the amplitude is zero
                                ! and, if so, skip the state.
                                if (tParentUnoccupied) cycle
                            end if

                            ! If this slot is unoccupied (and also not a core determinant) then add it to
                            ! the list of free slots and cycle.
                            if (tParentUnoccupied) then
                                iEndFreeSlot = iEndFreeSlot + 1
                                FreeSlot(iEndFreeSlot) = idet
                            end if

                            ! The current diagonal matrix element is stored persistently.
                            parent_hdiag = det_diagH(idet)
                            parent_hoffdiag = det_offdiagH(idet)

                            if (tTruncInitiator) call CalcParentFlag(idet, parent_flags)

                            call SumEContrib(nI_parent, ex_level_to_ref, parent_sign, ilut_parent, &
                                             parent_hdiag, parent_hoffdiag, 1.0_dp, tPairedReplicas, idet)

                            ! If we're on the Hartree-Fock, and all singles and
                            ! doubles are in the core space, then there will be
                            ! no stochastic spawning from this determinant, so
                            ! we can the rest of this loop.
                            if (ss_space_in%tDoubles .and. ex_level_to_hf == 0 .and. tDetermHFSpawning) cycle

                            if (tAllSymSectors) then
                                ms_parent = return_ms(ilut_parent)
                                nOccAlpha = (nel + ms_parent) / 2
                                nOccBeta = (nel - ms_parent) / 2
                            end if

                            ! If this condition is not met (if all electrons have spin up or all have spin down)
                            ! then there will be no determinants to spawn to, so don't attempt spawning.
                            if (abs(ms_parent) /= nbasis / 2) then

                                do ireplica = 1, lenof_sign

                                    call decide_num_to_spawn(parent_sign(ireplica), AvMCExcits, nspawn)

                                    do ispawn = 1, nspawn

                                        ! Zero the bit representation, to ensure no extraneous data gets through.
                                        ilut_child = 0_n_int

                                        call generate_excitation(nI_parent, ilut_parent, nI_child, &
                                                                 ilut_child, exFlag, ic, ex, tParity, prob, &
                                                                 HElGen, fcimc_excit_gen_store)

                                        ! If a valid excitation.
                                        if (.not. IsNullDet(nI_child)) then

                                            call encode_child(ilut_parent, ilut_child, ic, ex)
                                            ilut_child(IlutBits%ind_flag) = 0_n_int

                                            if (tSemiStochastic) then
                                                tChildIsDeterm = is_core_state(ilut_child, nI_child)

                                                ! Is the parent state in the core space?
                                                if (tParentIsDeterm) then
                                                    ! If spawning is both from and to the core space, cancel it.
                                                    if (tChildIsDeterm) cycle
                                                    call set_flag(ilut_child, flag_determ_parent)
                                                    if (tChildIsDeterm) call set_flag(ilut_child, flag_deterministic(core_run))
                                                end if
                                            end if

                                            child_sign = attempt_create(nI_parent, ilut_parent, parent_sign, &
                                                                        nI_child, ilut_child, prob, HElGen, ic, ex, tParity, &
                                                                        ex_level_to_ref, ireplica, unused_sign, &
                                                                        AvMCExcits, unused_rdm_real, precond_fac)

                                            child_sign = 0.0_dp
                                        end if

                                        ! If any (valid) children have been spawned.
                                        if (.not. (all(near_zero(child_sign) .or. ic == 0 .or. ic > 2))) then

                                            call new_child_stats_normal(iter_data_fciqmc, ilut_parent, &
                                                                 nI_child, ilut_child, ic, ex_level_to_ref, &
                                                                 child_sign, parent_flags, ireplica)

                                            call create_particle_with_hash_table(nI_child, ilut_child, child_sign, &
                                                                                 ireplica, ilut_parent, iter_data_fciqmc, ierr)

                                        end if ! If a child was spawned.

                                    end do ! Over mulitple particles on same determinant.

                                end do ! Over the replicas on the same determinant.

                            end if ! If connected determinants exist to spawn to.

                            ! If this is a core-space determinant then the death step is done in
                            ! determ_projection.
                            if (.not. tParentIsDeterm) then
                                call walker_death(iter_data_fciqmc, nI_parent, ilut_parent, parent_hdiag, &
                                                  parent_sign, idet, ex_level_to_ref)
                            end if

                        end do ! Over all determinants.

                        if (tSemiStochastic) call determ_projection()

                        TotWalkersNew = int(TotWalkers)
                        call end_iter_stats(TotWalkersNew)
                        call end_iteration_print_warn(TotWalkersNew)

                        call halt_timer(walker_time)

                        call set_timer(annihil_time)

                        call communicate_and_merge_spawns(MaxIndex, iter_data_fciqmc, .false.)

                        call DirectAnnihilation(TotWalkersNew, MaxIndex, iter_data_fciqmc, ierr)

                        TotWalkers = int(TotWalkersNew, int64)

                        call halt_timer(annihil_time)

                        if (tPrintReplicaOverlaps) call calc_replica_overlaps()

                        call update_iter_data(iter_data_fciqmc)

                        if (mod(iter, StepsSft) == 0) then
                            call set_timer(Stats_Comms_Time)
                            call calculate_new_shift_wrapper(iter_data_fciqmc, TotParts, tPairedReplicas)
                            call halt_timer(Stats_Comms_Time)

                            call ChangeVars(tSingBiasChange, tWritePopsFound)
                            if (tWritePopsFound) call WriteToPopsfileParOneArr(CurrentDets, TotWalkers)
                            if (tSingBiasChange) call CalcApproxpDoubles()
                            if (tSoftExitFound) exit outer_loop
                        end if

                    end do ! Over all iterations between Krylov vectors.

                end do ! Over all Krylov vectors.

                call calc_overlap_matrix(kp%nvecs, krylov_vecs, TotWalkersKP, overlap_matrix)

                if (tExactHamil) then
                    call calc_hamil_exact(kp%nvecs, krylov_vecs, TotWalkersKP, hamil_matrix, krylov_helems)
                    if (tSemiStochastic) then
                        call calc_projected_hamil(kp%nvecs, krylov_vecs, krylov_vecs_ht, TotWalkersKP, hamil_matrix, &
                                                  partial_determ_vecs_kp, full_determ_vecs_kp, krylov_helems)
                        call calc_projected_hamil(kp%nvecs, krylov_vecs, krylov_vecs_ht, TotWalkersKP, hamil_matrix, &
                    end if
                end if

                ! Sum the overlap and projected Hamiltonian matrices from the various processors.
                call communicate_kp_matrices(overlap_matrix, hamil_matrix)

                if (iProcIndex == root) call output_kp_matrices_wrapper(iconfig, overlap_matrices, hamil_matrices)

            end do ! Over all repeats for a fixed initial walker configuration.

            if (tOverlapPert) call average_and_comm_pert_overlaps(kp%nrepeats)

            if (iProcIndex == root) then
                call average_kp_matrices_wrapper(iconfig, kp%nrepeats, overlap_matrices, hamil_matrices, &
                                                 kp_overlap_mean, kp_hamil_mean, kp_overlap_se, kp_hamil_se)
                call find_and_output_lowdin_eigv(iconfig, kp%nvecs, kp_overlap_mean, kp_hamil_mean, nlowdin, lowdin_evals, .true.)

                ! Calculate data for the testsuite.
                s_sum = sum(kp_overlap_mean)
                h_sum = sum(kp_hamil_mean)
            end if

        end do outer_loop ! Over all initial walker configurations.


        if (tPopsFile) call WriteToPopsfileParOneArr(CurrentDets, TotWalkers)

        if (iProcIndex == root) call write_kpfciqmc_testsuite_data(s_sum, h_sum)

    end subroutine perform_kp_fciqmc