Problem with tungsten - DFT and HF (SCF) - unreasonable energies

I am trying to do some DFT calculations involving a tungsten complex, specifically the addition to H2 to trans-(iPr3P)2W(CO)3 – this is reaction #7 in the MOR41 database of Dohm et al. (JCTC, 2018, 14, 2596 - for which I have prepared a PSI4 database). I would like to be using some of the newer functionals available through libxc.

I am using Psi4 version 1.2.1 installed though miniconda.

I started with B3LYP for testing purposes. With the def2-SVP basis set I get reasonable results (see Table below). When I switch to the larger def2-TZVP basis set, the error is quite significant and repulsive, and with def2-TZVPP the error is completely ridiculous. I double-checked the basis set included against that in the EMSL basis set exchange and confirmed that they are the same. I did notice that Gaussian16 has a different ECP (notably lmax=5 instead of lmax=3 in PSI4) but using this ECP in PSI4 (I edited the basis set file) with the def2-TZVP basis set did not help. I also tried using Hartree-Fock but got the same problem. To me, this naïvely suggests a problem with the integral code; the problem with HF indicates that this cannot be a grid issue. However, what is weird is that I have only seen this problem with W - MOR41 includes many of the transition metals, but the problem only seems to occur with W.

The energies obtained from the calculations and the corresponding binding energies are:


                	reactant	 product	 H2	         ∆E
B3LYP/SVP	        -1800.179680	-1801.373734	-1.173627	-12.82
B3LYP/TZVP	        -1802.065479	-1803.088290	-1.179675	 98.43
B3LYP/TZVPP	        -1806.220179	-1805.715773	-1.180027      1057.00
B3LYP/TZVP(run with G16)-1801.629069	-1802.827269	-1.180004	-11.42
B3LYP/TZVP(G16ECP)      -1792.608765	-1793.612183	-1.179675	110.60
HF/TZVP	                -1792.860925	-1793.757451	-1.132519	148.09
HF/SVP	                -1791.168218	-1792.310048	-1.128942	 -8.09
HF/TZVPP	        -1796.679676	-1796.679676	-1.132944	710.93			0.00
				
reference (from Dohm et al.)		                                -16.2

Any help on resolving this issue would be appreciated. Attached are my input files (B3LYP/def2-TZVP).

Thanks in advance for your assistance.

Sincerely yours,

Mark.


Input file: reactant


# This is a psi4 input file auto-generated from the database() wrapper.

core.print_out('\n')
p4util.banner(' Database MOR41 Computation: Reagent MOR41-ED07-reagent \n    ')
core.print_out('\n')


