Obervations and question on the FCI module.

pwborthwick commented 1 year ago

Hi,

I have some observations on the code of the FCI/helper_CI modules. I know coding is a highly personal thing so these are only mentioned in the hope they may be of interest and not in any way as 'improvements' to the existing code.

https://github.com/psi4/psi4numpy/blob/68cd9fb8738eddb718e81e552bae8f481f0a6b24/Configuration-Interaction/helper_CI.py#L58 as an alternative to the existing code you could simply use
```
return bin(bits).count('1')
```
For a million integers the original took 16.6 sec and the alternative 2.7 sec.
https://github.com/psi4/psi4numpy/blob/cbef6ddcb32ccfbf773befea6dc4aaae2b428776/Configuration-Interaction/helper_CI.py#L84 as an alternative you could use list comprehension
```
return [n for n, bit in enumerate(bin(bits)[2:]]::-1]) if bit == '1']
```
For a million integers the original took 22.1 sec and the alternative 11.8 sec.
https://github.com/psi4/psi4numpy/blob/cbef6ddcb32ccfbf773befea6dc4aaae2b428776/Configuration-Interaction/helper_CI.py#L116 as an alternative you could use
```
bit = 0
return sum([obtList ^ (1 << i) for i in obtList])
```
this is only marginally faster than existing routine but saves 9 lines of code.
https://github.com/psi4/psi4numpy/blob/cbef6ddcb32ccfbf773befea6dc4aaae2b428776/Configuration-Interaction/helper_CI.py#L134 as an alternative you could use
```
return [n for n, i in enumerate(Determinant.obtBits2ObtIndexList(bits)) if i in orbitalIndexList]
```
You have two separate routines generateSingleExcitationsOfDet and generateDoubleExcitationsOfDet which are used for calculating CIS and CISD. A small change to the determinant loops in FCI generates CIS
```
ground_state = Determinant(alpha_indexlist=range(ndocc), beta_indexlist=range(ndocc))
detList = []
for alpha in combinations(range(nmo), ndocc):
for beta in combinations(range(nmo), ndocc):
    determinant = Determinant(alphaObtList=alpha, betaObtList=beta)
    if determinant.numberOfDiffOrbitals(ground_state) == 1:
        detList.append(determinant)
```
clearly for CISD the if statement becomes in range(0, 3) and so on for CISDT...Q...P etc. I'm not sure about the efficiency of this approach compared to the custom generateSingleExcitationsOfDet but it does give you CIS etc cheaply?
Are the eigenvectors produced by the FCI bitstring approach immeditaely usable? It seems that for a 'standard' CIS spin-adapted singlet Hamiltonian the ordering of an eigenvector is $_1^6$ $_1^7$ $_2^6$ $_2^7$ ... $_5^6$ $_5^7$ where $_i^a$ is an excitation from occupied orbital i to virtual orbital a. However itertools will produce, I think, an ordering of $_5^7$_5^6$ ... $_1^7$ $_2^6$ which will not be consistant with eg dipole component integrals. Can someone show eg an oscillator strength being computed using FCI bitstring generated eigenvectors.

Thank you for the interesting FCI module.

pwborthwick commented 1 year ago

Hi @tyzhang1993 ,

For HF in STO-3G first 12 roots

CIS                    0.64443155  0.64443155  0.64443155  0.64443155  0.64443155  0.64443155
                       0.71203465  0.71203465  0.79555934  0.79555934  0.79555934  1.05240884

spin-adapted singlets  0.71203465  0.71203465  1.05240884 

psi4numpy CIS(FCI)     0.644431546  0.644431546  0.678233098  0.678233098  0.795559340  1.05240884

The root 0.71203465 given by CIS and spin-adapted and confirmed by using Josh Goings McMurchie-Davidson program with both CIS and his version of FCI is missing and replaced by a root at 0.678233098 in the psi4numpy implementation. Direct eigensolvers have been used throughout. Any thoughts Best Peter

JonathonMisiewicz commented 1 year ago

There are a lot of issues here, and I'd like to get to them, but I have other obligations. I'll investigate on Saturday, at the latest.

I don't understand your last post: which codes generated these three sets of numbers? Links preferred.

pwborthwick commented 1 year ago

There are a lot of issues here, and I'd like to get to them, but I have other obligations. I'll investigate on Saturday, at the latest.

I don't understand your last post: which codes generated these three sets of numbers? Links preferred.

Thanks for your reply. The CIS results are from Josh Goings McMurchie-Davidson repo mmd calculated both from 'standard' CIS and by his bitstring version, also by my own code. The spin-adapted is from my own code and the last set of numbers from the CIS using bitstring from Psi4numpy. For water in a minimal STO-3G basis all sources agree. I'm not running psi4 at the moment so using another source of integrals for the psi4numpy code which I admit is not ideal. However, everything coincides for water, $H_2$ and LiH which I also tried. Thank you for your time, Peter

JonathonMisiewicz commented 1 year ago

Can I have the HF bond length for reproduction purposes?

pwborthwick commented 1 year ago

