Why am I getting strange results for transition metals in Psi4 1.5?

Hi everyone,

I’ve encountered a significant discrepancy in the atomic energies calculated for several transition metals (W, Re, Ir) using Psi4 at the wb97mv/def2-svpd level of theory. The issue seems to be specific to transition metals, as the calculations for organic molecules (e.g., H, C, O) yield consistent and accurate results when compared to other quantum chemistry packages like Jaguar.

Here’s a summary of the results for W:

Element Molecule Formula Jaguar Atomic Energy (Ha) Psi4 Atomic Energy (Ha) Difference (Ha)
H H2 H2 -0.5770035 -0.577114067 0.000110567
C CH4 CH4 -38.135503 -38.13510729 -0.000395709
O H2O H2O -75.193142 -75.19295249 -0.000189514
W isolated W -66.62167 -66.68160669 0.059936688
W REDNUK C7H4O5W -66.497504 -114.9241744 48.4266704
W FERDIP-CO C7H8O3W -66.575713 -81.59841358 15.02270058
W QAZWUJ-CO C14H12O6W -66.089077 -116.601408 50.51233104
W XODDAW-CO C34H34O5W -65.2532444 -215.5008829 150.2476385

For organic molecules (H, C, O), the results agree with Jaguar and show minimal discrepancies. However, when it comes to tungsten and other transition metals like Re and Ir, the Psi4 energies are highly inconsistent and sometimes deviate by as much as 150 Hartree!

Sample Input File:

molecule FERDIP_CO {
0 1
W           -0.05713407754315       -0.14864659018190       -0.15255709800078
O           -1.69442143536118       -1.63857098157205       -2.33149248076233
C           -1.07739031939630       -1.09829804926988       -1.52755770734550
C           -1.31553100534254        1.64601071076200        0.49991460882643
H           -0.82176552712375        2.58369265398967        0.24710740875750
H           -1.61331095047414        1.55627615525433        1.54407417449809
C           -2.08527618531366        1.01284435506803       -0.47146570867346
H           -2.16105518211603        1.45330887519575       -1.46021534574177
H           -2.94919653955078        0.43030814064546       -0.16873815447254
C           -0.89579336113555       -1.23579074600505        1.28418668071446
O           -1.29246802046766       -1.92572722812787        2.09833498026942
C            0.82025671115132        0.98986754269375       -1.52515066943175
O            1.40294885990693        1.57028068320966       -2.31239595904815
C            1.47043337808395       -1.65807330779779        0.57817639747092
H            1.22325885276241       -2.69061188659033        0.35051419646060
H            1.76002638587163       -1.48458376519588        1.61749901929082
C            2.07719007178722       -0.87156629156821       -0.41536041422695
H            2.28250492043074       -1.31787451836461       -1.38377338029728
H            2.80991017383053       -0.12417681214485       -0.10129088828741
}

set {
  basis def2-SVPD
}

energy('wb97mv')

The only variation between calculations was the molecular coordinates, yet the issue persists across all transition metal-containing molecules. I have tested this with multiple structures, and the results are consistently off when Psi4 is used for these transition metals.

Has anyone else encountered this issue with transition metals? Is there a known fix or a specific setting that I might need to adjust? Any insights would be greatly appreciated.

Thank you for your time and assistance.

Best regards,
Kareem Abdelhafez

What exactly are you calculating? You can calculate the energy of CH4, but how are you partitioning that into “the carbon contribution” and “the hydrogen contribution”? The same comment applies to your other molecules, but methane is a concrete example.

Basically, what I am trying to do is extract the atomic energies of different elements using “standard” molecules like H2 to extract the energy of H by dividing the total energy of H2 by 2. Using H atom’s energy, I extracted the energy of C from CH4 (subtracting 4*H from the total energy of the molecule). Same process for extracting the energy O from H2O.

Now having the “reference atomic energies” of C, H, and O, we can carry the same process to extract the energy of W of different organometallic molecules with the CSD_cod (i.e. mol_id): REDNUK, FERDIP, etc… What’s surprising is that the W atomic energy extracted from those different molecules is highly inconsistent and deviates by large amounts for the isolated W atomic energy.

I carried the same process using Jaguar software using the same set of molecules and the same level of theory, but the results were different, where W atomic energy stayed consistent and not deviating much from the isolated W atomic energy.

