UseTheForce Subroutine

    subroutine UseTheForce()

        ! This routine takes the old translation coefficients and Lambdas and moves them by a timestep in the direction
        ! of the calculated force.

        integer :: m, w, z, i, j, Symm, SymMin
        real(dp) :: NewCoeff, NewLambda

        DistCs = 0.0_dp

        do w = MinOccVirt, MaxOccVirt
            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

                ! Symmetry requirement that z must be from the same irrep as m.
                do z = SymLabelCounts2_rot(1, SymM + SymMin), &
                    (SymLabelCounts2_rot(1, SymM + SymMin) + &
                     SymLabelCounts2_rot(2, SymM + SymMin) - 1)

                    ! Only coeffs with sym of m and z the same have non-zero coeffs.
                    NewCoeff = 0.0_dp
                    NewCoeff = CoeffT1(z, m) - (TimeStep * DerivCoeff(z, m))
                    DistCs = DistCs + abs(TimeStep * DerivCoeff(z, m))
                    CoeffT1(z, m) = NewCoeff
                end do
            end do
        end do

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

        if (tLagrange) then

            DistLs = 0.0_dp
            LambdaMag = 0.0_dp
            do i = 1, NoOrbs
                do j = 1, NoOrbs
                    NewLambda = 0.0_dp
                    NewLambda = Lambdas(i, j) - (TimeStep * DerivLambda(i, j))  ! Timestep must be specified in the input file.
                    DistLs = DistLs + abs(TimeStep * DerivLambda(i, j))
                    Lambdas(i, j) = NewLambda
                    LambdaMag = LambdaMag + abs(NewLambda)
                end do
            end do
            DistLs = DistLs / (real(NoOrbs**2, dp))
            LambdaMag = LambdaMag / (real(NoOrbs**2, dp))

        end if

    end subroutine UseTheForce