The behavior of electrons during bond formation and breaking cannot commonly be accessed from experiments. Thus, bond perception is often based on chemical intuition or rule-based algorithms. Utilizing computational chemistry methods, we present intrinsic bond descriptors for the Diels–Alder reaction, allowing for an automatic bond perception. We show that these bond descriptors are available from localized orbitals and self-interaction correction calculations, e.g., from Fermi-orbital descriptors. The proposed descriptors allow a sparse, simple, and educational inspection of the Diels–Alder reaction from an electronic perspective. We demonstrate that bond descriptors deliver a simple visual representation of the concerted bond formation and bond breaking, which agrees with Lewis’ theory of bonding.

Bond type perception, the identification of bond types in a given molecule, is a longstanding problem in many modern-day applications, including visualization,1 force field calculations,2 or multiscale methods.3 While a human interpretation of bonding situations appears trivial in many cases, an ab initio description of bonding based on computational chemistry can be more challenging. Over the years, various simple descriptions of bonding have been developed, such as Lewis’ theory4 or Linnett’s double-quartet theory.5–8 Alternatively, heuristic criteria are often used for bond perception, e.g., the atomic distance being smaller than the sum of covalent radii. Such criteria can be used to determine bonded atoms, while additional rule-based algorithms can be employed to detect bonds to next-nearest neighbors in the local environment of an atom.2,9,10

Depending on the application, such criteria may be too simple to accurately describe the bond order for any molecule of interest. More sophisticated methods are, for example, Bader’s atoms in molecules analysis,11,12 where the Laplacian of the electronic density is analyzed, or the electron localization function13 that also accounts for the Pauli exclusion principle. Regardless, bond orders cannot directly be obtained from either method. Other methods provide bond order analysis, e.g., Wiberg-14 or Mayer-bond orders.15 Nonetheless, one often obtains fractional bond orders, making a classification into single, double, or triple bonds cumbersome.

An alternative description of bond orders can be found with the use of so-called Fermi-orbital descriptors (FODs).16 In recent research, it has been discussed that the FODs carry chemical bonding information and can be interpreted according to Lewis’ or Linnett’s double-quartet theory.16,17 Furthermore, it has been shown that the center of mass (centroids) of different localized orbitals can be used to get a reasonable initial guess for the FOD positions16 and to track bond rearrangements, e.g., by creating curly arrows as a representation of bond reorganization during a reaction.18,19

Sometimes, it may be desirable to track the bonding situation or even the bond formation and breaking of molecules in a reaction. This paper presents an analysis of bond descriptors for a chemical reaction, i.e., bond formation and bond breaking along a reaction path using automatic electron perception methods. As an example, we will model a well-known reaction—the Diels–Alder reaction.20–22 Despite often being used as a textbook example for a concerted reaction mechanism where all bond changes happen simultaneously,23–26 this reaction is still very important in modern applications, such as self-healing materials.27,28 In the simplest form of this reaction, 1,3-butadiene reacts with ethylene to form cyclohexene, as displayed in Fig. 1. In this work, we will provide insights into the bond formation for the Diels–Alder reaction from a bond perception and self-interaction perspective.

FIG. 1.

Schematic Diels–Alder reaction of 1,3-butadiene (C4H6) and ethylene (C2H4) to cyclohexene (C6H10). Two new carbon-carbon bonds will be created by breaking two double bonds and relocating another one. All bond changes happen simultaneously due to the concerted nature of the reaction.

FIG. 1.

Schematic Diels–Alder reaction of 1,3-butadiene (C4H6) and ethylene (C2H4) to cyclohexene (C6H10). Two new carbon-carbon bonds will be created by breaking two double bonds and relocating another one. All bond changes happen simultaneously due to the concerted nature of the reaction.

Close modal

This paper is structured as follows: In Sec. II, the theoretical background is introduced. The computational details and methodology to model the Diels–Alder reaction are described in Sec. III. The results are presented in Sec. IV. The calculated FODs and centroids will be shown and compared. Afterward, energies and other properties that have been calculated along the reaction path will be presented. The corresponding information on bond formation and breaking is discussed. Subsequently, the validity of the qualitative results is verified. The results are summarized in Sec. V.

Given its moderate computational effort and sufficient accuracy, density functional theory (DFT)29,30 has become one of the most commonly employed methods in computational chemistry and related research fields.31 The total energy of a system is given by the Kohn–Sham (KS) energy functional,
EKS[nα,nβ]=Ts[nα,nβ]+Vext[n]+EH[n]+EXC[nα,nβ],
(1)
with Ts being the kinetic energy, Vext being the external potential, EH being the Coulomb energy, and EXC being the exchange–correlation energy. In an open-shell description of a molecular system and related unrestricted calculations, there are two spin channels σ, i.e., the α and β spin channels. Accordingly, n denotes the total electronic density, while nσ denotes the spin density.
While being formally exact, for practical DFT calculations, additional assumptions regarding the exchange-correlation functional are needed. With those assumptions and further numerical approximations, reasonable results for a large variety of systems can be obtained.32 However, those density functional approximations can also introduce unwanted side effects, such as the artificial interaction of electrons with themselves.33 Since for one-electron systems, the Coulomb energy (EH) and exchange–correlation contributions (EXC) have to cancel each other out, violating this condition will result in the so-called self-interaction error (SIE) ESI,
ESI[n1σ]=EH[n1σ]+EXC[n1σ,0],
(2)
where n1σ is a single-electron density for an electron with spin σ. Recent research34 showed that the one-electron error,35,36 i.e., an error closely related to SIE, is still dominant in modern exchange–correlation functionals.
In the formulation of Perdew–Zunger (PZ),33 the SIE is removed from the KS energy EKS for all Nσ occupied orbitals, resulting in a self-interaction correction (SIC),
EPZ[nα,nβ]=EKS[nα,nβ]σiNσESI[niσ].
(3)
The FLO-SIC method37–40 is a novel form of PZ-SIC, where Fermi–Löwdin orbitals (FLOs) are used to calculate the one-electron densities. Fermi orbitals (FOs) can be obtained by applying the transformation matrix Rσ to the molecular orbital coefficients Cσ,
cFO,σ=CσRσ.
(4)
The transformation matrix Rσ can be obtained with
RijFO,σ=ψjσaiσnσ(aiσ).
(5)
Here, ψi are occupied orbitals, while aiσ denote position eigenstates localized at the FODs aiσ, which are reference electron positions to build the FOs. In FLO-SIC calculations, these FODs are optimized along with the electronic density. The FLOs are obtained by applying Löwdin’s symmetrical orthonormalization41 to the FOs. The resulting FLOs are well-localized, comparable to other localized orbitals, e.g., Foster–Boys (FB) orbitals42 where the orbital variances JFB are minimized,43,
JFB=iψir2ψiψirψi2.
(6)
The lower the value of JFB, the more localized is the state. It is known that SIC raises reaction energy barriers.44–46 For localized SIC orbitals, this can be explained by the increased noding of valence orbitals in the transition state (TS) that lowers the SIE (or raises the correction).47 The novel aspect of our contribution is the explicit analysis of how SIC and information derived from localized orbitals can be used to track the bond formation along a reaction path.

