solve_fermion returns incorrect energy when spin is contaminated

[kevinsung]

What should we add?

PySCF is not guaranteed to return a wavefunction with the requested spin squared. When the spin squared deviates from the value requested, then e_sci in this line is not the energy of the wavefunction.

qiskit-addon-sqd/qiskit_addon_sqd/fermion.py

Line 135 in 6a882f6

e_sci, sci_vec = fci.selected_ci.kernel_fixed_space(

Mario investigated and found this to be apparently the correct formula:

e_sci = e_sci-myci.ss_penalty*(spin_sq-spin_squared)**2

But there may be situations where it is not correct. In PySCF's normal FCI class, there is an e_tot attribute which is supposed to always return the energy, but SelectedCI doesn't have that. The surest thing to do is just to compute the energy directly from the RDMs:

dm1 = myci.make_rdm1(sci_vec, norb, (num_up, num_dn))
dm2 = myci.make_rdm2(sci_vec, norb, (num_up, num_dn))
e_sci = np.einsum('pr,pr->',dm1,hcore)+0.5*np.einsum('prqs,prqs->',dm2,eri)

(courtesy of Mario)

TODO construct a reproduction of this issue so it can be tested for

[caleb-johnson]

Thanks Kevin, the fermion module testing is light/nonexistent. I'll take a look at this and work up a unit test for it