FindandUsetheForce Subroutine

    subroutine FindandUsetheForce(TotForce, TotDiffCoeffs, CoeffT2)

! This takes the current lambdas with the derivatives of the constraints and calculates a force
! for each cm, with an orthonormalisation correction.
! This is then used to rotate the coefficients by a defined timestep.

        integer :: a, m, Symm, w, SymMin, TempMaxOccVirt
        real(dp) :: TotForce, TotDiffCoeffs
        HElement_t(dp) :: CoeffT2(NoOrbs, NoOrbs)

        call set_timer(findandusetheforce_time, 30)

        if (tSeparateOccVirt) then
            TempMaxOccVirt = 2
        else
            TempMaxOccVirt = 1
        end if

        do w = 1, TempMaxOccVirt
            ! The force will be zero on those coefficients not being mixed,
            ! but still want to run over all, so that the diagonal 1 values
            ! are maintained.
            if (w == 1) then
                SymMin = 1
                MinMZ = 1
                if (tSeparateOccVirt) then
                    MaxMZ = NoOcc
                else
                    MaxMZ = NoOrbs
                end if
            else
                SymMin = 9
                MinMZ = NoOcc + 1
                MaxMZ = NoOrbs
            end if

            do m = MinMZ, MaxMZ
                if (tStoreSpinOrbs) then
                    SymM = int(G1(SymLabelList2_rot(m))%sym%S)
                else
                    SymM = int(G1(SymLabelList2_rot(m) * 2)%sym%S)
                end if
                do a = SymLabelCounts2_rot(1, SymM + SymMin), &
                    (SymLabelCounts2_rot(1, SymM + SymMin) + &
                     SymLabelCounts2_rot(2, SymM + SymMin) - 1)
!
                    ! FIND THE FORCE
                    ! find the corrected force. (in the case where the uncorrected force
                    !is required, correction is set to 0.
                    ! DerivCoeff(m,a) is the derivative of the relevant potential energy w.r.t
                    !cm without any constraints (no lambda terms).
                    ! ForceCorrect is then the latest force on coefficients.  This is
                    !iteratively being corrected so that
                    ! it will finally move the coefficients so that they remain orthonormal.

                    ! use THE FORCE
                    ForceCorrect(a, m) = DerivCoeff(a, m) - Correction(a, m)
                    CoeffT2(a, m) = CoeffT1(a, m) - (TimeStep * ForceCorrect(a, m))
                    ! Using the force to calculate the coefficients at time T2
                    ! (hopefully more orthonomal than those calculated in the
                    ! previous iteration).

                    ! Calculate parameters for printing
                    TotForce = TotForce + ABS(ForceCorrect(a, m))
                    TotDiffCoeffs = TotDiffCoeffs + ABS(CoeffT2(a, m) - CoeffT1(a, m))
                end do
            end do
        end do

        TotForce = TotForce / (real(NoOrbs**2, dp))

        call halt_timer(findandusetheforce_time)

    end subroutine FindandUsetheForce