7.53. Population Analyses and Control of Output

At present ORCA knows three different ways of analyzing the computed SCF wavefunction that will be described below. All of these methods can produce a tremendous amount of output. However, this output can be precisely controlled by the user to his or her individual needs.

In general there is one compound key called PrintLevel which is there to choose reasonable amounts of output. All that PrintLevel does is to set certain flags in the array Print which holds the details about what to print and what not.

7.53.1. Controlling Output

The array Print allows the control of output. The general way of assigning elements of Print is:

%output
  PrintLevel     Normal
  Print[ Flag ]  0  # turn print off
                 1  # turn print on
                 n  # some flags are more sophisticated
end

The compound key PrintLevel can be used to select certain default settings for the print array. Specifying Print after PrintLevel can be used to modify these defaults.

%output
  PrintLevel     Nothing
                 Mini
                 Small
                 Normal
                 Maxi
                 Large
                 Huge
                 Debug
end

Print has presently the following elements that can be user controlled:

Flag

Action

P_InputFile

Echo the input file

P_Cartesian

Print the cartesian coordinates

P_Internal

Print the internal coordinates

P_Basis

\(=1:\) Print the basis set information

\(=2:\) Also print the primitives in input format

P_OneElec

Print of the one electron matrix

P_Overlap

Print the overlap matrix

P_KinEn

Print the kinetic energy matrix

P_S12

Print the S\(^{-1/2}\) matrix

P_GuessOrb

Print the initial guess orbitals

P_OrbEn

Print Orbital Energies

P_MOs

Print the MO coefficients on convergence

P_Density

Print the converged electron density

P_SpinDensity

Print the converged spin density

P_EHTDetails

Print initial guess extended Hückel details

P_SCFInfo

Print the SCF input flags

P_SCFMemInfo

Print the estimated SCF memory requirements

P_SCFIterInfo

\(=1 :\) print short iteration information

\(=2 :\) print longer iteration information

\(=3 :\) in a direct SCF also print integral progress

P_Fockian

Print Fockian matrix

P_DIISMat

Print DIIS matrix

P_DIISError

Print DIIS error

P_Iter_P

Print Density

P_Iter_C

Print MO coefficients

P_Iter_F

Print Fock matrix

P_Mayer

Print Mayer population analysis. Default \(=\) on.

P_NatPop

Print Natural population analysis. Default \(=\) off.

P_Hirshfeld

Print Hirshfeld population analysis. Default \(=\) off.

P_MBIS

Print MBIS population analysis. Default \(=\) off.

P_Mulliken

Print Mulliken population analysis. Default \(=\) on

P_AtCharges_M

Print Mulliken atomic charges

P_OrbCharges_M

Print Mulliken orbital charges

P_FragCharges_M

Print Mulliken fragment charges

P_FragBondOrder_M

Print Mulliken fragment bond orders

P_BondOrder_M

Print Mulliken bond orders

P_ReducedOrbPop_M

Print Mulliken reduced orb. charges

P_FragPopMO_M

Print Mulliken fragment population for each MO

P_FragOvlMO_M

Print Mulliken overlap populations per fragment pair

P_AtPopMO_M

Print Mulliken atomic charges in each MO

P_OrbPopMO_M

Print Mulliken orbital population for each MO

P_ReducedOrbPopMO_M

Print Mulliken reduced orbital population for each MO

P_Loewdin

Print Loewdin population analysis. Default \(=\) on.

P_AtCharges_L

Print Loewdin atomic charges

P_OrbCharges_L

Print Loewdin orbital charges

P_FragCharges_L

Print Loewdin fragment charges

P_FragBondOrder_L

Print Loewdin fragment bond orders

P_BondOrder_L

Print Loewdin bond orders

P_ReducedOrbPop_L

Print Loewdin reduced orb. charges

P_FragPopMO_L

Print Loewdin fragment population for each MO

P_FragOvlMO_L

Print Loewdin overlap populations per fragment pair

P_AtPopMO_L

Print Loewdin atomic charges in each MO

P_OrbPopMO_L

Print Loewdin orbital population for each MO

P_ReducedOrbPopMO_L

Print Loewdin reduced orbital population for each MO

P_NPA

Natural population analysis

P_NBO

Natural bond orbital analysis

P_Fragments

Print fragment information

P_GUESSPOP

Print initial guess populations

P_UNO_FragPopMO_M

Print Mulliken fragment population per UNO

P_UNO_OrbPopMO_M

Print Mulliken orbital pop. per UNO

P_UNO_AtPopMO_M

Print Mulliken atomic charges per UNO

P_UNO_ReducedOrbPopMO_M

Print Mulliken reduced orbital pop. per UNO

P_UNO_FragPopMO_L

Print Loewdin fragment population per UNO

P_UNO_OrbPopMO_L

Print Loewdin orbital pop. per UNO

P_UNO_AtPopMO_L

Print Loewdin atomic charges per UNO

P_UNO_ReducedOrbPopMO_L

Print Loewdin reduced orbital pop. per UNO

P_UNO_OccNum

Print occupation numbers per UNO

P_AtomExpVal

Print atomic expectation values

P_AtomBasis

Print atomic basis

P_AtomDensFit

Print electron density fit

P_Symmetry

Symmetry basic information

P_Sym_Salc

Symmetry process printing

P_SCFSTABANA

Information on progress, convergence, and results of the SCF stability analysis

P_DFTD

Print info on Grimme’s dispersion correction

