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



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](https://mathworld.wolfram.com/WeierstrassSigmaFunction.html) 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.


