Experiment 8: Stability Analysis and Transition Metals

Is your SCF solution really the ground state?

Excitation energies of Fe and Ru+ with and without SCF stability analysis.

At a glance. Five functionals on the Fe and Ru⁺ states, each run twice (without and with stability analysis), def2-TZVP: 24 jobs. You will learn what an unstable SCF solution is, how to detect and fix it, and why transition metals need extra care.

Warm-up

Loading the concept check…

Background

As you know, the coefficients of the basis functions in the Kohn–Sham orbitals come from an iterative procedure, the self-consistent field (SCF) procedure. Starting from an initial set of coefficients (the initial density matrix, or “guess”), it produces the final, optimized coefficients, which computational chemists call the converged SCF solution.

A converged SCF solution is a stationary point, but we can’t know in advance whether it is a minimum, a maximum or a saddle point. To be sure we have the lowest energy, we must characterize the stationary point and correct it if needed. This stability analysis requires the second derivatives of the energy with respect to orbital rotations (the electronic Hessian), and it can be repeated more than once. When an instability is detected, two constraints can be lifted to lower the energy  [1]: we can change the theoretical framework (from restricted to unrestricted), or we can allow complex instead of real orbitals. The wave function obtained after a stability analysis is stable and has a lower energy than the unstable one.

How can an SCF end on an unstable solution? For example, it can land on an excited state while you are computing the ground state. With a reasonable guess this is unlikely for molecules made mostly of main-group elements, but it is much more likely for systems with nearly degenerate low-lying states, such as transition-metal atoms. Without a stability analysis you would never notice. Good descriptions of the procedure are in  [1],  [2] and  [3].

In this experiment you will use five functionals with the def2-TZVP basis set  [4]: PBE  [5] for both atoms, BP86  [6–8] and PBE0  [9,10] for Ru⁺, and M06-2X  [11] and B97M-V  [12,13] for Fe. These are the combinations whose calculations converge reliably both with and without stability analysis. Other functionals (SPW92, B3LYP, SCAN, ωB97M-V) and Hartree–Fock were left out because their SCF often fails to converge for these atoms, itself a lesson about transition metals. Ru needs an effective core potential (ECP = def2-ECP).

You will compute the excitation energies of two atoms, Fe (reference 34.32 kcal/mol) and Ru⁺ (reference 21.96 kcal/mol). The data come from the 3dAEE8  [14–16] and 4dAEE5  [17] subsets of the Minnesota 2015B database  [18]. Run two sets of calculations:

1. Without stability analysis: compute the excitation energy E(excited) − E(ground).
2. With stability analysis: add INTERNAL_STABILITY TRUE to the $rem section, and check in the output whether an instability was found and the wave function re-optimized.

As usual, compute the errors with respect to the reference values.

The stability landscape

Before running anything, play with this toy model of what stability analysis does.

Loading the stability landscape…

The atoms

Loading the molecules…

Predict first

Loading…

Run the calculations

Each functional/atom pair appears twice: without and with stability analysis (the second set of inputs adds INTERNAL_STABILITY TRUE). In the output of the second run, look for the stability-analysis section, which reports whether an instability was found.

Loading the job list…

Enter your results

Record the total energy, SCF cycles and wall time for both runs, then enter the two excitation energies (excited − ground) for each run in the last tab.

Loading the results table…

Explore your data

Loading…

Questions

Loading the questions…

TipTip

This simple example shows that calculations involving transition metals are not easy. Even for “reactions” between atomic states the errors can be very large, and in extreme cases  [19] it isn’t clear which functionals give the right answer. As general advice, always run a stability analysis on the final solution. It is mandatory for systems with low-lying states, such as transition metals, especially with a small basis set. In these problematic cases, it is always a good idea to seek help from a theoretician.

Submit

Loading…

References

[1]

R. Seeger and J. A. Pople, J. Chem. Phys. 66, 3045 (1977).

[2]

R. Bauernschmitt and R. Ahlrichs, J. Chem. Phys. 104, 9047 (1996).

[3]

H. Fukutome, Int. J. Quantum Chem. 20, 955 (1981).

[4]

F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys. 7, 3297 (2005).

[5]

J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).

[6]

A. D. Becke, Phys. Rev. A 38, 3098 (1988).

[7]

J. P. Perdew, Phys. Rev. B 33, 8822 (1986).

[8]

J. P. Perdew, Phys. Rev. B 34, 7406 (1986).

[9]

C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999).

[10]

M. Ernzerhof and G. E. Scuseria, J. Chem. Phys. 110, 5029 (1999).

[11]

Y. Zhao and D. G. Truhlar, Theor. Chem. Acc. 120, 215 (2008).

[12]

N. Mardirossian and M. Head-Gordon, J. Chem. Phys. 142, 074111 (2015).

[13]

O. A. Vydrov and T. Van Voorhis, J. Chem. Phys. 133, 244103 (2010).

[14]

C. E. Hoyer, G. Li Manni, D. G. Truhlar, and L. Gagliardi, J. Chem. Phys. 141, 204309 (2014).

[15]

S. Luo, B. B. Averkiev, K. R. Yang, X. Xu, and D. G. Truhlar, J. Chem. Theory Comput. 10, 102 (2014).

[16]

H. Yu and D. G. Truhlar, J. Chem. Theory Comput. 10, 2291 (2014).

[17]

S. Luo and D. G. Truhlar, J. Chem. Theory Comput. 8, 4112 (2012).

[18]

H. S. Yu, X. He, S. L. Li, and D. G. Truhlar, Chem. Sci. 7, 5032 (2016).

[19]

L. Cao and U. Ryde, Phys. Chem. Chem. Phys. 21, 2480 (2019).