10 Molecular Properties and Analysis

10.16 Population of Effectively Unpaired Electrons

In a stretched hydrogen molecule the two electrons that are paired at equilibrium forming a bond become un-paired and localized on the individual H atoms. In singlet diradicals or doublet triradicals such a weak paring exists even at equilibrium. At a single-determinant SCF level of the theory the valence electrons of a singlet system like H2 remain perfectly paired, and one needs to include non-dynamical correlation to decouple the bond electron pair, giving rise to a population of effectively-unpaired (“odd”, radicalized) electrons.946, 84, 906 When the static correlation is strong, these electrons remain mostly unpaired and can be described as being localized on individual atoms.

These phenomena can be properly described within wave-function formalism. Within DFT, these effects can be described by broken-symmetry approach or by using SF-TDDFT (see Section 7.3.1). Below we describe how to derive this sort of information from pure DFT description of such low-spin open-shell systems without relying on spin-contaminated solutions.

The first-order reduced density matrix (1-RDM) corresponding to a single-determinant wave function (e.g., SCF or Kohn-Sham DFT) is idempotent:

ρσ(𝐫1)=γσSCF(1;2)γσSCF(2;1)𝑑𝐫2γσSCF(1;2)=ioccψiσKS(1)ψiσKS(2), (10.129)

where ρσ(1) is the electron density of spin σ at position 𝐫1, and γσSCF is the spin-resolved 1-RDM of a single Slater determinant. The cross product γσSCF(1;2)γσSCF(2;1) reflects the Hartree-Fock exchange (or Kohn-Sham exact-exchange) governed by the HF exchange hole:

γσSCF(1;2)γσSCF(2;1)=ρα(1)hXσσ(1,2)hXσσ(1,2)𝑑𝐫2=1. (10.130)

When 1-RDM includes electron correlation, it becomes non-idempotent:

Dσ(1)ρσ(1)-γσ(1;2)γσ(2;1)𝑑𝐫20. (10.131)

The function Dσ(1) measures the deviation from idempotency of the correlated 1-RDM and yields the density of effectively-unpaired (odd) electrons of spin σ at point 𝐫1.946, 791 The formation of effectively-unpaired electrons in singlet systems is therefore exclusively a correlation based phenomenon. Summing Dσ(1) over the spin components gives the total density of odd electrons, and integrating the latter over space gives the mean total number of odd electrons N¯u:

Du(1)=2σDσ(1)d𝐫1,N¯u=Du(1)𝑑𝐫1. (10.132)

The appearance of a factor of 2 in Eq. (10.132) above is required for reasons discussed in Ref. 791. In Kohn-Sham DFT, the SCF 1-RDM is always idempotent which impedes the analysis of odd electron formation at that level of the theory. Ref. 788 has proposed a remedy to this situation. It was noted that the correlated 1-RDM cross product entering Eq. (10.131) reflects an effective exchange, also known as cumulant exchange.84 The KS exact-exchange hole is itself artificially too delocalized. However, the total exchange-correlation interaction in a finite system with strong left-right (i.e., static) correlation is normally fairly localized, largely confined within a region of roughly atomic size.58 The effective exchange described with the correlated 1-RDM cross product should be fairly localized as well. With this in mind, the following form of the correlated 1-RDM cross product was proposed:788

γσ(1;2)γσ(2;1)=ρσ(1)h¯Xσσeff(1,2). (10.133)

The function h¯Xσσeff(1;2) is a model DFT exchange hole of Becke-Roussel (BR) form used in Becke’s B05 method.47 The latter describes left-right static correlation effects in terms of certain effective exchange-correlation hole.47 The extra delocalization of the HF exchange hole alone is compensated by certain physically motivated real-space corrections to it:47

h¯XCαα(1,2)=h¯Xααeff(1,2)+fc(1)h¯Xββeff(1,2). (10.134)

The BR exchange hole h¯Xσσeff is used in B05 as an auxiliary function, such that the potential from the relaxed BR hole equals that of the exact-exchange hole. This results in relaxed normalization of the auxiliary BR hole less than or equal to unity:

h¯Xσσeff(1;2)𝑑𝐫2=NXσeff(1)1. (10.135)

The expression of the relaxed normalization NXσeff(𝐫) is quite complicated, but it is possible to represent it in closed analytic form.790, 789 The smaller the relaxed normalization NXαeff(1), the more delocalized the corresponding exact-exchange hole.47 The α-α exchange hole is further deepened by a fraction of the β-β exchange hole, fc(1)h¯Xββeff(1,2), which gives rise to left-right static correlation. The local correlation factor fc in Eq.(10.134) governs this deepening and hence the strength of the static correlation at each point:47

