FCI Eigenvectors - Githubissues

Hi Josh, The eigenvectors produced by the FCI class do not seem to give correct results when used in eg transition dipole computations. For water in STO-3G using Crawfords geometry - I'm using this geometry as in his projects he gives reference Hamiltonians. Using CIS the first few printed energies and oscillator strengths are

Configuration Interaction Singles (CIS)
------------------------------
# Determinants:  40
nOcc * nVirt:    40
CIS state  1 (eV):       7.8166 (f=0.0000)
CIS state  2 (eV):       7.8166 (f=0.0000)
CIS state  3 (eV):       7.8166 (f=0.0000)
CIS state  4 (eV):       9.3723 (f=0.0000)
CIS state  5 (eV):       9.3723 (f=0.0000)
CIS state  6 (eV):       9.3723 (f=0.0000)
CIS state  7 (eV):       9.6998 (f=0.0023)

Using bitstrings via the commented out code (and making small modifications to the transition dipole computation line) you get

Configuration Interaction Singles (CIS)
------------------------------
# Determinants:  40
nOcc * nVirt:    40
CIS state  1 (eV):       7.8166 (f=0.0000)
CIS state  2 (eV):       7.8166 (f=0.0019)
CIS state  3 (eV):       7.8166 (f=0.0000)
CIS state  4 (eV):       9.3723 (f=0.2856)
CIS state  5 (eV):       9.3723 (f=0.0000)
CIS state  6 (eV):       9.3723 (f=0.0000)
CIS state  7 (eV):       9.6998 (f=0.0000)

