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.

<note 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. }}}