Dirac equation on general geometries? I have a numerical method for computing solutions to the Dirac equation for a spin 1/2 particle constrained to an arbitrary surface and am interested in finding applications where the configuration space has a complicated geometry, i.e., not just $R^2$ or the sphere, but a more general curved surface possibly with special boundary conditions.  On the sphere, for instance, one can take advantage of symmetry and simply use the spherical spinors, but for a more general configuration space (or one described only by measurements, for instance) it may not be possible to come up with a nice, closed-form solution.  However, I am not a physicist and would like to better understand where (or indeed, if) such problems arise.  Any pointers are greatly appreciated.  Thanks!
Edit: Note that I already have a numerical method for solving the Dirac equation -- I am not looking for information on how to solve it.  (And, interestingly enough, Clifford algebra / geometric algebra is already the starting point for what we do.)  Also, I am specifically interested in the case of surfaces, i.e., 2-manifolds, so (many) applications in GR probably don't apply.  Thanks!
 A: The lagrangian for a massless Dirac fermion in a general background is given by:
$$ L_D = i (\bar \psi e^{\mu}_I \gamma^I \mathcal{D}_{\mu} \psi + c.c.) $$
where $c.c.$ means complex conjugate. $\bar \psi = \gamma^0 \psi^+$, $e^{\mu}_I$ is the tetrad (or more generally n-bien) encoding the background metric via $g_{\mu\nu} = e^I_{\mu} e^J_{\mu} \eta_{IJ}$, where $\eta_{IJ}$ is the minkowski metric $diag(-1,1,1,1)$. $\mathcal{D}_{\mu}$ is the covariant derivative whose action on spinors is given by:
$$ \mathcal{D}_{\mu} \psi = \partial_{\mu} \psi + g A_{\mu}^{IJ}\gamma_I \gamma_J \psi $$
Here $\gamma_I$ are the Dirac matrices and $\gamma_I\gamma_J$ are the generators of the Lorentz lie algebra. For spinors which transform under the action of a general group the second term in the above expression can be generalized to $A_{\mu}^I T_I \psi$ where $T_I$ are the lie-algebra generators of the relevant group.
The n-bien encodes the background metric. The gauge connection encodes the background curvature via $F_{\mu\nu}^{IJ} = \partial_{[\mu}A_{\nu]}^{IJ} + g [A_{\mu}, A_{\nu}]^{IJ} $. The $[..]$ in the first term indicates antisymmetrization over the contained indices. The second term is the commutator of the matrices which constitute the connection and is non-zero only for a non-abelian group. The complete Lagrangian including that for the background geometry and the Dirac fields is:
$$ L_{GR + D} = \frac{1}{8\pi G} e \wedge e \wedge F + L_D $$
Of course, if the background geometry is or can be taken to be static then only the Dirac term matters. The corresponding equations of motion are easily determined via variation w.r.t $\psi$ and $\bar \psi$. The complete solution must take into account not only the bulk values of the n-bien and connection but also the boundary conditions imposed by the background geometry.
In the above, what is left unsaid is how the $e \wedge e \wedge F$ term corresponds to the Einstein-Hilbert lagrangian. This is essentially the connection formulation of GR for which an excellent pedagogical reference is the paper by Romano Geometrodynamics vs. Connection Dynamics.
Given all this formalism what is lacking is a concrete example. I have no doubt someone (@Lawrence ?) will soon supply one. But this should get you going. You can also go the Clifford/geometric algebra route @Carl mentioned and gain the benefits of a unified language at the expense of some time spent learning the framework. Hope this helps!

An explicit example of an exact solution can be found in this paper on Graphene Wormholes
A: To write down Dirac's equation over curved spacetime, first express the geometry in terms of vierbeins and spin connections. Spinor bundles are vector bundles over spacetime transforming locally under the local Lorentz gauge group. Using the spin connection, we can write down covariant derivatives for sections of the spinor bundle. To get the Dirac operator, we contract the covariant derivative with the inverse vierbein and contract that with the gamma matrices.
A: You might take a look at the version of general relativity (GR) done by the Cambridge Geometry Group. They translated GR into "geometric algebra", which amounts to the gamma matrices. So it's easy for them to do Dirac equation calculations on a black hole. Here's an example papers:
In P.G. Bergmann and V. de Sabbata eds, Advances in the Interplay Between Quantum and Gravity Physics, 251-283, Kluwer (2002), Anthony Lasenby and Chris Doran, Geometric Algebra, Dirac Wavefunctions and Black Holes
http://www.mrao.cam.ac.uk/~clifford/publications/abstracts/anl_erice_2001.html
Phys.Rev. D66 (2002) 024006, Chris Doran, Anthony Lasenby, Perturbation Theory Calculation of the Black Hole Elastic Scattering Cross Section
http://arxiv.org/abs/gr-qc/0106039v1
Dirac equation for the Kerr metric (rotating black hole):
Phys.Rev. D61 (2000) 067503, Chris Doran, A new form of the Kerr solution
http://arxiv.org/abs/gr-qc/9910099v3
These, and papers they cite or are cited by, should be sufficient.
A: I've never seen, and i think it would be fun to solve it inside the klein-bottle. Basically its a closed version of the Möebius strip (meaning that is also non-orientable); the boundary conditions are pretty simple
take a square, with vertices A,B,C,D and edges AB, BC,CD,DA. Now make the following identifications:
1) (A, AB, B) => (D, DC, C)
this gives you a periodic condition on x, making your square a cylinder
now if you would do (B, BC, C) => (A, AD, D) you would get back a torus $T^2$
if you do instead 
2) (B , BC, C) => (D, DA, A) you will get a surface homeomorphic to the klein-bottle
