PrintIntegrals Subroutine

    subroutine PrintIntegrals()

        integer :: i, j, k, l, io1, io2
        real(dp) :: DiagOneElPot, ERPot, ijVirtOneElPot, ijVirtCoulPot, ijVirtExchPot
        real(dp) :: singCoulijVirt, singExchijVirt, singCoulconHF, singExchconHF, ijklPot, ijklantisymPot
        real(dp) :: ijOccVirtOneElPot, ijOccVirtCoulPot, ijOccVirtExchPot

        io1 = 0
        io2 = 0
        if (tInitIntValues) then
            io1 = get_free_unit()
            open(io1, file='DiagIntegrals', status='unknown')
            write(io1, '(A10, 6A18)') "Iteration", "<i|h|i> ivirt", "<ii|ii> ivirt", "<ij|ij> iOccjVirt", "<ij|ji> iOccjVirt", &
                "<ij|ij> ijVirt", "<ij|ji> ijVirt"

            io2 = get_free_unit()
            open(io2, file='SingExcIntegrals', status='unknown')
            write(io2, '(A10, 6A18)') "Iteration", "<i|h|j> iOccjVirt", "<i|h|j> ijVirt", "<ik|jk> HFcon", "<ik|kj> HFcon", &
                "<ik|jk> ijVirt", "<ik|kj> ijVirt"

            DiagOneElPotInit = 0.0_dp
            ERPotInit = 0.0_dp
            ijVirtOneElPotInit = 0.0_dp
            ijVirtCoulPotInit = 0.0_dp
            ijVirtExchPotInit = 0.0_dp
            singCoulconHFInit = 0.0_dp
            singExchconHFInit = 0.0_dp
            singCoulijVirtInit = 0.0_dp
            singExchijVirtInit = 0.0_dp
            ijklPotInit = 0.0_dp
            ijklantisymPotInit = 0.0_dp
            ijOccVirtOneElPotInit = 0.0_dp
            ijOccVirtCoulPotInit = 0.0_dp
            ijOccVirtExchPotInit = 0.0_dp
            NoInts01 = 0
            NoInts02 = 0
            NoInts03 = 0
            NoInts04 = 0
            NoInts05 = 0
            NoInts06 = 0
            do i = 1, NoOrbs
                if (i > NoOcc) then
                    DiagOneElPotInit = DiagOneElPotInit + TMAT2DRot(i, i)
                    ERPotInit = ERPotInit + FourIndInts(i, i, i, i)
                    NoInts01 = NoInts01 + 1
                    do j = NoOcc + 1, NoOrbs
                        ! The i, j terms with i and j both virtual.
                        if (j > i) then
                            ijVirtOneElPotInit = ijVirtOneElPotInit + TMAT2DRot(i, j)
                            ijVirtCoulPotInit = ijVirtCoulPotInit + FourIndInts(i, j, i, j)
                            ijVirtExchPotInit = ijVirtExchPotInit + FourIndInts(i, j, j, i)
                            NoInts02 = NoInts02 + 1
                        end if
                        do k = 1, NoOrbs
                            if (k > (NoOcc + 1)) then
                                do l = NoOcc + 1, NoOrbs
                                    if (l == j) cycle
                                    ijklPotInit = ijklPotInit + FourIndInts(i, j, k, l)
                                    ijklantisymPotInit = ijklantisymPotInit + FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)
                                    NoInts04 = NoInts04 + 1
                                end do
                            else
                                if (i == j) cycle
                                singCoulijVirtInit = singCoulijVirtInit + FourIndInts(i, k, j, k)
                                singExchijVirtInit = singExchijVirtInit + FourIndInts(i, k, k, j)
                                NoInts03 = NoInts03 + 1
                            end if
                        end do
                    end do
                else
                    do j = NoOcc + 1, NoOrbs
                        do k = 1, NoOcc
                            singCoulconHFInit = singCoulconHFInit + FourIndInts(i, k, j, k)
                            singExchconHFInit = singExchconHFInit + FourIndInts(i, k, k, j)
                            NoInts06 = NoInts06 + 1
                        end do
                        ijOccVirtOneElPotInit = ijOccVirtOneElPotInit + TMAT2DRot(i, j)
                        ijOccVirtCoulPotInit = ijOccVirtCoulPotInit + FourIndInts(i, j, i, j)
                        ijOccVirtExchPotInit = ijOccVirtExchPotInit + FourIndInts(i, j, j, i)
                        NoInts05 = NoInts05 + 1
                    end do
                end if
            end do
        end if

        DiagOneElPot = 0.0_dp
        ERPot = 0.0_dp
        ijVirtOneElPot = 0.0_dp
        ijVirtCoulPot = 0.0_dp
        ijVirtExchPot = 0.0_dp
        singCoulconHF = 0.0_dp
        singExchconHF = 0.0_dp
        singCoulijVirt = 0.0_dp
        singExchijVirt = 0.0_dp
        ijklPot = 0.0_dp
        ijklantisymPot = 0.0_dp
        ijOccVirtOneElPot = 0.0_dp
        ijOccVirtCoulPot = 0.0_dp
        ijOccVirtExchPot = 0.0_dp
        do i = 1, NoOrbs
            if (i > NoOcc) then
                DiagOneElPot = DiagOneElPot + TMAT2DRot(i, i)
                ERPot = ERPot + FourIndInts(i, i, i, i)
                do j = NoOcc + 1, NoOrbs
                    ! The i, j terms with i and j both virtual.
                    if (j > i) then
                        ijVirtOneElPot = ijVirtOneElPot + TMAT2DRot(i, j)
                        ijVirtCoulPot = ijVirtCoulPot + FourIndInts(i, j, i, j)
                        ijVirtExchPot = ijVirtExchPot + FourIndInts(i, j, j, i)
                    end if
                    do k = 1, NoOrbs
                        if (k > (NoOcc + 1)) then
                            do l = NoOcc + 1, NoOrbs
                                if (l == j) cycle
                                ijklPot = ijklPot + FourIndInts(i, j, k, l)
                                ijklantisymPot = ijklantisymPot + FourIndInts(i, j, k, l) - FourIndInts(i, j, l, k)
                            end do
                        else
                            if (i == j) cycle
                            singCoulijVirt = singCoulijVirt + FourIndInts(i, k, j, k)
                            singExchijVirt = singExchijVirt + FourIndInts(i, k, k, j)
                        end if
                    end do
                end do
            else
                do j = NoOcc + 1, NoOrbs
                    do k = 1, NoOcc
                        singCoulconHF = singCoulconHF + FourIndInts(i, k, j, k)
                        singExchconHF = singExchconHF + FourIndInts(i, k, k, j)
                    end do
                    ijOccVirtOneElPot = ijOccVirtOneElPot + TMAT2DRot(i, j)
                    ijOccVirtCoulPot = ijOccVirtCoulPot + FourIndInts(i, j, i, j)
                    ijOccVirtExchPot = ijOccVirtExchPot + FourIndInts(i, j, j, i)
                end do
            end if
        end do

        DiagOneElPot = (DiagOneElPot - DiagOneElPotInit) / NoInts01
        ERPot = (ERPot - ERPotInit) / NoInts01
        ijVirtOneElPot = (ijVirtOneElPot - ijVirtOneElPotInit) / NoInts02
        ijVirtCoulPot = (ijVirtCoulPot - ijVirtCoulPotInit) / NoInts02
        ijVirtExchPot = (ijVirtExchPot - ijVirtExchPotInit) / NoInts02
        singCoulijVirt = (singCoulijVirt - singCoulijVirtInit) / NoInts03
        singExchijVirt = (singExchijVirt - singExchijVirtInit) / NoInts03
        singCoulconHF = (singCoulconHF - singCoulconHFInit) / NoInts06
        singExchconHF = (singExchconHF - singExchconHFInit) / NoInts06
        ijklPot = (ijklPot - ijklPotInit) / NoInts04
        ijklantisymPot = (ijklantisymPot - ijklantisymPotInit) / NoInts04
        ijOccVirtOneElPot = (ijOccVirtOneElPot - ijOccVirtOneElPotInit) / NoInts05
        ijOccVirtCoulPot = (ijOccVirtCoulPot - ijOccVirtCoulPotInit) / NoInts05
        ijOccVirtExchPot = (ijOccVirtExchPot - ijOccVirtExchPot) / NoInts05

        write(io1, '(I10, 6F18.10)') Iteration, DiagOneElPot, ERPot, ijOccVirtCoulPot, ijOccVirtExchPot, ijVirtCoulPot, &
            ijVirtExchPot
        write(io2, '(I10, 6F18.10)') Iteration, ijOccVirtOneElPot, ijVirtOneElPot, singCoulconHF, singExchconHF, singCoulijVirt, &
            singExchijVirt

        if ((.not. tNotConverged) .and. (.not. tInitIntValues)) then
            close(io1)
            close(io2)
        end if

    end subroutine PrintIntegrals