For transparent results, only open-source codes were used in this manuscript following the free and open-source software (FOSS) approach.48 The all-electron Gaussian-type orbital code PySCF49 was used for all DFT calculations, while the PyFLOSIC2 code50,51 was used for all FLO-SIC calculations. The computational parameters for these calculations follow Ref. 17. The local spin density approximation Slater exchange functional52,53 with the modified Perdew–Wang correlation functional54 (SPW92) has been used for all calculations, as implemented in the LIBXC library.55 The double-ζ polarization consistent pc-1 basis set56–58 was employed. A PySCF grid level of seven was used. This corresponds to a (90 974) grid for hydrogen and a (135 1202) grid for carbon in the multi-center quadrature scheme.59 Pruning has been disabled for all calculations. The self-consistent field (SCF) energy convergence threshold was set to 1 × 10−8 Eh. All calculations were performed spin-unrestricted.

To verify the results, DFT and FLO-SIC calculations have been performed with a self-written Julia code called chilli.jl.60 A comparison of the calculated energies can be found in the supplementary material. The nuclear geometry optimizations, TS optimization, and nudged elastic bands (NEB) calculations were carried out with the external optimizer pysisyphus61 using an interface to PySCF. The computational parameters and methodology used will be explained in detail below.

In the target reaction, the reactant state contains the separated molecules ethylene and 1,3-butadiene. Along a reaction path, the reactant state transforms via a transition state into the product state, cyclohexene. Nuclear geometry optimizations were started from experimental geometries: The reactant ethylene was taken from Ref. 62. For butadiene, the gauche-1,3 geometry was taken from Ref. 63. The structure for the product, cyclohexene, was taken from Ref. 64. The initial structures were optimized at the DFT level of theory using the rational function optimization (RFO) method.65 The threshold for the maximum force Fmax and the root-mean-square force Frms have been set to 1.5 × 10−5 Eh/a0 and 1 × 10−5 Eh/a0, respectively.

Subsequently, the optimized nuclear geometries for ethylene and butadiene were combined, with the initial fragment distance being similar to Ref. 66. The resulting reactant state was optimized utilizing translational and rotational intrinsic coordinates,67 using the same force thresholds as before. It has been verified that the resulting Hessian of all nuclear geometry optimizations has only positive eigenvalues to ensure that stable minima have been found.

To approximate the TS, a NEB calculation has been utilized at the DFT level of theory. The computational parameters were adapted from Refs. 68 and 69. Eleven images have been calculated using the limited-memory BFGS (L-BFGS) optimizer,70,71 as implemented in pysisyphus.61 A maximum force of 2 × 10−4 Eh/a0 (≈0.01 eV/Å) and a step size of 0.075  a0 (≈0.04 Å) were set. The climbing image was enabled, while the Kabsch algorithm72 to align images has been enabled. The splined highest energy image (HEI) from the NEB calculation was used to optimize the TS. This optimization uses the same thresholds as the previous nuclear geometry optimizations using a restricted-step RFO (RS-I-RFO) method.65,73 After the TS optimization, it has been checked that the resulting Hessian contains solely one imaginary frequency to ensure the result truly is a TS.

In addition, two NEB calculations were performed with the same parameters as the previous NEB calculation but with nine images and the climbing image disabled. One NEB samples the minimum energy path between the reactants and the TS, and the other samples between the TS and the product. This procedure resulted in 21 optimized geometries (one reactant, one TS, one product, and 2 × 9 NEB images).

It is known that Kohn–Sham orbitals (KSOs) tend to be delocalized over the whole extent of a given molecule.40 Thus, different localization schemes can be used to transform the KSOs into localized orbitals. In contrast to KSOs, localized orbitals can often be interpreted as orbitals carrying bond-related information.74 The investigated localized orbitals are FB, Edmiston–Ruedenberg (ER),75 Pipek–Mezey (PM) using Löwdin charges,76 generalized PM using Becke charges,77 and intrinsic bond orbitals (IBOs).78 All of these localization procedures were carried out using PySCF, where the co-iterative augmented Hessian (CIAH) method is used,79 except for the IBOs, where the implementation and minimization follow Ref. 78 with an additional symmetric orthogonalization. As a simple stability analysis for the minimization with the CIAH method, the eigenvalues of the resulting Hessian have been calculated, and it has been ensured that all eigenvalues are positive.

For the FLO-SIC calculations, the initial FOD configurations were taken from the PyCOM method,16 i.e., the centroids of localized FB orbitals. Note that the method was used with a stability analysis to ensure the local minima of the generated localized orbitals. The initial FOD configurations follow Lewis’ theory. With the need to optimize the density matrix and the FODs, the two-step SCF cycle following Ref. 80 was employed as implemented in PyFLOSIC2. The constrained L-BFGS-B method81 was used for the FOD optimization. The final maximum force Fmax acting on any FOD was below 1.5 × 10−4 Eh/a0.

In this section, we present an overview of the main results, focusing on the bond analysis during the Diels–Alder reaction, employing spatial descriptors and other selected properties.

