Sampling correlated truncated normal random variables

I am trying to do a simulation which requires generation of random variables. This simulation takes k>2 correlated truncated normal random variables. I have the covariance matrix and the means of those variables.

As a concrete example, I have three truncated random variables X, Y, Z mean is given by c(1,2,3). The truncation happens at +/-1.5*mean and covariance matrixis given bysigma=matrix(c(1,1,8.1,1,16,18,8.1,18,81),ncol=3)`.

Now, I have to sample random numbers using R such that I get the final correlation among the generated random variables given by the cov2cor(sigma). I understand that it will be an approximation and near match is also fine. I tried using tmvtnorm package but the final correlation is very different as the input covariance structure sigmais for the non-truncated normal distribution instead of the truncated one.

In a nutshell, I need to generate truncated normal random variables with a target correlation matrix using R.

Thanks!