Molecular Dynamics Simulation and Structure Changes of Polyester in Water and Non-Aqueous Solvents

Studying the changes in the microstructure of polyester (PET) in water and non-aqueous solvents is important to understand the swelling mechanism of PET, which can help to reduce water pollution during the dyeing process. This study uses molecular models of PET, water, and decamethyl-cyclopentasiloxane (D5) and employs molecular dynamics method to simulate the influence of solvents on the microstructure of PET. The results show that the glass transition temperature (Tg) of the pure PET system is close to the experimental value. The Tg of PET decreases with the addition of water and D5 solvents, and the free volume after adding D5 is higher compared to the free volume after adding water. Through molecular dynamics simulation of PET microstructure, it is found that D5 has a better SWELLING effect on PET than water.


Introduction
Polyethylene terephthalate (PET, polyester) fiber has a very tight structure and strong hydrophobicity. The cross-section of the PET fiber is circular under the microscope, its longitudinal thickness is uniform, and its surface is smooth [1][2][3]. During the PET dyeing process, the dispersed dye must enter the PET fiber through a narrow gap, which poses great difficulty. Therefore, to achieve the dyeing of PET [4][5][6], dyeing conditions of carrier, high temperature, or hot melt are generally required [7][8][9][10][11][12]. This kind of dyeing technology requires more equipment and higher cost. Moreover, the dyeing of polyester fabrics consumes a large amount of water, and the cost of waste treatment is high.
In order to solve the above problems, it is necessary to use other media with different properties from water. Hu [13] used 15 kinds of pure disperse red dyes to dye polyester/nylon microfiber fabrics for 60 min using self-developed supercritical dyeing equipment. Infrared spectral analysis showed that no new groups were introduced after the fabric was dyed, and the original properties of the fabric were maintained. Li et al. [14] reported that the non-aqueous solvent D5 has a certain swelling effect on PET. Wang et al. [15] added an accelerator to the D5 solvent dyeing system, and found that the accelerator increased the dye uptake rate of disperse dyes on PET. After measuring the fiber diameter with a confocal laser scanning microscope, they found that the accelerant improved the swelling of PET in silicone solvents. Therefore, it is necessary to study the influence of aqueous solvents and non-aqueous solvents on the microstructure of PET. The change in the microstructure of PET fiber is a necessary link in the dyeing process. To realize the dyeing of PET in solvent, the key lies in grasping the changes in PET microstructure in solvent. Research from a microscopic point of view to explore the variation of PET fiber structure under the action of water and non-aqueous solvents is conducive to deepening the knowledge and understanding of the dyeing mechanism.
Molecular dynamics (MD) simulation helps to understand the microstructure of polymers, analyze the changes in polymer structure at the molecular level, and provide new research methods [16,17]. MD simulation can explore the microscopic behavior of particles in the diffusion-limited aggregation process at the atomic level to reveal the aggregation mechanism of particles in aqueous solution [18]. MD simulation is also widely used to study the mechanism of interaction between different materials [19]. Sheng [20] studied the glass transition temperature of polyethylene/graphene nanocomposites through MD. The specific volumes of the three systems (polyethylene, polyethylene with small graphene flakes, and two small graphene flakes) were obtained as a function of temperature. It was found that the glass transition temperature decreases with the increase in content of graphene. Van der Waals energy changes significantly with the increase in content of graphene, and torsional energy also plays an important role in the glass transition of polymers. These results indicate that graphene can promote the movement of polymer segments and reduce the glass transition temperature of the polymer. Zhang et al. [21] used MD to study the swelling performance of polyvinyl alcohol (PVA) in water/ethanol solution, and found that pores containing water and ethanol were formed between the PVA chains. As the swelling degree increases, the number and size of the pores both increase. The diffusion coefficients of water and ethanol in swollen PVA increase linearly with the degree of swelling. Shanks and Pavel [22][23][24][25] used MD simulation to study the diffusion of methane and other gases and small molecule penetrants in aromatic PET, amorphous PET, and isomeric PET. They found that the diffusion coefficient and free volume and its distribution are related. Qiang et al. [26] used MD to study the cycling behavior of polyethylene (PE) under different loading conditions below the glass transition temperature at the nanoscale. The results show that van der Waals energy and dihedral angle energy are the primary factors affecting the cyclic behavior of PE. Wang et al. [27] simulated the glass transition behavior and mechanical properties of PET, PET/silica nanocomposites, and PET/hydroxylated silica nanocomposites. It was found that the content of silica had a certain influence on the T g of PET. Moreover, Guo et al. [28] simulated the adsorption of sulfamethazine (SMT) on six microplastics, such as PET by MD simulation. Their results showed that the main mechanism of adsorption was electrostatic and Van der Waals interaction.
In this study, MD was used to study the changes in the microstructure of PET in the three systems of pure PET, PET/H 2 O, and PET/D5 at different temperatures. The swelling mechanism of PET was studied from a microscopic point of view, which can provide theoretical support for further research on PET dyeing with non-aqueous solvents.

