(Probably) wrong results for sSAPT0, SAPT0 and SAPT2 with multiple basis sets

1

Hi there,

I’ve calculated sSAPT0 and SAPT0 for triiodide (fragmented as “I2…I-”) with the all-electron, relativistic contracted TZP-DKH basis set (retrieved from EMSL), but the results look very strange:

SAPT0/TZP-DKH (kcal/mol):

Bonding:        -511.92
                -------
Exchange:         77.18
Electrostatic:  -444.83 (75.51%)
Induction:      -136.44 (23.16%)
Dispersion:       -7.83 ( 1.32%)

sSAPT0/TZP-DKH (kcal/mol):

Bonding:        -419.99
                -------
Exchange:         77.18
Electrostatic:  -444.83 (89.47%)
Induction:       -46.56 ( 9.36%)
Dispersion:       -5.78 ( 1.16%)

This looks wrong, specially when I compare this interaction energy with EDA (ZORA-BP86-D3/TZ2P using ADF2016):

Bonding:         -49.75
                -------
Pauli Repulsion:  58.11
Electrostatic:   -45.17 (41.87%)
Orbital:         -59.82 (55.45%)
Dispersion:       -2.88 ( 2.66%)

Interaction energies differ by over an order of magnitude!
SAPT2/TZP-DKH does not help:

Bonding:        -515.43
                -------
Exchange:         77.61
Electrostatic:  -446.30 (75.25%)
Induction:      -138.91 (23.42%)
Dispersion:       -7.83 ( 1.32%)

Using a bigger/different basis set (ANO-RCC) does not help either:

sSAPT0/ANO-RCC (kcal/mol):

Bonding:        -422.54
                -------
Exchange:         85.27
Electrostatic:  -443.80 (87.39%)
Induction:       -49.12 ( 9.67%)
Dispersion:      -14.89 ( 2.93%)

I made a series of calculations changing a lot of things (slightly displacing atoms, changing the initial guess to GWH, changing the auxiliary basis set from Weigand to Ahlrich) and I could not get the two results (EDA and SAPT) to (minimally) agree with each other.

Finally, I would like to know if this could be a bug, specially when we look at the table below:

            Interation energies (kcal/mol)

 Structure  | sSAPT0/TZP-DKH @ Psi4 | EDA (ZORA-BP86-D3/TZ2P) @ ADF
   I-...I2  |   -511.92             |    -49.75
Me2S...ICl  |    -25.00             |    -28.71
Me2Se...ICl |    -20.98             |    -31.20
Me2Se...I2  |    -19.45             |    -26.62

(The above is only a piece of the whole table; the original contains at least three times more rows, all agreeing with each other in magnitude, except for I-…I2).

My questions are:

(a) Could this be correct? Could this be a bug?
(b) Can this be due to: (i) charge or (ii) symmetry?

All the best,
Felipe

P.S.: the following is the input I used for SAPT0/TZP-DKH:

memory 30 gb

molecule dimer {
 -1 1
 I          0.000          0.00000        3.069
 --
 0 1
 I          0.00000        0.00000        0.00000
 I          0.000          0.0000        -3.069
 units angstrom
}

set globals {
 maxiter 500
 relativistic DKH
}

basis {
 assign TZP-DKH
}

df_basis_sapt {
 assign weigand-jfit
}

df_basis_scf {
 assign weigand-jfit
}

energy('ssapt0')
2

You really should not use a JFIT set for DF-SAPT, you need an appropriate “-RI” basis (in psi4 language). The reason is that you need to fit arbitrary products of basis functions, which is much more sensitive than just Coulomb (“J”) integrals to approximate the density.
Also I think the Weigand set is meant for ECP calculations for most heavy atoms since they follow the Turbomole def2-XVP basis.

You need to turn off DF (possible?) or find a correct fitting basis set.

3

Thanks for the reply.

