Singular asymptotics of Gaussian integrals with periodic perturbations
Naive Approach: $$ \begin{align} \int_{-\infty}^\infty f(x)\,e^{-Ax^2}\,\mathrm{d}x &=\sum_{k=-\infty}^\infty\int_0^1f(x)\,e^{-A(x+k)^2}\,\mathrm{d}x\\ &=\int_0^1f(x)\sum_{k=-\infty}^\infty e^{-A(x+k)^2}\,\mathrm{d}x\tag{1} \end{align} $$ For $x\in[0,1]$ and $k\ge0$, $$ e^{-A(x+k+1)^2} \le\int_{x+k}^{x+k+1} e^{-At^2}\,\mathrm{d}t \le e^{-A(x+k)^2}\tag{2} $$ and $$ e^{-A(x-k-2)^2} \le\int_{x-k-2}^{x-k-1} e^{-At^2}\,\mathrm{d}t \le e^{-A(x-k-1)^2}\tag{3} $$ Summing $(2)$ and $(3)$ yields $$ -1+\int_{-\infty}^\infty e^{-At^2}\,\mathrm{d}t \le\sum_{k=-\infty}^\infty e^{-A(x+k)^2} \le2+\int_{-\infty}^\infty e^{-At^2}\,\mathrm{d}t\tag{4} $$ Since $\int_{-\infty}^\infty e^{-At^2}\,\mathrm{d}t=\sqrt{\frac\pi{A}}$, $(4)$ says $$ \sum_{k=-\infty}^\infty e^{-A(x+k)^2}=\sqrt{\frac\pi{A}}+O(1)\tag{5} $$ Combining $(1)$ and $(5)$ yields $$ \int_{-\infty}^\infty f(x)\,e^{-Ax^2}\,\mathrm{d}x=\left(\sqrt{\frac\pi{A}}+O(1)\right)\int_0^1f(x)\,\mathrm{d}x\tag{6} $$ For a given $f$, we can move the $O(1)$ outside to get $$ \int_{-\infty}^\infty f(x)\,e^{-Ax^2}\,\mathrm{d}x=\sqrt{\frac\pi{A}}\int_0^1f(x)\,\mathrm{d}x+O(1)\tag{7} $$
Fourier Series Approach:
Contour integration yields $$ \int_{-\infty}^\infty e^{2\pi inx}e^{-Ax^2}\,\mathrm{d}x =e^{-\pi^2n^2/A}\sqrt{\frac\pi A}\tag{8} $$ If $f=\sum\limits_na_ne^{2\pi inx}$, then $a_0=\int_0^1f(x)\,\mathrm{d}x$ $$ \int_{-\infty}^\infty f(x)e^{-Ax^2}\,\mathrm{d}x=\sqrt{\frac\pi{A}}\left(\int_0^1f(x)\,\mathrm{d}x+\sum_{n\ne0}a_ne^{-\pi^2n^2/A}\right)\tag{9} $$ and for $A\le1$, $$ \begin{align} \left|\sum_{n\ne0}a_ne^{-\pi^2n^2/A}\right| &\le\sup_{n\ne0}|a_n|\sum_{n\ne0}e^{-\pi^2n^2/A}\\ &\le\|f\|_{L^1[0,1]}2\sum_{n=1}^\infty e^{-\pi^2n/A}\\ &\le\|f\|_{L^1[0,1]}\frac{2e^{-\pi^2/A}}{1-e^{-\pi^2}}\tag{10} \end{align} $$ Combining $(9)$ and $(10)$ yields $$ \int_{-\infty}^\infty f(x)e^{-Ax^2}\,\mathrm{d}x =\sqrt{\frac{\pi}{A}}\int_0^1f(x)\,\mathrm{d}x+O\left(\frac1{\sqrt{A}}e^{-\pi^2/A}\right)\|f\|_{L^1[0,1]}\tag{11} $$
Convergence in the Fourier Series Approach:
Using the Fejér Kernel, we can show that trigonometric polynomials are dense in $L^1[0,1]$.
For trigonometric polynomials, there is no convergence problem in $(9)$ since all sums are finite. Likewise, $(11)$ also holds for trigonometric polynomials.
For any $\epsilon\gt0$, we can find a trigonometric polynomial $p$ so that $$ \|f-p\|_{L^1[0,1]}\le\epsilon\tag{12} $$ For $f-p$, both terms on the right side of $(11)$ are controlled by $(12)$. The only thing we need to control is the left hand side of $(11)$: $$ \begin{align} \int_{-\infty}^\infty|f(x)-p(x)|e^{-Ax^2}\,\mathrm{d}x &\le\|f-p\|_{L^1[0,1]}2\sum_{k=0}^\infty e^{-Ak^2}\\ &\le\|f-p\|_{L^1[0,1]}2\sum_{k=0}^\infty e^{-Ak}\\ &=\|f-p\|_{L^1[0,1]}\frac2{1-e^{-A}}\\ &\le\|f-p\|_{L^1[0,1]}\frac{2e}{e-1}\tag{13} \end{align} $$ Using $(12)$ and $(13)$ and the fact that $(11)$ holds for any trigonometric polynomial $p$, $(11)$ also holds for any $f\in L^1[0,1]$.
One way to look at this is as follows. Let $c = \int_0^1 f(x)\,dx$. Then $f(x) - c$ has integral $0$ over the period, and therefore is the derivative of a $C^1$ function $g(x)$ on ${\mathbb R}$ which also has period $1$. Then $$\int_{-\infty}^{\infty} f(x) e^{-Ax^2}\,dx = \int_{-\infty}^{\infty} c e^{-Ax^2}\,dx + \int_{-\infty}^{\infty} (f(x) - c) e^{-Ax^2}\,dx$$ The first term integrates to $c\sqrt{\frac{\pi}{A}} = \sqrt{\frac{\pi}{A}} \int_0^1 f(x)\,dx$, which is the main term. In the second term, if you integrate by parts you get $$\int_{-\infty}^{\infty} g(x)(2Ax e^{-Ax^2})\,dx$$ $g(x)$ is continuous and has period $1$, so $|g(x)| \leq M$ for some $M$. Thus we have $$\bigg|\int_{-\infty}^{\infty} g(x)(2Ax e^{-Ax^2})\,dx\bigg| \leq M \int_{-\infty}^{\infty} 2A|x| e^{-Ax^2}\,dx$$ $$=2M\int_{0}^{\infty} 2Ax e^{-Ax^2}\,dx$$ $$= 2M$$ This gives the error term.
Note that by the way $M$ is defined you have explicit bounds on $M$: For $0 \leq x \leq 1$ you have $$g(x) = \int_0^x f(y)\,dy - x\int_0^1f(y)\,dy = (1-x)\int_0^x f(y)\,dy - x\int_x^1 f(y)\,dy$$ So $$|g(x)| \leq (1-x)x\sup_{y \in \mathbb R}|f(y)| + x(1-x)\sup_{y \in \mathbb R}|f(y)|$$ $$=2x(1-x)\sup_{y \in \mathbb R}|f(y)|$$ Since the supremum of $x(1-x)$ on $[0,1]$ is ${1 \over 4}$, we get $$M = \sup_{x \in [0,1]}|g(x)| \leq {1 \over 2}\sup_{y \in \mathbb R}|f(y)|$$
I will also do it formally though maybe with a less number of steps to justify. It is clear that since $f$ is periodic, \begin{align} \int_{-\infty}^{\infty}f(x)e^{-Ax^2}dx&=\int_0^1 f(x)\sum_{n\in\mathbb{Z}}e^{-A(x+n)^2}dx=\\&= \int_0^1 f(x)e^{-A x^2}\vartheta_3\Bigl(iAx\Bigl|\Bigr.\frac{iA}{\pi}\Bigr)dx,\tag{1} \end{align} where $\vartheta_3(z|\tau)=\sum_{n\in\mathbb{Z}}e^{i\pi\tau n^2+2 i n z}$ denotes Jacobi theta function. Full asymptotic expansion of (1) in most easily obtained using
Jacobi imaginary transformation $$\vartheta_3(z|\tau)=(-i\tau)^{-\frac12}\exp\left\{-\frac{iz^2}{\pi\tau}\right\}\,\vartheta_3\Bigl(-\frac{z}{\tau}\Bigl|\Bigr.-\frac{1}{\tau}\Bigr)\tag{2}$$ after which (1) becomes $$ \sqrt{\frac{\pi}{A}}\int_0^1f(x)\,\vartheta_3\Bigl(\pi x\Bigl|\Bigr.\frac{i\pi}{A}\Bigr)dx.\tag{3}$$ Here the prefactor comes from $(-i\tau)^{-\frac12}$, the exponential in the integral is cancelled by the second factor in (2) and we also used that $\vartheta_3(z|\tau)$ is an even function of $z$. Let me stress that so far no approximations were made, the only assumption is that the summation can be interchanged with the integration at the very first step.
Note that the first argument of the theta function in (3) does not depend on $A$ and, especially, that the half-period ratio $\tau$ now tends to $i\infty$ instead of $i0^+$. This allows to extract the asymptotic expansion from the triple product formula $$\vartheta_3(z|\tau)=\prod_{n=1}^{\infty}(1-q^{2n})\prod_{n=1}^{\infty}(1+e^{2iz}q^{n-\frac12})\prod_{n=1}^{\infty}(1+e^{-2iz}q^{n-\frac12}),\tag{4}$$ where $q=e^{i\pi\tau}$. The point is that as $A\rightarrow 0^+$, the nome $q=e^{-\frac{\pi^2}{A}}$ in (3) rapidly goes to $0$, so that (3) becomes $$\sqrt{\frac{\pi}{A}}\int_0^1f(x)\,\Bigl[1+2\cos2\pi x\,e^{-\frac{\pi^2}{A}}+O(e^{-\frac{2\pi^2}{A}})\Bigr]dx.$$ It is clear that, in principle, one can derive from (4) any desired number of terms in the asymptotic expansion.
Added: It may seem that the above procedure works only because of the specific Gaussian form of the kernel. This is not so: if we had some other function $g(x)$ instead of $e^{-Ax^2}$, Jacobi imaginary transformation used above would be replaced by Poisson summation formula.
Your estimate of error term is correct. The following are just some supplementary details to make your argument more rigorous.
Let $(f_n)_{n\ge 1}$ be the sequence of partial sums of the Fourier series of $f$, i.e. $$ f_n(x) = \int_0^1 f(t)\,dt + \sum_{k=1}^n \big(a_k \cos(2\pi kx) + b_k \sin(2\pi kx)\big). $$ Note that $f_n$ converges to $f$ in $L^2([0,1])$, so by Cauchy-Schwarz's inequality, as $n\to\infty$, $$\int_0^1|f(t)-f_n(t)| dt\le \big(\int_0^1|f(t)-f_n(t)|^2dt\big)^{\frac{1}{2}}\to 0.\tag{1}$$ Also note that, for every $n\ge 1$,
$$|a_n|=2\cdot\big|\int_0^1 f(t)\cos(2\pi nt) dt\big|\le 2\int_0^1|f(t)|dt.\tag{2}$$
As you have shown, $$\int_{-\infty}^{+\infty}f_n(x)e^{-Ax^2}dx=\sqrt{\frac{\pi}{A}}(\int_0^1f(t)dt+\sum_{k=1}^na_ke^{\frac{-k^2\pi^2}{A}}).\tag{3}$$ From $(2)$ and $(3)$ we know, $$\big|\sqrt{\frac{A}{\pi}}\int_{-\infty}^{+\infty}f_n(x)e^{-Ax^2}dx-\int_0^1f(t)dt\big|\le 2\int_0^1|f(t)|dt\cdot\sum_{k=1}^\infty e^{\frac{-k^2\pi^2}{A}}\le M e^{\frac{-\pi^2}{A}},\tag{4}$$ where $M>0$ is independent of $n\ge 1$ and $0<A\le 1$. Since for every $m\in\mathbb{Z}$, $$\int_m^{m+1}|f_n(x)-f(x)|dx=\int_0^1|f_n(x)-f(x)|dx,$$
$$ \int_{-\infty}^{+\infty}|f_n(x)-f(x)|e^{-Ax^2}dx\le2 \int_0^1|f_n(x)-f(x)|dx\cdot\sum_{m=0}^\infty e^{-Am^2}<\infty.\tag{5} $$ Letting $n\to\infty$ in $(5)$, from $(1)$ we know that
$$\lim_{n\to\infty}\int_{-\infty}^{+\infty}|f_n(x)-f(x)|e^{-Ax^2}dx=0.\tag{6}$$ Combining $(4)$ and $(6)$, it follows that
$$|\int_{-\infty}^{+\infty}f(x)e^{-Ax^2}dx-\sqrt{\frac{\pi}{A}}\int_0^1f(t)dt|\le M\sqrt{\frac{\pi}{A}} e^{\frac{-\pi^2}{A}}.\tag{7}$$