Atomistic Molecular Dynamics Simulations of the Initial Crystallization Stage in an SWCNT-Polyetherimide Nanocomposite

Crystallization of all-aromatic heterocyclic polymers typically results in an improvement of their thermo-mechanical properties. Nucleation agents may be used to promote crystallization, and it is well known that the incorporation of nanoparticles, and in particular carbon-based nanofillers, may induce or accelerate crystallization through nucleation. The present study addresses the structural properties of polyetherimide-based nanocomposites and the initial stages of polyetherimide crystallization as a result of single-walled carbon nanotube (SWCNT) incorporation. We selected two amorphous thermoplastic polyetherimides ODPA-P3 and aBPDA-P3 based on 3,3′,4,4′-oxydiphthalic dianhydride (ODPA), 2,3′,3,4′-biphenyltetracarboxylic dianhydride (aBPDA) and diamine 1,4-bis[4-(4-aminophenoxy)phenoxy]benzene (P3) and simulated the onset of crystallization in the presence of SWCNTs using atomistic molecular dynamics. For ODPA-P3, we found that the planar phthalimide and phenylene moieties show pronounced ordering near the CNT (carbon nanotube) surface, which can be regarded as the initial stage of crystallization. We will discuss two possible mechanisms for ODPA-P3 crystallization in the presence of SWCNTs: the spatial confinement caused by the CNTs and π–π interactions at the CNT-polymer matrix interface. Based on our simulation results, we propose that ODPA-P3 crystallization is most likely initiated by favorable π–π interactions between the carbon nanofiller surface and the planar ODPA-P3 phthalimide and phenylene moieties.

