Introduction

Developing effective strategies to predict reactive site and chemical reactivity is one of the most important types of research directions in theoretical and computational chemistry1,2,3,4. Since the introduction of Mulliken atomic charge, numerous atomic charge calculation algorithms5 have been developed to forecast chemical reactivity at the atomic level. Noteworthy, nature population analysis (NPA) charge6 based on the nature bond orbital (NBO) method7 has gained widespread use. Subsequently, various chemical reactivity indexes derived from conceptual density functional theory methods8,9,10,11,12,13, including hardness, softness14,15, global and local indexes16,17,18, Fukui function19, and dual descriptor20, have been extensively developed and applied21,22 in chemical theory. These indices are also well-known as the conceptual density functional theory (CDFT)23,24,25 associated with analysis of the Euler-Lagrange equation derived from Hohenberg-Kohn theorem. They are invaluable for assessing nucleophilicity, electrophilicity, and other reactive properties of molecules and their constituent atoms. It should be mentioned that the condensed version of the density is without consideration of identifying specific reactivity directions. However, Fukui function and dual descriptor in an uncondensed way are the function of the position vector which give rise to plots in which the directional effects are prominent10, and a typical example is the Dunitz angle for the nucleophilic attack on a carbonyl group as retrieved from the Fukui function26. Moreover, Geerlings and De Proft group also explored the influence of external electric and magnetic field on the directionality of reactivity indices27,28. Furthermore, atomic forces29, force constants30, linear response functions31 and so on were analyzed and derived as more or less reactivity indicators. There are also some excellent theoretical methods that need to locate transition state to predict origin of the reaction chemoselectivity, and physical quantity changes could be tracked along the intrinsic reaction coordinate32,33,34. In the past two decades, these methods have played an important role in the interpretation and prediction of chemo- and regio-selectivities of a huge number of chemical reactions35, and wave function analysis programs such as Multiwfn36 further facilitates and promotes their applications for computing reactivity indices in the condensed version. To date, a lot of successful applications of the previous reactivity indexes have been contributed greatly to explain the role of catalysts and predict the chemoselectivity of chemical reactions37,38,39,40,41,42.

The previous reactivity indexes were generally computed with less consideration of identifying specific orbital overlap directions and modes, which should be corresponding to the possible reaction directions and thus can be used to predict the chemoselectivity and regioselectivity for the chemical reactions. As we know, the electron/charge transfer via the orbital overlap has a direction or several competing directions in the chemical reactions, so the ideal reactivity indexes should be computed with the identification of the orbital overlap direction, which is highly related to the reaction direction. As shown in Fig. 1a, if there is a type of C(sp2) atom in a reactant, that has only one possible orbital overlap direction, the computation of reactivity indexes without consideration of the direction should be proper. However, a considerable number of the substrates such as carbenes and linearly conjugated molecules need to consider the multiple reaction directions, which is due to the different possible orbital overlap directions (Fig. 1b). Therefore, developing a general method for computing reactivity indexes with the consideration of orbital overlap directions should be highly desirable in modern chemistry.

Fig. 1: The possible electron transfer directions associated with nucleophilic (Nu)/electronic (E) attack reaction directions on different types of molecules.
figure 1

a The molecules with one reaction direction. b The molecules with the atoms associated with multiple reaction directions. The red arrows represent nucleophilic (Nu) attack directions, while blue arrows represent electronic (E) attack directions.

Full size image

In this study, we introduce a directional orbital analysis method called Projection of Orbital Coefficient Vector (POCV) for calculating π electron properties and reactivity vectors. Our goal is to achieve quantitative and directional predictions of chemoselectivity in chemical reactions by consideration of valence orbital overlap modes. Although the purpose for predicting reactivity indexes is the same with CDFT, the concept of POCV is rather related to Hartree-Fock theory. Various substrates, including non-planar axially chiral molecules, carbenes, and linear conjugated molecules with multiple reaction directions, were selected as subjects of investigation. We compared the POCV method against traditional reactivity index methods and the Selected π Molecular Orbital (SMO) method to demonstrate its generality and superiority. Additionally, we developed a software program named Pywfn, implemented in Python 3.10, to apply these computational methods for computing the directional reactivity indexes (i.e., reactivity vector). We anticipate that this work will significantly impact quantum chemistry calculations43,44,45,46,47,48 related to reaction predictions and the π electron properties.

Results

Validation of reliability of POCV by computing π electronic properties

The π bond order changes of a specific axially chiral molecule, 2’-hydroxy-[1,1’-biphenyl]−2,5-dione depicted in Fig. 2a, were investigated by scanning through 36 steps at intervals of 10° for the dihedral angle ɸ(C12-C10-C3-C2). As shown in Fig. 2b, the π bond order change curves computed using the SMO method appears irregular compared to the widely accepted Mayer bond order (M)49, which does not differentiate between σ and π bonds. The SMO method struggles to identify certain π molecular orbitals in non-planar molecules50, as detailed in Supplementary Fig. S1.

Fig. 2: The π bond order changes along the flexible scanning of dihedral angle ɸ(C12-C10-C3-C2) of the axially chiral molecule 1.
figure 2