I couldn’t find much information in the manual for this, so here’s what I did:

  1. simply commenting the use of weigand-jfit gave me the following error:

    QcdbException BasisSetNotFound: BasisSet::construct: Unable to find a basis set for atom 1 for role JKFIT among:
  Shell Entries: ['I']
  Basis Sets: ['def2-qzvpp-jkfit']

  2. set scf_type to pk (since the default is df for SAPT), which asked me now for a new type of basis set (good):

    QcdbException BasisSetNotFound: BasisSet::construct: Unable to find a basis set for atom 1 for role RIFIT among:
  Shell Entries: ['I']
  Basis Sets: ['def2-qzvpp-ri']

  3. I’ve tried SAPT0 instead of SAPT2. No luck.

Due to the lack of ECPs (am I wrong?) in Psi4, I am stuck with all-electron, relativistic corrected basis sets. Furthermore, all-electron auxiliary basis sets for iodine seem scarce. And as it looks like it is easy to disable JKFIT, but I don’t know how to disable RIFIT.

The question now is then, how can I run a SAPT calculation without using a “-RI” basis set?

Thanks again and all the best,
Felipe

P.S.: the following input reproduces the error #2 above:

memory 30 gb

molecule dimer {
 -1 1
 I          0.000          0.00000        3.069
 --
 0 1
 I          0.00000        0.00000        0.00000
 I          0.000          0.0000        -3.069
 units angstrom
}

set globals {
 maxiter 500
 relativistic DKH
 scf_type pk
 df_scf_guess false
}

basis {
 assign TZP-DKH
}

energy('sapt2-ct')
4

I would need to think about it a bit more, but DKH could imbalance the Elst equations as they rely on large cancelations between the one-electron (attractive) and two-electron (repulsive) terms. As DKH is only applied to the one-electron part this could be an issue.

For Iodine there are no fitting basis set that are available and valid without ECP’s (at least AFAIK). You can compute the Elst, Exch, and Ind terms with PK using scf_type, but Disp requires density fitting. You can set SAPT0_E20DISP False to skip this step. It should be noted that the RIFIT basis will always be be built even if it is not used.

Work is under way to implement ECP’s, but it will probably be a month or two yet.

5

The potential DKH issue interests me as well, so I took an all-electron, relativist basis sets from the ORCA
program and used their new automatic generation of fitting basis sets (10.1021/acs.jctc.6b01041) feature.

Could you download these basis sets and try them? I don’t have DKH compiled.

Without DKH and scf_type pk I get at least the ballpark correct:
Total SAPT0 -64.38085856

You can just put them into the directory where you run the calculation and specify

set globals {
 basis oldTZVP
 DF_basis_sapt oldTZVP-RI
}
6

This sounds great! Here’s what I did:

First, the input I used with DKH was the following:

memory 30 gb

molecule dimer {
 -1 1
 I          0.000          0.00000        3.069
 --
 0 1
 I          0.00000        0.00000        0.00000
 I          0.000          0.0000        -3.069
 units angstrom
}

set globals {
 maxiter 500
 relativistic DKH
 scf_type pk
 df_scf_guess false
}

basis {
 assign oldTZVP
}

df_basis_sapt {
 assign oldTZVP-RI
}

energy('sapt2-ct')

Results:

# File: sapt2-ct-autoaux/i-ii0.out
# SAPT2 (kcal/mol)

Bonding:         -46.20
                -------
Exchange:         65.49
Electrostatic:   -51.68 (46.27%)
Induction:       -51.58 (46.18%)
Dispersion:       -8.42 ( 7.54%)

Then without DKH (just commented the relevant line in the input above and rerun):

# File: sapt2-ct-autoaux/i-ii0-nodkh.out
# SAPT2 (kcal/mol)

Bonding:         -45.29
                -------
Exchange:         71.00
Electrostatic:   -54.48 (46.85%)
Induction:       -52.98 (45.56%)
Dispersion:       -8.82 ( 7.58%)

As it seems, scf_type pk does change values (you said you got around -64.4 for SAPT0).
But DKH/no-DKH does not seem to change things much in this case (difference found = 0.91 kcal/mol).

Thanks,
Felipe

