# Relativistic Time-Dependent Density Functional Theory and Excited States Calculations  for the Zinc Dimer

Authors: Ossama Kullie
Journal: Journal of Atomic, Molecular, and Optical Physics
 (2012)

---

## Abstract

I present a time-dependent density functional study of the 20 low-lying excited states as well the ground states of the zinc dimerZn2, analyze its spectrum obtained from all electrons calculations performed using time-depended density functional with a relativistic 4-component and relativistic spin-free Hamiltonian as implemented in Dirac-Package, and show a comparison of the results obtained from different well-known and newly developed density functional approximations, a comparison with the literature and experimental values as far as available. The results are very encouraging, especially for the lowest excited states of this dimer. However, the results show that long-range corrected functionals such as CAMB3LYP gives the correct asymptotic behavior for the higher states, and for which the best result is obtained. A comparable result is obtained from PBE0 functional. Spin-free Hamiltonian is shown to be very efficient for relativistic systems such as Zn2.


---

## Body

## 1. Introduction

Zinc dimerZn2 is the first member of the group 12 (IIB) (Zn2, Cd2, Hg2, and Cn2) and has a representative character of these dimers. The interest in the dimers of the group IIB (12) is in part due to the possibility of laser applications in analogy with the rare gas dimers. A second point is the importance of the metallic complexes similar to the transition metal complexes [1–4] and some important application like the solar cell and renewable energy [5, 6] as well as electric battery for new cars technology [7, 8]. Zn2, Cd2, and Hg2 are exciter with a shallow, predominantly Van der Waals ground state and low-lying covalent bound excited states. They are also interesting from a theoretical point of view due to the different character of the ground and excited states and consequently the different methodological demands for an accurate theoretical description of the spectrum. The dimer of group 12 has been studied both experimentally and theoretically. Relevant reviews have been provided by Morse [9] and more recently by Koperski [10, 11]. The covalent contributions to the ground state bonding in the group 12 dimers have been investigated in [12], it was concluded that the bond is a mixture of 3/4  Van der Waals and 1/4 covalent interactions. Bucinisky et al. [13] provides spectroscopic constants using the coupled cluster method (CCSD(T)) and different level of the theory 4-component relativistic Hamiltonian, using Dirac-Coulomb Hamiltonian, relativistic spin-free Hamiltonian and nonrelativistic (NR) Hamiltonian. Furthermore, they investigated the relativistic effects and found to be about 5, 8, 19% of the binding energies for  Zn2, Cd2, and Hg2, respectively. Finally the last member of the group Cn2, copernicium, has an academic interest [14–16] due the chemical character of the bonding in comparison to Hg2 (and the lighter dimers of the group), and the influence of the relativistic effects on the atomic orbitals providing a change of the boding character in the dimer to more covalent or Van der Waals type.The paper presents all-electron calculations on the lowest-lying excited states as well as the ground state. The first 8 lowest exited states are discussed with a comparison to experimental and literature values, and several other higher excited states are presented and discussed. Earlier works investigated the lowest 8 excited states using different wave function methods. Ellingsen et al. [17] showed ab initio results for the ground and lowest 8 excited states of Zn2, they performed all electron calculations and present NR as well as relativistic spin-free Douglas-Kroll result, the spin-orbit coupling was accounted perturbatively. The ground state is studied at ACPTF (averaged coupled pair functional, CCSD(T) and CASPT2 (complete active space second-order perturbation theory) level, and the excited states are studied at MR-ACPF (multireference ACPF) and CASPT2 level. Czuchaj et al. [18–20] performed their computations for Zn2 (later for Cd2 and Hg2) using (NR) pseudopotential approach and MRCI (multireference configuration interaction), and the spin-orbit coupling was taken only approximately.In this work, we use a relativistic spin-free Hamiltonian (SFH), without spin-orbit coupling, with a comparison to a relativistic 4-component Dirac-Coulomb Hamiltonian (DCH), spin-orbit coupling included, in the framework of time-dependent density functional theory (TDDFT) and its linear-response approximation (LRA). The calculations are performed using Dirac-Package  (program for atomic and molecular direct iterative relativistic all-electron calculations) [21]. The relativistic effects for Zn2 (and even for Cd2) are small but visible and in some respects not negligible. To my experience, generally around zinc (Z=30) the relativistic effects started to become important for chemical properties. For Hg2, they are large enough (for Cn2 expected to be very large) to make it necessary to incorporate them into any properties that are sensitive to the potential [13]. This is predominantly due to the contraction of 6s orbital, a well-known and important relativistic effects in heavy atoms [22–25]. We will follow this issue in future works on the group 12 (IIB).The paper is organized as follows. Section2 is devoted to the theory and method. We briefly introduce in Section 2.1 the key concepts of the static density functional (DFT) and discuss its extension to the relativistic domain. In Section 2.2, we introduce the key concepts of time-dependent density functional (TDDFT) and the linear response approximation. Section 3 is devoted to the computational details and Section 4 to the result and discussion, and finally we give a conclusion in Section 5. Some useful (well-known) notations used in this paper are collected in Table 1.Table 1
Some of the acronyms used in this work.
HF
Hartree Fock method
NR
Nonrelativistic
DHF
Dirac or relativistic HF
DCH
Dirac-Coulomb Hamiltonian
MP2
Møller-Plesst 2nd-order perturbation theory
CCSD(T)
Coupled cluster singles-doubles (triples)
SFH
Relativistic spin-free Hamiltonian
(TD)DFT
(Time-depended) density functional theory
xc
Exchange-correlation
LR(A)
Linearresponse (approximation)
ALR
Adiabatic LR
srLDAMP2
Short-range LDA, long-range MP2


## 2. Theory and Methods

Time-dependent density functional theory (TDDFT) currently has a growing impact and intensive use in physics and chemistry of atoms, small and large molecules, biomolecules, finite systems, and solidstate. For excited states resulting from a single excitation that present a single jump from the ground state to an excited state, I used in this work the LRA as implemented in  Dirac-Package  [26–28] and well-known approximations of density functionals like LDA (SVWN5 correlation) [29, 30], PBE [31], PB86 [32–34], BPW91 (Becke exchange [32] and Perdew-Wang correlation [35]), long-range corrected PBE0 [36] and its gradient corrected functional GRAC-PBE0 [37, 38], BLYB and B3LYP [32, 39–41], or newly developed range-separated functionals such as CAMB3LYP [42]. Today's available DFT cannot describe the ground state of the group IIB dimers accurately due to a large contribution of dispersion in the bonding [12], despite this when calculating the covalently well-bound excited states the error is reduced considerably, quite possible accompanied with error cancellations.The ground state of the group 12 dimer has a (closed-shell) valence orbitals configuration:(ns2+ns2):σg2,σu2,n=4,5,6 for Zn2-Hg2. This configuration essentially arising from the interaction of atomic (ns) orbitals. It is weakly covalent and preponderantly dispersion interaction, well known especially in the rare gas dimers [43]. The potential curve displays a shallow van der Waals type of minimum. Exciting electrons from σg2 or σu2 to the lowest set of molecular orbitals spanned by the atomic orbitals Atom(ns2) + Atom(nsnp) or Atom(ns2) + Atom(ns(n+1)s), or Atom(ns2) + Atom(ns(n+1)p) gives rise to a manifold of states (see Table 2) among them states which strongly have covalent contributions as we will see in Section 4 Results and Discussion. This makes  TDDFT   using LRA and well-known functional approximations, adequate to describe these states [26].Table 2
Lowest excited states and the corresponding asymptotes.
Equation
(
1
)
(
(
n
s
2
)
1
S
+
(
n
s
n
p
)
3
P
)
:
3
Π
g
,
3
Π
u
,
3
Σ
g
+
,
3
Σ
u
+
Equation
(
2
)
(
(
n
s
2
)
1
S
+
(
n
s
n
p
)
1
P
)
:
1
Π
g
,
1
Π
u
,
1
Σ
g
+
,
1
Σ
u
+
Equation
(
3
)
(
(
n
s
2
)
1
S
+
(
n
s
(
n
+
1
)
s
)
3
S
)
:
3
Σ
g
+
,
3
Σ
u
+
Equation
(
4
)
(
(
n
s
2
)
1
S
+
(
n
s
(
n
+
1
)
s
)
1
S
)
:
1
Σ
g
+
,
1
Σ
u
+
Equation
(
5
)
(
(
n
s
2
)
1
S
+
(
n
s
(
n
+
1
)
p
)
3
P
)
:
3
Π
g
,
3
Π
u
,
3
Σ
g
+
,
3
Σ
u
+
Equation
(
6
)
(
(
n
s
2
)
1
S
+
(
n
s
(
n
+
1
)
p
)
1
P
)
:
1
Π
g
,
1
Π
u
,
1
Σ
g
+
,
1
Σ
u
+We will discuss the lowest 20 excited states dissociating to the atomic asymptotes (NR notation) given in Table2, resulting from exciting one electron from the ground state (4s21S+4s21S)1Σg+. The concern will be in the first place on the 8 lowest excited states corresponding to the asymptote Atom(ns2) + Atom(nsnp). States corresponding to the higher asymptotes Atom(ns2) + Atom(ns(n+1)s) and Atom(ns2) + Atom(ns(n+1)p) are computed and some of them are well-bound states, we will discuss their quality in view of the limit of the validity of the known DFT approximations yielding inaccurate potential curves and causing a disturbance near the avoiding crossing with states of the same symmetry (see Section 4). To my best knowledge, there is no experimental or theoretical values from DFT or wave function methods available for the higher states to compare with, this makes it difficult to judge the result of the present work. It is expected that the result of the lowest states will show an excellent agreement with the experimental data [10, 11] (and the references therein), whereas for the higher states a satisfactory result is expected showing the important features of these states. The comparison between spin-free and 4-component results shows clearly the capability of SFH to deal with the computation of the properties of the Zn2 dimer or similar systems. We also emphasize its importance for heavier relativistic systems [13], although spin-orbit effect is expected to be larger for Cd2, Hg2, and Cn2. Pyper et al. [22] pointed out that the relativistic ground-state potential well depth of Hg2  is 45% of the NR one and clearly it is stronger for Cn2.


### 2.1. Density Functional Theory

Density functional theory [44–46] has become recently a very large popularity as a good compromise between accuracy and computational expediency. The Hohenberg-Kohn theorem [44] proves the existence of an unique (up to an additive constant) external potential vext(r) for a given nondegenerate density n(r) of interacting Fermions. The key point behind this scheme is the very useful simplification, namely, the transformation of the many-body quantum problem to a set of equations of one-particle Schrödinger (or Dirac) type of a noninteracting reference system with the density as a central ingredient quantity to carry all the relevant information of the system under consideration, instead of the many-body quantum wave function in which all the information of the system is stored:(1)H^ϕi(r)=E[n(r)]ϕi(r),(2)H^=T^+Veff[n(r)]=∑i‍t^(ri)+veff(ri)[n(r)],(3)veff(ri)=vext(ri)+vH(ri)+vxc(ri)+vnn,(4)n(r)=∑i=1N‍‖ϕi(r)‖2,
where n(r) is the total density of the system and the sum is over N, that is, all occupied orbitals ϕi(r). t^(ri) is the one-particle kinetic energy operator, veff(ri) is the one-particle effective potential (also called Kohn-Sham potential veff(ri)≡vKS(ri)), with vext(ri) is the Coulombic interaction of the electron i with all the nuclei, called the external potential. vH(ri) is the Hartree and vxc(ri) exchange-correlation potential. And vnn is the classical Coulombic repulsion of the nuclei in the system. vH(ri) is given by the usual expression, but the crucial part vxc(ri) in this scheme is the explicitly unknown vxc(ri):(5)vH(ri)=∫‍d3rn(r)|ri-r|,vxc(ri)=∂Exc[n(r)]∂n(ri),
for which an appropriate good approximation must be found. Experiences in DFT (and  TDDFT) over the past decades shows that the density of atoms, molecules, finite systems, and solids have very complicated structures [47]. To find a good mathematical functionality form between the density (and its gradients) and an exchange-correlation potential with widely physical applications success is one of the most challenging problems in quantum physics and chemistry. Moreover, most of the problems arise when evaluating the results of the calculating systems can be tracked back to the limits of the validity of the today’s known and employed approximations specially the long-range behavior leaving quite a room for improvements. One should note that that in many applications the usual approximations are quite reliable and give good results and acceptable accuracies. The present work is not an exception as we will see when analyzing the results of the ground state and excited states of the Zn2 dimer.


#### 2.1.1. Density Functional Theory in the Relativistic Domain

In the relativistic Dirac theory in absence of electromagnetic field, the DCH has the same generic form as the NR Hamiltonian (for molecules) [26, 48]:(6)H^DC=∑iN‍hD(i)+12∑i≠jN‍g^Coul(i,j)+∑K≠K′M‍VK,K′nn,hD(i)=(c2β^+cα^⋅p^(i)-c2⋅I4)+I4⋅∑K=1M‍VKext(i),α^j=(0σjσj0),j=x,y,z;β^=(I200-I2),
where hD(i) is the one-particle DCH, and c is the speed of light in atomic units (atomic units are used throughout this work unless otherwise noted). Vnn is the classical nucleus-nucleus repulsion and VKext(i)=-ZK/riKis the external Coulombic interaction of the electron i with the nucleus K, and the sum is over all nuclei M. I2 and I4 are the 2×2- and 4×4-unity matrix and the term c2·I4 is a shift to align the relativistic and NR energy scales. β^ and α^=(αx,αy,αz) are the Dirac matrices, with the well-known Pauli matrices σ′s. The generic term
(7)g^Coul(i,j)=I4×I4rij
is the Coulombic instantaneous two-electron i,j interaction operator, it contains in the relativistic theory the spin-own orbit interaction. The DCH approximation reduces the density functional theory in the relativistic domain to the usual density functional theory with the density as the central ingredient, and there is no need to introduce the current density [48]. A density functional theory in the relativistic domain can be constructed on the the basis of (1)–(4) with the density is constructed from the relativistic 4-component wave function. The total energy of the system is given by
(8)E[n]=∑iN‍εi-EJ[n]+Exc[n]-∫‍d3rvxc(r,n)n(r)+Enn,
where εi are the electronic eigenvalues of the system and are calculated iteratively in a self-consistent manner (SCF iterations) in an effective many-body potential veff given in (3). Enn is the nuclear-nuclear repulsion energy, EJ[n] is the Hartree energy equation (9), and Exc[n] is the exchange-correlation energy, it can be further divided into exchange and correlation parts Exc[n]=Ex[n]+Ec[n]. At the (single, determinant) Hartree-Fock (HF) level, which in the relativistic calculations is usually called Dirac-Hartree-Fock (DHF), the two-particle interaction, the Hartree and exact exchange are given by (9) and (10) as follows:
(9)EJ[n]=12∫‍∫‍d3r1d3r2n(r1)n(r2)|r1-r2|,(10)Ex=-14∑i,jN‍∫‍∫‍d3r1d3r2ϕi†(r1)ϕj†(r2)ϕj(r1)ϕi(r2)|r1-r2|,
where r1 and r2 denote the coordinates of the electron one and two, respectively. EJ[n] is a classical interaction between two one-particle densities n(r1) and n(r2), whereas Ex is a quantum mechanical nonlocal part of many-particle interaction. The ϕ(r)s are the electronic one-particle HF-orbitals and the sum is over all the occupied orbitals N. A well-known approximation for the Hartree-Fock exchange energy is the (α-)Slater approximation [29] with remarkable performance for covalent bonding in covalently bound molecules with heavy atoms [49, 50]
(11)Exα[n]=-32αCx∫‍d3rn4/3(r),
where Cx=(3/4)(3/π)1/3 is a constant, in the Slater approximation the parameter α=0.7 is chosen. The missing of the correlation made the Slater approximation unpopular for chemical calculations. In the DFT, the exact Exc[n] is unknown as a functional of the density (and its gradients). Many approximations exist with different performance and accuracy depending on their application area. In LDA one assumes a slowly varying local density dependence; hence, the Dirac-formula [51] of the exchange energy for an uniform electronic gas equation (11) with α=2/3 is applied and the Vosko-Wilk-Nusair correlation formula [29, 30] for the correlation energy (we use SVWN5). LDA depends only on the density, whereas in the generalized gradient approximation (GGA) the density and its gradient are involved, meta GGAs [52] include higher gradients, this systematic improvements is known in the DFT community under the term “Jacob's ladder.” In hybrid functional, for example, BLYP and B3LYP [32, 39–41], one add a (fixed) suitable fraction of exact (Hartree-Fock) exchange (10) to the approximate x-energy part, which often improves the performance of the DFT approximation, whereas in the range-separated density functional [53] a parametric fraction of exchange (and possibly correlation) from wave function methods are added to the DFT exchange energy, with the parameter dictate the amount of exchange to be added, like CAMB3LYP [42], or of exchange-correlation like srLDAMP2 (see [43, 54–56] and the references therein), this improves the results considerably, unfortunately it is found that the optimum parameter value depends on the specific property of the system.


#### 2.1.2. The Relativistic 4-Component and SFH

The Dirac equation with the Dirac-Coulomb Hamiltonian (DCH) describes the important relativistic effects for chemical calculation, which become large for systems with largeZ. It is a firs-order differential equation(s), hence nonvariational “variational collapse” in contrast to the second-order differential Schrödinger equation in the NR case. The solutions to the Dirac equation describe both positrons (the “negative energy” states) and electrons (the “positive energy” states) as well as both spin orientations and a four-component wave function is involved called Dirac spinors:
(12)|ψ〉=(ΨLΨS),ΨL=(ϕ1ϕ2),ΨS=(ϕ3ϕ4),
where ΨL is called the large and ΨS the small component. This notation originally comes from the well-known kinetic balance approximation and is justified by the relation ~1/c between them, from which it follows the NR limit limc→∞ΨS=0 and one identify ΨL with the 2-component vector (spin up; down) of the Schrödinger equation. The full relativistic 4-component DCH is computationally demanding; therefore, it is desirable to reduce the computational effort in relativistic calculations by reducing the dimension of the involved quantities, normally by reducing or transforming the Hamiltonian to a new from, so that the calculations involving operators acting only on the large components and requiring a moderate computational effort by keeping the main physical features of the results. The relativistic SFH implemented in  Dirac-Package uses the Dyall's formulation [57] to obtain results without spin-orbit coupling for the four-component Hamiltonian in the default restricted kinetic balance scheme. In Section 4, we show that the results obtained for the excited states of Zn2 based on (relativistic) SFH are accurate similar and well comparable to those obtained from the 4-component DCH. For the deriving of this Hamiltonian, we kindly refer the reader to [57], see also [58] with advanced description in framework of second quantization formalism. The relativistic SFH permits factorization of the spin as in NR calculations so that standard NR post-SCF methods can be used for inclusion of electron correlation. The extension and implementation of relativistic SFH for many-body system or molecular calculation is straightforward see [21].


### 2.2. TDDFT  and Linear Response

In this section, we briefly introduce  TDDFT   formulation with a special emphasis on the linear density-response function and its connection to the electronic excitation spectrum, a more extensive derivations and wide discussions can be found in refs [47, 59–78] and the references therein.  TDDFT   was pioneered by a work of Zangwill and Soven [78], but the fundamental step was done later by Runge and Gross [60, 61], the Runge-Gross theorem is a rigorous foundation for the formally extension of the Hohenberg-Kohn theorem [44] to the time-dependent phenomena. It results in a time-dependent Kohn-Sham equation:
(13)[T^+vext,σ[n](rt)+vH[n](rt)+vxc,σ[n](rt)]ψjσ(rt)=i∂∂tψjσ(rt),
where T^ is the kinetic energy, vext,σ(rt),vH(rt), are vxc,σ(rt) are the time-dependent external, Hartree, and exchange-correlation potential respectively, and we adopt the notation (rt)≡(r,t). ψjσ(rt) is the wave function of a particle j with a spin σ. The external potential is unique determined via the total density:
(14)n(rt)=∑σ‍nσ(rt)=∑σ‍∑jNσ‍‖ψjσ(rt)‖
of the interacting system, where the sum is taken over all occupied spin-orbitals Nσ of a spin possibility σ.


#### 2.2.1. Linear Response

In the special case of the response of the ground-state density to a weak external field, that is, the case in the most optical applications, the slightly perturbed system, which can be written in a series expansionvext=vext0+vext1+⋯≈vext0+δvext, see [72], starts its evolution slowly from its ground-state density n0 corresponding to the ground-state external potential vext0. The xc can be expressed in terms of the states of (unperturbed) system, and thus as a functional of the ground-state density. The interacting real system and the Kohn-Sham fictitious system are connected via the same infinitesimal density change δn(rt). The infinitesimal change in the Hartree-xc-potential δvHxc=δvH+δvxc due to the infinitesimal change in the density can be expressed in its functional derivative:
(15)δvHxc(rt)=∫‍d3r′dt′fHxc(rr′,t-t′)δn(r′t′),
where fHxc is called the Hartree-xc-kernel and is given in LR regime by
(16)fHxc[n0](rr′;t-t′)=δ(t-t′)|r-r′|+δvxc[n](rt)δn(r′t′)|n=n0(r),
where δ(t-t′) is the Dirac-delta function. The first term in (16) is the Hartree contribution, it is instantaneous, or local in time. The second term in (16), fxc[n0], called the xc-kernel, is much simpler than vxc[n](rt) since it is a functional of the ground-state density n0, it is nonlocal in space and time [70].In the adiabatic approximation which is the most common in  TDDFT, one ignores all time-dependencies in the past and takes only the instantaneous densityn(t) being local in time. The adiabatic approach is a drastic simplification and a priori only justified for systems with a weak time-dependence, which are always locally close to equilibrium [72]. In practice, one takes a known ground-state functional approximation and insert n0(t) into it; thus, any ground-state approximation (LDA, GGA, …) provides an adiabatic approximation for the  TDDFT   xc-functional. The most common one is the ALDA.


## 2.1. Density Functional Theory

Density functional theory [44–46] has become recently a very large popularity as a good compromise between accuracy and computational expediency. The Hohenberg-Kohn theorem [44] proves the existence of an unique (up to an additive constant) external potential vext(r) for a given nondegenerate density n(r) of interacting Fermions. The key point behind this scheme is the very useful simplification, namely, the transformation of the many-body quantum problem to a set of equations of one-particle Schrödinger (or Dirac) type of a noninteracting reference system with the density as a central ingredient quantity to carry all the relevant information of the system under consideration, instead of the many-body quantum wave function in which all the information of the system is stored:(1)H^ϕi(r)=E[n(r)]ϕi(r),(2)H^=T^+Veff[n(r)]=∑i‍t^(ri)+veff(ri)[n(r)],(3)veff(ri)=vext(ri)+vH(ri)+vxc(ri)+vnn,(4)n(r)=∑i=1N‍‖ϕi(r)‖2,
where n(r) is the total density of the system and the sum is over N, that is, all occupied orbitals ϕi(r). t^(ri) is the one-particle kinetic energy operator, veff(ri) is the one-particle effective potential (also called Kohn-Sham potential veff(ri)≡vKS(ri)), with vext(ri) is the Coulombic interaction of the electron i with all the nuclei, called the external potential. vH(ri) is the Hartree and vxc(ri) exchange-correlation potential. And vnn is the classical Coulombic repulsion of the nuclei in the system. vH(ri) is given by the usual expression, but the crucial part vxc(ri) in this scheme is the explicitly unknown vxc(ri):(5)vH(ri)=∫‍d3rn(r)|ri-r|,vxc(ri)=∂Exc[n(r)]∂n(ri),
for which an appropriate good approximation must be found. Experiences in DFT (and  TDDFT) over the past decades shows that the density of atoms, molecules, finite systems, and solids have very complicated structures [47]. To find a good mathematical functionality form between the density (and its gradients) and an exchange-correlation potential with widely physical applications success is one of the most challenging problems in quantum physics and chemistry. Moreover, most of the problems arise when evaluating the results of the calculating systems can be tracked back to the limits of the validity of the today’s known and employed approximations specially the long-range behavior leaving quite a room for improvements. One should note that that in many applications the usual approximations are quite reliable and give good results and acceptable accuracies. The present work is not an exception as we will see when analyzing the results of the ground state and excited states of the Zn2 dimer.


### 2.1.1. Density Functional Theory in the Relativistic Domain

In the relativistic Dirac theory in absence of electromagnetic field, the DCH has the same generic form as the NR Hamiltonian (for molecules) [26, 48]:(6)H^DC=∑iN‍hD(i)+12∑i≠jN‍g^Coul(i,j)+∑K≠K′M‍VK,K′nn,hD(i)=(c2β^+cα^⋅p^(i)-c2⋅I4)+I4⋅∑K=1M‍VKext(i),α^j=(0σjσj0),j=x,y,z;β^=(I200-I2),
where hD(i) is the one-particle DCH, and c is the speed of light in atomic units (atomic units are used throughout this work unless otherwise noted). Vnn is the classical nucleus-nucleus repulsion and VKext(i)=-ZK/riKis the external Coulombic interaction of the electron i with the nucleus K, and the sum is over all nuclei M. I2 and I4 are the 2×2- and 4×4-unity matrix and the term c2·I4 is a shift to align the relativistic and NR energy scales. β^ and α^=(αx,αy,αz) are the Dirac matrices, with the well-known Pauli matrices σ′s. The generic term
(7)g^Coul(i,j)=I4×I4rij
is the Coulombic instantaneous two-electron i,j interaction operator, it contains in the relativistic theory the spin-own orbit interaction. The DCH approximation reduces the density functional theory in the relativistic domain to the usual density functional theory with the density as the central ingredient, and there is no need to introduce the current density [48]. A density functional theory in the relativistic domain can be constructed on the the basis of (1)–(4) with the density is constructed from the relativistic 4-component wave function. The total energy of the system is given by
(8)E[n]=∑iN‍εi-EJ[n]+Exc[n]-∫‍d3rvxc(r,n)n(r)+Enn,
where εi are the electronic eigenvalues of the system and are calculated iteratively in a self-consistent manner (SCF iterations) in an effective many-body potential veff given in (3). Enn is the nuclear-nuclear repulsion energy, EJ[n] is the Hartree energy equation (9), and Exc[n] is the exchange-correlation energy, it can be further divided into exchange and correlation parts Exc[n]=Ex[n]+Ec[n]. At the (single, determinant) Hartree-Fock (HF) level, which in the relativistic calculations is usually called Dirac-Hartree-Fock (DHF), the two-particle interaction, the Hartree and exact exchange are given by (9) and (10) as follows:
(9)EJ[n]=12∫‍∫‍d3r1d3r2n(r1)n(r2)|r1-r2|,(10)Ex=-14∑i,jN‍∫‍∫‍d3r1d3r2ϕi†(r1)ϕj†(r2)ϕj(r1)ϕi(r2)|r1-r2|,
where r1 and r2 denote the coordinates of the electron one and two, respectively. EJ[n] is a classical interaction between two one-particle densities n(r1) and n(r2), whereas Ex is a quantum mechanical nonlocal part of many-particle interaction. The ϕ(r)s are the electronic one-particle HF-orbitals and the sum is over all the occupied orbitals N. A well-known approximation for the Hartree-Fock exchange energy is the (α-)Slater approximation [29] with remarkable performance for covalent bonding in covalently bound molecules with heavy atoms [49, 50]
(11)Exα[n]=-32αCx∫‍d3rn4/3(r),
where Cx=(3/4)(3/π)1/3 is a constant, in the Slater approximation the parameter α=0.7 is chosen. The missing of the correlation made the Slater approximation unpopular for chemical calculations. In the DFT, the exact Exc[n] is unknown as a functional of the density (and its gradients). Many approximations exist with different performance and accuracy depending on their application area. In LDA one assumes a slowly varying local density dependence; hence, the Dirac-formula [51] of the exchange energy for an uniform electronic gas equation (11) with α=2/3 is applied and the Vosko-Wilk-Nusair correlation formula [29, 30] for the correlation energy (we use SVWN5). LDA depends only on the density, whereas in the generalized gradient approximation (GGA) the density and its gradient are involved, meta GGAs [52] include higher gradients, this systematic improvements is known in the DFT community under the term “Jacob's ladder.” In hybrid functional, for example, BLYP and B3LYP [32, 39–41], one add a (fixed) suitable fraction of exact (Hartree-Fock) exchange (10) to the approximate x-energy part, which often improves the performance of the DFT approximation, whereas in the range-separated density functional [53] a parametric fraction of exchange (and possibly correlation) from wave function methods are added to the DFT exchange energy, with the parameter dictate the amount of exchange to be added, like CAMB3LYP [42], or of exchange-correlation like srLDAMP2 (see [43, 54–56] and the references therein), this improves the results considerably, unfortunately it is found that the optimum parameter value depends on the specific property of the system.


### 2.1.2. The Relativistic 4-Component and SFH

The Dirac equation with the Dirac-Coulomb Hamiltonian (DCH) describes the important relativistic effects for chemical calculation, which become large for systems with largeZ. It is a firs-order differential equation(s), hence nonvariational “variational collapse” in contrast to the second-order differential Schrödinger equation in the NR case. The solutions to the Dirac equation describe both positrons (the “negative energy” states) and electrons (the “positive energy” states) as well as both spin orientations and a four-component wave function is involved called Dirac spinors:
(12)|ψ〉=(ΨLΨS),ΨL=(ϕ1ϕ2),ΨS=(ϕ3ϕ4),
where ΨL is called the large and ΨS the small component. This notation originally comes from the well-known kinetic balance approximation and is justified by the relation ~1/c between them, from which it follows the NR limit limc→∞ΨS=0 and one identify ΨL with the 2-component vector (spin up; down) of the Schrödinger equation. The full relativistic 4-component DCH is computationally demanding; therefore, it is desirable to reduce the computational effort in relativistic calculations by reducing the dimension of the involved quantities, normally by reducing or transforming the Hamiltonian to a new from, so that the calculations involving operators acting only on the large components and requiring a moderate computational effort by keeping the main physical features of the results. The relativistic SFH implemented in  Dirac-Package uses the Dyall's formulation [57] to obtain results without spin-orbit coupling for the four-component Hamiltonian in the default restricted kinetic balance scheme. In Section 4, we show that the results obtained for the excited states of Zn2 based on (relativistic) SFH are accurate similar and well comparable to those obtained from the 4-component DCH. For the deriving of this Hamiltonian, we kindly refer the reader to [57], see also [58] with advanced description in framework of second quantization formalism. The relativistic SFH permits factorization of the spin as in NR calculations so that standard NR post-SCF methods can be used for inclusion of electron correlation. The extension and implementation of relativistic SFH for many-body system or molecular calculation is straightforward see [21].


## 2.1.1. Density Functional Theory in the Relativistic Domain

In the relativistic Dirac theory in absence of electromagnetic field, the DCH has the same generic form as the NR Hamiltonian (for molecules) [26, 48]:(6)H^DC=∑iN‍hD(i)+12∑i≠jN‍g^Coul(i,j)+∑K≠K′M‍VK,K′nn,hD(i)=(c2β^+cα^⋅p^(i)-c2⋅I4)+I4⋅∑K=1M‍VKext(i),α^j=(0σjσj0),j=x,y,z;β^=(I200-I2),
where hD(i) is the one-particle DCH, and c is the speed of light in atomic units (atomic units are used throughout this work unless otherwise noted). Vnn is the classical nucleus-nucleus repulsion and VKext(i)=-ZK/riKis the external Coulombic interaction of the electron i with the nucleus K, and the sum is over all nuclei M. I2 and I4 are the 2×2- and 4×4-unity matrix and the term c2·I4 is a shift to align the relativistic and NR energy scales. β^ and α^=(αx,αy,αz) are the Dirac matrices, with the well-known Pauli matrices σ′s. The generic term
(7)g^Coul(i,j)=I4×I4rij
is the Coulombic instantaneous two-electron i,j interaction operator, it contains in the relativistic theory the spin-own orbit interaction. The DCH approximation reduces the density functional theory in the relativistic domain to the usual density functional theory with the density as the central ingredient, and there is no need to introduce the current density [48]. A density functional theory in the relativistic domain can be constructed on the the basis of (1)–(4) with the density is constructed from the relativistic 4-component wave function. The total energy of the system is given by
(8)E[n]=∑iN‍εi-EJ[n]+Exc[n]-∫‍d3rvxc(r,n)n(r)+Enn,
where εi are the electronic eigenvalues of the system and are calculated iteratively in a self-consistent manner (SCF iterations) in an effective many-body potential veff given in (3). Enn is the nuclear-nuclear repulsion energy, EJ[n] is the Hartree energy equation (9), and Exc[n] is the exchange-correlation energy, it can be further divided into exchange and correlation parts Exc[n]=Ex[n]+Ec[n]. At the (single, determinant) Hartree-Fock (HF) level, which in the relativistic calculations is usually called Dirac-Hartree-Fock (DHF), the two-particle interaction, the Hartree and exact exchange are given by (9) and (10) as follows:
(9)EJ[n]=12∫‍∫‍d3r1d3r2n(r1)n(r2)|r1-r2|,(10)Ex=-14∑i,jN‍∫‍∫‍d3r1d3r2ϕi†(r1)ϕj†(r2)ϕj(r1)ϕi(r2)|r1-r2|,
where r1 and r2 denote the coordinates of the electron one and two, respectively. EJ[n] is a classical interaction between two one-particle densities n(r1) and n(r2), whereas Ex is a quantum mechanical nonlocal part of many-particle interaction. The ϕ(r)s are the electronic one-particle HF-orbitals and the sum is over all the occupied orbitals N. A well-known approximation for the Hartree-Fock exchange energy is the (α-)Slater approximation [29] with remarkable performance for covalent bonding in covalently bound molecules with heavy atoms [49, 50]
(11)Exα[n]=-32αCx∫‍d3rn4/3(r),
where Cx=(3/4)(3/π)1/3 is a constant, in the Slater approximation the parameter α=0.7 is chosen. The missing of the correlation made the Slater approximation unpopular for chemical calculations. In the DFT, the exact Exc[n] is unknown as a functional of the density (and its gradients). Many approximations exist with different performance and accuracy depending on their application area. In LDA one assumes a slowly varying local density dependence; hence, the Dirac-formula [51] of the exchange energy for an uniform electronic gas equation (11) with α=2/3 is applied and the Vosko-Wilk-Nusair correlation formula [29, 30] for the correlation energy (we use SVWN5). LDA depends only on the density, whereas in the generalized gradient approximation (GGA) the density and its gradient are involved, meta GGAs [52] include higher gradients, this systematic improvements is known in the DFT community under the term “Jacob's ladder.” In hybrid functional, for example, BLYP and B3LYP [32, 39–41], one add a (fixed) suitable fraction of exact (Hartree-Fock) exchange (10) to the approximate x-energy part, which often improves the performance of the DFT approximation, whereas in the range-separated density functional [53] a parametric fraction of exchange (and possibly correlation) from wave function methods are added to the DFT exchange energy, with the parameter dictate the amount of exchange to be added, like CAMB3LYP [42], or of exchange-correlation like srLDAMP2 (see [43, 54–56] and the references therein), this improves the results considerably, unfortunately it is found that the optimum parameter value depends on the specific property of the system.


## 2.1.2. The Relativistic 4-Component and SFH

The Dirac equation with the Dirac-Coulomb Hamiltonian (DCH) describes the important relativistic effects for chemical calculation, which become large for systems with largeZ. It is a firs-order differential equation(s), hence nonvariational “variational collapse” in contrast to the second-order differential Schrödinger equation in the NR case. The solutions to the Dirac equation describe both positrons (the “negative energy” states) and electrons (the “positive energy” states) as well as both spin orientations and a four-component wave function is involved called Dirac spinors:
(12)|ψ〉=(ΨLΨS),ΨL=(ϕ1ϕ2),ΨS=(ϕ3ϕ4),
where ΨL is called the large and ΨS the small component. This notation originally comes from the well-known kinetic balance approximation and is justified by the relation ~1/c between them, from which it follows the NR limit limc→∞ΨS=0 and one identify ΨL with the 2-component vector (spin up; down) of the Schrödinger equation. The full relativistic 4-component DCH is computationally demanding; therefore, it is desirable to reduce the computational effort in relativistic calculations by reducing the dimension of the involved quantities, normally by reducing or transforming the Hamiltonian to a new from, so that the calculations involving operators acting only on the large components and requiring a moderate computational effort by keeping the main physical features of the results. The relativistic SFH implemented in  Dirac-Package uses the Dyall's formulation [57] to obtain results without spin-orbit coupling for the four-component Hamiltonian in the default restricted kinetic balance scheme. In Section 4, we show that the results obtained for the excited states of Zn2 based on (relativistic) SFH are accurate similar and well comparable to those obtained from the 4-component DCH. For the deriving of this Hamiltonian, we kindly refer the reader to [57], see also [58] with advanced description in framework of second quantization formalism. The relativistic SFH permits factorization of the spin as in NR calculations so that standard NR post-SCF methods can be used for inclusion of electron correlation. The extension and implementation of relativistic SFH for many-body system or molecular calculation is straightforward see [21].


## 2.2. TDDFT  and Linear Response

In this section, we briefly introduce  TDDFT   formulation with a special emphasis on the linear density-response function and its connection to the electronic excitation spectrum, a more extensive derivations and wide discussions can be found in refs [47, 59–78] and the references therein.  TDDFT   was pioneered by a work of Zangwill and Soven [78], but the fundamental step was done later by Runge and Gross [60, 61], the Runge-Gross theorem is a rigorous foundation for the formally extension of the Hohenberg-Kohn theorem [44] to the time-dependent phenomena. It results in a time-dependent Kohn-Sham equation:
(13)[T^+vext,σ[n](rt)+vH[n](rt)+vxc,σ[n](rt)]ψjσ(rt)=i∂∂tψjσ(rt),
where T^ is the kinetic energy, vext,σ(rt),vH(rt), are vxc,σ(rt) are the time-dependent external, Hartree, and exchange-correlation potential respectively, and we adopt the notation (rt)≡(r,t). ψjσ(rt) is the wave function of a particle j with a spin σ. The external potential is unique determined via the total density:
(14)n(rt)=∑σ‍nσ(rt)=∑σ‍∑jNσ‍‖ψjσ(rt)‖
of the interacting system, where the sum is taken over all occupied spin-orbitals Nσ of a spin possibility σ.


### 2.2.1. Linear Response

In the special case of the response of the ground-state density to a weak external field, that is, the case in the most optical applications, the slightly perturbed system, which can be written in a series expansionvext=vext0+vext1+⋯≈vext0+δvext, see [72], starts its evolution slowly from its ground-state density n0 corresponding to the ground-state external potential vext0. The xc can be expressed in terms of the states of (unperturbed) system, and thus as a functional of the ground-state density. The interacting real system and the Kohn-Sham fictitious system are connected via the same infinitesimal density change δn(rt). The infinitesimal change in the Hartree-xc-potential δvHxc=δvH+δvxc due to the infinitesimal change in the density can be expressed in its functional derivative:
(15)δvHxc(rt)=∫‍d3r′dt′fHxc(rr′,t-t′)δn(r′t′),
where fHxc is called the Hartree-xc-kernel and is given in LR regime by
(16)fHxc[n0](rr′;t-t′)=δ(t-t′)|r-r′|+δvxc[n](rt)δn(r′t′)|n=n0(r),
where δ(t-t′) is the Dirac-delta function. The first term in (16) is the Hartree contribution, it is instantaneous, or local in time. The second term in (16), fxc[n0], called the xc-kernel, is much simpler than vxc[n](rt) since it is a functional of the ground-state density n0, it is nonlocal in space and time [70].In the adiabatic approximation which is the most common in  TDDFT, one ignores all time-dependencies in the past and takes only the instantaneous densityn(t) being local in time. The adiabatic approach is a drastic simplification and a priori only justified for systems with a weak time-dependence, which are always locally close to equilibrium [72]. In practice, one takes a known ground-state functional approximation and insert n0(t) into it; thus, any ground-state approximation (LDA, GGA, …) provides an adiabatic approximation for the  TDDFT   xc-functional. The most common one is the ALDA.


## 2.2.1. Linear Response

In the special case of the response of the ground-state density to a weak external field, that is, the case in the most optical applications, the slightly perturbed system, which can be written in a series expansionvext=vext0+vext1+⋯≈vext0+δvext, see [72], starts its evolution slowly from its ground-state density n0 corresponding to the ground-state external potential vext0. The xc can be expressed in terms of the states of (unperturbed) system, and thus as a functional of the ground-state density. The interacting real system and the Kohn-Sham fictitious system are connected via the same infinitesimal density change δn(rt). The infinitesimal change in the Hartree-xc-potential δvHxc=δvH+δvxc due to the infinitesimal change in the density can be expressed in its functional derivative:
(15)δvHxc(rt)=∫‍d3r′dt′fHxc(rr′,t-t′)δn(r′t′),
where fHxc is called the Hartree-xc-kernel and is given in LR regime by
(16)fHxc[n0](rr′;t-t′)=δ(t-t′)|r-r′|+δvxc[n](rt)δn(r′t′)|n=n0(r),
where δ(t-t′) is the Dirac-delta function. The first term in (16) is the Hartree contribution, it is instantaneous, or local in time. The second term in (16), fxc[n0], called the xc-kernel, is much simpler than vxc[n](rt) since it is a functional of the ground-state density n0, it is nonlocal in space and time [70].In the adiabatic approximation which is the most common in  TDDFT, one ignores all time-dependencies in the past and takes only the instantaneous densityn(t) being local in time. The adiabatic approach is a drastic simplification and a priori only justified for systems with a weak time-dependence, which are always locally close to equilibrium [72]. In practice, one takes a known ground-state functional approximation and insert n0(t) into it; thus, any ground-state approximation (LDA, GGA, …) provides an adiabatic approximation for the  TDDFT   xc-functional. The most common one is the ALDA.


## 3. Computational Details

The reported results in this paper have been performed using a development version of the  Dirac10-Package [21] based on the 4-component relativistic DCH and SFH. We would like to stress, though, that the present implementation allows the use of all Hamiltonians implemented in the  Dirac-Package  such as the eXact 2-component relativistic Hamiltonian (X2C) [79] and the 4-component NR Lévy-Leblond Hamiltonian [80]. The nuclear charge distribution was described by a Gaussian model using the recommended values of [81].The values of the spectroscopic constantsRe, ωe, and De were extracted from a Morse potential fit based on at least ten equidistant points of step length 0.05 a.u. around the equilibrium distance a second fit using polynomial fit procedure available in  Dirac-Package  is used too, the comparison between the two fits show that 5-order polynomial fit is rather equivalent to a Morse potential fit, provided that Morse potential fit is performed for small region around the minimum which is done throughout this work, the agreement between the two fits gives us an additional criterion for the safety and correctness of the calculated spectroscopic constants reported in the present result.We employed the aug-cc-pVTZ (likewise aug-cc-pVQZ) Gaussian basis sets of Dunning and coworkers [82–84]. This basis set is widely used in the literature, thus simplifying the comparison between different works. The small components basis set for the 4-component relativistic calculations has been generated using restricted kinetic balance imposed in the canonical orthogonalization step [80]. All basis sets are used in uncontracted form. Test calculations with aug-cc-pVQZ basis sets indicate that the reported structures can be considered converged with respect to the chosen basis sets, see Section 4. The potential curves are generated with a bout 175 point densely chosen equidistant with of step length of 0.05 a.u. in the significant part of the potential curves 4.00–10.00 a.u. The asymptotic point is taken at 400 a.u., the value of this point is used to get the values (De(Ri)) at the point i.


## 4. Results and Discussion

In this section, we discuss our computational result based on our calculations with the linear response adiabatic TDDFT module in  Dirac-Package. Our main concern will be (beside the correctness of our computational result) to compare the behavior of different density functional approximations (and in comparison to other methods) to draw conclusions on the performance, the quality, and the validity of the different functional approximations, also in regard to applications to similar systems and possibly enlighten improvements of the DFT approximations in future works. The comparison with the literature values is accompanying our discussion, where works with different computational methods are available and with experimental values as far as available to judge the quality of our result.


### 4.1. Ground State

As already mentioned, the ground-state bond ofZn2 dimer is a mixture of 3/4  Van der Waals   and 1/4 covalent interactions [85] and the DFT can hardly deal with it as seen in Table 3, where the spectroscopic constants of the ground state are given for different density functional approximations. We note that the effect of the basis set size, typically by DFT, is very small clearly seen in Table 3 from PBE values calculated with aug-cc-pVTZ and aug-cc-pVQZ basis set. In Table 3, one sees that a comparable result is obtained by MP2 and srLDAMP2 as expected [43]. Similar to the rare-gas dimers [43], the range-separated DFT improves the DFT result (here LDA) for Zn2 and suitably cure the lack of correct long-range behavior known by pure DFT approximations because the long-range part of the exchange (and the long-range correlation in srLDAMP2) is treated by a wave function method (MP2). However, a crucial point is to determine a suitable value of the rage-separation parameter. Generally, a suitable range for this parameter is 0.2–0.5 a.u., for details and indepth discussion see [43] and the references therein. DFT approximations and CAMB3LYP, as well as srLDAMP2, do not yield a satisfactory result. Looking at the LDA, we see that the correction of the LDA by srLDA-MP2 is large; however, the improvement gives no advantage over the MP2 as they have similar computational coast. Dramatically behave the long-range corrected PBE0 and the hybrid functionals BLYP and B3LYP (contain a fixed fraction of exact HF-exchange only), they yield a dissociative ground state. BP86 is the only functional with accurate dissociation energy value, but its Re and ωe are not helpful. Although CAMB3LYP gives the best Re value comparison to experiment, this is not sufficient as the bond energy and vibrational frequency are not helpful. It is worthwhile to mention at this point that CAMB3LYP gives the correct asymptotic behavior for the excited states, see Figure 2, in contrast to pure (LDA, PBE, BPW91, BP86, …), long-range corrected (PBE0,GARC-PBE0) or hybrid (BLYP, B3LYP) DFTs, as seen in Figures 2 and 3. Whether this means that CAMB3LYP potential curves has a correct shape (in all regions) is difficult to say at the moment. The shape of the potential curve is an important feature for the DFT accuracy as noted by Grüning et al. [38].Table 3
Ground-state1Σg+ of  Zn2  dimer.
R
e(Å)
ω
e (cm−1)
D
e (eV)
exp1
25.7
0.034
exp2
4.19
25.9
0.035
HF-MP2Q
3.611
29
0.049
srLDAMP2Q
3.445
31
0.0459
PBEQ
3.157
48
0.678
PBE
3.156
49
0.683
PBE0
diss
diss
diss
BPW91
3.225
41
0.0154
BP86
3.181
46
0.036
BLYP
diss
diss
diss
B3LYP
diss
diss
diss
GRAC-PBE0
3.338
40.0
0.045
CAMB3LYP
4.219
11
0.001
LDA
2.846
85
0.225
a
3.959
22
0.024
b
3.96
22.5
0.030
c1
4.03
20.4
0.0205
c2
4.03
20.4
0.0205
pw using aug-cc-pVTZ basis set and SFH.Qaug-cc-pVQZ basis set, for PBE, HF-MP2 and srLDAMP2 (NR with parameter μ=0.5), see text:  1[86];  2[85].  a[12] using CCSD(T) in pseudopotential.  b[17] using NR-CCSD(T). c1[13] CCSD(T) with 4-comp. DCH.  c2[13] CCSD(T) with SFH.


### 4.2. Excited States

The excited states shown in the pw are given in Table2, where n=4 for Zn atom. The results are given in the Tables 5–8. We first discuss the lowest 8 states given in the Tables 5–8, then we proceed to discuss the higher states given in Table 8.At first we compare for PBE functional a 4-component and spin-free result for the four lowest states calculated in aug-cc-pVTZ basis set and demonstrate that SFH describes accurately the main relevant contributions of the relativistic effects. As seen in Table4, the difference between SFH and 4-components DCH is rather small. To see the difference and the splitting in the 4 component precisely De is given in cm-1. The splitting is very small or negligible clearly seen in Figure 1, where we compare visually the 8 lowest states of PBE functional using SFH and the corresponding 16 lowest excited states using 4-component DCH. We note that the CCSD(T) result of [13] for the ground state (see Table 3) using SFH and 4-components DCH confirms our result.Table 4
Comparison between SFH (NR state assignment) and 4-component DCH of the spectroscopic constant. AboveRe (Å), middle ωe (cm−1), and below De (cm−1), calculated with PBE functional and aug-cc-pVTZ basis set. For 4-component states assignment gerade, ungerade follow the symmetry of state in the first line.
3
Π
g
3
Σ
u
+
3
Π
u
3
Σ
g
+
SFH
2.347
2.534
4.795
4.79
4-c.0g-,0u-
2.345
—
4.874
—
4-c.0g+,0u+
2.345
—
4.480
—
4-c.0u,0g
—
2.534
—
4.553
4-c.1g(1u)
2.347
2.534
4.625
4.574
4-c.2g,2u
2.349
—
4.945
—
SFH
219
172
7
27
4-c.0g-,0u-
220
—
6
—
4-c.0g+,0u+
220
—
13
—
4-c.0u,0g
—
172
—
33
4-c.1g(1u)
219
172
13
34
4-c.2g,2u
219
—
8
—
SFH
13097
10870
52
405
4-c.0g-,0u-
12934
—
52
—
4-c.0g+,0u+
13130
—
417
—
4-c.0u,0g
—
10486
—
533
4-c.1g(1u)
12906
10680
235
550
4-c.2g,2u
13068
—
53
—Table 5
Bond lengthsRe (Å) of the lowest states corresponding to the lowest two asymptotes.
Method
3
Π
g
3
Σ
u
+
3
Π
u
3
Σ
g
+
1
Π
g
1
Σ
u
+
1
Π
u
1
Σ
g
+
PQ
2.345
2.532
4.254
4.765
2.350
2.596
4.735
2.573
PT
2.347
2.534
4.795
4.79
2.351
2.602
4.744
2.592
W91T
2.343
2.517
diss
4.546
2.347
2.621
dis
5.158
P0T
2.358
2.517
5.046
4.517
2.351
2.631
2.715
2.594
GP0T
2.356
2.522
diss
5.806
2.345
2.780
2.929
4.755
CB3LT
2.343
2.489
diss
diss
2.327
2.613
2.637
2.572
B3LT
2.371
2.566
diss
5.525
2.366
2.655
2.807
2.624
BLT
2.371
2.587
diss
4.882
2.376
2.648
diss
2.639
B86T
2.337
2.534
diss
4.583
2.341
2.611
4.647
5.370
LDAT
2.265
2.454
2.764
4.364
2.267
2.485
2.702
5.414
[17]a
2.33
2.48
3.99
diss
2.30
2.64
2.40
2.74
[17]b
2.35
2.50
4.11
diss
2.33
2.69
2.42
2.92
[87]c
2.41
2.70
diss
diss
2.33
3.22
2.40
3.05
[19]d
2.38
2.59
4.36
diss
2.38
2.64
2.65f
2.65f
[88]d
2.53
2.74
diss
—
2.51
2.97
2.64
3.07
[89]d
2.56
2.70
diss
diss
2.48
2.92
2.64
—
[90]e
2.372
2.53
exp
—
—
4.49g
—
—
3.0g
—
—
TPresent work calculated with aug-cc-pVTZ and  Qwith aug-cc-pVQZ basis set. P, W91, P0, GP0, B86, BL, B3L, and CB3L denote PBE, BPW91, PBE0, GRAC-PBE0, BP86, BLYP, B3LYP, and CAMB3LYP, respectively.  aWith DK-CASPT2.  bWith DK-MRACPF.  cWith CI. dWith MRCI. eWith CCSD(T). fValue are ca. gFrom [85], for 3Πu [91] gives the value 3.30.Table 6
Vibrational frequenciesωe (cm−1) of the lowest states corresponding to the lowest two asymptotes.
Method
3
Π
g
3
Σ
u
+
3
Π
u
3
Σ
g
+
1
Π
g
1
Σ
u
+
1
Π
u
1
Σ
g
+
PQ
220
173
9
28
219
136
12
142
PT
219
172
7
27
219
135
13
135
W91T
218
177
diss
33
220
129
diss
21
P0T
215
182
9
11
223
135
116
146
GP0T
216
181
diss
11
226
106
78
38
CB3LT
220
189
diss
diss
232
139
137
150
B3LT
211
167
diss
13
215
126
90
139
BLT
210
157
diss
27
207
122
diss
115
B86T
222
172
diss
37
222
131
14
29
LDAT
247
189
85
45
247
160
89
26
[17]a
231
200
23
diss
250
131
211
58
[17]b
220
208
32
diss
244
121
205
104
[87]c
211
169
diss
diss
212
77
175
112
[88]d
192
175
—
diss
210
134
178
—
[89]d
175
150
diss
diss
202
107
166
104
exp
223
±
5
e
161
±
5
f
20.3
±
0
.
2
g
—
—
122
±
1
0
h
148
±
6
i
—
For the acronyms, see Table5. T,Qas in Table 5.  aWith DK-CASPT2.  bWith DK-MRACPF.  cWith CI.  dWith MRCI.  eFrom [92].  fFrom [93].  gFrom [85].  hFrom [94].  iFrom [95].Table 7
Dissociation energies (eV) of the lowest states corresponding to the lowest two asymptotes.
Method
3
Π
g
3
Σ
u
+
3
Π
u
3
Σ
g
+
1
Π
g
1
Σ
u
+
1
Π
u
1
Σ
g
+
PQ
1.626
1.347
0.0065
0.050
1.703
0.579
0.0180
0.112*
PT
1.624
1.348
0.0065
0.050
1.698
0.572
0.0175
0.08*
W91T
1.423
1.23
diss
0.031
1.654
0.541
diss
0.027
P0T
1.481
1.332
0.0034
0.0031
2.387
1.247
0.413
0.10
GP0T
1.43
1.316
diss
0.0148
2.385
1.111
0.270
0.279
CB3LT
1.436
1.281
diss
diss
2.298
1.126
0.099
0.292*
B3LT
1.45
1.189
diss
0.0033
2.226
1.125
0.393
0.148*
BLT
1.514
1.181
diss
0.468
1.426
0.361
diss
0.542
B86T
1.593
1.312
diss
0.0673
1.688
0.546
0.008
0.058
LDAT
2.119
1.704
1.456
1.902
2.089
0.798
0.145
0.788
[17]a
1.502
1.225
0.026
diss
2.713
1.189
0.734
0.60
[17]b
1.457
1.204
0.110
diss
2.694
1.292
0.718
0.204
[87]c
0.91
0.90
diss
diss
2.35
0.71
—
—
[19]d
1.21
0.95
0.016
diss
2.26
1.12
0.63
0.32
[88]d
1.10
0.98
—
diss
2.43
1.13
0.66
—
[89]d
1.05
0.87
diss
diss
2.42
1.06
0.83
0.44
[90]e
1.41
1.21
—
—
—
—
—
—
exp
—
—
0.027
f
—
—
1.117
g
—
—
For the acronyms, see Table5.  T,QAs in Table 5. *See text. aWith DK-CASPT2. bWith DK-MRACPF. cWith CI.  d With MRCI.  eWith CCSD(T).  fFrom [96].  gFrom [94] (1.117±0.025), whereas [91] gives the value 1.30.Table 8
Higher states corresponding to higher asymptotes see Table2 and text.
State
R
e (Å)
ω
e (cm-1)
D
e (eV)
CB3L
P0
GP0
B3L
W91
B86
CB3L
P0
GP0
B3L
W91
B86
CB3L
P0
GP0
B3L
W91
B86
3
Σ
u
+
2.527
2.546
2.711
2.578
2.531
2.532
168
164
115
150
163
160
0.914
0.938
0.174
0.636*
0.555
0.644*
3
Σ
g
+
2.737
2.769
5.772
2.802
2.71
2.714
185
196
23
168
193
186
0.533
0.728
0.596
0.421
0.118
0.094
1
Σ
u
+
2.60
2.630
2.787
2.679
2.622
2.605
149
142
92
120
134
140
0.839
0.677
0.231
0.583
0.513
0.539
1
Σ
g
+
3.444
3.388
8.434
3.449
3.256
3.21
174
146
19
118
131
139
0.339
0.333
0.383
0.097
0.152
0.153
3
Π
u
2.919
3.080
3.162
3.352
3.323
3.451
99
82
72
59
51
45
1.416
0.95
0.90
0.646*
0.039
0.040*
3
Π
g
2.487
2.504
4.748
2.524
2.491
2.485
178
174
41
163
171
172
1.140
0.434
0.635
0.213*
0.143
0.20
3
Σ
u
+
2.519
2.532
diss
2.551
2.546
2.506
172
171
diss
158
166
164
0.905
0.270
diss
0.515
0.482
0.480
3
Σ
g
+
2.569
2.583
diss
2.603
2.513
2.563
153
150
diss
145
140
150
0.247
0.158*
diss
0.163*
0.150
0.157
1
Π
u
3.650
5.750
6.209
9.026
diss
diss
123
14
22
12
diss
diss
1.50
0.483
0.486
0.274
diss
diss
1
Π
g
2.459
2.482
6.317
2.495
2.472
2.465
190
184
26
174
180
182
1.417
0.344
0.482
0.46*
0.43
0.393
1
Σ
u
+
2.534
2.555
diss
2.585
2.537
2.533
169
167
diss
155
162
159
1.125
0.302*
diss
0.50*
0.561
0.560
1
Σ
g
+
2.704
2.682
diss
4.237
2.616
2.583
281
288
diss
296
244
210
0.517*
0.298
diss
0.218*
0.146
0.165
All values with SFH and aug-cc-pVTZ basis set. For the acronyms, see Table5. *See text.(a) Zn2 PBE functional, with SFH (left) ground state (lowest curve) and 8 lowest excited state (corresponding to the two asymptotes (4s2S1 + 4s4pP1) lower ones, and (4s2S1 + 4s4pPv) upper ones. And (b) accordingly the 16 excited states with the same asymptotes using 4-component DCH. Numbering in brackets shows the correspondence between states of (a) and (b).
(a)
(b)Zn
2, CAMB3LYP (a) and B3LYP (b) the 20 lowest states with SFH, corresponding to the asymptotes (from below) 4s2S1 + 4s4pP3, 4s2S1 + 4s4pP1, 4s2S1 + 4s5sS3, 4s2S1 + 4s5sS1, 4s2S1 + 4s5pP3, and 4s2S1 + 4s5pP1, respectively. Note some of the upper curves of B3LYP show incorrect asymptotes.
(a)
(b)Zn
2, PBE0 (a) and GRAC-PBE0 (b) the 20 lowest states with SFH, corresponding to the same asymptotes as in Figure 2. Note that some of the upper curves show an incorrect asymptotes, compare CAMB3LYP Figure 2 and see text.
(a)
(b)In Figure2, we show the 20 lowest excited states corresponding to the 6 asymptotes given in Table 2, for the CAMB3LYP and B3LYP functionals. The overall behavior in Figure 2 for CAMB3LYP is satisfactory, it shows a better behavior for all states, and the states follow (at least) qualitatively to the correct asymptotes. In contrast to the B3LYP as seen in Figure 2(b), where similar result is obtained for all other functionals used in this work. These functionals show an incorrect asymptotic limit and only for the lowest 8 states give the correct (two) asymptotes, whereas most of the higher states follow to a wrong asymptotic limit. This is somehow unexpected since B3LYP includes a (fixed) fraction of exact exchange.In Figure3, a second example is presented for PBE0 and GARC-PBE0. GARC-PBE0 is supposed to give a better result than PBE0, but for Zn2 dimer it does not show a correct description for the higher excited states. Indeed it is well known that pure DFT has incorrect long-range behavior which is the key point behind the range-separated DFT. It is clearly from this result that the separation of the two-electron interaction in short- and long-range parts as done in range-separated DFT like CAMB3LYP offers an advantage by treating the long-range part with a wave function method incorporating a suitable parametric amount of exact exchange. That only CAMB3LYP shows a better or a correct long-range behavior does not mean generally that a range-separated functional describes the excited states better in the short-range (or mid-range) region; however, its accuracy is satisfactory even it fails for the ground state (see Table 3) rather due to the lack of long-range correlation (in HF correlation is not present) important for dispersion interaction.Obviously, a crucial point in calculating the excited states in  TDDFT   is that the most of the DFT approximations are semilocal, the long-range interaction is incorrectly described, consequently a disturbed potential curves is obtained, especially near the avoiding crossing point where the disturbed curves show enhanced effects. This can be clearly seen for the1Σg+,3Σg+, and 1Πu+ in Figure 4. For CAMB3LYP, we see every two states of the same symmetry push each other away and later both follow to the correct limit. For PBE0, as an example, the avoiding crossing is clear for 1Σg+ and 3Σg+ states but not for 1Πu+, most likely because it is disturbed by the incorrect long-range behavior. Similar behavior to PBE0 was found in all other DFT approximations used in this work, that is, an incorrect long-range behavior, with (or leading to) an incorrect asymptotic limit (and a disturbed avoiding crossing) is responsible for incorrect description of the higher excited states. We will discuss the accuracies in detail in the next sections.Zn
2, spinfree Hamiltonian avoiding crossing, CAMB3LYP (a) and PBE0 (b) between the two Σg+1 corresponding to the asymptotes (4s2S1 + 4s4pP1) and (4s2S1 + 4s5sS1), the two πu1 corresponding to the asymptotes (4s2S1 + 4s4pP1) and (4s2S1 + 4s5pP1); the two Σg+3 states corresponding to the asymptotes (4s2S1 + 4s4pP3) and (4s2S1 + 4s5sS3). The highest Σg+1 is corresponding to the asymptotes (4s2S1 + 4s5pP1), see text.
(a)
(b)


#### 4.2.1. Lowest 8 Excited States

In Tables5–7, we give the evaluated spectroscopic constants for the lowest 8 excited states of Zn2 using  TDDFT, SFH, and aug-cc-pVTZ basis set. The lowest 8 excited states 3Πg; 3Πu; 3Σg+; 3Σu+ and 1Πg; 1Πu; 1Σg+; 1Σu+ are corresponding to the Atom((4s2)1S) + Atom((4s4p)3P) and Atom((4s2)1S) + Atom((4s4p)1P), respectively.First, we look at the PBE values using aug-cc-pVTZ basis set and aug-cc-pVQZ basis set. As we see from Tables3–5, the basis effect is small and only about 2*10-3 Å  for Re, about 1 unit for ωe and between 2–6 meV in De. Following this we conclude that the SFH (see Table 4) with aug-cc-pVTZ basis set enable us to calculate the excited states of zinc dimer accurately. Our result is sufficiently accurate to compare with experimental values, wave function methods and compare the behavior of different functional approximations with each other for this dimer.(a) The Lowest States3Πg, 3Πu, 3Σg+, 3Σu+
Looking at the Tables5–7, we see immediately that the best result is obtained for these states. For the lowest two state 3Πg, 3Σu+, all functionals give excellent agreement with wave function results giving in the literature, for example, [17] or the experimental value of ωe, although the agreement for the first excited state, 3Πg, is more pronounced. Recently, Determan et al. [90] have published accurate result for these two states using CCSD(T) and some density functional approximations, the excellent agreement with our values confirms our result. This is not surprising since these states are well bound and largely covalent in contrast to the ground state; moreover, the most known DFT approximations are more or less capable to describe (strong) covalent bonding due to its largely localized character in the bond region. It is also noticeable that all DFTs show for the eight lowest states asymptotically a correct behavior and the correct (two) asymptote, see Figures 2 and 3. For the lowest two states 3Πg, 3Σu+, only LDA strongly underestimates the dissociation energy and gives short bond lengths and large ωe's. PBE gives larger bond energy for both states, likewise BP86 for the first one. BLYP and PBE0 give smaller values for ωe. For Re all these approximations give a similar result. For the next lowest two states, 3Πu, 3Σg+, the situation is somehow complicated. For 3Πu the experimental value shows a weak bound state, whereas wave function methods show different results, likewise in the DFT. PBE and PBE0 describe it as a weak bound state, but apart from LDA all other DFTs give a dissociative state. Whereas for the 3Σu+ only CAMB3LYP shows a dissociative state in an agreement with the wave function methods. This is a first hint that CAMB3LYP gives a better long-range behavior and correct asymptotic limit for higher states than the other DFTs shown in the present work. This can be attributed to the fact that for high-quality response properties it is of primary importance for the potential curve to be accurate in the shape, rather than the condition to be met of being a functional derivative of a given density functional for the exchange-correlation energy [38]. For higher states, both the long-range behavior and the asymptotic limit in pure DFTs are incorrect and thus the shape of potential curves. BLYP gives De≈0.47eV for 3Σg+ which somehow large comparing to other functional. The state 3Σg+ (Atom(4s+4s)1S + Atom(4s+4p1)3P) shows a hump around 2.5 Å  clearly seen in Figure 4 due to an avoiding crossing with the higher state 3Σg+ (Atom(4s+4s)1S + Atom(4s+5p1)3P), the later is well bound (see Table 8) and shows a small hump around 2.2 Å  (hardly seen in Figure 4) presumably due to an avoiding crossing with a more higher state of the same symmetry.(b) The States1Πg, 1Πu, 1Σg+, 1Σu+
From Tables5–7, we again see a good agreement, especially for 1Πg and 1Σu+, between our result and the results of the wave function methods, where the agreement is less pronounced than the lowest two states. For 1Πg and 1Σu+ bond lengths, apart from LDA, all functionals give comparable results. For the vibrational frequencies, BLYP and B3LYP give smaller values, for 1Σu+ this is in excellent agreement with the experimental value of [94] or the value of [17]. CAMB3LYP gives the largest value of ωe. For the dissociation energy De, B3LYP, CAMB3LYP, PBE0, and GRAC-PBE0 give reasonable values with a good agreement with the experiment for 1Σu+. This remarkable result could be a hint that these three functionals have a correct mid-range behavior. From the agreement with the experiment and the wave function values, one concludes that the values of 1Πg of B3LYP, CAMB3LYP, PBE0, and GRAC-PBE0 should be close to the experiment. Next, we look at the two states 1Πu, 1Σg+, as mentioned above in Figure 4, these two states have avoided crossing with higher lying states of the same symmetry. From the tables, we now see a less agreement with the wave function method, and the lack of experimental values makes it more difficult to judge the result. If we take the values of [17], as a reference we see that reasonable DFTs values show larger bond lengths, smaller vibrational frequencies for 1Πu and for 1Σg+ vice versa for the most of the functionals. For 1Πu, the dissociation energies are smaller than the reference value. For 1Σg+, the obtained bond energy values for some functionals denoting the depth of the minimum (marked with “*”) relating to the shallowest point after the minimum, otherwise the incorrect asymptotic point will show a dissociative state, which of course an artifact of the (quantitatively) incorrect tail of the potential curve. We have seen in Figure 2 that CAMB3LYP has asymptotically a correct behavior specially for the higher states; however, it is quantitatively questionable and for some states seems to be inaccurate. In such cases, the spectroscopic constants are calculated relative to the shallowest point after the minimum, and not to the asymptotic point. This yields approximately the same Re and ωe, but the obtained value De will be definitely shallower than, or approximately equal to a value De relating to the “correct” asymptotic point. We note that all values marked with an “*” in Tables 7-8 are obtained this way. For the 1Σg+ state, we see from Table 7 that CAMB3LYP has a good agreement with [19], likewise B3LYP with [17], whereas BLYP shows an agreement with DK-CASPT2 value of [17], but to conclude we see that the result(s) of 1Σg+ are widely distributed; furthermore, the lack of any experimental value makes the situation more difficult.


#### 4.2.2. Higher Excited States

To deal with more higher excited states is difficult because of the above-mentioned reasons. Available approximations do not describe the long-range behavior correctly and/or fail to offer the correct asymptotic limit or predict it accurately [97]. We will discuss the higher molecular states given in Table 6 corresponding to the last four asymptotes (3–6) in Table 2. The result is given for the functionals BPW91 and BP86 (pure), B3LYP (hybrid), CAMB3LYP (range-separated), PBE0 (long-range corrected), and its gradient corrected one GRAC-PBE0. GRAC is an interpolation scheme, it is an asymptotic correction and supposed to be able to deal with higher excited states [37, 38]. The pw shows that the best result is obtained for CAMB3LYP and a comparable result is obtained for PBE0. Indeed strictly only CAMB3LYP was able to deal with higher excited states, it shows (at least qualitatively) the correct asymptotic as can be clearly seen in Figure 2. Other functionals do not show a correct asymptotic behavior as expected [37], including the ones for which no data shown in Table 8. B3LYP is given in Figure 2 as an example, it mixes the asymptotic for higher states with lower states. Our conclusion based on analyzing the data of all functionals and comparing them with each other. It is clear that lacking to the correct long-range behavior is primarily the origin of the problem, CAMB3LYP is able to cure this although not accurately, the question is why other corrections like GRAC does not have the expected improvement? At one side important is the nonlocal part of exact exchange which improves the situation considerably when the two-electron interaction is separated in short- and long-range part such as in CAMB3LYP, and we notice that there is no long-range correlation present in CAMB3LYP because HF offers only (nonlocal) exchange. Another point is the wrong long-range behavior of the response function [72, 77] caused by the incorrect long-range behavior of the density functional approximation is more crucial than it might be believed. This is supported by the fact that the spatial nonlocality of fxc is strongly frequency-dependent [98], in [98] Tokatly and Pankratov argued that not only any static approximation but also any LDA-based dynamic approximation (including any gradient corrections) for fxc cannot provide consistent result. To my best knowledge, there is no calculated or experimental result reported for any of the higher states given in Table 8, this makes the situation more difficult to analyze and be clarified. In Table 8 surprisingly we see that PBE0 gives a better result for higher excited states than its asymptotic corrected one GRAC-PBE0 and better than B3LYP, BP86, or BPW91. Furthermore, it gives for all states a comparable result to CAMB3LYP for Re and ωe. This supports our view and stress the importance of the long-range correction. It is a clear evident that PBE0 has a correct shape in inner part of the potential curve, and only its asymptotic part (tail of the potential curve) is incorrect, unfortunately the applied correction of GRAC is not good. As seen in Table 8 our next four states, 3Σu+, 3Σg+, and 1Σu+, 1Σg+ corresponding to Atom((4s2)1S + Atom((4s5s)3S), and Atom((4s2)1S + Atom((4s5s)1S), have more or less a similar result for all functionals, only GRAC-PBE0 shows unexplainable result, since it is supposed to show asymptotically a better behavior. We think that the CAMB3LYP result is the most correct one although it might be not satisfactory accurate. It is worthwhile to mention that states with avoiding crossing get a second shallow minimum after the avoiding crossing at large internuclear distances, this is not reported and only the first minimum is presented. Next, we look to the states 3πu, 3πg, 3Σu+, 3Σg+ corresponding to the Atom((4s2)1S + Atom((4s5p)3P). Here, we see that the result is distributed, BPW91, BP86, and B3LYP show similar results, whereas GRAC-PBE0 differs considerably from all approximations given in Table 8. PBE0 result is close to CAMB3LYP when looking to Re and ωe, but its De values are different clearly due to its incorrect asymptotic limit. The last states treated in this work 1πu, 1πg, 1Σu+, 1Σg+ are corresponding to the Atom((4s2)1S + Atom((4s5p)1P). The results of 1πu are puzzling and presumably only the values of CAMB3LYP are reasonable, whereas for 1πg all functional apart from GRAC-PBE0 give comparable values for ωe and Re, which could be a hint that these values are reasonable. 1Σu+, and 1Σg+ follow the general trend that PBE0 result is close to CAMB3LYP. BPW91, BP86, and B3LYP show a similar result, GRAC-PBE0 shows unexplainable result.The general conclusion of this section is that CAMB3LYP gives the best result due to its better treatment of the long-range part of the two-electron interaction and its asymptotically better behavior (tail of the potential curve) apparently due to including a suitable amount of exact exchange, PBE0 gives a comparable result, the main problem here is the tail of the potential curve. BPW91, BP86, and B3LYP are less satisfactory but still show acceptable result, whereas (most likely) the result of GRAC-PBE0 is not useful.


## 4.1. Ground State

As already mentioned, the ground-state bond ofZn2 dimer is a mixture of 3/4  Van der Waals   and 1/4 covalent interactions [85] and the DFT can hardly deal with it as seen in Table 3, where the spectroscopic constants of the ground state are given for different density functional approximations. We note that the effect of the basis set size, typically by DFT, is very small clearly seen in Table 3 from PBE values calculated with aug-cc-pVTZ and aug-cc-pVQZ basis set. In Table 3, one sees that a comparable result is obtained by MP2 and srLDAMP2 as expected [43]. Similar to the rare-gas dimers [43], the range-separated DFT improves the DFT result (here LDA) for Zn2 and suitably cure the lack of correct long-range behavior known by pure DFT approximations because the long-range part of the exchange (and the long-range correlation in srLDAMP2) is treated by a wave function method (MP2). However, a crucial point is to determine a suitable value of the rage-separation parameter. Generally, a suitable range for this parameter is 0.2–0.5 a.u., for details and indepth discussion see [43] and the references therein. DFT approximations and CAMB3LYP, as well as srLDAMP2, do not yield a satisfactory result. Looking at the LDA, we see that the correction of the LDA by srLDA-MP2 is large; however, the improvement gives no advantage over the MP2 as they have similar computational coast. Dramatically behave the long-range corrected PBE0 and the hybrid functionals BLYP and B3LYP (contain a fixed fraction of exact HF-exchange only), they yield a dissociative ground state. BP86 is the only functional with accurate dissociation energy value, but its Re and ωe are not helpful. Although CAMB3LYP gives the best Re value comparison to experiment, this is not sufficient as the bond energy and vibrational frequency are not helpful. It is worthwhile to mention at this point that CAMB3LYP gives the correct asymptotic behavior for the excited states, see Figure 2, in contrast to pure (LDA, PBE, BPW91, BP86, …), long-range corrected (PBE0,GARC-PBE0) or hybrid (BLYP, B3LYP) DFTs, as seen in Figures 2 and 3. Whether this means that CAMB3LYP potential curves has a correct shape (in all regions) is difficult to say at the moment. The shape of the potential curve is an important feature for the DFT accuracy as noted by Grüning et al. [38].Table 3
Ground-state1Σg+ of  Zn2  dimer.
R
e(Å)
ω
e (cm−1)
D
e (eV)
exp1
25.7
0.034
exp2
4.19
25.9
0.035
HF-MP2Q
3.611
29
0.049
srLDAMP2Q
3.445
31
0.0459
PBEQ
3.157
48
0.678
PBE
3.156
49
0.683
PBE0
diss
diss
diss
BPW91
3.225
41
0.0154
BP86
3.181
46
0.036
BLYP
diss
diss
diss
B3LYP
diss
diss
diss
GRAC-PBE0
3.338
40.0
0.045
CAMB3LYP
4.219
11
0.001
LDA
2.846
85
0.225
a
3.959
22
0.024
b
3.96
22.5
0.030
c1
4.03
20.4
0.0205
c2
4.03
20.4
0.0205
pw using aug-cc-pVTZ basis set and SFH.Qaug-cc-pVQZ basis set, for PBE, HF-MP2 and srLDAMP2 (NR with parameter μ=0.5), see text:  1[86];  2[85].  a[12] using CCSD(T) in pseudopotential.  b[17] using NR-CCSD(T). c1[13] CCSD(T) with 4-comp. DCH.  c2[13] CCSD(T) with SFH.


## 4.2. Excited States

The excited states shown in the pw are given in Table2, where n=4 for Zn atom. The results are given in the Tables 5–8. We first discuss the lowest 8 states given in the Tables 5–8, then we proceed to discuss the higher states given in Table 8.At first we compare for PBE functional a 4-component and spin-free result for the four lowest states calculated in aug-cc-pVTZ basis set and demonstrate that SFH describes accurately the main relevant contributions of the relativistic effects. As seen in Table4, the difference between SFH and 4-components DCH is rather small. To see the difference and the splitting in the 4 component precisely De is given in cm-1. The splitting is very small or negligible clearly seen in Figure 1, where we compare visually the 8 lowest states of PBE functional using SFH and the corresponding 16 lowest excited states using 4-component DCH. We note that the CCSD(T) result of [13] for the ground state (see Table 3) using SFH and 4-components DCH confirms our result.Table 4
Comparison between SFH (NR state assignment) and 4-component DCH of the spectroscopic constant. AboveRe (Å), middle ωe (cm−1), and below De (cm−1), calculated with PBE functional and aug-cc-pVTZ basis set. For 4-component states assignment gerade, ungerade follow the symmetry of state in the first line.
3
Π
g
3
Σ
u
+
3
Π
u
3
Σ
g
+
SFH
2.347
2.534
4.795
4.79
4-c.0g-,0u-
2.345
—
4.874
—
4-c.0g+,0u+
2.345
—
4.480
—
4-c.0u,0g
—
2.534
—
4.553
4-c.1g(1u)
2.347
2.534
4.625
4.574
4-c.2g,2u
2.349
—
4.945
—
SFH
219
172
7
27
4-c.0g-,0u-
220
—
6
—
4-c.0g+,0u+
220
—
13
—
4-c.0u,0g
—
172
—
33
4-c.1g(1u)
219
172
13
34
4-c.2g,2u
219
—
8
—
SFH
13097
10870
52
405
4-c.0g-,0u-
12934
—
52
—
4-c.0g+,0u+
13130
—
417
—
4-c.0u,0g
—
10486
—
533
4-c.1g(1u)
12906
10680
235
550
4-c.2g,2u
13068
—
53
—Table 5
Bond lengthsRe (Å) of the lowest states corresponding to the lowest two asymptotes.
Method
3
Π
g
3
Σ
u
+
3
Π
u
3
Σ
g
+
1
Π
g
1
Σ
u
+
1
Π
u
1
Σ
g
+
PQ
2.345
2.532
4.254
4.765
2.350
2.596
4.735
2.573
PT
2.347
2.534
4.795
4.79
2.351
2.602
4.744
2.592
W91T
2.343
2.517
diss
4.546
2.347
2.621
dis
5.158
P0T
2.358
2.517
5.046
4.517
2.351
2.631
2.715
2.594
GP0T
2.356
2.522
diss
5.806
2.345
2.780
2.929
4.755
CB3LT
2.343
2.489
diss
diss
2.327
2.613
2.637
2.572
B3LT
2.371
2.566
diss
5.525
2.366
2.655
2.807
2.624
BLT
2.371
2.587
diss
4.882
2.376
2.648
diss
2.639
B86T
2.337
2.534
diss
4.583
2.341
2.611
4.647
5.370
LDAT
2.265
2.454
2.764
4.364
2.267
2.485
2.702
5.414
[17]a
2.33
2.48
3.99
diss
2.30
2.64
2.40
2.74
[17]b
2.35
2.50
4.11
diss
2.33
2.69
2.42
2.92
[87]c
2.41
2.70
diss
diss
2.33
3.22
2.40
3.05
[19]d
2.38
2.59
4.36
diss
2.38
2.64
2.65f
2.65f
[88]d
2.53
2.74
diss
—
2.51
2.97
2.64
3.07
[89]d
2.56
2.70
diss
diss
2.48
2.92
2.64
—
[90]e
2.372
2.53
exp
—
—
4.49g
—
—
3.0g
—
—
TPresent work calculated with aug-cc-pVTZ and  Qwith aug-cc-pVQZ basis set. P, W91, P0, GP0, B86, BL, B3L, and CB3L denote PBE, BPW91, PBE0, GRAC-PBE0, BP86, BLYP, B3LYP, and CAMB3LYP, respectively.  aWith DK-CASPT2.  bWith DK-MRACPF.  cWith CI. dWith MRCI. eWith CCSD(T). fValue are ca. gFrom [85], for 3Πu [91] gives the value 3.30.Table 6
Vibrational frequenciesωe (cm−1) of the lowest states corresponding to the lowest two asymptotes.
Method
3
Π
g
3
Σ
u
+
3
Π
u
3
Σ
g
+
1
Π
g
1
Σ
u
+
1
Π
u
1
Σ
g
+
PQ
220
173
9
28
219
136
12
142
PT
219
172
7
27
219
135
13
135
W91T
218
177
diss
33
220
129
diss
21
P0T
215
182
9
11
223
135
116
146
GP0T
216
181
diss
11
226
106
78
38
CB3LT
220
189
diss
diss
232
139
137
150
B3LT
211
167
diss
13
215
126
90
139
BLT
210
157
diss
27
207
122
diss
115
B86T
222
172
diss
37
222
131
14
29
LDAT
247
189
85
45
247
160
89
26
[17]a
231
200
23
diss
250
131
211
58
[17]b
220
208
32
diss
244
121
205
104
[87]c
211
169
diss
diss
212
77
175
112
[88]d
192
175
—
diss
210
134
178
—
[89]d
175
150
diss
diss
202
107
166
104
exp
223
±
5
e
161
±
5
f
20.3
±
0
.
2
g
—
—
122
±
1
0
h
148
±
6
i
—
For the acronyms, see Table5. T,Qas in Table 5.  aWith DK-CASPT2.  bWith DK-MRACPF.  cWith CI.  dWith MRCI.  eFrom [92].  fFrom [93].  gFrom [85].  hFrom [94].  iFrom [95].Table 7
Dissociation energies (eV) of the lowest states corresponding to the lowest two asymptotes.
Method
3
Π
g
3
Σ
u
+
3
Π
u
3
Σ
g
+
1
Π
g
1
Σ
u
+
1
Π
u
1
Σ
g
+
PQ
1.626
1.347
0.0065
0.050
1.703
0.579
0.0180
0.112*
PT
1.624
1.348
0.0065
0.050
1.698
0.572
0.0175
0.08*
W91T
1.423
1.23
diss
0.031
1.654
0.541
diss
0.027
P0T
1.481
1.332
0.0034
0.0031
2.387
1.247
0.413
0.10
GP0T
1.43
1.316
diss
0.0148
2.385
1.111
0.270
0.279
CB3LT
1.436
1.281
diss
diss
2.298
1.126
0.099
0.292*
B3LT
1.45
1.189
diss
0.0033
2.226
1.125
0.393
0.148*
BLT
1.514
1.181
diss
0.468
1.426
0.361
diss
0.542
B86T
1.593
1.312
diss
0.0673
1.688
0.546
0.008
0.058
LDAT
2.119
1.704
1.456
1.902
2.089
0.798
0.145
0.788
[17]a
1.502
1.225
0.026
diss
2.713
1.189
0.734
0.60
[17]b
1.457
1.204
0.110
diss
2.694
1.292
0.718
0.204
[87]c
0.91
0.90
diss
diss
2.35
0.71
—
—
[19]d
1.21
0.95
0.016
diss
2.26
1.12
0.63
0.32
[88]d
1.10
0.98
—
diss
2.43
1.13
0.66
—
[89]d
1.05
0.87
diss
diss
2.42
1.06
0.83
0.44
[90]e
1.41
1.21
—
—
—
—
—
—
exp
—
—
0.027
f
—
—
1.117
g
—
—
For the acronyms, see Table5.  T,QAs in Table 5. *See text. aWith DK-CASPT2. bWith DK-MRACPF. cWith CI.  d With MRCI.  eWith CCSD(T).  fFrom [96].  gFrom [94] (1.117±0.025), whereas [91] gives the value 1.30.Table 8
Higher states corresponding to higher asymptotes see Table2 and text.
State
R
e (Å)
ω
e (cm-1)
D
e (eV)
CB3L
P0
GP0
B3L
W91
B86
CB3L
P0
GP0
B3L
W91
B86
CB3L
P0
GP0
B3L
W91
B86
3
Σ
u
+
2.527
2.546
2.711
2.578
2.531
2.532
168
164
115
150
163
160
0.914
0.938
0.174
0.636*
0.555
0.644*
3
Σ
g
+
2.737
2.769
5.772
2.802
2.71
2.714
185
196
23
168
193
186
0.533
0.728
0.596
0.421
0.118
0.094
1
Σ
u
+
2.60
2.630
2.787
2.679
2.622
2.605
149
142
92
120
134
140
0.839
0.677
0.231
0.583
0.513
0.539
1
Σ
g
+
3.444
3.388
8.434
3.449
3.256
3.21
174
146
19
118
131
139
0.339
0.333
0.383
0.097
0.152
0.153
3
Π
u
2.919
3.080
3.162
3.352
3.323
3.451
99
82
72
59
51
45
1.416
0.95
0.90
0.646*
0.039
0.040*
3
Π
g
2.487
2.504
4.748
2.524
2.491
2.485
178
174
41
163
171
172
1.140
0.434
0.635
0.213*
0.143
0.20
3
Σ
u
+
2.519
2.532
diss
2.551
2.546
2.506
172
171
diss
158
166
164
0.905
0.270
diss
0.515
0.482
0.480
3
Σ
g
+
2.569
2.583
diss
2.603
2.513
2.563
153
150
diss
145
140
150
0.247
0.158*
diss
0.163*
0.150
0.157
1
Π
u
3.650
5.750
6.209
9.026
diss
diss
123
14
22
12
diss
diss
1.50
0.483
0.486
0.274
diss
diss
1
Π
g
2.459
2.482
6.317
2.495
2.472
2.465
190
184
26
174
180
182
1.417
0.344
0.482
0.46*
0.43
0.393
1
Σ
u
+
2.534
2.555
diss
2.585
2.537
2.533
169
167
diss
155
162
159
1.125
0.302*
diss
0.50*
0.561
0.560
1
Σ
g
+
2.704
2.682
diss
4.237
2.616
2.583
281
288
diss
296
244
210
0.517*
0.298
diss
0.218*
0.146
0.165
All values with SFH and aug-cc-pVTZ basis set. For the acronyms, see Table5. *See text.(a) Zn2 PBE functional, with SFH (left) ground state (lowest curve) and 8 lowest excited state (corresponding to the two asymptotes (4s2S1 + 4s4pP1) lower ones, and (4s2S1 + 4s4pPv) upper ones. And (b) accordingly the 16 excited states with the same asymptotes using 4-component DCH. Numbering in brackets shows the correspondence between states of (a) and (b).
(a)
(b)Zn
2, CAMB3LYP (a) and B3LYP (b) the 20 lowest states with SFH, corresponding to the asymptotes (from below) 4s2S1 + 4s4pP3, 4s2S1 + 4s4pP1, 4s2S1 + 4s5sS3, 4s2S1 + 4s5sS1, 4s2S1 + 4s5pP3, and 4s2S1 + 4s5pP1, respectively. Note some of the upper curves of B3LYP show incorrect asymptotes.
(a)
(b)Zn
2, PBE0 (a) and GRAC-PBE0 (b) the 20 lowest states with SFH, corresponding to the same asymptotes as in Figure 2. Note that some of the upper curves show an incorrect asymptotes, compare CAMB3LYP Figure 2 and see text.
(a)
(b)In Figure2, we show the 20 lowest excited states corresponding to the 6 asymptotes given in Table 2, for the CAMB3LYP and B3LYP functionals. The overall behavior in Figure 2 for CAMB3LYP is satisfactory, it shows a better behavior for all states, and the states follow (at least) qualitatively to the correct asymptotes. In contrast to the B3LYP as seen in Figure 2(b), where similar result is obtained for all other functionals used in this work. These functionals show an incorrect asymptotic limit and only for the lowest 8 states give the correct (two) asymptotes, whereas most of the higher states follow to a wrong asymptotic limit. This is somehow unexpected since B3LYP includes a (fixed) fraction of exact exchange.In Figure3, a second example is presented for PBE0 and GARC-PBE0. GARC-PBE0 is supposed to give a better result than PBE0, but for Zn2 dimer it does not show a correct description for the higher excited states. Indeed it is well known that pure DFT has incorrect long-range behavior which is the key point behind the range-separated DFT. It is clearly from this result that the separation of the two-electron interaction in short- and long-range parts as done in range-separated DFT like CAMB3LYP offers an advantage by treating the long-range part with a wave function method incorporating a suitable parametric amount of exact exchange. That only CAMB3LYP shows a better or a correct long-range behavior does not mean generally that a range-separated functional describes the excited states better in the short-range (or mid-range) region; however, its accuracy is satisfactory even it fails for the ground state (see Table 3) rather due to the lack of long-range correlation (in HF correlation is not present) important for dispersion interaction.Obviously, a crucial point in calculating the excited states in  TDDFT   is that the most of the DFT approximations are semilocal, the long-range interaction is incorrectly described, consequently a disturbed potential curves is obtained, especially near the avoiding crossing point where the disturbed curves show enhanced effects. This can be clearly seen for the1Σg+,3Σg+, and 1Πu+ in Figure 4. For CAMB3LYP, we see every two states of the same symmetry push each other away and later both follow to the correct limit. For PBE0, as an example, the avoiding crossing is clear for 1Σg+ and 3Σg+ states but not for 1Πu+, most likely because it is disturbed by the incorrect long-range behavior. Similar behavior to PBE0 was found in all other DFT approximations used in this work, that is, an incorrect long-range behavior, with (or leading to) an incorrect asymptotic limit (and a disturbed avoiding crossing) is responsible for incorrect description of the higher excited states. We will discuss the accuracies in detail in the next sections.Zn
2, spinfree Hamiltonian avoiding crossing, CAMB3LYP (a) and PBE0 (b) between the two Σg+1 corresponding to the asymptotes (4s2S1 + 4s4pP1) and (4s2S1 + 4s5sS1), the two πu1 corresponding to the asymptotes (4s2S1 + 4s4pP1) and (4s2S1 + 4s5pP1); the two Σg+3 states corresponding to the asymptotes (4s2S1 + 4s4pP3) and (4s2S1 + 4s5sS3). The highest Σg+1 is corresponding to the asymptotes (4s2S1 + 4s5pP1), see text.
(a)
(b)


### 4.2.1. Lowest 8 Excited States

In Tables5–7, we give the evaluated spectroscopic constants for the lowest 8 excited states of Zn2 using  TDDFT, SFH, and aug-cc-pVTZ basis set. The lowest 8 excited states 3Πg; 3Πu; 3Σg+; 3Σu+ and 1Πg; 1Πu; 1Σg+; 1Σu+ are corresponding to the Atom((4s2)1S) + Atom((4s4p)3P) and Atom((4s2)1S) + Atom((4s4p)1P), respectively.First, we look at the PBE values using aug-cc-pVTZ basis set and aug-cc-pVQZ basis set. As we see from Tables3–5, the basis effect is small and only about 2*10-3 Å  for Re, about 1 unit for ωe and between 2–6 meV in De. Following this we conclude that the SFH (see Table 4) with aug-cc-pVTZ basis set enable us to calculate the excited states of zinc dimer accurately. Our result is sufficiently accurate to compare with experimental values, wave function methods and compare the behavior of different functional approximations with each other for this dimer.(a) The Lowest States3Πg, 3Πu, 3Σg+, 3Σu+
Looking at the Tables5–7, we see immediately that the best result is obtained for these states. For the lowest two state 3Πg, 3Σu+, all functionals give excellent agreement with wave function results giving in the literature, for example, [17] or the experimental value of ωe, although the agreement for the first excited state, 3Πg, is more pronounced. Recently, Determan et al. [90] have published accurate result for these two states using CCSD(T) and some density functional approximations, the excellent agreement with our values confirms our result. This is not surprising since these states are well bound and largely covalent in contrast to the ground state; moreover, the most known DFT approximations are more or less capable to describe (strong) covalent bonding due to its largely localized character in the bond region. It is also noticeable that all DFTs show for the eight lowest states asymptotically a correct behavior and the correct (two) asymptote, see Figures 2 and 3. For the lowest two states 3Πg, 3Σu+, only LDA strongly underestimates the dissociation energy and gives short bond lengths and large ωe's. PBE gives larger bond energy for both states, likewise BP86 for the first one. BLYP and PBE0 give smaller values for ωe. For Re all these approximations give a similar result. For the next lowest two states, 3Πu, 3Σg+, the situation is somehow complicated. For 3Πu the experimental value shows a weak bound state, whereas wave function methods show different results, likewise in the DFT. PBE and PBE0 describe it as a weak bound state, but apart from LDA all other DFTs give a dissociative state. Whereas for the 3Σu+ only CAMB3LYP shows a dissociative state in an agreement with the wave function methods. This is a first hint that CAMB3LYP gives a better long-range behavior and correct asymptotic limit for higher states than the other DFTs shown in the present work. This can be attributed to the fact that for high-quality response properties it is of primary importance for the potential curve to be accurate in the shape, rather than the condition to be met of being a functional derivative of a given density functional for the exchange-correlation energy [38]. For higher states, both the long-range behavior and the asymptotic limit in pure DFTs are incorrect and thus the shape of potential curves. BLYP gives De≈0.47eV for 3Σg+ which somehow large comparing to other functional. The state 3Σg+ (Atom(4s+4s)1S + Atom(4s+4p1)3P) shows a hump around 2.5 Å  clearly seen in Figure 4 due to an avoiding crossing with the higher state 3Σg+ (Atom(4s+4s)1S + Atom(4s+5p1)3P), the later is well bound (see Table 8) and shows a small hump around 2.2 Å  (hardly seen in Figure 4) presumably due to an avoiding crossing with a more higher state of the same symmetry.(b) The States1Πg, 1Πu, 1Σg+, 1Σu+
From Tables5–7, we again see a good agreement, especially for 1Πg and 1Σu+, between our result and the results of the wave function methods, where the agreement is less pronounced than the lowest two states. For 1Πg and 1Σu+ bond lengths, apart from LDA, all functionals give comparable results. For the vibrational frequencies, BLYP and B3LYP give smaller values, for 1Σu+ this is in excellent agreement with the experimental value of [94] or the value of [17]. CAMB3LYP gives the largest value of ωe. For the dissociation energy De, B3LYP, CAMB3LYP, PBE0, and GRAC-PBE0 give reasonable values with a good agreement with the experiment for 1Σu+. This remarkable result could be a hint that these three functionals have a correct mid-range behavior. From the agreement with the experiment and the wave function values, one concludes that the values of 1Πg of B3LYP, CAMB3LYP, PBE0, and GRAC-PBE0 should be close to the experiment. Next, we look at the two states 1Πu, 1Σg+, as mentioned above in Figure 4, these two states have avoided crossing with higher lying states of the same symmetry. From the tables, we now see a less agreement with the wave function method, and the lack of experimental values makes it more difficult to judge the result. If we take the values of [17], as a reference we see that reasonable DFTs values show larger bond lengths, smaller vibrational frequencies for 1Πu and for 1Σg+ vice versa for the most of the functionals. For 1Πu, the dissociation energies are smaller than the reference value. For 1Σg+, the obtained bond energy values for some functionals denoting the depth of the minimum (marked with “*”) relating to the shallowest point after the minimum, otherwise the incorrect asymptotic point will show a dissociative state, which of course an artifact of the (quantitatively) incorrect tail of the potential curve. We have seen in Figure 2 that CAMB3LYP has asymptotically a correct behavior specially for the higher states; however, it is quantitatively questionable and for some states seems to be inaccurate. In such cases, the spectroscopic constants are calculated relative to the shallowest point after the minimum, and not to the asymptotic point. This yields approximately the same Re and ωe, but the obtained value De will be definitely shallower than, or approximately equal to a value De relating to the “correct” asymptotic point. We note that all values marked with an “*” in Tables 7-8 are obtained this way. For the 1Σg+ state, we see from Table 7 that CAMB3LYP has a good agreement with [19], likewise B3LYP with [17], whereas BLYP shows an agreement with DK-CASPT2 value of [17], but to conclude we see that the result(s) of 1Σg+ are widely distributed; furthermore, the lack of any experimental value makes the situation more difficult.


### 4.2.2. Higher Excited States

To deal with more higher excited states is difficult because of the above-mentioned reasons. Available approximations do not describe the long-range behavior correctly and/or fail to offer the correct asymptotic limit or predict it accurately [97]. We will discuss the higher molecular states given in Table 6 corresponding to the last four asymptotes (3–6) in Table 2. The result is given for the functionals BPW91 and BP86 (pure), B3LYP (hybrid), CAMB3LYP (range-separated), PBE0 (long-range corrected), and its gradient corrected one GRAC-PBE0. GRAC is an interpolation scheme, it is an asymptotic correction and supposed to be able to deal with higher excited states [37, 38]. The pw shows that the best result is obtained for CAMB3LYP and a comparable result is obtained for PBE0. Indeed strictly only CAMB3LYP was able to deal with higher excited states, it shows (at least qualitatively) the correct asymptotic as can be clearly seen in Figure 2. Other functionals do not show a correct asymptotic behavior as expected [37], including the ones for which no data shown in Table 8. B3LYP is given in Figure 2 as an example, it mixes the asymptotic for higher states with lower states. Our conclusion based on analyzing the data of all functionals and comparing them with each other. It is clear that lacking to the correct long-range behavior is primarily the origin of the problem, CAMB3LYP is able to cure this although not accurately, the question is why other corrections like GRAC does not have the expected improvement? At one side important is the nonlocal part of exact exchange which improves the situation considerably when the two-electron interaction is separated in short- and long-range part such as in CAMB3LYP, and we notice that there is no long-range correlation present in CAMB3LYP because HF offers only (nonlocal) exchange. Another point is the wrong long-range behavior of the response function [72, 77] caused by the incorrect long-range behavior of the density functional approximation is more crucial than it might be believed. This is supported by the fact that the spatial nonlocality of fxc is strongly frequency-dependent [98], in [98] Tokatly and Pankratov argued that not only any static approximation but also any LDA-based dynamic approximation (including any gradient corrections) for fxc cannot provide consistent result. To my best knowledge, there is no calculated or experimental result reported for any of the higher states given in Table 8, this makes the situation more difficult to analyze and be clarified. In Table 8 surprisingly we see that PBE0 gives a better result for higher excited states than its asymptotic corrected one GRAC-PBE0 and better than B3LYP, BP86, or BPW91. Furthermore, it gives for all states a comparable result to CAMB3LYP for Re and ωe. This supports our view and stress the importance of the long-range correction. It is a clear evident that PBE0 has a correct shape in inner part of the potential curve, and only its asymptotic part (tail of the potential curve) is incorrect, unfortunately the applied correction of GRAC is not good. As seen in Table 8 our next four states, 3Σu+, 3Σg+, and 1Σu+, 1Σg+ corresponding to Atom((4s2)1S + Atom((4s5s)3S), and Atom((4s2)1S + Atom((4s5s)1S), have more or less a similar result for all functionals, only GRAC-PBE0 shows unexplainable result, since it is supposed to show asymptotically a better behavior. We think that the CAMB3LYP result is the most correct one although it might be not satisfactory accurate. It is worthwhile to mention that states with avoiding crossing get a second shallow minimum after the avoiding crossing at large internuclear distances, this is not reported and only the first minimum is presented. Next, we look to the states 3πu, 3πg, 3Σu+, 3Σg+ corresponding to the Atom((4s2)1S + Atom((4s5p)3P). Here, we see that the result is distributed, BPW91, BP86, and B3LYP show similar results, whereas GRAC-PBE0 differs considerably from all approximations given in Table 8. PBE0 result is close to CAMB3LYP when looking to Re and ωe, but its De values are different clearly due to its incorrect asymptotic limit. The last states treated in this work 1πu, 1πg, 1Σu+, 1Σg+ are corresponding to the Atom((4s2)1S + Atom((4s5p)1P). The results of 1πu are puzzling and presumably only the values of CAMB3LYP are reasonable, whereas for 1πg all functional apart from GRAC-PBE0 give comparable values for ωe and Re, which could be a hint that these values are reasonable. 1Σu+, and 1Σg+ follow the general trend that PBE0 result is close to CAMB3LYP. BPW91, BP86, and B3LYP show a similar result, GRAC-PBE0 shows unexplainable result.The general conclusion of this section is that CAMB3LYP gives the best result due to its better treatment of the long-range part of the two-electron interaction and its asymptotically better behavior (tail of the potential curve) apparently due to including a suitable amount of exact exchange, PBE0 gives a comparable result, the main problem here is the tail of the potential curve. BPW91, BP86, and B3LYP are less satisfactory but still show acceptable result, whereas (most likely) the result of GRAC-PBE0 is not useful.


## 4.2.1. Lowest 8 Excited States

In Tables5–7, we give the evaluated spectroscopic constants for the lowest 8 excited states of Zn2 using  TDDFT, SFH, and aug-cc-pVTZ basis set. The lowest 8 excited states 3Πg; 3Πu; 3Σg+; 3Σu+ and 1Πg; 1Πu; 1Σg+; 1Σu+ are corresponding to the Atom((4s2)1S) + Atom((4s4p)3P) and Atom((4s2)1S) + Atom((4s4p)1P), respectively.First, we look at the PBE values using aug-cc-pVTZ basis set and aug-cc-pVQZ basis set. As we see from Tables3–5, the basis effect is small and only about 2*10-3 Å  for Re, about 1 unit for ωe and between 2–6 meV in De. Following this we conclude that the SFH (see Table 4) with aug-cc-pVTZ basis set enable us to calculate the excited states of zinc dimer accurately. Our result is sufficiently accurate to compare with experimental values, wave function methods and compare the behavior of different functional approximations with each other for this dimer.(a) The Lowest States3Πg, 3Πu, 3Σg+, 3Σu+
Looking at the Tables5–7, we see immediately that the best result is obtained for these states. For the lowest two state 3Πg, 3Σu+, all functionals give excellent agreement with wave function results giving in the literature, for example, [17] or the experimental value of ωe, although the agreement for the first excited state, 3Πg, is more pronounced. Recently, Determan et al. [90] have published accurate result for these two states using CCSD(T) and some density functional approximations, the excellent agreement with our values confirms our result. This is not surprising since these states are well bound and largely covalent in contrast to the ground state; moreover, the most known DFT approximations are more or less capable to describe (strong) covalent bonding due to its largely localized character in the bond region. It is also noticeable that all DFTs show for the eight lowest states asymptotically a correct behavior and the correct (two) asymptote, see Figures 2 and 3. For the lowest two states 3Πg, 3Σu+, only LDA strongly underestimates the dissociation energy and gives short bond lengths and large ωe's. PBE gives larger bond energy for both states, likewise BP86 for the first one. BLYP and PBE0 give smaller values for ωe. For Re all these approximations give a similar result. For the next lowest two states, 3Πu, 3Σg+, the situation is somehow complicated. For 3Πu the experimental value shows a weak bound state, whereas wave function methods show different results, likewise in the DFT. PBE and PBE0 describe it as a weak bound state, but apart from LDA all other DFTs give a dissociative state. Whereas for the 3Σu+ only CAMB3LYP shows a dissociative state in an agreement with the wave function methods. This is a first hint that CAMB3LYP gives a better long-range behavior and correct asymptotic limit for higher states than the other DFTs shown in the present work. This can be attributed to the fact that for high-quality response properties it is of primary importance for the potential curve to be accurate in the shape, rather than the condition to be met of being a functional derivative of a given density functional for the exchange-correlation energy [38]. For higher states, both the long-range behavior and the asymptotic limit in pure DFTs are incorrect and thus the shape of potential curves. BLYP gives De≈0.47eV for 3Σg+ which somehow large comparing to other functional. The state 3Σg+ (Atom(4s+4s)1S + Atom(4s+4p1)3P) shows a hump around 2.5 Å  clearly seen in Figure 4 due to an avoiding crossing with the higher state 3Σg+ (Atom(4s+4s)1S + Atom(4s+5p1)3P), the later is well bound (see Table 8) and shows a small hump around 2.2 Å  (hardly seen in Figure 4) presumably due to an avoiding crossing with a more higher state of the same symmetry.(b) The States1Πg, 1Πu, 1Σg+, 1Σu+
From Tables5–7, we again see a good agreement, especially for 1Πg and 1Σu+, between our result and the results of the wave function methods, where the agreement is less pronounced than the lowest two states. For 1Πg and 1Σu+ bond lengths, apart from LDA, all functionals give comparable results. For the vibrational frequencies, BLYP and B3LYP give smaller values, for 1Σu+ this is in excellent agreement with the experimental value of [94] or the value of [17]. CAMB3LYP gives the largest value of ωe. For the dissociation energy De, B3LYP, CAMB3LYP, PBE0, and GRAC-PBE0 give reasonable values with a good agreement with the experiment for 1Σu+. This remarkable result could be a hint that these three functionals have a correct mid-range behavior. From the agreement with the experiment and the wave function values, one concludes that the values of 1Πg of B3LYP, CAMB3LYP, PBE0, and GRAC-PBE0 should be close to the experiment. Next, we look at the two states 1Πu, 1Σg+, as mentioned above in Figure 4, these two states have avoided crossing with higher lying states of the same symmetry. From the tables, we now see a less agreement with the wave function method, and the lack of experimental values makes it more difficult to judge the result. If we take the values of [17], as a reference we see that reasonable DFTs values show larger bond lengths, smaller vibrational frequencies for 1Πu and for 1Σg+ vice versa for the most of the functionals. For 1Πu, the dissociation energies are smaller than the reference value. For 1Σg+, the obtained bond energy values for some functionals denoting the depth of the minimum (marked with “*”) relating to the shallowest point after the minimum, otherwise the incorrect asymptotic point will show a dissociative state, which of course an artifact of the (quantitatively) incorrect tail of the potential curve. We have seen in Figure 2 that CAMB3LYP has asymptotically a correct behavior specially for the higher states; however, it is quantitatively questionable and for some states seems to be inaccurate. In such cases, the spectroscopic constants are calculated relative to the shallowest point after the minimum, and not to the asymptotic point. This yields approximately the same Re and ωe, but the obtained value De will be definitely shallower than, or approximately equal to a value De relating to the “correct” asymptotic point. We note that all values marked with an “*” in Tables 7-8 are obtained this way. For the 1Σg+ state, we see from Table 7 that CAMB3LYP has a good agreement with [19], likewise B3LYP with [17], whereas BLYP shows an agreement with DK-CASPT2 value of [17], but to conclude we see that the result(s) of 1Σg+ are widely distributed; furthermore, the lack of any experimental value makes the situation more difficult.


## 4.2.2. Higher Excited States

To deal with more higher excited states is difficult because of the above-mentioned reasons. Available approximations do not describe the long-range behavior correctly and/or fail to offer the correct asymptotic limit or predict it accurately [97]. We will discuss the higher molecular states given in Table 6 corresponding to the last four asymptotes (3–6) in Table 2. The result is given for the functionals BPW91 and BP86 (pure), B3LYP (hybrid), CAMB3LYP (range-separated), PBE0 (long-range corrected), and its gradient corrected one GRAC-PBE0. GRAC is an interpolation scheme, it is an asymptotic correction and supposed to be able to deal with higher excited states [37, 38]. The pw shows that the best result is obtained for CAMB3LYP and a comparable result is obtained for PBE0. Indeed strictly only CAMB3LYP was able to deal with higher excited states, it shows (at least qualitatively) the correct asymptotic as can be clearly seen in Figure 2. Other functionals do not show a correct asymptotic behavior as expected [37], including the ones for which no data shown in Table 8. B3LYP is given in Figure 2 as an example, it mixes the asymptotic for higher states with lower states. Our conclusion based on analyzing the data of all functionals and comparing them with each other. It is clear that lacking to the correct long-range behavior is primarily the origin of the problem, CAMB3LYP is able to cure this although not accurately, the question is why other corrections like GRAC does not have the expected improvement? At one side important is the nonlocal part of exact exchange which improves the situation considerably when the two-electron interaction is separated in short- and long-range part such as in CAMB3LYP, and we notice that there is no long-range correlation present in CAMB3LYP because HF offers only (nonlocal) exchange. Another point is the wrong long-range behavior of the response function [72, 77] caused by the incorrect long-range behavior of the density functional approximation is more crucial than it might be believed. This is supported by the fact that the spatial nonlocality of fxc is strongly frequency-dependent [98], in [98] Tokatly and Pankratov argued that not only any static approximation but also any LDA-based dynamic approximation (including any gradient corrections) for fxc cannot provide consistent result. To my best knowledge, there is no calculated or experimental result reported for any of the higher states given in Table 8, this makes the situation more difficult to analyze and be clarified. In Table 8 surprisingly we see that PBE0 gives a better result for higher excited states than its asymptotic corrected one GRAC-PBE0 and better than B3LYP, BP86, or BPW91. Furthermore, it gives for all states a comparable result to CAMB3LYP for Re and ωe. This supports our view and stress the importance of the long-range correction. It is a clear evident that PBE0 has a correct shape in inner part of the potential curve, and only its asymptotic part (tail of the potential curve) is incorrect, unfortunately the applied correction of GRAC is not good. As seen in Table 8 our next four states, 3Σu+, 3Σg+, and 1Σu+, 1Σg+ corresponding to Atom((4s2)1S + Atom((4s5s)3S), and Atom((4s2)1S + Atom((4s5s)1S), have more or less a similar result for all functionals, only GRAC-PBE0 shows unexplainable result, since it is supposed to show asymptotically a better behavior. We think that the CAMB3LYP result is the most correct one although it might be not satisfactory accurate. It is worthwhile to mention that states with avoiding crossing get a second shallow minimum after the avoiding crossing at large internuclear distances, this is not reported and only the first minimum is presented. Next, we look to the states 3πu, 3πg, 3Σu+, 3Σg+ corresponding to the Atom((4s2)1S + Atom((4s5p)3P). Here, we see that the result is distributed, BPW91, BP86, and B3LYP show similar results, whereas GRAC-PBE0 differs considerably from all approximations given in Table 8. PBE0 result is close to CAMB3LYP when looking to Re and ωe, but its De values are different clearly due to its incorrect asymptotic limit. The last states treated in this work 1πu, 1πg, 1Σu+, 1Σg+ are corresponding to the Atom((4s2)1S + Atom((4s5p)1P). The results of 1πu are puzzling and presumably only the values of CAMB3LYP are reasonable, whereas for 1πg all functional apart from GRAC-PBE0 give comparable values for ωe and Re, which could be a hint that these values are reasonable. 1Σu+, and 1Σg+ follow the general trend that PBE0 result is close to CAMB3LYP. BPW91, BP86, and B3LYP show a similar result, GRAC-PBE0 shows unexplainable result.The general conclusion of this section is that CAMB3LYP gives the best result due to its better treatment of the long-range part of the two-electron interaction and its asymptotically better behavior (tail of the potential curve) apparently due to including a suitable amount of exact exchange, PBE0 gives a comparable result, the main problem here is the tail of the potential curve. BPW91, BP86, and B3LYP are less satisfactory but still show acceptable result, whereas (most likely) the result of GRAC-PBE0 is not useful.


## 5. Conclusion

In the present work, we have studied the ground as well the 20 lowest exited states of the zinc dimer in the framework of DFT and TDDFT using well-known and newly developed functional approximations. We performed the calculations with  Dirac-Package  using relativistic 4-component DCH and SFH. First, we showed that SFH is capable to achieve the same accuracy as 4-components DCH and can describe quantitatively the main relevant contributions of the relativistic effects. In analyzing the results obtained from different functional approximations, comparing them with each other, with literature and experimental values as far as available, we drew some conclusions. The results show that the linear response in the adiabatic approximation with the known DFT approximations give good performance for the 8 lowest excited states ofZn2. For higher excited states, we found, somehow as expected, that most of DFT approximations used in the pw did not show a correct long-range behavior and the correct asymptotic limit to perform a fair accuracy for these states, where we have to stress that the lack of experimental or other theoretical results makes a judgment difficult. Nevertheless, we can say that the best result is obtained with the range-separated CAMB3LYP functional, which was the only one able (at least qualitatively) to show the correct asymptotic behavior. This can be led back to the separation of the two-electron interaction in a suitable manner, short- and long-range part, where the former is handled by the DFT and the later by HF. Showing that including a suitable (parametric) amount of the exact exchange improves the result considerably. Moreover, the (long-range corrected) PBE0 was able to give a comparable result to CAMB3LYP for the higher states although it fails to give the correct asymptotes. The comparison between CAMB3LY and other functionals allows us to conclude that for higher states the lack of a correct long-range and a suitable amount of exact exchange is responsible for incorrect result rather than the linear response approximation and the adiabatic limit. In addition, it causes a wrong long-range behavior of the response function a crucial point for the long-range behavior in TDDFT. In future works, we will be concerned with the heavier members of the group 12, Cd2, and Hg2, where relativistic effects are expected to be more important than in zinc dimer. Furthermore, the superheavy dimer Cn2 is under consideration, where the bonding character of its ground and excited states of academic interest due to the large relativistic effects and its influence on the atomic levels and hence on the molecular ground and excited states of the dimer.



---

*Source: 361947-2012-08-28.xml*