I think I’ve found and fixed a bug in the FISAPT module regarding external potentials. The monomer SCF calculations within the FISAPT module do not look for external potentials, but the dimer calculation in the SCF module does. If an external potential is present, it is accounted for in the dimer HF energy but not in the monomer HF energies such that the delta HF term is much larger than any of the other SAPT energy terms. Here is an example calculation with a nucleic acid base pair that has 268 electrons in the dimer. The external potential is point charges representing the flanking base pairs in a helix. Energies are in mEh.
Module | Extern? | Elst Exch Ind delta HF Disp Total SAPT0
SAPT | No | -47.7114 55.7805 -18.2558 -8.2544 -13.7969 -23.9836
SAPT | Yes | -47.6277 55.6425 -18.3326 -8.2151 -13.8167 -24.1345
FISAPT | No | -47.7114 55.7804 -18.2557 -8.2531 -13.7969 -23.9835
FISAPT | Yes | -47.7116 55.7805 -165.6920 -155.6894 -13.7969 -171.4200
The SAPT module predicts that the external potential makes the total SAPT0 energy more favorable by 0.15 mEh, while the FISAPT module predicts that the external potential stabilizes the dimer by an unreasonable 147 mEh.
I modified the FISAPT code to allow the monomer SCF calculations to see the external potential. As in the SCF module, I added the interaction of the external potential with the nuclei to the nuclear repulsion and added the interaction with the basis functions to the one-electron Hamiltonian. The latter interaction is stored as an additional Matrix object within the FISAPT class so that it can be added to the one-electron Hamiltonian during the FISAPT SCF calculations and later during the calculation of the delta HF term. I also had to add a method external_pot() to the Wavefunction class in order to access the protected external_pot_ object.
My modified code agrees with the SAPT module in how much the external potential changes the total SAPT0 energy of the dimer.
Module | Extern? | Elst Exch Ind delta HF Disp Total SAPT0
SAPT | No | -47.7114 55.7805 -18.2558 -8.2544 -13.7969 -23.9836
SAPT | Yes | -47.6277 55.6425 -18.3326 -8.2151 13.8167 -24.1345
FISAPT | No | -47.7114 55.7804 -18.2557 -8.2531 -13.7969 -23.9835
FISAPT | Yes | -47.6276 55.6424 -18.3326 -8.2138 -13.8167 -24.1345
I have also validated that my modified code passes the included smoke and sapt tests and that it produces reasonable energies for an I-SAPT calculation with and without an external potential, although I don’t know of a good way to compare the results with the SAPT module. Are there any other checks that I should make to validate that my changes didn’t screw up something in the code? Are the developers interested in merging this feature into the development branch?