How to get covariance matrix for random effects from LME4

So, I've fitted a LMM with two random intercepts in R: Y = Xbeta + Zb + e_i, where b ~ N (0, sigma).

I would like to get my hands on the underlying covariance matrix of b, which doesn't seem to be a trivial thing in lme4 package. You can get only the variances by VarCorr, not the actual correlation matrix.

According to this (page 2):

http://ift.tt/1V4R9rQ you can calculate the covariance of beta: e_i * lambda * t(lambda). And all those components you can extract from the output of lme4.

I was wondering if this is the way to go? Or do you have any other suggestions?