I would like to post a partial answer that I have for now. (sorry I am not sure how to type in equation here, this is copied from my answer in here https://mattermodeling.stackexchange.com/posts/1941/edit)
For the indices ijkl in the FCIDUMP file, it corresponds to the integral (see https://theochem.github.io/horton/2.0.2/user_hamiltonian_io.html)
$$
\tag{1}
\int dx_1dx_2\chi_i^(x_1)\chi_k^(x_2)\chi_j(x_1)\chi_l(x_2)
$$
while the indices in the paper, it is
$$
\tag{2}
\int dx_1dx_2\chi_i^(x_1)\chi_j^(x_2)\chi_k(x_2)\chi_l(x_1)
$$
I have ignored some other factors which are not relevant here.
We see that there are two differences, one is the order of the indices, and the other is the $x_{1,2}$ arguments. Thus we have the following mapping
ijkl in FCIDUMP file = iklj in the paper, for the spatial orb
Further, with the second equation above, we have the following symmetry property
$$
\tag{3}
ijkl = jilk = lkji
$$
where the second equality is due to reality requirement of the integral. With these, we are ready to understand the integrals in FCIDUMP file one by one. First, we need to subtract all the indices by 1, so that it matches with those in the paper. So we have
6.74493103326006093745E-01 0 0 0 0
6.63472044860555665302E-01 0 0 1 1
6.63472044860555665302E-01 1 1 0 0
6.97397949820408036281E-01 1 1 1 1
1.81287535812332034624E-01 1 0 1 0
-1.25247730398215462166E+00 0 0 -1 -1
-4.75934461144127241017E-01 1 1 -1 -1
7.43077168397780152276E-01 -1 -1 -1 -1
Next, the very last row is the nuclear repulsion energy, which is not of our interest here. The 2nd and 3rd to last rows are the 1-body integral, which are pretty easy to understand. So let’s focus on the 2-body integral.
$0000$. After using the mapping in the 3rd Equation, it is the same for the labeling of the spatial orb in the paper. The first two (from the right) means 2 0-th spatial orbs are occupied, and this can only be the case if both spin-up and down are occupied. Using the labeling of the spin-orbital, it would be 10 (or 01). Similarly for the next two indices. Thus we have 0110 in the spin-orbital basis (why it is not 1010, I am not sure, and maybe this is due to convention). By the symmetry in Eq. 4, we have also 1001.
$1111$. The argument is essentially the same as above, except that now we are dealing with 1st spatial orb, with the spin-up and down labeled as 23. Thus we have 3223 and 2332 for the spin-orb in the paper.
$0011$. Now with the mapping in Eq. 3, it is in fact 0110 for the spatial orb in the paper. The first two indices 10 means the 0th and 1st spatial orbs are occupied, and they could be either spin up and down. Thus we have $2\times2=4$ options: 20, 30, 21, 31 for the first two indices of the spin-orbs. Thus put it together, we have 0220, 0330, 1221, 1331 for the spin orbs. Again, I am not sure why we don’t have 2020, maybe due to convention.
$1100$. This is essentially the same as above, where we realize that it is 1001 for the spatial orb in the paper. With the same logic, we have 2002, 3003, 2112, and 3113. These can essentially be obtained with the symmetry property in Eq. 4.
$1010$. Ok, I am stuck here for now… Will update this after I figure it out.