The Hamiltonian in the bitstring case is substantially different from CIS because in CIS the 1st row is ordered (if $^a_i$ represents excitation from occupied orbital i to virtual orbital a and we have n occupied and u virtual spin-orbitals) $^1_1$ $^2_1$ ... $^u_1$, $^1_2$ $^2_2$... $^u_2$,..., $^1_n$ $^2_n$ ... $^u_n$ whereas itertools gives an order of $^1_n$ $^2_n$ ... $^u_n$, $^{n-1}_1$ $^{n-1}_2$ ... $^{n-1}_2$ ,..., $^1_1$ $^2_1$ ... $^u_1$ To get the bitstring Hamiltonian to the same form as the CIS one (which is reproduced in Crawford here you can reverse the whole determinant list and then reverse each 'number of virtual orbitals' sub-block. This leads to a Hamiltonian identical in absolute values to the CIS one but with some sign changes (the dominant diagonal is all positive - so this is not just arbitary phase of wavefunction). For example on the first row sixth column the bitstring Hamiltonian has -2.61966629e-02 whereas CIS has 2.61966629e-02 (and Crawford 0.0261967). I'm wondering now if the FCI energies are correct? Given the strong diagonal dominance of the sparse Hamiltonian some sign errors in off-diagonal terms may not show up except in higher precision than we're considering here. It would be good to get this resolved because not having the eigenvectors greatly reduces the usefulness of the bitstring FCI method.

I've used

det_list = list(zip(*(iter(det_list[::-1]),) * 4))
det_list = sum([list(i)[::-1] for i in det_list], [])
for i in det_list:
     print('{:6}     {:20}'.format(i, bin(i)[2:].zfill(nOrb)))

to make to order of the det_list same as CIS ordering. You might have to do some extra processing because your det_list is a list of arrays. What is reason for the array as everything seems to work fine without it? Peter

row col        det A            det B        excitation     degree    phase               value
37   34   00011111111110    00101111111101    [[ 2  2]        2         1         -0.02619666286190461
                                               [ 1  0]
                                               [10 11]]

Clearly the there are 2 new holes and 2 new particles [2, 2], occuring at 1 and 10 and 0 and 11. The level of excitation is clearly 2 and there is 1 permutation needed to move 1->0 and one needed to move 10->11 with no new holes between the excitations => $-1^{2}= +1$ . Can't see anything wrong there.

These are the eigenvectors for CIS singlet root 7 with energy 0.356461697 eV. from CIS

 0.35646169729593263
[0. ... 0.     0.     0.     0.     -0.7071  0.     0.     0.     0.     -0.7071   0.      0.  ]

For bitstring

[ 0.     -0.7071 -0.     0.      0.7071   0.     0. ... 0.  ]

For bitstring re-ordered before Hamiltonian generation

[ 0.  ... 0.     0.     0.     0.     0.7071  0.     0.     0.      0.    -0.7071   0.     0.   ]

Once the combinations have been re-ordered there is a sign change from the original CIS. I think this is due to the spin-scattered approach leading to an error in either the phase calculation or the ordering of the exc array. I'll close this for now and you can re-open if you want too.

I've reopened this having had another look at the issue. Moving to a simpler but non-trivial model - $H_2$ in 3-21g basis.

Again the eigenvalues look correct but oscillator strengths are wildly different. Firstly I've reordered the seed determinant list to give an orbital ordering the same as CIS. This means adding the lines

        det_list = list(zip(*(iter(det_list[::-1]),) * 6))
        det_list = sum([list(i)[::-1] for i in det_list], [])

before the combination loops in build_full_hamiltonian in the postscf module. This gives Hamiltonians

Bitstring
[[ 0.49  0.   -0.    0.   -0.06  0.    0.   -0.09  0.   -0.    0.   -0.08]
 [ 0.    0.4   0.   -0.    0.   -0.14  0.    0.    0.    0.    0.    0.  ]
 [-0.    0.    1.07  0.   -0.    0.    0.   -0.    0.   -0.11  0.   -0.  ]
 [ 0.   -0.    0.    0.95  0.   -0.    0.    0.    0.    0.    0.    0.  ]
 [-0.06  0.   -0.    0.    1.58  0.    0.   -0.08  0.   -0.    0.   -0.13]
 [ 0.   -0.14  0.   -0.    0.    1.45  0.    0.    0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.4   0.   -0.    0.   -0.14  0.  ]
 [-0.09  0.   -0.    0.   -0.08  0.    0.    0.49  0.   -0.    0.   -0.06]
 [ 0.    0.    0.    0.    0.    0.   -0.    0.    0.95  0.   -0.    0.  ]
 [-0.    0.   -0.11  0.   -0.    0.    0.   -0.    0.    1.07  0.   -0.  ]
 [ 0.    0.    0.    0.    0.    0.   -0.14  0.   -0.    0.    1.45  0.  ]
 [-0.08  0.   -0.    0.   -0.13  0.    0.   -0.06  0.   -0.    0.    1.58]]
CIS
[[ 0.49  0.   -0.    0.   -0.06  0.    0.    0.09  0.   -0.    0.    0.08]
 [ 0.    0.4   0.   -0.    0.   -0.14  0.    0.    0.    0.    0.    0.  ]
 [-0.    0.    1.07  0.   -0.    0.    0.    0.    0.    0.11  0.    0.  ]
 [ 0.   -0.    0.    0.95  0.   -0.    0.    0.    0.    0.    0.    0.  ]
 [-0.06  0.   -0.    0.    1.58  0.    0.    0.08  0.    0.    0.    0.13]
 [ 0.   -0.14  0.   -0.    0.    1.45  0.    0.    0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.4   0.   -0.    0.   -0.14  0.  ]
 [ 0.09  0.   -0.    0.    0.08  0.    0.    0.49  0.   -0.    0.   -0.06]
 [ 0.    0.    0.    0.    0.    0.   -0.    0.    0.95  0.   -0.    0.  ]
 [ 0.    0.    0.11  0.    0.    0.    0.   -0.    0.    1.07  0.   -0.  ]
 [ 0.    0.    0.    0.    0.    0.   -0.14  0.   -0.    0.    1.45  0.  ]
 [ 0.08  0.    0.    0.    0.13  0.    0.   -0.06  0.   -0.    0.    1.58]]

You can see [7, 0], [7, 4], [9, 2], [11, 0] and [11, 4] elements have sign differences. The data for these elements is given by

i  j           value            det1         det2                exc         degree    phase
7  0   -0.08898049469427943 000000001001 000000000110    [[2 2] [0 1] [3 2]]    2        1
7  4   -0.08008327234265367 000000001001 000001000010    [[2 2] [0 1] [3 6]]    2        1
9  2   -0.11428224892179743 000000100001 000000010010    [[2 2] [0 1] [5 4]]    2        1
11 0   -0.08008327234265365 000010000001 000000000110    [[2 2] [0 1] [7 2]]    2        1
11 4   -0.13123276680634516 000010000001 000001000010    [[2 2] [0 1] [7 6]]    2        1

It is immediately evident that all the sign differences occur for double excitations. All the correct signs are on single (or no) excitations. This suggests the problem occurs for double excitations. Looking at example [11, 0] the phase from 000000000110 -> 000000000101 -> 000000001001 -> 000000010001 -> 000000100001 -> 000001000001 -> 000010000001 takes 6 steps so $(-1)^6 = +1$ so phase seems right as there are no crossings. Let's look at the generation of exc matrix.

Scemama and Giner describe the algorithm in Figure 5 of their paper An efficient implementation of Slater-Condon rules. We can describe it as

loop over $\alpha$ and $\beta$ spins
find the particles in order
find the holes in order For [11, 0] splitting determinants in spin (000001)(001000) and (000010)(000001) for alpha particle appears at 0 and for beta at 3 giving [0,3] for alpha hole appears at 1 and for beta at 0 giving [1,0] correcting for spatial positions we get particles at [0,3] -> [0, 7] and for holes [1,0]->[2, 1] giving $\langle 2 1 \Vert 0 7 \rangle$ compared to $\langle 0 7 \Vert 1 2 \rangle$ which gives a sign difference!

I think where there are two $alpha$ or two $beta$ spin particles or holes there won't be a problem, however for mixed-spin you need to you need to ensure the even number ($alpha$ spin) comes first. This clearly won't affect the single excitations or no excitation cases.

As a quick and dirty test I added the lines

    if exc[1,1] % 2 == 1: exc[1,1], exc[2,1] = exc[2,1], exc[1,1]
    if exc[1,0] % 2 == 1: exc[1,0], exc[2,0] = exc[2,0], exc[1,0]

at the end of get_double_excitation in slater module, inserting lines before permutation calculations will cause a change of phase which I don't want. This gives

                       CIS                                        bitstring
CIS state  4 (eV):  15.75539775 (f=0.7983)       CIS state  4 (eV):  15.75539775 (f=0.7983)

So it looks like it was the spin-scattered approach that caused the problem but it's quite subtle. It looks also that if you look closely enough the eigenvalues must be a little off and it's only the strong diagonal dominance that is saving the day. I suppose rigorously one should adopt a spin-separated formalism within get_double_excitation everywhere else the spin-scattered method makes life much easier. The re-ordering of the determinant list is only valid to get correct eigenvectors for singles, it will cause errors at higher levels. The fix for spin in get_double_excitations is not the full answer, it leads to the wrong value for full CI in the sto-3g case for $H_2 O$ so something else is going on - but it does fix the oscillator strengths for CIS bitstring in that case. When I did my first calculations on $H_2$ I was getting weird results - diis looks converged but isn't, when I removed diis everything was normal. This seems to be a problem (we've encountered before I recall) with diis on small molecules. I've resorted to not switching diis on until its buffer is full which seems to stablize things, or you could look at the frontier orbitals which I think will look strange. Something for a rainy day maybe.

