Dear Psi4 Team,
While testing out some enhancements to third-order SAPT, Jonathan Waldrop and I seem to have found a bug in the E(30)ind implementation in psi4/psi4/src/psi4/libsapt_solver/ind30.cc.
The E(30)ind correction is given by Eq. (72) of JCP 125, 154107 (2006), which stems from evaluating the matrix elements in Eq. (75) using the second-order induction amplitudes in Eq. (79). Psi4 proceeds slightly differently, by actually evaluating the induction amplitudes from Eq. (79) (as far as I can see, the code in psi4/psi4/src/psi4/libsapt_solver/amplitudes.cc matches Eq. (79) exactly) and contracting them with the \omega electrostatic potential matrix elements in psi4/psi4/src/psi4/libsapt_solver/ind30.cc in SAPT2p3::ind30(). The problem is that this algorithm accounts for the first 2 matrix elements in Eq. (75) but not for the last one, <VS_AS_B>. Unless I understand something incorrectly, the Psi4 code is missing the <VS_AS_B> contribution in E(30)ind.
Going now back to Eq. (72), it is easy to prove that <VS_AS_B>=4s^r_a v^{ab}_{rs} s^s_b, that is, it contributes a quarter of the last term in Eq. (72). Our working hypothesis was that, because of the omission of this term, E(30)ind programmed in Psi4 has the last term in Eq. (72) multiplied by 12 instead of 16. Therefore, Jonathan modified his Psi4NumPy E(30)ind code to include this last term with the incorrect factor of 12, and now his values are matching Psi4 perfectly! In other words, we are quite certain that the Psi4 implementation includes the last term in Eq. (72) with an erroneous factor of 12 instead of 16. Fortunately, this doesn’t affect the SAPT2+(3) or SAPT2+3 results that also include delta HF, but it will (slightly) affect the third-order SAPT results computed without delta HF.
I think the E(30)ind,resp code is fine but the wrong E(30)ind will slightly affect the scaled approximated value of E(30)exch-ind,resp.
Please let me know if I need to explain anything better or if I can be of any further help.
All the best,
-Konrad Patkowski