FinalizeNewOrbs Subroutine

    subroutine FinalizeNewOrbs()

! At the end of the orbital rotation, have a set of coefficients CoeffT1 which transform
! the HF orbitals into a set of linear
! combinations ui which minimise |<ij|kl>|^2.  This is the final subroutine after
! all iterations (but before the memory deallocation)
! that calculates the final 4 index integrals to be used in the NECI calculation.

        use sym_mod, only: GenSymStatePairs

        integer :: i, a, j
        real(dp) :: TotGSConstraints, GSConstraint(TotNoConstraints)
        HElement_t(dp) :: CoeffTemp(SpatOrbs, SpatOrbs)

        ! First need to do a final explicit orthonormalisation. The orbitals
        ! are very close to being orthonormal, but not exactly. Need to make
        ! sure they are exact orthonormal using Gram Schmit.
        if (tStoreSpinOrbs .and. (.not. tMaxHLGap)) then
            CoeffTemp(:, :) = 0.0_dp
            do i = 1, SpatOrbs
                do j = 1, SpatOrbs
                    CoeffTemp(i, j) = CoeffT1(2 * i, 2 * j)
                end do
            end do

            call GRAMSCHMIDT_NECI(CoeffTemp, SpatOrbs)

            CoeffT1(:, :) = 0.0_dp
            do i = 1, SpatOrbs
                do j = 1, SpatOrbs
                    CoeffT1(2 * i, 2 * j) = CoeffTemp(i, j)
                    CoeffT1((2 * i) - 1, (2 * j) - 1) = CoeffTemp(i, j)
                end do
            end do
        else if (.not. tMaxHLGap) then
            call GRAMSCHMIDT_NECI(CoeffT1, NoOrbs)
        end if

! Put routine in here that takes this rotation matrix, CoeffT1, and forms raises it to the power of a small number, alpha.
! Changeing this number allows us to see the change in plateau level with various rotations.

! Write out some final results of interest, like values of the constraints, values of new coefficients.

        write(stdout, *) 'The final transformation coefficients after gram schmidt orthonormalisation'
        do i = 1, NoOrbs
            do a = 1, NoOrbs
                write(stdout, '(F10.4)', advance='no') CoeffT1(a, i)
            end do
            write(stdout, *) ''
        end do

        call WriteTransformMat()

        call CalcConstraints(CoeffT1, GSConstraint, TotGSConstraints)

        write(stdout, *) 'Final Potential Energy before orthogonalisation', PotEnergy

        call Transform2ElInts()

! Use these final coefficients to find the FourIndInts(i,j,k,l).
! These are now the <ij|kl> integrals we now want to use instead of the HF UMat.
! New potential energy is calculated in this routine using the orthogonalised coefficients.
! Compare to that before this, to make sure the orthogonalisation hasn't shifted them back to a non-minimal place.

        write(stdout, *) 'Final Potential Energy after orthogonalisation', PotEnergy

! Calculate the fock matrix, and print it out to see how much the off diagonal terms contribute.
! Also print out the sum of the diagonal elements to compare to the original value.
        call CalcFOCKMatrix()

        call RefillUMATandTMAT2D()
! UMat is the 4 index integral matrix (2 electron), whereas TMAT2D is the 2 index integral (1 el) matrix

! This is the keyword that tells the NECI calculation that the orbitals are not HF.  It means that contributions to
! the energy from walkers on singly occupied determinants are included in the values printed.
! Making it true here allows us to go directly from a Rotation into a spawn if required.
        tRotatedOrbs = .true.

        call GENSymStatePairs(SpatOrbs, .false.)

    end subroutine FinalizeNewOrbs