# Bugs

## Units in SAPT(DFT) output

The output correctly converts to mH and kcal/mol, but does not convert to kJ/mol. To fix this edit the file

  lib/psi4/driver/procrouting/sapt/sapt_util.py

to make the change to constants.hartree2kJmol:

def print_sapt_var(name, value, short=False, start_spacer="    "):
"""
Converts the incoming value as hartree to a correctly formatted Psi print format.
"""

vals = (name, value * 1000, value * constants.hartree2kcalmol, value * constants.hartree2kJmol)
if short:
return start_spacer + "%-20s % 15.8f [mEh]" % vals[:2]
else:
return start_spacer + "%-20s % 15.8f [mEh] % 15.8f [kcal/mol] % 15.8f [kJ/mol]" % vals

## Cannot use GGAs

The SAPT(DFT) module works correctly when the default, PBE0, functional is used. But when this is changed to a GGA like PBE it fails with the error message:

   => Auxiliary Basis Set <=

Basis Set: ANONYMOUS400E6A7D
Blend: AUG-CC-PVDZ-JKFIT + MB-SET-RI
Number of shells: 478
Number of basis function: 1452
Number of Cartesian functions: 1697
Spherical Harmonics?: true
Max angular momentum: 4

Traceback (most recent call last):
File "/home/alston/Psi4/install/1.2-gcc/bin/psi4", line 259, in <module>
exec(content)
File "<string>", line 71, in <module>
File "/home/alston/Psi4/install/1.2-gcc/lib//psi4/driver/driver.py", line 460, in energy
wfn = procedures['energy'][lowername](lowername, molecule=molecule, '''kwargs)
File "/home/alston/Psi4/install/1.2-gcc/lib//psi4/driver/procrouting/sapt/sapt_proc.py", line 234, in run_sapt_dft
cache = sapt_jk_terms.build_sapt_jk_cache(wfn_A, wfn_B, sapt_jk, True)
File "/home/alston/Psi4/install/1.2-gcc/lib//psi4/driver/procrouting/sapt/sapt_jk_terms.py", line 108, in build_sapt_jk_cache
cache["K_A"] = jk.K()[0].clone()

IndexError: list index out of range

'''* Psi4 encountered an error. Buy a developer more coffee!
'''* Resources and help at github.com/psi4/psi4.

The PBE functional can be selected using:

set {
scf_type  df
sapt_dft_functional     PBE
sapt_dft_do_dHF         False
sapt_dft_grac_shift_a   0.12589
sapt_dft_grac_shift_b   0.12589
}

The problem seems to be in Psi4/install/1.2-gcc/lib/psi4/driver/procrouting/sapt/sapt_jk_terms.py:

    # J and K matrices
jk.C_clear()

# Normal J/K for Monomer A

# Normal J/K for Monomer B

# K_O J/K
C_O_A = core.Matrix.triplet(
cache["D_B"], cache["S"], cache["Cocc_A"], False, False, False)

jk.compute()

# Clone them as the JK object will overwrite.
cache["J_A"] = jk.J()[0].clone()
cache["K_A"] = jk.K()[0].clone()

cache["J_B"] = jk.J()[1].clone()
cache["K_B"] = jk.K()[1].clone()

If a GGA is used, then the exchange matrices are not created so the clone operations on jk.K() fail.

Important

15-03-18

This bug has been fixed by Daniel. Version 1.1-938-gfc51d9a (current master on github) works with GGAs.

## SAPT(DFT) not possible without Asymptotic correction

Curiously, SAPT(DFT) calculations are not possible without the asymptotic correction. In Psi4/install/1.2-gcc/lib/psi4/driver/procrouting/sapt/sapt_proc.py we have the lines:

    if (sapt_dft_functional != "HF") and ((mon_a_shift == 0.0) or (mon_b_shift == 0.0)):
raise ValidationError('SAPT(DFT): must set both "SAPT_DFT_GRAC_SHIFT_A" and "B".')

if (core.get_option('SCF', 'REFERENCE') != 'RHF'):
raise ValidationError('SAPT(DFT) currently only supports restricted references.')

This needs to be changed as it should be possible to switch off the AC completely.

# Testing Psi4

## SAPT(DFT)

Currently, Psi4 can perform SAPT(DFT) calculations only with the ALDA kernel. This is the biggest limitation of the code. The following are the features of Psi4 1.2:

• ALDA kernel
• $E||{(2)}_{\rm ind,exch}$ is computed with the $S||2$ approximation.

• GRAC asymptotic correction.
• Range-separation also a possibility.
• DC or DC+ basis sets.
• Second-order terms only. Higher-order with the $\delta||{\rm HF}_{\rm int}$ approx.

• - What is the accuracy of SAPT(DFT) with the PBE0/GRAC functional with ALDA kernel? - Same for wPBE0(IP-tuned) with ALDA? - How fast is this code? What are the kinds of calculations we may expect to perform? - How well does the code parallelize?

### Accuracy: PBE0+(GRAC/LB94) + ALDA

#### Water dimers

Data from molecules/water/induction-damping/2017-ISA-ISApol-noS2/saptdft/aTZ-MC+/PBE0-GRLB/

• Ref(1) : SAPT(DFT) DALTON+CamCASP: PBE0/GRLB : ALDA+CHF kernel : aTZ MC+BS (Mb=3s2p1d1f)
• Psi4(1) : SAPT(DFT) PBE0/GRLB : ALDA kernel : aTZ DC+BS (Mb=3s2p1d1f). Here the $S||2$ approx is used to calculate the exch-ind energy. This is quite a poor approx for the 900 geom. Notice how the error is also present in the dHF term.

• Psi4(2) : SAPT(DFT) PBE/GRLB : AC-shift = 0.202938 a.u. Other details as with Psi4(1).
• UNITS: kJ/mol

• water2_900: $R_{O..O} = 2.2$ Angstrom

• water2_907: $R_{O..O} = 2.9924$ Angstrom

 water2_900 water2_907 Energy Ref(1) Psi4(1) Psi4(2) Ref(1) Psi4(1) Psi4(2) $E||{(1)}_{\rm elst}$ -160.995 -160.459 -162.198 -28.934 -28.814 -28.317 $E||{(1)}_{\rm exch}$ 378.326 376.891 392.954 25.311 24.687 26.413 $E||{(2)}_{\rm IND}$ -44.758 -71.218 -73.479 -4.474 -4.606 -4.654 $E||{(2)}_{\rm ind}$ -217.284 -215.744 -230.843 -11.054 -10.782 -11.348 $E||{(2)}_{\rm ind,exch}$ 172.526 144.526 157.364 6.580 6.176 6.693 $E||{(2)}_{\rm DISP}$ -44.617 -36.794 -38.893 -8.521 -7.409 -8.324 $E||{(2)}_{\rm int}$ 127.955 108.419 118.383 -16.619 -16.143 -14.882 $\delta||{\rm HF}_{\rm int}$ -56.708 -35.685 -35.685 -2.931 -2.864 -2.864 $E||{(2)}_{\rm int}+\delta||{\rm HF}_{\rm int}$ 71.246 72.734 82.698 -19.550 -19.008 -17.747
• The differences in $E|{(2)}_{\rm DISP}$ is mainly from the kernel. There is also a minor difference because of the basis (MC+BS vs DC+BS). The ALDA kernel results in too little dispersion.

• The differences in $E|{(2)}_{\rm IND}$ is from the $S|2$ approx and the kernel.

• The $S|2$ approx is the major source of error for the induction.

• What is the origin of the differences in the first-order exchange energies?

• The GRLB asymptotic correction in Psi4 uses the LB94+VWN at long-range. I would have thought they should have retained the PBEc functional entirely and corrected only the PBEx functional. What does DALTON do?

Q: How does PBE/ALDA fare (Psi4(2))?

• The second-order energies appear better, but any gains are offset by a too-repulsive first-order exchange energy.
• So there is insufficient binding with PBE.
• Almost want to combine PBE0 first-order with PBE second-order...

### Pyridine dimers

• Ref(1) : PBE0/CS00 ALDA+CHF aDZ+Mb3 MC+BS

• Ref(2) : PBE0/CS00 ALDA aDZ+Mb3 MC+BS

• Psi4-PBE0 : PBE0/GRLB ALDA aDZ+Mb3 DC+BS

• Psi4-PBE : PBE/GRLB ALDA
 Hb1_5 aDZ+Mb3 Energy Ref(1) Ref(2) Psi4-PBE0 Psi4-PBE $E||{(1)}_{\rm elst}$ -20.5598 -20.5598 -20.9766 -20.8048 $E||{(1)}_{\rm exch}$ 28.0175 28.0175 28.5536 31.0925 $E||{(2)}_{\rm IND}$ -3.5901 -3.1628 -3.7196 -3.6984 $E||{(2)}_{\rm ind}$ -8.5829 -7.8752 -8.9282 -9.5670 $E||{(2)}_{\rm ind,exch}$ 4.9928 4.7123 5.2085 5.8686 $E||{(2)}_{\rm DISP}$ -16.2333 -15.3627 -14.1152 -15.6091 $E||{(2)}_{\rm int}$ -12.3657 -11.0679 -10.2580 -9.0198 $\delta||{\rm HF}_{\rm int}$ -2.3259 -2.3259 -2.3422 -2.3422 $E||{(2)}_{\rm int}+\delta||{\rm HF}_{\rm int}$ -14.6917 -13.3939 -12.6002 -11.3620

Data from /home/alston/molecules/pyridine/scans/dimers_from_Model3_noS2/Min2-S1/aTZMb3-MC+/psi4

• Psi4-PBE: Notice how large the first-order exchange is for Psi4-PBE. This seems to be the case in general.
• Psi4-PBE: The second-order energies are closer to the reference values.
• Psi4-hybrid: first-order from Psi4-PBE0 and second from Psi4-PBE: Eint(2) = -11.7573 kJ/mol which becomes -14.0995 kJ/mol with deltaHF. This is much better. But does this hold for other configurations?

Similar data for the Hb1_5 dimer with the aTZ+Mb3 basis:

Cam-PBE0 : CamCASP PBE0/CS00 ALDA+CHF MC+BS basis
Cam-PBE  : CamCASP PBE/CS00 ALDA MC+BS basis
Psi-PBE0 : Psi4 PBE0/GRAC ALDA DC+BS
Psi-PBE  : Psi4 PBE/GRAC ALDA DC+BS
Psi-Hyb  : first-order from Psi-PBE0 and second-order from Psi-PBE

Psi-Cam-FUNC-KER : Psi4 MOs passed to CamCASP with functional FUNC and kernel KER.
Kernels are : ALDA or HYB=ALDA+CHF

kJ/mol                elst        exch         IND2        dHF        DISP2       Eint         CT2        POL2
Cam-PBE0          -20.45658    27.84655    -3.61371    -2.36352   -17.02206   -15.60932    -0.52554    -3.08817
Psi-PBE0          -20.83304    28.52921    -3.68395    -2.36648   -14.53605   -12.89031
Psi-Cam-PBE0-HYB  -20.84369    28.56603    -3.61448    -2.36352   -17.22447   -15.48013    -0.50929    -3.10519
Psi-Cam-PBE0-ALDA -20.84369    28.56603    -3.22050    -2.36352   -15.89760   -13.75929    -0.50848    -2.71203
---------------------------------------------------------------------------------------
Psi-PBE0-scaled   -20.83304    28.52921    -3.68395    -2.36648   -15.89755   -14.25181
Psi-Cam-PBE0-mix  -20.84369    28.56603    -3.61448    -2.36352   -15.89760   -14.15326
---------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------
Cam-PBE           -20.29438    31.92628    -3.44943    -2.36352   -18.02445   -12.20549    -0.51709    -2.93234
Psi-PBE           -20.69212    31.03235    -3.65223    -2.36648   -16.07874   -11.75723
Psi-PBE-scaled    -20.69212    31.03235    -3.65223    -2.36648   -18.06334   -13.74182
Psi-Cam-PBE-ALDA  -20.69365    31.03205    -3.56705    -2.36352   -17.92890   -13.52107    -0.54404    -3.02301
---------------------------------------------------------------------------------------
Psi-Hyb           -20.83304    28.52921    -3.65223    -2.36648   -16.07874   -14.40128
---------------------------------------------------------------------------------------
aug-cc-pVDZ+Mb3 DC+BS Psi4:                                           -15.16522
---------------------------------------------------------------------------------------

Points to note:

• Cam-PBE0 is very close to CCSD(T). The basis set differences are inconsequential here. I have data for Cam-PBE0 with the augA-Sadlej basis too.
• Both Psi4 results are off by a lot.
• I would have expected Cam-PBE to be similar to Psi-PBE, but there are substantial differences in the E(1)exch and E(2)DISP. This could be because the CS00 and GRAC schemes are substantially different at this dimer geometry. Both use the LB94 long-range functional, but differ in how this is spliced to the PBE functional. TO DO: use the Psi4 orbitals in CamCASP. We should get the same energies!

How does Psi4 calculate the second-order energies?

• Induction: CHF+ALDA kernel. This is possible as it uses the iterated solver for the coupled equations. S2 approx used for exch-ind.
• Dispersion: ALDA kernel. Exch-disp is calculated at un-coupled level. No scaling is applied. We have long recommended that scaling should be applied:

$$E||{(2)}_{\rm disp,exch}[C] \approx E||{(2)}_{\rm disp,exch}[UC] \times \frac{E||{(2)}_{\rm disp}[C]}{E||{(2)}_{\rm disp}[UC]}$$

• The Psi-FUNC-scaled results above use this expression and result in much better total interaction energies.

Important

Are CamCASP and Psi4 consistent?

I have now got CamCASP to work with Psi4 1.2 MOs. There were changes needed with the interface and a few other issues.

From the above discussion, we see that we should compare Psi-PBE0-scaled energies to CamCASP calculations computed using the first-order energies and the induction and delta-HF energies taken from Psi-Cam-PBE0-HYB and the second-order dispersion energy from Psi-Cam-PBE0-ALDA. This resulting theory is termed Psi-Cam-PBE0-mix above. And indeed, we see good agreement between the results from the two codes.

There are differences of a tenth of a kJ/mol, but these are due to density-fitting differences and the S||2 approx made by Psi4.

Important

Recommendation:

If we must use Psi4 only, then we should compute the interaction energy using PBE0 (as first-order is good) with scaled exchange-dispersion as shown above.

• This needs to be tested on different dimers.
• The error in the dispersion is still present, but is smaller.

AJMPublic/camcasp/psi4-testing (last edited 2021-04-14 12:53:59 by apw109)