WriteDoubHisttofile Subroutine

    subroutine WriteDoubHisttofile()

        integer :: i, j, k, l, BinNo, iunit
        real(dp) :: MaxFII, MinFII, BinIter, OnePartOrbEnValue, BinVal

        ! Histogramming all coulomb terms <ij|ij> where i<j, and i and j are
        ! both virtual. In reality we are looking at i = <j, but the ERhistograms
        ! will show the i = j terms.
        if (tROHistVirtCoulomb) then

            ROHistDCijOcklVir(:, :) = 0.0_dp
            MinFII = FourIndInts(1, 2, NoOcc + 1, NoOcc + 2)
            MaxFII = FourIndInts(1, 2, NoOcc + 1, NoOcc + 2)
            do i = 1, NoOcc
                do j = 1, NoOcc
                    do k = NoOcc + 1, NoOrbs
                        do l = NoOcc + 1, NoOrbs
                            if (FourIndInts(i, j, k, l) < MinFII) MinFII = FourIndInts(i, j, k, l)
                            if (FourIndInts(i, j, k, l) > MaxFII) MaxFII = FourIndInts(i, j, k, l)
                        end do
                    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
                ROHistDCijOcklVir(1, i) = BinVal
                BinVal = BinVal + BinIter
            end do
            do i = 1, NoOcc
                do j = 1, NoOcc
                    do k = NoOcc + 1, NoOrbs
                        do l = NoOcc + 1, NoOrbs
                            if (.not. near_zero(FourIndInts(i, j, k, l))) then
                                BinNo = CEILING((FourIndInts(i, j, k, l) - MinFII) * 4002 / (MaxFII - MinFII))
                                ROHistDCijOcklVir(2, BinNo) = ROHistDCijOcklVir(2, BinNo) + 1.0
                            end if
                        end do
                    end do
                end do
            end do

            ! Antisymmetric.
            ROHistASijOcklVir(:, :) = 0.0_dp
            MinFII = FourIndInts(1, 2, NoOcc + 1, NoOcc + 2) - FourIndInts(1, 2, NoOcc + 2, NoOcc + 1)
            MaxFII = FourIndInts(1, 2, NoOcc + 1, NoOcc + 2) - FourIndInts(1, 2, NoOcc + 2, NoOcc + 1)
            do i = 1, NoOcc
                do j = 1, NoOcc
                    do k = NoOcc + 1, NoOrbs
                        do l = NoOcc + 1, NoOrbs
                            if ((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) < MinFII) MinFII = (FourIndInts(i, j, k, l) - &
                                                                                                        FourIndInts(i, j, l, k))
                            if ((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) > MaxFII) MaxFII = (FourIndInts(i, j, k, l) - &
                                                                                                        FourIndInts(i, j, l, k))
                        end do
                    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
                ROHistASijOcklVir(1, i) = BinVal
                BinVal = BinVal + BinIter
            end do
            do i = 1, NoOcc
                do j = 1, NoOcc
                    do k = NoOcc + 1, NoOrbs
                        do l = NoOcc + 1, NoOrbs
                            if (.not. (FourIndInts(i, j, k, l) .isclose.FourIndInts(i, j, l, k))) then
                                BinNo = CEILING(((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) - MinFII) * 4002 / (MaxFII - MinFII))
                                ROHistASijOcklVir(2, BinNo) = ROHistASijOcklVir(2, BinNo) + 1.0
                            end if
                        end do
                    end do
                end do
            end do

            if (Iteration == 0) then
                iunit = get_free_unit()
                open(iunit, file='HistHFDoubijOcklVir', status='unknown')
                do j = 1, 4002
                    if (.not. (near_zero(ROHistDCijOcklVir(2, j)) .and. near_zero(ROHistASijOcklVir(2, j)))) then
                        write(iunit, '(4F20.10)') ROHistDCijOcklVir(1, j), ROHistDCijOcklVir(2, j), &
                            ROHistASijOcklVir(1, j), ROHistASijOcklVir(2, j)
                    end if
                end do
                close(iunit)
            end if
            if ((.not. tNotConverged) .and. (Iteration > 1)) then
                iunit = get_free_unit()
                open(iunit, file='HistRotDoubijOcklVir', status='unknown')
                do j = 1, 4002
                    if (.not. (near_zero(ROHistDCijOcklVir(2, j)) .and. near_zero(ROHistASijOcklVir(2, j)))) then
                        write(iunit, '(4F20.10)') ROHistDCijOcklVir(1, j), ROHistDCijOcklVir(2, j), &
                            ROHistASijOcklVir(1, j), ROHistASijOcklVir(2, j)
                    end if
                end do
                close(iunit)
            end if

            ROHistDCijklVir(:, :) = 0.0_dp
            MinFII = FourIndInts(NoOrbs - 1, NoOrbs, NoOrbs - 1, NoOrbs)
            MaxFII = FourIndInts(NoOrbs - 1, NoOrbs, NoOrbs - 1, NoOrbs)
            do i = NoOcc + 1, NoOrbs
                do j = NoOcc + 1, NoOrbs
                    do k = i + 1, NoOrbs
                        do l = j + 1, NoOrbs
                            if (FourIndInts(i, j, k, l) < MinFII) MinFII = FourIndInts(i, j, k, l)
                            if (FourIndInts(i, j, k, l) > MaxFII) MaxFII = FourIndInts(i, j, k, l)
                        end do
                    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
                ROHistDCijklVir(1, i) = BinVal
                BinVal = BinVal + BinIter
            end do
            do i = NoOcc + 1, NoOrbs
                do j = NoOcc + 1, NoOrbs
                    do k = i + 1, NoOrbs
                        do l = j + 1, NoOrbs
                            if (.not. near_zero(FourIndInts(i, j, k, l))) then
                                BinNo = CEILING((FourIndInts(i, j, k, l) - MinFII) * 4002 / (MaxFII - MinFII))
                                ROHistDCijklVir(2, BinNo) = ROHistDCijklVir(2, BinNo) + 1.0
                            end if
                        end do
                    end do
                end do
            end do

            ! Antisymmetric.
            ROHistASijklVir(:, :) = 0.0_dp
            MinFII = FourIndInts(NoOrbs - 3, NoOrbs - 2, NoOrbs - 1, NoOrbs) - FourIndInts(NoOrbs - 3, NoOrbs - 2, NoOrbs, NoOrbs - 1)
            MaxFII = FourIndInts(NoOrbs - 3, NoOrbs - 2, NoOrbs - 1, NoOrbs) - FourIndInts(NoOrbs - 3, NoOrbs - 2, NoOrbs, NoOrbs - 1)
            do i = NoOcc + 1, NoOrbs
                do j = NoOcc + 1, NoOrbs
                    do k = i + 1, NoOrbs
                        do l = j + 1, NoOrbs
                            if ((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) < MinFII) MinFII = (FourIndInts(i, j, k, l) - &
                                                                                                        FourIndInts(i, j, l, k))
                            if ((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) > MaxFII) MaxFII = (FourIndInts(i, j, k, l) - &
                                                                                                        FourIndInts(i, j, l, k))
                        end do
                    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
                ROHistASijklVir(1, i) = BinVal
                BinVal = BinVal + BinIter
            end do
            do i = NoOcc + 1, NoOrbs
                do j = NoOcc + 1, NoOrbs
                    do k = i + 1, NoOrbs
                        do l = j + 1, NoOrbs
                            if (.not. (FourIndInts(i, j, k, l) .isclose.FourIndInts(i, j, l, k))) then
                                BinNo = CEILING(((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) - MinFII) * 4002 / (MaxFII - MinFII))
                                ROHistASijklVir(2, BinNo) = ROHistASijklVir(2, BinNo) + 1.0
                            end if
                        end do
                    end do
                end do
            end do

            if (Iteration == 0) then
                iunit = get_free_unit()
                open(iunit, file='HistHFDoubijklVirt', status='unknown')
                do j = 1, 4002
                    if (.not. (near_zero(ROHistDCijklVir(2, j)) .and. near_zero(ROHistASijklVir(2, j)))) then
                        write(iunit, '(4F20.10)') ROHistDCijklVir(1, j), ROHistDCijklVir(2, j), &
                            ROHistASijklVir(1, j), ROHistASijklVir(2, j)
                    end if
                end do
                close(iunit)
            end if
            if ((.not. tNotConverged) .and. (Iteration > 1)) then
                iunit = get_free_unit()
                open(iunit, file='HistRotDoubijklVirt', status='unknown')
                do j = 1, 4002
                    if (.not. (near_zero(ROHistDCijklVir(2, j)) .and. near_zero(ROHistASijklVir(2, j)))) then
                        write(iunit, '(4F20.10)') ROHistDCijklVir(1, j), ROHistDCijklVir(2, j), ROHistASijklVir(1, j), &
                            ROHistASijklVir(2, j)
                    end if
                end do
                close(iunit)
            end if
        end if

        ! Histogramming all one particle orbital energies (occupied and
        ! virtual) even though we are not changing occupied.  Would like to
        ! see HOMO-LUMO gap etc.
        if (tROHistOneElInts) then

            ROHistHijVirt(:, :) = 0.0_dp
            MinFII = TMAT2DRot(NoOcc + 1, NoOcc + 2)
            MaxFII = TMAT2DRot(NoOcc + 1, NoOcc + 2)
            do i = NoOcc + 1, NoOrbs
                do j = i + 1, NoOrbs
                    if (TMAT2DRot(i, j) < MinFII) MinFII = TMAT2DRot(i, j)
                    if (TMAT2DRot(i, j) > MaxFII) MaxFII = TMAT2DRot(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
                ROHistHijVirt(1, i) = BinVal
                BinVal = BinVal + BinIter
            end do
            do i = NoOcc + 1, NoOrbs
                do j = i + 1, NoOrbs
                    if (.not. near_zero(TMAT2DRot(i, j))) then
                        BinNo = CEILING((TMAT2DRot(i, j) - MinFII) * 4002 / (MaxFII - MinFII))
                        ROHistHijVirt(2, BinNo) = ROHistHijVirt(2, BinNo) + 1.0
                    end if
                end do
            end do

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

            ROHistHijOccVirt(:, :) = 0.0_dp
            MinFII = TMAT2DRot(1, NoOcc + 1)
            MaxFII = TMAT2DRot(1, NoOcc + 1)
            do i = 1, NoOcc
                do j = NoOcc + 1, NoOrbs
                    if (TMAT2DRot(i, j) < MinFII) MinFII = TMAT2DRot(i, j)
                    if (TMAT2DRot(i, j) > MaxFII) MaxFII = TMAT2DRot(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
                ROHistHijOccVirt(1, i) = BinVal
                BinVal = BinVal + BinIter
            end do
            do i = 1, NoOcc
                do j = NoOcc + 1, NoOrbs
                    if (.not. near_zero(TMAT2DRot(i, j))) then
                        BinNo = CEILING((TMAT2DRot(i, j) - MinFII) * 4002 / (MaxFII - MinFII))
                        ROHistHijOccVirt(2, BinNo) = ROHistHijOccVirt(2, BinNo) + 1.0
                    end if
                end do
            end do

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

            ROHistHii(:, :) = 0.0_dp
            MinFII = TMAT2DRot(1, 1)
            MaxFII = TMAT2DRot(1, 1)
            do i = 1, NoOrbs
                if (TMAT2DRot(i, i) < MinFII) MinFII = TMAT2DRot(i, i)
                if (TMAT2DRot(i, i) > MaxFII) MaxFII = TMAT2DRot(i, i)
            end do
            BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
            MaxFII = MaxFII + BinIter
            MinFII = MinFII - BinIter
            BinVal = MinFII
            do i = 1, 4002
                ROHistHii(1, i) = BinVal
                BinVal = BinVal + BinIter
            end do
            do i = 1, NoOrbs
                BinNo = CEILING((TMAT2DRot(i, i) - MinFII) * 4002 / (MaxFII - MinFII))
                ROHistHii(2, BinNo) = ROHistHii(2, BinNo) + 1.0
            end do

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

        if (tROHistOnePartOrbEn) then
            ROHistOnePartOrbEn(:, :) = 0.0_dp
            MaxFII = 0.0_dp
            MinFII = 0.0_dp
            do i = 1, NoOrbs
                OnePartOrbEnValue = 0.0_dp
                OnePartOrbEnValue = OnePartOrbEnValue + TMAT2DRot(i, i)
                do j = 1, NoOcc
                    OnePartOrbEnValue = OnePartOrbEnValue + (2 * FourIndInts(i, j, i, j)) - FourIndInts(i, j, j, i)
                end do
                if (OnePartOrbEnValue > MaxFII) MaxFII = OnePartOrbEnValue
                if (OnePartOrbEnValue < MinFII) MinFII = OnePartOrbEnValue
            end do
            BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
            MaxFII = MaxFII + BinIter
            MinFII = MinFII - BinIter
            BinVal = MinFII
            do i = 1, 4002
                ROHistOnePartOrbEn(1, i) = BinVal
                BinVal = BinVal + BinIter
            end do
            do i = 1, NoOrbs
                OnePartOrbEnValue = 0.0_dp
                OnePartOrbEnValue = OnePartOrbEnValue + TMAT2DRot(i, i)
                do j = 1, NoOcc
                    OnePartOrbEnValue = OnePartOrbEnValue + (2 * FourIndInts(i, j, i, j)) - FourIndInts(i, j, j, i)
                end do
                BinNo = CEILING((OnePartOrbEnValue - MinFII) * 4002 / (MaxFII - MinFII))
                ROHistOnePartOrbEn(2, BinNo) = ROHistOnePartOrbEn(2, BinNo) + 1.0
            end do

            if (Iteration == 0) then
                iunit = get_free_unit()
                open(iunit, file='HistHFOnePartOrbEn', status='unknown')
                do j = 1, 4002
                    if (.not. near_zero(ROHistOnePartOrbEn(2, j))) then
                        do i = 1, 2
                            write(iunit, '(F20.10)', advance='no') ROHistOnePartOrbEn(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='HistRotOnePartOrbEn', status='unknown')
                do j = 1, 4002
                    if (.not. near_zero(ROHistOnePartOrbEn(2, j))) then
                        do i = 1, 2
                            write(iunit, '(F20.10)', advance='no') ROHistOnePartOrbEn(i, j)
                        end do
                        write(iunit, *) ''
                    end if
                end do
                close(iunit)
            end if
        end if

        if (tROHistDoubExc) then
            ROHistDoubExc(:, :) = 0.0_dp
            MaxFII = 0.0_dp
            MinFII = 0.0_dp
            do l = NoOcc + 1, NoOrbs
                do j = 1, NoOcc
                    do k = NoOcc + 1, NoOrbs
                        do i = 1, NoOcc
                            if ((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) > MaxFII) then
                                MaxFII = (FourIndInts(i, j, k, l)) - FourIndInts(i, j, l, k)
                            end if
                            if ((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) < MinFII) then
                                MinFII = (FourIndInts(i, j, k, l)) - FourIndInts(i, j, l, k)
                            end if
                        end do
                    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
                ROHistDoubExc(1, i) = BinVal
                BinVal = BinVal + BinIter
            end do
            do l = NoOcc + 1, NoOrbs
                do j = 1, NoOcc
                    do k = NoOcc + 1, NoOrbs
                        do i = 1, NoOcc
                            if (.not. (FourIndInts(i, j, k, l) .isclose.FourIndInts(i, j, l, k))) then
                                BinNo = CEILING(((FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)) - MinFII) * 4002 / (MaxFII - MinFII))
                                ROHistDoubExc(2, BinNo) = ROHistDoubExc(2, BinNo) + 1.0
                            end if
                        end do
                    end do
                end do
            end do

            if (Iteration == 0) then
                iunit = get_free_unit()
                open(iunit, file='HistHFDoubExc', status='unknown')
                do j = 1, 4002
                    if (.not. near_zero(ROHistDoubExc(2, j))) then
                        do i = 1, 2
                            write(iunit, '(F20.10)', advance='no') ROHistDoubExc(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='HistRotDoubExc', status='unknown')
                do j = 1, 4002
                    if (.not. near_zero(ROHistDoubExc(2, j))) then
                        do i = 1, 2
                            write(iunit, '(F20.10)', advance='no') ROHistDoubExc(i, j)
                        end do
                        write(iunit, *) ''
                    end if
                end do
                close(iunit)
            end if
        end if

        if (tROHistER) then
            ROHistER(:, :) = 0.0_dp
            MaxFII = 0.0_dp
            MinFII = 0.0_dp
            do i = 1, NoOrbs
                if (FourIndInts(i, i, i, i) > MaxFII) MaxFII = FourIndInts(i, i, i, i)
                if (FourIndInts(i, i, i, i) < MinFII) MinFII = FourIndInts(i, i, i, i)
            end do
            BinIter = ABS(MaxFII - MinFII) / 4000.0_dp
            MaxFII = MaxFII + BinIter
            MinFII = MinFII - BinIter
            BinVal = MinFII
            do i = 1, 4002
                ROHistER(1, i) = BinVal
                BinVal = BinVal + BinIter
            end do
            do i = 1, NoOrbs
                BinNo = CEILING((FourIndInts(i, i, i, i) - MinFII) * 4002 / (MaxFII - MinFII))
                ROHistER(2, BinNo) = ROHistER(2, BinNo) + 1.0
            end do

            if (Iteration == 0) then
                iunit = get_free_unit()
                open(iunit, file='HistHF-ER', status='unknown')
                do j = 1, 4002
                    if (.not. near_zero(ROHistER(2, j))) then
                        do i = 1, 2
                            write(iunit, '(F20.10)', advance='no') ROHistER(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='HistRot-ER', status='unknown')
                do j = 1, 4002
                    if (.not. near_zero(ROHistER(2, j))) then
                        do i = 1, 2
                            write(iunit, '(F20.10)', advance='no') ROHistER(i, j)
                        end do
                        write(iunit, *) ''
                    end if
                end do
                close(iunit)
            end if

        end if

    end subroutine WriteDoubHisttofile