a The green circle, blue square, orange triangle, and brown rhombus symbols represent the bond orders associated with the C10-C11, C10-C12, C10-C3, and C3-C2 bonds in the selected molecule 1. b Selected π molecular orbital (SMO)-computed π bond order P changes. c Projection of orbital coefficient vector (POCV)-computed π bond order Pπ changes. d Value changes of Mayer bond order (M) minus one. Source data are provided as a Source Data file.

Full size image

Subsequently, the π bond order changes for the same molecule have been computed by using the POCV method. The resulting curve shown in Fig. 2c, is considerably smoother than that obtained by using the SMO method. As depicted in Fig. 2d, it aligns more closely with the π bond order values derived from the Mayer bond order (M) minus one, which occasionally yields negative π bond order values. Moreover, π bond order change curves computed using different methods and basis sets were tested and shown in Supplementary Figs. S1–S6, demonstrating consistent results and the tiny influence of computational method choice on the POCV method. This analysis underscores the POCV method’s efficacy in accurately assessing π electronic properties such as π bond order in both planar and non-planar molecular structures.

The POCV method can be used to calculate the π electron population of atoms in conjugated molecules. Figure 3 depicts the π electron population of selected conjugated molecules with different delocalized π systems. The results computed at the B3LYP/6-31 G level are generally correct and almost the same to the theoretical values, indicating the POCV method should be reliable. In addition, the π electron population of molecule could also be used to compute the condensed Fukui function (CFF) and condensed dual descriptor (CDD) in future.

Fig. 3: Applications of the projection of orbital coefficient vector (POCV) method for computing π electron populations of the selected molecules, taking molecular symmetry into account.
figure 3

Theoretical values of delocalized π conjugated systems are marked in red, where the superscripts denote the number of π electrons and the subscripts indicate the atomic centers involved in the selected molecules. Computed π electron populations are shown in blue, while the computed contributions of individual atoms to the π electron population are highlighted in pink. All values are reported in unit of e, where e represents the atomic unit of charge.

Full size image

In theory, there are two π bond directions for a linear conjugated molecule. As shown in Fig. 4, the π bond orders in different directions of C1 = C2 and C2 = C3 bonds in a simple allene molecule were calculated using the POCV method at the B3LYP/6-31 G level. The maximum value direction of the π bond order of C1 = C2 is perpendicular to that of C2 = C3, which indicates that the π bond direction can be computed successfully by the POCV method, especially when there are multiple π bonds associated with one main group atom in a molecule. It can be concluded both quantitative and directional prediction on the π electron properties can be achieved by using the POCV method based on the QM calculations.

Fig. 4: Predicting π bond directions of the linear molecule 8.
figure 4

a The 3D structure of 8 in which the blue and yellow arrows are respectively corresponding to 0° and 90°. b The Radar chart indicates projection of orbital coefficient vector (POCV)-computed π bond order (Pπ) in different directions of C1 = C2 (blue line) and C2 = C3 (orange line) bonds in the selected molecule 8. Source data are provided as a Source Data file.

Full size image

Aromaticity is an important property of molecules, and there have been many related studies11,51,52,53. The POCV method’s reliability on computing π bond orders prompted us to propose a method for evaluating the aromaticity of the molecule, and the standard deviation (SD) of π bond orders within a conjugated molecule has been proposed for this purpose. Figure 5a illustrates that a larger SD value corresponds to a weaker degree of electronic delocalization, and the threads of SD value of each ring are consistent with the widely used Nucleus-Independent Chemical Shifts (NICS) calculations54,55. Notably, the POCV method accurately predicts the π bond orders and electronic delocalization degrees of planar molecules, underscoring its general applicability for computing π electronic properties in conjugated systems.

Fig. 5: Measuring aromaticity using the projection of orbital coefficient vector (POCV) method.
figure 5

a Standard deviation (SD) value of POCV-computed π bond orders (Pπ) of planar molecules. b SD value of POCV-computed π bond orders (Pπ) of non-planar molecules. The POCV-computed π bond orders (Pπ) are represented using a sequential color map, where blue indicates higher POCV-computed π bond orders and white indicates lower POCV-computed π bond orders. Bonds shown in black color have been excluded from the SD calculations. The nucleus-independent chemical shift (NICS) values of corresponding molecular rings with numbered circles are listed and also compared with their SD values in Supplementary Table S1.

Full size image

The POCV method could also be used to compute aromaticity for non-planar molecules, offering a different approach from the traditional methods like NICS. Figure 5b shows the computed aromaticity of the non-planar molecules including axial chiral molecule and the previously reported N-heterocyclic carbene (NHC)-stabilized diazoborane molecule56, which is a kind of challenging synthetic molecule that could be stabilized by aromaticity.

Prediction on reactivity vectors of carbenes and linear conjugated molecules

The POCV method uniquely incorporates directional information on orbital overlap, enabling precise calculations of chemical reactivity indices across different directions. This capability lets us quantitatively predict potential reaction directions for each active atom that is usually unsaturated main group atom with π electron in complex molecules.

For both nucleophilic and electrophilic carbenes shown in Fig. 6, it can be successfully predicted the max directional nucleophilic and electrophilic vectors \(\overrightarrow{{{{{\bf{F}}}}}_{{{{\bf{m}}}}{{{\bf{a}}}}{{{\bf{x}}}}}^{-}}\) and \(\overrightarrow{{{{{\bf{F}}}}}_{{{{\bf{m}}}}{{{\bf{a}}}}{{{\bf{x}}}}}^{+}}\) of carbene carbon by using POCV method. For the common nucleophilic carbene such as N-heterocyclic carbene57,58,59,60,61 depicted in Fig. 6a, the nucleophilic direction corresponds to the direction of lone pair σ2 electrons, which can be switched as the electrophilic direction in the electrophilic carbene (Fig. 6d) with the σ0π2 electronic configuration62,63.

