Different energies and much longer time for DKH2 comparing with molpro

Hi,

I am comparing DKH2 calculations from PSI4 with Molpro.
For the following inputs, I got from PSI4 (basis set obtained from https://www.basissetexchange.org/) and from Molpro energies above the convergence threshold:

molecule clbr {
0 1
Cl 0. 0. 0.
Br 0. 0. 2.136
}

basis mb {
[cc-pVTZ-DK-def]
spherical


Br 0
S 20 1.00
10639000.0000000 0.0001832
1593400.0000000 0.0005094
362610.0000000 0.0013353
102700.0000000 0.0032222
33501.0000000 0.0076707
12093.0000000 0.0182439
4715.9000000 0.0433826
1955.6000000 0.0991426
852.6100000 0.2010706
387.6700000 0.3191212
182.6800000 0.3179374
88.2450000 0.1406661
39.2630000 0.0138611
19.2340000 -0.0008414
9.4057000 0.0003982
4.1601000 -0.0002226
1.8995000 0.0000678
0.6047200 -0.0000338
0.3011400 0.0000221
0.1251500 -0.0000052
S 20 1.00
10639000.0000000 -0.0000591
1593400.0000000 -0.0001646
362610.0000000 -0.0004325
102700.0000000 -0.0010466
33501.0000000 -0.0025103
12093.0000000 -0.0060229
4715.9000000 -0.0146734
1955.6000000 -0.0347901
852.6100000 -0.0773387
387.6700000 -0.1442475
182.6800000 -0.1985461
88.2450000 -0.0726261
39.2630000 0.3809663
19.2340000 0.5761470
9.4057000 0.2033629
4.1601000 0.0112299
1.8995000 0.0010430
0.6047200 -0.0001708
0.3011400 0.0000729
0.1251500 -0.0000318
S 20 1.00
10639000.0000000 0.0000236
1593400.0000000 0.0000657
362610.0000000 0.0001726
102700.0000000 0.0004176
33501.0000000 0.0010031
12093.0000000 0.0024075
4715.9000000 0.0058916
1955.6000000 0.0140102
852.6100000 0.0316098
387.6700000 0.0600877
182.6800000 0.0872309
88.2450000 0.0324842
39.2630000 -0.2258144
19.2340000 -0.5093521
9.4057000 -0.1106266
4.1601000 0.7018218
1.8995000 0.5477988
0.6047200 0.0393571
0.3011400 -0.0089968
0.1251500 0.0019468
S 1 1.00
0.6047200 1.0000000
S 20 1.00
10639000.0000000 -0.0000072
1593400.0000000 -0.0000202
362610.0000000 -0.0000530
102700.0000000 -0.0001284
33501.0000000 -0.0003083
12093.0000000 -0.0007407
4715.9000000 -0.0018109
1955.6000000 -0.0043190
852.6100000 -0.0097329
387.6700000 -0.0186483
182.6800000 -0.0270679
88.2450000 -0.0105140
39.2630000 0.0747565
19.2340000 0.1753422
9.4057000 0.0399153
4.1601000 -0.3301097
1.8995000 -0.4463003
0.6047200 0.2957042
0.3011400 0.6810696
0.1251500 0.2785627
S 1 1.00
0.1251500 1.0000000
P 13 1.00
8676.5000000 0.0009606
2055.9000000 0.0050551
666.2300000 0.0234154
253.1000000 0.0841414
106.1200000 0.2229616
47.2420000 0.3880657
21.8250000 0.3534548
9.9684000 0.1085321
4.5171000 0.0042339
1.9982000 0.0016954
0.7098800 -0.0004284
0.2814500 0.0001980
0.1020400 -0.0000690
P 13 1.00
8676.5000000 -0.0003881
2055.9000000 -0.0020557
666.2300000 -0.0096085
253.1000000 -0.0355099
106.1200000 -0.0981287
47.2420000 -0.1828576
21.8250000 -0.1498313
9.9684000 0.1959655
4.5171000 0.5433698
1.9982000 0.3792448
0.7098800 0.0434121
0.2814500 -0.0040330
0.1020400 0.0017114
P 1 1.00
0.7098800 1.0000000
P 13 1.00
8676.5000000 0.0001001
2055.9000000 0.0005330
666.2300000 0.0024808
253.1000000 0.0092595
106.1200000 0.0255286
47.2420000 0.0486125
21.8250000 0.0376864
9.9684000 -0.0631036
4.5171000 -0.1878988
1.9982000 -0.1319307
0.7098800 0.2976662
0.2814500 0.5723946
0.1020400 0.3063896
P 1 1.00
0.1020400 1.0000000
D 9 1.00
403.8300000 0.0016248
121.1700000 0.0130068
46.3450000 0.0582881
19.7210000 0.1698577
8.8624000 0.3172628
3.9962000 0.3831152
1.7636000 0.2764135
0.7061900 0.0753432
0.2639000 0.0020696
D 1 1.00
0.7061900 1.0000000
D 1 1.00
0.2639000 1.0000000
F 1 1.00
0.5515000 1.0000000


Cl 0
S 15 1.00
456100.0000000 0.0002645
68330.0000000 0.0008905
15550.0000000 0.0029934
4405.0000000 0.0098871
1439.0000000 0.0315606
520.4000000 0.0903711
203.1000000 0.2134406
83.9600000 0.3646772
36.2000000 0.3377387
15.8300000 0.1006810
6.3340000 0.0030732
2.6940000 0.0010373
0.9768000 -0.0003695
0.4313000 0.0001520
0.1625000 -0.0000502
S 15 1.00
456100.0000000 -0.0000747
68330.0000000 -0.0002515
15550.0000000 -0.0008492
4405.0000000 -0.0028128
1439.0000000 -0.0091519
520.4000000 -0.0269903
203.1000000 -0.0699330
83.9600000 -0.1426906
36.2000000 -0.1978853
15.8300000 -0.0143283
6.3340000 0.5036507
2.6940000 0.5584234
0.9768000 0.0768656
0.4313000 -0.0080500
0.1625000 0.0022404
S 1 1.00
0.9768000 1.0000000
S 15 1.00
456100.0000000 0.0000226
68330.0000000 0.0000761
15550.0000000 0.0002572
4405.0000000 0.0008515
1439.0000000 0.0027786
520.4000000 0.0082052
203.1000000 0.0214933
83.9600000 0.0445276
36.2000000 0.0646680
15.8300000 0.0041911
6.3340000 -0.2092752
2.6940000 -0.4033247
0.9768000 0.0827071
0.4313000 0.7256082
0.1625000 0.3896193
S 1 1.00
0.1625000 1.0000000
P 9 1.00
663.3000000 0.0026741
156.8000000 0.0197715
49.9800000 0.0894991
18.4200000 0.2567502
7.2400000 0.4366374
2.9220000 0.3493206
1.0220000 0.0583589
0.3818000 -0.0045545
0.1301000 0.0022637
P 1 1.00
1.0220000 1.0000000
P 9 1.00
663.3000000 -0.0007253
156.8000000 -0.0053474
49.9800000 -0.0249761
18.4200000 -0.0730005
7.2400000 -0.1338923
2.9220000 -0.0940334
1.0220000 0.2627700
0.3818000 0.5640464
0.1301000 0.3413419
P 1 1.00
0.1301000 1.0000000
D 1 1.00
1.0460000 1.0000000
D 1 1.00
0.3440000 1.0000000
F 1 1.00
0.7060000 1.0000000


assign Cl cc-pVTZ-DK-def
assign Br cc-pVTZ-DK-def
}

basis mb2 {
[cc-pVTZ-DK-def-d]
spherical


Cl 0
S 1 1.00
456100.0000000 1.00000000D+00
S 1 1.00
68330.0000000 1.00000000D+00
S 1 1.00
15550.0000000 1.00000000D+00
S 1 1.00
4405.0000000 1.00000000D+00
S 1 1.00
1439.0000000 1.00000000D+00
S 1 1.00
520.4000000 1.00000000D+00
S 1 1.00
203.1000000 1.00000000D+00
S 1 1.00
83.9600000 1.00000000D+00
S 1 1.00
36.2000000 1.00000000D+00
S 1 1.00
15.8300000 1.00000000D+00
S 1 1.00
6.3340000 1.00000000D+00
S 1 1.00
2.6940000 1.00000000D+00
S 1 1.00
0.9768000 1.00000000D+00
S 1 1.00
0.4313000 1.00000000D+00
S 1 1.00
0.1625000 1.00000000D+00
P 1 1.00
663.3000000 1.00000000D+00
P 1 1.00
156.8000000 1.00000000D+00
P 1 1.00
49.9800000 1.00000000D+00
P 1 1.00
18.4200000 1.00000000D+00
P 1 1.00
7.2400000 1.00000000D+00
P 1 1.00
2.9220000 1.00000000D+00
P 1 1.00
1.0220000 1.00000000D+00
P 1 1.00
0.3818000 1.00000000D+00
P 1 1.00
0.1301000 1.00000000D+00
D 1 1.00
1.0460000 1.00000000D+00
D 1 1.00
0.3440000 1.00000000D+00
F 1 1.00
0.7060000 1.00000000D+00


Br 0
S 1 1.00
10639000.0000000 1.00000000D+00
S 1 1.00
1593400.0000000 1.00000000D+00
S 1 1.00
362610.0000000 1.00000000D+00
S 1 1.00
102700.0000000 1.00000000D+00
S 1 1.00
33501.0000000 1.00000000D+00
S 1 1.00
12093.0000000 1.00000000D+00
S 1 1.00
4715.9000000 1.00000000D+00
S 1 1.00
1955.6000000 1.00000000D+00
S 1 1.00
852.6100000 1.00000000D+00
S 1 1.00
387.6700000 1.00000000D+00
S 1 1.00
182.6800000 1.00000000D+00
S 1 1.00
88.2450000 1.00000000D+00
S 1 1.00
39.2630000 1.00000000D+00
S 1 1.00
19.2340000 1.00000000D+00
S 1 1.00
9.4057000 1.00000000D+00
S 1 1.00
4.1601000 1.00000000D+00
S 1 1.00
1.8995000 1.00000000D+00
S 1 1.00
0.6047200 1.00000000D+00
S 1 1.00
0.3011400 1.00000000D+00
S 1 1.00
0.1251500 1.00000000D+00
P 1 1.00
8676.5000000 1.00000000D+00
P 1 1.00
2055.9000000 1.00000000D+00
P 1 1.00
666.2300000 1.00000000D+00
P 1 1.00
253.1000000 1.00000000D+00
P 1 1.00
106.1200000 1.00000000D+00
P 1 1.00
47.2420000 1.00000000D+00
P 1 1.00
21.8250000 1.00000000D+00
P 1 1.00
9.9684000 1.00000000D+00
P 1 1.00
4.5171000 1.00000000D+00
P 1 1.00
1.9982000 1.00000000D+00
P 1 1.00
0.7098800 1.00000000D+00
P 1 1.00
0.2814500 1.00000000D+00
P 1 1.00
0.1020400 1.00000000D+00
D 1 1.00
403.8300000 1.00000000D+00
D 1 1.00
121.1700000 1.00000000D+00
D 1 1.00
46.3450000 1.00000000D+00
D 1 1.00
19.7210000 1.00000000D+00
D 1 1.00
8.8624000 1.00000000D+00
D 1 1.00
3.9962000 1.00000000D+00
D 1 1.00
1.7636000 1.00000000D+00
D 1 1.00
0.7061900 1.00000000D+00
D 1 1.00
0.2639000 1.00000000D+00
F 1 1.00
0.5515000 1.00000000D+00


assign Cl cc-pVTZ-DK-def-d
assign Br cc-pVTZ-DK-def-d
}

set {
reference rhf
relativistic dkh
dkh_order 2
scf_type pk
guess gwh
basis mb
basis_relativistic mb2
D_convergence 1e-10
E_convergence 1e-10
}

energy(โ€˜hfโ€™)

memory,120,m
gprint,basis
dkho=2
angstrom
geometry={
Cl 0. 0. 0.
Br 0. 0. 2.136
}

wf,charge=0,spin=0

gthresh,energy=1.d-10

basis,cl=cc-pVTZ-DK, br=cc-pVTZ-DK

{hf
}
{ccsd(t)}

psi-4 1.7 -3065.2229486335891124
molpro 2022 -3065.222943870092

the difference (5*10^-6) is above the convergence threshold (10^(-10))

My questions are:

  1. What is the reason for this energy difference above the convergence threshold?
  2. What is the reason that PSI4 spends much more time than molpro (about a factor of 20, for more electron systems, the difference is larger)? Is that due to handling of general contracted basis? Density-fitting would be faster. But, will the default, def2-universal-jkfit, suitable for basis set optimized for DKH? In Auxiliary Basis Sets, there is not an auxiliary basis for DK basis.
  3. Why in PSI4 one needs to prepare an uncontracted auxiliary basis for DKH as in Scalar relativistic Hamiltonians?

Thank you very much

P.S. seems four asterisks after a basis set disappeared in this post

  1. Is it the same basis set? It should be enough to specify cc-pVTZ-DK in psi4 and psi4 automatically decontracts the basis for the relativistic correction. Molpro should be doing the same.

  2. Most of it is likely due to the generally contacted basis for which psi4 has no special integral engine like molpro

  3. To increase the accuracy in the core region, also see 1)

