22.1 Molecular Properties, Wavefunction Analysis, and Localized Orbitals

Molecular properties (electrostatic moments, relativistic corrections, population analyses for densities and MOs, construction of localized MOs, NTOs, etc.) can be calculated with proper. This program is menu-driven and reads the input that determines which properties are evaluated from standard input (i.e. the terminal or an input file if started as proper < inputfile). The control file and files referenced therein are only used to determine the molecular structure, basis sets, and molecular orbitals and to read results computed before with other programs.

proper is a post-processing tool, mainly intended for interactive use which reads (almost all) its input from the terminal so that it is (usually) not necessary to modify the control file.

Several functionalities are also integrated in the programs that generate MOs or densities and can be invoked directly from the modules dscf, ridft, rimp2, mpgrad, ricc2 and egrad, if corresponding keywords are set in the control file. If one wants to skip the MO- or density generating step for dscf, ridft, rimp2, mpgrad, ricc2 it is possible to directly jump to the routine that carries out the analyses by starting the program with "<program> -proper". (For ricc2 it is, however, recommended to use instead ricctools -proper.) Currently, the respective keywords have to be inserted by hand (not with define) in the control file.

Here we briefly present the functionalities. A detailed description of the keywords that can be used in combination with the -proper flag is found Section 25.2.29.

22.1.1 Selection of densities

The proper program tries on start to read all densities that have been pre-calculated with any of the other programs of the TURBOMOLE package, prints a list with the densities that have been found and selects one of them for the calculation of properties. The default choice can be changed in the menu that is entered with the option dens. Here one can also list or edit additional attributes of the densities or build linear combinations of the available densities to form differences and superposition of densities.

The feature can thus be used to evaluate difference densities between the ground and excited electronic states or differences between densities calculated with different electronic structure methods.

The selected density can then be used in the subsequent menues to evaluate a variety of properties as expectation values or on a grid of points to generate interface files for visualization.

22.1.2 Electrostatic moments

Use the mtps option in the eval menu in proper, or add the $moments keyword control file when you start a program with the -proper flag to evaluate electrostatic moments. Up to quadrupole moments are calculated by default, on request also octupole moments are available. By default unnormalized traced Cartesian moments are calculated which are defined as \[Q^{(n)}_{\alpha\beta\ldots\nu} = \int d\mathbf{r} \rho(\mathbf{r}) r_\alpha r_\beta \ldots r_\nu ~,\] where \(\rho(\mathbf{r})\) is the charge density as position \(\mathbf{r}\). With the option Buckingham one can request in addition the computation of Cartesian traceless (Buckingham) multipole moments defined as: \[M^{(n)}_{\alpha\beta\ldots\nu} = \frac{(-1)^n}{n!} \int d\mathbf{r} \rho(\mathbf{r}) r^{2n} \frac{d^n}{dr_\alpha dr_\beta \ldots dr_\nu} \frac{1}{r}\]

22.1.3 Relativistic corrections

The option relcor in the eval menu of proper or the keyword $mvd (when starting a program with the -proper flag) initiate the calculation of relativistic corrections. With the -proper flag they are calculated for the SCF or DFT total density in case of dscf and ridft, for the SCF+MP2 density in case of rimp2 and mpgrad and for that of the calculated excited state in case of egrad. Quantities calculated are the expectation values \(\langle p^2\rangle\) , \(\langle p^4\rangle\) and the Darwin term \(\langle \sum_A 1/Z_A*\rho(R_A) \rangle\). Note, that at least the Darwin term requires an accurate description of the cusp in the wave function, thus the use of basis sets with uncontracted steep basis functions is recommended. Moreover note, that the results for these quantities are not really reasonable if ECPs are used (a respective warning is written to the output).

22.1.4 Population analyses

For population analyses enter the pop menu of proper. The available options and parameters that can be specified are the same as those for the $pop keyword (vide infra). If an electronic structure program is started with the -proper flag the population analyses is requested with the keyword $pop.

