I have been trying to estimate Heckman's two-step solution by hand. (for purposes of comparison with ML and application in a real sample of relatively large size). Therefore, I have reviewed the functions of the sampleSelection package. 'vcov' returns the variance-covariance matrix of the estimated coefficients, also the package documentation indicates that 'vcov' contains the correct standard errors.
On the other hand 'heckitVcov' calculates the asymptotic covariance matrix for the coefficients of a Heckit estimation, that's why I tried to use the 'heckitVcov' function, entering the arguments as indicated, but the result is 'NA' in all elements. I have cheated a bit as a shortcut by using values from the same 'selection' object, to compute the asymptotic covariance matrix.
Am I making a mistake in interpreting and calculating heckiVcov? Thanks

data("Mroz87")  
mroz<-as.data.frame(Mroz87)  
mroz$kids <- (mroz$kids5 + mroz$kids618 > 0)  
mod_TS2 <- selection(selection = lfp ~ age + I(age^2) + faminc + kids + educ, outcome = wage ~ exper + I(exper^2) + educ + city, data=mroz, method="2step")

Calculating 'heckitVcov'

heckitVcov(xMat=model.matrix(mod_TS2), wMat=model.matrix(mod_TS2$probit), vcovProbit=vcov(mod_TS2$probit), rho=mod_TS2$rho, delta=mod_TS2$imrDelta, sigma=mod_TS2$sigma)

its result

##                XO(Intercept) XOexper XOI(exper^2) XOeduc XOcity imrData$IMR1
##XO(Intercept)            NA      NA           NA     NA     NA           NA
##XOexper                  NA      NA           NA     NA     NA           NA
##XOI(exper^2)             NA      NA           NA     NA     NA           NA
##XOcity                   NA      NA           NA     NA     NA           NA
##imrData$IMR1             NA      NA           NA     NA     NA           NA

Note: I have already seen a similar question Heckman selection model in R manually but it is not conclusive.