print mini \(=\) 0

print small \(=\) 1

print normal \(=\) 1

print maxi \(=\) 2

print huge \(=\) 2

P_DFTD_GRAD

Print gradient info on Grimme’s dispersion correction

print mini \(=\) 0

print small \(=\) 0

print normal \(=\) 0

print maxi \(=\) 1

print huge \(=\) 2

P_G1EL2EL

Print one- and two-electron contributions to the g-tensor

The various choices for PrintLevel have the following defaults:

PrintLevel

Print settings

Mini

P_OrbEn

= 1

P_Cartesian

= 1

P_InputFile

= 1

P_SCFIterInfo

= 1

Small

all the previous plus

P_SCFInfo

= 1

P_Mayer

= 1

P_MULLIKEN

= 1

P_AtCharges_M

= 1

P_ReducedOrbPop_M

= 1

P_Loewdin

= 1

P_AtCharges_L

= 1

P_ReducedOrbPop_L

= 1

P_Fragments

= 1

P_FragCharges_M

= 1

P_FragBondOrder_M

= 1

P_FragCharges_L

= 1

P_FragBondOrder_L

= 1

Normal

all the previous plus

P_Internal

= 1

P_BondOrder_L

= 1

P_BondOrder_M

= 1

P_FragPopMO_L

= 1

P_ReducedOrbPopMO_L

= 1

P_SCFIterInfo

= 2

Maxi

all the previous plus

P_GuessOrb

= 1

P_MOs

= 1

P_Density

= 1

P_SpinDensity

= 1

P_Basis

= 1

P_FragOVLMO_M

= 1

P_OrbPopMO_M

= 1

P_OrbCharges_M

= 1

Huge

All the previous plus

P_OneElec

= 1

P_Overlap

= 1

P_S12

= 1

P_AtPopMO_M

= 1

P_OrbPopMO_M

= 1

P_AtPopMO_L

= 1

P_EHTDetails

= 1

Debug

print everything

7.53.2. Mulliken Population Analysis

The Mulliken population analysis [598] is, despite all its known considerable weaknesses, the standard in most quantum chemical programs. It partitions the total density using the assignment of basis functions to given atoms in the molecules and the basis function overlap. If the total charge density is written as \(\rho \left({ \vec{r} } \right)\) and the total number of electrons is \(N\) we have:

(7.409)\[\int{ \rho \left({ \vec{r} } \right)d\vec{r} } =N \]

and from the density matrix \(\mathrm{\mathbf{P} }\) and the basis functions \(\phi\):

(7.410)\[\rho \left({ \vec{r} } \right)=\sum\limits_{\mu \nu } { P_{\mu \nu } \phi _{\mu } \left({ \vec{r} } \right)\phi_{\nu } \left({ \vec{r} } \right)} \]

therefore:

(7.411)\[\int{ \rho \left({ \vec{{r} }} \right)d\vec{{r} }} =\sum\limits_{\mu \nu } {P_{\mu \nu } \underbrace{ \int{ \phi_{\mu } \left({ \vec{{r} }} \right)\phi _{\nu } \left({ \vec{{r} }} \right)d\vec{{r} }} }_{S_{\mu \nu } } } \]
(7.412)\[%\hspace{1cm} =\sum\limits_{\mu \nu } { P_{\mu \nu } S_{\mu \nu } } \]

After assigning each basis function to a given center this can be rewritten:

(7.413)\[\hspace{2.5cm} =\sum\limits_A { \sum\limits_B {\sum\limits_\mu{ ^{A}\sum\limits_\nu{ ^{B}P_{\mu \nu }^{AB} S_{\mu \nu }^{AB} } } } } \]
(7.414)\[\hspace{4cm} =\sum\limits_A { \sum\limits_\mu{^{A}\sum\limits_\nu{ ^{A}P_{\mu \nu }^{AA} S_{\mu \nu }^{AA} } } } +2\sum\limits_A { \sum\limits_{B<A} { \sum\limits_\mu{ ^{A}\sum\limits_\nu{^{B}P_{\mu \nu }^{AB} S_{\mu \nu }^{AB} } } } } \]

Mulliken proposed to divide the second term equally between each pair of atoms involved and define the number of electrons on center \(A\), \(N_{A}\), as:

(7.415)\[N_{A} =\sum\limits_\mu{ ^{A}\sum\limits_\nu{ ^{A}P_{\mu \nu }^{AA} S_{\mu \nu }^{AA} } } +\sum\limits_{B\ne A} { \sum\limits_\mu{ ^{A}\sum\limits_\nu{^{B}P_{\mu \nu }^{AB} S_{\mu \nu }^{AB} } } } \]

such that \(\sum\limits_A { N_{A} } =N\). The charge of an atom in the molecule is then:

(7.416)\[Q_{A} =Z_{A} -N_{A} \]

where \(Z_{A}\) is the core charge of atom \(A\). The cross terms between pairs of basis functions centered on different atoms is the overlap charge and is used in ORCA to define the Mulliken bond order:

(7.417)\[B_{AB} =2\sum\limits_\mu{ ^{A}\sum\limits_\nu{ ^{B}P_{\mu \nu }^{AB} S_{\mu \nu }^{AB} } } \]

The Mulliken population analysis is turned on by using:

%output
  Print[ P_Mulliken ] 1  # default = on
end

A number of additional options can be specified to control the details of the Mulliken population analysis. By default the Mulliken population analysis is turned on.

