calc_determ_hamil_sparse_hphf Subroutine

public subroutine calc_determ_hamil_sparse_hphf(rep)

Arguments

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

Contents


Source Code

    subroutine calc_determ_hamil_sparse_hphf(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, allocatable :: temp_store_nI(:, :)
        integer(TagIntType) :: HRTag, TempStoreTag
        HElement_t(dp), allocatable, dimension(:) :: hamiltonian_row
        HElement_t(dp) :: HOffDiag
        character(len=*), parameter :: t_r = "calc_determ_hamil_sparse_hphf"

        integer(n_int) :: tmp(0:NIfD)
        integer :: IC
        logical :: CS_I
        logical, allocatable :: cs(:)

        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, t_r, 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, t_r, TempStoreTag, ierr)
        allocate(temp_store_nI(nel, rep%determ_space_size), stat=ierr)
        allocate(cs(rep%determ_space_size), stat=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, rep%determ_sizes, rep%determ_displs)

        do i = 1, rep%determ_space_size
            call decode_bit_det(temp_store_nI(:, i), temp_store(:, i))
            cs(i) = TestClosedShellDet(temp_store(:, i))
        end do

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

            !call decode_bit_det(nI, SpawnedParts(:, i))
            nI = temp_store_nI(:, i + rep%determ_displs(iProcIndex))

            row_size = 0
            hamiltonian_row = 0.0_dp

            CS_I = cs(i + rep%determ_displs(iProcIndex))

            ! Loop over all deterministic states.
            do j = 1, rep%determ_space_size

                !call decode_bit_det(nJ, temp_store(:,j))
                nJ = temp_store_nI(:, j)

                ! If on the diagonal of the Hamiltonian.
                if (j == i + rep%determ_displs(iProcIndex)) then
                    hamiltonian_row(j) = get_diagonal_matel(nI, SpawnedParts(:, i)) - 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, SpawnedParts(:, i))
                        call set_det_offdiagH(i, HOffDiag)
                    end if
                    ! Always include the diagonal elements.
                    row_size = row_size + 1
                else
                    tmp = ieor(SpawnedParts(0:NIfD, i), temp_store(0:NIfD, j))
                    tmp = iand(SpawnedParts(0:NIfD, i), tmp)
                    IC = CountBits(tmp, NIfD)

                    if (IC <= maxExcit .or. ((.not. CS_I) .and. (.not. cs(j)))) then
                        if(t_evolve_adjoint(rep%first_run())) then
                            hamiltonian_row(j) = hphf_off_diag_helement_opt(nJ, temp_store(:, j), SpawnedParts(:, i), IC, cs(j), CS_I)
                        else
                            hamiltonian_row(j) = hphf_off_diag_helement_opt(nI, SpawnedParts(:, i), temp_store(:, j), IC, CS_I, cs(j))
                        endif
                        if (abs(hamiltonian_row(j)) > 0.0_dp) row_size = row_size + 1
                    end if
                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(t_r, TempStoreTag, ierr)
        deallocate(hamiltonian_row, stat=ierr)
        call LogMemDealloc(t_r, HRTag, ierr)
        deallocate(temp_store_nI, stat=ierr)
        deallocate(cs, stat=ierr)

    end subroutine calc_determ_hamil_sparse_hphf