This segment investigates the trends that can be found when analyzing the centroids of localized orbitals and FODs. For simplicity, the centroids of localized orbitals and FODs positions will be labeled as associated bond points (ABPs). As an example, the FODs have been displayed for selected structures in Fig. 2, along with the density. Corresponding figures using different centroids can be found in the supplementary material.

FIG. 2.

Diels–Alder reaction of 1,3-butadiene (C4H6) and ethylene (C2H4) to cyclohexene (C6H10). Electronic density (top) and FODs (bottom) for (a) the reactants, (b) the TS, (c) the image after the TS (TS + 1), and (d) the product. The isovalues have been chosen such that 85% of the density is contained. More information about the selection of the isovalue can be found in the supplementary material. Carbon atoms are colored in gray, hydrogen in white, and the FODs in red. Only one spin channel has been displayed since the FOD positions for these images coincide for both spin channels. The centroids of localized orbitals resemble the FODs (see the supplementary material).

FIG. 2.

Diels–Alder reaction of 1,3-butadiene (C4H6) and ethylene (C2H4) to cyclohexene (C6H10). Electronic density (top) and FODs (bottom) for (a) the reactants, (b) the TS, (c) the image after the TS (TS + 1), and (d) the product. The isovalues have been chosen such that 85% of the density is contained. More information about the selection of the isovalue can be found in the supplementary material. Carbon atoms are colored in gray, hydrogen in white, and the FODs in red. Only one spin channel has been displayed since the FOD positions for these images coincide for both spin channels. The centroids of localized orbitals resemble the FODs (see the supplementary material).

Close modal

For the reactant [Fig. 2(a)] and product state [Fig. 2(d)], the ABPs correctly describe the anticipated bonding situation (compare with Fig. 1). Counting the ABPs around each bond axis mimics the expected bond order16 (see Fig. 1). Note that the number of ABPs has to be divided by two since we are in the spin-unrestricted case. The TS [Fig. 2(b)] still has the bond order of the reactant state. Over the course of the reaction, the bond order will change once, i.e., from the reactant state bond order to the product state bond order. This transition takes place one image after the TS [Fig. 2(c)] for all bonds and can be associated with the formation and breaking of bonds. This imagery of all bonds forming concertedly in one step is also in line with other research.23–26 We note that this behavior remains when evoking other localization methods, with small deviations. While the FODs and centroids of FLO, FB, and ER orbitals are separated perpendicular to the bond axis, the centroids of PM orbitals, generalized PM orbitals, and IBOs mostly lie directly on the bond axis.

As a measure of the change of the ABPs during the reaction, one can use the distance dmax, where dmax is the maximum distance of the ABPs of bonding orbitals from the midpoint of the C–C bonds. As the formation of two new C–C bonds, alongside the breakup of two C–C double bonds over the course of the reaction, necessitates a spatial rearrangement of localized electrons pairs, this can be directly observed using dmax. A value of dmax ≈ 0 would resemble a single bond, where the FOD positions would directly lie on the C–C bond axis. Typical distances for dmax for double bonds, i.e., in ethylene, are about 0.4 Å (see the supplementary material for more details).

In Fig. 3(a), one can see the dmax for the FODs for all breaking double bonds, starting from the initial reactant state toward the product structure. The FOD distance dmax remains almost constant in the beginning at 0.4 Å, as typical for double bonds. Approaching the TS, dmax increases, reaching its peak at the TS (image 10). This means that FODs move farther away from the C–C bond axes. One can interpret the departing ABPs as the breaking of a bond. After the TS, the double bonds are broken into single bonds, with dmax being close to zero. This trend can be seen for all ABPs. For the PM orbitals and IBOs, one can see an oscillation of the distances for the centroids in ethylene. In addition, one can see that some orbitals behave differently depending on the spin channel, especially for PM orbitals. Thus, in these cases, the symmetry is broken between α and β spin channel. However, this behavior cannot be found in the generalized PM localization scheme. Since the trends between PM and generalized PM are otherwise the same, the oscillations and spin splitting are likely an artifact of the ill-defined Löwdin charges.77 

FIG. 3.

Maximum FOD distances dmax from the C–C bond midpoints (a) for the breaking double bonds and (b) for the forming double bond. The TS is indicated at image 10 via a dashed vertical line. For comparison, the long dashed horizontal line indicates the FOD distance in ethylene. The color-coded FOD positions in the molecule are indicated schematically, where carbon atoms are colored in gray and hydrogen in white.

FIG. 3.

Maximum FOD distances dmax from the C–C bond midpoints (a) for the breaking double bonds and (b) for the forming double bond. The TS is indicated at image 10 via a dashed vertical line. For comparison, the long dashed horizontal line indicates the FOD distance in ethylene. The color-coded FOD positions in the molecule are indicated schematically, where carbon atoms are colored in gray and hydrogen in white.

Close modal

As shown in the supplementary material, the values for dmax for the FB centroids closely resemble the ones of the FLOs. Thus, for the discussed reaction, PyCOM utilizing FB orbitals is suggested to result in suitable initial FODs. The second best approximation for initial FODs utilizing PyCOM would be achieved utilizing ER. PM, generalized PM, and IBOs perform worse for the investigated reaction.

The formation of the new double bond in the product can be tracked as well, as seen in Fig. 3(b). For the initial single bond, dmax is close to zero. The C–C single bond is present up to the TS. The peak of dmax is at the TS + 1 where the new double bond has been formed. The distance stabilizes afterward and reaches a value close to the previously observed 0.4 Å that is typical for a double bond.

In this segment, scalar properties, e.g., total energies, SIEs, and the absolute value of the dipole moment, are analyzed along the calculated reaction path. We restrict ourselves to scalar properties, aiming to find simple scalar descriptors for bond formation and bond breaking.