P.S.: I am now interested in redoing the basis set generation in ORCA. I generated the auxiliary basis sets using the following input in ORCA 4:

! BP86 D3 old-DKH-TZVPP DKH AutoAux LargePrint

Does the above look correct? The following was obtained for I (exponents look the same as yours):

Basis:

# Basis set for element : I
NewGTO I
S 8
  1  1386039.6978000000      0.0031588380
  2  207780.6500700000      0.0088650754

Auxiliary:

# Auxiliary basis set for element : I
NewAuxGTO I
S 1
  1   78957.7137313359      1.0000000000
S 1
  1   43865.3965174088      1.0000000000
S 1
  1   24369.6647318938      1.0000000000
7

I just used “old-TZVP” for a quick test, but your choice looks better. You can use !PRINTBASIS instead of Largeprint.
You can use my python script for the ORCA-PSI4 conversion.

orca2psi4.py (4.1 KB)

run something like: ./orca2psi4.py orca.out newbasis.gbs orca

There is a last problem with DKH! The relativistic basis set should be decontracted for the DKH transformations, and then re-contracted for the SCF (and everything else). That is how ORCA does it at least. I think this does not work correctly in a fully automated fashion. Part of the discussion is here: https://github.com/psi4/psi4/issues/662

You could try something like:
set {
basis_relativistic mybasis-decon
}

If you can reproduce ORCA results for HF (or get very close), that would be great.

8

Interesting! And thanks for your script, I ended up coding something similar before reading your reply.

Here is my attempt at producing the same energy values along both codes.
The following input for ORCA 4…

! HF old-DKH-TZVPP DKH AutoAux LargePrint

*xyz 0 1
 H 0.00 0.00 0.00
 F 0.00 0.00 0.92
*

…gave the following values of energy:

Total Energy:         -100.15187130 Eh
Nuclear Repulsion:       5.17673356 Eh

The following input for Psi4 (using the basis sets provided by ORCA 4)…

molecule {
 H 0.00 0.00 0.00
 F 0.00 0.00 0.92
}

set {
 scf_type pk
 basis hf
 relativistic dkh
}

energy('hf')

…produced the following (I am showing all significant digits each package gives me):

Total Energy:         -100.1460133852145304 Eh
Nuclear Repulsion:       5.1767335622934798 Eh

Not good!

I then tried to decontract the basis in ORCA (with Decontract). I got the following:

Total Energy:         -100.15275457 Eh
Nuclear Repulsion:       5.17673356 Eh

Using the decontracted basis set (from the last run with ORCA) generated the following with Psi4:

Total Energy:         -100.1527545735322491 Eh
Nuclear Repulsion:       5.1767335622934798 Eh

Pretty close! I will try to apply the same procedure to SAPT but I think the issue is settled!
Thank you again!

9

uhm…I meant Hartree-Fock not hydrogen fluoride :grinning:

So to sum up:
If you use fully decontracted basis sets in both ORCA and PSI4 you get agreement?

This is progress, but usually not what you want (expensive!). The basis set should only be decontracted during the DKH transformation, which the keyword basis_relativistic is meant for (If I understand correctly). If PSI4 automatically decontracts the basis I do not know. Maybe you need to set the basis_relativistic keyword explicitly.

set {
 scf_type pk
 basis MyContractedBasis
 basis_relativistic MyDecontractedBasis
 #or let PSI4 do it automatically with:
 #  basis_relativistic  MyContractedBasis-decon
 relativistic dkh
}
10

It turned out that hydrogen fluoride was the example in the link you’ve
sent :slight_smile:

Yes, agreement was reached with a fully contracted basis set.

I couldn’t make it to work with basis_relativistic (unknown keyword,
maybe because I’m using the stable version of Psi4?). And since all my
systems are small, I thought about using it this way.

11

Could be a version thing yes, I don’t know.

As far as I understand it, using a fully decontracted basis is a valid approach.
Good luck with your project!

12

Thank you, in particular for your interest and time!

All the best,
Felipe