Hope some of this helps. Peter

@jjgoings ,Hi Josh, I have looked carefully at the CIS by bitstring Hamiltonian and below is a dump of all 2 degree excitations. The listing shows idx, jdx, det1, det2 (in mixed spin), value bitstring, value CIS, np.isclose element, degree, phase and exc det1 alpha and beta det2 alpha and beta holes for alpha and beta and corrected back to mixed-spin particles for alpha and beta and corrected back to mixed-spin double-bar assignment for mixed-spin and separated spin

15 10  10001111110111  01001111111011  -0.06782   0.06782  False    2   1 [2 3] [13 12]
           10001111110111 0011111 1011101              01001111111011 1011101 0011111
              holes alpha  [1] holes beta      [6]                   adjusted alpha  [2]     adjusted beta  [13]
          particles alpha  [6] particles beta  [1]                   adjusted alpha  [12]    adjusted beta  [3]
                                 mixed spin <  2   13  ||  3   12  >          separated spin <  2   13  ||  12   3  >
--------------------------------------------------
16  2  00011111101111  01001111111110   0.00462   0.00462  True     2  -1 [0 4] [10 12]
           00011111101111 0111011 0011111              01001111111110 1011110 0011111
              holes alpha  [0, 5] holes beta      []                   adjusted alpha  [0, 10]     adjusted beta  []
          particles alpha  [2, 6] particles beta  []                   adjusted alpha  [4, 12]     adjusted beta  []
                                 mixed spin <  0   10  ||  4   12  >          separated spin <  0   10  ||  4   12  >