for this input psi4/input.dat at master ยท psi4/psi4 ยท GitHub we match with molpro

1 Like

Thank you very much for your reply. I modified the GitHub example to increase the convergence threshold as following

#! DKH comparison Psi4 vs. Molpro from Problems enabling DKH - #3 by Kirk

memory 450 mb

molecule Ne {
0 1
Ne
}

set {
reference rhf
basis cc-pvtz-dk
#rel_basis cc-pvtz-dk-decon
relativistic dkh
dkh_order 2
print 2
scf_type pk
D_convergence 1e-10
E_convergence 1e-10
}

e = energy(โ€˜scfโ€™)

print_variables()

and the molpro input is

memory,120,m
gprint,basis
dkho=2
angstrom
geometry={angstrom;ne}

wf,charge=0,spin=0
!basis=cc-pVDZ-DK

gthresh,energy=1.d-10
gthresh,orbital=1.d-10

basis,ne=cc-pVTZ-DK

{hf
}

I got -128.668916123700 from psi4 1.7 and -128.668916104512 from molpro 2022, about 1.92e-08 energy difference, looks excellent!

Thank you once again. For clbr, the input
molecule clbr {
0 1
Cl 0. 0. 0.
Br 0. 0. 2.136
}

set {
reference rhf
basis cc-pvtz-dk
relativistic dkh
dkh_order 2
guess gwh
D_convergence 1e-10
E_convergence 1e-10
}