What is the rationale for taking the energy of H to be 1/2 the energy of H2 as opposed to computing the energy of H atom?

yea I could have computed H, C, and O from using single atom calculations, but I used stable molecules as a surrogate of “standard molecules” that are stable, from which we can obtain the atomic energies. In any case, the obtained atomic energies are generally close to their isolated energies counterparts and they should be more or less consistent across different molecular environments (with ± 2 Ha), as shown in Jaguar computations for W-containing molecules.

The compelling issue here is with W, where the atomic energy range can go from 66 Ha up to 215 Ha using Psi4 for the same set of W-molecules evaluated with Jaguar!

I would argue that the atomic energies are not close to your derived values. The energy of a hydrogen atom is -0.5 Eh, which means your energy of -0.577114067 Eh is off by more than 200 kJ/mol.

I can’t speak to why Jaguar is giving different energies for the molecules, but there are lots of potential parameters to adjust. If both codes are working correctly, you should be able to obtain identical results between the two of them, starting with H2.

Yes, I agree that the method of evaluating the atomic energies of elements could give different results, whether it is derived from a molecule or an isolated atom. The key is that if we stick to a certain method, we should at least expect consistent values (up to a certain threshold). I would be ok if the W atomic energies range from -64 to -67 Ha. As in the case of Jaguar, the range is from -65.25 to -66.6 Ha. However, the variation I see with Psi4 goes from -66 to -215!

As for the input parameters, I kept the input files consistent keeping everything set to default, only specifying the basis set and the density functional. Again, this is also how I did it with Jaguar.

Thanks for the feedback!

What you’re telling us is that various linear combinations of single-point energies aren’t agreeing. This isn’t enough information for us to go on.

Give us the energies between Psi4 and Jaguar for all your single-points. If the disagreement is very specifically because Psi4 computes high energies for molecules containing tungsten, post your SCF iterations for those molecules.

mol_id formula Jaguar sp energy Psi4 sp energy
H2 H2 -1.154007 -1.154228134
CH4 CH4 -40.443517 -40.44356356
H2O H2O -76.347149 -76.34718062
isolated W W -66.62167 -66.68160669
REDNUK C7H4O5W -711.719749 -760.1431441
FERDIP-CO C7H8O3W -563.719688 -578.7399346
QAZWUJ-CO C14H12O6W -1058.069013 -1108.575994
XODDAW-CO C34H34O5W -1757.444175 -1907.70332

Here is the output file for XODDAW-CO to view scf iterations
XODDAW-CO.txt (49.7 KB)

As you can see in the scf iterations below, the initial guess for the electronic density is very close to the converged energy predicted by Jaguar (-1757.08282934246540), it then jumps to -1905.547 in the first iteration:

@DF-RKS iter SAD: -1757.08282934246540 -1.75708e+03 0.00000e+00
@DF-RKS iter 1: -1905.54724228132773 -1.48464e+02 1.72795e-03 DIIS
@DF-RKS iter 2: -1840.99210992340477 6.45551e+01 9.32183e-03 DIIS
@DF-RKS iter 3: -1906.86456210147776 -6.58725e+01 1.12896e-03 DIIS
@DF-RKS iter 4: -1905.43782761383159 1.42673e+00 1.41857e-03 DIIS
@DF-RKS iter 5: -1907.14387073606690 -1.70604e+00 7.92111e-04 DIIS
@DF-RKS iter 6: -1907.59623588246905 -4.52365e-01 3.02615e-04 DIIS
@DF-RKS iter 7: -1907.65609799190906 -5.98621e-02 1.81589e-04 DIIS
@DF-RKS iter 8: -1907.68117150869693 -2.50735e-02 9.28408e-05 DIIS
@DF-RKS iter 9: -1907.69392754935552 -1.27560e-02 3.66317e-05 DIIS
@DF-RKS iter 10: -1907.69826762971411 -4.34008e-03 3.47514e-05 DIIS
@DF-RKS iter 11: -1907.69946302303788 -1.19539e-03 3.33364e-05 DIIS
@DF-RKS iter 12: -1907.70069366383132 -1.23064e-03 2.91722e-05 DIIS
@DF-RKS iter 13: -1907.70192496526352 -1.23130e-03 1.98540e-05 DIIS
@DF-RKS iter 14: -1907.70281339133498 -8.88426e-04 7.21430e-06 DIIS
@DF-RKS iter 15: -1907.70312104075879 -3.07649e-04 4.69758e-06 DIIS
@DF-RKS iter 16: -1907.70322082777989 -9.97870e-05 3.88449e-06 DIIS
@DF-RKS iter 17: -1907.70326147746437 -4.06497e-05 2.70179e-06 DIIS
@DF-RKS iter 18: -1907.70329555696367 -3.40795e-05 1.45785e-06 DIIS
@DF-RKS iter 19: -1907.70330947615298 -1.39192e-05 1.02536e-06 DIIS
@DF-RKS iter 20: -1907.70331574636748 -6.27021e-06 7.48315e-07 DIIS
@DF-RKS iter 21: -1907.70331777531715 -2.02895e-06 4.16886e-07 DIIS
@DF-RKS iter 22: -1907.70331933362831 -1.55831e-06 2.68930e-07 DIIS
@DF-RKS iter 23: -1907.70331995425590 -6.20628e-07 1.73163e-07 DIIS