%output
  Print[ P_AtCharges_M  ] 1   # Print atomic charges
  Print[ P_OrbCharges_M ] 1   # Print orbital charges
  Print[ P_FragCharges_M] 1   # Print fragment charges
  Print[ P_BondOrder_M  ] 1   # Print bond orders
  Print[ P_FragBondOrder_M ] 1# Print fragment b.o.
  Print[ P_ReducedOrbPop_M ] 1# Print reduced orb. Charges
  Print[ P_AtPopMO_M    ] 1   # Print atomic charges in
                              # each MO
  Print[ P_OrbPopMO_M   ] 1   # Print orbital populaiton
                              # for each MO
  Print[ P_ReducedOrbPopMO_M] 1 # Print reduced orbital
                                # pop for each MO
  Print[ P_FragPopMO_M  ] 1     # Print the fragment
                                # population for for each MO
end

These options allow to get very detailed information about the computed wavefunctions and is much more convenient than to look at the MOs directly. A “reduced orbital population” is a population per angular momentum type. For example the sum of populations of each p\(_{z}\) orbital at a given atom is the reduced orbital population of the p\(_{z}\) function.

Note that for finite temperature HF or KS-DFT calculations (SmearTemp \(>\) 0 K, fractional occupation numbers or FOD analysis, see Finite Temperature HF/KS-DFT), only the Mulliken reduced orbital charges based on \(\rho^{FOD}\) will be printed. They can be used to get a first impression about the localization of hot electrons in the molecule without generating the corresponding FOD plot (see FOD plots). The following example shows the corresponding printout for the first carbon atom of p-benzyne based on a FOD analysis with default settings (see Fractional Occupation Number Weighted Electron Density (FOD)).

------------------------------------------
FOD BASED MULLIKEN REDUCED ORBITAL CHARGES
------------------------------------------
  0 C s       :     0.006371  s :     0.006371
      pz      :     0.016375  p :     0.030785
      px      :     0.009893
      py      :     0.004516
      dz2     :     0.004248  d :     0.010308
      dxz     :     0.000254
      dyz     :     0.004855
      dx2y2   :     0.000860
      dxy     :     0.000091
      f0      :     0.000006  f :     0.000378
      f+1     :     0.000014
      f-1     :     0.000309
      f+2     :     0.000002
      f-2     :     0.000006
      f+3     :     0.000010
      f-3     :     0.000032

If other population analysis printouts are wanted the user is referred to the Löwdin analysis (Löwdin Population Analysis) which is turned on by default using the total SCF density of the calculation, also in the case of finite electronic temperature.

7.53.3. Löwdin Population Analysis

The Löwdin analysis [839] is somewhat more straightforward than the Mulliken analysis. In the Löwdin method one changes to a basis where all overlap integrals vanish. This is accomplished via Löwdins symmetric orthogonalization matrix \(\mathrm{\mathbf{S} }^{-1/2}\). Using this transformation matrix the new basis functions are multicentered but are in a least square sense as similar as possible to the original, strictly localized, atomic basis functions. The similarity of the transformed functions and original functions is explored in the population analysis. The density matrix transforms as:

(7.418)\[\mathrm{\mathbf{P} }^{L}=\mathrm{\mathbf{S} }^{1/2}\mathrm{\mathbf{PS} }^{1/2} \]

Then the atomic populations are:

(7.419)\[N_{A} =\sum\limits_\mu{ ^{A}P_{\mu \mu }^{L} } \]

The bond order is defined from the Wiberg index [888] that was first used in the context of semiempirical methods (that are formulated in the Löwdin basis right from the start):

(7.420)\[B_{AB} =\sum\limits_\mu{ ^{A}\sum\limits_\nu{ ^{B}\left({ P_{\mu \nu }^{L} } \right)^{2} } } \]

The output for the Löwdin population analysis (that I personally prefer over the Mulliken analysis) is closely similar. By default the Löwdin population analysis is turned on and provides some more detail than the Mulliken analysis.

%output
  Print[ P_Loewdin ] 1  # default = on
end

The flags to regulate the details are almost identical:

%output
  Print[ P_AtCharges_L  ] 1     # Print atomic charges
  Print[ P_OrbCharges_L ] 1     # Print orbital charges
  Print[ P_FragCharges_L] 1     # Print fragment charges
  Print[ P_BondOrder_L  ] 1     # Print bond orders
  Print[ P_FragBondOrder_L ] 1  # Print fragment b.o.
  Print[ P_ReducedOrbPop_L ] 1  # Print reduced orb. Charges
  Print[ P_AtPopMO_L    ] 1     # Print atomic charges in
                                # each MO
  Print[ P_OrbPopMO_L   ] 1     # Print orbital population
                                # for each MO
  Print[ P_ReducedOrbPopMO_L] 1 # Print reduced orbital
                                # pop for each MO
  Print[ P_FragPopMO_L  ] 1     # Print the fragment
                                # population for each MO
end

In addition one can set, in the method block, the threshold for the printing of the bond order.

%method
  LOEWDIN_BONDORDERTHRESH 0.05
end

7.53.4. Mayer Population Analysis

Mayers bonding analysis [567, 568, 569, 570] is another creative attempt to define chemically useful indices. The Mayer atomic charge is identical to the Mulliken charge. The Mayer bond order is defined as:

