Summation of $\sum_{n=0}^{\infty}a^nq^{n^2}$

A possible approach could be to search an approximated estimate of the summation. Applying the Euler-Maclaurin formula to $f(n)=a^nq^{n^2}$ in the interval between $0$ and $N$, we have

$$\displaystyle \sum_0^N a^nq^{n^2}=\int_0^N a^nq^{n^2}+ \frac{a^Nq^{N^2}+1}{2} \\ + \sum_{k=1}^{\infty }\,{\frac {B_{2k}}{(2k)!}}\left(f^{(2k-1)}(N)-f^{(2k-1)}(0)\right)$$

where $B_{2k}$ denotes the $2k^{th}$ Bernoulli number and $f^{(2k-1)}$ is the $(2k-1)^{th}$ derivative.

To estimate the integral, we can make the substitution $$n=\frac{t}{\sqrt{\log q}}-\frac{\log a}{2 \log q}$$ We have

$$\displaystyle a^nq^{n^2}=a^{t/\sqrt{\log q} -\log a/(2\log q) } q^ {[t/\sqrt{\log q} -\log a/(2\log q) ]^2} \\ = e^{-\log^2 a/(4 \log q)}e^{t^2} $$

Also, we have

$$t=\frac{\log a +2n \log q}{2\sqrt{\log q}}$$ $$dt=\sqrt{\log q } \, dn$$ $$dn=\frac{1}{\sqrt{\log q}} dt$$

Then our integral can be rewritten as

$$\displaystyle \int_0^N a^nq^{n^2}dn\\= \frac{e^{-\log^2 a/(4 \log q)}}{\sqrt{\log q}} \int_{\frac{\log a}{2\sqrt{\log q}}}^{\frac{\log a +2N\log q}{2\sqrt{\log q}}} e^{t^2}$$

Setting by simplicity $$K_1=\log a/(2\sqrt{\log q})$$ and $$K_2=(\log a +2N\log q)/(2\sqrt{\log q})$$ we can further rewrite the integral as

$$\displaystyle \int_0^N a^nq^{n^2}dn =\frac{e^{-K_1^2 }}{\sqrt{\log q}} \int_{K_1}^{K_2} q^{t^2}\\ =\frac{e^{-K_1^2}}{\sqrt{\log q}} \left(\int_{0}^{K_2} q^{t^2} -\int_{0}^{K_1} q^{t^2} \right)\\ = \frac{\sqrt{\pi}\,e^{-K_1^2 }}{2\sqrt{\log q}} \bigg[\operatorname{erfi}(K_2) - \operatorname{erfi}(K_1)\bigg]\\$$

where $\operatorname{erfi}$ indicates the imaginary error function, defined as

$$\displaystyle \operatorname {erfi} (x)={\frac {2}{\sqrt {\pi }}}\int_{0}^{x}e^{t^{2}}\,dt$$

Alternatively, our integral could be written in terms of the Dawson function. In fact

$$\displaystyle \int_0^N a^nq^{n^2}dn\\ =\frac{e^{-\frac{\log^2 a}{4\log q} }}{\sqrt{\log q}} \left(\int_{0}^{K_2} q^{t^2} -\int_{0}^{K_1} q^{t^2} \right)\\ = \frac{e^{-K_2^2 +n \log a +2n^2\log q }}{\sqrt{\log q}} \int_{0}^{K_2} q^{t^2} \\ - \frac{e^{-K_1^2 }}{\sqrt{\log q}} \int_{0}^{K_1} q^{t^2} \\ =\frac{a^Nq^{N^2}}{\sqrt{\log q}} \, e^{-K_2^2} \int_{0}^{K_2} e^{t^2}\\ - \frac{e^{-K_1^1}}{ \sqrt{\log q}} \int_{0}^{K_1} e^{t^2} \\ =\frac{{a^Nq^{N^2} \,F_+(K_2) }- {F_+(K_1)}}{\sqrt{\log q}} $$

where $F_+(K_2)$ and $F_+(K_1)$ indicate the Dawson function, i.e. the one-sided Fourier–Laplace sine transform of the Gaussian function, usually expressed in the following two forms:

$${\displaystyle F_{+}(x)=e^{-x^{2}}\int _{0}^{x}e^{t^{2}}\,dt}$$

$${\displaystyle F_{-}(x)=e^{x^{2}}\int _{0}^{x}e^{-t^{2}}\,dt.\!}$$

Note that the Dawson function is linked to the error function and the imaginary error function by

$${\displaystyle \operatorname {erf} (x)={\frac {2}{\sqrt {\pi }}e^{-x^{2}}F_-(x)}}$$

$${\displaystyle \operatorname {erfi} (x)={\frac {2}{\sqrt {\pi }}e^{-x^{2}}F_+(x)}}$$

---

To proceed further, we have to note that the initial sum $\sum_0^\infty a^n q^{n^2}$ converges for $q<1$ and diverges for $q>1$ (the case $q=1$ reduces to the trivial sum $\sum_0^\infty a^n$, whose convergence or divergence depends on $a$). These two cases determine the presence of real or imaginary terms in the last equations. Since the OP asks for the result of infinite sums, I will focus only on the first case. However, the same following approach could be used to estimate partial sums in the second case.