The total KS energies EKS for the sampled reaction are displayed in Fig. 4(a), shifted by the total energy of the product. As expected, one can indicate the TS as the HEI. Calculating the PZ SIE from the FLO-SIC calculations [see Eq. (3)], one obtains Fig. 4(b). First, one can see that the SIE is larger for the product than for the reactant. Interestingly, the SIE gets minimal for the TS. However, the smaller SIE for the TS is in line with the expectation from Ref. 47. Applying this correction to the KS energies, one obtains the FLO-SIC energies EPZ, as seen in Fig. 4(a). The reaction barrier increases for the corrected energies, agreeing with other research.45,46

FIG. 4.

Energies along the reaction images. The TS is indicated at image 10 via a dashed line.(a) Total KS energies EKS and FLO-SIC energies EPZ, shiftedby the respective total energy of the product. Note that theFLO-SIC energies are in any case lower than the KS energies.See the supplementary material for the unshifted energies.(b) FLO-SIC SIEs ESI. The respective structures for the reactants, TS, and product are sketched in (b) accordingly, with carbon atoms colored in gray, hydrogen in white, and the FODs in red. Only one spin channel has been displayed since the FOD positions for these images coincide for both spin channels.

FIG. 4.

Energies along the reaction images. The TS is indicated at image 10 via a dashed line.(a) Total KS energies EKS and FLO-SIC energies EPZ, shiftedby the respective total energy of the product. Note that theFLO-SIC energies are in any case lower than the KS energies.See the supplementary material for the unshifted energies.(b) FLO-SIC SIEs ESI. The respective structures for the reactants, TS, and product are sketched in (b) accordingly, with carbon atoms colored in gray, hydrogen in white, and the FODs in red. Only one spin channel has been displayed since the FOD positions for these images coincide for both spin channels.

Close modal

Figure 5 displays the orbital variances JFB [see Eq. (6)] of the occupied FB orbitals and FLOs. The FLOs are more localized than the FB orbitals but both show a similar trend. The figure roughly resembles the trend of the negative SIE in Fig. 4(b). Accordingly, as the SIE gets minimal at the TS, the TS is also the most delocalized state as indicated by the maximal value of the FB cost function (Fig. 5). The bond breaking, as indicated previously by the ABPs [see Fig. 3(a)], can also be seen in the orbital variances of all localized orbitals, which provide a measure of delocalization [see Eq. (6)]. The orbital variances of the remaining localized orbitals can be found in the supplementary material.

FIG. 5.

Orbital variances JFB of occupied FB orbitals and FLOs along the reaction images. The TS is indicated at image 10 via a dashed line.

FIG. 5.

Orbital variances JFB of occupied FB orbitals and FLOs along the reaction images. The TS is indicated at image 10 via a dashed line.

Close modal

Similar behavior can be seen in other properties, such as the absolute electric dipole moment μ for both the KS and the FLO-SIC densities, shown in Fig. 6, or the ionization potential (see the supplementary material). As proposed in Ref. 17, the absolute value of the dipole moment is a simple descriptor of the density. We now find that this simple descriptor is also able to correctly monitor density changes leading to bond breaking and formation in this reaction. The dipole moment can indicate the re-arrangement of the density due to the relocation of the electrons, which is the actual bond breaking.

FIG. 6.

Absolute electric dipole moments μ of the systems along the reaction images, both for the DFT and the FLO-SIC densities. The TS is indicated at image 10 via a dashed line.

FIG. 6.

Absolute electric dipole moments μ of the systems along the reaction images, both for the DFT and the FLO-SIC densities. The TS is indicated at image 10 via a dashed line.

Close modal

Additionally, we investigated the possibility to monitor bond changes using universal force field (UFF) energies.82 In Fig. 7, one can see the calculated UFF energies EUFF using PyFLOSIC2 and Open Babel.83  PyFLOSIC2 utilizes a novel bond perception based on the FOD positions. Using nearest-neighbor relations between nuclei and FODs, each FOD gets classified as core, bond, or lone FOD. This information enables the calculation of the bond order matrix and the determination of the local chemical environment of each atom in a molecule. Since the centroids give the same bond order as the FODs, they will result in the same UFF energies.

FIG. 7.

UFF energies EUFF of the systems along the reaction images using PyFLOSIC2 and Open Babel. The TS is indicated at image 10 via a dashed line.

FIG. 7.

UFF energies EUFF of the systems along the reaction images using PyFLOSIC2 and Open Babel. The TS is indicated at image 10 via a dashed line.

Close modal

Interestingly, the bond assessment for the UFF based on FODs delivers energies resembling the trend of the dipole moment μ (compare to Fig. 6). While for the reactant and product state, the ABP-based energies agree with the Open Babel derived ones, there are significant differences around the TS. The ABP-based bond perception describes the concerted change of the bond order better than the derived bonds from Open Babel. Further investigations are necessary for a generalization of these findings. Despite this, the ABP-based UFF seems to be a promising candidate for an efficient scalar function tracking bond formation and breaking for similar systems.

To verify the qualitative trend of the found results, additional DFT calculations have been performed on the optimized structures, similar to the verification segment in Ref. 17. To analyze the effect of the exchange–correlation functional, calculations have been performed utilizing the generalized gradient approximation (GGA) functional PBEsol84 and the meta-GGA functional r2SCAN.85 The main result is that independent of the functional, one can still see the bond formation in the dipole moment and the ionization potential. Moreover, the orbital variances of FB orbitals display the same trends, while the dmax of FB orbitals are mostly unaffected by the choice of the functional.

To analyze the effect of the basis set, DFT calculations using the family of polarization consistent basis sets pc-n56–58 have been performed. In particular, the split-valence, double-ζ, triple-ζ, quadruple-ζ, and quintuple-ζ basis sets pc-0, pc-1, pc-2, pc-3, and pc-4 have been investigated. One finds a suitable agreement for pc-1 with the larger basis sets. Similar to the analysis of the functional, the dmax of FB orbitals remain largely unaffected by the basis set choice. Accordingly, all basis sets recover the same trend for the dipole moment, ionization potential, and orbital variances.

Additionally, to confirm that there is no bias in the values of dmax due to the starting guess of the FOD optimization, an exemplary calculation using ER orbitals has been performed. As expected, both configurations result in the same optimized FOD configuration and, therefore, in the same values of dmax, while the FB guess needs fewer iterations to converge.

