Nothing is wrong with the details of the computation, but you are treating the Euler-Maclaurian formula as an equality, when it isn't.
Before continuing, I will switch the sum a little bit: You are working with
$$S:=\sum_{k=0}^{\infty} e^{-k^2}.$$
The formulas are a little cleaner in terms of
$$T:=\sum_{k=-\infty}^{\infty} e^{-k^2}.$$
These are related by $T=2S-1$, so it is easy to switch between them.
The Euler-Macluarin formula with remainder term, combined with your correct computation that $f^{(k)}(x)$ goes to $0$ as $x \to \pm \infty$, tells us that we have
$$\sum_{k=-\infty}^{\infty} e^{-k^2} = \int_{x=-\infty}^{\infty} e^{-x^2} dx + \int_{-\infty}^{\infty} \frac{B_N(x-\lfloor x \rfloor)}{N!} \frac{d^N e^{-x^2}}{(dx)^N} dx$$
for any positive integer $N$. Here $B_N$ is the $N$-th Bernouli polynomial and $\lfloor x \rfloor$ is $x$ rounded down to an integer. For some functions, the remainder $\int \frac{B_N(x-\lfloor x \rfloor)}{N!} \frac{d^N f}{(dx)^N} dx$ goes to $0$ as $N \to \infty$, but it often doesn't, and it doesn't in your case.
To give some other examples, if you compute the asymptotics of $\sum_{k \leq M} \frac{1}{k}$ using Euler-Maclaurin, the Euler-Mascheroni constant occurs as the limit of the remainder term. If you derive Stirling's formula by Euler-Maclaurin summation of $\sum \log k$, then the $\log (2 \pi)$ occurs as the limit of the remainder term.
There is a really good discussion of this in Chapter 9 of Concrete Mathematics, by Graham, Knuth and Patashnik. You might particularly enjoy the fourth example in Chapter 9.6, where they use Euler Maclaurin summation to show that, for any $N$, we have
$$\sum_{k=-\infty}^{\infty} e^{-k^2/t} = \sqrt{\pi t} + O(t^{-N}) \ \mbox{as}\ t \to \infty.$$
Both the Euler-Maclaurin formula and Poisson summation are basic examples of deeper ideas. But at their core, the underlying ideas are very different.
The Euler-Maclaurin formula can be thought of as coming from repeated integration by parts, regarding the initial sum as a Rieman-Stieltjes integral
$$ \sum f(n) = \int_a^b f(t) d \lfloor t \rfloor, $$
where the constant terms in the antiderivatives appearing in the integrations by parts are chosen to minimize some cross error. This turns out to be the same view as estimating the expected error in an iterated composite trapezoidal rule for estimating an integral.
At its core, Euler-Maclaurin concerns estimating the difference between the discrete and continuous sums of a nice function $f(t)$. To that end, there exist generalizations to higher dimensions, and in particular to functions on polytopes. See this MO question for more on this connection, but the theme remains that one estimates a discretized sum (which is often quite hard to understand, but easier to actually compute) with an integral (which is often much easier to understand properties of, but can be hard to actually compute).
Many examples from elementary and analytic number theory of successes of Euler-Maclaurin formulas arise from computing high accuracy estimates of discrete sums of extremely well-behaved functions, such as computing
$$ \sum \log n, \quad \sum \frac{1}{x}.$$
In these cases, the smoothness properties of the function $f$ are more easily viewed from the integral.
Conversely, Euler-Maclaurin formulas are used to give high accuracy estimates of integrals in standard numerical analysis methods, as discrete sums with understandable error terms are actually computable and estimable. Note that when approximating an integral $\int_1^b f(x) dx$, the error comes in the form of high derivatives of $f$, so this still leverages smoothness properties of $f$, but in a different way.
Poisson summation is very different in nature. Poisson summation concerns Fourier series, or perhaps stated more fundamentally, the Fourier transform. The underlying big picture concerns the structure of $\mathbb{R}$ and $\mathbb{Z}$, and how $\mathbb{Z}$ behaves inside $\mathbb{R}$ as a discrete subgroup.
One can study Poisson summation for $\mathbb{Z}^n$ sitting within $\mathbb{R}^n$, or more general group quotients. Carrying this to one extreme leads to the Selberg Trace Formula (and in fact, one can directly view regular Poisson summation as a "trivial" trace formula. This point of view is apparent in the introductory chapter of Bump's Automorphic Forms and Representations, as a way of motivating eventual progress towards discussion trace formulas). From this point of view, it is clear that Poisson summation is fundamentally a consequence of representation-theoretic information.
Although over $\mathbb{R}/\mathbb{Z}$, this looks superficially similar to the integrals in Euler-Maclaurin summation, I would say that the underlying ideas are substantively different.
Best Answer
Finally got a solution. Posting it in case anyone's interested.
First, notice that because $ \int_0^{\infty} f(x) \, dx $ converges, $ x \to \infty \implies f^{(n)}(x) \to 0 $ for $ n = 0,1,\cdots, \infty $. Let us also remember the sine and cosine functions are bounded
Let's do the last integral by parts:
$ \begin{align} \int_{\mathbb{R}} f(x) e^{-2 \pi i x n} \, dx &= 2 \int_0^{\infty} f(x) \cos{(2 \pi x n)} \, dx \tag{because of parity} \\ &= 2 \left[ \left( \frac{f(x) \sin{(2 \pi n x)}}{2 \pi n} \right)_0^{\infty} - \int_0^{\infty} \frac{f'(x) \sin{(2 \pi n x)}}{2 \pi n} \, dx \right] \\ &= \cdots \\ &\approx 2 \left[ - \frac{f'(0)}{(2 \pi n)^2} + \frac{f'''(0)}{(2 \pi n)^4} \right] \end{align} $
The rest is easy.