Navigation:

• CamCASP

• AJM main

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:

sapt_util.py (fragment)

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:

sapt_jk_terms.py

    # J and K matrices
    jk.C_clear()

    # Normal J/K for Monomer A
    jk.C_left_add(wfn_A.Ca_subset("SO", "OCC"))
    jk.C_right_add(wfn_A.Ca_subset("SO", "OCC"))

    # Normal J/K for Monomer B
    jk.C_left_add(wfn_B.Ca_subset("SO", "OCC"))
    jk.C_right_add(wfn_B.Ca_subset("SO", "OCC"))

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

    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.

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:

sapt_proc.py

    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.

Questions to be answered:

• - 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

• 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

Data from: /home/alston/molecules/pyridine/scans/dimers_from_Model3_noS2/Min1-Hb1/aDZMb3-MC+/psi4

• 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

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
---------------------------------------------------------------------------------------
CCSD(T) augA-Sadlej+Mb3 DC+BS Psi4:                                           -15.84991
        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.