Finally, to ensure the reproducibility of our findings, reference DFT and FLO-SIC calculations were performed with chilli.jl. The calculated DFT and FLO-SIC energies along the reaction path show the same trend as obtained by PyFLOSIC2. This accentuates that the presented FLO-SIC trends are reproducible. The full analysis of all reference calculations can be found in the supplementary material.

In this article, we investigated bond formation and breaking exemplary for the well-known Diels–Alder reaction from a self-interaction and local orbital perspective, using nudged elastic bands (NEB) calculations and optimizing the corresponding transition state (TS) using density functional theory (DFT). It has been shown that the Fermi-orbital descriptors (FODs) and centroids of different localized orbitals, i.e., Fermi–Löwdin orbitals (FLOs), Foster–Boys (FB), and Edmiston–Ruedenberg (ER), can be used to describe and display the bonding situation not only in the reactant and product state but along the reaction as well. All associated bond points (ABPs) show a similar trend and indicate the bond formation and breaking at the same reaction image. More precisely, we find that the bond order changes exactly once along the reaction path while going from the TS to image TS + 1. This clearly shows that the concerted mechanism of the Diels–Alder reaction can be described using ABPs. It has been shown that the (maximum) ABPs distance from the middle point of the bond axis can be used to indicate the bond formation and breaking. We showed that for this reaction, FB orbitals are the best starting point from the tested centroids to generate initial FODs, utilizing the PyCOM procedure.

DFT and FLO-SIC calculations have been performed along the same reaction. In agreement with other research,45,46 it has been shown that the reaction barrier gets raised when applying the self-interaction correction (SIC) and the self-interaction error (SIE) gets minimal for the TS. The bond formation can be observed in different properties, such as the absolute value of the dipole moment or the universal force field (UFF) energy derived from ABPs.

For the investigated reaction, a Lewis configuration for the FODs was found. For any reaction having a Lewis configuration, one can expect similar behavior for the proposed measures of tracking bond changes. However, there exist cases where instead of Lewis configurations, Linnett double-quartet structures are preferred in FLO-SIC.17,51 For these cases, the FLOs and contained bonding information are vastly different from the localized orbitals of DFT calculations. Thus, it could be interesting to apply the proposed measures to investigate reactions involving Linnett double-quartet structures.

The main finding of this work is that fully optimized FODs and centroids of localized orbitals are able to describe the bond formation and bond breaking for the analyzed Diels–Alder reaction. While other research has hinted at the bond information contained in FOD positions, the centroids are often more accessible. Such descriptors can be helpful to categorize the bonding situation in various applications, e.g., when analyzing self-healing materials. Therefore, it is promising that they show similar trends. Whether these trends continue for more complex reactions should be investigated in future research.

Figures that display all ABPs for selected structures can be found in the supplementary material. The summarized results of the verification calculations using chilli.jl and the calculations using different functionals and basis sets are discussed. Additionally, figures for the values of dmax for all ABPs, the ionization potentials, and UFF contributions from PyFLOSIC2 and Open Babel are listed. A reference implementation to calculate the isovalues from Fig. 2 is included. More information about the symmetry of the optimized nuclear geometries, the localization procedure, bond perception, and the definition of the electric dipole moment are appended.

The authors thank Dr. Johannes Steinmetzer for helpful discussions and technical support regarding pysisyphus. The authors thank the Universitätsrechenzentrum of the Friedrich Schiller University Jena for computational time and support. W.T. Schulze and S. Gräfe highly acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)–Research unit FuncHeal, project ID 455748945–FOR 5301 (project P5). S. Schwalbe has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)–Project ID 421663657–KO 1924/9-2. S. Schwalbe thanks Dr. Sebastian Borrmann for technical support and the HPC Freiberg for computational time. The authors thank Dr. Susi Lehtola for comments on the original manuscript. The authors thank an anonymous reviewer whose comments guided us to fix the FB cost function in PySCF.

The authors have no conflicts to disclose.

Wanja Timm Schulze: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Software (supporting); Visualization (lead); Writing – original draft (lead). Sebastian Schwalbe: Conceptualization (equal); Data curation (supporting); Formal analysis (supporting); Investigation (supporting); Methodology (equal); Software (lead); Supervision (equal); Validation (lead); Visualization (supporting); Writing – review & editing (equal). Kai Trepte: Conceptualization (equal); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Supervision (supporting); Visualization (supporting); Writing – review & editing (equal). Alexander Croy: Conceptualization (equal); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Project administration (equal); Supervision (supporting); Visualization (supporting); Writing – review & editing (equal). Jens Kortus: Conceptualization (equal); Funding acquisition (lead); Visualization (supporting); Writing – review & editing (equal). Stefanie Gräfe: Conceptualization (equal); Funding acquisition (lead); Project administration (equal); Resources (lead); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material. The initial and optimized structures, NEB images, and images with ABPs included can be found openly available under https://gitlab.com/wangenau/bond_formation_supplementary.

