FCI Eigenvectors about mcmurchie-davidson

[pwborthwick]

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.

[pwborthwick]

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

[pwborthwick]

@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

[pwborthwick]

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?

[pwborthwick]

@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

[jjgoings]

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!