init_semi_stochastic Subroutine

public subroutine init_semi_stochastic(core_in, tStartedFromCoreGround)


Type IntentOptional Attributes Name
type(subspace_in) :: core_in
logical, intent(out) :: tStartedFromCoreGround


Source Code

Source Code

    subroutine init_semi_stochastic(core_in, tStartedFromCoreGround)

        ! Initialise the semi-stochastic information. This includes enumerating a list of all
        ! determinants or CSFs in the deterministic space and calculating and storing the resulting
        ! Hamiltonian matrix elements. The lists which will store the walker amplitude vectors in
        ! the deterministic space are also allocated.

        type(subspace_in) :: core_in
        logical, intent(out) :: tStartedFromCoreGround

        integer :: i, ierr, run, num_core_runs
        integer :: nI(nel)
        integer(MPIArg) :: mpi_temp
        character(len=*), parameter :: t_r = "init_semi_stochastic"
#ifndef CMPLX_
        real(dp), allocatable :: e_values(:)
        HElement_t(dp), allocatable :: e_vectors(:, :), gs_vector(:)
        real(dp) :: gs_energy
        ! If we are load balancing, this gets disabled once semi stochastic
        ! has been initialised. Therefore we should do a last-gasp load
        ! adjustment at this point.
        if (tLoadBalanceBlocks .and. .not. tFillingStochRDMOnFly) then
            call adjust_load_balance(iter_data_fciqmc)
        end if

        call MPIBarrier(ierr, tTimeIn=.false.)

        call set_timer(SemiStoch_Init_Time)

        if (t_global_core_space) then
            num_core_runs = 1
            num_core_runs = inum_runs
        end if
        ! Allocate the corespace replicas

        write(stdout, '(/,12("="),1x,a30,1x,12("="))') "Semi-stochastic initialisation"
        call neci_flush(stdout)

        do run = 1, size(cs_replicas)
            associate(rep => cs_replicas(run))
                allocate(rep%determ_sizes(0:nProcessors - 1))
                allocate(rep%determ_displs(0:nProcessors - 1))
                allocate(rep%determ_last(0:nProcessors - 1))
                call rep%associate_run(t_global_core_space, run)
                rep%determ_sizes = 0_MPIArg
                rep%determ_last = 0_MPIArg

                if (.not. (tStartCAS .or. core_in%tPops .or. core_in%tDoubles .or. core_in%tCAS .or. core_in%tRAS .or. &
                           core_in%tOptimised .or. core_in%tLowE .or. core_in%tRead .or. core_in%tMP1 .or. &
                           core_in%tFCI .or. core_in%tHeisenbergFCI .or. core_in%tHF .or. &
                           core_in%tPopsAuto .or. core_in%tPopsProportion)) then
                    call stop_all("init_semi_stochastic", "You have not selected a semi-stochastic core space to use.")
                end if
                if (.not. tUseRealCoeffs) call stop_all(t_r, "To use semi-stochastic you must also use real coefficients.")

                ! Call the enumerating subroutines to create all excitations and add these states to
                ! SpawnedParts on the correct processor. As they do this, they count the size of the
                ! deterministic space (on their own processor only).
                write(stdout, '("Generating the deterministic space...")'); call neci_flush(stdout)
                if (core_in%tPopsAuto) then
                    write(stdout, '("Choosing 10% of initiator space as core space, if this number is larger than 50000, then use 50000")')
                    call neci_flush(stdout)
                    ! from my understanding npops refers to the total core space size
                    write(stdout, '("Estimated size of core space:",1X,i5)') int(AllNoInitDets(run) * 0.1_dp)
                    call neci_flush(stdout)
                    if (int(AllNoInitDets(run) * 0.1_dp) > 50000) then
                        core_in%npops = 50000
                        core_in%npops = int(AllNoInitDets(run) * 0.1_dp)
                    end if
                    ! might also need to check if the total space is too large so that
                    ! tApproxSpace should be used instead...
                end if

                if (core_in%tPopsProportion) then
                    write(stdout, '("Choosing ",F7.2,"% of initiator space as core space")') 100 * core_in%npops_proportion
                    write(stdout, *) "Estimated size of core space: ", int(AllNoInitDets(run) * core_in%npops_proportion)
                    core_in%npops = max(1, int(AllNoInitDets(run) * core_in%npops_proportion))
                end if

                if (core_in%tApproxSpace) write(stdout, '(" ... approximately using the factor of",1X,i5)') core_in%nApproxSpace
                call generate_space(core_in, run)

                ! So that all procs store the size of the deterministic spaces on all procs.
                mpi_temp = rep%determ_sizes(iProcIndex)
                call MPIAllGather(mpi_temp, rep%determ_sizes, ierr)

                rep%determ_space_size = sum(rep%determ_sizes)
                rep%determ_space_size_int = int(rep%determ_space_size)

                ! This is now the total size on the replica with the largest space
                ! Typically, all replicas will have either similar or the same space size
                write(stdout, '("Total size of deterministic space:",1X,i8)') rep%determ_space_size
                if (rep%determ_space_size > (AllNoInitDets(run) .div. 2_int64) .and. tTruncInitiator) then
                    write(stdout, *)"WARNING: Total size of deterministic space is greater than&
                        & 50% of the initiator space."
                    write(stdout, *)"         Reducing the size of the deterministic space is&
                        & encouraged."
                    if (iProcIndex == 0) then
                        write(stderr, *)"WARNING: Total size of deterministic space is greater than&
                            & 50% of the initiator space."
                        write(stderr, *)"         Reducing the size of the deterministic space is&
                            & encouraged."
                    end if
                end if
                write(stdout, '("Size of deterministic space on this processor:",1X,i8)') rep%determ_sizes(iProcIndex)
                call neci_flush(stdout)

                call rep%alloc_wf()

                ! This array will hold the positions of the deterministic states in CurrentDets.
                allocate(rep%indices_of_determ_states(rep%determ_sizes(iProcIndex)), stat=ierr)
                call LogMemAlloc('indices_of_determ_states', int(rep%determ_sizes(iProcIndex)), &
                                  sizeof_int, t_r, rep%IDetermTag, ierr)

                ! Calculate the indices in the full vector at which the various processors end.
                rep%determ_last(0) = rep%determ_sizes(0)
                do i = 1, nProcessors - 1
                    rep%determ_last(i) = rep%determ_last(i - 1) + rep%determ_sizes(i)
                end do

                call sort(spawnedparts(0:NIfTot, 1:rep%determ_sizes(iprocindex)), ilut_lt, ilut_gt)

                ! Do a check that no states are in the deterministic space twice. The list is sorted
                ! already so simply check states next to each other in the list.
                do i = 2, rep%determ_sizes(iProcIndex)
                    if (all(SpawnedParts(0:nifd, i - 1) == SpawnedParts(0:nifd, i))) then
                        call decode_bit_det(nI, SpawnedParts(:, i))
                        write(stdout, '("State found twice:")')
                        write(stdout, *) SpawnedParts(:, i)
                        call write_det(stdout, nI, .true.)
                        call stop_all("init_semi_stochastic", &
                                      "The same state has been found twice in the deterministic space.")
                    end if
                end do

                ! Store every core determinant from all processors on all processors, in core_space.
                call store_whole_core_space(rep)
                ! Create the hash table to address the core determinants.
                call initialise_shared_rht(rep%core_space, rep%determ_space_size_int, rep%core_ht)

                if (tWriteCore) call write_core_space(rep)

                write(stdout, '("Generating the Hamiltonian in the deterministic space...")'); call neci_flush(stdout)
                if (tAllSymSectors .or. tReltvy .or. nOccAlpha <= 1 .or. nOccBeta <= 1 &
                    .or. tGUGA .or. t_mol_3_body .or. t_non_hermitian_2_body) then
                    ! In the above cases the faster generation is not implemented, so
                    ! use the original algorithm.
                    call set_timer(SemiStoch_Hamil_Time)
                    if (tHPHF) then
                        call calc_determ_hamil_sparse_hphf(rep)
                        call calc_determ_hamil_sparse(rep)
                    end if
                    call halt_timer(SemiStoch_Hamil_Time)
                    write(stdout, '("Total time (seconds) taken for Hamiltonian generation:", f9.3)') &
                    if (tHPHF) then
                        call calc_determ_hamil_opt_hphf(rep)
                        call calc_determ_hamil_opt(rep)
                    end if
                end if