(7.421)\[B_{AB} =\sum\limits_\mu{ ^{A}\sum\limits_\nu{ ^{B}\left({ \mathrm{\mathbf{PS} }} \right)_{\mu \nu } \left({ \mathrm{\mathbf{PS} }} \right)_{\nu \mu } +\left({ \mathrm{\mathbf{RS} }} \right)_{\mu \nu } \left({ \mathrm{\mathbf{RS} }} \right)_{\nu \mu } } } \]

Here \(\mathrm{\mathbf{P} }\) is the total electron density matrix and \(\mathrm{\mathbf{R} }\) is the spin-density matrix. These Mayer bond orders are very useful. Mayer’s total valence for atom \(A\) is defined as:

(7.422)\[V_{A} =2N_{A} -\sum\limits_\mu{ ^{A}\sum\limits_\nu{ ^{A}\left({ \mathrm{\mathbf{PS} }} \right)_{\mu \nu } \left({ \mathrm{\mathbf{PS} }} \right)_{\nu \mu } } } \]

In normal bonding situations and with normal basis sets \(V_{A}\) should be reasonably close to the valence of atom \(A\) in a chemical sense (i.e. close to four for a carbon atom). The bonded valence is given by:

(7.423)\[X_{A} =V_{A} -\sum\limits_{B\ne A} { B_{AB} } \]

and finally the free valence (a measure of the ability to form further bonds) is given by:

(7.424)\[F_{A} =V_{A} -X_{A} \]

The Mayer population analysis is turned on by:

%output
  Print[ P_Mayer ] 1  # default = on
end

The output is rather simple and short and can not be further controlled. By default the Mayer population analysis is turned on. In addition one can set, in the method block, the threshold for the printing of the bond order.

%method
  MAYER_BONDORDERTHRESH 0.1
end

7.53.5. Natural Population Analysis

A popular and useful method for population analysis is the natural population analysis due to Weinhold and co-workers. It is implemented in the NBO interface.

7.53.6. Local Spin Analysis

It is common practice in various areas of chemistry to think about the interaction of open-shell systems in terms of local spin states. For example, in dimeric or oligomeric transition metal clusters, the ‘exchange coupling’ between open shell ions that exist locally in high-spin states is a much studied phenomenon. Diradicals would be typical systems in organic chemistry that show this phenomenon. In quantum mechanics, however, the total spin is not a local property, but instead a property of the system as a whole. The total spin squared, \(S^2\), and its projection onto the z-axis, \(S_z\), commute with the non-relativistic Hamiltonian and hence, the eigenfunctions of the non-relativistic Hamiltonian can be classified according to good quantum numbers \(S\) and \(M\) according to:

\[{\mathbf{S} }_{}^2\left|{ \Psi _{}^{SM} } \right\rangle = S(S + 1)\left|{ \Psi _{}^{SM} } \right\rangle\]
\[S_z^{}\left|{ \Psi _{}^{SM} } \right\rangle = M\left|{ \Psi _{}^{SM} } \right\rangle\]

where \(\left|{ \Psi _{}^{SM} } \right\rangle\) is an exact eigenfunction of the non-relativistic Hamiltonian or an approximation to it that conserves the total spin as a good quantum number. The total spin itself is given by the sum over the individual electron spins as:

\[{\mathbf{S} } = \sum\limits_i { {\mathbf{s} }(i) }\]

And hence,

\[{\mathbf{S} }_{}^2 = \sum\nolimits_{i,j} { {\mathbf{s} }(i){\mathbf{s} }(j) }\]

is a two-electron property of the system. It is obviously not trivial to relate the chemically very meaningful concept of local spin to a rigorous quantum mechanical treatment. While there are various proposals of how to deal with this problem, we follow here a proposal of Clark and Davidson (Clark, A.E.; Davidson, E.R., J. Chem. Phys. 2001, 115, 7382-7392). The following equations are implemented in the SCF and CASSCF modules of Orca.

Clark and Davidson define fragment projection operators with the property:

\[P_A^{}P_B^{} = \delta _{AB}^{}P_A^{}\]

and:

\[\sum\limits_A { P_A^{} } = 1\]

Then using this identity:

\[{\mathbf{S} } = \sum\limits_i { \sum\limits_A { {\mathbf{s} }(i)P_{A(i) }^{} } }\]
\[{\mathbf{S} } = \sum\limits_A { \sum\limits_i { {\mathbf{s} }(i)P_A^{}(i) } }\]
\[\,\,\,\, = \sum\limits_A { {\mathbf{S} }_A^{}(i) }\]

they show that the local spin operators obey the standard relations for spin operators:

\[{\mathbf{S} }_A^{} = { \mathbf{S} }_A^\dagger\]
\[{\mathbf{S} }_A^{} \times{ \mathbf{S} }_A^{} = i\hbar{ \mathbf{S} }_A^{}\]

Hence

\[{\mathbf{S} }_{}^2 = \sum\limits_A { \sum\limits_B { {\mathbf{S} }_A^{}{\mathbf{S} }_B^{} } }\]

But then importantly:

\[{\mathbf{S} }_A^{}{\mathbf{S} }_B^{} = \sum\limits_i { \sum\limits_j { {\mathbf{s} }(i){\mathbf{s} }(j) } P_A^{}(i)P_B^{}(j) }\]
\[\,\,\,\,\,\,\,\,\,\,\,\, = { \textstyle{3 \over 4} }\delta _{AB}^{}\sum\limits_A { P_A^{}(i) } + \sum\limits_i { \sum\limits_{j > i} { {\mathbf{s} }(i){\mathbf{s} }(j)\{ P_A^{}(i)P_B^{}(j) + P_A^{}(j)P_B^{}(i)\} } }\]

