Navigation:

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.

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:

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:

Questions to be answered:

Accuracy: PBE0+(GRAC/LB94) + ALDA

Water dimers

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

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

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

Pyridine dimers

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

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

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:

How does Psi4 calculate the second-order energies?

$$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]}$$

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.