perform_subspace_fciqmc Subroutine

public subroutine perform_subspace_fciqmc(kp)


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


Source Code

    subroutine perform_subspace_fciqmc(kp)

        use fcimc_helper, only: create_particle_with_hash_table
        use FciMCData, only: HashIndex, nWalkerHashes
        use orthogonalise, only: orthogonalise_replicas, orthogonalise_replica_pairs
        use orthogonalise, only: calc_replica_overlaps

        type(kp_fciqmc_data), intent(inout) :: kp

        integer :: iiter, idet, ireplica, ispawn, ierr, err
        integer :: iconfig, irepeat, ireport, nlowdin
        integer :: nspawn, parent_flags, unused_flags
        integer :: ex_level_to_ref, ex_level_to_hf
        integer :: TotWalkersNew, determ_ind, ic, ex(2, maxExcit), MaxIndex
        integer :: nI_parent(nel), nI_child(nel), unused_vecslot
        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(:, :), lowdin_spin(:, :)
        logical :: tChildIsDeterm, tParentIsDeterm, tParentUnoccupied
        logical :: tParity, tSingBiasChange, tWritePopsFound
        HElement_t(dp) :: HElGen, parent_hoffdiag

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

        type(ll_node), pointer :: temp_node

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

        ! Unused factor
        precond_fac = 1.0_dp

        call init_kp_fciqmc(kp)

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

        if (tCalcSpin) then
            allocate(spin_matrices(kp%nvecs, kp%nvecs, kp%nrepeats, kp%nreports), stat=ierr)
            allocate(lowdin_spin(kp%nvecs, kp%nvecs), stat=ierr)
            spin_matrices = 0.0_dp
            lowdin_spin = 0.0_dp
        end if

        outer_loop: do irepeat = 1, kp%nrepeats

            call init_kp_fciqmc_repeat(iconfig, irepeat, kp%nrepeats, kp%nvecs, iter_data_fciqmc)
            if (iProcIndex == root) call write_ex_state_header(kp%nvecs, irepeat)

            do ireport = 1, kp%nreports

                ! 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, ireport)
                hamil_matrix => hamil_matrices(:, :, irepeat, ireport)
                overlap_matrix(:, :) = 0.0_dp
                hamil_matrix(:, :) = 0.0_dp
                if (tCalcSpin) then
                    spin_matrix => spin_matrices(:, :, irepeat, ireport)
                    spin_matrix(:, :) = 0.0_dp
                end if

                call calc_overlap_matrix(kp%nvecs, CurrentDets, int(TotWalkers), overlap_matrix)

                if (tExactHamil) then
                    call calc_hamil_exact(kp%nvecs, CurrentDets, int(TotWalkers), hamil_matrix)
                    if (tSemiStochastic) then
                        associate(rep => cs_replicas(core_run))
                            call calc_projected_hamil(kp%nvecs, CurrentDets, HashIndex, int(TotWalkers), &
                                                      hamil_matrix, rep%partial_determ_vecs, rep%full_determ_vecs)
                        end associate
                        call calc_projected_hamil(kp%nvecs, CurrentDets, HashIndex, int(TotWalkers), &
                    end if
                end if

                !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, 3(f18.7, 4x), l1)') idet, CurrentDets(0,idet), parent_sign, &
                !            test_flag(CurrentDets(:,idet), flag_deterministic)
                !    else
                !        write(stdout,'(i7, i12, 3(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

                ! Sum the overlap and projected Hamiltonian matrices from the various processors.
                if (tCalcSpin) then
                    ! Calculate the spin squared projected into the subspace.
                    call calc_projected_spin(kp%nvecs, CurrentDets, HashIndex, int(TotWalkers), spin_matrix)
                    call communicate_kp_matrices(overlap_matrix, hamil_matrix, spin_matrix)
                    call communicate_kp_matrices(overlap_matrix, hamil_matrix)
                end if

                if (iProcIndex == root) then
                    call output_kp_matrices_wrapper(iter, overlap_matrices(:, :, 1:irepeat, ireport), &
                                                    hamil_matrices(:, :, 1:irepeat, ireport))
                    if (tCalcSpin) then
                        call find_and_output_lowdin_eigv(iter, kp%nvecs, overlap_matrix, hamil_matrix, nlowdin, &
                                                         lowdin_evals, .false., spin_matrix, lowdin_spin)
                        call write_ex_state_data(iter, nlowdin, lowdin_evals, hamil_matrix, overlap_matrix, &
                                                 spin_matrix, lowdin_spin)
                        call find_and_output_lowdin_eigv(iter, kp%nvecs, overlap_matrix, hamil_matrix, nlowdin, &
                                                         lowdin_evals, .false.)
                        call write_ex_state_data(iter, nlowdin, lowdin_evals, hamil_matrix, overlap_matrix)
                    end if
                end if

                do iiter = 1, kp%niters(ireport)

                    call set_timer(walker_time)

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

                    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

                        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.

                        ! 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 (tOrthogKPReplicas .and. iter > orthog_kp_iter) then
                        call orthogonalise_replicas(iter_data_fciqmc)
                    else if (tPrintReplicaOverlaps) then
                        call calc_replica_overlaps()
                    end if

                    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 per report cycle.

            end do ! Over all report cycles.

        end do outer_loop ! Over all repeats of the whole calculation.

        ! Output the Lowdin estimates for the final *averaged* matrices.
        if (iProcIndex == root) then
            iter = 0
            do ireport = 1, kp%nreports
                call average_kp_matrices_wrapper(iter, kp%nrepeats, overlap_matrices(:, :, 1:kp%nrepeats, ireport), &
                                                 hamil_matrices(:, :, 1:kp%nrepeats, ireport), kp_overlap_mean, &
                                                 kp_hamil_mean, kp_overlap_se, kp_hamil_se)
                if (tCalcSpin) then
                    call find_and_output_lowdin_eigv(iter, kp%nvecs, kp_overlap_mean, kp_hamil_mean, nlowdin, &
                                                     lowdin_evals, .true., spin_matrix, lowdin_spin)
                    call find_and_output_lowdin_eigv(iter, kp%nvecs, kp_overlap_mean, kp_hamil_mean, nlowdin, &
                                                     lowdin_evals, .true.)
                end if
                ! Update the iteration label.
                iter = iter + kp%niters(ireport)
            end do
        end if

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


        if (tPopsFile) call WriteToPopsfileParOneArr(CurrentDets, TotWalkers)

        call write_kpfciqmc_testsuite_data(s_sum, h_sum)

    end subroutine perform_subspace_fciqmc