I have a quantum mechanical system of two interacting electrons in one spatial dimension. The Hamiltonian of the system is of the form $H = h + \frac{1}{|x_1 - x_2|}$, where $h$ is a one-electron part irrelevant to the current question, and $x_1, x_2 \in [0, a],\, a > 0$ are spatial coordinates of the electrons in the system. I have decided to use configuration interaction method(Ritz method if you will) to solve the ground state energy of this system numerically. This involves expanding the exact two-electron wave function in a basis of Slater determinants $$\frac{1}{\sqrt{2}}\begin{vmatrix} \chi_i(\mathbf{x}_1) & \chi_j(\mathbf{x}_1) \\ \chi_i(\mathbf{x}_2) & \chi_j(\mathbf{x}_2) \end{vmatrix}.$$ $\chi_i$ is a spin-"orbital", which is simply a product of a one-electron wave function $\psi_i$ and either spin eigenfunction $\alpha$ or $\beta$ (corresponding to spin-up and spin-down, this notation is from Modern Quantum Chemistry by Szabo and Ostlund). $\alpha$ and $\beta$ are functions of a spin-variable $\omega$ and they are orthonormal wrt an inner product which is denoted by an integral, so e.g. $\int \alpha^*(\omega)\beta(\omega) \,d\omega = 0 $ (bit of an abuse of notation). In addition, $\mathbf{x}_i = (x_i, \omega_i)$. Below I use overbar to denote the spin-part of a spin-orbital so that $\chi_i = \psi_i\alpha$ and $\bar{\chi}_i = \psi_i\beta$.
The problem I have lies in the evaluation of the integrals needed to calculate the Hamiltonian matrix elements for the configuration interaction method. One type of these integrals is $$I(\bar{\chi}_i, \chi_j, \bar{\chi}_k, \chi_l) \equiv \int d\mathbf{x}_1d\mathbf{x}_2 \big[\bar{\chi}_i(\mathbf{x}_1)\chi_j(\mathbf{x}_2) - \chi_j(\mathbf{x}_1)\bar{\chi}_i(\mathbf{x}_2)\big] \frac{\bar{\chi}_k(\mathbf{x}_1)\chi_l(\mathbf{x}_2) - \chi_l(\mathbf{x}_1)\bar{\chi}_k(\mathbf{x}_2)}{|x_1-x_2|} \tag{1},$$ which involve Slater determinants formed by both a spin-orbital with spin-up and a spin-orbital with spin-down. Now, if I were to directly use the orthonormality of $\alpha$ and $\beta$, I would get $$I(\bar{\chi}_i, \chi_j, \bar{\chi}_k, \chi_l) = \int \frac{dx_1dx_2}{|x_1-x_2|} \big[\psi_i(x_1)\psi_j(x_2)\psi_k(x_1)\psi_l(x_2) + \psi_j(x_1)\psi_i(x_2)\psi_l(x_1)\psi_k(x_2)\big] \tag{2},$$ which very much looks like a diverging integral to me since $\frac{1}{|x_1-x_2|}$ has a singularity at $x_1 = x_2$, while the nominator of the integrand does not vanish at $x_1 = x_2$. Integrals such as $I(\chi_i, \chi_j, \chi_k, \chi_l)$ and $I(\bar{\chi}_i, \bar{\chi}_j, \bar{\chi}_k, \bar{\chi}_l)$ can be simplified to $$ \int dx_1dx_2 \big[\psi_i(x_1)\psi_j(x_2) - \psi_j(x_1)\psi_i(x_2)\big] \frac{\psi_k(x_1)\psi_l(x_2) - \psi_l(x_1)\psi_k(x_2)}{|x_1-x_2|}.$$ which is finite thanks to the antisymmetrized wavefunctions, which vanish at $x_1=x_2$, in the integrand.
Assuming steps between (1) and (2) aren't justified, I have tried other approaches for calculating $(1)$ and found that $$ \int dx_1dx_2 \big[\psi_i(x_1)\psi_j(x_2) - \psi_j(x_1)\psi_i(x_2)\big] \frac{\psi_k(x_1)\psi_l(x_2) - \psi_l(x_1)\psi_k(x_2)}{|x_1-x_2|} = I(\bar{\chi}_i, \chi_j, \bar{\chi}_k, \chi_l) + I(\bar{\chi}_i, \chi_j, {\chi}_k, \bar{\chi}_l) $$ but this doesn't seem to lead anywhere.
Is there a way to express $I(\bar{\chi}_i, \chi_j, \bar{\chi}_k, \chi_l)$ as a convergent integral involving only spatial wave functions $\psi_i$ or is this model/formalism unfounded applied to this system?
Edit: To give a bit more background, I am trying to get results similar to the ones obtained in this paper. One of the methods they use in that paper is actually Ritz method, but with a bit different basis set and and without spin.
I am also aware of the close relationship to 1D hydrogen atom which was pointed out by Vadim in their answer. I believe I should be able to get finite results in this case(see the earlier paper). And after thinking about it more, I think my best bet is to use a different set of basis functions.