Materials and Methods
The MD simulation was carried out using the Forcite and Amorphous Cell modules of Material Studio 8.0 software (Accelrys, San Diego, CA, USA) with the condensed-phase optimized molecular potentials for atomistic simulation studies (COMPASS, Center for Molecslar Science, Institute of Chemistry, Beijing, China) [29] force field, which is widely adopted to predict the structure of polymers. In the COMPASS force field, the total energy of the system is expressed as: where E total is the total energy, E b is bond stretching energy, and E θ is angle bending energy. E ϕ is dihedral torsion energy, E χ is out-of-plane energy, E cross is cross term interaction energy, E ele is electrostatic interaction energy, and E vdw is van der Waals interaction energy. The electrostatic and van der Waals forces were calculated by using the atom-based summation method with a cut off value of 12.5 Å. Figure 1 shows the construction process of the molecular simulation. First, a PET molecular chain containing 50 repeating units was constructed (Figure 1c), and the cell energy of each cell was minimized using the Smart Minimizer method. After the energy of the periodic cells was reduced through the geometric optimization process to obtain a more stable structure, the NPT (constant molecule, pressure, and temperature) ensemble was selected for annealing. The annealing process was performed for 20 cycles, and a total of 600 ps was run. After annealing, the system equilibrium was completed. First, 200 ps NPT ensemble simulation was conducted and then 500 ps NVT (constant molecule number, volume and temperature) ensemble simulation was carried out. Andersen [30] thermostat and Berendsen [31] barostat were used for temperature control and pressure control, respectively. System 2 and System 3 were subjected to the same simulation process. Finally, equilibrated simulation boxes were used to analyze the fraction of free volume (FFV) and other characteristic parameters (Figure 1g).  Figure 1 shows the construction process of the molecular simulation. First, a PET molecular chain containing 50 repeating units was constructed (Figure 1c), and the cell energy of each cell was minimized using the Smart Minimizer method. After the energy of the periodic cells was reduced through the geometric optimization process to obtain a more stable structure, the NPT (constant molecule, pressure, and temperature) ensemble was selected for annealing. The annealing process was performed for 20 cycles, and a total of 600 ps was run. After annealing, the system equilibrium was completed. First, 200 ps NPT ensemble simulation was conducted and then 500 ps NVT (constant molecule number, volume and temperature) ensemble simulation was carried out. Andersen [30] thermostat and Berendsen [31] barostat were used for temperature control and pressure control, respectively. System 2 and System 3 were subjected to the same simulation process. Finally, equilibrated simulation boxes were used to analyze the fraction of free volume (FFV) and other characteristic parameters (Figure 1g). The basic parameters of PET and the combined system with H2O and D5 solvents are listed in Table 1. The volume and density of the polymer after optimization are also shown. The final density of the ensemble simulated PET system is 1.241 × 10 3 kg/m 3 , which is close to the experimentally measured density of amorphous PET polymer (1.330 × 10 3 kg/m 3 ).  The basic parameters of PET and the combined system with H 2 O and D5 solvents are listed in Table 1. The volume and density of the polymer after optimization are also shown. The final density of the ensemble simulated PET system is 1.241 × 10 3 kg/m 3 , which is close to the experimentally measured density of amorphous PET polymer (1.330 × 10 3 kg/m 3 ).

