1. Initialize \( H(\alpha,\beta)=0.0 \)
2. Loop over \( \alpha = 1, NSD \)
3. For each \( \alpha \), loop over \( a=1,ntbme \) and fetch \( V(a) \) and the single-particle indices \( p,q,r,s \)
4. If \( V(a) = 0 \) skip. Otherwise, apply \( \hat{a}^\dagger_p\hat{a}^\dagger_q \hat{a}_s \hat{a}_r \) to the Slater determinant labeled by \( \alpha \).
5. Find, if any, the label \( \beta \) of the resulting Slater determinant and the phase (which is 0, +1, -1).
6. If phase \( \neq 0 \), then update \( H(\alpha,\beta) \) as \( H(\alpha,\beta) + phase*V(a) \). The sum is important because multiple operators might contribute to the same matrix element.
7. Continue loop over \( a \)
8. Continue loop over \( \alpha \).

You should force the resulting matrix \( H \) to be symmetric. To do this, when updating \( H(\alpha,\beta) \), if \( \alpha \neq \beta \), also update \( H(\beta,\alpha) \).