molecule dbmol {
    units Angstrom
    no_com
    no_reorient
    0 1
    W                0.090640800000     0.066260000000    -0.088000300000
    C                0.106952100000    -1.928814200000    -0.466549400000
    O                0.096133700000    -3.065375600000    -0.740415000000
    C                0.028274400000     0.379980900000    -1.997828000000
    O                0.028029300000     0.578417300000    -3.159275700000
    C               -0.093307900000     2.080442700000     0.142785900000
    O               -0.257037700000     3.235495700000     0.206221600000
    P               -2.383484100000    -0.093746500000     0.193528000000
    C               -2.920792500000    -0.437281000000     1.956908800000
    H               -4.016616300000    -0.455293800000     1.975816700000
    C               -3.316389800000     1.491151500000    -0.179503100000
    H               -2.732679400000     2.213265400000     0.405760500000
    C               -3.171527500000    -1.529706200000    -0.720131900000
    H               -2.611149500000    -2.369346900000    -0.287224300000
    C                3.519099200000    -1.266466700000    -0.796919500000
    P                2.532536100000    -0.033802100000     0.199724100000
    H                3.444524000000    -0.211399300000     2.511588000000
    H                3.107501500000    -2.219452200000    -0.441725700000
    C                3.585414000000     1.508545400000     0.239418500000
    H                4.620761700000     1.185809300000     0.407929600000
    C               -3.188174200000     1.901699700000    -1.653699400000
    H               -3.484484300000     2.950588100000    -1.768427500000
    H               -3.845391300000     1.299978000000    -2.288545900000
    H               -2.164782300000     1.791235700000    -2.018112000000
    C               -4.776111600000     1.586022800000     0.287039300000
    H               -5.430390800000     0.926602700000    -0.288175000000
    H               -5.131429100000     2.613182000000     0.138362200000
    H               -4.897503100000     1.346733800000     1.347270300000
    C               -4.665911000000    -1.776671700000    -0.470122300000
    H               -4.929291900000    -1.762487500000     0.591439900000
    H               -4.941115900000    -2.760855300000    -0.868429500000
    H               -5.281061700000    -1.033761700000    -0.985240800000
    C               -2.870200200000    -1.526003600000    -2.225596100000
    H               -3.144218800000    -2.498142700000    -2.651782600000
    H               -1.812675800000    -1.348054600000    -2.425012500000
    H               -3.448170400000    -0.758446700000    -2.747043500000
    C               -2.408396600000    -1.803345600000     2.432906700000
    H               -2.872536800000    -2.627943400000     1.884708800000
    H               -2.631951200000    -1.938346500000     3.497793800000
    H               -1.323907900000    -1.887034700000     2.300676200000
    C               -2.438924400000     0.690855900000     2.879924400000
    H               -1.350131300000     0.797763100000     2.830651500000
    H               -2.713842600000     0.472510200000     3.918633000000
    H               -2.876995800000     1.656057300000     2.609917300000
    C                5.027356300000    -1.245450500000    -0.518440400000
    H                5.252864000000    -1.270709500000     0.553401900000
    H                5.493823500000    -0.349511700000    -0.941425500000
    H                5.505229900000    -2.115642400000    -0.983606400000
    C                3.211772300000    -1.184836700000    -2.299117800000
    H                3.611414100000    -0.271091500000    -2.747925100000
    C                1.275062600000    -0.181216800000     2.574089300000
    H                2.135619700000    -1.206862800000    -2.482717000000
    H                3.670616200000    -2.038173500000    -2.812060800000
    C                2.705816400000    -2.193295000000     2.086512200000
    H                0.415776100000    -0.627664700000     2.006295600000
    H                1.117135600000    -0.557046900000     3.593282600000
    H                1.188968300000     0.905789500000     2.590191600000
    C                2.605867800000    -0.667968300000     1.970815600000
    H                2.602648800000    -2.490168400000     3.137199300000
    H                3.666061900000    -2.571011900000     1.726042600000
    H                1.908387600000    -2.684812300000     1.519271800000
    C                3.176687800000     2.445829000000     1.383263000000
    H                2.128772900000     2.744492100000     1.290794200000
    H                3.786333900000     3.356035300000     1.342193600000
    H                3.326985000000     1.990449400000     2.366969500000
    C                3.486197100000     2.235599400000    -1.109636000000
    H                2.440861800000     2.462784300000    -1.343426600000
    H                3.892085500000     1.637683200000    -1.929626800000
    H                4.044777600000     3.178009400000    -1.068198700000


}

core.set_memory_bytes(15000000000)

core.set_global_option('BASIS', 'DEF2-TZVP')
core.set_global_option('GUESS', 'GWH')
core.set_global_option('MAXITER', 1000)

core.set_global_option('WRITER_FILE_LABEL', 'MOR41-ED07-reagent')

pickle_kw = ("""(dp0
S'db_func'
p1
cpsi4.driver.driver
energy
p2
sS'db_name'
p3
S'MOR41'
p4
sS'db_mode'
p5
S'sow'
p6
sS'name'
p7
S'B3LYP'
p8
s.""")

kwargs = pickle.loads(pickle_kw)
electronic_energy = energy(**kwargs)

core.print_variables()
core.print_out('\nDATABASE RESULT: computation 11258 for reagent MOR41-ED07-reagent yields electronic energy %20.12f\n' % (electronic_energy))

core.set_variable('NATOM', dbmol.natom())

input file: product


# This is a psi4 input file auto-generated from the database() wrapper.

core.print_out('\n')
p4util.banner(' Database MOR41 Computation: Reagent MOR41-PR07-reagent \n    ')
core.print_out('\n')


