DFT and ONIOM Simulation of 1,3-Butadiene Polymerization Catalyzed by Neodymium-Based Ziegler–Natta System

Using modern methods of quantum chemistry, a theoretical substantiation of the high cis-stereospecificity of 1,3-butadiene polymerization catalyzed by the neodymium-based Ziegler–Natta system was carried out. For DFT and ONIOM simulation, the most cis-stereospecific active site of the catalytic system was used. By analyzing the total energy, as well as the enthalpy and Gibbs free energy of the simulated catalytically active centers, it was found that the coordination of 1,3-butadiene in the trans-form was more favorable than in the cis-form by 11 kJ/mol. However, as a result of π-allylic insertion mechanism modeling, it was found that the activation energy of cis-1,3-butadiene insertion into the π-allylic neodymium–carbon bond of the terminal group on the reactive growing chain was 10–15 kJ/mol lower than the activation energy of trans-1,3-butadiene insertion. The activation energies did not change when both trans-1,4-butadiene and cis-1,4-butadiene were used for modeling. That is, 1,4-cis-regulation was due not to the primary coordination of 1,3-butadiene in its cis-configuration, but to its lower energy of attachment to the active site. The obtained results allowed us to clarify the mechanism of the high cis-stereospecificity of 1,3-butadiene polymerization by the neodymium-based Ziegler–Natta system.


Introduction
For more than 75 years, Ziegler-Natta catalysts have continued to be the main type of catalytic system used for the industrial production of polymers with a regular structure [1][2][3][4]. Here, it is necessary to note that the design of Ziegler-Natta catalysts was developed for the polymerization of conjugated dienes and the production of synthetic rubbers [5][6][7][8][9]. The large-scale industrial production of stereoregular rubbers for mass use, including cis-1,4-polyisoprene (an analogue of natural rubber) and cis-1,4-polybutadiene, is one of the outstanding achievements of global science and technology [10]. Over the past three decades, technologies for the production of synthetic rubbers using lanthanide, primarily neodymium-containing catalytic systems, have been actively developed [11][12][13][14][15][16]. Despite the large number of works, research in this area continues to be relevant. This is due to the unique capabilities of neodymium-containing catalytic systems: the high activity, stereoregularity, and controlled molecular parameters of the resulting polymers, and the ability to use various monomers in their synthesis [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34]. In industrial practice, the so-called "binary" [35][36][37][38][39] and "ternary" [40][41][42][43] neodymium-containing catalytic systems have received the greatest development. Despite the differences in the nature of the initial components, we agree that the structure of the active centers, under the conditions of a comparable set of atoms, is similar. This is confirmed by the similar characteristics of