mulliken or $pop without any extension start a Mulliken population analysis (MPA). For -proper the analysis is carried out for all densities present in the respective program, e.g. total (and spin) densities leading to Mulliken charges (and unpaired electrons) per atom in RHF(UHF)-type calculations in dscf or ridft, SCF+MP2 densities in rimp2 or mpgrad, excited state densities in egrad. Suboptions (see Section 25.2.29) also allow for the calculation of Mulliken contributions of selectable atoms to selectable MOs including provision of data for graphical output (simulated density of states).

With $pop nbo a Natural Population Analysis (NPA) [383] is done. Currently only the resulting charges are calculated.

With $pop paboon a population analyses based on occupation numbers [384] is performed yielding "shared electron numbers (SENs)" and multicenter contributions. For this method always the total density is used, i.e. the sum of alpha and beta densities in case of UHF, the SCF+MP2-density in case of MP2 and the GHF total density for (two-component-)GHF. Note that the results of such an analysis may depend on the choice of the number of modified atomic orbitals ("MAOs"). By default, the number of MAOs is chosen such that they are reasonable in most cases (see Section 25.2.29). Nevertheless it is recommended to read carefully the information concerning MAOs given in the output before looking at the results for atomic charges and shared electron numbers. For different ways of selecting MAOs see Section 25.2.29.

With $pop wiberg Wiberg bond indices (WBI) [385] are calculated.

22.1.5 Generation of localized MOs

The option lmos in the mos menu of proper and the keyword $localize with the -proper flag trigger the calculation of localized molecular orbitals. Per default a Boys localization including all occupied MOs is carried out (i.e. the squared distance of charge centers of different LMOs is maximized). Alternative localization methods are Pipek-Mezey, Intrinsic Bond Orbitals (IBOs), and minimizations of (powers) of the second or fourth orbital moment.

pm

Pipek-Mezey localization, maximizes the sum of the squared orbital charges on the nuclei: \[L_{\texttt{PM}} = \sum_i \sum_A^{atoms} \big(q_i^A\big)^2\] It has the lowest operation count of all the implemented localization procedures, but gives only well-localized orbitals for basis sets without diffuse functions. In difference to Foster-Boys (see below) it conserves the \(\pi\sigma\) separation.

boys

Foster-Boys localization, minimizes the average orbital spread: \[L_{\texttt{FB}} = \sum_i \langle i | (\mathbf{r} - \mathbf{r}_i)^2 | i \rangle\] where \(\mathbf{r}_i = \langle i |\mathbf{r}|i\rangle\) is the center of the LMO. It is almost as fast as Pipek-Mezey, but stable with the basis set. Its disadvantage relative to Pipek-Mezey and IBO localization is that it breaks the \(\pi-\sigma\) separation for double bonds.

sm

second moment localization, minimizes a power of the sum of orbital spreads: \[L_{\texttt{SM},n} = \sum_i \langle i | (\mathbf{r} - \mathbf{r}_i)^2 | i \rangle^n\] The exponent is defined with the additional option exp=n where n must be an integer \(>0\). For \(n=1\) SM localization is identical to Foster-Boys. Values \(n>1\) cause a larger penalty for LMOs with larger spreads and thereby reduce the variance of the LMO spreads at the price of a small increase of the average spread and a small increase in the computational costs.

fm

fourth moment localization[386], mimimizes a power of the sum of the fourth central moments: \[L_{\texttt{FM},n} = \sum_i \langle i | (\mathbf{r} - \mathbf{r}_i)^4 | i \rangle^n\] The exponent is again defined with the additional option exp=n where n must be an integer \(>0\). Fourth moment localization produces in contrast to PM, IBO, FB, and SM localization LMOs with smaller tails Values \(n>1\) cause a larger penalty for LMOs with a larger Kurtosis and thereby reduce the variance of the LMO Kurtosis at the price of a small increase of the average Kurtosis. The computational costs for FM localization are 3–5 times larger than for SM localization.

