WriteSingHisttofile Subroutine

    subroutine WriteSingHisttofile()

        integer :: i, j, k, BinNo, a, b, iunit
        real(dp) :: MaxFII, MinFII, BinIter, BinVal, SingExcit(NoOrbs, NoOrbs)

        ! <ik|jk> terms where all i, j and k are virtual
        ! Coulomb.
        ROHistSCijkVir(:, :) = 0.0_dp
        MaxFII = FourIndInts(NoOcc + 1, NoOcc + 1, NoOcc + 2, NoOcc + 1)
        MinFII = FourIndInts(NoOcc + 1, NoOcc + 1, NoOcc + 2, NoOcc + 1)
        do i = NoOcc + 1, NoOrbs
            do k = NoOcc + 1, NoOrbs
                do j = i + 1, NoOrbs
                    if (FourIndInts(i, k, j, k) > MaxFII) MaxFII = FourIndInts(i, k, j, k)
                    if (FourIndInts(i, k, j, k) < MinFII) MinFII = FourIndInts(i, k, j, k)
                end do
            end do
        end do
        BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
        MaxFII = MaxFII + BinIter
        MinFII = MinFII - BinIter
        BinVal = MinFII
        do i = 1, 4002
            ROHistSCijkVir(1, i) = BinVal
            BinVal = BinVal + BinIter
        end do
        do i = NoOcc + 1, NoOrbs
            do k = NoOcc + 1, NoOrbs
                do j = i + 1, NoOrbs
                    if (.not. near_zero(FourIndInts(i, k, j, k))) then
                        BinNo = CEILING((FourIndInts(i, k, j, k) - MinFII) * 4002 / (MaxFII - MinFII))
                        ROHistSCijkVir(2, BinNo) = ROHistSCijkVir(2, BinNo) + 1.0
                    end if
                end do
            end do
        end do

        ! Exchange.
        ROHistSEijkVir(:, :) = 0.0_dp
        MaxFII = FourIndInts(NoOcc + 1, NoOcc + 1, NoOcc + 1, NoOcc + 2)
        MinFII = FourIndInts(NoOcc + 1, NoOcc + 1, NoOcc + 1, NoOcc + 2)
        do i = NoOcc + 1, NoOrbs
            do k = NoOcc + 1, NoOrbs
                do j = i + 1, NoOrbs
                    if (FourIndInts(i, k, k, j) > MaxFII) MaxFII = FourIndInts(i, k, k, j)
                    if (FourIndInts(i, k, k, j) < MinFII) MinFII = FourIndInts(i, k, k, j)
                end do
            end do
        end do
        BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
        MaxFII = MaxFII + BinIter
        MinFII = MinFII - BinIter
        BinVal = MinFII
        do i = 1, 4002
            ROHistSEijkVir(1, i) = BinVal
            BinVal = BinVal + BinIter
        end do
        do i = NoOcc + 1, NoOrbs
            do k = NoOcc + 1, NoOrbs
                do j = i + 1, NoOrbs
                    if (.not. near_zero(FourIndInts(i, k, k, j))) then
                        BinNo = CEILING((FourIndInts(i, k, k, j) - MinFII) * 4002 / (MaxFII - MinFII))
                        ROHistSEijkVir(2, BinNo) = ROHistSEijkVir(2, BinNo) + 1.0
                    end if
                end do
            end do
        end do

        ! Antisymmetric.
        ROHistSASijkVir(:, :) = 0.0_dp
        MaxFII = FourIndInts(NoOcc + 1, NoOcc + 1, NoOcc + 2, NoOcc + 1) - FourIndInts(NoOcc + 1, NoOcc + 1, NoOcc + 1, NoOcc + 2)
        MinFII = FourIndInts(NoOcc + 1, NoOcc + 1, NoOcc + 2, NoOcc + 1) - FourIndInts(NoOcc + 1, NoOcc + 1, NoOcc + 1, NoOcc + 2)
        do i = NoOcc + 1, NoOrbs
            do k = NoOcc + 1, NoOrbs
                do j = i + 1, NoOrbs
                if ((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) > MaxFII) MaxFII = FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)
                if ((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) < MinFII) MinFII = FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)
                end do
            end do
        end do
        BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
        MaxFII = MaxFII + BinIter
        MinFII = MinFII - BinIter
        BinVal = MinFII
        do i = 1, 4002
            ROHistSASijkVir(1, i) = BinVal
            BinVal = BinVal + BinIter
        end do
        do i = NoOcc + 1, NoOrbs
            do k = NoOcc + 1, NoOrbs
                do j = i + 1, NoOrbs
                    if (.not. (FourIndInts(i, k, j, k) .isclose.FourIndInts(i, k, k, j))) then
                        BinNo = CEILING(((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) - MinFII) * 4002 / (MaxFII - MinFII))
                        ROHistSASijkVir(2, BinNo) = ROHistSASijkVir(2, BinNo) + 1.0
                    end if
                end do
            end do
        end do

        if (Iteration == 0) then
            iunit = get_free_unit()
            open(iunit, file='HistHFSingijkVir', status='unknown')
            do j = 1, 4002
                if (any(.not. near_zero([ROHistSCijkVir(2, j), ROHistSEijkVir(2, j), ROHistSASijkVir(2, j)]))) then
                    write(iunit, '(6F20.10)') ROHistSCijkVir(1, j), ROHistSCijkVir(2, j), ROHistSEijkVir(1, j), ROHistSEijkVir(2, j),&
                                                        ROHistSASijkVir(1, j), ROHistSASijkVir(2, j)
                end if
            end do
            close(iunit)
        end if
        if ((Iteration > 1) .and. (.not. tNotConverged)) then
            iunit = get_free_unit()
            open(iunit, file='HistRotSingijkVir', status='unknown')
            do j = 1, 4002
                if (any(.not. near_zero([ROHistSCijkVir(2, j), ROHistSEijkVir(2, j), ROHistSASijkVir(2, j)]))) then
                    write(iunit, '(6F20.10)') ROHistSCijkVir(1, j), ROHistSCijkVir(2, j), ROHistSEijkVir(1, j), ROHistSEijkVir(2, j),&
                                                        ROHistSASijkVir(1, j), ROHistSASijkVir(2, j)
                end if
            end do
            close(iunit)
        end if

        ! <ik|jk> where k is occupied, and i and j are both virtual.
        ! Coulomb.
        ROHistSCkOcijVir(:, :) = 0.0_dp
        MaxFII = FourIndInts(NoOcc + 1, 1, NoOcc + 2, 1)
        MinFII = FourIndInts(NoOcc + 1, 1, NoOcc + 2, 1)
        do i = NoOcc + 1, NoOrbs
            do k = 1, NoOcc
                do j = i + 1, NoOrbs
                    if (FourIndInts(i, k, j, k) > MaxFII) MaxFII = FourIndInts(i, k, j, k)
                    if (FourIndInts(i, k, j, k) < MinFII) MinFII = FourIndInts(i, k, j, k)
                end do
            end do
        end do
        BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
        MaxFII = MaxFII + BinIter
        MinFII = MinFII - BinIter
        BinVal = MinFII
        do i = 1, 4002
            ROHistSCkOcijVir(1, i) = BinVal
            BinVal = BinVal + BinIter
        end do
        do i = NoOcc + 1, NoOrbs
            do k = 1, NoOcc
                do j = i + 1, NoOrbs
                    if (.not. near_zero(FourIndInts(i, k, j, k))) then
                        BinNo = CEILING((FourIndInts(i, k, j, k) - MinFII) * 4002 / (MaxFII - MinFII))
                        ROHistSCkOcijVir(2, BinNo) = ROHistSCkOcijVir(2, BinNo) + 1.0
                    end if
                end do
            end do
        end do

        ! Exchange.
        ROHistSEkOcijVir(:, :) = 0.0_dp
        MaxFII = FourIndInts(NoOcc + 1, 1, 1, NoOcc + 2)
        MinFII = FourIndInts(NoOcc + 1, 1, 1, NoOcc + 2)
        do i = NoOcc + 1, NoOrbs
            do k = 1, NoOcc
                do j = i + 1, NoOrbs
                    if (FourIndInts(i, k, k, j) > MaxFII) MaxFII = FourIndInts(i, k, k, j)
                    if (FourIndInts(i, k, k, j) < MinFII) MinFII = FourIndInts(i, k, k, j)
                end do
            end do
        end do
        BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
        MaxFII = MaxFII + BinIter
        MinFII = MinFII - BinIter
        BinVal = MinFII
        do i = 1, 4002
            ROHistSEkOcijVir(1, i) = BinVal
            BinVal = BinVal + BinIter
        end do
        do i = NoOcc + 1, NoOrbs
            do k = 1, NoOcc
                do j = i + 1, NoOrbs
                    if (.not. near_zero(FourIndInts(i, k, k, j))) then
                        BinNo = CEILING((FourIndInts(i, k, k, j) - MinFII) * 4002 / (MaxFII - MinFII))
                        ROHistSEkOcijVir(2, BinNo) = ROHistSEkOcijVir(2, BinNo) + 1.0
                    end if
                end do
            end do
        end do

        !antisymmetric
        ROHistSASkOcijVir(:, :) = 0.0_dp
        MaxFII = FourIndInts(NoOcc + 1, 1, NoOcc + 2, 1) - FourIndInts(NoOcc + 1, 1, 1, NoOcc + 2)
        MinFII = FourIndInts(NoOcc + 1, 1, NoOcc + 2, 1) - FourIndInts(NoOcc + 1, 1, 1, NoOcc + 2)
        do i = NoOcc + 1, NoOrbs
            do k = 1, NoOcc
                do j = i + 1, NoOrbs
                if ((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) > MaxFII) MaxFII = FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)
                if ((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) < MinFII) MinFII = FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)
                end do
            end do
        end do
        BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
        MaxFII = MaxFII + BinIter
        MinFII = MinFII - BinIter
        BinVal = MinFII
        do i = 1, 4002
            ROHistSASkOcijVir(1, i) = BinVal
            BinVal = BinVal + BinIter
        end do
        do i = NoOcc + 1, NoOrbs
            do k = 1, NoOcc
                do j = i + 1, NoOrbs
                    if (.not. (FourIndInts(i, k, j, k) .isclose.FourIndInts(i, k, k, j))) then
                        BinNo = CEILING(((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) - MinFII) * 4002 / (MaxFII - MinFII))
                        ROHistSASkOcijVir(2, BinNo) = ROHistSASkOcijVir(2, BinNo) + 1.0
                    end if
                end do
            end do
        end do

        if (Iteration == 0) then
            iunit = get_free_unit()
            open(iunit, file='HistHFSingkOcijVir', status='unknown')
            do j = 1, 4002
                if (any(.not. near_zero([ROHistSCkOcijVir(2, j), ROHistSEkOcijVir(2, j), ROHistSASkOcijVir(2, j)]))) then
                    write(iunit, '(6F20.10)') ROHistSCkOcijVir(1, j), ROHistSCkOcijVir(2, j), ROHistSEkOcijVir(1, j), &
                        ROHistSEkOcijVir(2, j), ROHistSASkOcijVir(1, j), ROHistSASkOcijVir(2, j)
                end if
            end do
            close(iunit)
        end if
        if ((Iteration > 1) .and. (.not. tNotConverged)) then
            iunit = get_free_unit()
            open(iunit, file='HistRotSingkOcijVir', status='unknown')
            do j = 1, 4002
                if (any(.not. near_zero([ROHistSCkOcijVir(2, j), ROHistSEkOcijVir(2, j), ROHistSASkOcijVir(2, j)]))) then
                    write(iunit, '(6F20.10)') ROHistSCkOcijVir(1, j), ROHistSCkOcijVir(2, j), ROHistSEkOcijVir(1, j), &
                        ROHistSEkOcijVir(2, j), ROHistSASkOcijVir(1, j), ROHistSASkOcijVir(2, j)
                end if
            end do
            close(iunit)
        end if

        ! <ik|jk> where i and k are both occupied, and j virtual.
        ! Coulomb.
        ROHistSCikOcjVir(:, :) = 0.0_dp
        MaxFII = FourIndInts(1, 1, NoOcc + 1, 1)
        MinFII = FourIndInts(1, 1, NoOcc + 1, 1)
        do i = 1, NoOcc
            do k = 1, NoOcc
                do j = NoOcc + 1, NoOrbs
                    if (FourIndInts(i, k, j, k) > MaxFII) MaxFII = FourIndInts(i, k, j, k)
                    if (FourIndInts(i, k, j, k) < MinFII) MinFII = FourIndInts(i, k, j, k)
                end do
            end do
        end do
        BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
        MaxFII = MaxFII + BinIter
        MinFII = MinFII - BinIter
        BinVal = MinFII
        do i = 1, 4002
            ROHistSCikOcjVir(1, i) = BinVal
            BinVal = BinVal + BinIter
        end do
        do i = 1, NoOcc
            do k = 1, NoOcc
                do j = NoOcc + 1, NoOrbs
                    if (.not. near_zero(FourIndInts(i, k, j, k))) then
                        BinNo = CEILING((FourIndInts(i, k, j, k) - MinFII) * 4002 / (MaxFII - MinFII))
                        ROHistSCikOcjVir(2, BinNo) = ROHistSCikOcjVir(2, BinNo) + 1.0
                    end if
                end do
            end do
        end do

        ! Exchange.
        ROHistSEikOcjVir(:, :) = 0.0_dp
        MaxFII = FourIndInts(1, 1, 1, NoOcc + 1)
        MinFII = FourIndInts(1, 1, 1, NoOcc + 1)
        do i = 1, NoOcc
            do k = 1, NoOcc
                do j = NoOcc + 1, NoOrbs
                    if (FourIndInts(i, k, k, j) > MaxFII) MaxFII = FourIndInts(i, k, k, j)
                    if (FourIndInts(i, k, k, j) < MinFII) MinFII = FourIndInts(i, k, k, j)
                end do
            end do
        end do
        BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
        MaxFII = MaxFII + BinIter
        MinFII = MinFII - BinIter
        BinVal = MinFII
        do i = 1, 4002
            ROHistSEikOcjVir(1, i) = BinVal
            BinVal = BinVal + BinIter
        end do
        do i = 1, NoOcc
            do k = 1, NoOcc
                do j = NoOcc + 1, NoOrbs
                    if (.not. near_zero(FourIndInts(i, k, k, j))) then
                        BinNo = CEILING((FourIndInts(i, k, k, j) - MinFII) * 4002 / (MaxFII - MinFII))
                        ROHistSEikOcjVir(2, BinNo) = ROHistSEikOcjVir(2, BinNo) + 1.0
                    end if
                end do
            end do
        end do

        ! Antisymmetrised.
        ROHistSASikOcjVir(:, :) = 0.0_dp
        MaxFII = FourIndInts(1, 1, NoOcc + 1, 1) - FourIndInts(1, 1, 1, NoOcc + 1)
        MinFII = FourIndInts(1, 1, NoOcc + 1, 1) - FourIndInts(1, 1, 1, NoOcc + 1)
        do i = 1, NoOcc
            do k = 1, NoOcc
                do j = NoOcc + 1, NoOrbs
                if ((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) > MaxFII) MaxFII = FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)
                if ((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) < MinFII) MinFII = FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)
                end do
            end do
        end do
        BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
        MaxFII = MaxFII + BinIter
        MinFII = MinFII - BinIter
        BinVal = MinFII
        do i = 1, 4002
            ROHistSASikOcjVir(1, i) = BinVal
            BinVal = BinVal + BinIter
        end do
        do i = 1, NoOcc
            do k = 1, NoOcc
                do j = NoOcc + 1, NoOrbs
                    if (.not. (FourIndInts(i, k, j, k) .isclose.FourIndInts(i, k, k, j))) then
                        BinNo = CEILING(((FourIndInts(i, k, j, k) - FourIndInts(i, k, k, j)) - MinFII) * 4002 / (MaxFII - MinFII))
                        ROHistSASikOcjVir(2, BinNo) = ROHistSASikOcjVir(2, BinNo) + 1.0
                    end if
                end do
            end do
        end do

        if (Iteration == 0) then
            iunit = get_free_unit()
            open(iunit, file='HistHFSingikOcjVir', status='unknown')
            do j = 1, 4002
                if (any(.not. near_zero([ROHistSCikOcjVir(2, j), ROHistSEikOcjVir(2, j), ROHistSASikOcjVir(2, j)]))) then
                    write(iunit, '(6F20.10)') ROHistSCikOcjVir(1, j), ROHistSCikOcjVir(2, j), ROHistSEikOcjVir(1, j), &
                        ROHistSEikOcjVir(2, j), ROHistSASikOcjVir(1, j), ROHistSASikOcjVir(2, j)
                end if
            end do
            close(iunit)
        end if
        if ((Iteration > 1) .and. (.not. tNotConverged)) then
            iunit = get_free_unit()
            open(iunit, file='HistRotSingikOcjVir', status='unknown')
            do j = 1, 4002
                if (any(.not. near_zero([ROHistSCikOcjVir(2, j), ROHistSEikOcjVir(2, j), ROHistSASikOcjVir(2, j)]))) then
                    write(iunit, '(6F20.10)') ROHistSCikOcjVir(1, j), ROHistSCikOcjVir(2, j), ROHistSEikOcjVir(1, j), &
                        ROHistSEikOcjVir(2, j), ROHistSASikOcjVir(1, j), ROHistSASikOcjVir(2, j)
                end if
            end do
            close(iunit)
        end if

        ! Single excitations connected to the HF determinant.
        ROHistSing(:, :) = 0.0_dp
        MaxFII = 0.0_dp
        MinFII = 0.0_dp
        do j = NoOcc + 1, NoOrbs
            do i = 1, NoOcc
                SingExcit(i, j) = 0.0_dp
                if (i == j) cycle
                a = SymLabelList2_rot(i)
                b = SymLabelList2_rot(j)
                do k = 1, NoOcc + 1
              SingExcit(i, j) = SingExcit(i, j) + real(TMAT2D(2 * a, 2 * b), dp) + ((2 * FourIndInts(i, k, j, k)) - FourIndInts(i, k, k, j))
                end do
                if (SingExcit(i, j) > MaxFII) MaxFII = SingExcit(i, j)
                if (SingExcit(i, j) < MinFII) MinFII = SingExcit(i, j)
            end do
        end do
        BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
        MaxFII = MaxFII + BinIter
        MinFII = MinFII - BinIter
        BinVal = MinFII
        do i = 1, 4002
            ROHistSing(1, i) = BinVal
            BinVal = BinVal + BinIter
        end do
        do j = NoOcc + 1, NoOrbs
            do i = 1, NoOcc
                if (i == j) cycle
                if (.not. near_zero(SingExcit(i, j))) then
                    BinNo = CEILING((SingExcit(i, j) - MinFII) * 4002 / (MaxFII - MinFII))
                    ROHistSing(2, BinNo) = ROHistSing(2, BinNo) + 1.0
                end if
            end do
        end do

        if (Iteration == 0) then
            iunit = get_free_unit()
            open(iunit, file='HistHFSingExcHF', status='unknown')
            do j = 1, 4002
                if (.not. near_zero(ROHistSing(2, j))) then
                    do i = 1, 2
                        write(iunit, '(F20.10)', advance='no') ROHistSing(i, j)
                    end do
                    write(iunit, *) ''
                end if
            end do
            close(iunit)
        end if
        if ((Iteration > 1) .and. (.not. tNotConverged)) then
            iunit = get_free_unit()
            open(iunit, file='HistRotSingExcHF', status='unknown')
            do j = 1, 4002
                if (.not. near_zero(ROHistSing(2, j))) then
                    do i = 1, 2
                        write(iunit, '(F20.10)', advance='no') ROHistSing(i, j)
                    end do
                    write(iunit, *) ''
                end if
            end do
            close(iunit)
        end if

    end subroutine WriteSingHisttofile