Computational Details
Quantum-chemical calculations were carried out using the Orca 5.0 program package [66,67] in the framework of the density functional theory. The hybrid density functional B3LYP [68,69] was used in the calculations in combination with the second-generation splitvalence triple-ζ atomic basis set def2-TZVP presented by Ahlrichs et al. [70,71]. The core electrons of the neodymium atom were described using the def2-ECP pseudopotential [72] recommended for use with def2 basis sets. Since the Nd(III) ion has three unpaired electrons, and the total multiplicity of all the complexes under study was four (confirmed by a series of additional calculations), the spin-unrestricted method was used for all calculations.
A feature of the systems under study was that weak Van der Waals interactions played an important role in them, which were taken into account in the framework of the D3 semiempirical model [73] presented by Grimme et al. All calculations were performed with the full optimization of molecular geometry without any symmetry constrains. Geometry optimization was carried out taking into account the influence of the solvent (hexane) in the self-consistent reaction field (SCRF) conductor-like polarizable continuum model (C-PCM) [74,75].
The calculations of vibrational spectra performed after geometry optimization did not contain imaginary modes. This meant that the structures found corresponded to minima on the total potential energy surface. These calculations were also used to estimate the thermal corrections needed to calculate the total Gibbs free energy of particles (at 298.15 K and 1 atm).
Since it was found in the first part of the work that the results obtained did not allow us to draw unambiguous conclusions (see conclusions in Section 3.1 below), it became clear that to continue the work, larger catalyst systems would need to be considered. Therefore, it was decided to carry out further studies (all results provided in Section 3.2 below) within the framework of the ONIOM two-layer (QM1/QM2) model [76], in which a part was allocated to the system under study, calculated using the high-level quantum chemical method QM1; the system as a whole was calculated using the low-level QM2 method. Then, the results of these two calculations were combined using a special technique.
As a high-level (QM1) method, we used the same approximation as in the first part of the work, B3LYP/def2-TZVP. Based on the results of preliminary calculations, the tightbinding DFT (semi-empirical DFT) method XTB1 presented by Grimme et al. [77,78] was chosen as a low-level QM2 method.
Van der Waals interactions were also taken into account within the semi-empirical D3 model designed by Grimme et al. [73]. All calculations were performed with the full optimization of molecular geometry without any symmetry constrains. Geometry optimization was carried out taking into account the influence of the solvent (hexane) in the continuum model ALPB-a method for the analytical linearization of the Poisson-Boltzmann equation [79].
To calculate the activation energies of the process, we used the standard procedure for searching for the transition state. The calculations of vibrational spectra performed after geometry optimization did not contain imaginary modes for minima and contained one imaginary mode for transition states. Each transition state identified was checked for compliance with the considered chemical process using the internal reaction coordinate (IRC) procedure.
The ONIOM technique involved the selection of a fragment in the system under study, which was calculated within the framework of the high-level QM1 method. This fragment was the same for all systems that were considered in the framework of this study. The active side (AS) fragment chosen for QM1 calculations consisted of 43 atoms (see Figure S1 in Supplementary Materials). Dangling chemical bonds within the framework of the ONIOM model were modeled using hydrogen atoms, which were located along the broken bonds. These hydrogen atoms were completed by the program itself. All other atoms, including the growing polymer chain, were calculated using the low-level QM2 method.

Structural and Thermochemical Analysis of AS
It is known that Nd(III) is characterized by coordination numbers from six to nine. According to numerous studies, Nd(III) is capable of coordinating diene monomers. As a result, the generally accepted mechanism of interaction becomes possible, which includes the stage of the coordination of the monomer at the active site of polymerization. It was shown in [55] that the polymerization of butadiene initiated by the neodymium catalytic system can lead to the formation of six types of active sites (ASs). In this work, the active site AS+C 1 H 2 C 2 HC 3 HC 4 H 2 (Figure 1), which exhibited the highest cis-stereospecificity of action, was chosen as a model structure for calculations.
It is known that the trans-configuration is more favorable for 1,3-butadiene [80]. Our calculations confirmed that the trans-configuration was the most favorable for the existence of 1,3-butadiene both in the gas phase and in hexane. Thus, the difference in Gibbs free energy between the cisand trans-configurations of 1,3-butadiene was 12.2 kJ/mol (∆H 0 298 is 14.5 kJ/mol). Due to the fact that the AS+C 1 H 2 C 2 HC 3 HC 4 H 2 catalytic system considered in this work exhibited high cis-stereospecificity, when calculating the active center, its complexes with both cisand trans-forms of 1,3-butadiene were considered ( Figure 1).
As a result of minimizing the total energy, the optimal structures of the AS+C 1 H 2 C 2 HC 3 HC 4 H 2 catalytic complex were found for both configurations of 1,3-butadiene (hereinafter AS+η-cis-C 4 H 6 and AS+η-trans-C 4 H 6 ) ( Figure 2). The most important geometric characteristics of these structures are presented in Tables 1 and 2. According to Figure 2, the observed free space near Nd(III) made it possible to judge the absence of any steric hindrances for the coordination interaction of Nd(III) with 1,3butadiene in both cisand trans-configurations. In the presented structures, the Nd(III) ion was connected to the isobutyl carbanion and aluminum ions by chloride bridges. In addition, according to Figure 2, the Nd(III) ion in a given ligand environment created coordination structures of distorted octahedral symmetry.
According to Table 1, neodymium, chlorine, and aluminum ions in the common molecular structure lay almost in the same plane. At the same time, the planes of two Al-Cl-Nd-Cl quadrilaterals were almost perpendicular. Scheme 1 shows the possible modes of the coordination interaction of 1,3-butadiene with the active site of polymerization. As a result of minimizing the total energy, the optimal structures of the AS+C 1 H2C 2 HC 3 HC 4 H2 catalytic complex were found for both configurations of 1,3-butadiene (hereinafter AS+η-cis-C4H6 and AS+η-trans-C4H6) ( Figure 2). The most important geometric characteristics of these structures are presented in Table 1 and Table 2.
The orientation of cisand trans-1,3-butadiene relative to i-C 4 H 9 in the active center had to be considered separately, since the further incorporation of the monomer occurred as a result of the interaction of these two particles. Table 1 shows the shortest distances between the terminal atoms of these particles, that is, those atoms that should interact with each other in the future. As can be seen from Table 1, in the case of cis-1,3-butadiene, this distance was~0.4 Å shorter, which indicated a greater possibility for the interaction of this   Table 2. Some geometric characteristics of i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS and i-C 4 H 9 -π-syn-H 2 CHC γ HC β H 2 C α -AS structures optimized at the B3LYP/def2-TZVP level.  According to Figure 2, the observed free space near Nd(III) made it possible to judge the absence of any steric hindrances for the coordination interaction of Nd(III) with 1,3butadiene in both cis-and trans-configurations. In the presented structures, the Nd(III) ion was connected to the isobutyl carbanion and aluminum ions by chloride bridges. In addition, according to Figure 2, the Nd(III) ion in a given ligand environment created coordination structures of distorted octahedral symmetry.