ibo

intrinsic bond orbital localization[387], maximizes the sum of the fourth power of the orbital charges on the nuclei: \[L_{\texttt{IBO}} = \sum_i \sum_A^{atoms} \big(q_i^A\big)^4\] The IBO localization is similar to Pipek-Mezey localization, but the localization is carried out in in a minimal basis of intrinsic atomic orbitals (IAOs, see below) which make the procedure stable with respect to basis basis extension. The exponent of 4 reduces the problem of (nearly) degenerate minima of the PM functional. The cost of the IBO localization is for large systems and/or basis sets dominated by the costs for the construction of the IAOs.

As output one gets localized MOs (written to files lmos or lalp/lbet in UHF cases). In addition information about the dominant contributions of canonical MOs to the LMOs and a population analysis and, if requested, also about the centers, the spread and the Kurtosis of the LMOs is written to standard output.

22.1.6 Intrinsic Bond Orbitals Analysis

The option ibos in the mos menu of proper initiates the computation of intrinsic bond orbitals[387] (IBOs). Instead of a Mulliken PA in the AO basis the contributions a population analysis in the basis of the intrinsic atomic orbitals[387] (IAOs) and an IAO atomic charge analysis is done. The start guess for the IBO optimization is slightly different then the one used with the option lmos, which in the case of degenerate solutions (core orbitals, symmetry-degenerate LMOs located at the same centers) can lead to slightly different IBOs. As proposed in Ref. [387] the IAOs are obtained by projection of the occupied orbitals onto pre-computed atomic orbitals from Hartree-Fock calculations on the isolated atoms in the correlation consistent triple-\(\zeta\) basis (cc-pVTZ and cc-pVTZ-PP for pseudo-potential basis sets). These are available for most main group and transition metal atoms, for atoms beyond Kr with the standard pseudo-potential cores. If additional definitions are needed, they can be added to the basis set library.

Available reference orbitals for IAOs:

For transition (earth) alkali and transition metal atoms the inclusion of the highest \(s\) orbital (\(2s\) for Li–Be, \(3s\) for Na–Mg, and \(4s\) for K–Zn) can cause problems for systems that contain the respect cations. In the cations these \(s\) orbitals are (essentially) unoccupied and lead the unphysical delocalized IAOs. To avoid this problem it is possible to restrict the reference basis for the construction of the IAOs by adding to the control file data group $iaoopts:

$iaoopts
   subset=3s2p1d <list of atom indes>

22.1.7 Natural transition orbitals

For excited states calculated at the CIS (or CCS) level the transition density between the ground and an excited state \[\begin{align} E_{ia} = \langle \Psi_{ex}|a^\dagger_i a_a |\Psi_{ex} \rangle \end{align}\](22.1) can be brought to a diagonal form through a singular value decomposition (SVD) of the excitation amplitudes \(E_{ia}\): \[\begin{align} [\mathbf{O}^\dagger \mathbf{E} \mathbf{V}]_{ij} = \delta_{ij} \sqrt{\lambda}_i \end{align}\](22.2) The columns of the matrices \(\mathbf{O}\) and \(\mathbf{V}\) belonging to a certain singular value \(\lambda_i\) can be interpreted as pairs of occupied and virtual natural transition orbitals[388, 389] and the singular values \(\lambda_i\) are the weights with which this occupied-virtual pair contributes to the excitation. Usually electronic excitations are dominated by one or at least just a few NTO transitions and often the NTOs provide an easier understanding of transition than the excitation amplitudes \(E_{ia}\) in the canonical molecular orbital basis.

From excitation amplitudes computed with the ricc2 program NTOs and their weights (the singular values) can be calculated with the ntos option in the mos menu of proper or with ricctools. E.g. using the right eigenvectors for the second singlet excited state in irrep 1 with:

