This is an updated version of my initial answer. The OP has raised several little issues in the comments that are now corrected. I also provide more details that are supposed to make it easier to follow/generalise the argument. I apologise for the length.
As pointed out in the comments such a formula is not explicitly in the literature. However, it can be derived from a general Voronoi type formula due to A. Corbett. Let me try to shed some light on how to do this in the following. For simplicity I will treat the modified sum $$S(X) = \sum_{n\in\mathbb{Z}} \frac{A_F(n,1)}{\vert n\vert}e\left(n\frac{a}{c}\right)W\left(\frac{n}{X}\right),$$ for a smooth function $W$ with compact support in $\mathbb{R}_+$. To derive this formula it will be useful to switch to the language of automorphic representations. So let us assume that $f$ and $F$ are newforms and let $\pi_f$ and $\Pi_F=\text{sym}^2(\pi_f)$ be the associated automorphic representations.
Before we can come to the summation formula in question we gather some preliminaries. Note that, because $N$ is square free and $f$ has trivial nebentypus, we know that $\pi_{f,p}=\text{St}$ is the Steinberg representation for all $p\mid N$. From this we can derive several useful facts. At the archimedean place we find that $\Pi_{F,\infty}$ has representation parameters $$\lambda=(2\nu,-2\nu,0) \text{ and } \delta=\begin{cases} (0,0,0) &\text{ in the Maa\ss\ case,}\\ (1,0,1) &\text{in the holomorphic case.}\end{cases} $$ In particular, $$ L_{\infty}(s,\Pi_{F,}) = \begin{cases}
\Gamma_{\mathbb{R}}(s+2\nu)\Gamma_{\mathbb{R}}(s)\Gamma_{\mathbb{R}}(s-2\nu) &\text{ in the Maa\ss\ case,} \\ \Gamma_{\mathbb{C}}(s+2k-1)\Gamma_{\mathbb{R}}(s+1) &\text{ in the holomorphic case.}
\end{cases}$$ See for example Proposition~5.12 in ["Summation formulas, from Poisson and Voronoi to the present"](https://arxiv.org/abs/math/0304187). The associated gamma-quotient then is
\begin{equation}
\gamma(s,\Pi_{F,\infty},\psi_{\infty}) = \epsilon(s,\Pi_{F,\infty},\psi_{\infty})\frac{L_{\infty}(1-s,\tilde{\Pi}_F)}{L_{\infty}(s,\Pi_F)}.\nonumber
\end{equation}
The computation at the finite places is slightly more involved. We start by mentioning that the Steinberg representation corresponds under the Local Langlands Correspondence to $\Vert \cdot \Vert^{-\frac{1}{2}}\otimes \text{sp}(2)$, where $\text{sp}(2)$ is the $2$-dimensional special representation of $W_{\mathbb{Q}_p}$. Further, using Artin reciprocity we can identify a character $\chi$ of $\mathbb{Q}_p^{\times}$ with a character of $W_{\mathbb{Q}_p}$. The upshot is that, we can compute the $L$- and $\epsilon$-factors of $\chi \Pi_{F,p}$ on the Weil-Deligne side. Indeed we find that
\begin{equation}
\left[\Vert \cdot \Vert^{-\frac{1}{2}}\otimes \text{sp}(2)\right]\otimes \left[\Vert \cdot \Vert^{-\frac{1}{2}}\otimes\text{sp}(2)\right] \cong 1 \oplus \left[ \Vert \cdot \Vert^{-1}\otimes\text{sp}(3)\right]. \nonumber
\end{equation}
From this we conclude that instead of computing the local factors of $\chi \Pi_{F,p}$ we can compute those of $\chi\Vert \cdot \Vert^{-1}\otimes \text{sp}(3)$ on the Weil-Deligne side of the Local Langlands Correspondence. According to ["Elliptic Curves and the Weil-Deligne group"](https://math.berkeley.edu/~dyott/Elliptic Curves and the Weil-Deligne Group.pdf) by D. E. Rohrlich we have that
\begin{multline}
a(\chi \Pi_{F,p}) = \begin{cases}
2 &\text{ if $\chi$ is unramified,}\\
3a(\chi) &\text{ if $\chi$ is ramified,}
\end{cases}, L(s,\Pi_{F,p}) = L(s+1,\chi) \\ \text{ and } \epsilon(\frac{1}{2},\chi\Pi_{F,p}) = \begin{cases}
1 &\text{ if $\chi=1$,}\\ \epsilon(\frac{1}{2},\chi)^3 &\text{ if $\chi$ is unitary and ramified.}
\end{cases} \nonumber
\end{multline}
(Here some clarification concerning the epsilon factors is necessary. We are always following Langlands convention and consider the canonical additive character with the corresponding self dual Haar measure as fixed and drop them from the notation.)
Our local computation implies the well known fact that the conductor of $\Pi_F$ is $N^2$. Further, we can write
\begin{multline}
A_F(m,n)=A_F\left(\frac{m}{(m,N^{\infty})},\frac{n}{(n,N^{\infty})}\right) (mn,N^{\infty})\\ \cdot \prod_{p\mid N} W_{F,p}\left(\left(\begin{matrix} (mn,p^{\infty})&0&0\\0&(n,p^{\infty})&0\\0&0&1\end{matrix}\right)\right), \nonumber
\end{multline}
where $W_{F,p}$ denotes the (suitably normalised) Whittaker new vector of $\Pi_{F,p}$. The Fourier coefficients $A_F(m,d)$ depend only on unramified data of $\pi_F$. In particular we have the following well known relations:
\begin{equation*}
A_F(m,1) = A_F(1,m) = \sum_{d^2\mid n} a_f(\frac{n^2}{d^4}) \text{ and } A_F(m,n) = \sum_{d\mid (m,n)}\mu(d)A_F(\frac{n}{d},1)A_F(\frac{m}{d},1),
\end{equation*}
for $(md,N)=1$. Knowing the local $L$-factor we can compute
\begin{multline}
W_{F,p}\left(\left(\begin{matrix} p^{f_1+f_2} & 0 & 0\\ 0 & p^{f_2} & 0 \\ 0 & 0 &1\end{matrix}\right)\right) = p^{-f_1-f_2}s_{(f_1+f_2,f_2)}(p^{-1},0)W_{F,p}(1) \\= \begin{cases} p^{-2f_1}W_{F,p}(1) &\text{ if }f_2=0, \\
0& \text{ if } f_2\neq 0.
\end{cases} \nonumber
\end{multline}
This can be deduced from Theorem 4.1 in ["Whittaker functions associated to newforms of GL(n)"](https://projecteuclid.org/euclid.jmsj/1390600834) by M. Miyauchi.
These preliminaries show that the case of symmetric square lifts (of forms of square free level) is very special. We now turn towards developing the requested Voronoi formula. Sine our method is flexible enough we will start by working in greater generality. Let $\Pi$ be an automorphic representation for $\text{GL}_3(\mathbb{A})$ of level $Q$ containing the newform $G$ of character $\chi$. (In the setting above we have $\Pi=\Pi_F$, $G=F$, $\chi=1$ and $Q=N^2$.)
We now make the first step towards our summation formula. By applying Theorem~1.1 from A. Corbett's paper ["Voronoi summation for $\text{GL}_n$: collusion between level and modulus"](https://arxiv.org/abs/1807.00716) with $\chi=1$, $c_2=1$, $\phi_{\infty}=W(\frac{\cdot}{X})$, $M=1$, $l=1$, $q=c$ we get
\begin{multline}
S_G(X) = \sum_{n\in\mathbb{Z}} \frac{A_G(n,1)}{\vert n\vert}e\left(n\frac{a}{c}\right)W\left(\frac{n}{X}\right)\\
=c\sum_{\substack{m,r\in\mathbb{Z}_{\neq 0},\\ (m,N)=1,\\ r\mid N^{\infty}}}\sum_{d\mid c}\text{KL}(\overline{aQL^3}r;m;c,1,d)\chi(\frac{md}{c})\frac{A_G(d,m)}{\vert md\vert} \\
\cdot \mathcal{B}_{\Pi_{\infty}, \phi_{\infty}}\left(\frac{rmd^2}{c^3QL^3}\right)\cdot \prod_{p\mid N} \mathcal{B}_{\Pi_{p}}\left(\frac{rmd^2}{c^3QL^3}\right), \nonumber
\end{multline}
where $L$ is the largest square free integer dividing $Q$. Here $\text{KL}$ is a 2-dimensional Kloosterman sum defined by $$\text{KL}_2(x,y;\frac{c}{d})\colon =\text{KL}(x;y;c,1,d) = \sum_{\alpha\in \left(\mathbb{Z}/\frac{c}{d}\mathbb{Z}\right)^{\times}}e\left(-\frac{xd\alpha}{c}+\frac{y\overline{\alpha}}{c/d}\right)$$ and the function $\mathcal{B}_{\Pi_{\infty},\phi_{\infty}}(\cdot)$ is a Bessel-transform of $\phi_{\infty}$ given by
\begin{multline}
\mathcal{B}_{\Pi_{\infty},\phi_{\infty}\infty}(y) = \frac{1}{4\pi i}\sum_{r=0,1} \text{sgn}(y)^r\int_{(\sigma)}\gamma(1-s,\text{sgn}^r\Pi_{\infty},\psi_{\infty})\vert y\vert^{1-s}\int_0^{\infty}W(\frac{x}{X})\vert x\vert^{-1-s}dxds \\
= \frac{1}{X} [\mathcal{H}_{\Pi_{\infty}} W] \left(\frac{y}{X}\right). \nonumber
\end{multline}
Here we have rewritten $\mathcal{B}_{\Pi_{\infty},\phi_{\infty}\infty}$ in terms of the transform \begin{equation}[\mathcal{H}_{\Pi_{\infty}} W](y) = \frac{1}{4\pi i}\sum_{r=0,1} \text{sgn}(y)^r\int_{(\sigma)}\gamma(1-s,\text{sgn}^r\Pi_{\infty},\psi_{\infty})\vert y\vert^{1-s}[\mathfrak{M}W](-s) ds, \label{C}\end{equation} where $[\mathfrak{M}W]$ is just the normal Mellin-transform of $W$.
The transforms $\mathcal{B}_{\Pi_{p},\phi_{p}}(\cdot)$ are $p$-adic versions of this Bessel transform. They are given by
\begin{multline}
\mathcal{B}_{\Pi_{p},\phi_{p}}(y)= \frac{\log(p)}{2\pi} \sum_{\substack{\xi\colon F^{\times}\to S^1,\\ \xi(p)=1}}\xi(y)\int_{\sigma-i\frac{\pi}{\log(p)}}^{\sigma+i\frac{\pi}{\log(p)}}\epsilon(\frac{1}{2},\xi\Pi_{p})p^{a(\xi \Pi_{p})(s-\frac{1}{2})}\frac{L(s,\xi^{-1}\tilde{\Pi}_{p})}{L(1-s,\chi\Pi_{p})}\vert y \vert_p^{1-s} \\ \cdot \int_{\mathbb{Q}_p^{\times}} \xi(x)\psi_p(x\frac{a}{c})W_{G,p}\left(\left(\begin{matrix}
x&0&0\\0&1&0\\0&0&1\end{matrix}\right)\right)\vert x\vert_p^{-s}d^{\times}x ds. \nonumber
\end{multline}
We claim that
\begin{equation}
\int_{\mathbb{Q}_p^{\times}} \xi(x)\psi_p(x\frac{a}{c})W_{G,p}\left(\left(\begin{matrix}
x&0&0\\0&1&0\\0&0&1\end{matrix}\right)\right)\vert x\vert_p^{-s}d^{\times}x ds = \begin{cases}
L(1-s,\Pi_p)W_{G,p}(1) &\text{ if }\xi=1,\\
0 &\text{ else.}
\end{cases} \nonumber
\end{equation}
To see this we first note that the support of $W_{G,p}$ on the matrices in question allows us to integrate over $\mathbb{Z}_p$ instead of $\mathbb{Q}_p$. Next, since $\psi_p$ is unramified and $(c,p)=1$ we have $\psi_p(x\frac{a}{c})=1$ for $x\in \mathbb{Z}_p$. Now if $\xi$ is non-trivial the integral must vanish for obvious reasons. On the other hand, if $\xi=1$, the claim follows from Theorem~5.1 in ["Whittaker functions associated to newforms of GL(n)"](https://projecteuclid.org/euclid.jmsj/1390600834).
With this claim being established the Bessel-transform takes the following nice form:
\begin{equation}
\mathcal{B}_{\Pi_{p},\phi_{p}}(y)= \epsilon(\frac{1}{2},\Pi_p)\frac{\log(p)}{2\pi}\int_{\sigma-i\frac{\pi}{\log(p)}}^{\sigma+i\frac{\pi}{\log(p)}}p^{a( \Pi_{p})(s-\frac{1}{2})}L(s,\tilde{\Pi}_{p})\vert y \vert_p^{1-s}dsW_{G,p}(1). \nonumber
\end{equation}
Writing $L(s,\tilde{\Pi}_p) = \sum_{k\geq 0} s_{(k,0)}(\tilde{\alpha})p^{-ks}$ we find
\begin{align}
\mathcal{B}_{\Pi_{p},\phi_{p}}(y) &= \delta_{v_p(y)\geq -a(\Pi_p)}\epsilon(\frac{1}{2},\Pi_p)\vert y \vert_p p^{-\frac{a(\Pi_p)}{2}} s_{(a(\Pi_p)+v_p(y))}(\tilde{\alpha})W_{G,p}(1) \nonumber \\
&= p^{\frac{a(\Pi_p)}{2}} \epsilon(\frac{1}{2},\Pi_p) W_{\tilde{G},p}\left(\left(\begin{matrix} p^{a(\Pi_p)+v_p(y)} & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1\end{matrix} \right)\right)\frac{W_{G,p}(1)}{W_{\tilde{G},p}(1)}. \nonumber
\end{align}
Here we applied Theorem~4.1 from ["Whittaker functions associated to newforms of GL(n)"](https://projecteuclid.org/euclid.jmsj/1390600834) to express the symmetric Polynomial in terms of the Whittaker vector for the contragredient newform $\tilde{G}$. (If $\Pi_p$ is self-dual we can replace $\tilde{G}$ by $G$ in the formula above. This applies in particular to $\Pi_p=\Pi_{F,p}$, where this can be checked more directly using the explicit computations from the beginning.)
With this at hand we can explicate our formula as follows:
\begin{multline}
S_G(X) = \frac{cQ^{\frac{1}{2}}}{X} \cdot \left[\prod_{p\mid Q}\epsilon(\frac{1}{2},\Pi_p)\right]\cdot \sum_{\substack{(m,Q)=1,\\ L^3\mid r\mid Q^{\infty}}} \chi(\frac{md}{c})\frac{A_G(d,m)}{\vert md\vert}\text{KL}_2(\overline{aQL^3}r,m;\frac{c}{d})\\
\cdot \left[\prod_{p\mid Q} W_{\tilde{G},p}\left(\left(\begin{matrix} p^{v_p(r/L^3)} & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1\end{matrix} \right)\right)\frac{W_{G,p}(1)}{W_{\tilde{G},p}(1)} \right] \cdot [\mathcal{H}_{\Pi_{\infty}}W]\left(\frac{rmd^2}{c^3QL^3X}\right). \nonumber
\end{multline}
The condition $L^3\mid r$ comes from the support of the local Whittaker functions. Indeed the denominator $L^3$ appears in a the first place due to a rough estimate of the $p$-adic Bessel-transforms. Our refined analysis of the case at hand has now revealed the true support. This allows us to shift the $r$-sum accordingly to get:
\begin{multline}
S_G(X) = \frac{cQ^{\frac{1}{2}}}{X} \cdot \left[\prod_{p\mid Q}\epsilon(\frac{1}{2},\Pi_p)\right]\cdot \sum_{\substack{(m,Q)=1,\\ r\mid Q^{\infty}}} \chi(\frac{md}{c})\frac{A_G(d,m)}{\vert md\vert}\text{KL}_2(\overline{aQ},rm;\frac{c}{d})\\
\cdot \left[\prod_{p\mid Q} W_{\tilde{G},p}\left(\left(\begin{matrix} p^{v_p(r)} & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1\end{matrix} \right)\right)\frac{W_{G,p}(1)}{W_{\tilde{G},p}(1)} \right] \cdot [\mathcal{H}_{\Pi_{\infty}}W]\left(\frac{rmd^2}{c^3QX}\right). \nonumber
\end{multline}
We now specialise to the case $\Pi=\Pi_F$. (Actually we only use that in this case $\Pi_p$ is self-dual it is easy to obtain a more general formula from the discussion above.) First we write $Q=N^2$, $\epsilon(\frac{1}{2},\Pi_{F,p})=1$ and $\chi=1$ in the formula above. Further, self-duality allows us to write
\begin{equation}
\frac{A_F(d,m)}{\vert md\vert }\left[\prod_{p\mid Q} W_{F,p}\left(\left(\begin{matrix} p^{v_p(r)} & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1\end{matrix} \right)\right) \right] = \frac{A_F(d,mr)}{\vert mdr\vert }.\nonumber
\end{equation}
Thus the $m$-sum and the $r$-sum can be combined. Thus, we have proved the following result.
$\mathbf{Theorem}$: Let $F$ be a newform for $\text{SL}_3(\mathbb{R})$ underlying the symmetric square lift of a newform $f$ of square free level $N$ and trivial character. Further fix a smooth function $W$ of compact support and take $(c,Na)=1$ as well as $X\in \mathbb{R}_{>0}$. Then we have
$$ \sum_{n\in\mathbb{Z}} \frac{A_F(n,1)}{\vert n\vert}e\left(n\frac{a}{c}\right)W\left(\frac{n}{X}\right) = \frac{cQ^{\frac{1}{2}}}{X} \cdot \sum_{m\in\mathbb{Z}_{\neq 0}} \frac{A_F(d,m)}{\vert md\vert }\text{KL}_2(\overline{aQ},m;\frac{c}{d})
\cdot [\mathcal{H}_{\Pi_{\infty}}W]\left(\frac{md^2}{c^3QX}\right).$$
Here the Bessel transform $[\mathcal{H}_{F}W]$ is defined above.
Let us mention that the formula provided by A. Corbett is very general and very flexible. Thus, along those lines it is possible to derive many more explicit Voronoi type formulae for $\text{GL}_n$ in particular one can generalise the approach presented above in several directions. However, in general the computational overhead might be far bigger.
Let us finish off by making some remarks concerning a follow up question by the OP in the comments. Looking at the $\text{GL}_2$-case one observes that very convenient Voronoi formulae for sums $$\sum_{n\in\mathbb{Z}} e(n\frac{a}{c})a_f(n)W(\frac{n}{X})$$ are available whenever $(c,N)=1$ or $N\mid c$. A for analytic purposes useful completely general formula was only recently found in ["Subconvexity for modular form $L$-functions in the $t$ aspect"](https://arxiv.org/abs/1707.01576). In particular the case $N=p^{2k}$ and $(N,c)=p^k$ is very hard. Returning to the $\text{GL}_3$-world we can say the following. The Voronoi formula for the sum $S_G(X)$ if $(c,Q)=1$ is discussed in this post (with a little more work one can derive an explicit formula for general $G$ dropping the self-duality assumption). The case $Q\mid c$ is covered in ["The Voronoi formula on $\text{GL}(3)$ with ramification"](https://arxiv.org/abs/1806.10786) (along the lines of this answer one can also derive the corresponding formula from A.Corbett's work). The case for general $c$ will be very difficult. It is covered by A. Corbett's work but the compute the local Bessel transforms promises to be very hard.
To illustrate this let us look at the following toy situation. Let $F$ be the symmetric-square lift of a newform $f$ of level $p$ and trivial character. Then $F$ has level $p^2$. We now take $c=c'p$ for $(c',p)=1$. (We are brief in details and translate everything straight away in the explicit situation at hand!) This hits exactly the hardest case as discussed above. However, our advantage is that we have full information on $\Pi_{F,p}$ as discussed above. Applying A. Corbett's theorem once again we find
\begin{equation}
S(X) = \frac{c'}{X} \sum_{\substack{(m,p)=1,\\ r\mid p^{\infty}}}\sum_{d\mid c'} \text{KL}_2(\overline{ap^4}r,m;\frac{c'}{d}) \frac{A_F(d,m)}{\vert md\vert} [\mathcal{H}_FW] \left(\frac{rmd^2}{Xp^5c'^3} \right)\cdot \mathcal{B}_{\Pi_{F,p},\phi_p}\left(\frac{rmd^2}{p^5c'^3}\right).\label{A}
\end{equation}
The job is now to evaluate the $p$-adic Bessel transform. For simplicity we assume that $W_{F,p}(1)=1$. In this case we have
\begin{multline}
\mathcal{B}_{\Pi_{p},\phi_{p}}(y)= \frac{\log(p)}{2\pi} \sum_{\substack{\xi\colon F^{\times}\to S^1,\\ \xi(p)=1}}\xi(y)\int_{\sigma-i\frac{\pi}{\log(p)}}^{\sigma+i\frac{\pi}{\log(p)}}\epsilon(\frac{1}{2},\xi\Pi_{p})p^{a(\xi \Pi_{p})(s-\frac{1}{2})}\frac{L(s,\xi^{-1}\tilde{\Pi}_{p})}{L(1-s,\chi\Pi_{p})}\vert y \vert_p^{1-s} \\ \cdot \int_{\mathbb{Z}_p\setminus \{0\}} \xi(x)\psi_p(x\frac{a}{c'p})\vert x\vert_p^{2-s}d^{\times}x ds. \nonumber
\end{multline}
The $x$-integral can be evaluated using Gau\ss\ sums. Indeed one gets
\begin{equation}
\int_{\mathbb{Z}_p\setminus \{0\}} \xi(x)\psi_p(x\frac{a}{c'p})\vert x\vert_p^{2-s}d^{\times}x = \begin{cases}
-\frac{p}{p-1}+L(1-s,\Pi_{F,p}) &\text{ if } \xi=1,\\
(1-p^{-1})^{-1}p^{-\frac{1}{2}}\xi(\frac{c'}{a})\epsilon(\frac{1}{2},\xi^{-1}) &\text{ if }a(\xi)=1,\\
0&\text{ else.}
\end{cases} \nonumber
\end{equation}
We conclude that in this case the Bessel transform reads
\begin{multline}
\mathcal{B}_{\Pi_{p},\phi_{p}}(y)=\frac{\log(p)}{2\pi} \int_{\sigma-i\frac{\pi}{\log(p)}}^{\sigma+i\frac{\pi}{\log(p)}} \bigg[p^{2s-1}L(s,\tilde{\Pi}_{F,p})-\frac{p^{2s}}{p-1}\frac{L(s,\tilde{\Pi}_{F,p})}{L(1-s,\Pi_{F,p})} \\ +\frac{p^{3s-1}}{p-1}\sum_{a(\xi)=1} \xi(-y\frac{c'}{a})\epsilon(\frac{1}{2},\xi)^2 \bigg] \vert y\vert^{1-s}ds. \label{B}
\end{multline}
With enough will-power one can certainly go further in the evaluation of this expression but we stop at this point. Our goal was to highlight some of the difficulties that arise when the additive twist 'colludes' with the level and we hopefully succeeded.