i-C 4 H 9 -π-anti-
According to Table 1, neodymium, chlorine, and aluminum ions in the common molecular structure lay almost in the same plane. At the same time, the planes of two Al-Cl-Nd-Cl quadrilaterals were almost perpendicular. Scheme 1 shows the possible modes of the coordination interaction of 1,3-butadiene with the active site of polymerization. It follows from Table 1 that the distance between Nd(III) and the carbon atoms in the 1,3-butadiene in both the cis-and trans-forms was ~3 Å. This indicated that the main role in the coordination bonding of Nd(III) with 1,3-butadiene was played by π-bonds, as well as electrostatic (dispersion and Van der Waals) interactions. The η 4 -cis/trans-coordination Scheme 1. Structurally different modes of butadiene coordination at the metal M in butadiene π-complexes [56]. An analysis of the total energy, as well as the enthalpy and Gibbs free energy, of the AS+η-cis-C 4 H 6 and AS+η-trans-C 4 H 6 complexes showed that the coordination of trans-1,3butadiene was more favorable than the coordination binding of cis-1,3-butadiene by 11.0 kJ/mol. The changes in enthalpy and Gibbs free energy under standard conditions for transcis transformation were 13.7 and 7.8 kJ/mol, respectively. The fact that these values were somewhat lower (in particular, the change in the Gibbs free energy was lower by 5 kJ/mol) in comparison with 1,3-butadiene, which was outside the zone of the polymerization process, was not a sufficient justification for the formation of cis-1,4-polybutadiene.
The result of the next stage, the polymerization initiation reaction involving AS and 1,3-butadiene, was the formation of the i-C 4 H 9 -H 2 CHC γ HC β H 2 C α -AS structure. The possible variations in the terminal group structures of the reactive growing chain are presented in Scheme 2: Since the π-allylic neodymium-carbon bond was an integral part of the terminal group on the reactive growing chain, the i-C4H9-π-anti-H2CHC γ HC β H2C α -AS and i-C4H9π-syn-H2CHC γ HC β H2C α -AS complexes were used for further calculations.
As can be seen from Figure 3, there were noticeable differences for the terminal allylic end unit, which existed as anti-and syn-π-allylic complexes. First of all, in the case of the anti-π-allylic complexes, the distance from three atoms of trans-1,3-butadiene to Nd(III) practically coincided. This indicated the predominance of π-and dispersion interactions when this form was bound to Nd(III). In the case of the syn-π-allylic complexes, monodentate metal-ligand bonds and π-interactions with the double bond of the terminal unit of the polymer chain were formed. Since the π-allylic neodymium-carbon bond was an integral part of the terminal group on the reactive growing chain, the i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS and i-C 4 H 9 -π-syn-H 2 CHC γ HC β H 2 C α -AS complexes were used for further calculations.
As can be seen from Figure 3, there were noticeable differences for the terminal allylic end unit, which existed as antiand syn-π-allylic complexes. First of all, in the case of the anti-π-allylic complexes, the distance from three atoms of trans-1,3-butadiene to Nd(III) practically coincided. This indicated the predominance of πand dispersion interactions when this form was bound to Nd(III). In the case of the syn-π-allylic complexes, monodentate metal-ligand bonds and π-interactions with the double bond of the terminal unit of the polymer chain were formed.
As a result, the ligand environment of the Nd(III) ion in the terminal group on the reactive growing chain underwent changes when the π-allylic structure was repositioned from the anti-π-allylic complex to the syn-π-allylic complex. The ligand environment of Nd(III) in the active syn-π-allylic complex did not change much compared to the AS+trans-C 4 H 9 complex described above. In the case of the anti-π-allylic complex, the chlorine atoms were located in the same plane and formed a rectangle 3.5 × 4 Å in size. The Nd(III) ion was located above the center of this rectangle, and above it was a semiring of cis-1,3-butadiene.
The Gibbs free energy of the formation of i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS and i-C 4 H 9 -π-syn-H 2 CHC γ HC β H 2 C α -AS complexes from the more favorable and, accordingly, more probable complex AS+trans-C 4 H 6 was −84.5 and −84.9 kJ/mol, respectively. That is, we could conclude that the formation of both forms of the complex was almost equally likely, but this fact also could not explain the predominant formation of cis-1,4-polybutadiene.
Next, we studied the complexes simulating the coordination of 1,3-butadiene in the terminal group on the reactive growing chain. Based on the assumption of the cisstereospecificity of the AS from the whole variety of structures, the complexes that lead to the formation of cis-1,4-polybutadiene were considered: i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS+η-cis-C 4 H 6 (A) and i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS+η-trans-C 4 H 6 (B) (Figure 3). The coordination of 1,3-butadiene from the volume of the reaction mixture was accompanied by the displacement of the growing polymer chain relative to the active center and the freeing of space for the monomer. According to classical concepts (Scheme 2), at this moment, a σ-cis-allylic (σ-C1-cis) structure in the terminal group should form, in which a covalent neodymium-carbon bond of the reactive growing chain is formed. However, our calculations showed that at this stage, too, the π-allylic structure of the anti-configuration was energetically more favorable than the σ-C1-cis structure. At the same time, the coordination of cis-1,3-butadiene by the Nd(III) ion mainly occurred due to π-interaction with the terminal CH-CH 2 group of 1,3-butadiene (the smallest distance to Nd(III) was presented by the C 1 and C 2 atoms of 1,3-butadiene), while the coordination of trans-1,3-butadiene occurred due to the π-interaction with the central CH-CH group (the smallest distance to Nd(III) was presented by the C 2 and C 3 atoms of 1,3-butadiene). It should also be noted that the trans-form of 1,3-butadiene remained practically planar, while in the cis-form some distortion was observed, and the terminal CH 2 group deviated from the plane formed by the other three carbon atoms by about 15 • . As for the coordination of the polymer fragment relative to AS, it was almost identical in these two cases (Table 3). Table 3. Some geometric characteristics of i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS+η-cis-C 4 H 6 and i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS+η-trans-C 4 H 6 structures optimized at the B3LYP/def2-TZVP level.  It should be noted that the coordination of 1,3-butadiene in the cisand trans-forms significantly changed the geometry of the active site. If in i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS all four chlorine atoms lay almost in the same plane (Figure 3a), then when 1,3-butadiene was coordinated, the Cl 1 -Nd-Cl 2 plane was almost perpendicular to the Cl 3 -Nd-Cl 4 plane (the Cl 1 and Cl 2 atoms were coordinated with one aluminum atom, and Cl 3 and Cl 4 with the second). Due to the fact that both the terminal π-allylic structure and 1,3-butadiene in the cisand trans-forms were linked to Nd(III) via π-bonds, it was not possible to reliably determine the coordination number of Nd(III).