Fig. 6: Prediction on reactivity vectors of carbenes 17 and 18.
figure 6

a The 2D graph of nucleophilic carbene (σ2π0). b Predicting the nucleophilic vectors \(\overrightarrow{{{{{\bf{F}}}}}^{-}}\) (red arrows) and the corresponding maximum value \(\overrightarrow{{{{{\bf{F}}}}}_{{{{\bf{m}}}}{{{\bf{a}}}}{{{\bf{x}}}}}^{-}}\) (yellow arrow). c The distribution of reactivity vectors \(\overrightarrow{{{{\bf{F}}}}}\) (gray circle), \(\overrightarrow{{{{{\bf{F}}}}}^{+}}\) (blue square), and \(\overrightarrow{{{{{\bf{F}}}}}^{-}}\) (red triangle) associated with the directions in the semicircle area selected by the same angle with the two C-N bonds. d The 2D graph of electronic carbene (σ0π2). e The electrophilic vector \(\overrightarrow{{{{{\bf{F}}}}}^{+}}\) (blue arrows) and the corresponding maximum value \(\overrightarrow{{{{{\bf{F}}}}}_{{{{\bf{m}}}}{{{\bf{a}}}}{{{\bf{x}}}}}^{+}}\) (yellow arrow). f The distribution of reactivity vectors \(\overrightarrow{{{{\bf{F}}}}}\) (gray circle), \(\overrightarrow{{{{{\bf{F}}}}}^{+}}\) (blue square), and \(\overrightarrow{{{{{\bf{F}}}}}^{-}}\) (red triangle) associated with the directions in the semicircle area selected by the same angle with the two C-P bonds. Source data are provided as a Source Data file.

Full size image

To further validate the generality of the POCV method, we selected an example of an intermediate with the presence of a linearly conjugate molecule based on the previous work64, and performed the atomic reactivity vector calculations. As shown in Fig. 7a, there are two possible nucleophilic sites including α- and γ-carbon atoms with the non-directional nucleophilic Parr function and Fukui function are 0.315/0.417 and 0.042/0.058, but the corresponding energy barriers of α- and γ-addition are 14.8 kcal·mol⁻¹ and 19.9 kcal·mol⁻¹, respectively, indicating the non-directional nucleophilicity should be not used for predicting the active sites for the linear conjugated molecules. Subsequently, as shown in Fig. 7b, we have additionally calculated the maximum value of atomic reactivity vector (\(\overrightarrow{{{{{\bf{F}}}}}_{{{{\bf{m}}}}{{{\bf{a}}}}{{{\bf{x}}}}}}\)) of α- and γ-carbon atoms are 1.27 and 1.17, respectively, which is consistent with the above energy barrier trend. That is to say, the greater atomic reactivity vector is associated with the lower energy barrier, which confirms the accuracy of the POCV method in these kinds of systems.

Fig. 7: Predicting reactivity vectors of linearly conjugated molecules with multiple reaction sites.
figure 7

a The energy barriers of α- and γ-addition reaction pathways at different reaction sites. b The reactivity vectors of different reaction sites of intermediate molecule 19 in panel a, and the corresponding maximum values \((\overrightarrow{{{{{\bf{F}}}}}_{{{{\bf{m}}}}{{{\bf{a}}}}{{{\bf{x}}}}}})\) are marked beside. Source data are provided as a Source Data file.

Full size image

Discussion

In contrast to previous computational methods that neglect molecular orbital overlap directions and modes, our study introduces the POCV method for calculating π bond orders and reactivity vector indices in selected conjugated molecules. To validate the efficacy of POCV, we conducted π electron population and π bond order calculations on conjugated molecules, and assessed electronic delocalization degree by computing the standard deviation of π bond orders across planar and non-planar structures. Additionally, we proposed the reactivity vector can be used to more accurately predict reaction sites in the carbenes and linear conjugated molecules with multiple orbital overlap directions. Our findings affirm that the POCV method offers robust rationality, accuracy, and versatility in chemical predictions. It could also be proposed that the POCV method can be connected to and used in the CDFT calculations. POCV might even be used to predict the orbital interaction of two fragments, and this is inspired by the analysis method of Natural Orbitals for Chemical Valence (NOCV)65,66,67,68 originated by Nalewajski and Mrozek, which is allowed to look into the details of the orbital deformation density in terms of NOCV orbitals or densities to understand the nature of the orbital interaction better. In summary, POCV is a method for analyzing directional orbital interactions. At present, the projection method of p-type orbital has been proven effective in this work, and the directional analysis method for interactions of d-type and other type orbitals is still being studied in our laboratory. Future work will focus on further methodological exploration, software development, and establishing a reactivity vectors database to facilitate precise directional chemical reaction design by using POCV.

Methods

All molecular structures used in this study were optimized using the Gaussian16 program69 at the B3LYP70/6-31 G(d) level. Relevant information such as Cartesian coordinates, orbital coefficients, and overlap matrices of the molecules were extracted for computing directional reactivity indices and π electron properties using the POCV method implemented in the Pywfn software. The computational results using different computational methods and basis sets show consistent results, and more details can be found in Supporting Information.

