$\newcommand{\De}{\Delta} \newcommand{\ep}{\epsilon} \newcommand{\ga}{\gamma} \newcommand{\Ga}{\Gamma} \newcommand{\la}{\lambda} \newcommand{\Si}{\Sigma} \renewcommand{\th}{\theta} \newcommand{\R}{\mathbb{R}} \newcommand{\F}{\mathcal{F}} \newcommand{\E}{\operatorname{\mathsf E}} \newcommand{\PP}{\operatorname{\mathsf P}} \newcommand{\ii}[1]{\operatorname{\mathsf I}\{#1\}}$

Knowing only the 2nd and 3rd moments of the (presumable centered) summand random variables (r.v.'s) is hardly enough for good approximation.

Also, it is somewhat unclear what you mean by "3 moments". In the Berry--Esseen bound, the absolute 3rd moments are used, whereas the Edgeworth expansion (which is in a sense more accurate) is given in terms of initial moments.

So, given only the limited information you have, here is a cheap and fast way to approximate and simulate the sum $S_n$ of independent centered r.v.'s $X_1,\dots,X_n$ with given $\E X_i^2$ and $\E X_i^3$: just match the 2nd and 3rd (initial) moments of $S_n$ with the corresponding moments of a r.v. of the form $cZ$, where $c$ is a real number, $Z=Y-\E Y=Y-\la$, and $Y$ has the Poisson distribution with some parameter $\la>0$. Note that $\E S_n^2=B_2:=\sum_1^n \E X_i^2$, $\E S_n^3=B_3:=\sum_1^n \E X_i^3$, and $\E Z^2=\la=\E Z^3$. So, the "matching" system of equations is \begin{equation} c^2\la=B_2,\quad c^3\la=B_3. \end{equation} Solving this system, we get \begin{equation} c=\frac{B_3}{B_2},\quad \la=\frac{B_2^3}{B_3^2}. \end{equation} With these values of $c$ and $\la$, we can approximately simulate $S_n$ as $c(Y-\la)$ with $Y\sim\text{Poisson}(\la)$.

$\newcommand{\De}{\Delta} \newcommand{\ep}{\epsilon} \newcommand{\ga}{\gamma} \newcommand{\Ga}{\Gamma} \newcommand{\la}{\lambda} \newcommand{\Si}{\Sigma} \renewcommand{\th}{\theta} \newcommand{\R}{\mathbb{R}} \newcommand{\F}{\mathcal{F}} \newcommand{\E}{\operatorname{\mathsf E}} \newcommand{\PP}{\operatorname{\mathsf P}} \newcommand{\ii}[1]{\operatorname{\mathsf I}\{#1\}}$

Knowing only the 2nd and 3rd moments of the (presumable centered) summand random variables (r.v.'s) is hardly enough for very good approximation.

Also, it is somewhat unclear what you mean by "3 moments". In the Berry--Esseen bound, the absolute 3rd moments are used, whereas the Edgeworth expansion (which is in a sense more accurate) is given in terms of initial moments.

So, given only the limited information you have, here is a fast way to approximate and simulate the sum $S_n$ of independent centered r.v.'s $X_1,\dots,X_n$ with given $\E X_i^2$ and $\E X_i^3$: just match the 2nd and 3rd (initial) moments of $S_n$ with the corresponding moments of a r.v. of the form $cZ$, where $c$ is a real number, $Z=Y-\E Y=Y-\la$, and $Y$ has the Poisson distribution with some parameter $\la>0$. Note that $\E S_n^2=B_2:=\sum_1^n \E X_i^2$, $\E S_n^3=B_3:=\sum_1^n \E X_i^3$, and $\E Z^2=\la=\E Z^3$. So, the "matching" system of equations is \begin{equation} c^2\la=B_2,\quad c^3\la=B_3. \end{equation} Solving this system, we get \begin{equation} c=c_*:=\frac{B_3}{B_2},\quad \la=\la_*:=\frac{B_2^3}{B_3^2}. \end{equation} So, we can approximately simulate $S_n$ as $c_*(Y-\la_*)$ with $Y\sim\text{Poisson}(\la_*)$.

One may note that in the case when the $X_i$'s are identically distributed, $\la_*$ is inversely proportional to the coefficient in the leading term of the Edgeworth expansion. More generally, if the $\E X_i^2$'s are of the same order of magnitude and if the $\E X_i^3$'s are also of the same order of magnitude, then $\la_*$ is of the order of $n$; so, $\la_*$ is then large if $n$ is large, and then the distributions of $S_n$ and its Poisson approximation $c_*(Y-\la_*)$ are both close to $N(0,B_2)$ and hence to each other.