i-C 4 H 9 -π-anti-
The thermochemical parameters for the formation of i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS+η-cis-C 4 H 6 and i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS+η-trans-C 4 H 6 complexes ( Figure 4) were evaluated by the reaction: i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS+trans-C 4 H 6 i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS+η-trans/cis-C 4 H 6 (1) As can be seen from Figure 3, there were noticeable differences for the terminal allylic end unit, which existed as anti-and syn-π-allylic complexes. First of all, in the case of the anti-π-allylic complexes, the distance from three atoms of trans-1,3-butadiene to Nd(III) practically coincided. This indicated the predominance of π-and dispersion interactions when this form was bound to Nd(III). In the case of the syn-π-allylic complexes, monodentate metal-ligand bonds and π-interactions with the double bond of the terminal unit of the polymer chain were formed. As a result, the ligand environment of the Nd(III) ion in the terminal group on the reactive growing chain underwent changes when the π-allylic structure was repositioned from the anti-π-allylic complex to the syn-π-allylic complex. The ligand environment of Nd(III) in the active syn-π-allylic complex did not change much compared to the AS+trans-C4H9 complex described above. In the case of the anti-π-allylic complex, the chlorine atoms were located in the same plane and formed a rectangle 3.5 × 4 Å in size. The Nd(III) ion was located above the center of this rectangle, and above it was a semiring of cis-1,3-butadiene. i-C4H9-π-anti-H2CHC γ HC β H2C α -AS+η-cis-C4H6 i-C4H9-π-anti-H2CHC γ HC β H2C α -AS+η-trans-C4H6 Figure 4. Structures of i-C4H9-π-anti-H2CHC γ HC β H2C α -AS+ η-cis-C4H6 and i-C4H9-π-anti-H2CHC γ HC β H2C α -AS+ η-trans-C4H6 optimized at the B3LYP/def2-TZVP level.

