read_spinfree_2rdm_files Subroutine

public subroutine read_spinfree_2rdm_files(rdm_defs, rdm, spawn)

Arguments

Type IntentOptional Attributes Name
type(rdm_definitions_t), intent(in) :: rdm_defs
type(rdm_list_t), intent(inout) :: rdm
type(rdm_spawn_t), intent(inout) :: spawn

Contents


Source Code

    subroutine read_spinfree_2rdm_files(rdm_defs, rdm, spawn)

        ! Read spinfree 2-RDMs from standard NECI spinfree_TwoRDM files.

        ! The resulting object will be in the form of a spinfree 2-RDM, *not*
        ! of a standard NECI RDM, so care should be taken - one can not expect
        ! most routines to act correctly on this resulting object. For example,
        ! only spatial orbital labels or stored, symmetry is not taken account
        ! of (so i<j and k<l is not enforced) and orbital labels are in a
        ! different order to that of traditional NECI 2-RDMs.

        use hash, only: clear_hash_table, FindWalkerHash, add_hash_table_entry
        use Parallel_neci, only: iProcIndex, nProcessors
        use rdm_data, only: rdm_list_t, rdm_spawn_t
        use rdm_data, only: rdm_definitions_t
        use rdm_data_utils, only: calc_separate_rdm_labels, extract_sign_rdm, add_to_rdm_spawn_t
        use rdm_data_utils, only: communicate_rdm_spawn_t_wrapper, annihilate_rdm_list
        use util_mod, only: get_free_unit, int_fmt

        type(rdm_definitions_t), intent(in) :: rdm_defs
        type(rdm_list_t), intent(inout) :: rdm
        type(rdm_spawn_t), intent(inout) :: spawn

        integer(int_rdm) :: ijkl, rdm_entry(0:rdm%sign_length)
        integer :: ij, kl, i, j, k, l, ij_proc_row, sign_length_old
        integer :: irdm, ielem, iunit(rdm%sign_length), ierr, hash_val
        real(dp) :: rdm_sign(rdm%sign_length)
        logical :: file_exists, nearly_full, finished, all_finished
        logical :: file_end(rdm%sign_length)
        character(30) :: rdm_filename(rdm%sign_length)
        character(len=*), parameter :: t_r = 'read_2rdm_popsfile'

        write(stdout, '(1X,"Reading in the spinfree 2-RDMs...")')

        ! If we're about to fill up the spawn list, perform a communication.
        nearly_full = .false.
        ! Have we finished adding RDM elements to the spawned list?
        finished = .false.
        ! Have we reached the end of the various files?
        file_end = .false.

        associate(state_labels => rdm_defs%state_labels, repeat_label => rdm_defs%repeat_label)
            do irdm = 1, rdm%sign_length
                if (state_labels(1, irdm) == state_labels(2, irdm)) then
                    write(rdm_filename(irdm), '("spinfree_TwoRDM.",'//int_fmt(state_labels(1, irdm), 0)//')') irdm
                else
                    write(rdm_filename(irdm), '("spinfree_TwoRDM.",'//int_fmt(state_labels(1, irdm), 0)//',"_",' &
                           //int_fmt(state_labels(2, irdm), 0)//',".",i1)') &
                        state_labels(1, irdm), state_labels(2, irdm), repeat_label(irdm)
                end if
                inquire (file=trim(rdm_filename(irdm)), exist=file_exists)

                if (.not. file_exists) then
                    call stop_all(t_r, "Attempting to read in a spinfree 2-RDM from "//trim(rdm_filename(irdm))//", &
                                       &but this file does not exist.")
                end if
            end do
        end associate

        ! Only let the root processor do the reading in.
        if (iProcIndex == 0) then
            do irdm = 1, rdm%sign_length
                iunit(irdm) = get_free_unit()
                open(iunit(irdm), file=trim(rdm_filename(irdm)), status='old')
            end do

            do
                do irdm = 1, rdm%sign_length
                    if (.not. file_end(irdm)) then
                        rdm_sign = 0.0_dp
                        read(iunit(irdm), "(4I15, F30.20)", iostat=ierr) i, j, k, l, rdm_sign(irdm)
                        if (ierr > 0) call stop_all(t_r, "Error reading "//trim(rdm_filename(irdm))//".")
                        ! The final line is printed as follows, by convention.
                        if (i == -1 .and. j == -1 .and. k == -1 .and. l == -1) then
                            file_end(irdm) = .true.
                            exit
                        end if

                        ! If the spawned list is nearly full, perform a communication.
                        if (nearly_full) then
                            call communicate_rdm_spawn_t_wrapper(spawn, rdm, finished, all_finished)
                            nearly_full = .false.
                        end if

                        call add_to_rdm_spawn_t(spawn, i, j, k, l, rdm_sign, .true., nearly_full)
                    end if
                end do
                if (all(file_end)) exit
            end do

            do irdm = 1, rdm%sign_length
                close(iunit(irdm))
            end do
        end if

        finished = .true.
        ! Keep performing communications until all RDM spawnings on every
        ! processor have been communicated.
        do
            call communicate_rdm_spawn_t_wrapper(spawn, rdm, finished, all_finished)
            if (all_finished) exit
        end do

        call annihilate_rdm_list(rdm)

        ! Fill in the hash table to the RDM. Clear it first, just in case.
        call clear_hash_table(rdm%hash_table)
        do ielem = 1, rdm%nelements
            call calc_separate_rdm_labels(rdm%elements(0, ielem), ij, kl, i, j, k, l)
            hash_val = FindWalkerHash((/i, j, k, l/), size(rdm%hash_table))
            call add_hash_table_entry(rdm%hash_table, ielem, hash_val)
        end do

        ! Clear the spawn object.
        spawn%free_slots = spawn%init_free_slots(0:nProcessors - 1)
        call clear_hash_table(spawn%rdm_send%hash_table)

    end subroutine read_spinfree_2rdm_files