Both programs store the results for the occupied and virtual NTOs in files named, respectively, ntos_occ and ntos_vir. The option nto in the grid menu of the proper program can used to evaluate NTOs for visualization on a grid of points.

Note that the NTO analysis ignores for the correlated methods (CIS(D), ADC(2), CC2, CCSD, etc.) the double excitation contributions and correlation contributions to the ground state. This is no problem for single excitation dominated transition out of a “good” single reference ground state, in particular if only a qualitative picture is wanted, but one has to be aware of these omissions when using NTOs for states with large double excitation contributions or when they are used for quantitative comparisons.

Difference densities based on natural transition orbitals

If the excitation vectors have been obtained starting from a GHF reference, the NTOs are complex and contain contributions from both spin function. Moreover, the transitions are usually dominated by two NTOs at least. Thus, the interpretation of 2c-NTOs may become difficult. To get a simple picture of the transition at hand still, approximate difference densities can be computed according to \[\begin{equation} \rho {({\bf r})}_n = \Re\Big(\sum_{ab}^{N_\text{vir}} {\phi_a ({\bf r})}^\dagger \phi_b ({\bf r}) \sum_i^{N_\text{occ}} {C_i^a}^* C_i^b - \sum_{ij}^{N_\text{occ}} {\phi_i ({\bf r})}^\dagger \phi_j ({\bf r}) \sum_a^{N_\text{vir}} {C_i^a}^* C_j^a \Big) \, . \end{equation}\](22.3) The first term corresponds to the increase of the occupation of the virtual NTOs, while the second term corresponds to the decrease of the occupation of the occupied NTOs.
This approximate difference density is available for excitation vectors obtained with the following methods: CCS/CIS, CIS(D\(\infty\)), ADC(2) and CC2. Symmetry other than C\(_1\) is currently not supported. Note that the approximate difference densities are based on the same approximations as the NTOs, namely ignoring correlation and double excitation contributions.
From excitation amplitudes computed with the ricc2 program the approximate difference densities are computed with ricctools. E.g. using the right eigenvectors for the second singlet excited state in irrep 1:

This resulting density file can be visualized using the analysis mode of the ricc2 program as described in Section 10.3.3, e.g. by adding the following lines to the control file

 
$anadens
 calc my_approx_diffden from
 1d0 cc2-1a-002-approxdiffden.cao
$pointval

and running

  ricc2 -fanal

22.1.8 Corresponding Spin Orbitals

The analysis of spin-unrestricted open-shell calculations (UHF or UKS) are often hampered by the fact that the spatial parts of canonical \(\alpha\)- and \(\beta\)-spin orbitals can differ a lot. This makes it difficult to identify singly occupied molecular orbitals (SOMOs) and to distinguish almost doubly occupied orbitals from SOMOs and strongly spin-polarized “magnetic orbitals” that might be present in multireference situations as e.g. during the dissociation of covalent bonds.

Corresponding spin orbitals (CSOs) are defined through unitary transformations of the \(\alpha\)- and the \(\beta\)-spin orbitals that maximize the similarity (or overlap) of the spatials parts of \(\alpha\)- and \(\beta\)-spin orbitals with same indezes and thereby rotate strongly-polarized “magnetic” parts and SOMOs into a small set of orbitals. This is achieved by a singular value decomposition (SVD) of the overlap matrix between the occupied \(\alpha\)- and \(\beta\)-spin orbitals: \[\begin{align} \mathbf{U}^T \mathbf{S}^{\alpha\beta,oo} \mathbf{V} = \mathbf{s} \end{align}\](22.4) where \(s_{ij} = \delta_{ij} s_i\) and \(s_i\) are the singular values. The occupied CSOs \(\tilde{\phi}^\sigma_k\) are then obtained by transforming the canonical MOs \(\phi^\sigma_k\) with the unitary matrices \(\mathbf{U}\) and \(\mathbf{V}\): \[\begin{align} \tilde{\phi}^\alpha_k & = \sum_j \phi^\alpha_j U_{jk} \\ \tilde{\phi}^\beta_k & = \sum_j \phi^\beta_j V_{jk} \end{align}\](22.5–22.6) The CSOs are sorted according to decreasing singular values which range from 1.0 to 0.0 and can be classified as follows:

