which integers take the form $x^2 + xy + y^2$?

Let $r_2(n)$ be the number of representations of $n$ as a sum of two squares, and let $l(n)$ be the number of ways to write $n$ as $x^2+xy+y^2$. Then as you mentioned we have that $$\sum_{y,x\in \mathbb{Z}} q^{x^2+y^2}=\sum_{n=0}^\infty r_2(n) q^n \ \text{and} \ \sum_{y,x\in \mathbb{Z}} q^{x^2+xy+y^2}=\sum_{n=0}^\infty l(n) q^n.$$

Exact number of representations: We can explicitely write down $r_2(n)$, the number of representations of $n$ as a sum of two squares. See Sum of Squares Function.

Also, in a similar manner we can explicitely write down $l(n)$, the number of representations of $n$ as $x^2+xy+y^2$. See the answer on this math stack exchange post.

Theta Functions: We can evaluate the infinite series in terms of Jacobi theta functions without much difficulty. To deal with the sum of squares, notice that

$$\sum_{y,x\in \mathbb{Z}} q^{x^2+y^2} =\left( \sum_{n=-\infty}^\infty q^{n^2}\right)^2= \vartheta_3(q)^2.$$ Next, we can transform $x^2 +xy+y^2$ into $\frac{n^2+3m^2}{4}$ where $m,n$ must have the same parity. Then $$\sum_{x,y}q^{x^{2}+xy+y^{2}}=\sum_{\begin{array}{c} x,y\\ both\ odd \end{array}}q^{\frac{x^{2}+3y^{2}}{4}}+\sum_{\begin{array}{c} x,y\\ both\ even \end{array}}q^{\frac{x^{2}+3y^{2}}{4}} $$

$$=q\sum_{n,m\in\mathbb{Z}}q^{n(n+1)+3m(m+1)}+\sum_{n,m\in\mathbb{Z}}q^{n^{2}+3m^{2}}.$$This can be rewritten in terms of the jacobi theta functions as $$\vartheta_2(q)\vartheta_2(q^3)+\vartheta_3(q)\vartheta_3(q^3).$$

Hope that helps,

Just to point out that $x^2 + x y + y^2$ represents exactly the same numbers as $x^2 + 3 y^2,$ just with different frequency. The similar fact for indefinite forms is that $x^2 + x y - y^2$ represents exactly the same numbers as $x^2 - 5 y^2.$ In both cases $x^2 + xy \pm y^2$ the sum is not even unless $x,y$ are both even, which says that the order with which the prime 2 divides the result must be even. For odd primes, the analogous statement for primes dividing the result can be found using the Legendre symbol, in your case $(-3 | p) = (p | 3).$ It is easy to prove that $x^2 + x y + k y^2$ represents a superset of $x^2 + (4k-1) y^2,$ for almost all $k$ a strict superset.

In his recent book, Number Theory in the Spirit of Ramanujan, Bruce Berndt discusses representation problem for (1, 0, 1), (1, 0, 2), (1, 1, 1), (1, 0, 3).

This is a quote from the arXiv version of Alexander Berkovich and Hamza Yesilyurt which has:

Comments: 26 pages, no figures, fun to read

On page 149 of Rational Quadratic Forms by Cassels, Lemma 6.3 gives a count for the number of primitive representations of a number $n = x^2 + y^2,$ and this lemma can be applied to find the total number of representations, even if there are none primitive. Something very similar should work for $x^2 + x y + y^2,$ perhaps it is in the Berndt book.

As far as "closed form," the best we can expect is in the world of modular forms and, in particular, "eta quotients." You might just look at all Alex Berkovich's manuscripts on the arXiv, he plays with these infinite sums all day.

EDIT, 24 January 2012: Exercise 2 on page 80 of Introduction to the Theory of Numbers by Leonard Eugene Dickson: The number of representations of positive integer $n$ by $x^2 + x y + y^2$ is $6 E(n),$ where $E(n)$ is the excess of the number of divisors that are $1 \pmod 3$ over the number of divisors that are $2 \pmod 3.$

One can start showing the following:

The integers $n$ which are of the form $x^2+xy+y^2$, for two relatively prime integers $x,y$ are precisely those positive integers occurring as divisors of $m^2+m+1$, for some integer $m$.

In other words, the polynomial $f(x)=x^2+x+1$ has the property that the positive divisors of the integers it represents are precisely those integers that can be properly represented by its homogenization $F(x,y)=x^2+xy+y^2$ ("properly" here refers to the condition $(x,y)=1$).

The proof uses the fact that the imaginary quadratic order $\mathbf{Z}[x]/f(x)$ has class number one, as anticipated by Elkies. It goes as follows:

Let $n$ be a positive divisor of $m^2+m+1$, for some integer $m$. Consider the quadratic form in $x,y$ given by:

$Q(x,y)=\frac{m^2+m+1}{n}x^2-(2m+1)xy+ny^2$;

it has integer coefficients, positive definite, and has discriminant equal to $-3$ (in particular it is primitive).

Since $h(-3)=1$, there is only one reduced, positive definite quadratic form of discriminant $-3$. This is $E(x,y)=x^2+xy+y^2$. Therefore $Q$ and $E$ are properly equivalent (that is there is a determinant-one change of variables taking one into the other), and since $Q$ certainly properly represents $n$, so does $E$.

The converse is similar and uses the fact that if a quadratic form $Q$ properly represents an integer $n$, then $Q$ is properly equivalent to a form of the type $nx^2+bxy+cy^2$ (this is lemma 2.3 of Cox's wonderful book "Primes of the form x^2+ny^2").

Once you understand the positive integers $n$ that are properly represented by $E$, then you can get them all, after scaling by squares.

The original problem is then reduced to understanding those integers $n$ for which $x^2+x+1$ has a zero in $\mathbf{Z}/(n)$. Using the Chinese Remainder Theorem, this can be reduced to the case where $n$ is a prime power $p^s$. Then for $p\neq 3$ Hensel's Lemma tells you that your equation has a solution mod $p^s$ if and only if it has a solution mod $p$. With quadratic reciprocity you can conclude that any prime divisor of $n$ has to be congruent to $1$ mod $3$. I am at the moment missing how you can solve the equation $x^2+x+1=0$ in $\mathbf{Z}/3^s$, but I think that a version of Hensel's Lemma applies. I'll think about it.

[EDIT: The equation $x^2+x+1=0$ has a solution in $\mathbf{Z}/(3^s)$, with $s\geq 1$, if and only if $s=1$. (Therefore $x^2+xy+y^2$ represents properly only $3$ and $1$ as powers of $3$.) One can see this by checking that there is no solution for $s=2$, and therefore for $s>2$. More generally, if $p>2$ then $x^{p-1}+x^{p-2}+\ldots+x+1=0$ has no solutions in $\mathbf{Z}/(p^s)$ with $s>1$. The $p$--adic valuation of an integer of the form $x^{p-1}+x^{p-2}+\ldots+x+1=(x^p-1)/(x-1)$ is either zero or one.]

(Suggested reading. Cox's book quoted above and Serre's paper: $\Delta=b^2-4ac$)