POCV concept

According to the quantum calculations performed in Gaussian, the self-consistent field equation of a molecule in Eq. (1) can be solved to obtain the wave functions (i.e., molecular orbitals) with the form of the linear combination of atomic orbitals (LCAO) in Eq. (2). For each of the molecular orbital, the linear combination of px, py, pz atomic orbitals can be equal to a coefficient vector \(\overrightarrow{{{{\bf{p}}}}}({C}_{px},{C}_{px},{C}_{px})\) multiplied by the corresponding basis functions in Eq. (3). The expressions of px, py, pz atomic orbitals have been provided according to the basis sets in Supporting Information, and a theoretical derivation has been contributed as following to answer why the three coefficients of p-type atomic orbitals can be proposed as a vector.

$$FC=SCE$$
(1)
$${\psi }_{n}={\sum }_{i=1}^{{N}_{b}}{C}_{i,n}{\varphi }_{i}$$
(2)
$${\psi }_{n}=\cdots+C{}_{s}\varphi _{s}+\overrightarrow{p}({C}_{px},{C}_{py},{C}_{pz})\cdot \left[\begin{array}{c}{\varphi }_{px}\\ {\varphi }_{py}\\ {\varphi }_{pz}\end{array}\right]+\cdots$$
(3)

The possible reaction directions should correspond to possible overlap directions of the valence p atomic orbitals of the main group atoms, which are selected as the projection directions area \(\{\overrightarrow{n}\}\) of \(\overrightarrow{{{{\bf{p}}}}}\) vector and have to be ≥90° with the existed σ bonds, avoiding the contribution of the projection of the \(\overrightarrow{{{{\bf{p}}}}}\) vector that is belong to the existed σ bonds. The flow chart of Pywfn is depicted in Fig. 8a. As shown in Fig. 8b, the defined only one projection direction \(\overrightarrow{{{{\bf{n}}}}}\) for a planar C(sp2) atom maintains equal angles that are 90° with the existing σ bonds, a linear C(sp) atom has projection directions \(\{\overrightarrow{n}\}\) spread across a circular area, and carbene-type C(sp2) atom has their projection directions area \(\{\overrightarrow{n}\}\) distributed within the spherical wedge. The POCV method is supposed to be mainly used for handling the unsaturated main group atoms with π electron, and more methods for selecting orbital projection directions for other types of main group atoms will be supplemented in the future with updates to the pywfn software.

Fig. 8: Details of the projection of orbital coefficient vector (POCV) methodology.
figure 8

a The Pywfn flowchart. Quantum mechanics (QM)-derived Cartesian coordinates and Canonical molecular orbital coefficient matrices (C) can be used to define the projection directions \(\{\overrightarrow{{{{\bf{n}}}}}\}\) for the atoms selected by the user. Based on these inputs, the POCV algorithm is applied to replace each original vector \(\overrightarrow{{{{\bf{p}}}}}\) by the projected vector \(\overrightarrow{{{{\bf{p}}}}^{\prime}}\) of user-selected atoms in all the occupied molecule orbital. This process generates two types of projected molecular orbital coefficient matrices (Cπ and Cf), enabling the calculation of physical quantities including π electron properties, π electron based condensed Fukui function (CFF-π) and condensed dual descriptor (CDD-π), reactivity vectors, and so on. b The possible project directions \(\{\overrightarrow{{{{\bf{n}}}}}\}\) of different types of atoms including planar, linear, and carbene-type conjugated atoms, in which blue and green arrows represent initial projection direction \(\overrightarrow{{{{{\bf{n}}}}}_{{{{\bf{0}}}}}}\) and the possible reaction direction \(\overrightarrow{{{{\bf{n}}}}}\), respectively. c Schematic representation of the modification of molecular orbital coefficient matrices: The p atomic orbital coefficient vectors of all occupied molecular orbitals would be updated from \(\overrightarrow{{{{\bf{p}}}}}\) to \(\overrightarrow{{{{\bf{p}}}}^{\prime}}\) by using POCV method. d Projection diagrams for POCV-computed π bond order between A and B atoms \(({P}_{AB}^{\pi })\). e Projection diagrams for computing bound bond order between A and B atoms \(({P}_{AB}^{f})\).

Full size image

Since the direction of the orbital coefficient vector \(\overrightarrow{{{{\bf{p}}}}}({C}_{px},{C}_{px},{C}_{px})\) is not always the same as the projection direction \(\overrightarrow{{\mathbf{n}} }\), we have to compute the projection of the orbital coefficient vector \(\overrightarrow{{{{\bf{p}}}}}\) in the direction \(\overrightarrow{{{{\bf{n}}}}}\) to represent the contribution of atomic orbital in the molecule. As shown in Fig. 8c, the \(\overrightarrow{{{{\bf{p}}}}}\) of the related atoms could be projected in the selected direction \(\overrightarrow{{\mathbf{n}} }\) as a vector \(\overrightarrow{{{{\bf{p}}}}^{\prime}}(C^{\prime}_{px},C^{\prime}_{px},C^{\prime}_{px})\) by using Eq. (4), and a theoretical derivation of the POCV concept and the projection Eq. (4) can be found as following.

