read_fcimcstats Subroutine

public subroutine read_fcimcstats(iShiftVary, tFailRead)

Uses

Arguments

Type IntentOptional Attributes Name
integer, intent(inout) :: iShiftVary
logical, intent(out) :: tFailRead

Contents

Source Code


Source Code

    subroutine read_fcimcstats(iShiftVary, tFailRead)
        use SystemData, only: tMolpro, MolproID, tMolproMimic
        integer, intent(inout) :: iShiftVary
        logical, intent(out) :: tFailRead
        character(len=1) :: readline
        character(len=*), parameter :: t_r = "read_fcimcstats"
        character(len=24) :: filename
        logical :: exists, tRefToZero
        integer :: eof, comments, i, ierr
        integer :: iunit, WalkersDiffProc
        real(dp) :: doubs, change, Ann, Died, Born, rewalkers, imwalkers
        real(dp) :: shift, rate, reproje, improje, reinstproje, iminstproje
        real(dp) :: AccRat, IterTime, FracFromSing, TotImagTime, HFShift, InstShift
        real(dp) :: denom, renum, imnum, normhf, norm_psi, curr_S2, Avshift, dud, tote
        real(dp) :: curr_S2_init, AbsProjE
        integer(int64) :: iters, validdata, datapoints, totdets
        real(dp), dimension(lenof_sign) :: insthf

        block
            ! We close the 'FCIMCStats' file that is opened in appending mode
            ! and reopen in reading mode.
            ! If this breaks other code on the technical level, this is good,
            ! because there is a logical error if data is appended to
            ! 'FCIMCStats' after doing the error analysis.
            integer :: file_id
            logical :: connected
            inquire (file='FCIMCStats', number=file_id, opened=connected)
            if (connected) close(file_id)
        end block

        !Open file (FCIMCStats or FCIQMCStats)
        iunit = get_free_unit()
        if (tMolpro .and. .not. tMolproMimic) then
            filename = 'FCIQMCStats_'//adjustl(MolproID)
            inquire (file=filename, exist=exists)
            if (.not. exists) then
                write(stdout, *) 'No FCIQMCStats file found for error analysis'
                tFailRead = .true.
                return
            end if
            open(iunit, file=filename, status='old', action='read', position='rewind')
            write(stdout, "(A)") "Reading back in FCIQMCStats datafile..."
        else
            filename = 'FCIMCStats'
            inquire (file='FCIMCStats', exist=exists)
            if (.not. exists) then
                write(stdout, *) 'No FCIMCStats file found for error analysis'
                tFailRead = .true.
                return
            end if
            open(iunit, file='FCIMCStats', status='old', action='read', position='rewind')
            write(stdout, "(A)") "Reading back in FCIMCStats datafile..."
            if (inum_runs == 2 .and. exists) then
                write(stdout, "(A)") " This is a DNECI run! But we just analyse the first FCIMCStats!"
            end if
        end if

        comments = 0
        datapoints = 0
        validdata = 0
        tRefToZero = .false.
        do while (.true.)

            read(iunit, "(A)", advance='no', iostat=eof) readline
            if (eof < 0) then
                exit
            else if (eof > 0) then
                call stop_all(t_r, "Error reading FCIMCStats file")
            end if

            if (readline /= '#') then
                !Valid data on line

                if (lenof_sign == 2 .and. inum_runs == 1) then
                    !complex fcimcstats
                    read(iunit, *, iostat=eof) &
                        iters, &                !1.
                        shift, &                              !2.
                        change, &   !3.
                        rate, &   !4.
                        rewalkers, &                       !5.
                        imwalkers, &                       !6.
                        reproje, &                !7.     real \sum[ nj H0j / n0 ]
                        improje, &                   !8.     Im   \sum[ nj H0j / n0 ]
                        reinstproje, &                 !9.
                        iminstproje, &                    !10.
                        tote, &       ! Tot.ProjE.iter (Re)
                        insthf(1), &                         !11.
                        insthf(lenof_sign), &                         !12.
                        doubs, &                         !13.
                        AccRat, &                               !14.
                        TotDets, &                        !15.
                        IterTime, &                             !16.
                        FracFromSing, &     !17.
                        WalkersDiffProc, &                           !18.
                        TotImagTime, &                               !19.
                        HFShift, &                                   !20.
                        InstShift, &                                 !21.
                        denom     !24     |n0|^2  This is the denominator for both calcs
                else
                    read(iunit, *, iostat=eof) &
                        Iters, &
                        shift, &
                        change, &
                        rate, &
                        rewalkers, &
                        Ann, &
                        Died, &
                        Born, &
                        reproje, &
                        Avshift, &
                        reinstproje, &
                        insthf(1), &
                        doubs, &
                        AccRat, &
                        totdets, &
                        IterTime, &
                        FracFromSing, &
                        WalkersDiffProc, &
                        TotImagTime, &
                        dud, &
                        HFShift, &
                        InstShift, &
                        tote, &
                        denom
                end if

                if (eof < 0) then
                    call stop_all(t_r, "Should not be at end of file")
                else if (eof > 0) then
                    ! This is normally due to a difficulty reading NaN or
                    ! Infinity. Assume that this occurs when NoAtRef --> 0.
                    ! Therefore we can safely wipe the stats.
                    if (iters > iShiftVary) then
                        tRefToZero = .true.
                        validdata = 0
                        iShiftVary = int(Iters) + 1
                    end if
                else
                    if (iters > iShiftVary) then
                        if (abs(denom) < 1.0e-5_dp) then
                            !Denominator gone to zero - wipe the stats
                            tRefToZero = .true.
                            validdata = 0
                            iShiftVary = int(Iters) + 1
                        else
                            validdata = validdata + 1
                        end if
                    end if
                end if
                datapoints = datapoints + 1

            else
                !Just read it again without the advance=no to move onto the next line
                read(iunit, "(A)", iostat=eof) readline
                if (eof < 0) then
                    call stop_all(t_r, "Should not be at end of file")
                else if (eof > 0) then
                    call stop_all(t_r, "Error reading FCIMCStats file")
                end if
                comments = comments + 1
            end if
        end do

        if (tRefToZero) then
            write(stdout, "(A)") "After shift varies, reference population goes to zero"
            write(stdout, "(A,I14)") "Blocking will only start after non-zero population, at iteration ", iShiftVary - 1
        end if

        write(stdout, "(A,I12)") "Number of comment lines found in file: ", comments
        write(stdout, "(A,I12)") "Number of data points found in file: ", datapoints
        write(stdout, "(A,I12)") "Number of useable data points: ", validdata

        if (validdata <= 1) then
            write(stdout, "(A)") "No valid datapoints found in file. Exiting error analysis."
            tFailRead = .true.
            return
        else
            tFailRead = .false.
        end if

        !Allocate arrays
        allocate(numerator_data(validdata), stat=ierr)
        if (lenof_sign == 2 .and. inum_runs == 1) then
            allocate(imnumerator_data(validdata), stat=ierr)
        end if
        allocate(pophf_data(validdata), stat=ierr)
        allocate(shift_data(validdata), stat=ierr)
        if (ierr /= 0) call stop_all(t_r, "Memory allocation error")
        numerator_data(:) = 0.0_dp
        pophf_data(:) = 0.0_dp
        shift_data(:) = 0.0_dp
        rewind (iunit)

        !Now read all the data in again
        i = 0
        do while (.true.)

            read(iunit, "(A)", advance='no', iostat=eof) readline
            if (eof < 0) then
                exit
            else if (eof > 0) then
                call stop_all(t_r, "Error reading FCIMCStats file")
            end if

            if (readline /= '#') then
                !Valid data on line

                if (lenof_sign == 2 .and. inum_runs == 1) then
                    ! complex fcimcstats
                    read(iunit, *, iostat=eof) &
                        iters, &                !1.
                        shift, &                              !2.
                        change, &   !3.
                        rate, &   !4.
                        rewalkers, &                       !5.
                        imwalkers, &                       !6.
                        reproje, &                !7.     real \sum[ nj H0j / n0 ]
                        improje, &                   !8.     Im   \sum[ nj H0j / n0 ]
                        reinstproje, &                 !9.
                        iminstproje, &                    !10.
                        tote, &     ! Tot.PorjE.iter (Rm)
                        insthf(1), &                         !11.
                        insthf(lenof_sign), &                         !12.
                        doubs, &                         !13.
                        AccRat, &                               !14.
                        TotDets, &                        !15.
                        IterTime, &                             !16.
                        FracFromSing, &     !17.
                        WalkersDiffProc, &                           !18.
                        TotImagTime, &                               !19.
                        HFShift, &                                   !20.
                        InstShift, &                                 !21.
                        denom, &     !24     |n0|^2  This is the denominator for both calcs
                        renum, &   !22.    Re[\sum njH0j] x Re[n0] + Im[\sum njH0j] x Im[n0]   No div by StepsSft
                        imnum, &       !23.    Im[\sum njH0j] x Re[n0] - Re[\sum njH0j] x Im[n0]   since no physicality
                        normhf, &
                        norm_psi, &
                        curr_S2
                else
                    read(iunit, *, iostat=eof) &
                        Iters, &
                        shift, &
                        change, &
                        rate, &
                        rewalkers, &
                        Ann, &
                        Died, &
                        Born, &
                        reproje, &
                        Avshift, &
                        reinstproje, &
                        insthf(1), &
                        doubs, &
                        AccRat, &
                        totdets, &
                        IterTime, &
                        FracFromSing, &
                        WalkersDiffProc, &
                        TotImagTime, &
                        dud, &
                        HFShift, &
                        InstShift, &
                        tote, &
                        denom, &
                        renum, &
                        normhf, &
                        norm_psi, &
                        curr_S2, curr_S2_init, &
                        AbsProjE
                end if

                if (eof < 0) then
                    call stop_all(t_r, "Should not be at end of file")
                else if (eof > 0) then
                    ! This is normally due to a difficulty reading NaN or
                    ! Infinity. Assume that this occurs when NoAtRef --> 0.
                    ! Therefore, we should not be trying to store the stats,
                    ! and can ignore the error.
                    if (iters > iShiftVary) &
                        call stop_all(t_r, "This is not a valid data line - &
                                            &should not be trying to store")
                end if
                if (iters > iShiftVary) then
                    i = i + 1
                    if (tGivenEquilibrium) then
                        !Count the update cycles we are to ignore
                        if (iters < Shift_Start) then
                            relaxation_time_shift = relaxation_time_shift + 1
                        end if
                        if (iters < ProjE_Start) then
                            relaxation_time_proje = relaxation_time_proje + 1
                        end if
                    end if
                    numerator_data(i) = renum
                    if (lenof_sign == 2 .and. inum_runs == 1) then
                        imnumerator_data(i) = imnum
                    end if
                    pophf_data(i) = denom
                    shift_data(i) = shift
                end if
            else
                !Just read it again without the advance=no to move onto the next line
                read(iunit, "(A)", iostat=eof) readline
                if (eof < 0) then
                    call stop_all(t_r, "Should not be at end of file")
                else if (eof > 0) then
                    call stop_all(t_r, "Error reading FCIMCStats file")
                end if
            end if
        end do

        if (i /= validdata) then
            call stop_all(t_r, "Data arrays not filled correctly 1")
        end if

        if (abs(shift_data(validdata)) < 1.0e-10_dp .or. abs(pophf_data(validdata)) < 1.0e-10_dp .or. &
            abs(numerator_data(validdata)) < 1.0e-10_dp) then
            write(stdout, *) shift_data(validdata), pophf_data(validdata), numerator_data(validdata)
            call stop_all(t_r, "Data arrays not filled correctly 2")
        end if

        close(iunit)

    end subroutine read_fcimcstats