With the first- and second-order density matrix:

\[\gamma ({\mathbf{x} },{\mathbf{x'} }) = N\int{ \Psi ({\mathbf{x} },{\mathbf{x} }_2^{},...,{\mathbf{x} }_N^{})\Psi _{}^*({\mathbf{x'} },{\mathbf{x} }_2^{},...,{\mathbf{x} }_N^{}) } d{\mathbf{x} }_2^{}...d{\mathbf{x} }_N^{}\]
\[\Gamma ({\mathbf{x} }_1^{},{\mathbf{x'} }_1^{};{\mathbf{x} }_2^{},{\mathbf{x'} }_2^{}) = \left({ \begin{matrix} N \cr 2 \cr \end{matrix} } \right)\int{ \Psi ({\mathbf{x} }_1^{},{\mathbf{x} }_2^{},...,{\mathbf{x} }_N^{})\Psi _{}^*({\mathbf{x'} }_1^{},{\mathbf{x'} }_2^{},...,{\mathbf{x} }_N^{}) } d{\mathbf{x} }_3^{}...d{\mathbf{x} }_N^{}\]

(with \({\textstyle{N \choose 2} }= { \textstyle{1 \over 2} }N(N - 1)\)). Then:

\[\left\langle { {\mathbf{S} }_A^{}{\mathbf{S} }_B^{} } \right\rangle = { \textstyle{3 \over 4} }\delta _{AB}^{}tr(\gamma P_A^{}) + 2tr(P_A^{}(1)P_B^{}(2){\mathbf{s} }(1){\mathbf{s} }(2)\Gamma (1,1;2,2))\]

In terms of the number of electrons on site ‘A’and the expectation value of \(S_z^A\)

\[\left\langle { S_z^A} \right\rangle = { \textstyle{1 \over 2} }tr(\gamma _{}^{\alpha - \beta }P_A^{})\]
\[\left\langle { N_{}^A} \right\rangle = tr(\gamma _{}^{\alpha + \beta }P_A^{})\]

in terms of molecular orbitals:

\[\left\langle { S_z^A} \right\rangle = { \textstyle{1 \over 2} }\sum\limits_{p,q} { \gamma _{pq}^{\alpha - \beta }\left\langle { p|P_A^{}|q} \right\rangle}\]
\[\left\langle { N_{}^A} \right\rangle = \sum\limits_{p,q} { \gamma _{pq}^{\alpha + \beta }\left\langle { p|P_A^{}|q} \right\rangle}\]

McWeeny and Kutzelnigg (McWeeny, R.; Kutzelnigg, W. Int. J. Quant. Chem. 1968, 11, 187-203) show that for the expectation value of s(1)s(2), the relevant irreducible part of the two-body density can be expressed in terms of the spinless density matrix of second order:

\[R_0^{(0) }(1,1';2,2') = - { \textstyle{1 \over 3} }\Gamma (1,1';2,2') - { \textstyle{2 \over 3} }\Gamma (2,1';1,2')\]
\[\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, = - { \textstyle{1 \over 3} }\sum\limits_{pqrs} { \Gamma _{rs}^{pq}p(1)q(1')r(2)s(2') } + 2\Gamma _{rs}^{pq}p(2)q(1')r(1)s(2')\]
\[\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, = - { \textstyle{1 \over 3} }\sum\limits_{pqrs} { (\Gamma _{rs}^{pq} + 2\Gamma _{ps}^{rq})p(1)q(1')r(2)s(2') }\]

with a normalization factor of \(3 \over 4\) after spin integration. Hence using this:

\[\left\langle { {\mathbf{S} }_A^{}{\mathbf{S} }_B^{} } \right\rangle = { \textstyle{3 \over 4} }\delta _{AB}^{}tr(\gamma P_A^{}) + { \textstyle{6 \over 4} }tr(P_A^{}(1)P_B^{}(2)R_0^{(0) }(1,1;2,2))\]

And then performing the integral:

\[\left\langle { {\mathbf{S} }_A^{}{\mathbf{S} }_B^{} } \right\rangle = { \textstyle{3 \over 4} }\delta _{AB}^{}tr(\gamma P_A^{}) - \underbrace{ {\textstyle{6 \over 4} }{\textstyle{1 \over 3} }}_{{\textstyle{1 \over 2} }}\sum\limits_{pqrs} { (\Gamma _{rs}^{pq} + 2\Gamma _{ps}^{rq})P_{pq}^AP_{rs}^B}\]

This is the final and perhaps most compact equation. The projection operator can be defined in very many different ways. The easiest is to Löwdin orthogonalize the basis set:

\[\left|{ \mu_L^A} \right\rangle = \sum\limits_{\nu _{}^A} { \left|{ \nu _{}^A} \right\rangle S_{\mu \nu }^{ - 1/2} }\]

where \(L\) denotes the Löwdin basis. This means that molecular orbitals are expressed in the orthogonal basis as:

\[{\mathbf{c} }_L^{} = { \mathbf{S} }_{}^{ + 1/2}{\mathbf{c} }\]

and the density as:

\[{\mathbf{P} }_L^{} = { \mathbf{S} }_{}^{ + 1/2}{\mathbf{PS} }_{}^{ + 1/2}\]

The fragment projector is defined as:

\[P_A^{} = \sum\limits_{\mu _L^{} \in A} { \left|{ \mu _L^{} } \right\rangle\left\langle { \mu _L^{} } \right|}\]

Clark and Davidson suggest a slightly more elaborate projector in which first, the intra-fragment overlap is eliminated. This happens with a matrix U that for two fragments takes form:

\[{\mathbf{U} } = \left( \begin{matrix}{ {\mathbf{S} }_A^{ - 1/2} } & { \mathbf{0} } \cr{ \mathbf{0} } & { {\mathbf{S} }_B^{ - 1/2} } \cr \end{matrix} \right)\]

where is the block of basis functions belonging to fragment A. Likewise:

\[ \begin{align}\begin{aligned}{\mathbf{U} }_{}^{ - 1} = \left( \begin{matrix}{ {\mathbf{S} }_A^{ + 1/2} } & { \mathbf{0} } \cr{ \mathbf{0} } & { {\mathbf{S} }_B^{ + 1/2} } \cr \\ \end{matrix} \right)\end{aligned}\end{align} \]

Then the ‘pre-overlap’is:

\[{\mathbf{\bar S} } = { \mathbf{U} }_{}^\dagger{ \mathbf{SU} }\]

This contains the unit matrix in the intra-fragment blocks and non-zero elements elsewhere. This overlap matrix is the finally orthogonalized to obtain the globally orthogonal Löwdin basis. We finally transform the MO coefficients by the following transformation:

\[{\mathbf{c} }_L^{} = { \mathbf{S} }_{}^{ + 1/2}{\mathbf{U} }_{}^{ - 1}{\mathbf{c} }\]

For the projectors, operating with the two MOs i and j gives:

\[\left\langle { i|P_A^{}|j} \right\rangle = \sum\limits_{\mu _L^{} \in A} { \sum\limits_{\kappa _L^B\tau _L^C} { \left\langle { \kappa _L^B|\mu _L^A} \right\rangle\left\langle { \mu _L^A|\tau _L^C} \right\rangle c_{\kappa i}^Lc_{\tau j}^L} }\]
\[\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, = \sum\limits_{\mu _L^{} \in A} { \sum\limits_{\kappa _L^B\tau _L^C} { \delta _{AB}^{}\delta _{AC}^{}\delta _{\kappa \mu }^{}\delta _{\tau \mu }^{}c_{\kappa i}^Lc_{\tau j}^L} }\]
\[\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, = \sum\limits_{\mu _L^{} \in A} { c_{\mu i}^Lc_{\mu j}^L}\]

Herrmann et al. (Herrmann, C.; Reiher, M.; Hess, B.A. J. Chem. Phys. 2005, 122, 34102) give the correct expression of the expectation values for a single spin-unrestricted determinant

\[\begin{aligned} & \left\langle { {\mathbf{S} }_A^{}{\mathbf{S} }_B^{} } \right\rangle = { \textstyle{3 \over 4} }\delta _{AB}^{}\left\{{ \sum\limits_i { P_{ii}^A} + \sum\limits_{\bar i} { P_{\bar i\bar i}^A} } \right\} \cr & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, + { \textstyle{1 \over 4} }\left\{{ \sum\limits_{ij} { P_{ii}^AP_{jj}^B} + \sum\limits_{\bar i\bar j} { P_{\bar i\bar i}^AP_{\bar j\bar j}^B} - \sum\limits_{ij} { P_{ij}^AP_{ij}^B} - \sum\limits_{\bar i\bar j} { P_{\bar i\bar j}^AP_{\bar i\bar j}^B} - \sum\limits_{\bar ij} { P_{\bar i\bar i}^AP_{jj}^B} - \sum\limits_{i\bar j} { P_{ii}^AP_{\bar j\bar j}^B} } \right\} \cr & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, - \sum\limits_{i\bar j} { P_{i\bar j}^AP_{i\bar j}^B} \cr\end{aligned}\]

Which is used in the Orca implementation.

The use of the Local spin-implementation is very easy. All that is required is to divide the molecule into fragments. The rest happens automatically. For example, let us consider two nitrogen atoms at the dissociation limit. While the total spin state is S=0, the tow nitrogen atoms local exist in high-spin states (S=3/2). Consider the following test job:

! HF def2-SVP UHF TightSCF PModel 

%scf brokensym 3,3 end 

* xyz 0 1
N(1) 0 0 0
N(2) 0 0 1094
*

and the output:

-------------------
LOCAL SPIN ANALYSIS (Loewdin* projector)
-------------------


(1) A.E. Clark; E.R. Davison J. Chem. Phys. (2001), 115(16), pp 7382-7392 
(2)  C. Herrmann, M. Reiher, B.A. Hess J. Chem. Phys. (2005) 122, art 034102-1

Number of fragments       = 2
Number of basis functions = 28
Number of atoms           = 2

   ... Fragment AO indices were mapped
   ... intra-fragment orthogonalization completed
   ... Global Loewdin orthogonalizer constructed
   ... Loewdin orthogonalized occupied orbitals constructed


 <SA*SB>        1            2     
----------------------------------
   1    :     3.7568  
   2    :    -2.2500       3.7568  


               <SzA>      Seff(A)
        --------------------------
   1    :     1.5000      1.5017
   2    :    -1.5000      1.5017

thus perfectly corresponding to the expectations. The same can be done at the CASSCF level:

! HF def2-SVP UHF TightSCF PModel

%casscf nel 6 norb 6 nroots 1 end 