As shown in Eq. (5) and Fig. 9, the p-type atomic orbital (\({C}_{px}^{B}{\varphi }_{px}(r)+{C}_{py}^{B}{\varphi }_{py}(r)+{C}_{pz}^{B}{\varphi }_{pz}(r)\)) under the base coordinate system B(\(\overrightarrow{{{{\bf{x}}}}}\),\(\overrightarrow{{{{\bf{y}}}}}\),\(\overrightarrow{{{{\bf{z}}}}}\)) should be equal to that (\({C}_{px}^{E}{\varphi }_{px}(rT)+{C}_{py}^{E}{\varphi }_{py}(rT)+{C}_{pz}^{E}{\varphi }_{pz}(rT)\)) under the base coordinate system E(\(\overrightarrow{{{{\boldsymbol{\sigma }}}}}\),\(\overrightarrow{{{{\boldsymbol{\gamma }}}}}\),\(\overrightarrow{{{{\bf{n}}}}}\)). The \(r(x,y,z)\) represents the coordinates of any point in space, and it can be transformed to \(r^{\prime}(x^{\prime},y^{\prime},z^{\prime})\) by multiplying a customized matrix \(T=[\overrightarrow{{{{\boldsymbol{\sigma }}}}},\overrightarrow{{{{\boldsymbol{\gamma }}}}},\overrightarrow{{{{\bf{n}}}}}]\) in Eq. (6). Then, we need to find the transformation relationship between the orbital coefficients \({C}_{px}^{B},{C}_{py}^{B},{C}_{pz}^{B}\) and the updated coefficients \({C}_{px}^{E},{C}_{py}^{E},{C}_{pz}^{E}\). The x′, y′, z′ can be written as those in Eq. (7). The base functions φ(r) in three directions of the p-type atomic orbital are shown in Eq. (8), and the r and r′ can be substituted into the Eq. (8), whereas the orbital angular momentums (l,m,n) are (1,0,0), (0,1,0), and (0,0,1) for \({\varphi }_{px}\), \({\varphi }_{py}\), and \({\varphi }_{pz}\) respectively to obtain Eqs. (9) and (10). Since base coordinate system transformation does not change the length of r vector, i.e., \({x}^{2}+{y}^{2}+{z}^{2}={x^{\prime} }^{2}+{y^{\prime} }^{2}+{z^{\prime} }^{2}\), so we can ignore the changes of exponential functions and normalization constant N, and the following theoretical derivation can be continued with the simplified basis function in Eqs. (11, 12). Then, we let the simplified basis functions (Eqs. (11, 12)) substitute into Eq. (5) to obtain Eq. (13), and we can also obtain Eq. (14), in which it can be concluded that there is a multiplication relationship between the transform matrix T and the two sets of coefficients. Therefore, the two sets of coefficients can be reasonably proposed as two vectors in mathematics (Eq. (15)). Due to the orthogonality among these three unit vectors (\(\overrightarrow{{{{\boldsymbol{\sigma }}}}},\overrightarrow{{{{\boldsymbol{\gamma }}}}},{{{\rm{and}}}}\,\overrightarrow{{{{\bf{n}}}}}\)), we can conclude \({T}^{T}={T}^{-1}\), and obtain the coefficient vector in base coordinate system E in Eq. (16). Since we only retain the \(\overrightarrow{{{{\bf{n}}}}}\) directional p-type orbital, the p-type orbital coefficient vector \(\overrightarrow{{{{{\bf{P}}}}}^{{{{\bf{E}}}}}}({C}_{px}^{E},{C}_{py}^{E},{C}_{pz}^{E})\) in base coordinate system E can be reduced to the only one in \(\overrightarrow{{{{\bf{n}}}}}\) direction, which can be further transformed to the wave function in the base coordinate system B through the following Eq. (17). Subsequently, the projection equation of \({{{{\bf{p}}}}}^{{{{\bf{B}}}}}\) (i.e., the \(\overrightarrow{{{{\bf{p}}}}}\) vector) in \(\overrightarrow{{{{\bf{n}}}}}\) direction can be reasonably proved in Eq. (4), and the reasonability of POCV method could be confirmed by the above theoretical derivation. At last, \(\overrightarrow{{{{\bf{p}}}}^{\prime}}(C^{\prime}_{px},C^{\prime}_{px},C^{\prime}_{px})\) (i.e., \(\overrightarrow{{{{{\bf{P}}}}}^{{{{\bf{E}}}}}}({C}_{px}^{E},{C}_{py}^{E},{C}_{pz}^{E})\)) can be used to replace the original p-type atomic orbital coefficients vector \(\overrightarrow{{{{\bf{p}}}}}({C}_{px},{C}_{px},{C}_{px})\) to obtain the updated orbital coefficients matrixes (\({C}^{\pi }/{C}^{f}\)) and density matrixes \(({D}^{\pi }/{D}^{f})\) computed by Eq. (18).