Evaluation of System Balance
In order to obtain reasonable parameters of the system, the system must be simulated for a long enough time for it to reach an equilibrium state. There are two main criteria for judging whether the system reaches equilibrium: one is temperature equilibrium, which requires that the fluctuation of the equilibrium value should not exceed 10 K. The second is energy balance, which requires the energy to be constant or fluctuate up and down along a constant value. Figure 2 presents the temperature and energy balance diagrams of PET under the NPT ensemble of 360 K and 150 ps, respectively. The standard deviation of the calculated temperature is 4.99 K. The standard deviations of the total energy, potential energy, non-bonding energy, and kinetic energy are 107.63 kcal/mol, 69.92 kcal/mol, 28.65 kcal/mol, and 65.51 kcal/mol, respectively, indicating that the energy of the simulated system has reached an equilibrium state.

Evaluation of System Balance
In order to obtain reasonable parameters of the system, the system must be simulated for a long enough time for it to reach an equilibrium state. There are two main criteria for judging whether the system reaches equilibrium: one is temperature equilibrium, which requires that the fluctuation of the equilibrium value should not exceed 10 K. The second is energy balance, which requires the energy to be constant or fluctuate up and down along a constant value. Figure 2 presents the temperature and energy balance diagrams of PET under the NPT ensemble of 360 K and 150 ps, respectively. The standard deviation of the calculated temperature is 4.99 K. The standard deviations of the total energy, potential energy, non-bonding energy, and kinetic energy are 107.63 kcal/mol, 69.92 kcal/mol, 28.65 kcal/mol, and 65.51 kcal/mol, respectively, indicating that the energy of the simulated system has reached an equilibrium state.

Simulation Analysis of Glass Transition Temperature of PET
In this study, the glass transition temperature (Tg) of Systems 1-3 was simulated and measured. MD simulation was used to determine the change curve of specific volume with temperature under constant pressure in the range of 200-500 K.
The initial stage temperature was set to 500 K, and the temperature was reduced by 20 K after each stage of the experiment was completed. The kinetic final equilibrium conformation of the previous stage was used as the initial conformation of the later stage MD simulation. At each temperature point, a 50 ps NVT ensemble dynamics simulation was performed, and then a 150 ps NPT ensemble dynamics simulation was performed. After the system reached equilibrium, data acquisition and analysis were performed.
The simulated specific volume versus temperature plot of PET is shown in Figure 3. The intersection of the two straight lines obtained by fitting is the glass transition temperature of the system [20]. Figure 3a shows that the glass transition temperature of PET is 370 K, which is close to the corresponding experimentally measured value of 342-350 K [32]. This indicates that the PET model established in this study can reflect the properties of the polymer, which further verifies its rationality for experiments. It can be seen from Figure 3b that the Tg of PET, after adding H2O solvent, is 368 K. Compared with the pure PET system, the temperature drop is smaller, indicating that the effect of water on reducing the Tg of PET is not obvious. It can be seen from Figure 3c that the Tg of PET and D5 system is 346 K. The simulation results show that the Tg of PET decreases significantly after adding D5, and the decrease of Tg is beneficial to reduce the dyeing temperature of PET and accelerate the dyeing of PET. This result indicates that the PET polymer model

Simulation Analysis of Glass Transition Temperature of PET
In this study, the glass transition temperature (T g ) of Systems 1-3 was simulated and measured. MD simulation was used to determine the change curve of specific volume with temperature under constant pressure in the range of 200-500 K.
The initial stage temperature was set to 500 K, and the temperature was reduced by 20 K after each stage of the experiment was completed. The kinetic final equilibrium conformation of the previous stage was used as the initial conformation of the later stage MD simulation. At each temperature point, a 50 ps NVT ensemble dynamics simulation was performed, and then a 150 ps NPT ensemble dynamics simulation was performed. After the system reached equilibrium, data acquisition and analysis were performed.
The simulated specific volume versus temperature plot of PET is shown in Figure 3. The intersection of the two straight lines obtained by fitting is the glass transition temperature of the system [20]. Figure 3a shows that the glass transition temperature of PET is 370 K, which is close to the corresponding experimentally measured value of 342-350 K [32]. This indicates that the PET model established in this study can reflect the properties of the polymer, which further verifies its rationality for experiments. It can be seen from Figure 3b that the T g of PET, after adding H 2 O solvent, is 368 K. Compared with the pure PET system, the temperature drop is smaller, indicating that the effect of water on reducing the T g of PET is not obvious. It can be seen from Figure 3c that the T g of PET and D5 system is 346 K. The simulation results show that the T g of PET decreases significantly after adding D5, and the decrease of T g is beneficial to reduce the dyeing temperature of PET and accelerate the dyeing of PET. This result indicates that the PET polymer model established in this paper can reflect the characteristics of the polymer, and the parameters used are correct, which is of instructive significance for conducting the related experiments. This also verified the rationality of the proposed model. On this basis, 413 K [33] was determined to be the suitable kinetic simulation temperature for analyzing the microstructure of PET in the dyeing process. established in this paper can reflect the characteristics of the polymer, and the parameters used are correct, which is of instructive significance for conducting the related experiments. This also verified the rationality of the proposed model. On this basis, 413 K [33] was determined to be the suitable kinetic simulation temperature for analyzing the microstructure of PET in the dyeing process.

