attempt_detect_equil_time Subroutine

public subroutine attempt_detect_equil_time(this, equil_time, tFail)

Arguments

Type IntentOptional Attributes Name
real(kind=dp), allocatable :: this(:)
integer :: equil_time
logical, intent(out) :: tFail

Contents


Source Code

    subroutine attempt_detect_equil_time(this, equil_time, tFail)
        ! General routine, does not require global data
        ! Tries to remove the equilibration time from the start of the array

        real(dp), allocatable :: this(:)
        integer :: length
        real(dp), allocatable :: tmp(:)
        integer, allocatable :: i_tmp(:)
        logical :: swapped
        logical, intent(out) :: tFail
        real(dp) :: store
        integer :: i, i_store
        integer :: equil_time
        character(len=*), parameter :: t_r = 'attempt_detect_equil_time'

        equil_time = 0
        if (Errordebug > 0) write(stdout, *) "Attempting to detect equilibration time"

        if (.not. allocated(this)) call stop_all(t_r, "Error, array not allocated on entry into " &
            & //"attempt_remove_equil_time")
        length = size(this, 1)
        allocate(tmp(length))
        allocate(i_tmp(length))
        tmp = this

        do i = 1, length
            i_tmp(i) = i
        end do

        ! bubble sort - yay
        do
            swapped = .false.
            do i = 2, length
                if (tmp(i) < tmp(i - 1)) then
                    store = tmp(i - 1)
                    tmp(i - 1) = tmp(i)
                    tmp(i) = store
                    i_store = i_tmp(i - 1)
                    i_tmp(i - 1) = i_tmp(i)
                    i_tmp(i) = i_store
                    swapped = .true.
                end if
            end do
            if (.not. swapped) exit
        end do

!        do i=1,length
        !write(stdout,*) tmp(i),i_tmp(i) ! this just prints the sorted data
        ! useful to check routine does what is expected
!        end do

        ! A crude way to remove equilibration time automatically
        ! is to remove the first value (labelled 1st in i_tmp) in the list
        ! if it's ranked first or last in the sorted list
        ! Then to repeat this process.
        ! This can be done without resorting the list because now the 2nd element
        ! must lie at the (n-1)th position or the 2nd position (depending on where
        ! the first element was).

        tFail = .false.

        if (i_tmp(1) == 1) then
            if (Errordebug > 0) write(stdout, *) "Equilib time detected at lowest value"
            do i = 1, length
                if (i_tmp(i) /= i) then
                    if (Errordebug > 0) write(stdout, *) "Equilib time detected to be", i - 1, "values"
                    equil_time = i - 1
                    exit
                end if
                if (i == length) then
                    if (Errordebug > 0) then
                        write(stdout, *) "ERROR: the whole data file would be removed by automatic equilibrium detection"
                    end if
                    tFail = .true.
                    return
                end if
            end do
        end if
        if (i_tmp(length) == 1) then
            if (Errordebug > 0) write(stdout, *) "Equilib time detected at the highest value"
            do i = 1, length
                if (i_tmp(length - i + 1) /= i) then ! Yuck, integer adding, prone to error
                    if (Errordebug > 0) write(stdout, *) "Equilib time detected to be", i - 1, "values"
                    equil_time = i - 1
                    exit
                end if
                if (i == length) then
                    if (Errordebug > 0) then
                        write(stdout, *) "ERROR: the whole data file would be removed by automatic equilibrium detection"
                    end if
                    tFail = .true.
                    return
                end if
            end do
        end if

    end subroutine attempt_detect_equil_time