Study of Polymerization Mechanism
At the initial stage, it was decided to calculate the transition state for the trans-cis transition of 1,3-butadiene in the reaction system, taking into account its interaction with the AS. The resulting transition state structures are presented in the Supplementary Ma- Figure 4. Structures of i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS+η-cis-C 4 H 6 and i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS+η-trans-C 4 H 6 optimized at the B3LYP/def2-TZVP level.
The change in the Gibbs free energy under standard conditions for the formation of the i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS+η-cis-C 4 H 6 complex was 11.1 kJ/mol, and for the i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS+η-trans-C 4 H 6 formation it was 18.8 kJ/mol. This meant that under standard conditions, this reaction did not proceed. However, both reactions were exothermic: the ∆H 0 298 value in the first case was -46.5 kJ/mol, and in the second -43.2 kJ/mol. Accordingly, the ∆S 0 298 value for i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS+η-cis-C 4 H 6 was -193.3 J/K·mol, and for i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS+η-trans-C 4 H 6 it was -208.1 J/K·mol. Due to the fact that the ∆H and ∆S had the same signs, these reactions were reversible, and an increase in the concentration of reagents shifted the equilibrium towards the formation of products. Taking into account that the ∆G 0 298 in both cases had a rather small positive value, under the conditions of a real reaction mixture, where the concentration of 1,3-butadiene is significantly higher than the concentration of active sites, both reactions should proceed in the forward direction. It should also be noted that the Gibbs free energy of the trans-cis transition of 1.3-butadiene upon interaction with the active center remained equal to 7.7 kJ/mol, as in Figure 2. The values of ∆G 0 298 for these reactions were quite close. This meant that approximately equal concentrations of both the antiand syn-forms of terminal π-allyl structures should be expected in the reaction system.
Thus, the analysis of the structural characteristics and thermochemical features of the formation of AS complexes with 1,3-butadiene, corresponding to stages such as monomer coordination and the initiation and growth of the polymer chain, was insufficient to justify the predominant formation of cis-1,4-polybutadiene by the neodymium-based Ziegler-Natta system. This required a more detailed study of the process mechanism.