molecule dbmol {
    units Angstrom
    no_com
    no_reorient
    0 1
    W               -0.000070400000    -0.017989600000    -0.146457100000
    C                0.000157200000     1.955377300000    -0.631545200000
    O                0.000411800000     3.068043100000    -0.986288700000
    C                0.000004600000     0.271128000000     1.825131000000
    O                0.000003100000     0.436870000000     2.983755900000
    C                0.000079000000    -2.055167300000    -0.056719200000
    O                0.000369700000    -3.220382300000    -0.103637600000
    H                0.422788400000    -0.287137400000    -1.984721800000
    H               -0.422970800000    -0.287404200000    -1.985319300000
    P               -2.498968200000     0.054091800000    -0.052469100000
    C               -3.398360300000    -0.535031900000    -1.584174500000
    H               -4.468441100000    -0.526892600000    -1.348265500000
    C               -3.356878400000    -0.865573300000     1.355714100000
    H               -3.371868200000    -0.108729700000     2.151443200000
    C               -3.157712400000     1.802287600000     0.144667000000
    H               -2.704252500000     2.317272900000    -0.711479200000
    C                3.158107500000     1.802481500000     0.146809300000
    P                2.499417800000     0.054538100000    -0.052471700000
    H                4.468463700000    -0.525230000000    -1.349702400000
    H                2.705790900000     2.318402900000    -0.709381300000
    C                3.358098800000    -0.866609200000     1.354368100000
    H                3.374620800000    -0.110180300000     2.150474500000
    C               -4.805511900000    -1.297741000000     1.073085200000
    H               -5.254728700000    -1.674152000000     1.999881200000
    H               -4.835629400000    -2.112614900000     0.342511600000
    H               -5.434310400000    -0.486090300000     0.702882400000
    C               -2.560531500000    -2.067317500000     1.879344900000
    H               -3.061648100000    -2.475072300000     2.765696100000
    H               -1.542735600000    -1.790111500000     2.156345200000
    H               -2.501752200000    -2.863967000000     1.131819400000
    C               -4.680578800000     1.979997900000     0.065888300000
    H               -4.916432600000     3.048909800000    -0.001308900000
    H               -5.169896500000     1.596146400000     0.966264000000
    H               -5.123495400000     1.487333900000    -0.804742400000
    C               -2.618029800000     2.465249900000     1.419594100000
    H               -2.906610500000     3.522573500000     1.432737600000
    H               -1.529918400000     2.405362800000     1.476662200000
    H               -3.033061600000     1.997312800000     2.319029600000
    C               -3.156120400000     0.402360000000    -2.773590000000
    H               -3.624805500000    -0.011995200000    -3.673909400000
    H               -2.084672200000     0.514366300000    -2.975522500000
    H               -3.570938900000     1.399990900000    -2.605324600000
    C               -2.983408700000    -1.970962500000    -1.930114600000
    H               -3.139122600000    -2.660074400000    -1.095275000000
    H               -1.922494900000    -2.013238900000    -2.196445300000
    H               -3.565671300000    -2.333115000000    -2.785725800000
    C                4.681143700000     1.979452700000     0.069932400000
    H                5.169200400000     1.594107000000     0.970364900000
    H                4.917733700000     3.048301200000     0.004377900000
    H                5.124722600000     1.487610300000    -0.800835900000
    C                2.617304900000     2.464485800000     1.421754200000
    H                3.030128100000     1.994753700000     2.321273700000
    C                2.982560800000    -1.967591100000    -1.933470700000
    H                1.529026700000     2.406081300000     1.476946600000
    H                2.907394000000     3.521374000000     1.436786500000
    C                3.156164200000     0.407188300000    -2.773025600000
    H                1.921746200000    -2.008793700000    -2.200354200000
    H                3.564969200000    -2.328960800000    -2.789307000000
    H                3.137352200000    -2.657957400000    -1.099490900000
    C                3.398322700000    -0.532419600000    -1.585326900000
    H                3.622375400000    -0.006874700000    -3.674767900000
    H                3.573575400000     1.403617700000    -2.603965500000
    H                2.084603000000     0.521976000000    -2.972940400000
    C                4.806061200000    -1.299837300000     1.069941600000
    H                5.255832800000    -1.677434300000     1.995985100000
    H                5.435230800000    -0.488422700000     0.699837500000
    H                4.834740800000    -2.114113900000     0.338643400000
    C                2.561315100000    -2.067960800000     1.878190500000
    H                2.501285900000    -2.864212200000     1.130337600000
    H                1.543963400000    -1.790105100000     2.156196500000
    H                3.062895500000    -2.476510500000     2.763911500000


}

core.set_memory_bytes(15000000000)

core.set_global_option('BASIS', 'DEF2-TZVP')
core.set_global_option('BASIS_GUESS', 'DEF2-SVP')

core.set_global_option('WRITER_FILE_LABEL', 'MOR41-PR07-reagent')

pickle_kw = ("""(dp0
S'db_func'
p1
cpsi4.driver.driver
energy
p2
sS'db_name'
p3
S'MOR41'
p4
sS'db_mode'
p5
S'sow'
p6
sS'name'
p7
S'B3LYP'
p8
s.""")

