How would Hamiltonian for several fermions with spin look? All discussions of Pauli exclusion principle I read usually talked about antisymmetric wavefunctions, from which the princinple appears. But I would like to see a Hamiltonian for multiple fermions, which would have Pauli exclusion principle already incorporated into its eigenstates. I've looked into Pauli equation, but it seems to only describe one fermion. The same would hold for the more complete Dirac equation.
So, my question is: what equation describes multiple (more than two) fermions with spin in non-relativistic regime, so that its eigenstates would automatically obey Pauli exclusion principle? Do I need the full machinery of QFT for this?
 A: Let's speak about 1D particles for simplicity.
What should be understood first of all is that for indistinguishable particles configuration space isn't the same as for distinguishable ones. For two distinguishable spinless 1D particles configuration space is a square: one side is for $x_1$, another for $x_2$. But if the particles appear indistinguishable, then half of the space is redundant: states obtained by exchange of particle coordinates are identical. So, for two spinless particles configuration space is really a triangle:

Here lower purple triangle is the config space. Now, Schrödinger equation for two indistinguishable spinless 1D particles becomes trivially to obtain from that for distinguishable particles: we just have to impose boundary conditions for $x_1=x_2$ line. For bosons wavefunction must be symmetric, and this implies homogeneous Neumann conditions. For fermions wavefunction must be antisymmetric, so it has a node at the line of interparticle collision, and this means homogeneous Dirichlet conditions.
Now for fermions this automatically gives us the correct eigenstates, which indeed obey Pauli exclusion principle — they just can't have nonzero values on $x_1=x_2$ line. Here're first several1:

If we now want to return to our full (square) configuration space, we should just append zeroes to missing triangle, and then subtract the same wavefunction but with exchanged particles from it. It will automatically appear differentiable at $x_1=x_2$ line:

Now let's include spin in our picture. Thanks to spin being a discrete degree of freedom with finite range of values, for single particle we can just concatenate its configuration space parts corresponding to different spin values. For example, a state of a spin-$\frac12$ particle in infinite box in state $\left|1\uparrow\right\rangle+\left|3\downarrow\right\rangle$ could look like this:

Here left part corresponds to state of spin-up component of wavefunction, and right one is for spin-down component. Note that the Hamiltonian shouldn't try to differentiate the wavefunction at the joint — as far as we neglect spin-dependent interactions, its matrix should just be a diagonal block matrix of two spinless Hamiltonians. Now for two indistinguishable particles nothing changes except we now use 4 times bigger full configuration space (than in case without spin), or 2 times bigger correct triangle-shaped space. Here's how full configuration space would look:

Here spin states are denoted as $\left|s_1\right\rangle\left|s_2\right\rangle$ where $s_i$ is spin of particle $i$. Lines of interparticle collisions are denoted by slight dark color. Now the truncated symmetrized configuration space:

The upper-left rectangle is now merged with lower-right one. How — it will depend on the state. Also note that there's still possibility of interparticle collision — inside this lower-right rectangle — which we haven't symmetrized, and in fact must not do (there's no reason to do it in this case).
Now it's easy to see that the states with identical spins are forced to have antisymmetric orbital — their orbital parts of configuration space are triangular. The states with non-equal spins may be symmetric or antisymmetric (or asymmetric at all). This corresponds to known classification of two-particle states into spin-singlets and spin-triplets. Just looking at the wavefunction, one could immediately say what type of state we have: triplet states will have a node along interparticle collision line — this means the orbital is antisymmetric.
Let's now look at first several states:


Here we can see: spin-singlet states — numbers 1,2,6, and triplets — numbers 3,4,5.
For more than two particles it seems quite straightforward to generalize: just truncate the configuration space so that it is no longer possible to exchange particles, and impose the correct boundary conditions on the hypersurfaces which appear after these cuts. For many particles this may even give some computational resource savings: hypervolume of configuration space would reduce by $n!$ where $n$ is number of particles, as compared to space for distinguishable particles.

1: I used simplified way of computing these states, just putting very large potential barrier in the upper-left triangle. To decouple different spin parts of wavefunctions I also added thin barriers. In real computation one should, of course, use the opportunity to remove extra data from processing and more correctly define the Hamiltonian matrix.
