Spin Contamination Appearance after Truncating Wavefunction

[tbhuber]

I am using CIPSI for all-electron energy calculations.

When I build a wavefunction for a CH3 radical with a 'fixed' number of determinants, I obtain an S^2 value of 0.75. I expected this value, since I am forcing the wavefunction to be an eigenfunction of S^2. But once I truncate my wavefunction, the S^2 value becomes 1.00 (qp_run truncate_wf_spin based on ci_threshhold). I observe the same results for a plethora of wavefunctions with different quantity of determinants and different molecules. See below for example of fci output and pt2 output respectively.

Please correct me if I am wrong, but there shouldn't be any triplet (or higher) states during the truncating step since there were none generated in the previous step of building the wave function. Is it possible that the S^2 value is being miscalculated during the pt2 calculation?

N_det = 4626334
N_states = 1
N_sop = 311088

• State 1
  < S^2 > = 0.75000002079000971
  E = -39.811606632198909
  Variance = 7.1084719944476002E-002
  PT norm = 1.4633135370513599E-002
  PT2 = -2.7270756336503156E-003
  rPT2 = -2.7264918136369214E-003
  E+PT2 = -39.814333707832560 +/- 5.4679372548341915E-006
  E+rPT2 = -39.814333124012549 +/- 5.4667666634644601E-006

N_det = 1855164
N_states = 1
N_sop = 263186

• State 1
  < S^2 > = 1.0000241195616941
  E = -39.811260478586149
  Variance = 9.6033003388409588E-002
  PT norm = 1.5253829569668725E-002
  PT2 = -3.2844583094855159E-003
  rPT2 = -3.2836942617489638E-003
  E+PT2 = -39.814544936895636 +/- 6.4219870258990843E-006
  E+rPT2 = -39.814544172847896 +/- 6.4204931099504075E-006

[kgasperich]

Thanks for pointing this out, we'll look into it
@anbenali
@TApplencourt

[kgasperich]

@tbhuber do you have a small example so we can be sure we're reproducing the same issue?
thanks!

[kgasperich]

Here is a small example that shows this.
The <S^2> in the fci and pt2 output don't match.
pt2-s2-issue.tar.gz

[eginer]

By truncating the wave function you clearly destroy the eigenfunction of S2 because you throw away a lot of "small" determinants which are useless for the energy but very important for S2. But I guess that if you would reminimize the energy within the truncated set of determinants, without the constraint of being an eigenfunction of S2, you will obtain something much closer to 0.75 . So as a test, once that you have truncated the wave function, you can unset the eigenfunction of S2 (check the keywords only_expected_s2 in davidson and s2_eig in determinants, set them to false). Then, just run the executable diagonalize_h and it will minimize the energy by optimizing the CI coefficients. Then, I expecy that the value of S2 will be closer to 0.75. If it is not the case, then set to True the keuwords I mentioned above, but the Davidson might not converge. Best, Manu Le jeu. 20 août 2020 à 23:33, Kevin Gasperich <notifications@github.com> a écrit :

…

Here is a small example that shows this. The <S^2> in the fci and pt2 output don't match. pt2-s2-issue.tar.gz <https://github.com/QuantumPackage/qp2/files/5105583/pt2-s2-issue.tar.gz> — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#126 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABWLQJZRXGPX25ZSB3ZI3R3SBWJCDANCNFSM4QBYRQGQ> .

[kgasperich]

We have identified the issue, and it should be fixed soon.
It looks like it's just a problem with what's being printed, so your wavefunction and energy are all still correct.

When pt2 is called, it builds psi_s2 from the provider, and this sets it to <S^2> - S_z*(S_z-1) (or <S_+ S_-> in terms of ladder operators).

When fci is called, it sets psi_s2 equal to CI_S2, and this is shifted by S_z*(S_z-1)
The shift is applied here for full diagonalization, and here for the davidson.

So the <S^2> in the fci output is actually <S^2>, and the value from the pt2 output is <S^2> - S_z*(S_z-1).

[kgasperich]

for the high-spin cases (i.e. max Sz for that multiplicity), you should see:

fci   pt2
0.75  1
2     2
3.75  3
6     4
8.75  5

[tbhuber]

Thank you gentlemen. I appreciate the time and effort you have invested in investigating/resolving the issue.

[kgasperich]

@scemama
One way to make everything consistent would be to add S_z2_Sz in diag_S_mat_elem
so line 28 would be diag_S_mat_elem = dble(nup) + S_z2_Sz.

It might make the function a bit heavier (I think it will add some if(.not.s_z2_sz_is_built) ...), but if this is a problem then we could pass the S_z2_Sz as an additional argument.
This might be easier than remembering that diag_S_mat_elem actually returns <S^2> - S_z*(S_z-1) and applying the shift anywhere it is used.

If we make this change, we should also remove the shift in u_0_S2_u_0, davidson_diag_hjj_sjj (1) and (2), and the provider for CI_s2.

[scemama]

OK. Maybe the best would be to remove the variable S_z2_Sz and to compute it directly in the diagonal elements of S^2. I will do it when I come back from vacation.
Thanks!

[TApplencourt]

It's should be fixed after merging #127.
Can you try to git pull and see if it fixed your problem @tbhuber ?
Thanks!

[tbhuber]

@TApplencourt Unfortunately, I have qp2 installed on a cluster owned by a university and do not have privilege to git pull. They do not support using head versions, but I am currently seeking a work-a-round.

[kgasperich]

you could try wget https://github.com/QuantumPackage/qp2/archive/dev.zip