fc(𝐫)=min(fα(𝐫),fβ(𝐫),1) (10.136a)
0fc(𝐫)1 (10.136b)
fα(𝐫)=1-NXαeff(𝐫)NXβeff(𝐫). (10.136c)

Using Eqs. (10.136), (10.132), and (10.133), the density of odd electrons becomes:

Dα(1)=ρα(1)(1-NXαeff(1))=ρα(1)fc(1)NXβeff(1). (10.137)

The final formulas for the spin-summed odd electron density and the total mean number of odd electrons read:

Du(1)=4andopfc(1)[ρα(1)NXβeff(1)+ρβ(1)NXαeff(1)]N¯u=Du(𝐫𝟏)𝑑𝐫𝟏. (10.138)

Here acnd-opp=0.526 is the SCF-optimized linear coefficient of the opposite-spin static correlation energy term of the B05 functional.47, 789

It is informative to decompose the total mean number of odd electrons into atomic contributions. Partitioning in real space the mean total number of odd electrons N¯u as a sum of atomic contributions, we obtain the atomic population of odd electrons (FAr) as:

FAr=ΩADu(𝐫𝟏)𝑑𝐫1. (10.139)

Here ΩA is a subregion assigned to atom A in the system. To define these atomic regions in a simple way, we use the partitioning of the grid space into atomic subgroups within Becke’s grid-integration scheme.52 Since the present method does not require symmetry breaking, singlet states are calculated in restricted Kohn-Sham (RKS) manner even at strongly stretched bonds. This way one avoids the destructive effects that the spin contamination has on FAr and on the Kohn-Sham orbitals. The calculation of FAr can be done fully self-consistently only with the RI-B05 and RI-mB05 functionals. In these cases no special keywords are needed, just the corresponding EXCHANGE rem line for these functionals. Atomic population of odd electron can be estimated also with any other functional in two steps: first obtaining a converged SCF calculation with the chosen functional, then performing one single post-SCF iteration with RI-B05 or RI-mB05 functionals reading the guess from a preceding calculation, as shown on the input example below:

Example 10.46  To calculate the odd-electron atomic population and the correlated bond order in stretched H2, with B3LYP/RI-mB05, and with fully SCF RI-mB05

$comment
   Stretched H2: example of B3LYP calculation of the atomic population of odd electrons
   with post-SCF RI-BM05 extra iteration.
$end

$molecule
   0 1
   H   0.  0.   0.0
   H   0.  0.   1.5000
$end

$rem
   SCF_GUESS         CORE
   METHOD            B3LYP
   BASIS             G3LARGE
   PURECART          222
   THRESH            14
   MAX_SCF_CYCLES    80
   PRINT_INPUT       TRUE
   SCF_FINAL_PRINT   1
   INCDFT            FALSE
   XC_GRID           000128000302
   SYM_IGNORE        TRUE
   SYMMETRY          FALSE
   SCF_CONVERGENCE   9
$end

@@@

$comment
   Now one RI-B05 extra-iteration after B3LYP to generate the odd-electron atomic
   population and the correlated bond order.
$end

$molecule
   read
$end

$rem
   SCF_GUESS         READ
   EXCHANGE          BM05
   PURECART          22222
   BASIS             G3LARGE
   AUX_BASIS         riB05-cc-pvtz
   THRESH            14
   PRINT_INPUT       TRUE
   INCDFT            FALSE
   XC_GRID           000128000302
   SYM_IGNORE        TRUE
   SYMMETRY          FALSE
   MAX_SCF_CYCLES    0
   SCF_CONVERGENCE   9
   DFT_CUTOFFS       0
$end

@@@

$comment
   Finally, a fully SCF run RI-B05 using the previous output as a guess.
   The following input lines are obligatory here:
   PURECART     22222
   AUX_BASIS    riB05-cc-pvtz
   DFT_CUTOFFS  0
$end

$molecule
   read
$end

$rem
   SCF_GUESS        READ
   EXCHANGE         BM05
   PURECART         22222
   BASIS            G3LARGE
   AUX_BASIS        riB05-cc-pvtz
   THRESH           14
   PRINT_INPUT      TRUE
   INCDFT           FALSE
   IPRINT           3
   XC_GRID          000128000302
   SYM_IGNORE       TRUE
   SCF_FINAL_PRINT  1
   SYMMETRY         FALSE
   MAX_SCF_CYCLES   80
   SCF_CONVERGENCE  8
   DFT_CUTOFFS      0
$end

Once the atomic population of odd electrons is obtained, a calculation of the corresponding correlated bond order of Mayer’s type follows in the code, using certain exact relationships between FAr, FBr, and the correlated bond order of Mayer type BAB. Both new properties are printed at the end of the output, right after the multipoles section. It is useful to compare the correlated bond order with Mayer’s SCF bond order. To print the latter, use SCF_FINAL_PRINT = 1.