(sec:pop.detailed)= # 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. (sec:pop.output.detailed)= ## Controlling Output The array `Print` allows the control of output. The general way of assigning elements of `Print` is: ```orca %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. ```orca %output PrintLevel Nothing Mini Small Normal Maxi Large Huge Debug end ``` `Print` has presently the following elements that can be user controlled: :::{flat-table} :header-rows: 1 * - Flag - Action * - `P_InputFile` - Echo the input file * - `P_Cartesian` - Print the cartesian coordinates * - `P_Internal` - Print the internal coordinates * - {rspan}`1` `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 * - {rspan}`2` `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 * - {rspan}`5` `P_DFTD` - Print info on Grimme's dispersion correction * - print mini $=$ 0 * - print small $=$ 1 * - print normal $=$ 1 * - print maxi $=$ 2 * - print huge $=$ 2 * - {rspan}`5` `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` | | (sec:pop.mulliken.detailed)= ## Mulliken Population Analysis The Mulliken population analysis {cite}`mulliken1955pop` 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: $$\int{ \rho \left({ \vec{r} } \right)d\vec{r} } =N $$ (eqn:299) and from the density matrix $\mathrm{\mathbf{P} }$ and the basis functions $\phi$: $$\rho \left({ \vec{r} } \right)=\sum\limits_{\mu \nu } { P_{\mu \nu } \phi _{\mu } \left({ \vec{r} } \right)\phi_{\nu } \left({ \vec{r} } \right)} $$ (eqn:300) therefore: $$\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 } } } $$ (eqn:301) $$%\hspace{1cm} =\sum\limits_{\mu \nu } { P_{\mu \nu } S_{\mu \nu } } $$ (eqn:302) After assigning each basis function to a given center this can be rewritten: $$\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} } } } } $$ (eqn:303) $$\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$ 0 K, fractional occupation numbers or FOD analysis, see {ref}`sec:scfconv.temp.detailed`), 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 {ref}`sec:plots.surface.fod.detailed`). The following example shows the corresponding printout for the first carbon atom of *p*-benzyne based on a FOD analysis with default settings (see {ref}`sec:scfconv.temp.fod.detailed`). ``` ------------------------------------------ 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 ({ref}`sec:pop.loewdin.detailed`) which is turned on by default using the total SCF density of the calculation, also in the case of finite electronic temperature. (sec:pop.loewdin.detailed)= ## Löwdin Population Analysis The Löwdin analysis {cite}`szabo1989modern` 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: $$\mathrm{\mathbf{P} }^{L}=\mathrm{\mathbf{S} }^{1/2}\mathrm{\mathbf{PS} }^{1/2} $$ (eqn:308) Then the atomic populations are: $$N_{A} =\sum\limits_\mu{ ^{A}P_{\mu \mu }^{L} } $$ (eqn:309) The bond order is defined from the Wiberg index {cite}`wiberg1968tetrahedron` that was first used in the context of semiempirical methods (that are formulated in the Löwdin basis right from the start): $$B_{AB} =\sum\limits_\mu{ ^{A}\sum\limits_\nu{ ^{B}\left({ P_{\mu \nu }^{L} } \right)^{2} } } $$ (eqn:310) 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. ```orca %output Print[ P_Loewdin ] 1 # default = on end ``` The flags to regulate the details are almost identical: ```orca %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. ```orca %method LOEWDIN_BONDORDERTHRESH 0.05 end ``` (sec:pop.mayer.detailed)= ## Mayer Population Analysis Mayers bonding analysis {cite}`mayer1983chem,mayer1984int,mayer1985theor,mayer1987modelling` 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: $$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 } } } $$ (eqn:311) 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: $$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 } } } $$ (eqn:312) 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: $$X_{A} =V_{A} -\sum\limits_{B\ne A} { B_{AB} } $$ (eqn:313) and finally the free valence (a measure of the ability to form further bonds) is given by: $$F_{A} =V_{A} -X_{A} $$ (eqn:314) The Mayer population analysis is turned on by: ```orca %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. ```orca %method MAYER_BONDORDERTHRESH 0.1 end ``` (sec:pop.natural.detailed)= ## 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. ## 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: $${\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)$$ 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: ```orca ! 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 1 2 ---------------------------------- 1 : 3.7568 2 : -2.2500 3.7568 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: ```orca ! 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: ``` 1 2 ---------------------------------- 1 : 3.7500 2 : -3.7500 3.7500 * Seff(A) -------------------------- 1 : n.a. 1.5000 2 : n.a. 1.5000 * = for a singlet state all 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: ```orca ! 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. ``` 1 2 ---------------------------------- 1 : 7.7009 2 : -5.3721 7.7012 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. (sec:pop.UNO.detailed)= ## UNO Orbital Printing The analysis of UNO's can be controlled similarly. The flags together with their default values are shown below: ```orca %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 ``` (sec:pop.hirshfeld.detailed)= ## Hirshfeld Charges The partitioning method by Hirshfeld is one of the most used approaches in the so-called atoms in molecules (AIM) methods.{cite}`hirshfeld1977` In this case, the AIM density of atom A, $\rho_A(\vec{r})$ is written as: $$ \rho_A(\vec{r}) = \rho(\vec{r})w_A(\vec{r}) $$ (eq_aim_charge) 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: $$ w_A(\vec{r}) = \frac{\rho_A^0(\vec{r})}{\rho^0(\vec{r})} $$ (eq_w_hirshfeld) 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. {eq}`eq_aim_charge` is known as _stockholder_. From eqs. {eq}`eq_aim_charge` and {eq}`eq_w_hirshfeld` one can calculate the Hirshfeld charges as: $$ Q_A^{\mathrm{Hirsh.}} = Z_A -\int\rho_A(\vec{r})d\vec{r} $$ (eq_q_hirshfeld) 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 ```orca ! Hirshfeld ``` in the ORCA input file, or alternatively via the `%output` block: ```orca %output Print[ P_Hirshfeld ] 1 # default = off end ``` For instance, if we request the Hirshfeld charges for a water molecule: ```orca !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: ```orca ------------------ 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 ``` (sec:pop.mbis.detailed)= ## MBIS Charges The Minimal Basis Iterative Stockholder (MBIS) method is a variant of the Hirshfeld method.{cite}`mbis2016` 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})$: $$ \rho_A^0(\vec{r}) = \sum_{i=1}^{m_A}\rho_{Ai}^0(\vec{r}) $$ (eq_sum_mbis) with $\rho_{Ai}^0(\vec{r})$ equal to: $$ \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) $$ (eq_rai_mbis) Here, $m_A$ is the number of shells of atom $A$. The populations $N_{Ai}$, and the widths $\sigma_{Ai}$ can be written as: $$ N_{Ai} = \int\rho(\vec{r})\frac{\rho_{Ai}^0(\vec{r})}{\rho^0(\vec{r})}d\vec{r} $$ (eq_nai_mbis) $$ \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} $$ (eq_sai_mbis) 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. {eq}`eq_sum_mbis` and {eq}`eq_rai_mbis`, (2) the new set of $N_{Ai}$ and $\sigma_{Ai}$ are obtained via eqs. {eq}`eq_nai_mbis` and {eq}`eq_sai_mbis`, (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: $$ Q_A^{\mathrm{MBIS.}} = Z_A -\int\rho_A(\vec{r})d\vec{r} $$ (eq_q_mbis) The calculation of the MBIS charges in ORCA is requested by writing ```orca ! MBIS ``` in the ORCA input file, or alternatively via the `%output` block: ```orca %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: ```orca !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: ```orca ------------------ 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: ```orca %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: ```orca %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: ```orca %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 ```