as far as i can remember, GM[X] = exp(mu) and GSD[X] = exp(sd), where mu = E[log X] and sd = sqrt(Var[log X]), i.e. the moments of the normally distributed Y = log X which are plugged in as meanlog and sdlog in your R-function.

so just call rlnorm with log GM[X] and log GSD[X] instead of GM[X] and GSD[X]:

> X  gm_mean = function(x, na.rm=TRUE){exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))}

example:

> X gm_mean(X)

[1] 5.109866

re-check before use! ;-)

Thank you Thomas.