Find an element with given periodicity

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.

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$).