Simulation Analysis of Molecular Ability and Diffusion Characteristics of Polyester Molecular Chains
The mean square displacement (MSD) represents the mobility of particles in the system and is used to characterize the mobility of molecules. The calculation formula of MSD is as follows: where ri(0) and ri(t) are the positions of particle i at the initial time and time t, respectively; the symbol indicates the average value for all particles.
The particle diffusion coefficient D can be calculated from the slope of MSD using the Einstein relationship as follows:

Simulation Analysis of Molecular Ability and Diffusion Characteristics of Polyester Molecular Chains
The mean square displacement (MSD) represents the mobility of particles in the system and is used to characterize the mobility of molecules. The calculation formula of MSD is as follows: where r i (0) and r i (t) are the positions of particle i at the initial time and time t, respectively; the symbol indicates the average value for all particles. The particle diffusion coefficient D can be calculated from the slope of MSD using the Einstein relationship as follows: (3) Figure 4 shows the MSD results of PET molecular chain in each of the three systems at a temperature of 413 K. It can be seen that the chain movement of the swollen PET macromolecule is greater than that of the PET molecular chain. In order to analyze molecular diffusion behavior more clearly, the MSD curves of three systems are calculated by Log10 and fitted by segments (as shown in Figure 5). In Figure 5, a is the slope of the fitting straight line, which is anomalous (non-Einstein) behavior when a < 1, and Einstein behavior when m is close to 1 [34]. As can be seen from the figure, at the back end of the simulation time, three systems transform to Einstein diffusion. During this time, the diffusion coefficient can be calculated by Einstein Equation (3). The diffusion coefficients of the three systems are shown in Table 2.
Log10 and fitted by segments (as shown in Figure 5). In Figure 5, a is the slope of the fitting straight line, which is anomalous (non-Einstein) behavior when a < 1, and Einstein behav ior when m is close to 1 [34]. As can be seen from the figure, at the back end of the simu lation time, three systems transform to Einstein diffusion. During this time, the diffusion coefficient can be calculated by Einstein Equation (3). The diffusion coefficients of th three systems are shown in Table 2.  straight line, which is anomalous (non-Einstein) behavior when a < 1, and Einstein behavior when m is close to 1 [34]. As can be seen from the figure, at the back end of the simulation time, three systems transform to Einstein diffusion. During this time, the diffusion coefficient can be calculated by Einstein Equation (3). The diffusion coefficients of the three systems are shown in Table 2.   It can be seen from Table 2 that as the temperature increases, the molecular chain obtains more energy, so the movement is more active. At the same time, the higher the temperature, the softer the PET molecular chain. Consequently, more free volume is generated, and more space or path is provided for the diffusion of solvent molecules. In water and D5, the diffusion coefficient of D5 is higher than that of water.