Study of Polymerization Mechanism
At the initial stage, it was decided to calculate the transition state for the trans-cis transition of 1,3-butadiene in the reaction system, taking into account its interaction with the AS. The resulting transition state structures are presented in the Supplementary Materials (Tables S7-S57). Further, within the framework of this section, we consider only the change of the Gibbs free energy; however, the Supplementary Materials also provide tables similar to Table 4 but containing the values of the change in the total energy E 0 (Table S58) and enthalpy (Table S59). The activation energy (∆G = ) determined for the trans-cis transition for 1,3-butadiene in the free state, associated with the rotation of the CH 2 group around the C=C bond, was 27.0 kJ/mol. A similar transition was calculated for 1,3-butadiene interacting with the active site, and the activation energy in this case increased to 43 kJ/mol. That is, the assumption that the active center facilitated the trans-cis transition for 1,3butadiene was not confirmed. Thus, the initial coordination of 1,3-butadiene at the active center in the trans-form was legitimate for subsequent calculations.
Next, the transition states associated with the insertion of 1,3-butadiene into the polymer chain were simulated. At the stage of initiation, this was ensured by the rotation of 1,3-butadiene in such a way that the terminal CH 2 group was oriented toward the i-C 4 H 9 group bound to Nd(III), as well as with some distance between i-C 4 H 9 and Nd(III) and approaching the terminal CH 2 group of 1,3-butadiene. At the stage of polymer chain growth, everything should have proceeded in exactly the same way, only the position of i-C 4 H 9 was occupied by a reactive growing chain. Transition states were successfully identified for transand cis-1,3-butadiene. The correspondence of the transition states of the studied reaction was confirmed by the calculation of "descents" to the reactants and products in the IRC procedure. The calculated activation energy (∆G = ) of the polymerization process for trans-1,3-butadiene was 67.0 kJ/mol, and for cis-1,3-butadiene it was 54.3 kJ/mol. It should be noted that the activation energy of depolymerization was twice the activation energy of the direct reaction and was 122.2 kJ/mol for trans-1,3-butadiene and 112.8 kJ/mol for cis-1,3-butadiene. This indicated an extremely low probability of this process occurring. If we assume that the initial form in this reaction was the active center with coordinated trans-1,3-butadiene (i-C 4 H 9 -AS-trans-C 4 H 6 ), then the insertion of the first trans-unit into the reactive growing chain proceeded in one stage, the activation energy of which was the activation energy of the insertion of 1,3-butadiene in the trans-configuration into the terminal group on the reactive growing chain, that is, ∆G = trans = 67.0 kJ/mol. Since the insertion of 1,3-butadiene in the cis-configuration occurred in two stages, the activation energy had to be measured from the structure with the lowest energy ( Figure 5), and in this case it was ∆G = cis = 54.3 + 6.6 = 60.9 kJ/mol. Thus, it was shown that at the stage of initiation, the formation of the terminal πallylic structure of the anti-configuration proceeded with a slightly lower activation energy compared to the formation of the syn-π-allylic complexes. At the same time, the energy effect of the insertion of trans-1,3-butadiene into the active center from the volume of the reaction mixture for both forms was almost the same, amounting to ∆G cis = 51.7 kJ/mol and ∆G trans = 55.1 kJ/mol.
At the next stage of the reactive polymer chain growth, there were two almost thermodynamically equiprobable structures, i-C 4 H 9 -π-anti-H 2 CHC γ HC β H 2 C α -AS+η-trans-C 4 H 6 and i-C 4 H 9 -π-syn-H 2 CHC γ HC β H 2 C α -AS+η-trans-C 4 H 6 , with somewhat lower energy activation during the formation of the first structure. Further, the formation of four different structures was possible, which for brevity we denote as i-C 4 H 9 -CC-AS, i-C 4 H 9 -CT-AS, i-C 4 H 9 -TC-AS, and i-C 4 H 9 -TT-AS (where "T" indicates trans-C 4 H 6 or π-syn-H 2 CHC γ HC β H 2 C α fragments and "C" indicates cis-C 4 H 6 or π-anti-H 2 CHC γ HC β H 2 C α fragments). To calculate the activation energy of i-C 4 H 9 -TC-AS and i-C 4 H 9 -TT-AS formation, we considered i-C 4 H 9 -C-AS-T and i-C 4 H 9 -T-AS-T, respectively. However, to form the other two structures, it was first necessary to calculate the activation energy of the trans-cis transition of 1,3-butadiene. The activation energy (∆G = ) of the transition i-C 4 H 9 -C-AS-T i-C 4 H 9 -C-AS-C (2) practically did not change in comparison with the initiating stage and was equal to 41.0 kJ/mol. The activation energy (∆G = ) of the transition was also almost identical and equal to 41.9 kJ/mol. That is, this activation energy was constant and was the same in all cases.
Polymers 2023, 15, x FOR PEER REVIEW 11 of 19 studied reaction was confirmed by the calculation of "descents" to the reactants and products in the IRC procedure. The calculated activation energy (ΔG ≠ ) of the polymerization process for trans-1,3-butadiene was 67.0 kJ/mol, and for cis-1,3-butadiene it was 54.3 kJ/mol. It should be noted that the activation energy of depolymerization was twice the activation energy of the direct reaction and was 122.2 kJ/mol for trans-1,3-butadiene and 112.8 kJ/mol for cis-1,3-butadiene. This indicated an extremely low probability of this process occurring. If we assume that the initial form in this reaction was the active center with coordinated trans-1,3-butadiene (i-C4H9-AS-trans-C4H6), then the insertion of the first trans-unit into the reactive growing chain proceeded in one stage, the activation energy of which was the activation energy of the insertion of 1,3-butadiene in the trans-configuration into the terminal group on the reactive growing chain, that is, ΔG ≠ trans = 67.0 kJ/mol. Since the insertion of 1,3-butadiene in the cis-configuration occurred in two stages, the activation energy had to be measured from the structure with the lowest energy ( Figure 5), and in this case it was ΔG ≠ cis = 54.3 + 6.6 = 60.9 kJ/mol.  Next, the transition states of the insertion of the second 1,3-butadiene molecule into the reactive polymer chain were calculated in all four reactions. Generalized detailed information regarding the change in activation energy and Gibbs free energy during this process is presented in Table 4.
From the analysis of the data in Table 4, it can be seen that the activation energy increased in all cases when the second unit was inserted, and when a cis-unit was inserted into the reactive growing polymer chain, it remained lower than when a trans-unit was inserted. The lowest activation energy was observed in the case of the formation of pure cis-1,4-polybutadiene. In this case, the difference in activation energies compared to the insertion of a trans-unit already exceeded 10 kJ/mol. Also interesting is the fact that the insertion of a cis-unit into a trans-unit still presented a lower activation energy, although the difference was not large. The least likely was the insertion of a trans-unit into a cis-unit.
Next, we consider the reaction of the insertion of the third and subsequent links into the reactive growing polymer chain. Here, we assumed that trans-1,3-butadiene was also inserted into the structures obtained at the previous stage from the bulk of the reaction mixture. However, a serious problem then arose, related to which combinations of cisand trans-segments of the polymer should be taken into account in the calculations. Even if at the second stage it was still not difficult to determine this, since there were only four options, further on the number of possible combinations increased exponentially. Therefore, it was decided to consider only two options-the insertion of a cisand trans-1,4-butadiene unit into the most probable polymer structure obtained at the previous stage, and the formation of trans-1,4-polybutadiene for comparison, as practice shows that this is possible with a slow process. It was also decided to carry out calculations before the insertion of five 1,3-butadiene molecules into the reactive growing polymer chain. Since, in this case, for cis-1,4-polybutadiene the number of monomer units exceeded the size of the Kuhn segment, which was 2-3 units, it became possible to extrapolate the results obtained to real objects, the polymer chain length of which reached several thousand monomer units.
At the third stage, taking into account the above, we considered (trans-C 4 H 6 ) 2 -AS-trans-C 4 H 6 and (cis-C 4 H 6 ) 2 -AS-trans-C 4 H 6 as initial structures, and as reaction products (trans-C 4 H 6 ) 3 -AS, (cis-C 4 H 6 ) 3 -AS and trans-C 4 H 6 (cis-C 4 H 6 ) 2 -AS were obtained. The formation of (cis-C 4 H 6 ) 3 -AS again required a preliminary trans-cis transition step in butadiene according to the reaction: (cis-C 4 H 6 ) 2 -AS-trans-C 4 H 6 (cis-C 4 H 6 ) 2 -AS-cis-C 4 H 6 (4) The calculated activation energy of this reaction was significantly lower than the reactions of polymerization initiation and the beginning stage of the reactive polymer chain growth, at 24.6 kJ/mol. This value turned out to be even lower than the activation energy of the trans-cis transition of 1,3-butadiene outside the zone of the polymerization process. The change in the Gibbs free energy in the polymerization reaction zone according to Equation (4) was 14.6 kJ/mol, which was somewhat higher than the value obtained outside the polymerization process zone. Changes in the activation energy and Gibbs free energy were calculated in a similar way for the next stages of the 1,3-butadiene polymerization process.
Thus, the analysis of the data presented in Table 4 allowed us to conclude that the activation energy of the trans-cis transition of 1,3-butadiene interacting with the active center practically stabilized at a value approximately 3 kJ/mol lower compared to the activation energy for 1,3-butadiene, which was outside the zone of the polymerization process. It is expected that this value will be maintained in the future. The activation energy of the insertion of cis-1,3-butadiene into the reactive growing chain was 10-15 kJ/mol lower than the activation energy of the insertion of trans-1,3butadiene. The activation energy values did not change when both transand cis-1,4-polybutadiene terminal units were used for calculations. The almost coinciding activation energy values for the insertion of cis-1,3butadiene into the reactive growing chain at the last two stages indicated that the value of 61 kJ/mol could be taken as the activation energy of the formation of cis-1,4-polybutadiene at the AS. The obtained value agreed with the experimental data on the activation energy of the polymerization of dienes on lanthanide catalysts equal 30-60 kJ/mol [10,31]. If there are some macroscopic factors that cannot be taken into account in the framework of a quantum chemical study, then the activation energy of the insertion of cis-1,3-butadiene into the reactive growing chain will decrease by another~10 kJ/mol and will be 51 kJ/mol, while the activation energy of the insertion of trans-1,3-butadiene will not change. Thus, in this case, the difference in activation energies will already be about 20-25 kJ/mol.

