Bug in the E(30)ind SAPT correction?

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


It would be nice if you could open up in issue on github. There it will get a lot more exposure to the devs and it is easier to ping those who might know the part of the code.

You (either of you) are also very welcome to directly make a pull-request on github with the changes and description of the problem.

edit: On github you can also find the psi4 dev slack channel for more direct feedback.

No problem, the new issue has been opened!