1.
F.
Lazzari
,
A.
Salvadori
,
G.
Mancini
, and
V.
Barone
, “
Molecular perception for visualization and computation: The Proxima library
,”
J. Chem. Inf. Model.
60
,
2668
(
2020
).
2.
Q.
Zhang
,
W.
Zhang
,
Y.
Li
,
J.
Wang
,
L.
Zhang
, and
T.
Hou
, “
A rule-based algorithm for automatic bond type perception
,”
J. Cheminf.
4
,
26
(
2012
).
3.
P.
Seeber
,
S.
Seidenath
,
J.
Steinmetzer
, and
S.
Gräfe
, “
Growing Spicy ONIOM and many body expansions
,”
WIREs Comput. Mol. Sci.
(published online)
.
4.
G. N.
Lewis
, “
The atom and the molecule
,”
J. Am. Chem. Soc.
38
,
762
(
1916
).
5.
J. W.
Linnett
, “
Valence-bond structures: A new proposal
,”
Nature
187
,
859
(
1960
).
6.
J. W.
Linnett
, “
A modification of the Lewis-Langmuir octet rule
,”
J. Am. Chem. Soc.
83
,
2643
(
1961
).
7.
W. F.
Luder
, “
Electronic structure of molecules (Linnett, J. W.)
,”
J. Chem. Educ.
43
,
55
(
1966
).
8.
P. B.
Empedocles
and
J. W.
Linnett
, “
The electronic structure of benzene
,”
Proc. R. Soc. London, Ser. A
282
,
166
(
1964
).
9.
J.
Wang
,
W.
Wang
,
P. A.
Kollman
, and
D. A.
Case
, “
Automatic atom type and bond type perception in molecular mechanical calculations
,”
J. Mol. Graphics Modell.
25
,
247
(
2006
).
10.
S.
Artemova
,
L.
Jaillet
, and
S.
Redon
, “
Automatic molecular structure perception for the universal force field
,”
J. Comput. Chem.
37
,
1191
(
2016
).
11.
R. F. W.
Bader
, “
A quantum theory of molecular structure and its applications
,”
Chem. Rev.
91
,
893
(
1991
).
12.
R. F. W.
Bader
,
Atoms in Molecules, International Series of Monographs on Chemistry
(
Clarendon Press
,
1994
).
13.
A. D.
Becke
and
K. E.
Edgecombe
, “
A simple measure of electron localization in atomic and molecular systems
,”
J. Chem. Phys.
92
,
5397
(
1990
).
14.
K. B.
Wiberg
, “
Application of the pople-santry-segal CNDO method to the cyclopropylcarbinyl and cyclobutyl cation and to bicyclobutane
,”
Tetrahedron
24
,
1083
(
1968
).
15.
I.
Mayer
, “
Bond order and valence indices: A personal account
,”
J. Comput. Chem.
28
,
204
(
2007
).
16.
S.
Schwalbe
,
K.
Trepte
,
L.
Fiedler
,
A. I.
Johnson
,
J.
Kraus
,
T.
Hahn
,
J. E.
Peralta
,
K. A.
Jackson
, and
J.
Kortus
, “
Interpretation and automatic generation of Fermi-orbital descriptors
,”
J. Comput. Chem.
40
,
2843
(
2019
).
17.
K.
Trepte
,
S.
Schwalbe
,
S.
Liebing
,
W. T.
Schulze
,
J.
Kortus
,
H.
Myneni
,
A. V.
Ivanov
, and
S.
Lehtola
, “
Chemical bonding theories as guides for self-interaction corrected solutions: Multiple local minima and symmetry breaking
,”
J. Chem. Phys.
155
,
224109
(
2021
).
18.
J. E. M. N.
Klein
,
G.
Knizia
, and
H. S.
Rzepa
, “
Epoxidation of alkenes by peracids: From textbook mechanisms to a quantum mechanically derived curly-arrow depiction
,”
ChemistryOpen
8
,
1244
(
2019
).
19.
G.
Sciortino
,
A.
Lledós
, and
P.
Vidossich
, “
Bonding rearrangements in organometallic reactions: From orbitals to curly arrows
,”
Dalton Trans.
48
,
15740
(
2019
).
20.
O.
Diels
and
K.
Alder
, “
Synthesen in der hydroaromatischen Reihe
,”
Justus Liebigs Ann. Chem.
460
,
98
(
1928
).
21.
K.
Alder
, “
Neuere Methoden der präparativen organischen Chemie. 16. Die Methode der Diensynthese
,”
Angew. Chem.
55
,
53
(
1942
).
22.
E. B.
Hershberg
and
J. R.
Ruhoff
, “
1,3-Butadiene
,”
Org. Synth.
17
,
25
(
1937
).
23.
R. E.
Townshend
,
G.
Ramunni
,
G.
Segal
,
W. J.
Hehre
, and
L.
Salem
, “
Organic transition states. V. The Diels-Alder reaction
,”
J. Am. Chem. Soc.
98
,
2190
(
1976
).
24.
K. N.
Houk
,
Y. T.
Lin
, and
F. K.
Brown
, “
Evidence for the concerted mechanism of the Diels-Alder reaction of butadiene with ethylene
,”
J. Am. Chem. Soc.
108
,
554
(
1986
).
25.
F.
Bernardi
,
A.
Bottoni
,
M. J.
Field
,
M. F.
Guest
,
I. H.
Hillier
,
M. A.
Robb
, and
A.
Venturini
, “
MC-SCF study of the Diels-Alder reaction between ethylene and butadiene
,”
J. Am. Chem. Soc.
110
,
3050
(
1988
).
26.
E.
Goldstein
,
B.
Beno
, and
K. N.
Houk
, “
Density functional theory prediction of the relative energies and isotope effects for the concerted and stepwise mechanisms of the Diels–Alder reaction of butadiene and ethylene
,”
J. Am. Chem. Soc.
118
,
6036
(
1996
).
27.
R.
Geitner
,
J.
Kötteritzsch
,
M.
Siegmann
,
T. W.
Bocklitz
,
M. D.
Hager
,
U. S.
Schubert
,
S.
Gräfe
,
B.
Dietzek
,
M.
Schmitt
, and
J.
Popp
, “
Two-dimensional Raman correlation spectroscopy reveals molecular structural changes during temperature-induced self-healing in polymers based on the Diels–Alder reaction
,”
Phys. Chem. Chem. Phys.
17
,
22587
(
2015
).
28.
C. R.
Ratwani
,
A. R.
Kamali
, and
A. M.
Abdelkader
, “
Self-healing by Diels-Alder cycloaddition in advanced functional polymers: A review
,”
Prog. Mater. Sci.
131
,
101001
(
2023
).
29.
P.
Hohenberg
and
W.
Kohn
, “
Inhomogeneous electron gas
,”
Phys. Rev.
136
,
B864
(
1964
).
30.
W.
Kohn
and
L. J.
Sham
, “
Self-consistent equations including exchange and correlation effects
,”
Phys. Rev.
140
,
A1133
(
1965
).
31.
T.
van Mourik
,
M.
Bühl
, and
M.-P.
Gaigeot
, “
Density functional theory across chemistry, physics and biology
,”
Philos. Trans. R. Soc., A
372
,
20120488
(
2014
).
32.
K.
Burke
, “
Perspective on density functional theory
,”
J. Chem. Phys.
136
,
150901
(
2012
).
33.
J. P.
Perdew
and
A.
Zunger
, “
Self-interaction correction to density-functional approximations for many-electron systems
,”
Phys. Rev. B
23
,
5048
(
1981
).
34.
S.
Schwalbe
,
K.
Trepte
, and
S.
Lehtola
, “
How good are recent density functionals for ground and excited states of one-electron systems?
,”
J. Chem. Phys.
157
,
174113
(
2022
).
35.
D. R.
Lonsdale
and
L.
Goerigk
, “
The one-electron self-interaction error in 74 density functional approximations: A case study on hydrogenic mono- and dinuclear systems
,”
Phys. Chem. Chem. Phys.
22
,
15805
(
2020
).
36.
D. R.
Lonsdale
and
L.
Goerigk
, “
One-electron self-interaction error and its relationship to geometry and higher orbital occupation
,”
J. Chem. Phys.
158
,
044102
(
2023
).
37.
M. R.
Pederson
,
A.
Ruzsinszky
, and
J. P.
Perdew
, “
Communication: Self-interaction correction with unitary invariance in density functional theory
,”
J. Chem. Phys.
140
,
121103
(
2014
).
38.
M. R.
Pederson
, “
Fermi orbital derivatives in self-interaction corrected density functional theory: Applications to closed shell atoms
,”
J. Chem. Phys.
142
,
064112
(
2015
).
39.
M. R.
Pederson
and
T.
Baruah
, “
Self-interaction corrections within the Fermi-orbital-based formalism
,” in
Advances in Atomic, Molecular, and Optical Physics
(
Elsevier
,
2015
), pp.
153
180
.
40.
Z. h.
Yang
,
M. R.
Pederson
, and
J. P.
Perdew
, “
Full self-consistency in the Fermi-orbital self-interaction correction
,”
Phys. Rev. A
95
,
052505
(
2017
).
41.
P. O.
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.
18
,
365
(
1950
).
42.
J. M.
Foster
and
S. F.
Boys
, “
Canonical configurational interaction procedure
,”
Rev. Mod. Phys.
32
,
300
(
1960
).
43.
D. A.
Kleier
,
T. A.
Halgren
,
J. H.
Hall
, and
W. N.
Lipscomb
, “
Localized molecular orbitals for polyatomic molecules. I. A comparison of the Edmiston-Ruedenberg and Boys localization methods
,”
J. Chem. Phys.
61
,
3905
(
1974
).
44.
B. G.
Johnson
,
C. A.
Gonzales
,
P. M. W.
Gill
, and
J. A.
Pople
, “
A density functional study of the simplest hydrogen abstraction reaction. Effect of self-interaction correction
,”
Chem. Phys. Lett.
221
,
100
(
1994
).
45.
S.
Patchkovskii
and
T.
Ziegler
, “
Improving “difficult” reaction barriers with self-interaction corrected density functional theory
,”
J. Chem. Phys.
116
,
7806
(
2002
).
46.
A. J.
Johansson
,
M. R. A.
Blomberg
, and
P. E. M.
Siegbahn
, “
Quantifying the effects of the self-interaction error in density functional theory: When do the delocalized states appear? II. Iron-oxo complexes and closed-shell substrate molecules
,”
J. Chem. Phys.
129
,
154301
(
2008
).
47.
C.
Shahi
et al, “
Stretched or noded orbital densities and self-interaction correction in density functional theory
,”
J. Chem. Phys.
150
,
174102
(
2019
).
48.
S.
Lehtola
and
A. J.
Karttunen
, “
Free and open source software for computational chemistry education
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
12
,
e1610
(
2022
).
49.
Q.
Sun
et al, “
Recent developments in the PySCF program package
,”
J. Chem. Phys.
153
,
024109
(
2020
).
50.
S.
Schwalbe
,
L.
Fiedler
,
J.
Kraus
,
J.
Kortus
,
K.
Trepte
, and
S.
Lehtola
, “
PyFLOSIC: Python-based Fermi–Löwdin orbital self-interaction correction
,”
J. Chem. Phys.
153
,
084104
(
2020
).
51.
S.
Liebing
,
K.
Trepte
, and
S.
Schwalbe
, “
Effect of molecular and electronic geometries on the electronic density in FLO-SIC
,” in
Springer Proceedings in Physics
(
Springer International Publishing
,
2022
), pp.
167
186
.
52.
F.
Bloch
, “
Bemerkung zur Elektronentheorie des Ferromagnetismus und der elektrischen Leitfähigkeit
,”
Z. Phys.
57
,
545
(
1929
).
53.
P. A. M.
Dirac
, “
Note on exchange phenomena in the Thomas atom
,”
Math. Proc. Cambridge Philos. Soc.
26
,
376
(
1930
).
54.
J. P.
Perdew
and
Y.
Wang
, “
Accurate and simple analytic representation of the electron-gas correlation energy
,”
Phys. Rev. B
45
,
13244
(
1992
).
55.
S.
Lehtola
,
C.
Steigemann
,
M. J. T.
Oliveira
, and
M. A. L.
Marques
, “
Recent developments in libxc—A comprehensive library of functionals for density functional theory
,”
SoftwareX
7
,
1
(
2018
).
56.
F.
Jensen
, “
Polarization consistent basis sets: Principles
,”
J. Chem. Phys.
115
,
9113
(
2001
).
57.
F.
Jensen
, “
Polarization consistent basis sets. II. Estimating the Kohn–Sham basis set limit
,”
J. Chem. Phys.
116
,
7372
(
2002
).
58.
F.
Jensen
, “
Polarization consistent basis sets. III. the importance of diffuse functions
,”
J. Chem. Phys.
117
,
9234
(
2002
).
59.
A. D.
Becke
, “
A multicenter numerical integration scheme for polyatomic molecules
,”
J. Chem. Phys.
88
,
2547
(
1988
).
60.
Dataset: S.
Schwalbe
,
K.
Trepte
, and
W. T.
Schulze
,
(2022)
“chilli.jl [chilli_jl],” Zenodo https://zenodo.org/record/7596230.
61.
J.
Steinmetzer
,
S.
Kupfer
, and
S.
Gräfe
, “
pysisyphus: Exploring potential energy surfaces in ground and excited states
,”
Int. J. Quantum Chem.
121
,
e26390
(
2021
).
62.
P.
Linstrom
, NIST chemistry WebBook, NIST standard reference database 69,
1997
.
63.
J. H.
Baraban
,
M. A.
Martin‐Drumel
,
P. B.
Changala
,
S.
Eibenberger
,
M.
Nava
,
D.
Patterson
,
J. F.
Stanton
,
G. B.
Ellison
, and
M. C.
McCarthy
, “
The molecular structure of gauche-1,3-butadiene: Experimental establishment of non-planarity
,”
Angew. Chem., Int. Ed.
57
,
1821
(
2018
).
64.
J. F.
Chiang
and
S. H.
Bauer
, “
Molecular structure of cyclohexene
,”
J. Am. Chem. Soc.
91
,
1898
(
1969
).
65.
A.
Banerjee
,
N.
Adams
,
J.
Simons
, and
R.
Shepard
, “
Search for stationary points on surfaces
,”
J. Phys. Chem.
89
,
52
(
1985
).
66.
B.
Ramirez
,
T.
Cordova
,
F.
Ruette
, and
G.
Chuchani
, “
Inquiry of the reaction paths in thermal retro-Diels–Alder reactions in the gas phase: Theoretical study on the concerted and stepwise elimination mechanisms of cyclohexenes
,”
Comput. Theor. Chem.
1067
,
103
(
2015
).
67.
L.-P.
Wang
and
C.
Song
, “
Geometry optimization made simple with translation and rotation coordinates
,”
J. Chem. Phys.
144
,
214108
(
2016
).
68.
H. C.
Herbol
,
J.
Stevenson
, and
P.
Clancy
, “
Computational implementation of nudged elastic band, rigid rotation, and corresponding force optimization
,”
J. Chem. Theory Comput.
13
,
3250
(
2017
).
69.
A. W.
Ruttinger
,
D.
Sharma
, and
P.
Clancy
, “
Protocol for directing nudged elastic band calculations to the minimum energy pathway: Nurturing errant calculations back to convergence
,”
J. Chem. Theory Comput.
18
,
2993
(
2022
).
70.
J.
Nocedal
, “
Updating quasi-Newton matrices with limited storage
,”
Math. Comput.
35
,
773
(
1980
).
71.
D. C.
Liu
and
J.
Nocedal
, “
On the limited memory BFGS method for large scale optimization
,”
Math. Program.
45
,
503
(
1989
).
72.
W.
Kabsch
, “
A solution for the best rotation to relate two sets of vectors
,”
Acta Crystallogr., Sect. A
32
,
922
(
1976
).
73.
E.
Besalú
and
J. M.
Bofill
, “
On the automatic restricted-step rational-function-optimization method
,”
Theor. Chem. Acc.
100
,
265
(
1998
).
74.
J. J. P.
Stewart
, “
An examination of the nature of localized molecular orbitals and their value in understanding various phenomena that occur in organic chemistry
,”
J. Mol. Model.
25
,
7
(
2019
).
75.
C.
Edmiston
and
K.
Ruedenberg
, “
Localized atomic and molecular orbitals
,”
Rev. Mod. Phys.
35
,
457
(
1963
).
76.
J.
Pipek
and
P. G.
Mezey
, “
A fast intrinsic localization procedure applicable for ab initio and semiempirical linear combination of atomic orbital wave functions
,”
J. Chem. Phys.
90
,
4916
(
1989
).
77.
S.
Lehtola
and
H.
Jónsson
, “
Pipek–Mezey orbital localization using various partial charge estimates
,”
J. Chem. Theory Comput.
10
,
642
(
2014
).
78.
G.
Knizia
, “
Intrinsic atomic orbitals: An unbiased bridge between quantum theory and chemical concepts
,”
J. Chem. Theory Comput.
9
,
4834
(
2013
).
79.
Q.
Sun
, “
Co-iterative augmented Hessian method for orbital optimization
,” arXiv.1610.08423.
80.
A.
Karanovich
,
Y.
Yamamoto
,
K. A.
Jackson
, and
K.
Park
, “
Electronic structure of mononuclear Cu-based molecule from density-functional theory with self-interaction correction
,”
J. Chem. Phys.
155
,
014106
(
2021
).
81.
R. H.
Byrd
,
P.
Lu
,
J.
Nocedal
, and
C.
Zhu
, “
A limited memory algorithm for bound constrained optimization
,”
SIAM J. Sci. Comput.
16
,
1190
(
1995
).
82.
A. K.
Rappe
,
C. J.
Casewit
,
K. S.
Colwell
,
W. A.
Goddard
, and
W. M.
Skiff
, “
UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations
,”
J. Am. Chem. Soc.
114
,
10024
(
1992
).
83.
N. M.
O’Boyle
,
M.
Banck
,
C. A.
James
,
C.
Morley
,
T.
Vandermeersch
, and
G. R.
Hutchison
, “
Open Babel: An open chemical toolbox
,”
J. Cheminf.
3
,
33
(
2011
).
84.
J. P.
Perdew
,
A.
Ruzsinszky
,
G. I.
Csonka
,
O. A.
Vydrov
,
G. E.
Scuseria
,
L. A.
Constantin
,
X.
Zhou
, and
K.
Burke
, “
Restoring the density-gradient expansion for exchange in solids and surfaces
,”
Phys. Rev. Lett.
100
,
136406
(
2008
).
85.
J. W.
Furness
,
A. D.
Kaplan
,
J.
Ning
,
J. P.
Perdew
, and
J.
Sun
, “
Accurate and numerically efficient r2SCAN meta-generalized gradient approximation
,”
J. Phys. Chem. Lett.
11
,
8208
(
2020
).
Published open access through an agreement with Friedrich-Schiller-Universität Jena

Supplementary Material