FindNewOrbs Subroutine

    subroutine FindNewOrbs()

        if (tERLocalization .and. (.not. tStoreSpinOrbs)) then
            call Transform2ElIntsERlocal()
        else
! Find the partially (and completely) transformed 4 index integrals to be used in further calcs.
            call Transform2ElInts()
        end if

!Find derivatives of the c and lambda matrices and print the sum of off-diagonal matrix elements.
        call FindTheForce()
        ! This finds the unconstrained force (unless the lagrange keyword is present).

!Update coefficents by moving them in direction of force. Print sum of squared changes in coefficients.
        if (tShake) then
            call ShakeConstraints()
            ! Find the force that moves the coefficients while keeping them orthonormal, and use it
            ! to get these new coefficients.
        else
            call UseTheForce()
            ! This can be either completely unconstrained, or have the lagrange constraints imposed.
        end if
!The coefficients coefft1(a,m) are now those that have been shifted by the time step.

!Test these for orthonomaility and then convergence.
!If they do not meet the convergence criteria, they go back into the previous step to produce another set of coefficients.

        call set_timer(testorthoconver_time, 30)

        call TestOrthonormality()
!Force should go to zero as we end in minimum - test for this

        call TestForConvergence()

        call halt_timer(testorthoconver_time)

    end subroutine FindNewOrbs