--------------------------------------------------
16  7  00011111101111  10001111111101   0.00304  -0.00304  False    2   1 [1 4] [10 13]
           00011111101111 0111011 0011111              10001111111101 0011111 1011110
              holes alpha  [5] holes beta      [0]                   adjusted alpha  [10]     adjusted beta  [1]
          particles alpha  [2] particles beta  [6]                   adjusted alpha  [4]      adjusted beta  [13]
                                 mixed spin <  1   10  ||  4   13  >          separated spin <  10   1  ||  4   13  >
--------------------------------------------------
16 10  00011111101111  01001111111011  -0.05646  -0.05646  True     2  -1 [2 4] [10 12]
           00011111101111 0111011 0011111              01001111111011 1011101 0011111
              holes alpha  [1, 5] holes beta      []                   adjusted alpha  [2, 10]     adjusted beta  []
          particles alpha  [2, 6] particles beta  []                   adjusted alpha  [4, 12]     adjusted beta  []
                                 mixed spin <  2   10  ||  4   12  >          separated spin <  2   10  ||  4   12  >
--------------------------------------------------

As I observed before, the order of the mixed-spin orbitals are incorrect in some of the cases, but correct in others. Examination of the different cases shows that the wrong sign occurs where a beta-spin orbital is mixed with an alpha-spin on the same side of the ||. This is a consequence of the mixed-spin as Scemama and Giner's algorithms are based on alpha before beta. This makes no difference to the single excitations or to the calculation of phase. Within the mixed-spin approach the solution is simply to ensure if $\langle \beta \alpha \vert$ or $\vert \beta \alpha \rangle$ occurs to swap the $\beta$ and $\alpha$. In get_double_excitation function I have added

    if (exc[1,0] % 2 ) == 1 and (exc[2,0] % 2) != 1: exc[1,0], exc[2,0] = exc[2,0], exc[1,0]
    if (exc[1,1] % 2 ) == 1 and (exc[2,1] % 2) != 1: exc[1,1], exc[2,1] = exc[2,1], exc[1,1]

before the return statement.

Running for $H_2O$ in STO-3G I now get

CIS state  7 (eV):       9.6998 (f=0.0023)
CIS state 15 (eV):      13.7589 (f=0.0649)
CIS state 19 (eV):      15.1075 (f=0.0155)
CIS state 23 (eV):      17.8321 (f=1.2519)
CIS state 24 (eV):       0.9101 (f=0.8488)

compared to from CIS without bitstring

CIS state  7 (eV):       9.6998 (f=0.0023)
CIS state 15 (eV):      13.7589 (f=0.0649)
CIS state 19 (eV):      15.1075 (f=0.0155)
CIS state 23 (eV):      17.8321 (f=1.2519)
CIS state 24 (eV):      24.7657 (f=0.8488)

For HF in STO-3G for bit-string

CIS state  7 (eV):      29.3428 (f=0.0064)
CIS state  8 (eV):      29.3428 (f=0.0064)
CIS state 12 (eV):      34.4850 (f=0.2290)
CIS state 16 (eV):      64.9118 (f=0.2003)
CIS state 20 (eV):     719.0314 (f=0.0339)