• For CSO pairs with singular values close to 1.0 the spatials parts of the \(\alpha\)- and \(\beta\)-spin orbitals are almost the same and they correspond thus to closed-shell MOs in spin-restricted calculations.

• For the majority spin there will be \(|n_\alpha - n_\beta|\) CSOs with singular values of 0.0: These are the SOMOs.

• For CSO pairs with singular values significantly lower than 1.0 (but non-zero) the spatial parts of \(\alpha\)- and \(\beta\)-spin orbitals are significantly different: these are strongly spin-polarized “magnetic orbitals”. Their presence indicates multireference situations.

In the implementation in TURBOMOLE similar transformations are applied to the virtual \(\alpha\)- and \(\beta\)-spin orbitals, but in this case sorted according to increasing singular values for \(\mathbf{S}^{\alpha\beta,vv}\). Thus, in the output CSO sets the unoccupied counterparts of the SOMOs (in the CSOs for the minority spin) have the same spatials parts as the SOMOs. They are then followed by strongly spin-polarized unoccupied CSOs and the non-polarized CSO pairs that correspond to the virtuals orbitals of spin-restricted calculations have the highest indezes.

By default the coefficients of the CSOs are written to file in the cartesian AO basis (as it is done for NTOs). With the option lsymao one can request that the output is in the symmetry-adapted spherical AO basis (as for canonical MOs).

For single-reference cases without strongly spin-polarized “magnetic orbitals” one might want to isolate only the SOMOs, but keep the remaining MOs that correspond to the closed-shell and virtuals MOs of spin-restricted calculations close to canonical MOs. The can be done obtained wit the option match. If switched on, the occupied orbitals for the minority spin are kept in the canonical basis. For the majortiy spin, the CSOs are first calculated as described above and then the \(|n_\alpha - n_\beta\) CSOs with the largest singular values are transformed with the inverse of the transformation matrix for other spin to make as similar as possible to the canonical MOs of the minority spin. Again, similar transformations are applied to the virtuals MOs.

22.1.9 Orbitals for weakly interacting fragments

The frag option in the mos menu of properallows to extract from a supermolecular HF or DFT calculation on weakly interaction fragments MOs and occupation numbers for the individual fragments. This is can be in particular usefull for supermolecular UHF or UKS calculations to obtain information on the spin states of the fragments. The option requires as prerequisite that the control file contains a valid input for the $frag data group assigning all atoms to fragments. Furthermore, the interaction between the fragments must be weak so that after localization each occupied LMO can be assigned to one of the fragments. All further input options for frag are passed as input to the calculation of LMOs and have the same meaning as for the LMO option.

For each fragment the program will generate files containing the orbitals and coordinates and simple control files containing the orbital occupations, basis set information and references to the coordinate and MO data groups that can be used to run calculations for the individual fragments.

22.1.10 Fit of charges due to the electrostatic potential:

$esp_fit fits point charges at the positions of nuclei to electrostatic potential arising from electric charge distribution (for UHF cases also for spin density, also possible in combination with $soghf). For this purpose the ("real") electrostatic potential is calculated at spherical shells of grid points around the atoms. By default, Bragg-Slater radii, \(r_{BS}\), are taken as shell radii.

A parametrization very close to that suggested by Kollman (a multiple-shell model with shells of radii ranging from 1.4*\(r_{vdW}\) to 2.0*\(r_{vdW}\), \(r_{vdW}\) is the van-der-Waals radius; U.C. Singh, P.A. Kollman, J. Comput. Chem. 5(2), 129-145 (1984)) is used if the keyword is extended:

$esp_fit kolman