Here, C(x k ) is the constraint (Eq. 39); x k is a generic variable in the orbital rotation space (Eq. 26) for the k-th optimization cycle; p k is the walking direction for x k , i.e. x k+1 = x k +αp k and the step length α is obtained from a line search; λ k+1 is the lagrange multiplier for the next (k + 1)-th optimization cycle. To reduce the computational cost, the hessian ∇ 2 xx L is approximated and updated by a BFGS scheme.