Cluster-Robust Standard Errors in Maximum Likelihood Estimation

I am estimating a model on pooled panel data by Maximum Likelihood using fminunc. I want to compute the cluster-robust standard errors after the estimation. This is a sandwich estimator, where the "bread" is given by the inverse Hessian and the "meat" involves the contribution of the k-th group to the score vector. I get the Hessian straight out of the fminunc algorithm.

But: How do I use the fminunc output to compute the contribution of the k-th group to the score vector?

Or, how would I alternatively compute the cluster-robust standard errors after minimizing the minus log-Likelihood function?