Theoretical background

Hilbert-space partitioning

In order to obtain information of the atomic contributions in a given chemical system (for instance atomic populations and electron sharing indices) it is crucial to define an atom in a molecule (AIM), which can either be real-space partition (allocating each point of the 3D space fully or partially to a specific atom) or Hilbert-space partition (separating the atomic basis functions belonging to a certain atom). The ESI-3D code [1] developed by Dr. Eduard Matito mainly used Bader’s Quantum Theory of Atoms in Molecules (QTAIM, real-space scheme) [2] as the AIM for the calculations. However, in this program we propose the use of Hilbert-space schemes (Mulliken [3], Löwdin [4], Meta-Löwdin [5], NAO [6], and IAO [7]) available in the PySCF [8, 9, 10, 11] framework as the partition of the system. QTAIM relies on numerical integrations, so the unavoidable errors associated to them make some of these aromaticity descriptors unviable in large systems. This newer approach, however, does not require numerical integration, but rather relies on the separation of the molecule by using their atom-centered functions, leading to an exact partition of the system. The following formulas will be expressed for single-determinant, closed-shell wavefunctions. The most fundamental magnitude is the Atomic Overlap Matrix ( AOM, \(\mathbf{S}^{\text{A}}\) ) in the Molecular Orbitals (MO, \(\mathbf{\phi}\) ) basis, with elements

\[S_{ij}^\text{A}=\int_{\text{A}}\phi_i^*(\textbf{r})\phi_j(\textbf{r})\text{d}\textbf{r}.\]

The average number of electrons in a given atom (\(N_\text{A}\)) can be expressed in terms of an Atomic Orbitals (AO, \(\chi\)) basis as

\[N_{\text{A}} = \sum_{\nu\in\text{A}}^\text{M} \sum_\mu^\text{M} P_{\nu\mu}S_{\mu\nu}^\text{AO} = \sum_{\nu\in\text{A}}^\text{M} (PS^\text{AO})_{\nu\nu}\]

where we can introduce the elements of the overlap matrix in an AO basis, \(S_{\mu\nu}^\text{AO}=\int\chi_\mu^{*}(\textbf{r}){\chi_\nu}(\textbf{r})d\textbf{r}\). The elements of the P-matrix, \(P_{\nu\mu} = 2 \sum_i ^{M} c_{\nu i} c_{i\mu}^+\), showcase the orbital occupancies. In the simplest case of a single-determinant wavefunction, it s a unit matrix over all occupied for \(\text{M}>nocc\) and 0s in the rest. For multi-determinant wavefunctions, the diagonalization of the P-matrix in the MO representation gives the natural orbital occupancies (\(n_i\)) and the transformation matrix to the new basis (\(\Gamma\)), which are not constrained to occupied orbitals anymore, by performing a unitary transformation:

\[\phi^{NO} = \Gamma^{+}C^{+}\chi^{\text{AO}}C\;\Gamma = (\Gamma^{'})^+\chi^{\text{AO}}\;\Gamma^{'}\]

In this sense, \(\Gamma\) is the diagonal representation of the MO basis, \(C\) the transformation matrix from MOs into AOs, and \(\Gamma^{'} = \Gamma\, C\) the transformation matrix from AOs into NOs.

In the simplest case of a single-determinant wavefunction, Mulliken’s approach lets us obtain information from a specific atom by only taking into account its atomic basis functions. Moreover, the Delocalization Index (DI, \(\delta\)), also referred to as Bond Order (BO) [12], measures the average number of electrons shared between two atoms A and B, as

\[\delta(\text{A,B})=\sum^\text{M}_{\mu\in\text{A}}\sum^\text{M}_{\nu\in\text{B}}(PS^\text{AO})_{\nu\mu}(PS^\text{AO})_{\mu\nu}.\]

In order to mimic the expression of the AOM as that of QTAIM, one can introduce a new auxiliary matrix, \(\mathbf{\eta}^{\text{A}}\), which is a bock-truncated unit matrix with all elements being zero except \(\eta_{\mu\mu}^\text{A}=1\) for \(\mu\in\text{A}\). Hence, the general expression for Mulliken’s approach is the following:

\[\mathbf{S}^\text{A,Mull}=\mathbf{c}^{+}\mathbf{S}^{AO}\mathbf{\eta}^\text{A}\mathbf{c}.\]

The resulting matrix is non-symmetric due to the underlying AO basis being non-orthogonal. To overcome these issues, chemists have explored alternative Hilbert-space methods that rely on orthogonalized AO bases, mainly obtained through a unitary transformation of the original AO basis used in calculations. Löwdin first proposed the symmetric orthogonalization procedure by using \(T_{\mu\nu}=S_{\mu\nu}^{-1/2}\). Following his steps, several different approaches have been reported to find more robust schemes of basis set orthogonalization, being the ones applied in this article the meta-Löwdin and the Natural Atomic Orbitals (NAO). Alternatively, Knizia proposed an ingenious scheme to express in an exact number the occupied MOs of a calculation in an orthogonal basis of reduced rank, the so-called Intrinsic Atomic Orbitals (IAO) approach. In all cases, the mapping from real-space to Hilbert-space can be performed as follows:

\[\mathbf{S}^\text{A,X}=\mathbf{c}^{+}({\mathbf{T}}^{-1})^{+}\mathbf{\eta}^\text{A}\mathbf{T}^{-1}\mathbf{c}.\]

Electron-Sharing Indices

The Electron Sharing Indices (ESI) present in this program rely on the atomic overlap matrices. The following aromaticity indicators will be expressed in the set of connected atoms in ring connectivity \(\mathscr{A}=\{\text{A}_1, \text{A}_2, \cdot\cdot\cdot, \text{A}_n\}\).

Para-delocalization index (PDI)

Fulton [13] reported that the delocalization indices in a given aromatic 6-membered ring in the para position were larger than that in the meta position. From that idea, Poater and coworkers proposed to average the DIs in the para position in a 6-membered ring, so the Para-Delocalization Index (PDI) [14] reads as:

\[\text{PDI}(\mathscr{A}) = \frac{\delta_{\text{A}_1\text{A}_4}+\delta_{\text{A}_2\text{A}_5}+\delta_{\text{A}_3\text{A}_6}}{3},\]

A larger PDI value indicates a more aromatic character. The index can only be calculated for rings of \(n=6\), so it will not be computed for rings of different sizes.

Iring

Giambiagi and coworkers proposed to express an index in terms of the generalized bond order along the ring, the \(\textbf{I}_{\textbf{ring}}\) [15]. That is, to account for the delocalization along the ring, following the specified connectivity:

\[\text{I}_{\text{ring}}(\mathscr{A})= 2^{n} \sum_{i_1,i_2\ldots i_n} S_{i_1i_2}^{\text{A}_{1}} S_{i_2i_3}^{\text{A}_{2}} \cdot \cdot \cdot S_{i_ni_1}^{\text{A}_{n}}\]

This index relies on the multicenter character of a molecule. A larger I\(_{\text{ring}}\) value indicates larger aromaticity along the ring.

Multicenter index (MCI)

As an aim to improve the I\(_{\text{ring}}\), Bultinck and coworkers proposed the Multicenter Index (MCI) [16] by not only taking into account the Kekulé structure of the system but rather all the \(n!\) possible ring connectivities generated by permuting the position of all atoms in the ring. That way, the delocalization is measured throughout the system, rather than along the ring. Denoting the different permutations as \(\mathscr{P}(\mathscr{A})\):

\[\text{MCI}(\mathscr{A}) = \frac{1}{2n} \sum_{\mathscr{P}(\mathscr{A})} \text{I}_{\text{ring}}(\mathscr{A})\]

As well as the previous indices, a larger MCI value denotes a more aromatic character. Due to the exponential growth of the calculation, we do not suggest computing the MCI for rings larger than \(n=12\) for single-core processes and \(n=14\) for multi-core processes. See MCI TIMINGS for details and timings of the algorithms.

AV1245 and AVmin

When using real-space schemes as the atomic partition, their inherent numerical integration errors made the multicenter indices in large rings non-viable. Matito proposed an index that contained the multicenter character as those of I\(_{\text{ring}}\) and MCI, but without the size-extensivity problem. Therefore, he suggested to average all the 4c-MCI values along the ring that keep the positional relationship of 1,2,4,5, so designing the new index AV1245 [17] as follows:

\[\text{AV1245}(\mathscr{A}) = \frac{1000}{3} \sum_{i=1}^n\text{MCI}(\{\text{A}_i, \text{A}_{i+1}, \text{A}_{i+3}, \text{A}_{i+4}\})\]

where if \(i>n\) \(\text{A}_i\) should be replaced by \(\text{A}_{i-n}\). In addition, Matito defined the AVmin index as the minimum (absolute) value of all the 4-MR MCI indices that enter the AV1245 expression. A higher AV1245 and AVmin values indicates more aromaticity in the system, and the index can not be computed for rings smaller than 6 centers.

Fluctuation Index (FLU)

The Fluctuation Index (FLU) [18] measures the resemblance of a series of reference \(\delta\) to some typical aromatic molecules:

\[\text{FLU}(\mathscr{A}) = \frac{1}{n} \sum_{i=1}^{n} \left[\left(\frac{V(A_i)}{V(A_{i-1})} \right)^\alpha \frac{\delta(A_i, A_{i-1}) - \delta_{ref}(A_i, A_{i-1})}{\delta_{ref}(A_i, A_{i-1})} \right]^2\]

Where one can separate it into two parts: the polarizability of the bond and the comparison to some reference \(\delta\) ( for instance, the “CC”, “CN”, “BN”, “NN” and “CS” bonds). \(\alpha\) is a simple function to make sure the first term is always greater or equal to 1. The index is close to zero for aromatic molecules and greater than zero in non-aromatic or antiaromatic molecules, and should not be used to study reactivity as they measure the similarity with respect to some molecule.

Bond Order Alternation (BOA)

The Bond Order Alternation (BOA) reflects the alternation of the delocalization indices along a conjugated circuit and is built upon the BLA premise (see below in the Geometrical Indicators section):

\[\text{BOA}(\mathscr{A}) = \frac{1}{n_1} \sum_{i=1}^{n_1} \delta(A_{2i-1},A_{2i}) - \frac{1}{n_2} \sum_{i=1}^{n_2} \delta(A_{2i},A_{2i+1})\]

where \(n_1 = \lfloor (n+1)/2 \rfloor\) and \(n_2 = \lfloor n/2 \rfloor\), being \(\lfloor x \rfloor\) the floor function of \(x\) returning the largest integer less or equal than \(x\). As well as for the BLA index, for odd-centered closed circuits this index may fail, so instead the \(\text{BOA}_c\) index should be used as the comparison of \(\delta(A_i, A_{i+1}) - \delta(A_{i+1}, A_{i+2})\):

\[\text{BOA}_c(\mathscr{A}) = \frac{1}{N} \sum_{i=1}^{N} \left| \delta(A_{i},A_{i+1}) - \delta(A_{i+1},A_{i+2}) \right|\]

Geometrical Indicators

HOMA and HOMER

The Harmonic Oscillator Model of Aromaticity (HOMA) [19] was defined by Kruszewski and Krygowski and relies only on geometrical data.

\[\text{HOMA}(\mathscr{A}) = 1 - 257.7\frac{1}{n} \cdot \sum_i^n (R_{opt} - R_{A_i,A_{i+1}})^2 = 1 - 257.7\frac{1}{n} \cdot ((R_{opt} - \bar{R})^2 + \sum_i^n (R_{A_i,A_{i+1}} - \bar{R})^2) = 1 - (EN + GEO)\]

The formula depends on a series of tabulated optimal bond distances, \(R_{opt}\), as well as the normalization factor \(\alpha\) for each bond to make the index 1 for benzene and 0 and negative values for non-aromatic or antiaromatic molecules, which makes it a good option for most organic molecules but fails for newer systems. The HOMA index is separated into the EN and GEO subparts, which measure the deviation of the interatomic distance into some tabulated numbers and the variance of this interatomic distance, respectively, and are close to zero for aromatic molecules. The implemented version of this index is [19]. The HOMER aromaticity index is a reparametrization of the HOMA for the lowest lying triplet excited state, T1 [20]. Different parameters can be introduced using the homarefs and homerrefs attributes.

Bond-Length Alternation (BLA)

The Bond-Length Alternation (BLA) index measures the average of the bond lengths of consecutive bonds in the ring

\[\text{BLA}(\mathscr{A}) = \frac{1}{n_1} \sum_{i=1}^{n_1} r_{A_{2i-1},A_{2i}} - \frac{1}{n_2} \sum_{i=1}^{n_2} r_{A_{2i},A_{2i+1}}\]

where \(n_1 = \lfloor (n+1)/2 \rfloor\) and \(n_2 = \lfloor n/2 \rfloor\), being \(\lfloor x \rfloor\) the floor function of \(x\) returning the largest integer less or equal than \(x\). This index was designed for open chains, and thus does not provide reliable results for closed circuits with and odd number of centers, so in those cases this index should be dismissed. Instead, one can use its closed-circuits counterpart, \(\text{BLA}_c\):

\[\text{BLA}_c(\mathscr{A}) = \frac{1}{N} \sum_{i=1}^{N} \vert r_{A_{i},A_{i+1}} - r_{A_{i+1},A_{i+2}} \vert\]

This new definition can indeed be used for closed rings, but produces numbers that even if qualitatively agree with BLA, they do not match completely.

References

1

Eduard Matito. ESI-3D: Electron Sharing Indexes Program for 3D Molecular Space Partitioning. 2006.

2

Richard F. W. Bader. Atoms in molecules: a quantum theory. Number 22 in The International series of monographs on chemistry. Clarendon Press ; Oxford University Press, 1994. ISBN 9780198558651.

3

R. S. Mulliken. Electronic Population Analysis on LCAO–MO Molecular Wave Functions. I. J. Chem. Phys., 1955.

4

Per-Ölov Löwdin. On the Non-Orthogonality Problem Connected with the Use of Atomic Wave Functions in the Theory of Molecules and Crystals. J. Chem. Phys., 1950.

5

Qiming Sun and Garnet Kin-Lic Chan. Exact and Optimal Quantum Mechanics/Molecular Mechanics Boundaries. J. Chem. Theory Comput., 2014.

6

Alan E. Reed, Robert B. Weinstock, and Frank Weinhold. Natural population analysis. J. Chem. Phys, 1985.

7

Gerald Knizia. Intrinsic Atomic Orbitals: An Unbiased Bridge between Quantum Theory and Chemical Concepts. J. Chem. Theory Comput., 2013.

8

Qiming Sun, Timothy C Berkelbach, Nick S Blunt, George H Booth, Sheng Guo, Zhendong Li, Junzi Liu, James D McClain, Elvira R Sayfutyarova, Sandeep Sharma, and others. Pyscf: the python-based simulations of chemistry framework. WIREs, Comput. Mol. Sci., 8(1):e1340, 2018.

9

Qiming Sun, Xing Zhang, Samragni Banerjee, Peng Bao, Marc Barbry, Nick S. Blunt, Nikolay A. Bogdanov, George H. Booth, Jia Chen, Zhi-Hao Cui, Janus J. Eriksen, Yang Gao, Sheng Guo, Jan Hermann, Matthew R. Hermes, Kevin Koh, Peter Koval, Susi Lehtola, Zhendong Li, Junzi Liu, Narbe Mardirossian, James D. McClain, Mario Motta, Bastien Mussard, Hung Q. Pham, Wirawan Purwanto Artem Pulkin, Paul J. Robinson, Enrico Ronca, Elvira R. Sayfutyarova, Maximilian Scheurer, Henry F. Schurkus, James E. T. Smith, Chong Sun, Shi-Ning Sun, Shiv Upadhyay, Lucas K. Wagner, Xiao Wang, Alec White, James Daniel Whitfield, Mark J. Williamson, Sebastian Wouters, Jun Yang, Jason M. Yu, Tianyu Zhu, Timothy C. Berkelbach, Sandeep Sharma, Alexander Yu. Sokolov, and Garnet Kin-Lic Chan. Recent developments in the pyscf program package. J. Chem. Phys., 2020.

10

Q. Sun, T.C. Berkelbach, N.S. Blunt, G.H. Booth, S. Guo, Z. Li, J. Liu, J.D. McClain, E.R. Sayfutyarova, S. Sharma, S. Wouters, and G.K.-L. Chan. Pyscf: the python-based simulations of chemistry framework. Comput. Mol. Sci., 2018.

11

Qiming Sun. Libcint: an efficient general integral library for gaussian basis functions. J. Comput. Chem., 2015.

12

I. Mayer. Charge, bond order and valence in the AB initio SCF theory. Chem. Phys. Lett., 1983.

13

Robert L. Fulton. Sharing of electrons in molecules. J. Phys. Chem., 97:7516–7529, 1993.

14

Jordi Poater, Xavier Fradera, Miquel Duran, and Miquel Solà. The delocalization index as an electronic aromaticity criterion: application to a series of planar polycyclic aromatic hydrocarbons. Chem. Eur. J., 9:400–406, 2003.

15

Mario Giambiagi, Myriam Segre De Giambiagi, Cassia D. Dos Santos Silva, and Aloysio Paiva De Figueiredo. Multicenter bond indices as a measure of aromaticity. Phys. Chem. Chem. Phys., 2(15):3381–3392, 2000. doi:10.1039/b002009p.

16

Patrick Bultinck, Robert Ponec, and Sofie Van Damme. Multicenter bond indices as a new measure of aromaticity in polycyclic aromatic hydrocarbons. J. Phys. Org. Chem., 18(8):706–718, August 2005. doi:10.1002/poc.922.

17

Eduard Matito. An electronic aromaticity index for large rings. Phys. Chem. Chem. Phys., 18(17):11839–11846, 2016. doi:10.1039/C6CP00636A.

18

E. Matito, M. Duran, and M. Solà. The aromatic fluctuation index (flu): a new aromaticity index based on electron delocalization. J. Chem. Phys., 122:014109, 2005.

19(1,2)

J Kruszewski and T. M. Krygowski. Definition of aromaticity basing on the harmonic oscillator model. Tetrahedron Letters, 13:3839–3842, 1972.

20

Enrique M. Arpa and Bo Durbeej. HOMER: a reparameterization of the harmonic oscillator model of aromaticity (HOMA) for excited states. Phys. Chem. Chem. Phys., 25(25):16763–16771, 2023. doi:10.1039/D3CP00842H.