Conclusions
According to the structural and thermochemical study of the polymerization of 1,3butadiene initiated by the most cis-stereospecific active site of the NdCl 3 ·3TBF-Al(i-C 4 H 9 ) 3 catalytic system, the coordination of 1,3-butadiene by Nd(III) was equally probable in its cisand trans-configurations. An analysis of the total energy, as well as the enthalpy and Gibbs free energy, of the AS+η-cis-C 4 H 6 and AS+η-trans-C 4 H 6 complexes showed that the coordination of 1,3-butadiene in the trans-form was more favorable than in the cis-form by 11 kJ/mol.
The transition states associated with the insertion of 1,3-butadiene into the reactive growing polymer chain were simulated. It was found that the energy of the insertion of cis-1,3-butadiene into the reactive growing polymer chain was 10-15 kJ/mol lower than the activation energy of the insertion of trans-1,3-butadiene. The activation energy values did not change when both transand cis-1,4-polybutadiene terminal units were used for calculations. The almost coinciding values of the activation energy for the insertion of cis-1,3-butadiene into the reactive growing polymer chain at the last two stages indicated that the value of 61 kJ/mol could be taken as the activation energy of the formation of cis-1,4-polybutadiene at the AS.
Based on the results of the study, the following mechanism of 1,3-butadiene polymerization by the neodymium-based Ziegler-Natta system could be proposed: the coordination of trans-1,3-butadiene at the active site; the trans-cis isomerization of 1,3-butadiene; and the insertion of cis-1,3-butadiene into the reactive growing polymer chain.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.