init_spectral_lanczos Subroutine

    subroutine init_spectral_lanczos()

        use bit_rep_data, only: nifd, extract_sign
        use DetBitOps, only: DetBitEq, EncodeBitDet, ilut_lt, ilut_gt
        use gndts_mod, only: gndts
        use load_balance_calcnodes, only: DetermineDetNode
        use sort_mod, only: sort
        use SystemData, only: nbasis, nel, BRR, nBasisMax, G1, tSpn, LMS, tParity, SymRestrict

        integer :: ndets_this_proc, ndets_tot
        integer(MPIArg) :: mpi_temp
        integer, allocatable :: nI_list(:, :)

        integer(n_int) :: ilut(0:NIfTot)
        character(len=*), parameter :: t_r = 'init_spectral_lanczos'
        integer :: i, j, ndets, proc, ierr, temp(1, 1), hf_ind
        real(dp) :: real_sign(lenof_sign)
        real(dp) :: norm_pert

        write(stdout, '(/,1x,a39,/)') "Beginning spectral Lanczos calculation."
        call neci_flush(stdout)

        allocate(ndets_sl(0:nProcessors - 1))
        allocate(disps_sl(0:nProcessors - 1))

        write(stdout, '(1x,a56)', advance='yes') "Enumerating and storing all determinants in the space..."
        call neci_flush(stdout)

        ! Determine the total number of determinants.
        call gndts(nel, nbasis, BRR, nBasisMax, temp, .true., G1, tSpn, lms, tParity, SymRestrict, ndets, hf_ind)
        allocate(nI_list(nel, ndets))
        ! Now generate them again and store them this time.
        call gndts(nel, nbasis, BRR, nBasisMax, nI_list, .false., G1, tSpn, lms, tParity, SymRestrict, ndets, hf_ind)

        ! Count the number of determinants belonging to this processor.
        ndets_this_proc = 0
        do i = 1, ndets
            proc = DetermineDetNode(nel, nI_list(:, i), 0)
            if (proc == iProcIndex) ndets_this_proc = ndets_this_proc + 1
        end do

        allocate(sl_ilut_list(0:NIfTot, ndets_this_proc))
        sl_ilut_list = 0_n_int

        j = 0
        do i = 1, ndets
            proc = DetermineDetNode(nel, nI_list(:, i), 0)
            if (proc == iProcIndex) then
                call EncodeBitDet(nI_list(:, i), ilut)
                j = j + 1
                sl_ilut_list(0:nifd, j) = ilut(0:nifd)
            end if
        end do

        call sort(sl_ilut_list, ilut_lt, ilut_gt)

        write(stdout, '(1x,a9)') "Complete."
        call neci_flush(stdout)

        mpi_temp = int(ndets_this_proc, MPIArg)
        call MPIAllGather(mpi_temp, ndets_sl, ierr)

        disps_sl(0) = 0
        do i = 1, nProcessors - 1
            disps_sl(i) = disps_sl(i - 1) + ndets_sl(i - 1)
        end do

        ndets_tot = int(sum(ndets_sl))

        write(stdout, '(1x,a44)', advance='no') "Allocating arrays to hold Lanczos vectors..."
        call neci_flush(stdout)
        allocate(sl_vecs(ndets_this_proc, n_lanc_vecs_sl))
        sl_vecs = 0.0_dp
        allocate(full_vec_sl(ndets_tot))
        write(stdout, '(1x,a9)') "Complete."
        call neci_flush(stdout)

        allocate(pert_ground_left(ndets_this_proc), stat=ierr)
        if (ierr /= 0) then
            write(stdout, '(1x,a11,1x,i5)') "Error code:", ierr
            call stop_all(t_r, "Error allocating array to hold left perturbed ground state components.")
        end if

        allocate(sl_hamil(n_lanc_vecs_sl, n_lanc_vecs_sl))
        allocate(sl_h_eigv(n_lanc_vecs_sl))
        allocate(sl_overlaps(n_lanc_vecs_sl))
        allocate(trans_amps_left(n_lanc_vecs_sl))
        allocate(trans_amps_right(n_lanc_vecs_sl))
        sl_hamil = 0.0_dp
        sl_h_eigv = 0.0_dp
        sl_overlaps = 0.0_dp
        trans_amps_left = 0.0_dp
        trans_amps_right = 0.0_dp

        write(stdout, '(1x,a48)') "Allocating and calculating Hamiltonian matrix..."
        call neci_flush(stdout)
        call calculate_sparse_ham_par(ndets_sl, sl_ilut_list, .true.)
        write(stdout, '(1x,a48,/)') "Hamiltonian allocation and calculation complete."
        call neci_flush(stdout)

        call return_perturbed_ground_spec(left_perturb_spectral, sl_ilut_list, pert_ground_left, left_pert_norm)
        call return_perturbed_ground_spec(right_perturb_spectral, sl_ilut_list, sl_vecs(:, 1), right_pert_norm)

        ! Normalise the perturbed wave function, currently stored in sl_vecs(:,1),
        ! so that it can be used as the first lanczos vector.
        sl_vecs(:, 1) = sl_vecs(:, 1) * sqrt(pops_norm) / right_pert_norm

    end subroutine init_spectral_lanczos