Bonded ALMO-EDA jobs failing at last step

I’m trying to reproduce some results from the PNAS paper titled ’ Energy decomposition analysis of single bonds within Kohn–Sham density functional theory

Two part question:

1) As a start I have used example input 12.23 given in the Q-Chem 5.3 manual for the fluorine diatomic (changing the method and basis-set to wB97M-V/aug-cc-pVTZ as is specified in the PNAS paper).This job runs fine, however there is a significant difference in my results to the reported values, mainly coming from the spin-coupling, and charge-transfer terms, resulting in a difference in the total interaction energy of 10.4 kcal/mol:

Bonded-EDA (wB97M-V/aug-cc-pVTZ) - fluorine diatomic (All values in kcal/mol)

My results
E(Prep) = +9.2
E(FRZ) = +186.3
E(SC) = -129.6
E(POL) = -38.8
E(CT) = -77.6
E(Total) = -50.5

PNAS Results (Table 3)
E(Prep) = +9.2
E(FRZ) = +186.3
E(SC) = -124.1
E(POL) = -37.2
E(CT) = -74.3
E(Total) = -40.1

This difference is quite large. Am I missing something from the paper? No other details are provided in the methods section.

This brings me to my second question.

2) If I try to increase the SCF convergence threshold I get hung up on the SCF convergence for the charge-transfer term on the super system. The job fails - Segmentation fault - Error in the serial run.

This happens for almost every molecule I have tried. For those jobs that do finish successfully if I use the same input and simply change the molecule (replace F2 for Cl2) the job fails at the same SCF step for charge-transfer of the super-system. This mainly happens when is use basis sets of triple-zeta quality. Maybe I’m missing something but i cant understand why a job that finishes successfully for F2 would fail for Cl2.

This may be naive, but it seems to me that the bonded EDA is a lot less ‘black-box’ than the non-bonded ALMO-EDA model. What sort of things should i be looking-out for when preparing the input for these jobs? Any advice from the ALMO-EDA developers would help.

Thank You

I looked into the inconsistency between the results you obtained (I was also able to reproduce the same data) and the numbers in the PNAS paper. It seems that the raw data given by q-chem has not taken the scaling factor (defined by eq. 9 in that paper) into account. In the SI of that paper, you can find that the scaling for F-F with wB97M-V/aTZ is 0.958, and multiplying the E_SC, E_POL, and E_CT terms you got with this factor will reproduce the result in the PNAS paper. Note that according to eq. 9, this scaling factor is not empirical but determined by several single-point energies.

I was also able to reproduce the segmentation fault that occasionally occurs in the CT step, which seems to be a bug in the code. Further investigation is needed.

Indeed, the bonded EDA feature is less mature than the non-bonded one, which is in part because the main developer of that method is no longer active in the Q-Chem developer community. Some necessary code improvement and maintenance are definitely needed. Eventually we will have to modernize this code, making it more robust and easier to use.

Will this be fixed (anytime soon) in an upcoming version? I reckon that the major problem is SCF convergence. Will there be other SCF algorithms rather than gdm_ls become available?

I forgot to mention that it is not fixed in the latest version (5.4.2)