generate_optimised_space Subroutine

public subroutine generate_optimised_space(opt_data, tLimitSpace, ilut_list, space_size, max_space_size)


Type IntentOptional Attributes Name
type(opt_space_data), intent(in) :: opt_data
logical :: tLimitSpace
integer(kind=n_int), intent(inout) :: ilut_list(0:,:)
integer, intent(inout) :: space_size
integer, intent(in), optional :: max_space_size


Source Code

    subroutine generate_optimised_space(opt_data, tLimitSpace, ilut_list, space_size, max_space_size)

        ! This routine generates a deterministic space by diagonalising a small fraction
        ! of the whole space, and choosing the basis states with the largest weights in
        ! the ground state for the deterministic states. This therefore aims to find
        ! some of the basis states with the largest weights in the true ground state.

        ! More specifically, we start with the Hartree-Fock state, and generate a first
        ! space by finding all states connected to it. We then find the ground state of
        ! the Hamiltonian in this space. Using this ground state, we pick out the basis
        ! states with the largest amplitudes (up to some user-specified criteria), and
        ! define these basis states as our new space. We then find all states connected
        ! to the states in this space, and diagonalise the Hamiltonian in this new space
        ! and again pick out the basis states with the largest weights. This process is
        ! iterated for as many time as the user requests.

        ! In: opt_data: Derived type containing the parameters for the algorithm to
        !     generate the space.
        ! In: tLimitSpace: If true then, every iteration of generating algorithm, after
        !     all connections have generated, the space is cut off to have a space with
        !     a maximum size of max_space_size. This is done by removing the highest
        !     energy determinants.
        ! In/Out: ilut_list - List of determinants generated.
        ! In/Out: space_size - Number of determinants in the generated space.
        !             If ilut_list is not empty on input and you want to keep
        !             the states already in it, then on input space_size should
        !             be equal to the number of states to be kept in ilut_list,
        !             and new states will be added in from space_size+1.
        !             Otherwise, space_size must equal 0 on input.
        !             On output space_size will equal the total number of
        !             generated plus what space_size was on input.
        ! In (optional): max_space_size - Only used if tLimitSpace is true. See
        !     tLimitSpace for an explanation of use.

        type(opt_space_data), intent(in) :: opt_data
        logical :: tLimitSpace
        integer(n_int), intent(inout) :: ilut_list(0:, :)
        integer, intent(inout) :: space_size
        integer, optional, intent(in) :: max_space_size

        integer(n_int), allocatable, dimension(:, :) :: ilut_store, temp_space
        integer :: counter, i, j, ierr
        integer :: old_num_states, new_num_states
        integer(MPIArg) :: proc_space_sizes(0:nProcessors - 1), disps(0:nProcessors - 1), &
                           sendcounts(0:nProcessors - 1), recvcount, this_proc_size
        integer(TagIntType) :: IlutTag, TempTag
        character(len=*), parameter :: t_r = "generate_optimised_space"

        type(DavidsonCalcType) :: davidsonCalc

        if (iProcIndex /= root) then
            ! Allocate some space so that the MPIScatterV call does not crash.
            allocate(ilut_store(0:NIfTot, 1), stat=ierr)
            call LogMemAlloc("ilut_store", (NIfTot + 1), size_n_int, t_r, IlutTag, ierr)
        else if (iProcIndex == root) then
            ! Allocate the stores of ilut's that will hold these deterministic states.
            ! For now, assume that we won't have a deterministic space of more than one
            ! million states. Could make this user-specified later.
            allocate(ilut_store(0:NIfTot, 1000000), stat=ierr)
            call LogMemAlloc("ilut_store", 1000000 * (NIfTot + 1), size_n_int, t_r, IlutTag, ierr)
            allocate(temp_space(0:NIfTot, 1000000), stat=ierr)
            call LogMemAlloc("temp_store", 1000000 * (NIfTot + 1), size_n_int, t_r, TempTag, ierr)
            ilut_store = 0_n_int
            temp_space = 0_n_int

            ! Put the Hartree-Fock state in the list first.
            ilut_store(0:NIfTot, 1) = ilutHF(0:NIfTot)

            ! old_num_states will hold the number of deterministic states in the current
            ! space. This is just 1 for now, with only the Hartree-Fock.
            old_num_states = 1

            ! Now we start the iterating loop. Find all states which are either a single or
            ! double excitation from each state in the old ilut store, and then see if they
            ! have a non-zero Hamiltonian matrix element with any state in the old ilut store:

            ! Over the total number of iterations.
            do i = 1, opt_data%ngen_loops

                write(stdout, '(a37,1X,i2)') "Optimised space generation: Iteration", i
                call neci_flush(stdout)

                ! Find all states connected to the states currently in ilut_store.
                write(stdout, '(a29)') "Generating connected space..."
                call neci_flush(stdout)
                ! Allow for up to 1 million connected states to be created.
                new_num_states = 1000000
                call generate_connected_space(old_num_states, ilut_store(:, 1:old_num_states), &
                                              new_num_states, temp_space(:, 1:1000000))
                write(stdout, '(a26)') "Connected space generated."
                call neci_flush(stdout)

                ! Add these states to the ones already in the ilut stores.
                ilut_store(:, old_num_states + 1:old_num_states + new_num_states) = &
                    temp_space(:, 1:new_num_states)

                new_num_states = new_num_states + old_num_states

                call remove_repeated_states(ilut_store, new_num_states)

                write(stdout, '(i8,1X,a13)') new_num_states, "states found."
                call neci_flush(stdout)

                if (tLimitSpace) call remove_high_energy_orbs(ilut_store(:, 1:new_num_states), &
                                                              new_num_states, max_space_size, .false.)

                write(stdout, '(a27)') "Constructing Hamiltonian..."
                call neci_flush(stdout)

                if (t_non_hermitian_2_body) then
                    call calculate_sparse_hamiltonian_non_hermitian(new_num_states, ilut_store(:, 1:new_num_states))
                    call calculate_sparse_hamiltonian(new_num_states, ilut_store(:, 1:new_num_states))
                end if

                write(stdout, '(a29)') "Performing diagonalisation..."
                call neci_flush(stdout)

                ! Now that the Hamiltonian is generated, we can finally find the ground state of it:
                if (t_non_hermitian_2_body) then
                    call stop_all(t_r, &
                                  "perform_davidson not adapted for non-hermitian Hamiltonians!")
                end if
                call perform_davidson(davidsonCalc, sparse_hamil_type, .false.)

                associate ( &
                    davidson_eigenvector => davidsonCalc%davidson_eigenvector &

                    ! davidson_eigenvector now stores the ground state eigenvector. We want to use the
                    ! vector whose components are the absolute values of this state:
                    davidson_eigenvector = abs(davidson_eigenvector)
                    ! Multiply by -1.0_dp so that the sort operation later is the right way around.
                    davidson_eigenvector = -1.0_dp * davidson_eigenvector

                    ! Now decide which states to keep for the next iteration. There are two ways of
                    ! doing this, as decided by the user. Either all basis states with an amplitude
                    ! in the ground state greater than a given value are kept (tAmpCutoff = .true.),
                    ! or a number of states to keep is specified and we pick the states with the
                    ! biggest amplitudes (tAmpCutoff = .false.).
                    if (opt_data%tAmpCutoff) then
                        counter = 0
                        do j = 1, new_num_states
                            if (abs(davidson_eigenvector(j)) > opt_data%cutoff_amps(i)) then
                                counter = counter + 1
                                ilut_store(:, counter) = ilut_store(:, j)
                            end if
                        end do
                        old_num_states = counter
                        ! Sort the list in order of the amplitude of the states in the ground state,
                        ! from largest to smallest.
                        call sort(davidson_eigenvector(:), ilut_store(:, 1:new_num_states))

                        ! Set old_num_states to specify the number of states which should be kept.
                        old_num_states = opt_data%cutoff_nums(i)
                        if (old_num_states > new_num_states) old_num_states = new_num_states
                    end if

                    write(stdout, '(i8,1X,a12)') old_num_states, "states kept."
                    call neci_flush(stdout)

                    call deallocate_sparse_ham(sparse_ham, SparseHamilTags)
                    deallocate(hamil_diag, stat=ierr)
                    call LogMemDealloc(t_r, HDiagTag, ierr)

                end associate
            end do

        end if ! If on root.

        ! At this point, the space has been generated on the root processor. The rest is
        ! just sending the right info to the right processors...

        if (iProcIndex == root) then
            ! Find which core each state belongs to and sort accordingly.
            call sort_space_by_proc(ilut_store(:, 1:old_num_states), old_num_states, proc_space_sizes)

            ! Create displacement and sendcount arrays for MPIScatterV later:
            sendcounts = int(proc_space_sizes * (NIfTot + 1), MPIArg)
            disps(0) = 0_MPIArg
            do i = 1, nProcessors - 1
                disps(i) = int(disps(i - 1) + proc_space_sizes(i - 1) * (NIfTot + 1), MPIArg)
            end do
        end if

        ! Send the number of states on each processor to the corresponding processor.
        call MPIScatter(proc_space_sizes, this_proc_size, ierr)
        recvcount = int(this_proc_size * (NIfTot + 1), MPIArg)

        ! Finally send the actual determinants to the ilut_list array.
        call MPIScatterV(ilut_store, sendcounts, disps, ilut_list(:, space_size + 1:space_size + this_proc_size), recvcount, ierr)

        space_size = space_size + int(this_proc_size)

        ! Finally, deallocate arrays.
        if (allocated(ilut_store)) then
            deallocate(ilut_store, stat=ierr)
            call LogMemDealloc(t_r, IlutTag, ierr)
        end if
        if (iProcIndex == root) then
            deallocate(temp_space, stat=ierr)
            call LogMemDealloc(t_r, TempTag, ierr)
        end if

    end subroutine generate_optimised_space