energy(โ€˜hfโ€™)

works and gives the energy -3065.2229486336927948, which differs molpro as ~ 5*10^-6. Perhaps would be the difference of the basis set as you mentioned

https://www.molpro.net/info/basis.php?search=1&element=Br&l=s&basis=cc-pVTZ-DK&print=1
in psi4, Br, contraction coefficient of s function, 0.0001832 , molpro Br, 0.000183

In addition, if I do not use pk type scf, the default will use def2-universal-jkfit.gbs for density fitting. It is much faster. May I know how accurate def2-universal-jkfit is with DK type basis set?

The def2-JK basis sets expect ECPs for heavy elements. Though I think Br is still fine. Missing core auxiliarly functions would be bad.

Maybe you can try to compare other energy components from the output?

1 Like

I am sorry for spamming one more reply. To look into the effect on the last few digits in basis set parameters

I use the https://www.basissetexchange.org/ basis set as molpro input (the same for psi4) in the first post.

memory,120,m
gprint,basis
dkho=2
angstrom
geometry={
Cl 0. 0. 0.
Br 0. 0. 2.136
}

wf,charge=0,spin=0
!basis=cc-pVDZ-DK

gthresh,energy=1.d-10

!basis,cl=cc-pVTZ-DK, br=cc-pVTZ-DK