$$\overrightarrow{{{{\bf{p}}}}^{\prime}}(C^{\prime}_{px},C^{\prime}_{py},C^{\prime}_{pz})=\overrightarrow{{{{\bf{p}}}}}\cdot \overrightarrow{{{{\bf{n}}}}}\cdot \overrightarrow{{{{\bf{n}}}}}$$
(4)
$${C}_{px}^{B}{\varphi }_{px}(r)+{C}_{py}^{B}{\varphi }_{py}(r)+{C}_{pz}^{B}{\varphi }_{pz}(r)={C}_{px}^{E}{\varphi }_{px}(rT)+{C}_{py}^{E}{\varphi }_{py}(rT)+{C}_{pz}^{E}{\varphi }_{pz}(rT)$$
(5)
$$r\cdot T=[\begin{array}{ccc}x & y & z\end{array}] \left[\begin{array}{ccc}{T}_{11} & {T}_{12} & {T}_{13}\\ {T}_{21} & {T}_{22} & {T}_{23}\\ {T}_{31} & {T}_{31} & {T}_{33}\end{array} \right]=[\begin{array}{ccc}{x}^{\prime} & {y}^{\prime} & {z}^{\prime} \end{array}]={r}^{\prime}$$
(6)
$$\begin{array}{c}x^{\prime}={T}_{11}x+{T}_{21}y+{T}_{31}z\\ y^{\prime}={T}_{12}x+{T}_{22}y+{T}_{32}z\\ z^{\prime}={T}_{13}x+{T}_{23}y+{T}_{33}z\end{array}$$
(7)
$$\varphi (r)=N\cdot {x}^{l}\cdot {y}^{m}\cdot {z}^{n}\cdot {e}^{a({x}^{2}+{y}^{2}+{z}^{2})}$$
(8)
$$\begin{array}{c}{\varphi }_{px}(r)=N\cdot x\cdot {e}^{a({x}^{2}+{y}^{2}+{z}^{2})}\\ {\varphi }_{py}(r)=N\cdot y\cdot {e}^{a({x}^{2}+{y}^{2}+{z}^{2})}\\ {\varphi }_{pz}(r)=N\cdot z\cdot {e}^{a({x}^{2}+{y}^{2}+{z}^{2})}\end{array}$$
(9)
$${\varphi }_{px}(r ^{\prime})=N\cdot x^{\prime}\cdot {e}^{a(x{^{\prime 2}}+y{ ^{\prime 2}}+z{^{\prime 2}} )}\\ {\varphi }_{py}(r^{\prime})=N\cdot {y}^{\prime}\cdot {e}^{a(x{^{\prime 2}}+{y}{^{\prime 2}}+z{^{\prime 2}} )}\\ {\varphi }_{pz}(r^{\prime})=N\cdot z^{\prime}\cdot {e}^{a(x{^{\prime 2}}+y{^{\prime 2}}+z{^{\prime 2}} )}$$
(10)
$$\begin{array}{c}{\varphi }_{px}(r)=x\\ {\varphi }_{py}(r)=y\\ {\varphi }_{pz}(r)=z\end{array}$$
(11)
$$\begin{array}{c}{\varphi }_{px}(r{\prime} )=x^{\prime}={T}_{11}x+{T}_{21}y+{T}_{31}z\\ {\varphi }_{py}(r^{\prime} )=y^{\prime}={T}_{12}x+{T}_{22}y+{T}_{32}z\\ {\varphi }_{pz}(r^{\prime} )=z^{\prime}={T}_{13}x+{T}_{23}y+{T}_{33}z\end{array}$$
(12)
$${C}_{px}^{B}x+{C}_{py}^{B}y+{C}_{pz}^{B}z= {C}_{px}^{E}({T}_{11}x+{T}_{21}y+{T}_{31}z)+{C}_{py}^{E}({T}_{12}x+{T}_{22}y+{T}_{32}z) \\ +{C}_{pz}^{E}({T}_{13}x+{T}_{23}y+{T}_{33}z)=({C}_{px}^{E}{T}_{11}+{C}_{py}^{E}{T}_{12}+{C}_{pz}^{E}{T}_{13})x \\ +({C}_{px}^{E}{T}_{21}+{C}_{py}^{E}{T}_{22}+{C}_{pz}^{E}{T}_{23})y+({C}_{px}^{E}{T}_{31}+{C}_{py}^{E}{T}_{32}+{C}_{pz}^{E}{T}_{33})z$$
(13)
$$\left[\begin{array}{ccc}{T}_{11} & {T}_{12} & {T}_{13}\\ {T}_{21} & {T}_{22} & {T}_{23}\\ {T}_{31} & {T}_{32} & {T}_{33}\end{array}\right] \left[\begin{array}{c}{C}_{px}^{E}\\ {C}_{py}^{E}\\ {C}_{pz}^{E}\end{array}\right]=\left[\begin{array}{c}{C}_{px}^{B}\\ {C}_{py}^{B}\\ {C}_{pz}^{B}\end{array}\right]$$
(14)
$$T\cdot \overrightarrow{{{{{\bf{P}}}}}^{{{{\bf{E}}}}}}=\overrightarrow{{{{{\bf{P}}}}}^{{{{\bf{B}}}}}}$$
(15)
$$\overrightarrow{{{{{\bf{P}}}}}^{{{{\bf{E}}}}}}={T}^{-1}\cdot \overrightarrow{{{{{\bf{P}}}}}^{{{{\bf{B}}}}}}={T}^{T}\cdot \overrightarrow{{{{{\bf{P}}}}}^{{{{\bf{B}}}}}}$$
(16)
$${C}_{pz}^{E}{\varphi }_{pz}( {r}^{\prime} )= [({T}_{13},{T}_{23},{T}_{33})\cdot ({C}_{px}^{B},{C}_{py}^{B},{C}_{pz}^{B})]\cdot [({T}_{13},{T}_{23},{T}_{33}) \\ \cdot ({\varphi }_{px}(r),{\varphi }_{py}(r),{\varphi }_{pz}(r))]=[\overrightarrow{{{{\bf{n}}}}}\cdot \overrightarrow{{P}^{B}}]\cdot \overrightarrow{{{{\bf{n}}}}}\cdot [{\varphi }_{px}(r),{\varphi }_{py}(r),{\varphi }_{pz}(r)] \\= \overrightarrow{{{{\bf{p}}}} ^{\prime}} (C^{\prime}_{x},C^{\prime}_{y},C^{\prime}_{z}) \cdot [{\varphi }_{px}(r),{\varphi }_{py}(r),{\varphi }_{pz}(r)]$$
(17)
$${D}_{\mu \nu }^{\pi /f}={\eta }_{i}{\sum}_{i}{C}_{\mu i}^{\pi /f}{C}_{\nu i}^{\pi /f}$$
(18)
Fig. 9
figure 9