If $q<1$, then $K_2$, $K_1$, and the term $\sqrt{ \log q}$ in the denominators of the RHS of the equations above are all imaginary. So if we write

$$K_2=\frac{\log a +2N\log q}{2\sqrt{\log q}}= \frac{\log a +2N\log q}{2\sqrt{-\log 1/q}}= -i\frac{\log a +2N\log q}{2\sqrt{\log 1/q}} =i\hat{K_2} $$

and similarly

$$K_1=\frac{\log a }{2\sqrt{\log q}}= \frac{\log a }{2\sqrt{-\log 1/q}}= -i\frac{\log a }{2\sqrt{\log 1/q}} =i\hat{K_1} $$

reminding that $\operatorname{erfi}(z i)= i\,\operatorname{erf}(z)$, we have

$$\displaystyle \int_0^N a^nq^{n^2}dn\\ = \frac{\sqrt{\pi} \, e^{\hat{K_1}^2} }{2\sqrt{\log q}}\left[ { \operatorname{erfi}(i\hat{K_2}) }- { \operatorname{erfi}(i\hat{K_1}) }\right] \\ = \frac{\sqrt{\pi} \, e^{\hat{K_1}^2} } {2\sqrt{\log (1/q)}}\left[ { \operatorname{erf}(\hat{K_2}) }- { \operatorname{erf}(\hat{K_1}) }\right] $$

Now we can consider that, as $N$ increases, the error function shows a fast convergence to $1$, satisfying $\operatorname{erf}(z)=1+O(e^{-z^2}z^{-1})$. Therefore, making the inverse substitution $\hat{K_1}= - (\log a )/(2\sqrt{\log 1/q}) $, we get

$$\displaystyle \int_0^\infty a^nq^{n^2}dn\\ = \frac{\sqrt{\pi} \, e^{ \frac{\log^2 a}{4\log 1/q}} } {2\sqrt{\log (1/q)}}\left[ 1- { \operatorname{erf}\left(-\frac{ \log a }{2\sqrt{\log 1/q}}\right) }\right] $$

---

Turning back to the initial Euler-Maclaurin expansion, the degree of approximation for the sum $\sum_0^\infty a^nq^{n^2}$ depends on the number of terms considered after the integral. Among these, it is not difficult to note that the leading terms result from the values of the function and the derivatives in $n=0$. For most values of $a$ and $q$, a good level of approximation can be obtained considering these terms up to the third derivative. Because $$f^{(1)}(0)=\log a$$ $$f^{(3)}(0)=\log^2 a + 6\log a \log q$$

and $$B_2=\frac{1}{6}$$ $$B_4=-\frac{1}{30}$$

we have

$$\displaystyle \sum_0^N a^nq^{n^2}\approx \frac{\sqrt{\pi} \, e^{ \frac{\log^2 a}{4\log 1/q}} } {2\sqrt{\log (1/q)}}\left[ 1- { \operatorname{erf}\left(-\frac{ \log a }{2\sqrt{\log 1/q}}\right) }\right] + \frac{1}{2} -\frac{\log a}{12}+\frac{\log a(\log^2 a +6 \log q)}{720} $$

For example, setting $a=1/2$ and $q=3/4$, we have

$$\sum_0^\infty a^nq^{n^2}=1.46413758...$$

as shown by WA here. Our approximated formula yields

$$\displaystyle \sum_0^N a^nq^{n^2}\approx \frac{\sqrt{\pi} \, e^{ \frac{\log^2 1/2}{4\log 4/3}} } {2\sqrt{\log (4/3)}}\left[ 1- { \operatorname{erf}\left(-\frac{ \log 1/2 }{2\sqrt{\log 4/3}}\right) }\right] + \frac{1}{2} -\frac{\log 1/2}{12}+\frac{\log 1/2 (\log^2 1/2 +6 \log 4/3)}{720} \\ =1.46407819...$$

as confirmed by WA here. The relative error, expressed in percentage, is $-0.004\%$.

Including further terms of the Euler-Maclaurin espansion, the approximation can be improved. For example, considering the terms up to the fifth derivative, taking into account that

$$f^{(5)}(0)\\=\log^5 a + 20 \log^3 a \log q \\ + 60 \log a \log^2 q $$

and $$B_6=-\frac{1}{42}$$

the approximated formula above would include the additional term

$$-\frac{1}{30240} \bigg[\log^5 (1/2) + 20 \log^3 (1/2) \log (3/4) + 60 \log(1/2) \log^2 (3/4) \bigg]=0.00005574... $$

as shown here. In this way, the estimated sum would become

$$=1.46407819 + 0.00005574 = 1.46413393 $$

with a relative error of $-0.00025\%$.

Lastly, to facilitate the computation of the approximated formula (e.g. in Python), the value of the $\operatorname{erf}$ function can be calculated using several numerical approximations. These allow to calculate the value of $$\operatorname{erf}\left(-\frac{ \log a }{2\sqrt{\log 1/q}}\right) $$ with good precision using elementary functions.