kwargs = pickle.loads(pickle_kw)
electronic_energy = energy(**kwargs)

core.print_variables()
core.print_out('\nDATABASE RESULT: computation 11258 for reagent MOR41-PR07-reagent yields electronic energy %20.12f\n' % (electronic_energy))

core.set_variable('NATOM', dbmol.natom())

# This is a psi4 input file auto-generated from the database() wrapper.

core.print_out('\n')
p4util.banner(' Database MOR41 Computation: Reagent MOR41-H2-reagent \n    ')
core.print_out('\n')


molecule dbmol {
    units Angstrom
    no_com
    no_reorient
    0 1
    H                0.798743700000     0.046507300000     0.075739400000
    H                1.071087200000     0.408924800000     0.665066700000


}

core.set_memory_bytes(15000000000)

core.set_global_option('BASIS', 'DEF2-TZVP')
core.set_global_option('BASIS_GUESS', 'DEF2-SVP')

core.set_global_option('WRITER_FILE_LABEL', 'MOR41-H2-reagent')

pickle_kw = ("""(dp0
S'db_func'
p1
cpsi4.driver.driver
energy
p2
sS'db_name'
p3
S'MOR41'
p4
sS'db_mode'
p5
S'sow'
p6
sS'name'
p7
S'B3LYP'
p8
s.""")

kwargs = pickle.loads(pickle_kw)
electronic_energy = energy(**kwargs)

core.print_variables()
core.print_out('\nDATABASE RESULT: computation 11258 for reagent MOR41-H2-reagent yields electronic energy %20.12f\n' % (electronic_energy))

core.set_variable('NATOM', dbmol.natom())

I tried this input file on the developer version of Psi, and I can’t even get the product SCF to converge with TZVP, only the small basis guess. std::cout gets hit with over 500 “Failed to converge:” messages indicating the ECP is acting up, and Psi segfaults and dies.

I strongly suspect this problem is with the ECPs. I’m not familiar with the ECP code, but I don’t believe there’s a remedy for this until we rewrite that code completely. I’ll bring this to the attention of some developers who might know more.

(Technical note: I removed the pickled part and replaced **kwargs with b3lyp. The next version of Psi switches from Python2 to Python3, which is going to change the handling of pickled objects.)

Thanks for the quick reply.

If the SCF won’t converge, there is always the option of using “basis_guess def2-svp” since this basis set converges.

If it is an issue with the ECP, then why do I get reasonable results with def2-SVP, which uses the same ECP, and why do I get reasonable results with the other transition metals (inter alia Pd, Pt, Ru, Os, Ir, Rh)?

Thanks for your help. Is there any hope in resolving this issue or should I consider other options?

Sincerely,

Mark.

Your problem was put into github: https://github.com/psi4/psi4/issues/1338
Best to wait a few days to let the code author comment if possible. Especially what those messages mean.

If possible, it would be best to rule out that PSI4’s SCF converged to a wrong state.
The SCF could be hard to converge correctly even if the ECP is correct.

Actually, the option of using basis_guess def2-svp is taken! If you look at the input file you posted, there is already a line to do precisely that.

My biggest concern is that we’re not even seeing the same error. The error you see is nonsense energies for the large basis. The error I’m seeing is a segfault that prevents me from even seeing a large basis energy. While I have grounds to suspect an ECP related bug behind my segfault, as Holger alluded to, your nonsense energies may be more complicated. It’ll take time to investigate this more fully, and the first step on my end is going to be to see if I can reproduce your bug. It would be worth trying the computation in serial and also trying Psi 1.2.1 instead of the dev version I was using. Once I can reproduce your bug, I should be able to run a stability analysis on your Hartree-Fock computation to check for excited electronic states. If those come up clean, then the prime suspect is ECPs again.

I can’t guarantee a fix for the bug in the near future. I’ll need to properly reproduce your bug before I can say much for sure. Fortunately, the Psi developer’s conference is this weekend, and the state of the ECPs is likely to be on the discussion list, so I should have some news for you. If I were in your position, I’d wait a few days for more news from the developers, but be thinking in the back of my mind about alternatives if Psi can’t get this.

I’ve attempted to reproduce this computation with modern Psi and been able to get much more reasonable numbers:

Reactant: -1801.61184030879 E_h
H2: -1.17967496502912 E_h
Product: -1802.80957303765 E_h
Reaction Energy: -11.3 kcal/mol

No special tricks needed - issue closed at long last!