Transformation of p-type atomic orbitals that is the combination of basis functions in different coordinate systems. The same p-type atomic orbital in orange color formed by the basic functions in red, gray, and blue colors under the two base coordinate systems \({{{\bf{B}}}}(\overrightarrow{{{{\bf{x}}}}},\overrightarrow{{{{\bf{y}}}}},\overrightarrow{{{{\bf{z}}}}})\) and \({{{\bf{E}}}}(\overrightarrow{{{{\boldsymbol{\sigma }}}}},\overrightarrow{{{{\boldsymbol{\gamma }}}}},\overrightarrow{{{{\bf{n}}}}})\) transformed by using linearity transformation matrix T[\(\overrightarrow{{{{\boldsymbol{\sigma }}}}}\),\(\overrightarrow{{{{\boldsymbol{\gamma }}}}}\),\(\overrightarrow{{{{\bf{n}}}}}\)]. The small black balls represent two bonded atoms A and B, \(\overrightarrow{{{{\boldsymbol{\sigma }}}}}\) vector represents the unit vector in the bond axis direction, whereas \(\overrightarrow{{{{\boldsymbol{\gamma }}}}}\) and \(\overrightarrow{{{{\bf{n}}}}}\) vectors are also the unit vectors in corresponding directions, and the three vectors are orthogonal to each other. The Cpx, Cpy, and Cpz denote the coefficients of basis functions in the global coordinate system B, whereas C′px, C′py, and C′pz represent the corresponding coefficients of basis functions in the local coordinate system E.

Full size image

π electronic properties

Firstly, all the unsaturated main group atoms in a molecule would be identified, then the possible pπ electron population direction for each unsaturated atom was selected as projection direction \(\overrightarrow{{{{\bf{n}}}}}\). Noteworthy, the unsaturated atoms of non-planar molecules could have the different \(\overrightarrow{{{{\bf{n}}}}}\) directions that are always ≥90° with the existed σ bonds, and this is why we can handle the non-planar conjugated molecules. When there is sp3-type hybridized main group atom without π electron, it cannot be treated as the center A atom shown in Fig. 8d, and its contribution is ignored for computing the π electron properties. Then, all the \(\vec{n}\) projection direction vectors of the unsaturated main group atoms including \(\overrightarrow{{{{\bf{p}}}} ^{\prime}}\)(A) and \(\overrightarrow{{{{\bf{p}}}} ^{\prime}}\)(B) in Fig. 8d have been obtained and only retained in the corresponding projected molecule orbital coefficient matrix \({C}^{\pi }\). At last, π bond order (\({P}_{AB}^{\pi }\)) and π electron population \({Q}_{A}^{\pi }\) could be further computed according to Eq. (19) and Eq. (20). Combining the π electron population with some concepts from CDFT (such as the CDD) yields the CDD-π according to Eq. (21).

$${P}_{AB}^{\pi }={\sum}_{\mu \in A}{\sum}_{\nu \in B}{({D}^{\pi }S)}_{\mu \nu }{({D}^{\pi }S)}_{\nu \mu }$$
(19)
$${Q}_{A}^{\pi }={Z}_{A}-{\sum}_{\mu \in A}{({D}^{\pi }S)}_{\mu \mu }$$
(20)
$$CD{D}^{\pi }={Q}_{A}^{\pi (N-1)}+{Q}_{A}^{\pi (N+1)}-{Q}_{A}^{\pi (N)}$$
(21)

Reactivity vector indexes

Generally, the reactivity vector direction (i.e., the most possible reaction direction) of an unsaturated main group atom is highly related to the upcoming formation direction of σ bond through the overlap of s- and p-orbitals, which should be already included in the above projection directions \(\{\overrightarrow{{{{\bf{n}}}}}\}\). In addition, the reactivity vector value should be associated with the bond order values contributed by s-orbital and projected p-orbital of the selected A atom with all the adjacent atoms (Fig. 8e), which has been defined as the total bound bond order and can be computed as follows. It should be mentioned that main group atoms have very few d-type and other type orbital contributions, which can be directly retained in the bound bond order calculations.

