ShakeApproximation Subroutine

    subroutine ShakeApproximation()

        ! This is an approximation in which only the diagonal elements are
        ! considered in the matrix of the derivative of the constraints
        ! DerivConstrT1T2.

        integer :: m, l

        ! Use 'shake' algorithm in which the iterative scheme is applied to
        ! each constraint in succession.
        write(stdout, *) 'DerivConstrT1T2Diag calculated from the shake approx'

        DerivConstrT1T2Diag(:) = 0.0_dp
        do l = 1, TotNoConstraints
            do m = 1, NoOrbs
                DerivConstrT1T2Diag(l) = DerivConstrT1T2Diag(l) + Dot_Product(DerivConstrT2(:, m, l), DerivConstrT1(:, m, l))
            end do
            ShakeLambdaNew(l) = Constraint(l) / ((-1) * TimeStep * DerivConstrT1T2Diag(l))
            write(stdout, *) DerivConstrT1T2Diag(l)
        end do

    end subroutine ShakeApproximation