calc_determ_hamil_sparse Subroutine

public subroutine calc_determ_hamil_sparse(rep)

Arguments

Type IntentOptional Attributes Name
type(core_space_t), intent(inout) :: rep

Contents


Source Code

    subroutine calc_determ_hamil_sparse(rep)

        type(core_space_t), intent(inout) :: rep
        integer :: i, j, row_size, counter, ierr
        integer :: nI(nel), nJ(nel)
        integer(n_int), allocatable, dimension(:, :) :: temp_store
        integer(TagIntType) :: HRTag, TempStoreTag
        HElement_t(dp), allocatable, dimension(:) :: hamiltonian_row

        integer :: pos, nExcit
        integer(n_int) :: ilutG(0:nifguga)
        integer(n_int), pointer :: excitations(:, :)
        type(ExcitationInformation_t) :: excitInfo

        character(len=*), parameter :: this_routine = "calc_determ_hamil_sparse"

        integer(n_int) :: tmp(0:NIfD), ilutI_tmp(0:NIfTot), ilutJ_tmp(0:niftot)
        integer :: IC, nI_tmp(nel), nJ_tmp(nel)
        integer(n_int) :: ilutI(0:niftot), ilutJ(0:niftot)
        type(CSF_Info_t) :: csf_i, csf_j
        HElement_t(dp) :: tmp_mat, tmp_mat_2, HOffDiag

        allocate(rep%sparse_core_ham(rep%determ_sizes(iProcIndex)), stat=ierr)
        allocate(rep%SparseCoreHamilTags(2, rep%determ_sizes(iProcIndex)))
        allocate(hamiltonian_row(rep%determ_space_size), stat=ierr)
        call LogMemAlloc('hamiltonian_row', int(rep%determ_space_size), 8, this_routine, HRTag, ierr)
        allocate(rep%core_ham_diag(rep%determ_sizes(iProcIndex)), stat=ierr)
        allocate(temp_store(0:NIfTot, rep%determ_space_size), stat=ierr)
        call LogMemAlloc('temp_store', rep%determ_space_size * (NIfTot + 1), 8, this_routine, TempStoreTag, ierr)

        ! Stick together the deterministic states from all processors, on
        ! all processors.
        ! n.b. Explicitly use 0:NIfTot, as NIfTot may not equal NIfBCast
        call MPIAllGatherV(SpawnedParts(0:NIfTot, 1:rep%determ_sizes(iProcIndex)), &
                           temp_store(0:niftot, 1:), rep%determ_sizes, rep%determ_displs)


        ! Loop over all deterministic states on this processor.
        do i = 1, rep%determ_sizes(iProcIndex)

            ilutI = SpawnedParts(0:niftot, i)
            call decode_bit_det(nI, IlutI)
            if (tGUGA) csf_i = CSF_Info_t(ilutI)

            row_size = 0
            hamiltonian_row = 0.0_dp
            ! Loop over all deterministic states.
            do j = 1, rep%determ_space_size

                ilutJ = temp_store(:, j)
                call decode_bit_det(nJ, ilutJ)
                if (tGUGA) csf_j = CSF_Info_t(ilutJ)

                if(t_evolve_adjoint(rep%first_run())) then
                    nI_tmp = nJ
                    nJ_tmp = nI
                    ilutI_tmp = ilutJ
                    ilutJ_tmp = ilutI
                else
                    nI_tmp = nI
                    nJ_tmp = nJ
                    ilutI_tmp = ilutI
                    ilutJ_tmp = ilutJ
                end if

                ! If on the diagonal of the Hamiltonian.
                if (DetBitEq(IlutI, ilutJ_tmp, nifd)) then
                    hamiltonian_row(j) = get_diagonal_matel(nI_tmp, IlutI_tmp) - Hii
                    rep%core_ham_diag(i) = hamiltonian_row(j)
                    ! We calculate and store the diagonal matrix element at
                    ! this point for later access.
                    if (.not. tReadPops) then
                        call set_det_diagH(i, Real(hamiltonian_row(j), dp))
                        HOffDiag = get_off_diagonal_matel(nI_tmp, ilutI_tmp)
                        call set_det_offdiagH(i, HOffDiag)
                    end if
                    ! Always include the diagonal elements.
                    row_size = row_size + 1
                else
                    if (tHPHF) then
                        hamiltonian_row(j) = hphf_off_diag_helement(nI_tmp, nJ_tmp, IlutI_tmp, ilutJ_tmp)
                    else if (tGUGA) then
                        call calc_guga_matrix_element(&
                                ilutI_tmp, csf_i, ilutJ_tmp, csf_j, excitInfo, tmp_mat, .true.)
#ifdef DEBUG_
                        call calc_guga_matrix_element(&
                                ilutI_tmp,  csf_i, ilutJ_tmp, csf_j, excitInfo, tmp_mat_2, .true.)
                        if (.not. near_zero(tmp_mat - tmp_mat_2)) then
                            call stop_all(this_routine, "type 1 and 2 do not agree!")
                        end if
                        call calc_guga_matrix_element(&
                                ilutJ_tmp, csf_j, ilutI_tmp, csf_i, excitInfo, tmp_mat_2, .true.)
                        if (.not. near_zero(tmp_mat - conjgt(tmp_mat_2))) then
                            call stop_all(this_routine, "not hermititan!")
                        end if
#endif

#ifdef CMPLX_
                        hamiltonian_row(j) = conjg(tmp_mat)
#else
                        hamiltonian_row(j) = tmp_mat
#endif
                    else

                        tmp = ieor(IlutI_tmp(0:NIfD), ilutJ_tmp(0:NIfD))
                        tmp = iand(IlutI_tmp(0:NIfD), tmp)
                        IC = CountBits(tmp, NIfD)

                        if (IC <= maxExcit) then
!                             hamiltonian_row(j) = get_helement(nI_tmp, nJ_tmp, IC, ilutI_tmp, ilutJ_tmp)
                            hamiltonian_row(j) = get_helement(nJ_tmp, nI_tmp, IC, ilutJ_tmp, ilutI_tmp)
                        end if
                    end if
                    if (abs(hamiltonian_row(j)) > 0.0_dp) row_size = row_size + 1
                end if
            end do
            ! Now we know the number of non-zero elements in this row of the Hamiltonian, so allocate it.
            call allocate_sparse_ham_row(rep%sparse_core_ham, i, row_size, "sparse_core_ham", rep%SparseCoreHamilTags(:, i))

            rep%sparse_core_ham(i)%elements = 0.0_dp
            rep%sparse_core_ham(i)%positions = 0
            rep%sparse_core_ham(i)%num_elements = row_size

            counter = 1
            do j = 1, rep%determ_space_size
                ! If non-zero or a diagonal element.
                if (abs(hamiltonian_row(j)) > 0.0_dp .or. (j == i + rep%determ_displs(iProcIndex))) then
                    rep%sparse_core_ham(i)%positions(counter) = j
                    rep%sparse_core_ham(i)%elements(counter) = hamiltonian_row(j)
                    counter = counter + 1
                end if
                if (counter == row_size + 1) exit
            end do

        end do

        ! Don't time the mpi_barrier call, because if we did then we wouldn't
        ! be able separate out some of the core Hamiltonian creation time from
        ! the MPIBarrier calls in the main loop.
        call MPIBarrier(ierr, tTimeIn=.false.)

        deallocate(temp_store, stat=ierr)
        call LogMemDealloc(this_routine, TempStoreTag, ierr)
        deallocate(hamiltonian_row, stat=ierr)
        call LogMemDealloc(this_routine, HRTag, ierr)

    end subroutine calc_determ_hamil_sparse