CalcFOCKMatrix Subroutine

    subroutine CalcFOCKMatrix()

        use SystemData, only: nBasis
        use LoggingData, only: tRDMonfly

        integer :: i, j, k, l, a, b, ierr
        real(dp) :: FOCKDiagSumHF, FOCKDiagSumNew
        character(len=*), parameter :: this_routine = 'CalcFOCKMatrix'
        !NEED TO FIX THIS!

! This subroutine calculates and writes out the fock matrix for the transformed orbitals.
! ARR is originally the fock matrix in the HF basis.
! ARR(:,1) - ordered by energy, ARR(:,2) - ordered by spin-orbital index.

! When transforming the orbitals into approximate natural orbitals, we want to save memory, so don't bother
! calculating the whole matrix, just the diagonal elements that we actually need.

        if (tUseMP2VarDenMat .or. tFindCINatOrbs .or. tUseHFOrbs .or. tRDMonfly) then
            allocate(ArrDiagNew(NoOrbs), stat=ierr)
            call LogMemAlloc('ArrDiagNew', NoOrbs, 8, this_routine, ArrDiagNewTag, ierr)
            ArrDiagNew(:) = 0.0_dp
        else
            allocate(ArrNew(NoOrbs, NoOrbs), stat=ierr)
            call LogMemAlloc('ArrNew', NoOrbs**2, 8, this_routine, ArrNewTag, ierr)
            ArrNew(:, :) = 0.0_dp
        end if

        ! First calculate the sum of the diagonal elements, ARR.
        ! Check if this is already being done.
        FOCKDiagSumHF = 0.0_dp
        do a = 1, nBasis
            FOCKDiagSumHF = FOCKDiagSumHF + Arr(a, 2)
        end do

        write(stdout, *) 'Sum of the fock matrix diagonal elements in the HF basis set = ', FOCKDiagSumHF

        FOCKDiagSumNew = 0.0_dp
        do j = 1, NoRotOrbs
            l = SymLabelList3_rot(j)
            if (tUseMP2VarDenMat .or. tFindCINatOrbs .or. tUseHFOrbs .or. tRDMonfly) then
                do a = 1, NoOrbs
                    b = SymLabelList2_rot(a)
                    if (tStoreSpinOrbs .or. tTurnStoreSpinOff) then
                        ArrDiagNew(l) = ArrDiagNew(l) + (CoeffT1(a, j) * ARR(b, 2) * CoeffT1(a, j))
                    else
                        ArrDiagNew(l) = ArrDiagNew(l) + (CoeffT1(a, j) * ARR(2 * b, 2) * CoeffT1(a, j))
                    end if
                end do
                if (tStoreSpinOrbs .or. tTurnStoreSpinOff) then
                    FOCKDiagSumNew = FOCKDiagSumNew + (ArrDiagNew(l))
                else
                    FOCKDiagSumNew = FOCKDiagSumNew + (ArrDiagNew(l) * 2)
                end if
            else
                do i = 1, NoRotOrbs
                    k = SymLabelList2_rot(i)
                    ArrNew(k, l) = 0.0_dp
                    do a = 1, NoOrbs
                        b = SymLabelList2_rot(a)
                        if (tStoreSpinOrbs .or. tTurnStoreSpinOff) then
                            ArrNew(k, l) = ArrNew(k, l) + (CoeffT1(a, i) * Arr(b, 2) * CoeffT1(a, j))
                        else
                            ArrNew(k, l) = ArrNew(k, l) + (CoeffT1(a, i) * Arr(2 * b, 2) * CoeffT1(a, j))
                        end if
                    end do
                end do
                if (tStoreSpinOrbs .or. tTurnStoreSpinOff) then
                    FOCKDiagSumNew = FOCKDiagSumNew + (ArrNew(l, l))
                else
                    FOCKDiagSumNew = FOCKDiagSumNew + (ArrNew(l, l) * 2)
                end if
                ! Only running through spat orbitals, count each twice to
                ! compare to above.
            end if
        end do
        ! If we are truncation the virtual space, only the unfrozen entries
        ! will be transformed.

        write(stdout, *) 'Sum of the fock matrix diagonal elements in the transformed basis set = ', FOCKDiagSumNew

! Refill ARR(:,1) (ordered in terms of energies), and ARR(:,2) (ordered in terms of orbital number).
! ARR(:,2) needs to be ordered in terms of symmetry and then energy (like SymLabelList), so currently this ordering will not be
! correct when reading in qchem intDUMPS as the orbital number ordering is by energy.

        ! If we are only writing out 1 ROFCIDUMP or we are not truncating at
        ! all - can refill ARR etc.
        if (NoDumpTruncs <= 1) then

            if (tUseMP2VarDenMat .or. tFindCINatOrbs .or. tUseHFOrbs .or. tRDMonfly) then
                if (tStoreSpinOrbs .or. tTurnStoreSpinOff) then
                    do j = 1, NoOrbs
                        ARR(j, 2) = ArrDiagNew(j)
                        ARR(j, 1) = ArrDiagNew(BRR(j))
                    end do
                else
                    do j = 1, NoOrbs
                        ARR(2 * j, 2) = ArrDiagNew(j)
                        ARR(2 * j - 1, 2) = ArrDiagNew(j)
                        ARR(2 * j, 1) = ArrDiagNew(BRR(2 * j) / 2)
                        ARR(2 * j - 1, 1) = ArrDiagNew(BRR(2 * j) / 2)
                    end do
                end if
            else
                if (tStoreSpinOrbs .or. tTurnStoreSpinOff) then
                    do j = 1, NoRotOrbs
                        ARR(j, 2) = ArrNew(j, j)
                        ARR(j, 1) = ArrNew(BRR(j), BRR(j))
                    end do
                else
                    do j = 1, NoRotOrbs
                        ARR(2 * j, 2) = ArrNew(j, j)
                        ARR(2 * j - 1, 2) = ArrNew(j, j)
                        ARR(2 * j, 1) = ArrNew(BRR(2 * j) / 2, BRR(2 * j) / 2)
                        ARR(2 * j - 1, 1) = ArrNew(BRR(2 * j) / 2, BRR(2 * j) / 2)
                    end do
                end if
            end if

        end if

        if ((tUseMP2VarDenMat .or. tFindCINatOrbs .or. tUseHFOrbs .or. tRDMonfly) .and. (NoDumpTruncs <= 1)) then
            deallocate(ArrDiagNew)
            call LogMemDealloc(this_routine, ArrDiagNewTag)
        else if (NoDumpTruncs <= 1) then
            deallocate(ArrNew)
            call LogMemDealloc(this_routine, ArrNewTag)
        end if

        write(stdout, *) 'end of calcfockmatrix'
        call neci_flush(stdout)

    end subroutine CalcFOCKMatrix