Since the s-type orbitals that cannot be separated and also contribute to the σ bonds, both the s-type orbitals and the projected vector \(\overrightarrow{{{{\bf{p}}}}^{\prime}}\) for the A atom depicted in Fig. 8e were retained in the corresponding projected molecule orbital coefficient matrix \({C}^{f}\), in which all the orbital coefficients of the adjacent atoms were also retained. Subsequently, the bound bond order (\({P}_{AB}^{f}\)) in Fig. 8e can be computed according to the Eq. (22) using matrix \({D}^{f}\). The atomic reactivity vector \(\overrightarrow{{{{\bf{F}}}}}\) can be computed based on Eq. (23), in which the \({P}_{\max }\) is defined the max bond order of carbon atom A (\({P}_{\max }\)(C) = 3) in a single direction. Hereinafter, \(\overrightarrow{{{{{\bf{F}}}}}^{+}}\) and \(\overrightarrow{{{{{\bf{F}}}}}^{-}}\) are defined as the electrophilic and nucleophilic reactivity vectors, and can be obtained in single-point energy calculations on the optimized molecule that additionally adds or subtracts one electron, respectively. For the sp3-type hybridized main group atom, its four bond axis directions could be selected as the four possible reaction directions \(\{\overrightarrow{{{{\bf{n}}}}}\}\), and the bound bond orders could still be computed to further obtain the reactivity vector.

$${P}_{AB}^{f}={\sum}_{\mu \in A}{\sum}_{\nu \in B}{({D}^{f}S)}_{\mu \nu }{({D}^{f}S)}_{\nu \mu }$$
(22)
$$\overrightarrow{{F}_{A}}=[{M}_{\max}-{\sum}_{B}{P}_{AB}^{f}]\overrightarrow{\cdot n}$$
(23)

For the function \(F(\overrightarrow{{{{\bf{n}}}}})\), we need to solve iteratively using the gradient descent method to find the \(\overrightarrow{{{{{\bf{F}}}}}_{{{{\bf{m}}}}{{{\bf{a}}}}{{{\bf{x}}}}}}\) vector associated with the possible reaction direction. The iteration stops when \(F(\overrightarrow{{{{\bf{n}}}}})\) reaches the maximum value or when the condition is not satisfied, at which point the final \(\overrightarrow{{{{{\bf{F}}}}}_{{{{\bf{m}}}}{{{\bf{a}}}}{{{\bf{x}}}}}}\) is obtained. Noteworthy, there are three initial guesses for the projection directions n (i.e., \(\overrightarrow{{{{{\bf{n}}}}}_{{{{\bf{0}}}}}}\)) as follows. As shown in Fig. 8b, when an atom A has three neighboring atoms (i.e., B1, B2, and B3), \(\overrightarrow{{{{{\bf{n}}}}}_{{{{\bf{0}}}}}}\) is the two normal vectors of the plane formed by the three neighboring atoms with\(\overrightarrow{{{{{\bf{n}}}}}_{{{{\bf{0}}}}}}={\mbox{Norm}}(\overrightarrow{{{{{\bf{B}}}}}_{{{{\bf{1}}}}}{{{{\bf{B}}}}}_{{{{\bf{2}}}}}}\times \overrightarrow{{{{\bf{B}}}}{}_{{{{\bf{1}}}}}{{{\bf{B}}}}_{{{{\bf{2}}}}}})\); When an atom A has two neighboring atoms (B1, B2) and in the same line, \(\overrightarrow{{{{{\bf{n}}}}}_{{{{\bf{0}}}}}}\) is any direction perpendicular to the line with \(\overrightarrow{{{{{\bf{n}}}}}_{{{{\bf{0}}}}}}={\mbox{Norm}}(\overrightarrow{{{{{\bf{B}}}}}_{{{{\bf{1}}}}}{{{{\bf{B}}}}}_{{{{\bf{2}}}}}}\times \overrightarrow{{{{\bf{N}}}}})\), in which \(\overrightarrow{{{{\bf{N}}}}}\) is a random unit vector; When an atom A has two neighboring atoms (B1, B2) and not in the same line, \(\overrightarrow{{{{{\bf{n}}}}}_{{{{\bf{0}}}}}}\) is the sum of the two vectors from atom A to (B1, B2) with\(\overrightarrow{{{{{\bf{n}}}}}_{{{{\bf{0}}}}}}={\mbox{Norm}}(\overrightarrow{{{{\bf{A}}}}{{{{\bf{B}}}}}_{{{{\bf{1}}}}}}\times \overrightarrow{{{{\bf{A}}}}{{{{\bf{B}}}}}_{{{{\bf{2}}}}}})\).

Selecting π molecular orbital (SMO) method

In order to compare with the π bond order results computed by using the selecting π molecular orbital (SMO) and POCV methods, we have devised a method to realize the selection of π molecular orbital and calculated the π bond order as following. The initial step of the SMO method is to specify π electron population direction, which is the normal vector direction of the plane where the C(sp2)-type atom is located. To facilitate the selection, we have set some rules so that the computer can automatically select π molecular orbitals, and the results are almost the same to those picked by human. The two rules of SMO method are (1) the contribution of 1 s orbital is much smaller and the contribution of all p orbital is larger than that of s orbital, and (2) the angle between p orbital maximum value direction of the C(sp2)-type atom and the specified normal vector direction is smaller than 0.2 π. All the occupied molecular orbitals that satisfy above two conditions would be selected as π molecular orbitals. Finally, we have calculated the π bond order for the bond between atoms A and B based on the Mayer bond order equation by using all the selected occupied π molecular orbitals.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.