Find an element with given periodicity Sorry for all the confusion. I think what I am actually asking is: Can we find an explicit smooth non-zero function on $\mathbb R^2$ that satisfies
$$f(x_1,x_2) =e^{-i\pi x_2} f(x_1+1,x_2) \text{ and } f(x_1,x_2) =e^{i\pi x_1} f(x_1,x_2+1).$$
 A: Your space can be considered as sections of a complex line bundle over the torus. Note that the usual partial derivatives $\partial_1,\partial_2$ do not preserve it, but the operators
$$
D_1 = \partial_1 - i\pi x_2,D_2 = \partial_2 + i\pi x_1
$$
do — they define a (unitary) connection on your line bundle, which essentially means that $D_i(fg) = (\partial_i f)g + f(D_i g)$ for $f$ periodic and $g$ in your subspace. Now we have
$$
[D_1,D_2] = 2i\pi
$$
so these "connection partial derivatives" no longer commute. This calculation immediately shows that there can be no smooth function $g$ which is nowhere vanishing in your subspace; if that were the case, we could set $D_ig = A_ig$ for some periodic functions $A_i$, which yields
\begin{align*}
(2\pi i)g &= D_1D_2 g - D_2D_1 g\\
&= D_1(A_2 g) - D_2(A_1 g)\\
&= (\partial_1 A_2 + A_1A_2 - \partial_2 A_1 -A_1A_2)g\\
\implies 2\pi i &= \partial_1 A_2 - \partial_2 A_1
\end{align*}
But this is impossible since the integral of a partial derivative of a periodic function vanishes. (This proof essentially uses that $[D_1,D_2]\mathrm dx_1\mathrm dx_2$, the curvature of the connection, represents a nontrivial de Rham cohomology class on the torus, so that the line bundle is nontrivial, while such a function $g$ would determine a trivialization.)
The main advantage of the Fourier basis is that it diagonalizes the partial derivative operators. This doesn't work for the $D_i$ since they satisfy the canonical commutation relations (their commutator is a non-zero multiple of the identity); this is essentially the Heisenberg uncertainty relation. However, the CCR have a unique (irreducible) representation, and from this we obtain the existence of an isometry from your space to $L^2(\mathbb R)$ which carries $D_1$ to multiplication by $2\pi i x$ and $D_2$ to $-\partial_x$.
In fact, one can write this isomorphism down quite explicitly: the operator $-D_1^2 - D_2^2 = (D_1 -iD_2)^\dagger(D_1 -iD_2) +2\pi$ has smallest eigenvalue $2\pi$ with one-dimensional eigenspace spanned by a function $f_0$ such that $D_1 f_0 = iD_2 f_0$. One can then act by the "creation operator" $(D_1 -iD_2)$; setting $f_n = (D_1 -iD_2)^n f_0$, these functions form an orthogonal basis in which the action of the $D_i$ is tridiagonal (the resulting matrix vanishes except for the two diagonals around the main diagonal).
It remains to solve the equation $D_1 f = iD_2 f$ with your boundary conditions. Without boundary conditions, the general solution on $\mathbb R^2\cong \mathbb C$ is of the form $e^{-\pi\lvert z\rvert^2/2}g(z)$ with $g$ holomorphic. The boundary conditions become
\begin{align*}
g(z+1) & {}= e^{-\pi z +\pi/2}g(z), \\
g(z+i) & {}= e^{i\pi z +\pi/2}g(z).
\end{align*}
This is the functional equation of the Weierstrass sigma function of the square lattice.
Of course, this basis is probably not that useful for numerics. I would try to discretize the annihilation and creation operators, find the kernel of the former and act on it by the latter; in the resulting basis the operators $D_i$ should have small off-diagonal terms.
A: Fix a Schwartz function $g_0 : {\bf R} \to {\bf C}$,
say $g_0(x) = e^{-\pi x^2}$ (which looks like a natural choice in this context).
For $x_1,x_2 \in \bf R$ define $g(x_1,x_2)$ by the absolutely convergent sum
$$
g(x_1,x_2) = \sum_{n=-\infty}^\infty g_0(x_2+n) e^{2\pi i n x_1}.
$$
Clearly $g(x_1,x_2) = g(x_1+1,x_2)$ for all $x_1,x_2$.  Also
$$
g(x_1,x_2) = \sum_{n=-\infty}^\infty g_0(x_2+n+1) e^{2\pi i {n+1} x_1}
= e^{2\pi i x_1} \sum_{n=-\infty}^\infty g_0(x_2+n+1) e^{2\pi i n x_1}
= e^{2\pi i x_1} g(x_1, x_2 + 1).
$$
Finally define
$$
f(x_1,x_2) = e^{\pi i x_1 x_2} g(x_1,x_2).
$$
Then $f$ is a smooth function satisfying the required identities
$$
f(x_1,x_2) = e^{-i\pi x_2} f(x_1+1,x_2), \quad 
f(x_1,x_2) = e^{i\pi x_1} f(x_1,x_2+1)
\tag{$\star$}
$$
for all real $x_1,x_2$.
This construction was surmised by working backwards, observing that
if $f$ satisfies the required quasiperiodicity then $g$ is periodic in $x_1$
and thus has a Fourier series
$\sum_{n=-\infty}^\infty g_n(x_2) e^{2\pi i n x_1}$,
and then the identity $g(x_1,x_2) = e^{2\pi i x_1} g(x_1, x_2 + 1)$
yields $g_n(x_2+1) = g_{n+1}(x_2)$, whence $g_n(x_2) = g_0(x_2+n)$ for each $n$.
The $L^2$ norm of $g$ on the unit square is the square root of
$$
\sum_{n=-\infty}^\infty \int_0^1 |g_n(x_2)|^2 \, dx_2
= \sum_{n=-\infty}^\infty \int_n^{n+1} |g_0(x_2)|^2 \, dx_2
= \int_{-\infty}^\infty |g_0(x_2)|^2 \, dx_2,
$$
which is the $L^2$ norm of $g$ on $\bf R$.  So we have an injection
(and probably an isomorphism) of $L^2({\bf R}, {\bf C})$ into the
Hilbert space of functions satisfying ($\star$) which takes
Schwartz functions to smooth functions.
The construction breaks the (anti)symmetry between $x_1$ and $x_2$.
However, if we make $f$ periodic in $x_2$ by multiplying by
$e^{\pi i x_1 x_2}$ instead of $e^{-\pi i x_1 x_2}$,
we end up with much the same formula but with $x_1,x_2$ switched and
$g_0$ replaced by its
Fourier transform!  (This calculation is similar to the one that gives
the Poisson summation formula.)  In particular, for our choice
$g_0(x) = e^{-\pi x^2}$ it is the same function.  This suggests
using the Hermite functions (orthogonal polynomials times $e^{-\pi x^2}$)
for $g_0$ to obtain an orthogonal basis for the space of functions
satisfying ($\star$).