#ifndef CMPLX_
                if (t_print_core_info) then
                    ! i think i also want information, like the energy and the
                    ! eigenvectors of the core-space
                    if (t_print_core_hamil) then
                        call print_basis(rep)
                        call print_hamiltonian(rep)
                    end if
                    root_print "I am before the diagonalization step with", t_non_hermitian_2_body
                    if (t_non_hermitian_2_body) then
                        call diagonalize_core_non_hermitian(e_values, e_vectors, rep)
                        if (t_choose_trial_state) then
                            gs_energy = e_values(trial_excit_choice(1))
                            gs_energy = e_values(1)
                        end if
                        call diagonalize_core(gs_energy, gs_vector, rep)
                    end if
                    root_print "semi-stochastic space GS energy: ", gs_energy

                    if (t_print_core_vec) then
                        root_print "semi-stochastic gs vector:"
                        call print_vec(gs_vector, "semistoch-vec")
                    end if
                end if

                if (tRDMonFly) call generate_core_connections(rep)

                if (tDetermProjApproxHamil) call init_var_space(rep)

                ! Move the states to CurrentDets.
                call add_core_states_currentdet_hash(run)

                ! If using a trial wavefunction, and that initialisation has already
                ! been performed, then the current_trial_amps array needs correcting
                ! after the core states were added and sorted into CurrentDets.
                call reinit_current_trial_amps()

                ! If starting from a popsfile then global_determinant_data will not
                ! have been initialised, or if in the middle of a calculation then new
                ! determinants may have been added.
                if (tReadPops .or. (semistoch_shift_iter /= 0)) call fill_in_diag_helements()

                SpawnedParts = 0_n_int
                TotWalkersOld = TotWalkers

                tStartedFromCoreGround = .false.
                if (tStartCoreGroundState .and. (.not. tReadPops) .and. tStaticCore .and. (.not. tTrialInit)) then
                    if (t_non_hermitian_2_body) then
                        call set_timer(SemiStoch_nonhermit_Time)
                        call start_walkers_from_core_ground_nonhermit(tPrintInfo=.true., run=run)
                        call halt_timer(SemiStoch_nonhermit_Time)
                        tStartedFromCoreGround = .true.
                        write(stdout, '("Total time (seconds) taken for non-hermitian diagonalization:", f9.3)') &
                        call set_timer(SemiStoch_Davidson_Time)
                        call start_walkers_from_core_ground(tPrintInfo=.true., run=run)
                        call halt_timer(SemiStoch_Davidson_Time)
                        tStartedFromCoreGround = .true.
                        write(stdout, '("Total time (seconds) taken for Davidson calculation:", f9.3)') &
                    end if
                end if
            end associate
        end do
        ! Call MPIBarrier here so that Semistoch_Init_Time will give the
        ! initialisation time for all processors to finish.
        call MPIBarrier(ierr, tTimeIn=.false.)

        call halt_timer(SemiStoch_Init_Time)

        write(stdout, '("Semi-stochastic initialisation complete.")')
        write(stdout, '("Time (seconds) taken for semi-stochastic initialisation:", f9.3, /)') &
        call neci_flush(stdout)

    end subroutine init_semi_stochastic