* xyz 0 1
N(1) 0 0 0
N(2) 0 0 1094
*

With the result:

<SA*SB>        1            2     
----------------------------------
   1    :     3.7500  
   2    :    -3.7500       3.7500  


               <SzA>*      Seff(A)
        --------------------------
   1    :       n.a.      1.5000
   2    :       n.a.      1.5000

* = for a singlet state all <SzA> values are zero by definition

Thus, cleanly confirming the expectations.

As a less trivial example, consider a typical Fe(III) antiferromatically coupled transition metal dimer. An appropriate input may be:

! pbe def2-sv(p) tightscf kdiis soscf pmodel

%scf
  brokensym 5,5
end


* xyz -2 1
Fe(1)	-1.93818	 0.53739	-0.00010
Fe(2)	 1.06735	 0.47031	 0.00029
S(3)	-0.38935	 2.59862	-0.00983
S(3)	-0.48170	-1.59050	 0.01091
S(1)	 2.68798	 0.43924	 1.99710
S(1)	 2.68692	 0.42704	-1.99712
S(2)	-3.55594	 0.56407	-1.99889
S(2)	-3.55850	 0.58107	 1.99646
H(1)	 3.91984	 0.39462	 1.47608
H(1)	 3.91940	 0.39536	-1.47662
H(2)	-4.78410	 0.69179	-1.48280
H(2)	-4.78991	 0.49249	 1.47983
*

Where one of the bridging sulfurs was assigned to each site respectively.

<SA*SB>        1            2     
----------------------------------
   1    :     7.7009  
   2    :    -5.3721       7.7012  


               <SzA>      Seff(A)
        --------------------------
   1    :     1.7579      2.3197
   2    :    -1.7579      2.3198

Nice shows the expected results with the local site spins being close to their ideal value 2.5 which would hold for a high-spin Fe(III) ion.

7.53.7. UNO Orbital Printing

The analysis of UNO’s can be controlled similarly. The flags together with their default values are shown below:

%output
  Print[ P_UNO_OccNum     ] = 1; # Occupation numbers
  Print[ P_UNO_AtPopMO_M  ] = 0; # Mulliken atom pop.
                                 # per UNO
  Print[ P_UNO_OrbPopMO_M] = 0;  # Mulliken orbital pop.
                                 # per UNO
  Print[ P_UNO_ReducedOrbPopMO_M] = 0;
                                 # Mulliken reduced orbital
                                 # pop. per UNO
  Print[ P_UNO_AtPopMO_L  ] = 0; # Loewdin atom pop.
                                 # per UNO
  Print[ P_UNO_OrbPopMO_L] = 0;  # Loewdin orbital pop.
                                 # per UNO
  Print[ P_UNO_ReducedOrbPopMO_L] = 0;
                                 # Loewdin reduced orbital
                                 # pop. per UNO
end

7.53.8. Hirshfeld Charges

The partitioning method by Hirshfeld is one of the most used approaches in the so-called atoms in molecules (AIM) methods.[394] In this case, the AIM density of atom A, \(\rho_A(\vec{r})\) is written as:

(7.425)\[ \rho_A(\vec{r}) = \rho(\vec{r})w_A(\vec{r}) \]

Here, \(\rho(\vec{r})\) is the total charge density at position \(\vec{r}\), and \(w_A(\vec{r})\) a weighting function, that within the Hirshfeld method is equal to:

(7.426)\[ w_A(\vec{r}) = \frac{\rho_A^0(\vec{r})}{\rho^0(\vec{r})} \]

where \(\rho_A^0(\vec{r})\) is the pro-atomic density of atom \(A\) and \(\rho^0(\vec{r})=\sum_A\rho_A^0(\vec{r})\) the pro-molecular density. The ratio in eq. (7.425) is known as stockholder. From eqs. (7.425) and (7.426) one can calculate the Hirshfeld charges as:

(7.427)\[ Q_A^{\mathrm{Hirsh.}} = Z_A -\int\rho_A(\vec{r})d\vec{r} \]

In ORCA, the pro-atomic density within the Hirshfeld method is calculated via density fitting with a set of Gaussian s-functions per element.

The calculation of the Hirshfeld charges in ORCA is requested by writing

! Hirshfeld

in the ORCA input file, or alternatively via the %output block:

%output
  Print[ P_Hirshfeld ] 1  # default = off
end

For instance, if we request the Hirshfeld charges for a water molecule:

!HF cc-pvdz tightscf Hirshfeld

%maxcore 4000

* xyz 0 1
  O    0.00000006589375       0.00157184228646       0.00000000004493
  H    0.77316868532439      -0.58666889665624      -0.00000000000005
  H   -0.77316876182122      -0.58666895650640      -0.00000000000005
*

ORCA prints the following information at the end of the output file:

------------------
HIRSHFELD ANALYSIS
------------------

Total integrated alpha density =      4.999998580
Total integrated beta density  =      4.999998580

  ATOM     CHARGE      SPIN
   0 O   -0.333756    0.000000
   1 H    0.166879    0.000000
   2 H    0.166879    0.000000

  TOTAL   0.000003    0.000000

7.53.9. MBIS Charges

The Minimal Basis Iterative Stockholder (MBIS) method is a variant of the Hirshfeld method.[869] The idea behind this approach is that the pro-atomic density \(\rho_A^0(\vec{r})\) is expanded in a minimal set of atom-centered s-type Slater functions \(\rho_{Ai}^0(\vec{r})\):