spherical
basis={
!
! chlorine (15s,9p,2d,1f) โ†’ [5s,4p,2d,1f]
s, CL , 456100.0000000, 68330.0000000, 15550.0000000, 4405.0000000, 1439.0000000, 520.4000000, 203.1000000, 83.9600000, 36.2000000, 15.8300000, 6.3340000, 2.6940000, 0.9768000, 0.4313000, 0.1625000
c, 1.15, 0.0002645, 0.0008905, 0.0029934, 0.0098871, 0.0315606, 0.0903711, 0.2134406, 0.3646772, 0.3377387, 0.1006810, 0.0030732, 0.0010373, -0.0003695, 0.0001520, -0.0000502
c, 1.15, -0.0000747, -0.0002515, -0.0008492, -0.0028128, -0.0091519, -0.0269903, -0.0699330, -0.1426906, -0.1978853, -0.0143283, 0.5036507, 0.5584234, 0.0768656, -0.0080500, 0.0022404
c, 13.13, 1.0000000
c, 1.15, 0.0000226, 0.0000761, 0.0002572, 0.0008515, 0.0027786, 0.0082052, 0.0214933, 0.0445276, 0.0646680, 0.0041911, -0.2092752, -0.4033247, 0.0827071, 0.7256082, 0.3896193
c, 15.15, 1.0000000
p, CL , 663.3000000, 156.8000000, 49.9800000, 18.4200000, 7.2400000, 2.9220000, 1.0220000, 0.3818000, 0.1301000
c, 1.9, 0.0026741, 0.0197715, 0.0894991, 0.2567502, 0.4366374, 0.3493206, 0.0583589, -0.0045545, 0.0022637
c, 7.7, 1.0000000
c, 1.9, -0.0007253, -0.0053474, -0.0249761, -0.0730005, -0.1338923, -0.0940334, 0.2627700, 0.5640464, 0.3413419
c, 9.9, 1.0000000
d, CL , 1.0460000, 0.3440000
c, 1.1, 1.0000000
c, 2.2, 1.0000000
f, CL , 0.7060000
c, 1.1, 1.0000000
!
! bromine (20s,13p,9d,1f) โ†’ [6s,5p,3d,1f]
s, BR , 10639000.0000000, 1593400.0000000, 362610.0000000, 102700.0000000, 33501.0000000, 12093.0000000, 4715.9000000, 1955.6000000, 852.6100000, 387.6700000, 182.6800000, 88.2450000, 39.2630000, 19.2340000, 9.4057000, 4.1601000, 1.8995000, 0.6047200, 0.3011400, 0.1251500
c, 1.20, 0.0001832, 0.0005094, 0.0013353, 0.0032222, 0.0076707, 0.0182439, 0.0433826, 0.0991426, 0.2010706, 0.3191212, 0.3179374, 0.1406661, 0.0138611, -0.0008414, 0.0003982, -0.0002226, 0.0000678, -0.0000338, 0.0000221, -0.0000052
c, 1.20, -0.0000591, -0.0001646, -0.0004325, -0.0010466, -0.0025103, -0.0060229, -0.0146734, -0.0347901, -0.0773387, -0.1442475, -0.1985461, -0.0726261, 0.3809663, 0.5761470, 0.2033629, 0.0112299, 0.0010430, -0.0001708, 0.0000729, -0.0000318
c, 1.20, 0.0000236, 0.0000657, 0.0001726, 0.0004176, 0.0010031, 0.0024075, 0.0058916, 0.0140102, 0.0316098, 0.0600877, 0.0872309, 0.0324842, -0.2258144, -0.5093521, -0.1106266, 0.7018218, 0.5477988, 0.0393571, -0.0089968, 0.0019468
c, 18.18, 1.0000000
c, 1.20, -0.0000072, -0.0000202, -0.0000530, -0.0001284, -0.0003083, -0.0007407, -0.0018109, -0.0043190, -0.0097329, -0.0186483, -0.0270679, -0.0105140, 0.0747565, 0.1753422, 0.0399153, -0.3301097, -0.4463003, 0.2957042, 0.6810696, 0.2785627
c, 20.20, 1.0000000
p, BR , 8676.5000000, 2055.9000000, 666.2300000, 253.1000000, 106.1200000, 47.2420000, 21.8250000, 9.9684000, 4.5171000, 1.9982000, 0.7098800, 0.2814500, 0.1020400
c, 1.13, 0.0009606, 0.0050551, 0.0234154, 0.0841414, 0.2229616, 0.3880657, 0.3534548, 0.1085321, 0.0042339, 0.0016954, -0.0004284, 0.0001980, -0.0000690
c, 1.13, -0.0003881, -0.0020557, -0.0096085, -0.0355099, -0.0981287, -0.1828576, -0.1498313, 0.1959655, 0.5433698, 0.3792448, 0.0434121, -0.0040330, 0.0017114
c, 11.11, 1.0000000
c, 1.13, 0.0001001, 0.0005330, 0.0024808, 0.0092595, 0.0255286, 0.0486125, 0.0376864, -0.0631036, -0.1878988, -0.1319307, 0.2976662, 0.5723946, 0.3063896
c, 13.13, 1.0000000
d, BR , 403.8300000, 121.1700000, 46.3450000, 19.7210000, 8.8624000, 3.9962000, 1.7636000, 0.7061900, 0.2639000
c, 1.9, 0.0016248, 0.0130068, 0.0582881, 0.1698577, 0.3172628, 0.3831152, 0.2764135, 0.0753432, 0.0020696
c, 8.8, 1.0000000
c, 9.9, 1.0000000
f, BR , 0.5515000
c, 1.1, 1.0000000
}

{hf
}
{ccsd(t)}

The molpro output includes
!RHF STATE 1.1 Energy -3065.222943870093
still differs with PSI4 -3065.2229486335891124, by ~5*10^-6. I think the basis sets are the same now, since both are from https://www.basissetexchange.org/ (I can print out psi4 build-in basis, then convert to molpro format, that will be an additional check)