Manually calculating robust standard errors for pglm using gradient and hessian matrix

Question

I want to manually calculate robust standard error for a fixed effect poisson model produced using the pglm function that, unlike the plm function, does support sandwich error matrices. This make impossible to use the standard vcovHC() function to calculate standard errors, since it does know how to extract terms in a format it understands from the result of pglm (see Robust Standard Errors: Poisson Panel Regression (pglm, lmtest) and https://stats.stackexchange.com/questions/273152/vcovhc-heteroskedasticity-in-pooled-and-panel-probit).

Looking on the web, I found the following way to calculate cluster robust standard errors:

library(readstata13)
library(pglm)
library(lmtest)
library(MASS)

ships<-readstata13::read.dta13("http://www.stata-press.com/data/r13/ships.dta")
ships$lnservice=log(ships$service)
res1 <- pglm(accident ~ op_75_79+co_65_69+co_70_74+co_75_79+lnservice,family = poisson, data = ships, effect = "individual", model="within", index = "ship")
summary(res1)

standard_se<-ginv(-res1$hessian)
coeftest(res1,standard_se)

# Similar to e(sample) in STATA
esample<-as.numeric(rownames(as.matrix(res1$gradientObs)))
fc <- ships[esample,]$ship #isolates the groups used in estimation

# Calculates the new Meat portion of our covariance matrix
m <- length(unique(fc))
k <- 5
u <- res1$gradientObs
u.clust <- matrix(NA, nrow=m, ncol=k)
for(j in 1:k){
  u.clust[,j] <- tapply(u[,j], fc, sum)
}
cl.vcov <-ginv(-res1$hessian)%*%( t(u.clust) %*% (u.clust))%*%ginv(-res1$hessian)
coeftest(res1,cl.vcov)

However, what I need are robust and not cluster robust standard errors. Could anyone shed some light on how to do that manually using gradient and hessian matrix?

Any help would be much appreciated. Thank you!