and CIS without bitstring

CIS state  7 (eV):      29.3428 (f=0.0064)
CIS state  8 (eV):      29.3428 (f=0.0064)
CIS state 12 (eV):      34.4850 (f=0.2290)
CIS state 16 (eV):      64.9118 (f=0.2003)
CIS state 20 (eV):     719.0314 (f=0.0339)

! This still needs a reordering of the determinant list which generates the Hamiltonian, I've used

        det_list = list(zip(*(iter(det_list[::-1]),) * nVir))
        det_list = sum([list(i)[::-1] for i in det_list], [])

where nVir is the number of spin-virtual orbitals.

BTW A neat one-line for trailz is

(n & -n).bit_length() - 1

Summary: To get CIS by bitstring working I've

Re-ordered the generating determinant list so spin-orbitals will match ordering of MO spin dipoles generated in the usual way, This involves reversing the determinant list and then reversing each 'number of spin virtual orbital' partition of the list.
For degree 2 excitations if a hole or particle pair have the $\beta$ (odd) orbital before the $\alpha$ (even) orbital then reverse the order.

all this is valid only for getting singles

I think I understand why I couldn't use the bit-string CIS eigenvectors to compute transition dipoles and oscillator strengths. I think the bit-string eigenvectors which are calculated in a second-quantisation formalism ie 'true' vacuum cannot be expected to be used with Hartree-Fock ie Fermi vacuum quantities. The normal-order is different for the vacuums so phases (and orbital positions) will be different?

@jjgoings , Hi Josh, I'm closing this now as I understand the problem I had. For MMD the resolution is that the eigenvectors produced via the determinant basis must be re-ordered and the the phase corrected before used with Hartree-Fock (Fermi reference vacuum) quantities. This can for C1 amplitudes be accomplished by adding these lines in the 'bitstring' specific code

            idx = list(zip(*(iter(list(range(nOV))[::-1]),) * nVir))
            idx = np.array(sum([list(i)[::-1] for i in idx], []))
            sign = np.array([1*((-1)**j)  for j in range(nOcc) for i in range(nVir)])

and then adding

            if construction == 'bitstring':
                transition_density = sign * transition_density[idx]

before calculation of the transition dipoles. This gives agreement with CIS (plain) for the oscillator strengths. I've tested on water in STO-3G, 3-21G and DZ basis and $H_2$ in STO-3G and 3-21G and oscillator strengths are in complete agreement.

The FCI is only used to get the ground-state correction so it would acceptable to use a spin-adapted determinant basis. You could determine the $\alpha$ and $\beta$ excitations and reject any $\alpha \rightarrow \beta$ or $\beta \rightarrow \alpha$ excitations, something like

 def  is_spin_adapted(det, nOrb):
            '''
            get the alpha and beta particle numbers from determinant - 
            if no alpha->beta or beta->alpha then spin-adapted
            '''
            occupations = [(int(bool(det & (1 << 2*i))), int(bool(det & (2 << 2*i)))) for i in range(nOrb//2)]
            nalpha, nbeta = [i.count(0) for i in zip(*occupations)]

            return nalpha == nbeta

then simply add the condition

                if is_spin_adapted(det[0], nOrb): det_list.append(det)

This will speed up things eg $H_2 O$ in STO-3G 1001 determinants in basis to 441 and of course the Hamiltonian is half the dimension.

I noticed you have a FIXME in ao2mo for tiling the 2-electron repulsions. I've used

        spin_block = np.kron(np.kron(self.mol.single_bar, spin).transpose(), spin).real
        self.mol.double_bar = spin_block.transpose(0,2,1,3) - spin_block.transpose(0,2,3,1)

which seems to work fine. I've looked at psi4numpys FCI and decided that the spin-scattered treatment is much better and far more concise. Best Peter

Great, thank you for following up on this! Feel free to make a pull request, and I would be happy to merge in these fixes!

jjgoings / McMurchie-Davidson

FCI Eigenvectors #24