Simulation Analysis of PET Microvoids and Their Distribution
Free volume (FV) refers to the volume in the system that is not occupied by molecules. It is an important parameter used to characterize the structure of polymers. The molecular chains can only move due to the existence of free volume. FFV is the percentage of free volume to the total volume of the system. Figure 6 shows the schematic of van der Waals surface and Connolly surface. The van der Waals surface intersects with the van der Waals radii of the atoms in the structure and Connolly surface is at the boundary between the Connolly probe and the atoms. Free volume is defined as the volume circled by the Connolly surface. FFV was calculated by the ratio of free volume to the total volume of the model. In this study, the Atom Volumes and Surfaces tool in MS software was used to obtain the accessible FV of the probe by the Connolly surface method. The free volume was calculated when the probe radius was r and r + ∆r in turn. The difference between the two is the free volume of the cavity with radius r~r + ∆r.
FFV(∆r) = FFV(r) − FFV(r + ∆r) In the above formula, FFV(r) and FFV(r + r), respectively, represent the percentage of the free volume that can be contacted when the probe radius is r and r + r to the total volume of the system; FFV( r) represents the free volume fraction of holes in the range of r~r + ∆r. Figure 7 shows the volume morphology with different diameters at the same temperature. The blue and white areas represent the free volume created by transient gaps caused by inefficient chain packing or chain rearrangement. The free volume decreases with increase in probe diameter. Figure 8 shows the free volume fraction of PET in the three systems at different temperatures. It can be seen that, as the temperature rises over the range from 400 K to 440 K, the free volume fraction of PET increases. In this temperature range, the free volume fraction of PET increases after adding water, indicating that the addition of water has a certain swelling effect on PET. However, the increase in free volume after adding D5 is larger.   Figure 8 shows the free volume fraction of PET in the three systems at different tem peratures. It can be seen that, as the temperature rises over the range from 400 K to 440 K the free volume fraction of PET increases. In this temperature range, the free volum    Figure 8 shows the free volume fraction of PET in the three systems at different temperatures. It can be seen that, as the temperature rises over the range from 400 K to 440 K, the free volume fraction of PET increases. In this temperature range, the free volume the presence of water and non-aqueous solvent D5 increases the free volume. The presence of D5 significantly increases the total free volume of PET. It can also be seen from Figure 8b that the major changes in the free volume fractions of the three systems are mainly due to the cavities in the range of 1.0 Å to 4.0 Å in diameter. The contribution of the increase in the total free volume fraction is also mainly from the cavities in the range of 1.0 Å to 4.0 Å . The MD simulation results show that the non-aqueous solvent D5 has a better swelling effect on PET than water.

Interfacial Interaction Model and MD Simulation
The binding energy (EBinding) is an effective parameter to quantitatively characterize the strength of interactions between different substances. The total energy (Etotal) of the simulation system includes non-bond interaction energy (Enon-bond), valence energy (Evalence), and valence energy (cross terms) (Ecrossterm). The Enon-bond is the most important part of the total energy of the simulation system. The Enon-bond is composed of three important parts, including Van der Waals energy, electrostatic energy, and hydrogen bonding energy. However, the hydrogen bond energy calculated by the COMPASS force field is not listed separately, and it is included with other non-bond interaction terms [35]. The Etotal of the simulation system can be expressed by the following Equation (5) where ETotal represents the total energy of PET and H2O, EPET represents the total energy of PET after removing H2O, and EH 2 O represents the total energy of H2O after removing PET. In order to further explore the influence of H 2 O and D5 on the change in free volume of PET, the distribution of cavities in PET was studied. Figure 8a shows the distribution of free volume fractions in the three systems at a temperature of 413 K. The total free volume fractions FFV of PET are 37.459%, 37.560%, and 37.892%, respectively, indicating that the presence of water and non-aqueous solvent D5 increases the free volume. The presence of D5 significantly increases the total free volume of PET. It can also be seen from Figure 8b that the major changes in the free volume fractions of the three systems are mainly due to the cavities in the range of 1.0 Å to 4.0 Å in diameter. The contribution of the increase in the total free volume fraction is also mainly from the cavities in the range of 1.0 Å to 4.0 Å. The MD simulation results show that the non-aqueous solvent D5 has a better swelling effect on PET than water.

