# FinalizeNewOrbs Subroutine

None

## Source Code

    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
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