Start with $f\in \mathcal{S}(\mathbb{R}^n)$. For every $L>0$ let $\chi_1\in C_c^\infty(\mathbb{R})$ be the standard cut-off function such that $\varphi(t) =0$ for $\vert t \vert \geq 1$, $\varphi(0)=1$ and $\varphi^{(k)}(1)=0$ for all $k\geq 1$. Then we define for every $L\geq 0$ the cut-off function
$$ \varphi_L: \mathbb{R} \rightarrow \mathbb{R}, \varphi_L(t) = \begin{cases} 0,& \vert t \vert\geq L+1, \\ 1,& \vert t \vert \leq L,\\ \varphi(t-\mathrm{sgn(t)L}),& L<\vert t \vert<L+1. \end{cases} $$
One readily checks that $\varphi_L\in C_c^\infty(\mathbb{R})$. We now define the cut-off function
$$ \chi_L:\mathbb{R}^n \rightarrow \mathbb{R}, \chi_L(x) = \prod_{j=1}^n \varphi_L(x_j). $$
It is easy to check that for every $\alpha \in \mathbb{N}^n$ and all $L\geq 1$ we have
$$ \sup_{x\in \mathbb{R}^n} \vert \partial^\alpha \chi_L(x) \vert = \sup_{x\in \mathbb{R}^n} \vert\partial^\alpha \chi_1(x) \vert < \infty. $$
The problem boils down to find a positive $g\in \mathcal{S}(\mathbb{R}^n)$ such that $\sqrt{g}\in \mathcal{S}(\mathbb{R}^n)$. Once we have found such a function, then we can define
$$ f_\varepsilon(x):=\sqrt{\chi_{1/\varepsilon}(x)f(x)+ \varepsilon g(x)}. $$
Define the semi-norms
$$ \Vert f \Vert_{\alpha, \beta} := \sup_{x \in \mathbb{R}^n} \vert x^\beta \partial^\alpha f(x) \vert. $$
Since
\begin{align*}
\sup_{x\in \mathbb{R}^n} \vert \partial^\alpha \chi_{1/\epsilon}(x) \vert <\infty
\end{align*}
as noted above, we have
\begin{align*}
\lVert f_\epsilon^2 \rVert_{\alpha. \beta} <\infty
\end{align*}
for all $\alpha, \beta$, which implies that $f_\epsilon^2 \in \mathcal{S}(\mathbb{R}^n)$.
Now, consider another Schwartz function
\begin{align*}
F_\epsilon:=(1-\chi_{1/\epsilon}) f
\end{align*}
Since $F_\epsilon$ is supported on $[\frac{1}{\epsilon}, \infty)^n \sqcup (-\infty, -\frac{1}{\epsilon}]^n$, the estimate
\begin{align*}
\frac{1}{\epsilon^n} \lvert x^{\beta}\partial^\alpha F_\epsilon(x) \rvert \leq \Bigl \lvert \Bigl(\prod_{i=1}^n x_i \Bigr) x^{\beta}\partial^\alpha F_\epsilon(x) \Bigr \rvert=\lvert x^{\beta+(1,\cdots,1)}\partial^\alpha F_\epsilon(x) \rvert
\end{align*}
holds for $x \in \mathbb{R}^n$ so that
\begin{align*}
\frac{1}{\epsilon^n} \Vert F_\epsilon \Vert_{\alpha, \beta} \leq \Vert F_\epsilon \Vert_{\alpha, \beta+(1, \dots, 1)}
\end{align*}
We also observed above that the supremum of $\partial^\alpha \chi_L$ does not depend on $L$. Hence $\Vert F_\epsilon \Vert_{\alpha, \beta+(1, \dots, 1)}$ does not depend on $\epsilon$, implying
\begin{align*}
\Vert f_\varepsilon^2 - f \Vert_{\alpha, \beta}
\leq \Vert F_\epsilon \Vert_{\alpha, \beta} + \varepsilon \Vert g \Vert_{\alpha, \beta}
\leq \epsilon^n \Vert F_\epsilon \Vert_{\alpha, \beta+(1, \dots, 1)} + \varepsilon \Vert g \Vert_{\alpha, \beta} \to 0^+ \text{ as } \epsilon \to 0^+
\end{align*}
for all $\alpha,\beta$. That is, $f_\varepsilon^2 \to f$ in the Fréchet topology of $\mathcal{S}(\mathbb{R}^n)$.
Lastly, we need to check if $f_\epsilon \in \mathcal{S}(\mathbb{R}^n)$. However, because $g$ is assumed to be positive, it is clear that $f_\epsilon$ is smooth on whole $\mathbb{R}^n$. Moreover, $f_\varepsilon (x)=\sqrt{\varepsilon g(x)}$ for $\vert x \vert > 1/\varepsilon$ so that
\begin{equation}
\lvert x^\beta \partial^\alpha f_\epsilon(x) \rvert \leq \sup_{ \lvert x \rvert \leq \frac{1}{\epsilon} } \lvert x^\beta \partial^\alpha f_\epsilon(x) \rvert + \sqrt{\epsilon}\sup_{ \lvert x \rvert > \frac{1}{\epsilon} } \lvert x^\beta \partial^\alpha \sqrt{g}(x) \rvert
\end{equation}
for all $x \in \mathbb{R}^n$. Therefore,
\begin{equation}
\lVert f_\epsilon \rVert_{\alpha, \beta} < \infty
\end{equation}
for all $\alpha, \beta$, implying $f_\epsilon \in \mathcal{S}(\mathbb{R}^n)$.
Note that $g(x)=e^{-x^2}$ does the job as $\sqrt{g(x)}= e^{-x^2/2}$ is again Schwartz.
Best Answer
$F$ is not necessarily Schwartz. If it was, then its Fourier transform will also be Schwartz. Now since $F = f * |\cdot|^{-1}$, and the Fourier transform of $|\cdot|^{-1}$ as a tempered distribution is $c_3 |\xi|^{-2}$ (where $c_3$ is a dimensional constant), we see that $\hat{F} = c_3\hat f |\cdot|^{-2}$ as convolution corresponds to multiplication on the Fourier side. We could design $\hat f $ to be a smooth compactly supported bump function equal to $1$ on $|\xi| \le 1$ and consequently, $$ \sup_{\xi \in \mathbb{R^3}} c_3|\xi|^{-2} |\hat f(\xi)| = +\infty, $$ i.e. $\hat F$ is not Schwartz.