The actual bond length for the quoted results was part of a dissociation study so was short for HF at 0.74 angstrom.

pwborthwick commented 1 year ago

I've re-run this and psi4numpy is in agreement with other results, No idea what happened first time. So sorry for wasting your time, sincere apologies. Peter

JonathonMisiewicz commented 1 year ago

Thanks. I'll hope to respond to your other points (and clear some of the other Psi4Numpy issues) over the weekend.

pwborthwick commented 1 year ago

Hi @JonathonMisiewicz , Originally I thought the problem was just the ordering of the CI amplitudes produced by the psi4numpy FCI code but now I've got issues with the values of the eigenvectors themselves. I've done the following using psi4

using Crawford Projects geometry, written the CIS spin-adapted singlets method and verified that I get roots in agreement with Crawford Project 12 for CIS.
From mints got AO dipoles and transformed to MO basis and re-shaped to [ndocc, nvirt]
Constructed transition dipole for an arbitary root (0.5056282996 - 3rd singlet) and computed the oscillator strength which all agree with values from my program, Josh Goings McMurchie-Davidson program (which in turn was checked with Gaussian).
Imported CI_helper from psi4numpy FCI module and copied code from it's CIS module.
calculated energies and eigenvectors using FCI/CIS code.
for the root with value 0.5056282996 the eigenvectors can be easily reshaped (nb CIS has different initial determinant order from FCI). There are two contiguous identical blocks [ndocc, nvirt] long. Taking one block and re-normalizing you can then re-shape into a [ndocc, nvirt] matrix which can be compared with the one obtained from CIS spin-adapted method
the matrices are identical in absolute values, but there are differences in sign - I reproduce them below (spin-adapted first)
```
[[ 0.000784  0.      ]
[ 0.064045  0.      ]
[ 0.       -0.344767]
[ 0.936501  0.      ]
[ 0.        0.      ]]
```

[[ 0.000784 0. ] [-0.064045 0. ] [ 0. -0.344767] [-0.936501 0. ] [ 0. 0. ]]


The oscillator strengths do not of course agree. This is the situation for all roots I've tried - right absolute values and some wrong signs. I attach the script I used...[psi.txt](https://github.com/psi4/psi4numpy/files/10103222/psi.txt)  (there might some minor typos in file as I had to un-install psi4 before I'd tidied the file file up as the psi4 installation destroyed my existing Python installation). 
     Peter

JonathonMisiewicz commented 1 year ago

I imagine that the two codes you're using produce the same eigenvector in determinant space but change the signs of their basis vectors and thereby have different expansion coefficients. I strongly suggest you double-check what the basis vectors in those two codes are.

pwborthwick commented 1 year ago

I imagine that the two codes you're using produce the same eigenvector in determinant space but change the signs of their basis vectors and thereby have different expansion coefficients. I strongly suggest you double-check what the basis vectors in those two codes are.

Yes I'm aware that the wavefunction is only defined upto phase so a sign change is possible. However as you will have seen from the script I have used the same set of mo coefficients throughout from a single SCF calculation so all signs should be compatible. All quantities are taken from psi4 and should be compatible in both spin-adapted and FCI/CIS computations. It is not the overall sign of the transition dipole that is the issue, it is that the numerical value of the dipole and of the resultant (signless) oscillator strengths are not in agreement. Thank you for taking an interest.

JonathonMisiewicz commented 1 year ago

It is not the overall sign of the transition dipole that is the issue, it is that the numerical value of the dipole and of the resultant (signless) oscillator strengths are not in agreement.

I was talking about the sign in the CI coefficients.

I imagine that the two codes you're using produce the same eigenvector in determinant space but change the signs of their basis vectors and thereby have different expansion coefficients. I strongly suggest you double-check what the basis vectors in those two codes are.

Yes I'm aware that the wavefunction is only defined upto phase so a sign change is possible. However as you will have seen from the script I have used the same set of mo coefficients throughout from a single SCF calculation so all signs should be compatible.

Why are you assuming that the MO coefficients are the only place for a sign disagreement?

For all I know, the third basis vector in the first program is 1b ^ 2a, and the third basis vector in the second program is 2a ^ 1b. Because these are determinants/wedge-products, these disagree by a sign. Even if the two CI solvers are input the same MO coefficients, the basis vectors disagree by a sign, so the same eigenvector will have different coefficients in different bases.

To be explicit, I'm raising this as the first thing I would check, rather than as the definitive cause.

pwborthwick commented 1 year ago

Thank you for your suggestion. I'm not sure what you mean by 'first' and 'second' programs though? There is only one program with one basis set and one set of mo coefficents. The bottom line is can you, using the CIS module in psi4numpy for singlet root 0.505628 get an oscillator strength of 0.002341 (psi4 values) for water in Crawford geometry and STO-3G basis? I can't get the right value by this method although can by a plethora of other methods :) Anyway thank you for your suggestions and help and giving some of your time to this. I suspect your feeling is that there is something wrong with what I am doing and not with the psi4numpy code so I will close this now.

psi4 / psi4numpy

Obervations and question on the FCI module. #131