(7.428)\[ \rho_A^0(\vec{r}) = \sum_{i=1}^{m_A}\rho_{Ai}^0(\vec{r}) \]

with \(\rho_{Ai}^0(\vec{r})\) equal to:

(7.429)\[ \rho_{Ai}^0(\vec{r}) = \frac{N_{Ai}}{\sigma_{Ai}^38\pi}\mathrm{exp}\left(-\frac{\left|\vec{r}-\vec{R_A}\right|}{\sigma_{Ai}}\right) \]

Here, \(m_A\) is the number of shells of atom \(A\). The populations \(N_{Ai}\), and the widths \(\sigma_{Ai}\) can be written as:

(7.430)\[ N_{Ai} = \int\rho(\vec{r})\frac{\rho_{Ai}^0(\vec{r})}{\rho^0(\vec{r})}d\vec{r} \]
(7.431)\[ \sigma_{Ai}=\frac{1}{3N_{Ai}}\int\rho(\vec{r})\frac{\rho_{Ai}^0(\vec{r})}{\rho^0(\vec{r})}\left|\vec{r}-\vec{R_A}\right| d\vec{r} \]

In order to compute the AIM densities \(\rho_A(\vec{r})\), the MBIS method uses an iterative algorithm where: (1) an initial guess is generated for the set of \(N_{Ai}\) and \(\sigma_{Ai}\) and the pro-atomic densities are calculated through eqs. (7.428) and (7.429), (2) the new set of \(N_{Ai}\) and \(\sigma_{Ai}\) are obtained via eqs. (7.430) and (7.431), (3) if convergence is reached for \(\rho_A(\vec{r})\), the iterative process stops, otherwise we go back to (1) but now one uses the last estimates for \(N_{Ai}\) and \(\sigma_{Ai}\).

Once, the MBIS iterative process stops, the MBIS charges are calculated as:

(7.432)\[ Q_A^{\mathrm{MBIS.}} = Z_A -\int\rho_A(\vec{r})d\vec{r} \]

The calculation of the MBIS charges in ORCA is requested by writing

! MBIS

in the ORCA input file, or alternatively via the %output block:

%output
  Print[ P_MBIS ] 1  # default = off
end

If we request the MBIS charges for a HF calculation at the cc-pvdz level of a chloroform molecule:

!HF cc-pvdz tightscf MBIS

* xyz 0 1
  C   -0.00000997794639     -0.00091664148112      0.45499807439812
  H    0.00000069467312      0.00031189002174      1.53703126401237
  Cl   0.00003188789531      1.69433732001280     -0.08420513240263
  Cl   1.46635420502892     -0.84684178730039     -0.08421103795485
  Cl  -1.46637680965097     -0.84689178125304     -0.08420916805301
*

ORCA prints the following information at the end of the output file:

------------------
MBIS ANALYSIS
------------------

Convergence threshold (charges)      ...        1.000e-06
Number of iterations                 ...        46

Total integrated alpha density       ...     29.000001385
Total integrated beta density        ...     29.000001385

  ATOM     CHARGE    POPULATION     SPIN
   0 C    0.208633    5.791367    0.000000
   1 H    0.169417    0.830583    0.000000
   2 Cl  -0.126877   17.126877    0.000000
   3 Cl  -0.125586   17.125586    0.000000
   4 Cl  -0.125590   17.125590    0.000000

  TOTAL  -0.000003   58.000003    0.000000

 MBIS VALENCE-SHELL DATA:
  ATOM   POPULATION   WIDTH(A.U.)
   0 C    4.122213    0.508675
   1 H    0.830583    0.358785
   2 Cl   8.532439    0.524031
   3 Cl   8.531380    0.523959
   4 Cl   8.531381    0.523959

The second block corresponds to the valence Slater function, which is caracterized by its population \(N_{A,v}\) and width \(\sigma_{A,v}\).

The convergence threshold for the MBIS charges is set to \(10^{-6}\). However, it can be changed via the tag MBIS_CHARGETHRESH in the %method block:

%method
  MBIS_CHARGETHRESH 0.0001
end

ORCA can also print the following MBIS-related quantities: (1) atomic dipole moments, (2) atomic quadrupole moments, (3) atomic octupole moments, and (4) third radial moment of the MBIS density \(\left(\langle r^3\rangle_A = \int \left|\vec{r}-\vec{R_A}\right|^3\rho_A(\vec{r})d\vec{r}\right)\). The printing of these properties is controlled by the tag MBIS_LARGEPRINT, to be specified in the %method block:

%method
  MBIS_LARGEPRINT TRUE  # default = FALSE
end

If this option is activated, an extra iteration is performed after reaching the convergence threshold for the charges.

The origin for the calculation of the atomic dipole, quadrupole and octupole moments is the center of each atom (default). However, the user can also define a global origin (independent of the atom) through the tag MBIS_ORIGIN_MULT in the %method block:

%method
  MBIS_ORIGIN_MULT  CenterOfCoords     # origin of coordinate system (0,0,0)
                    CenterOfMass       # center of mass
                    CenterOfNucCharge  # center of nuclear charge
                    CenterXYZ          # arbitrary position, set coordinates with MBIS_ORIMULT_XYZ
                    CenterOfEachAtom   # center of each atom (default)

  MBIS_ORIMULT_XYZ  x,y,z # set the coordinates, otherwise 0,0,0 (unit: Angstrom)
end