The nanofiller surface may have an impact on the structure of neighboring polymer chains [26]. When the nanofiller size is comparable to the size of a polymer chain, the nanofiller particles may create a geometric confinement, which is known to result in the orientation of adjacent polymer chains [16]. Furthermore, the nanofiller surface may confine and slow down the segmental relaxation of polymers [27,28]. A significant structural ordering near the nanofiller surface may result via specific interactions at the polymer-nanofiller interface, especially when the size of the nanofiller particle exceeds that of the polymer chains, similar to the size-dependent soft epitaxy [16,26]. The structural properties of the polymer near the nanofiller are also determined by the chemical composition of the participating components [10,[29][30][31][32][33][34]. For aliphatic polymers, researchers have shown that there are specific interactions between CH-groups and CNT atoms [35]. In the case of conjugated polymers, there is evidence that due to π-π interactions, chain fragments may orient along the nanofiller surface [36,37].
However, questions of why crystallization of the polymer matrix takes place near the nanofiller surface and what is the role of geometric confinement or other specific interactions (for example, π-π) remain largely unanswered, and this is in particular true for aromatic heterocyclic polymers [38].
A detailed investigation of the nanocomposite structure at the polymer-nanofiller interface allows for a more in-depth understanding of the mechanisms of crystallization initiation of heterocyclic polymers as a result of carbon nanofiller incorporation. This is essential for the development of new nanocomposite materials with improved mechanical properties. Atomistically detailed computer-aided simulations make it possible to study the structure of the polymer-nanofiller interface on the scale of single atoms [39][40][41][42] and, therefore, help to provide insights into the structural ordering near the carbon nanofiller surface and why it improves the nanocomposite mechanical properties.
Simulation studies of the structural properties of polymers near the nanofiller surface have been performed using coarse-grained [43] and atomistic [44] models. Rather simple polymers in terms of chemical structure have been simulated so far, such as polyethylene [43,45] or poly(trimethylene terephthalate) [44], and the effects of single-walled nanotubes (SWCNTs) and graphene on the structural properties of alkane polymer melts have been considered [46]. Researchers have demonstrated that the structural ordering of polymer chains near the CNT surface may be the result of geometric confinements due to the presence of the nanofiller. The graphene surface may have a different impact; according to the results of Yang et al. [46], the polymer-graphene interface structure is determined by specific molecular interactions between the participating components. Simulations have been performed to study the structural ordering of conjugated polymers near the carbon nanofiller surface. Asadinezhad et al. have reported that the planar aromatic rings of poly(trimethylene terephthalate) are oriented parallel with the CNT axis, which may be attributed to π-π interactions [44]. More recent research has focused on fully-aromatic high-performance polymers with high temperature stability, as for example, thermoplastic polyimides [7,47,48]. Typical commercially available thermoplastic PEIs do not crystallize in the presence of nanofillers. However, at the Institute of Macromolecular Compounds of the Russian Academy of Sciences, a crystallizable thermoplastic polyetherimide R-BAPB based on 1,3-bis-(3 ,4-dicarboxyphenoxy)-benzene (dianhydride R) and 4,4 -bis-(4"-aminophenoxy)-diphenyl (diamine BAPB), was synthesized [49], which was used as a polymer binder to produce nanocomposites by melt processing techniques. In our previous studies, we used atomistic molecular dynamics (MD) simulations to investigate the structural ordering of planar moieties of the R-BAPB chain near the surface of SWCNTs [50,51] and graphene [52,53]. It is important to note that the structural ordering of the planar moieties of R-BAPB leads to an improvement in mechanical properties at the polymer-nanofiller interface. Our simulations [52] showed that the average values of the elastic modulus and yield stress of crystallized R-BAPB near the nanofiller surface increased in comparison to R-BAPB in the amorphous state. We were able to clearly demonstrate that the structural ordering of the R-BAPB planar moieties near the nanofiller surface resulted in an improvement in the mechanical properties of the polymer-CNT nanocomposite [52]. The model developed by us previously made it possible to qualitatively describe the initial stages of R-BAPB crystallization near the CNT surface, which was confirmed by experimental results [24,25]. Based on our previous work, we want to stress that special care should be taken to choose the force field, in particular for the parameterization of electrostatic interactions in the case of heterocyclic polyetherimides with polar groups [50][51][52][53][54][55][56][57][58][59][60][61][62].
Quite new is the fact that the incorporation of SWCNTs may nucleate crystallization even in amorphous polyetherimides. In fact, the structural, thermal and mechanical properties of two amorphous thermoplastic polyetherimides ODPA-P3 and aBPDA-P3 based on 3,3 ,4,4 -oxydiphthalic dianhydride (ODPA), 2,3 ,3,4 -biphenyltetracarboxylic dianhydride (aBPDA) and diamine 1,4-bis[4-(4-aminophenoxy)phenoxy]benzene (P3) were investigated [63,64], and the authors showed that incorporation of SWCNTs resulted in the crystallization of the ODPA-P3 matrix, as opposed to aBPDA-P3, which remained amorphous after the addition of the SWCNTs. The mechanism of SWCNT-induced crystallization of ODPA-P3 remains unclear; in fact, it is not clear whether, for example, crystallization is initiated due to spatial confinement, because of the carbon nanofiller presence, or due to π-π interactions between the aromatic cycles of the CNTs and the PEI phenylene and/or heterocyclic (phthalimide) moieties.
To understand better how SWCNTs nucleate crystallization in amorphous thermoplastic polyetherimides, we studied SWCNT-filled nanocomposites based on ODPA-P3 and aBPDA-P3. The backbone of both polymers is almost identical with the exception of the dianhydrides used. 3,3′,4,4′oxydiphthalic dianhydride (ODPA) induces a slight kink in the polymer backbone, whereas 2,3′,3,4′biphenyltetracarboxylic dianhydride (aBPDA) induces a local 90° kink in the polymer backbone, which severely compromises the backbone linearity; see Figure 1. Chemical structure of the thermoplastic polyetherimide (PEI) repeat units: 3,3′,4,4′-oxydiphthalic dianhydride (ODPA)-P3 and 2,3′,3,4′-biphenyltetracarboxylic dianhydride (aBPDA)-P3. Arrows mark the vectors D and P aligned along phenylene and phthalimide planar moieties of the polyetherimide chains under study, for which the orientation to the nanotube axis was investigated in this paper. Phenylene rings in the moiety designated by the D vector can rotate out of a planar structure; however, such deviations from a planar structure are small, especially in the vicinity of the CNT surface, and phenylene rings are mostly coplanar.
In order to study the intermolecular structure at the polymer-SWCNT interface, we analyzed the orientation-related ordering of the PEI chain planar moieties shown by the vectors P or D in Figure 1. It should be noted that the phenylene rings designated by the vector D could form non-planar conformations due to the rotation of the rings. However, such deviations are rather small especially in the vicinity of the CNT surface, and both rings are mostly in a coplanar conformation. Chemical structure of the thermoplastic polyetherimide (PEI) repeat units: 3,3 ,4,4 -oxydiphthalic dianhydride (ODPA)-P3 and 2,3 ,3,4 -biphenyltetracarboxylic dianhydride (aBPDA)-P3. Arrows mark the vectors D and P aligned along phenylene and phthalimide planar moieties of the polyetherimide chains under study, for which the orientation to the nanotube axis was investigated in this paper. Phenylene rings in the moiety designated by the D vector can rotate out of a planar structure; however, such deviations from a planar structure are small, especially in the vicinity of the CNT surface, and phenylene rings are mostly coplanar.
In order to study the intermolecular structure at the polymer-SWCNT interface, we analyzed the orientation-related ordering of the PEI chain planar moieties shown by the vectors P or D in Figure 1. It should be noted that the phenylene rings designated by the vector D could form non-planar conformations due to the rotation of the rings. However, such deviations are rather small especially in the vicinity of the CNT surface, and both rings are mostly in a coplanar conformation.
Previously, we used atomistic molecular dynamics to explore the structural properties of nanocomposites based on amorphous R-BAPS containing polar sulfone groups (-SO 2 -) and partly crystallizable R-BAPB [50]. We showed that composites based on R-BAPB build up an ordered structure near the filler surface due to the orientation of the polymer chain planar fragments primarily along the nanoparticle surface [50,51,53]. At the same time, polyetherimide R-BAPS showed no orientation-related ordering in the corresponding composites. The obtained results correlate with the published experimental data on nanocomposites based on these PEIs [24,25]. However, the presence of a polar sulfone group in the repeat unit of the polyetherimide R-BAPS does not give a definitive answer to the question whether the lack of crystallization of this polyetherimide is related to the particular non-conformation behavior of the chains of this polymer or to specific dipole-dipole interactions between the polar groups. Therefore, in the present study, we simulated the structural properties of two PEIs, ODPA-P3 and aBPDA-P3, with a simpler molecular backbone composition, i.e., only phenyl rings, ether linkages and phthalimide units are present. The similarity of their chemical structure and the absence of the sulfone groups will enable us to provide modeling insights into the initial PEI SWCNT-induced crystallization.
The main goal of the present study is to elucidate the underlying molecular mechanism responsible for the initial stages of SWCNT-ODPA-P3 crystallization based on microsecond simulations. In studies of polymer structural properties, typical simulation times do not exceed 100 ns; however, this is not enough for heterocyclic polymers, which require much longer simulation times [55]. To validate the simulation approach and the chosen force field, we started with reproducing the known experimental differences between the structural properties of ODPA-P3 and aBPDA-P3 in the presence of SWCNTs. Then, with the help of excluded volume interaction parameters, we compared the influence of spatial restrictions, due to the nanotube presence and π-π stacking with the initial stage of SWCNT-ODPA-P3 crystallization.

Models and Simulation Method
The initial configurations of ODPA-P3 and aBPDA-P3 were created using an approach previously developed by us [50,51]. To start, 27 polymer chains with a polymerization degree of N p = 8 were randomly placed into a periodic cell. After that, a SWCNT with a chirality of (5,5), a diameter of 0.7 nm, and a length of 4.7 nm comprised of 400 carbon atoms was placed into the center of the cell. In order to maintain the proper carbon atom valency, hydrogen atoms were attached at the carbon atoms positioned at both carbon nanotube termini. No partial charges were assigned to atoms in the CNT. The CNT model is similar to the one used in our previous studies [50,51].
The initially-generated nanocomposite was compressed for 26 ns. To relieve residual stresses after the compression step, an annealing process was performed, consisting of three cycles of temperature drop and increase, within the range of 800 to 300 K. After the annealing step, the systems were instantly cooled down from T = 800 K to 600 K followed by a simulation for 2 µs at a temperature exceeding the glass transition point (T g ) of the polyetherimides. This simulation timescale is close to the slowest relaxation processes related to moving a polymer coil over distances comparable to its own size, as was shown for the R-BAPB and R-BAPS aromatic polyetherimides with polymerization degree N p = 8 at temperature 600 K [55,62]. During equilibration, the average sizes of polymer chains reached their equilibrium values, which are in good agreement with the theoretical estimates; see Figure S1 in the Supplementary Materials. To investigate the polymer structure at the polymer-CNT interface, the 1 µs-long production runs were performed at T = 600 K and T = 580 K. The choice of said temperatures was motivated by previous experimental results [63]. In fact, experimental results reported by Hegde et al. [63] showed that the melting point of the nanocomposite based on ODPA-P3 with the minimum SWCNT loading (≈0.1 vol %) is about 580 K. Therefore, the simulations were performed in the temperature range close to the experimental melting point.

Results and Discussion
At the initial simulation stage, as in our previous study [50], a test was carried out to estimate the influence of the partial charges on the structural properties of thermoplastic PEIs. Assigning partial charges substantially increases the demand on computational resources [50,55,58]. Thus, in our study, we performed only 100 ns-long simulations with partial charges after a 1.5 µs-long equilibration procedure.
The values of partial charges were calculated using the Hartree-Fock quantum chemical method in the basic set of wave functions 6-31G* [58]. The electrostatic interactions during MD simulations were calculated using the Particle-mesh Ewald (PME) scheme [68] with parameters similar to our previous study [54].
As for the R-BAPB [50], the results obtained confirmed only a weak impact of the partial charges on the polymer layer structure at the interface; see Figure S3 in the Supplementary Materials. This may be attributed to the fact that the values of partial charges on the carbon nanotube atoms equal zero, and the explicit partial charges in the polymer have practically no influence on the general pattern of the structural ordering of the PEI planar fragments near the CNT surface. Therefore, it was decided to ignore the polymer partial charges, which allowed for a significant reduction in simulation time.
To investigate the PEI structure, the pair distribution functions were calculated for different groups of atoms where N A and N B are numbers of type A and B atoms, respectively, ρ B is the mean density of type B atoms located inside the sphere of radius r from the atom A, r ij is the distance between two atoms, and δ is the Dirac delta function. The distant-dependent density distribution of the PEIs monomers has been calculated: where m(r) is the polymer mass in the cylinder layer ∆r at a distance r from the CNT axis. The cylinder layer height is determined by the CNT length, h CNT ≈ 4.7 nm. The pair distribution functions between all carbon atoms in the polyetherimide aromatic rings and carbon nanotube g ar−CNT (r), all carbon atoms in the aromatic rings g ar−ar (r), as well as the dependence of the PEI density on the distance from the CNT axis ρ(r) have been calculated; see Figure 2.
The pair distribution functions g ar−CNT (r) and g ar−ar (r) of the two PEIs do not differ much from each other; see Figure 2a,b. The pair distribution functions g ar−CNT (r), as usual [50,69,70], demonstrate two shoulders responsible for interactions between the nanotube and the polymer chains [50]. The vertices of the curves for both polyetherimides are located at the same distance; see Figure 2a. The pair distribution functions g ar−ar (r) (Figure 2b) show a similar trend for both PEIs.
Two local peaks on g ar−ar (r) at r ≈ 0.45 nm and r ≈ 0.58 nm may be related to interactions between the carbon atoms of the heterocyclic rings of the PEI and the CNT [57]. The dependences of ρ(r) near the CNT surface for the two PEIs are quite different; see Figure 2c. The density of the ODPA-P3 chains at a distance of the first peak (r = 0.75 nm) near the CNT surface is much higher than the density of the aBPDA-P3 chains. Besides, the dependence ρ(r) of ODPA-P3 has a more pronounced second maximum, which shows that the ODPA-P3 chains adsorb more strongly onto the CNT surface than the aBPDA-P3 chains. Two local peaks on ( ) − ar ar g r at r ≈ 0.45 nm and r ≈ 0.58 nm may be related to interactions between the carbon atoms of the heterocyclic rings of the PEI and the CNT [57]. The dependences of ( ) ρ r near the CNT surface for the two PEIs are quite different; see Figure 2c. The density of the ODPA-P3 chains at a distance of the first peak (r = 0.75 nm) near the CNT surface is much higher than the density of the aBPDA-P3 chains. Besides, the dependence ( ) ρ r of ODPA-P3 has a more pronounced second maximum, which shows that the ODPA-P3 chains adsorb more strongly onto the CNT surface than the aBPDA-P3 chains. (c) the dependence of the aBPDA-P3 and ODPA-P3 density on the distance from the CNT axis. T = 600 К in all cases.
In order to study the intermolecular structure at the polymer-SWCNT interface, we analyzed the orientation-related ordering of the phenylene and phthalimide planar moieties in the PEI chains, represented by the D and P vectors in Figure 1 correspondingly, and calculated the order parameters, In order to study the intermolecular structure at the polymer-SWCNT interface, we analyzed the orientation-related ordering of the phenylene and phthalimide planar moieties in the PEI chains, represented by the D and P vectors in Figure 1 correspondingly, and calculated the order parameters, S(r) where r is the distance from the nanotube axis to the polymer chain planar moiety and θ is the angle between the vector P or D directed along the planar moiety of the PEI chains and the carbon nanotube axis, Figure 3. The order parameters S(r) calculated for the P and D vectors showed that at T = 600 K, the planar fragments of ODPA-P3 orient themselves more strongly along the CNT surface. The first peak in the S(r) dependence for both polyetherimides is attributed to the orientation of the PEI planar fragments directly at the polymer-SWCNT interface. The order parameter for the P vector in ODPA-P3 shows the second peak at a distance of about 1.5 nm, which is absent in the S(r) dependence of aBPDA-P3. This fact demonstrates a rather strong structural ordering of the ODPA-P3 planar fragments even far from the nanotube surface. As opposed to the S(r) for ODPA-P3, the dependence S(r) for aBPDA-P3 has one maximum only and rapidly decreases to zero at increasing distances from the nanotube surface.
The emerging structural ordering of the ODPA-P3 polymer chain fragments may be interpreted as the initial crystallization stage of this polyetherimide [71,72]. In fact, previous results allowed the interpretation of polymer crystallization through several stages. First of all, polymer chain fragments orient themselves along the nanofiller surface. Further on, polymer chains orient and order themselves on larger spatial scales, which are accompanied by the order parameter growth in the system. This is followed by the entire ordering of polymer chains in the system, accompanied by an increase in the system local density and the formation of a structure similar to a crystalline one. Therefore, the structural ordering of the ODPA-P3 planar fragments near the CNT surface may be regarded as the initial polymer crystallization stage on local spatial scales. Similar results were obtained in our previous study for PEI R-BAPB [50].
where r is the distance from the nanotube axis to the polymer chain planar moiety and θ is the angle between the vector P or D directed along the planar moiety of the PEI chains and the carbon nanotube axis, Figure 3.
The order parameters ( ) S r calculated for the P and D vectors showed that at T = 600 К, the planar fragments of ODPA-P3 orient themselves more strongly along the CNT surface. The first peak in the ( ) S r dependence for both polyetherimides is attributed to the orientation of the PEI planar fragments directly at the polymer-SWCNT interface. The order parameter for the P vector in ODPA-P3 shows the second peak at a distance of about 1.5 nm, which is absent in the ( ) S r dependence of aBPDA-P3. This fact demonstrates a rather strong structural ordering of the ODPA-P3 planar fragments even far from the nanotube surface. As opposed to the ( ) S r for ODPA-P3, the dependence ( ) S r for aBPDA-P3 has one maximum only and rapidly decreases to zero at increasing distances from the nanotube surface.
The emerging structural ordering of the ODPA-P3 polymer chain fragments may be interpreted as the initial crystallization stage of this polyetherimide [71,72]. In fact, previous results allowed the interpretation of polymer crystallization through several stages. First of all, polymer chain fragments orient themselves along the nanofiller surface. Further on, polymer chains orient and order themselves on larger spatial scales, which are accompanied by the order parameter growth in the system. This is followed by the entire ordering of polymer chains in the system, accompanied by an increase in the system local density and the formation of a structure similar to a crystalline one. Therefore, the structural ordering of the ODPA-P3 planar fragments near the CNT surface may be regarded as the initial polymer crystallization stage on local spatial scales. Similar results were obtained in our previous study for PEI R-BAPB [50]. For the systems under the study, we also analyzed the distributions of the orientation angles θ of the vectors P and D relative to the CNT axis; the results are shown in Figure 4.
The orientation angle distribution ( Figure 4) demonstrates that in composites based on ODPA-P3, several layers in the vicinity of the CNT surface are forming where the PEI chain planar moieties are aligned along the CNT axis. Furthermore, the maximum in the distribution diagram is reached at the values of the angles θ within five to 20 degrees, and the structural ordering of the ODPA-P3 planar moieties extends to a distance of about 3 nm from the CNT axis (Figure 4a,b), i.e., considering the periodic boundary conditions, the orientation of the planar moieties takes place almost across the entire volume of the simulation cell. At the same time, the diagrams of the angle distribution between the carbon nanotube axis and vectors P and D for aBPDA-P3 do not confirm the For the systems under the study, we also analyzed the distributions of the orientation angles θ of the vectors P and D relative to the CNT axis; the results are shown in Figure 4.
The orientation angle distribution ( Figure 4) demonstrates that in composites based on ODPA-P3, several layers in the vicinity of the CNT surface are forming where the PEI chain planar moieties are aligned along the CNT axis. Furthermore, the maximum in the distribution diagram is reached at the values of the angles θ within five to 20 degrees, and the structural ordering of the ODPA-P3 planar moieties extends to a distance of about 3 nm from the CNT axis (Figure 4a,b), i.e., considering the periodic boundary conditions, the orientation of the planar moieties takes place almost across the entire volume of the simulation cell. At the same time, the diagrams of the angle distribution between the carbon nanotube axis and vectors P and D for aBPDA-P3 do not confirm the orientation of planar fragments near the CNT surface at a simulation time of~1 µs. From Figure 4c,d, it is clear that the  The absence of structural ordering of the aBPDA-P3 planar moieties near the CNT surface may be caused by the high temperatures used. Therefore, we studied the nanocomposite intermolecular structure at a lower temperature of Т = 580 К. To the end of the equilibration run at Т = 600 К, the temperature was instantly reduced to 580 К, and additional simulations of the system have been carried out for 1 μs. After that, the order parameter ( ) S r was calculated for the vectors P and D again; see Figure 5. aBPDA-P3 (c,d) and at temperatures of 580 K and 600 K subject to the distance from the CNT axis after a 1-μs production run.
The results obtained demonstrate that in the case of aBPDA-P3, the temperature drop from The absence of structural ordering of the aBPDA-P3 planar moieties near the CNT surface may be caused by the high temperatures used. Therefore, we studied the nanocomposite intermolecular structure at a lower temperature of T = 580 K. To the end of the equilibration run at T = 600 K, the temperature was instantly reduced to 580 K, and additional simulations of the system have been carried out for 1 µs. After that, the order parameter S(r) was calculated for the vectors P and D again; see Figure 5.  The absence of structural ordering of the aBPDA-P3 planar moieties near the CNT surface may be caused by the high temperatures used. Therefore, we studied the nanocomposite intermolecular structure at a lower temperature of Т = 580 К. To the end of the equilibration run at Т = 600 К, the temperature was instantly reduced to 580 К, and additional simulations of the system have been carried out for 1 μs. After that, the order parameter ( ) S r was calculated for the vectors P and D again; see Figure 5. aBPDA-P3 (c,d) and at temperatures of 580 K and 600 K subject to the distance from the CNT axis after a 1-μs production run.
The results obtained demonstrate that in the case of aBPDA-P3, the temperature drop from 600 to 580 К does not lead to any structural ordering near the CNT surface; see Figure 5. However, the temperature drop impacts the molecular structure of the ODPA-P3-based composite significantly. For this polymer, the temperature decrease leads to the increase in the second maximum value in the The results obtained demonstrate that in the case of aBPDA-P3, the temperature drop from 600 to 580 K does not lead to any structural ordering near the CNT surface; see Figure 5. However, the temperature drop impacts the molecular structure of the ODPA-P3-based composite significantly. For this polymer, the temperature decrease leads to the increase in the second maximum value in the S(r) dependence and higher values of S(r) at larger distances r from the CNT surface.
At a temperature of 580 K, we also calculated the distributions of the orientation angles θ between the planar moieties of the repeat unit of ODPA-P3 and aBPDA-P3 after a 1-µs production run; see Figure 6.
The analysis of the results obtained shows that the maximum distributions of orientation angle for ODPA-P3 composites shift to lower θ values after the temperature drop (compare Figures 6a and 4a). Furthermore, the formation of more pronounced second layer at a distance of ≈1.25-1.5 nm from the nanotube axis with planar moieties oriented largely at an angle of 25 degrees is observed, as well as the formation of oriented layers at larger distances. At the same time, the analysis of the distributions of orientation angles θ between the aBPDA-P3 planar fragments and the CNT axis at a temperature of 580 K failed to confirm the structural ordering of this polymer near the CNT surface; see Figure 6b. At a temperature of 580 К, we also calculated the distributions of the orientation angles θ between the planar moieties of the repeat unit of ODPA-P3 and aBPDA-P3 after a 1-μs production run; see Figure 6.
The analysis of the results obtained shows that the maximum distributions of orientation angle for ODPA-P3 composites shift to lower θ values after the temperature drop (compare Figures 6a   and 4a). Furthermore, the formation of more pronounced second layer at a distance of ≈1.25-1.5 nm from the nanotube axis with planar moieties oriented largely at an angle of 25 degrees is observed, as well as the formation of oriented layers at larger distances. At the same time, the analysis of the distributions of orientation angles θ between the aBPDA-P3 planar fragments and the CNT axis at a temperature of 580 К failed to confirm the structural ordering of this polymer near the CNT surface; see Figure 6b. Notably, different structural properties of nanocomposites based on thermoplastic PEIs ODPA-P3 and aBPDA-P3 may also be observed during the analysis of their instant configurations. Figure 7 shows snapshots of the PEIs under study after a 1-μs simulation at a temperature of 580 K. Notably, different structural properties of nanocomposites based on thermoplastic PEIs ODPA-P3 and aBPDA-P3 may also be observed during the analysis of their instant configurations. Figure 7 shows snapshots of the PEIs under study after a 1-µs simulation at a temperature of 580 K.    The analysis of Figure 7a suggests that chains of crystallizable polyetherimide ODPA-P3 orient themselves near the surface of the CNT alongside the nanofiller surface. At the same time, there is no orientation of planar moieties in the amorphous polyetherimide aBPDA-P3 near the CNT surface; see Figure 7b.
As opposed to aBPDA-P3, structural ordering of the ODPA-P3 planar moieties is due to the different chemical structures of their dianhydride moieties. In this study, we calculated the distributions P(ϕ) of the dihedral angles of internal rotation (ω, χ) in ODPA-P3 and aBPDA-P3; see Figure 8a. The analysis of Figure 7a suggests that chains of crystallizable polyetherimide ODPA-P3 orient themselves near the surface of the CNT alongside the nanofiller surface. At the same time, there is no orientation of planar moieties in the amorphous polyetherimide aBPDA-P3 near the CNT surface; see Figure 7b.
As opposed to aBPDA-P3, structural ordering of the ODPA-P3 planar moieties is due to the different chemical structures of their dianhydride moieties. In this study, we calculated the distributions   The analysis of the results obtained shows that distributions of dihedral angles ω of ODPA-P3 and aBPDA-P3 almost fully coincide. The most probable values of the angles ω in both polyetherimides amount to approximately 10 • and 160 • . However, distributions of dihedral angles χ of the two polyetherimides differ from each other. For the aBPDA-P3, the most probable value of the dihedral angle χ is equal to 90 • , and the most probable value of this angle for the ODPA-P3 amounts to approximately 160 • ; see Figure 8b. Therefore, the phthalimide moieties in aBPDA are unfolded perpendicularly in relation to each other, which prevents the structural ordering of the planar moieties near the CNT surface. To this end, this spatial confinement both influences the adsorption of the polymer chain to the CNT surface in the dianhydride fragment and prevents the structural ordering of phenylene rings of the diamine fragment of aBPDA-P3. At the same time, the dianhydride fragment of ODPA-P3 has practically a planar structure facilitating the orientation of phthalimide moieties and phenylene rings along the CNT surface.
Therefore, we succeeded in simulating the experimental effect reported previously on the ODPA-P3 crystallization initiation after the CNT incorporation [63,64]. Thus, results obtained here provide at least the correct choice of the force field to predict differences in the PEIs' structural properties.
As mentioned in the Introduction, the ordering of planar moieties of the polymer repeat unit may be regarded as the initial crystallization stage of the thermoplastic polyetherimide ODPA-P3 [71,72]. It is assumed that the molecular mechanism of the initial crystallization stage of the ODPA-P3 is determined by either the CNT presence or influence of π-π interactions between CNT and ODPA-P3 atoms. Then, we will try to answer the question as to what the main factor is in the ordering process of ODPA-P3 planar moieties towards the CNT surface.
Notably, the correct consideration of π-π interactions in molecular dynamics is a non-trivial task when classical mechanical force fields are used, which do not account for quantum effects. Nevertheless, since in reality the molecular dynamics force fields are semi-empirical, π-π interactions are efficiently accounted for by excluded volume interactions. Moreover, it was shown that simulations with classical force fields allow reproducing structural properties of molecular systems determined by π-π interactions [73]. In our previous studies [50,51,58], we have investigated the influence of electrostatic interactions on the structure of nanocomposites based on the polyetherimide R-BAPB. Ordered structures in these systems were formed during the simulation both with and without partial charges. Since, as it was shown in these studies, the ordering is related to the orientation of aromatic fragments of the polyetherimide chains towards the CNT axis and there are π-π interactions between them in real systems, it may be acceptable to state that π-π are quite accurately accounted for in the Gromos53a5 force field through excluded volume interactions between specific atom types.
Thus, changes in the form and parameters of the excluded volume interactions between the CNT atoms and polymer chain atoms in the simulation allow for investigating the influence of each of the two effects in question on the polymer-CNT interface structure. Therefore, we studied how changes in the form and parameters of the excluded volume interactions between the carbon atoms of the polymer and the CNT influence the structural ordering of the ODPA-P3 planar fragments near the CNT surface. In the force field Gromos53a5 [67,74], the excluded volume interactions are set using the Lennard-Jones (LJ) potential: where ε is the potential well depth, σ is the distance at which U L−J r ij becomes equal to zero (van der Waals radius), and r ij is the distance between the i-th and j-th atoms. At the beginning, to verify the impact of the geometric confinement through the CNT in the nanocomposites, during the initial crystallization stage of thermoplastic polyetherimides, we performed a simulation of composites based on ODPA-P3 using the modified LJ potential written in the form of the pure repulsive Weeks-Chandler-Anderson (WCA) potential, W r ij [75]: where r 0 = 6 √ 2σ is the distance corresponding to the minimum of the function U L−J r ij . LJ and WCA potentials are given in Figure 9.
Polymers 2017, 9,548 12 of 20 where 0 r = σ 6 2 is the distance corresponding to the minimum of the function ( ) − L J ij U r . LJ and WCA potentials are given in Figure 9. The use of the WCA potential to describe interactions between the carbon atoms of the CNT and heterocyclic rings of ODPA-P3 enables us to exclude the attraction of planar fragments to the CNT surface due to π-π interactions. In this case, the orientation of the ODPA-P3 planar moieties could be the result of spatial confinements promoted by the CNT presence.
The change of the form of excluded volume interaction potential was followed by additional simulation for 1 μs. Then, we calculated the dependences of the order parameter ( ) S r and density ( ) ρ r of the ODPA-P3 on the distance from the CNT axis and distributions of the orientation angles θ between the ODPA-P3 planar moieties and CNT axis; see Figure 10.
The analysis of the data presented in Figure 10 allows us to conclude that 1-μs simulation using the WCA potential instead of the LJ potential leads to the destruction of the structural ordering of the ODPA-P3 planar segments near the CNT surface. The dependence of the order parameter ( )  The use of the WCA potential to describe interactions between the carbon atoms of the CNT and heterocyclic rings of ODPA-P3 enables us to exclude the attraction of planar fragments to the CNT surface due to π-π interactions. In this case, the orientation of the ODPA-P3 planar moieties could be the result of spatial confinements promoted by the CNT presence.
The change of the form of excluded volume interaction potential was followed by additional simulation for 1 µs. Then, we calculated the dependences of the order parameter S(r) and density ρ(r) of the ODPA-P3 on the distance from the CNT axis and distributions of the orientation angles θ between the ODPA-P3 planar moieties and CNT axis; see Figure 10.
The analysis of the data presented in Figure 10 allows us to conclude that 1-µs simulation using the WCA potential instead of the LJ potential leads to the destruction of the structural ordering of the ODPA-P3 planar segments near the CNT surface. The dependence of the order parameter S(r) on the distance features a drop of the second and third maximum values and the shift of their position for a larger distance away from the CNT axis. Similar conclusions can be formulated from the analysis of the dependence ρ(r) of ODPA-P3. The dependence ρ(r) obtained using the WCA potential shows the position of the first and the second maxima significantly lower and shifted to larger r values relative to the position of the maxima in the dependence ρ(r) calculated using the LJ potential, which provides evidence that the ODPA-P3 polymer chains cease to be attracted to the CNT surface.
Comparative analysis of the distributions of the orientation angles θ of the ODPA-P3 planar moieties near the CNT surface for the WCA potential ( Figure 10c) and LJ potential (Figure 6a) shows a practically complete disappearance of the structural ordering of the ODPA-P3 planar fragments in the case of the WCA potential. The analysis of the dependence S(r) and distributions of the orientation angles θ of the ODPA-P3 planar fragments relative to the CNT also confirms the disappearance of all structural ordering of the ODPA-P3 planar moieties near the CNT surface; see Figures S4 and S5 in the Supplementary Materials. Notably, when the simulation was extended to 2 µs, there was still no structural ordering of the ODPA-P3 planar moieties near the CNT surface; see Figure S5 in the Supplementary Materials. This confirms the lack of interaction between the carbon atoms of the PEI and CNT, i.e., insufficient strength of π-π interactions at the polymer-CNT phase interface. In the simulation, this leads to the disappearance of the structural ordering of the polymer planar moieties near the CNT surface. Thus, spatial confinements are not the determining factor for initiation of ODPA-P3 crystallization.  Figure S5 in the Supplementary Materials. This confirms the lack of interaction between the carbon atoms of the PEI and CNT, i.e., insufficient strength of π-π interactions at the polymer-CNT phase For an additional check of the influence of π-π interactions on the polymer-CNT interphase structure, we performed simulations with varying depth of the potential well ε of the interaction energy (4) between the carbon atoms of the CNT and the ODPA-P3 heterocyclic rings. As was shown in the present study, the orientation of planar moieties near the CNT surface has been found in computer simulation both with and without taking into account electrostatic interactions. Since, the ordering is related to the orientation of aromatic fragments of the PEI chains towards the CNT axis and there are π-π interactions between them in real systems, it may be acceptable to state that π-π interactions are quite accurately accounted for in the Gromos53a5 force field through volume interactions between specific atom types in agreement with results of our previous investigations of the structural properties of R-BAPB-based nanocomposites [50,51,58].
Varying the value of the parameter ε, we sought to efficiently relax or enhance the π-π interactions at the polymer-CNT phase interface. Figure 11 shows dependences U L−J r ij for various values of ε used in our study. At that, the van der Waals radius value σ of the carbon atoms remained constant. energy (4) between the carbon atoms of the CNT and the ODPA-P3 heterocyclic rings. As was shown in the present study, the orientation of planar moieties near the CNT surface has been found in computer simulation both with and without taking into account electrostatic interactions. Since, the ordering is related to the orientation of aromatic fragments of the PEI chains towards the CNT axis and there are π-π interactions between them in real systems, it may be acceptable to state that π-π interactions are quite accurately accounted for in the Gromos53a5 force field through volume interactions between specific atom types in agreement with results of our previous investigations of the structural properties of R-BAPB-based nanocomposites [50,51,58].
Varying the value of the parameter ε, we sought to efficiently relax or enhance the π-π interactions at the polymer-CNT phase interface. Figure 11 shows dependences ( ) − L J ij U r for various values of ε used in our study. At that, the van der Waals radius value σ of the carbon atoms remained constant. Changes in the potential well depth ε of the excluded volume interactions between the carbon atoms of the CNT and ODPA-P3 were followed by the additional simulation for 1 μs. Further, we investigated the structural properties of ODPA-P3. To this end, we calculated ( ) S r , ( ) ρ r and distributions of the orientation angles θ between the PEI planar fragments and the CNT axis; see Figure 12.
The analysis of the dependence of the order parameter ( ) S r for the vector P shows that the 10-times reduction of ε leads to ODPA-P3 planar fragments moving away from the CNT surface. The dependence ( ) S r has only one maximum, which slowly decreases when the distance is growing; see Figure 12a. This may be attributed to the fact that after a 10-times reduction of the potential well depth ε , the ODPA-P3 atoms cannot adsorb on the CNT surface due to thermal fluctuations. When the parameter ε is increased by five-times, the dependence ( ) S r demonstrates the rise of the first maximum height, which confirms the fact that the ODPA-P3 heterocyclic rings find it more energetically favorable to orient themselves near the CNT surface. The analysis of the dependence ( ) ρ r for varying ε values shows that the increase in ε also causes the height of the first and second maxima to grow in the dependence ( ) ρ r of ODPA-P3, as well; see Figure 12b. With ε increasing, this polyetherimide starts to be more strongly attracted to the CNT surface. The increase or reduction of the potential well depth ε between the carbon atoms of the CNT and ODPA-P3 influences both the structural ordering of the ODPA-P3 planar moieties near the CNT surface and the polymer density near the CNT, i.e., on the adsorption of polymer chains on the CNT surface. These conclusions are Figure 11. Van der Waals energy U L−J r ij for the varying potential well depth ε multiplied by 0.1-and five-times as compared to the initial parameters of the Gromos53a5 force field (ε = 1 ).
Changes in the potential well depth ε of the excluded volume interactions between the carbon atoms of the CNT and ODPA-P3 were followed by the additional simulation for 1 µs. Further, we investigated the structural properties of ODPA-P3. To this end, we calculated S(r), ρ(r) and distributions of the orientation angles θ between the PEI planar fragments and the CNT axis; see Figure 12.
The analysis of the dependence of the order parameter S(r) for the vector P shows that the 10-times reduction of ε leads to ODPA-P3 planar fragments moving away from the CNT surface. The dependence S(r) has only one maximum, which slowly decreases when the distance is growing; see Figure 12a. This may be attributed to the fact that after a 10-times reduction of the potential well depth ε, the ODPA-P3 atoms cannot adsorb on the CNT surface due to thermal fluctuations. When the parameter ε is increased by five-times, the dependence S(r) demonstrates the rise of the first maximum height, which confirms the fact that the ODPA-P3 heterocyclic rings find it more energetically favorable to orient themselves near the CNT surface. The analysis of the dependence ρ(r) for varying ε values shows that the increase in ε also causes the height of the first and second maxima to grow in the dependence ρ(r) of ODPA-P3, as well; see Figure 12b. With ε increasing, this polyetherimide starts to be more strongly attracted to the CNT surface. The increase or reduction of the potential well depth ε between the carbon atoms of the CNT and ODPA-P3 influences both the structural ordering of the ODPA-P3 planar moieties near the CNT surface and the polymer density near the CNT, i.e., on the adsorption of polymer chains on the CNT surface. These conclusions are substantiated also by the analysis of the distribution of orientation angles θ of the PEI planar fragments relative to the CNT axis. The results obtained have demonstrated that the reduction of the parameter ε by 10-times in the distribution of the orientation angles θ (Figure 12c) causes the second and third ordered polymer layers near the CNT surface to disappear, whereas these layers were observed earlier using the standard parameters of the force field; see Figure 6a. In the case of the five-times increase in the parameter ε, distributions of the orientation angles θ feature the first peak at a distance of ≈0.75-0.8 nm, where the planar fragments of the ODPA-P3 repeat unit are oriented at an angle of~10 • , and the second peak at a distance of ≈1.25-1.5 nm, where the planar moieties are oriented at an angle of ≈20 • ; see Figure 12d. The analysis of the results obtained shows a strong influence of specific π-π interactions of the polymer-CNT interface structure. Similar conclusions are suggested by the analysis of the dependences S(r) and distribution angles θ between the vector D and the CNT axis; see Figure S6 in the Supplementary Materials. Notably, if a simulation accounts for the attraction between the carbon atoms of the CNT and polymer, even in the case of the 10-times reduction of the parameter S(r) as compared to its initial value in the force field, weak structural ordering of the ODPA-P3 planar fragments still remains near the CNT surface after 1 µs.
istance of ≈0.75-0.8 nm, where the planar fragments of the ODPA-P3 repeat unit are oriented angle of ~10°, and the second peak at a distance of ≈1.25-1.5 nm, where the planar moieties a iented at an angle of ≈20°; see Figure 12d. The analysis of the results obtained shows a stron fluence of specific π-π interactions of the polymer-CNT interface structure. Similar conclusions ar ggested by the analysis of the dependences ( ) S r and distribution angles θ between the vector d the CNT axis; see Figure S6 in the Supplementary Materials. Notably, if a simulation accoun r the attraction between the carbon atoms of the CNT and polymer, even in the case of the 10-time duction of the parameter ( ) S r as compared to its initial value in the force field, weak structur dering of the ODPA-P3 planar fragments still remains near the CNT surface after 1 μs. five-times as compared to the initial parameters of the Gromos53a5 force field after a 1-μs production run.
Therefore, structural ordering of the ODPA-P3 planar fragments near the CNT surface is th sult of π-π interactions between the carbon atoms of the CNT aromatic rings and polym terocyclic rings. Regarding this structural ordering as the initial stage of the thermoplast PA-P3 crystallization due to the nucleating action of the carbon nanofiller, it can be affirmed th e presence of π-π interactions between the heterocyclic polymer and nanofiller surface initiate PA-P3 crystallization. We believe that this mechanism may be valid for other heterocycl lymers, as well. Notably, these conclusions on the influence of π-π interactions on the initial stag the thermoplastic polyetherimide ODPA-P3 crystallization correlate very well with th perimental findings on the adsorption of low molecular compounds with (locally) conjugate omatic cycles on the CNT surface [76,77].

Conclusions
In the present study, we used microsecond-long atomistic molecular-dynamics simulations t vestigate the molecular mechanisms responsible for the initial stage of thermoplast Figure 12. Dependences of the (a) order parameter S(r) for the vector P and (b) density ρ(r) on the distance from the CNT axis for varying values of the parameter ε. Distributions of angles θ between the vector P and CNT axis at the increase in the potential well depth ε by (c) 0.1-and (d) five-times as compared to the initial parameters of the Gromos53a5 force field after a 1-µs production run. Therefore, structural ordering of the ODPA-P3 planar fragments near the CNT surface is the result of π-π interactions between the carbon atoms of the CNT aromatic rings and polymer heterocyclic rings. Regarding this structural ordering as the initial stage of the thermoplastic ODPA-P3 crystallization due to the nucleating action of the carbon nanofiller, it can be affirmed that the presence of π-π interactions between the heterocyclic polymer and nanofiller surface initiates ODPA-P3 crystallization. We believe that this mechanism may be valid for other heterocyclic polymers, as well. Notably, these conclusions on the influence of π-π interactions on the initial stage of the thermoplastic polyetherimide ODPA-P3 crystallization correlate very well with the experimental findings on the adsorption of low molecular compounds with (locally) conjugated aromatic cycles on the CNT surface [76,77].

Conclusions
In the present study, we used microsecond-long atomistic molecular-dynamics simulations to investigate the molecular mechanisms responsible for the initial stage of thermoplastic polyetherimide crystallization. To this end, we simulated the structural properties of two amorphous thermoplastic polyetherimides, ODPA-P3 and aBPDA-P3, in the presence of SWCNTs. We observed that the planar phthalimide and phenyl rings of ODPA-P3 show an increase in order when they approach the CNT surface. This, in our opinion, can be considered as the initial stage of ODPA-P3 crystallization [63]. At the same time, no orientation ordering of the aBPDA-P3 planar fragments could be observed. The difference in structural properties of the PEIs investigated is mainly dictated by the difference in chemical composition of the dianhydride moiety.
To provide computational insights into the molecular mechanisms of the initial stage of PEI crystallization, we simulated the interface structure of the ODPA-P3 nanocomposite reinforced with CNTs. We explored two possible scenarios that could lead to ordering: the polymer orients itself near the CNT surface due to the spatial confinements or because of the strong influence of the π-π interactions at the polymer-CNT interface. It has been established here that in the absence of the attraction between the CNT surface and the PEI matrix, no structural ordering of the ODPA-P3 planar fragments near the CNT surface is observed. We have shown that by varying the strength of the van der Waals interactions between the CNT atoms and the polymer segments, it is possible to weaken or enhance the structural ordering of the ODPA-P3 planar fragments near the CNT surface. Therefore, we propose that the molecular mechanism of the thermoplastic PEI crystallization initiation, as a result of CNT incorporation, is the result of favorable π-π interactions between the heterocyclic polymer rings and the CNT surface where all carbon atoms are in sp 2 hybridization state.
In the work presented in this paper, we did not consider CNT agglomeration effects, as this would have complicated our initial modeling efforts. CNT agglomeration (bundling) will most definitely affect polymer-CNT interactions and complicate PEI crystallization. We have already started to look at PEI crystallization in the presence of different CNT packing motives, e.g., crystallization at intersecting CNTs, and this work will be the basis for our next publication.