The Psi4 code is using density-fitting (DF), but we don’t know if Jaguar is doing the same. Furthermore, we don’t know the definition(s) of Jaguar’s integration grids to compare with those of Psi4. There are multiple ways that the codes could differ, and we don’t have any knowledge here of Jaguar’s internals. (At least, I don’t.)

Psi4 having lower single points than Jaguar is not what I was expecting.

Can you repeat the tungsten single points using Hartree-Fock using both Jaguar and Psi4? My intuition is that you’re having problems with multiple SCF solutions, and that’s easier to analyze at the Hartree-Fock level than DFT. I’d also strongly advise you to use a more recent version of Psi4.

I’m using the latest Psi4 1.9.1 version.

Your output file was “Psi4 1.5 release”.

Oh, my mistake. I have 1.9.1 on my PC but the run was done on an HPC, which is apparently on an older version. I will rerun on my PC and update u.

It seems that the issue is version related! Using 1.9.1

==> Iterations <==

                       Total Energy        Delta E     RMS |[F,P]|

@DF-RKS iter SAD: -1757.27861902064615 -1.75728e+03 0.00000e+00
@DF-RKS iter 1: -1755.94011605362584 1.33850e+00 1.50490e-03 DIIS/ADIIS
@DF-RKS iter 2: -1754.72469095614497 1.21543e+00 2.31154e-03 DIIS/ADIIS
@DF-RKS iter 3: -1757.22800302415476 -2.50331e+00 6.05038e-04 DIIS/ADIIS
@DF-RKS iter 4: -1757.42491368503761 -1.96911e-01 2.44053e-04 DIIS/ADIIS
@DF-RKS iter 5: -1757.45899314327494 -3.40795e-02 9.00216e-05 DIIS
@DF-RKS iter 6: -1757.46521984824199 -6.22670e-03 3.43606e-05 DIIS
@DF-RKS iter 7: -1757.46719574808276 -1.97590e-03 1.63219e-05 DIIS
@DF-RKS iter 8: -1757.46805534752593 -8.59599e-04 7.44219e-06 DIIS
@DF-RKS iter 9: -1757.46830380563756 -2.48458e-04 4.11558e-06 DIIS
@DF-RKS iter 10: -1757.46837859210768 -7.47865e-05 2.21891e-06 DIIS
@DF-RKS iter 11: -1757.46839432329261 -1.57312e-05 7.90736e-07 DIIS
@DF-RKS iter 12: -1757.46839729881049 -2.97552e-06 4.91272e-07 DIIS
@DF-RKS iter 13: -1757.46839787870158 -5.79891e-07 1.88300e-07 DIIS
Energy and wave function converged.

Notice here the energy didn’t jump to -1905.5 after the first scf iteration as it did using psi4/1.5.

I noticed this issue with multiple other molecules with psi4/1.5, like Ir and Re.

Is it generally known that version 1.5 might have issues with specific heavy metals in certain environments, which is solved in version 1.9?

1 Like

We updated our engine to handle ECPs in 1.6, which is probably the cause of the difference.

1 Like