EquateDiagFock Subroutine

    subroutine EquateDiagFock()

        integer :: irr, NumInSym, Orbi, Orbj, w, i, j, k, ConjInd, OrbjConj
        real(dp) :: Angle, AngleConj, Check, Norm

        CoeffT1(:, :) = 0.0_dp

        ! Do virtual and occupied orbitals seperately.
        do w = MinOccVirt, MaxOccVirt

            ! Loop over irreps.
            do irr = 1, 8

                NumInSym = SymLabelCounts2_rot(2, (w - 1) * 8 + irr)

                ! Loop over the j-orthogonal vectors to create in this
                ! symmetry block.
                do j = 1, NumInSym

                    Orbj = SymLabelList2_rot(SymLabelCounts2_rot(1, (w - 1) * 8 + irr) - 1 + j)

                    ! See if this vector has already been done.
                    Check = 0.0_dp
                    do i = 1, NoOrbs
                        Check = Check + CoeffT1(i, Orbj)
                    end do
                    if (.not. near_zero(Check)) then
                        ! This vector is a conjugate pair of another vector and
                        ! has already been worked out...
                        cycle
                    end if

                    ! Find out if we this vector will be complex. It will be
                    ! real if j = N or j = N/2
                    if (j == NumInSym) then
                        !i The vector will be the normalized 1, 1, 1 vector.

                        do i = 1, NumInSym
                            Orbi = SymLabelList2_rot(SymLabelCounts2_rot(1, (w - 1) * 8 + irr) - 1 + i)
                            CoeffT1(Orbi, Orbj) = 1 / SQRT(real(NumInSym, dp))
                        end do

                    else if ((mod(NumInSym, 2) == 0) .and. (j == (NumInSym / 2))) then

                        do i = 1, NumInSym
                            Orbi = SymLabelList2_rot(SymLabelCounts2_rot(1, (w - 1) * 8 + irr) - 1 + i)
                            if (mod(i, 2) == 1) then
                                CoeffT1(Orbi, Orbj) = -1 / SQRT(real(NumInSym, dp))
                            else
                                CoeffT1(Orbi, Orbj) = 1 / SQRT(real(NumInSym, dp))
                            end if
                        end do

                    else
                        ! Vector is complex - find its conjugate vector - do
                        ! these at the same time.
                        ConjInd = NumInSym - j
                        OrbjConj = SymLabelList2_rot(SymLabelCounts2_rot(1, (w - 1) * 8 + irr) - 1 + ConjInd)

                        do i = 1, NumInSym

                            Orbi = SymLabelList2_rot(SymLabelCounts2_rot(1, (w - 1) * 8 + irr) - 1 + i)

                            Angle = real(i * j * 2, dp) * PI / real(NumInSym, dp)
                            AngleConj = real(i * ConjInd * 2, dp) * PI / real(NumInSym, dp)

                            CoeffT1(Orbi, Orbj) = (1 / SQRT(real(2 * NumInSym, dp))) * (COS(Angle) + COS(AngleConj))
                            CoeffT1(Orbi, OrbjConj) = (1 / SQRT(real(2 * NumInSym, dp))) * (SIN(Angle) - SIN(AngleConj))

                        end do

                    end if

                end do

            end do
        end do

        do j = 1, NoOrbs
            Norm = 0.0_dp
            do i = 1, NoOrbs
                Norm = Norm + (CoeffT1(i, j)**2)
            end do
            if (near_zero(Norm)) then
                CoeffT1(j, j) = 1.0_dp
            end if
        end do

        do j = 1, NoOrbs
            do i = 1, NoOrbs
                write(stdout, "(G13.5)", advance='no') CoeffT1(j, i)
            end do
            write(stdout, *) ""
        end do

        !Check normalization.
        do j = 1, NoOrbs
            Norm = 0.0_dp
            do i = 1, NoOrbs
                Norm = Norm + (CoeffT1(i, j)**2)
            end do
            if (abs(Norm - 1.0_dp) > 1.0e-7_dp) then
                call Stop_All("EquateDiagFock", "Rotation Coefficients not normalized")
            end if
        end do

        ! Check orthogonality.
        do j = 1, NoOrbs
            do i = 1, NoOrbs
                if (i == j) cycle
                Norm = 0.0_dp
                do k = 1, NoOrbs
                    Norm = Norm + (CoeffT1(k, j) * CoeffT1(k, i))
                end do
                if (abs(Norm) > 1.0e-7_dp) then
                    write(stdout, *) "COLUMNS: ", j, i
                    call Stop_All("EquateDiagFock", "RotationCoefficients not orthogonal")
                end if
            end do
        end do

    end subroutine EquateDiagFock