read_dense Subroutine

private subroutine read_dense(this, filename)

Read the 6-index integrals from a file to dense format @param[in] filename name of the file to read from

Type Bound



Type IntentOptional Attributes Name
class(dense_lMat_t), intent(inout) :: this
character(len=*), intent(in) :: filename


Source Code

Source Code

    subroutine read_dense(this, filename)
        class(dense_lMat_t), intent(inout) :: this
        character(*), intent(in) :: filename
        integer :: iunit, ierr
        integer(int64) :: indices(6)
        HElement_t(dp) :: matel
        character(*), parameter :: t_r = "readLMat"
        integer(int64) :: counter, index
        character(1024) :: char1024

        call this%alloc(this%lMat_size())

        if (tHDF5LMat) then
#ifdef USE_HDF_
            call this%read_hdf5_dense(filename)
            call stop_all(t_r, "HDF5 integral files disabled at compile time")
            if (iProcIndex_intra == 0) then
                iunit = get_free_unit()
                open(iunit, file=filename, status='old')
                ! If first line cannot be parsed as matrix element we skip it
                ! (it would contain the number of orbitals then).
                read(iunit, '(a)', iostat=ierr) char1024
                if (ierr/=0) call stop_all(t_r, "Error reading TCDUMP file")
                ! Only assume first two indices are on same line.
                read(char1024, *, iostat=ierr) matel, indices(1:2)
                if (ierr==0) rewind(iunit)
                counter = 0
                    read(iunit,*, iostat = ierr) matel, indices
                    ! end of file reached?
                    if (ierr < 0) then
                    else if (ierr > 0) then
                        ! error while reading?
                        call stop_all(t_r, "Error reading TCDUMP file")
                        ! permutational factor of 3 (not accounted for in integral calculation)
                        matel = 3.0_dp * matel
                        call freeze_lmat(matel,indices)
                        ! else assign the matrix element

                        if(abs(matel) > LMatEps) then
                            counter = counter + 1
                            index = this%indexFunc(&
                                indices(1),indices(2),indices(3),indices(4),indices(5), indices(6))
                            if(index > this%lMat_size()) then
                                counter = index
                                write(stdout, *) "Warning, exceeding size"
                            this%lMat_vals%ptr(index) = matel

                end do

                write(stdout, *) "Sparsity of LMat", real(counter) / real(this%lMat_size())
                write(stdout, *) "Nonzero elements in LMat", counter
                write(stdout, *) "Allocated size of LMat", this%lMat_size()
            end if
            call MPIBcast(counter)
        end if

        ! Synchronize the shared resource
        call this%lMat_vals%sync()
    end subroutine read_dense