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))
else
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
else
! 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