Interfacial Interaction Model and MD Simulation
The binding energy (E Binding ) is an effective parameter to quantitatively characterize the strength of interactions between different substances. The total energy (E total ) of the simulation system includes non-bond interaction energy (E non-bond ), valence energy (E valence ), and valence energy (cross terms) (E crossterm ). The E non-bond is the most important part of the total energy of the simulation system. The E non-bond is composed of three important parts, including Van der Waals energy, electrostatic energy, and hydrogen bonding energy. However, the hydrogen bond energy calculated by the COMPASS force field is not listed separately, and it is included with other non-bond interaction terms [35]. The E total of the simulation system can be expressed by the following Equation (5): The binding energy between PET and H 2 O can be expressed by the following Equation (6) [36]: where E Total represents the total energy of PET and H 2 O, E PET represents the total energy of PET after removing H 2 O, and E H2O represents the total energy of H 2 O after removing PET. The negative value of E Interaction is E Binding . The larger the value of E Binding , the stronger the interaction. Table 3 shows the binding energies of PET/D5 and PET/H 2 O at 413 K. The binding energy of PET and D5 is significantly larger than that of PET and H 2 O, which demonstrates that the interfacial interaction force between PET and D5 is greater than that between PET and H 2 O. This conclusion can also be obtained by observing their interfacial interaction models after MD simulation ( Figure 9). As shown in Figure 9a, before annealing, H 2 O molecules are distributed around the PET molecular chains, and D5 molecules also surround the PET molecular chains. After annealing, very few H 2 O molecules can enter the PET molecular chain. In contrast, more D5 molecules (Figure 9b) can enter the PET molecular chain and they completely interact with each other. The negative value of EInteraction is EBinding. The larger the value of EBinding, the stronger the interaction. Table 3 shows the binding energies of PET/D5 and PET/H2O at 413 K. The binding energy of PET and D5 is significantly larger than that of PET and H2O, which demonstrates that the interfacial interaction force between PET and D5 is greater than that between PET and H2O. This conclusion can also be obtained by observing their interfacial interaction models after MD simulation ( Figure 9). As shown in Figure 9a, before annealing, H2O molecules are distributed around the PET molecular chains, and D5 molecules also surround the PET molecular chains. After annealing, very few H2O molecules can enter the PET molecular chain. In contrast, more D5 molecules (Figure 9b) can enter the PET molecular chain and they completely interact with each other.

Simulation Analysis of Solubility of PET Molecular Chains
Solubility parameters and cohesive energy density are important for evaluating the compatibility of different materials. Solubility parameter (δ), which is defined as the square root of the cohesive energy density (CED), has been widely used to predict the compatibility of polymer blends.

Simulation Analysis of Solubility of PET Molecular Chains
Solubility parameters and cohesive energy density are important for evaluating the compatibility of different materials. Solubility parameter (δ), which is defined as the square root of the cohesive energy density (CED), has been widely used to predict the compatibility of polymer blends. δ = √ CED (7) When the δ of the two substances is relatively close, that is, when ∆δ is less than 2.05 (Cal 1/2 /cm 3/2 ), the compatibility of the two substances is better. When ∆δ is greater than 10.02 (Cal 1/2 /cm 3/2 ), the two substances are incompatible [37,38].
In this simulation study, using the Forcite module in the Material studio simulation software and the CED function, the solubility parameter value of each component was simulated and calculated. Table 4 shows the simulated solubility parameters of PET, H 2 O, and D5. The MD results are in good agreement with the experimental values [39], suggesting the reliability of the MD simulation. PET and D5 have similar solubility parameters, indicating their good compatibility.

Conclusions
In this paper, MD simulation was used to study the changes in PET microstructure in H 2 O and D5. The study found that after adding H 2 O, the glass transition temperature of PET decreases, the diffusion coefficient increases, and the free volume increases. After adding D5, the glass transition temperature is further reduced, and it is easier to dye PET. Moreover, the diffusion coefficient and free volume are larger than that with water, indicating that the swelling effect of D5 on PET is more ideal than that of water, which can be explained by the results of the binding energy and the solubility parameters. The binding energy of PET/D5 is significantly larger than that of PET/H 2 O, which demonstrates that the interfacial interaction force between PET and D5 is greater than that between PET and H 2 O. PET and D5 have similar solubility parameters, indicating their good compatibility. This study of the microstructure of PET is helpful to deepen the understanding of the mechanism involved during the process of dyeing PET with water and non-aqueous solvents. It can provide strong support for further research on dyeing PET with non-aqueous solvents.