Molecular Dynamics Study of the Swelling of Poly(methyl methacrylate) in Supercritical Carbon Dioxide

The swelling of a poly (methyl methacrylate) in supercritical carbon dioxide was studied by means of full atomistic classical molecular dynamics simulation. In order to characterize the polymer swelling, we calculated various properties related to the density, structure, and dynamics of polymer chains as a function of the simulation time, temperature, and pressure. In addition, we compared the properties of the macromolecular chains in supercritical CO2 with the properties of the corresponding bulk system at the same temperature and atmospheric pressure. It was shown that diffusion of CO2 molecules into the polymer led to a significant increase in the chain mobility and distances between them. Analysis of diffusion coefficients of CO2 molecules inside and outside the poly(methyl methacrylate) sample has shown that carbon dioxide actively interacts with the functional groups of poly (methyl methacrylate). Joint analysis of the radial distribution functions obtained from classical molecular dynamics and of the averaging interatomic distances from Car-Parrinello molecular dynamics allows us to make a conclusion about the possibility of formation of weak hydrogen bonds between the carbon dioxide oxygen atom and the hydrogen atoms of the polymer methyl groups.


Introduction
In recent years, much attention has been paid to the study of polymer-supercritical fluid systems. This is eloquently echoed by a large number of reviews devoted to the use of supercritical carbon dioxide (sc-CO 2 ) as a solvent for the physical processing of polymeric materials [1,2], as well as in polymer modification, formation of polymer composites, polymer blending, microcellular foaming, polymerization [3], particle formation [4], and applications of sc-CO 2 in the fabrication of polymer systems for drug delivery [5]. Being a non-toxic and environmentally friendly solvent, supercritical CO 2 can replace toxic organic compounds in a number of chemical processes [6][7][8], which is of great interest to the pharmaceutical industry, for example, when creating prolonged dosage forms of drugs by polymer doping [2,9]. The process of polymer impregnation is based on plasticization: when treated in a supercritical solvent, the polymer adsorbs it, which leads to an increase in the mobility of the segments and chains, and an increase in the average distance between the polymer monomers, thereby facilitating the penetration of various additives into the polymer and desorption of undesirable impurities [10,11]. It has been shown previously that CO 2 is a good plasticizer for amorphous polymers like poly (methyl methacrylate) (PMMA) [12,13]. PMMA/CO 2 systems were investigated by means of different experimental techniques in a wide range of conditions [14][15][16][17][18][19][20][21][22][23][24]. In particular, sorption and swelling of different polymers including PMMA in the presence of sub-and supercritical carbon dioxide have been measured at pressures up to 30 MPa at 308.2 and 323.2 K in reference [15]. Handa Y.P. and coauthors have studied glass transition in the system PMMA/compressed gas as a function of the gas pressure using high-pressure calorimetric measurements [16]. Sorption of CO 2 in PMMA at 308−573 K and concurrent dilation of the polymer at 308-358 K over a pressure range up to 5 MPa were studied in reference [18]. Nikitin L.N. and colleges have performed the study of sorption by poly (methyl methacrylate) and poly (butyl methacrylate) in sc-CO 2 conditions using technique of direct optical observation [19]. In reference [20], a purely gravimetric approach based on the use of a magnetic suspension balance was proposed to simultaneously measure swelling and sorption of supercritical fluids in PMMA in a commercially available setup. Pantoula M. and coauthors used the quartz crystal microbalance and the mass-loss analysis to investigate the sorption of sc-CO 2 onto PMMA and polystyrene [21], and used the magnetic suspension balance and the optical determination of the volume change to study the swelling process of this polymers [22] under pressures up to 40 MPa and temperatures T = 308-405 K. The authors of reference [25] have experimentally investigated the behavior of PMMA-based systems not only in pure sc-CO 2 , but also in sc-CO 2 modified by acetone, ethanol, and methylene chloride. They have found that processing PMMA-based polymers with pure sc-CO 2 leads to polymer swelling, and addition of a liquid cosolvent to CO 2 enhances polymer dissolution. However, selection of a cosolvent and its concentration is crucial for optimizing solubility. In particular, the polymer solubility in CO 2 acetone, as a function of cosolvent concentration, reaches a maximum at about 15 wt % of the cosolvent, while for CO 2 -ethanol, the solubility is low and practically not affected by the increase in the cosolvent percentage. In reference [26], the swelling of poly (methyl methacrylate) and poly (butyl methacrylate) bulk samples in supercritical carbon dioxide was studied in situ. The kinetics of swelling, the diffusion coefficients of CO 2 in the polymers were calculated, the effects of temperature and pressure on the obtained values were analyzed and it was concluded that the degree of PMMA swelling at a fixed exposure temperature (311 K) increased as the pressure grew and showed a nonmonotonic dependence on the temperature.
Along with experimental studies, simulation approaches are useful and powerful tools for investigating systems containing a polymer at the molecular level. Van der Vegt et al. using molecular dynamics (MD) simulations have made calculations of the sorption thermodynamics of CO 2 in a model glassy polymeric membrane [27] and studied the temperature dependence of carbon dioxide transport in an amorphous polyethylene melt [28]. CO 2 sorption and swelling in glassy matrices of atactic polystyrene obtained by coarse-graining, equilibration, and reverse-mapping were simulated over the temperatures ranging from 308 to 405 K at pressures of up to 30 MPa [29]. Combining experimental and computational techniques, Zhang et al [30] studied CO 2 -induced plasticization in a polyimide membrane. They discussed the effect of the CO 2 interactions with the ether groups on the mobility of polyimide chains, calculated the glass transition temperature of the polyimide at different CO 2 content values in the polymer matrix. In reference [31], the molecular dynamics simulation was used to investigate the adsorption of PMMA and polyvinyl acetate on anα-quartz surface and to understand the interactions between the quartz surface and the polymers. An interesting study of P. Xue et al. [32] focused on the study of the mechanism of a supercritical CO 2 thickener using MD simulations. The authors examined poly (vinyl acetate-covinyl ether) used as a thickener and showed that adding the polymer reduced the diffusion of supercritical CO 2 indicating an interaction between the solvent molecules and the polymer functional groups. In particular, it was found that the ester group had better ability to bind a CO 2 molecule than the ether group.
Although experimental studies of the PMMA swelling process in sc-CO 2 have already been reported, this phenomenon has not yet been investigated at the microscopic level. To the best of our knowledge, PMMA swelling in sc-CO 2 has never been studied in detail at the molecular level using MD simulations. It should be noted that the time and length scales of MD simulations are limited to a few nanoseconds and nanometers, and therefore it is not possible to achieve the same system size or observation time as in the experiment. Nevertheless, we can estimate the structural and dynamical changes of the polymer matrix at the molecular level at the very beginning of the swelling process.
We present here a molecular dynamics study of two types of systems: an atactic PMMA in bulk at atmospheric pressure in the temperature interval 273 K-533 K, and PMMA in sc-CO 2 at temperatures of 333 K and 353 K and pressures 10-25 MPa. We have also calculated and discussed the structural and dynamic characteristics including polymer density, specific volume, coefficient of thermal expansion, glass transition temperature, mean squared displacement, radial distribution functions, gyration radius of the chain, and end-to-end distance. Figure 1 illustrates the monomer unit of PMMA. The polymer chain was built from 100 monomers; then 27 chains were used to construct a three-dimensional structure with an initial density of 1000.0 kg/m 3 by Materials Studio [33]. As is well known, force field selection is determinant of the accuracy of simulation results of a specific system. In this study, for MD simulations of PMMA in bulk and in sc-CO 2 , the OPLSAA (Optimized Potentials for Liquid Simulations All Atom) force field was used [34,35]. Validation of this force field by comparing the PMMA structure and dynamics with neutron experiments was carried out by C. Chen and co-authors [36]. It was shown that the simulation model provides a fair description of real PMMA samples. The Lennard-Jones (LJ) and partial atomic charge parameters of the OPLSAA force field used for PMMA are given in Table 1. For carbon dioxide, we used the model developed by Z. Zhang and Z. Duan [37]. As for cross site-site interactions, we applied the geometric mean mixing rule for both LJ parameters. The temperature and pressure were controlled by a Nose'-Hoover thermostat [38,39] and a Parrinello-Rahman barostat [40], respectively. The leap-frog integrator was adopted to integrate the equations of motion [41]. The cutoff radius was set of 1.5 nm for all interactions. For the long-range electrostatic interactions, a particle mesh Ewald [42,43] with a grid spacing of 0.25 nm and an interpolation order of four was used. The constraints were implemented using the LINCS algorithm [44]. The time step was 1 fs for all the simulations. In order to gain a homogenous sample of PMMA, the polymer was equilibrated in accordance with the heating and cooling scheme as the one used in a few previous studies [45,46]. All the MD simulations were conducted using GPU-accelerated GROMACS v5.0.7 [47]. Processes of energy minimization, canonical ensemble (NVT, 0.5 ns), and isothermal-isobaric (NPT, 5 ns) ensemble were applied to the PMMA sample at 533 K and 0.1 MPa. After that, the sample was subjected to cooling (for t = 0.5 ns) starting from 533 K and ending at 273 K in ∆T = 10-20 K steps. The equilibrium configurations are performed in NPT ensemble with a pressure of 0.1 MPa and temperature of 273K-533K with 10K-20K steps. At the end of this stage, the densities of the samples and the coefficient of thermal expansion were calculated and compared with the experimental ones [48,49]. We present here a molecular dynamics study of two types of systems: an atactic PMMA in bulk at atmospheric pressure in the temperature interval 273 K-533 K, and PMMA in sc-CO2 at temperatures of 333 K and 353 K and pressures 10-25 MPa. We have also calculated and discussed the structural and dynamic characteristics including polymer density, specific volume, coefficient of thermal expansion, glass transition temperature, mean squared displacement, radial distribution functions, gyration radius of the chain, and end-to-end distance. Figure 1 illustrates the monomer unit of PMMA. The polymer chain was built from 100 monomers; then 27 chains were used to construct a three-dimensional structure with an initial density of 1000.0 kg/m 3 by Materials Studio [33]. As is well known, force field selection is determinant of the accuracy of simulation results of a specific system. In this study, for MD simulations of PMMA in bulk and in sc-CO2, the OPLSAA (Optimized Potentials for Liquid Simulations All Atom) force field was used [34,35]. Validation of this force field by comparing the PMMA structure and dynamics with neutron experiments was carried out by C. Chen and co-authors [36]. It was shown that the simulation model provides a fair description of real PMMA samples. The Lennard-Jones (LJ) and partial atomic charge parameters of the OPLSAA force field used for PMMA are given in Table 1. For carbon dioxide, we used the model developed by Z. Zhang and Z. Duan [37]. As for cross site-site interactions, we applied the geometric mean mixing rule for both LJ parameters. The temperature and pressure were controlled by a Nose'-Hoover thermostat [38,39] and a Parrinello-Rahman barostat [40], respectively. The leap-frog integrator was adopted to integrate the equations of motion [41]. The cutoff radius was set of 1.5 nm for all interactions. For the long-range electrostatic interactions, a particle mesh Ewald [42,43] with a grid spacing of 0.25 nm and an interpolation order of four was used. The constraints were implemented using the LINCS algorithm [44]. The time step was 1 fs for all the simulations. In order to gain a homogenous sample of PMMA, the polymer was equilibrated in accordance with the heating and cooling scheme as the one used in a few previous studies [45,46]. All the MD simulations were conducted using GPU-accelerated GROMACS v5.0.7 [47]. Processes of energy minimization, canonical ensemble (NVT, 0.5 ns), and isothermal-isobaric (NPT, 5 ns) ensemble were applied to the PMMA sample at 533 K and 0.1 MPa. After that, the sample was subjected to cooling (for t = 0.5 ns) starting from 533 K and ending at 273 K in ΔT = 10-20 K steps. The equilibrium configurations are performed in NPT ensemble with a pressure of 0.1 MPa and temperature of 273K-533K with 10K-20K steps. At the end of this stage, the densities of the samples and the coefficient of thermal expansion were calculated and compared with the experimental ones [48,49].    In order to investigate the PMMA swelling process, the samples of the polymer equilibrated at 333 K, 353 K, and 0.1 MPa were placed in the center of a cubic cell with periodic boundary conditions and "embedded" by previously equilibrated sc-CO 2 (85086 molecules) at certain thermodynamic parameters ( Table 2). The production run simulations were performed for 10 ns with a time step of 1 fs. The data were collected for analysis every 0.1 ps. Table 2. Stated parameters of the simulated systems: Temperature T, Pressure P, experimental Density of CO 2 ρ(CO 2 ) EXP [50], total Density of the systems at the end of simulation ρ MD , and size of cubic cell L. In addition to the classical molecular dynamics, a Car-Parrinello (CPMD) simulation of a small PMMA fragment (consisting of 3 monomer units) in sc-CO 2 has been carried out at T = 333 K and ρ = 725 kg/m 3 in the CPMD-3.13.2 [51] program package. The total Car-Parrinello dynamical system consists of two adiabatically decoupled subsystems: the cold electronic degrees of freedom and the nuclear degrees of freedom at the relevant physical temperature. The computations have been performed using the gradient-corrected BLYP functional [52,53]. The orbitals of the valence electrons were expanded in a plane wave basis set to a 25 Ry cutoff. The interaction between the core and the valence electrons was described by the ultrasoft pseudopotential in the Vanderbilt form [54]. The Brillouin zone was sampled at the Γ-point only. The fictitious electronic mass and integration step were set up to 600 a.u. and 5 a.u., respectively. The initial configuration consisting of a PMMA fragment surrounded by 58 CO 2 molecules was simulated by the classical molecular dynamics. The obtained classical trajectory was then equilibrated for 10 ps by means of the CPMD simulation. The production run length was 10 ps; the statistics were gathered every 1.2 fs. The simulation was performed in the NVT ensemble with a Nose-Hoover chain thermostat [38,39].

Bulk PMMA
In order to demonstrate that the force field and simulation parameters were chosen correctly, the polymer physical properties must be realistically represented in the simulation. As it was mentioned in the previous section, the densities of the samples and the coefficient of thermal expansion were calculated and compared with the experimental ones. The results with the theoretical and experimental values, available in the literature [48,49], are reported in Table 3. The densities and coefficients of thermal expansion, obtained for PMMA at different temperatures by means of MD simulations, are in good agreement with the experimental ones. Table 3. Comparison of properties (density ρ and coefficient of thermal expansion β) at atmospheric pressure in the temperature range 273 -533 K obtained by using MD and experimentally.
T, K ρ, kg/m 3 ρ exp , kg/m 3 β, 10 −4 1/K β exp , 10 − [49] Glass transition temperature (T g ) is a unique property of polymers. Below T g , polymers behave like glass that is hard and brittle, and above T g , polymers act like rubber that is soft and viscous [55,56]. Although polymers possess two completely different states below or above T g , glass transition is a second order phase transition, so that many first order properties, for example, volume, change gradually when the temperature increases. Therefore, T g can be derived from fitting the intercept of the two linear trend lines at a low and a high temperature, respectively. There are several ways to estimate T g from MD simulations. In particular, in reference [57] T g was determined through tracing the variations in the macroscopic (thermal conductivity, volume, thermal expansion and Young's modulus) and microscopic properties (radial distribution functions, mean squared displacement, non-bonded energy) of the polymer during temperature cooling scans. The authors [57] found that the density and volume method was less time consuming in determining the T g than the thermal conductivity and Young's modulus method. Therefore, we calculated the specific volume, i.e., the inverse density, at each temperature during the cooling process. The temperature dependence of the specific volume of the bulk PMMA is shown in Figure 2.
The temperature gradient of specific volume has a discontinuity at Tg, therefore, the intersection between the lines obtained as interpolation of specific volume values below and above Tg is the estimation of the glass transition temperature. We used this method to obtain the simulated glass transition temperature of the polymer (T g = 417 K). The glass transition temperature value obtained from the simulation was found to be a little higher than the corresponding experimental values (*T g = 378 K [58], T g = 363-387 [59]). These deviations could be caused by the extremely fast cooling rate of the MD simulation (≈10 9 K/s) relative to the real experiment (≈0.3 K/s) [59]. Moreover, as it was mentioned in reference [57], the T g obtained from MD simulations, is dependent on the polymerization degree, i.e., the higher the polymerization degree is, the higher the glass transition temperature is (T g = 450 K and T g = 381 K for polymerization degree of PMMA 100 and 10, respectively, at a cooling rate 20 K/ns). the two linear trend lines at a low and a high temperature, respectively. There are several ways to estimate Tg from MD simulations. In particular, in reference [57] Tg was determined through tracing the variations in the macroscopic (thermal conductivity, volume, thermal expansion and Young's modulus) and microscopic properties (radial distribution functions, mean squared displacement, non-bonded energy) of the polymer during temperature cooling scans. The authors [57] found that the density and volume method was less time consuming in determining the Tg than the thermal conductivity and Young's modulus method. Therefore, we calculated the specific volume, i.e., the inverse density, at each temperature during the cooling process. The temperature dependence of the specific volume of the bulk PMMA is shown in Figure 2. Glass transition of polymers is well correlated with mobility of polymer chains, which can be affected by chemical constituents and intermolecular interactions. When a polymer is heated above T g, the intermolecular interactions become weaker and, as a result, the mobility and flexibility of the chains increase. In order to examine the mobility of PMMA during the glass transition process, the mean squared displacement (MSD) was calculated. The MSD curves were calculated by the relation: where r i (t) is the position vector of atom i at time t; the symbol < . . . > denotes the average for all the atoms as well as for all the time origins.
The steeper slope of MSD indicates higher mobility of the polymer chain. As one can see in Figure 3, the slopes of MSD above T g are much higher than those below T g , showing higher mobility of the polymer chains above the glass transition temperature (namely 417 K). The temperature gradient of specific volume has a discontinuity at Tg, therefore, the intersection between the lines obtained as interpolation of specific volume values below and above Tg is the estimation of the glass transition temperature. We used this method to obtain the simulated glass transition temperature of the polymer (Tg = 417 K). The glass transition temperature value obtained from the simulation was found to be a little higher than the corresponding experimental values (*Tg = 378 K [58], Tg = 363-387 [59]). These deviations could be caused by the extremely fast cooling rate of the MD simulation (≈10 9 K/s) relative to the real experiment (≈0.3 K/s) [59]. Moreover, as it was mentioned in reference [57], the Tg obtained from MD simulations, is dependent on the polymerization degree, i.e., the higher the polymerization degree is, the higher the glass transition temperature is (Tg = 450 K and Tg = 381 K for polymerization degree of PMMA 100 and 10, respectively, at a cooling rate 20 K/ns).
Glass transition of polymers is well correlated with mobility of polymer chains, which can be affected by chemical constituents and intermolecular interactions. When a polymer is heated above Tg, the intermolecular interactions become weaker and, as a result, the mobility and flexibility of the chains increase. In order to examine the mobility of PMMA during the glass transition process, the mean squared displacement (MSD) was calculated. The MSD curves were calculated by the relation: (1) where ri(t) is the position vector of atom i at time t; the symbol <…> denotes the average for all the atoms as well as for all the time origins.
The steeper slope of MSD indicates higher mobility of the polymer chain. As one can see in Figure 3, the slopes of MSD above Tg are much higher than those below Tg, showing higher mobility of the polymer chains above the glass transition temperature (namely 417 K). In reference [57] the same behavior of temperature dependence of MSD curves of PMMA is observed. Namely, the MSD curves remain constant with variation in temperature below 470 K and MSD values notably increased with temperatures raised above 470 K. The difference between our results and data presented in reference [57] is probably due to a small number of PMMA molecules in the MD cell (only 3) which was simulated by M. Mohammadi and coauthors compared with our simulation box containing 27 molecules.
In order to understand how the structural behavior of the polymer samples depends on temperature, we calculated the radius of gyration (Rg), which is one of the most important quantities In reference [57] the same behavior of temperature dependence of MSD curves of PMMA is observed. Namely, the MSD curves remain constant with variation in temperature below 470 K and MSD values notably increased with temperatures raised above 470 K. The difference between our results and data presented in reference [57] is probably due to a small number of PMMA molecules In order to understand how the structural behavior of the polymer samples depends on temperature, we calculated the radius of gyration (R g ), which is one of the most important quantities in conformational statistics of polymer chains [60][61][62]. The radius of gyration was calculated by the following Equation: where m i is the mass of site i and r i is the position of site i relative to the center of mass of the molecule. As can be seen from Figure 4, there is an abrupt change in the slope of R g as a function of temperature around 417 K. Thus, at temperatures above 417 K, significant changes in the structure of the polymer are observed, i.e., the polymer transitions from a glassy state to a highly elastic state. It should be noted that the Tg value obtained from the temperature dependence of the specific volume of the bulk PMMA coincides with the value obtained from R g (T).
where mi is the mass of site i and ri is the position of site i relative to the center of mass of the molecule.

213
As can be seen from Figure 4, there is an abrupt change in the slope of Rg as a function of temperature 214 around 417 K. Thus, at temperatures above 417 K, significant changes in the structure of the polymer

PMMA in sc-CO 2
The first stage of dissolution of any polymer is its swelling. Swelling is the process of absorption of a low molecular weight solvent by a polymer, accompanied by an increase in the mass and volume of the polymer and a change in the conformation of its macromolecules. Because the PMMA solubility is negligible in sc-CO 2 (it is less than 0.001 percentile weight of the extracted amount per unit mass of CO 2 at 333 K and 20 MPa) [25], the dissolution process stops at the stage of swelling, and, thus, one can discuss limited swelling of PMMA in supercritical carbon dioxide. PMMA swelling in sc-CO 2 leads to an increase in the end-to-end distance R ete and radius of gyration R g ( Figure 5) in comparison with the bulk PMMA where R ete and R g fluctuate around a constant value during the entire simulation time. In order to investigate the mobility of PMMA segments, we calculated the self-diffusion coefficients of the PMMA from the slope of the MSD according to the Einstein's relation [63] as follows: As Table 4 shows, the PMMA swelling in CO2 changes the mobility of the chains. The selfdiffusion coefficient of the polymer in the supercritical solvent is three orders of magnitude higher than in the bulk.  12 20 0.85 ± 0. 12 25 0.78 ± 0.15 PMMA in sc-CO2 has comparatively higher mobility, and thus requires a lower temperature compared with the bulk PMMA to obtain the same segment mobility. The increase in mobility of the polymer molecules after adding CO2 can be understood as a plasticization effect due to the increased In order to investigate the mobility of PMMA segments, we calculated the self-diffusion coefficients of the PMMA from the slope of the MSD according to the Einstein's relation [63] as follows: As Table 4 shows, the PMMA swelling in CO 2 changes the mobility of the chains. The self-diffusion coefficient of the polymer in the supercritical solvent is three orders of magnitude higher than in the bulk. PMMA in sc-CO 2 has comparatively higher mobility, and thus requires a lower temperature compared with the bulk PMMA to obtain the same segment mobility. The increase in mobility of the polymer molecules after adding CO 2 can be understood as a plasticization effect due to the increased Materials 2019, 12, 3315 9 of 18 space between the chain segments. Figure 6 illustrates the swelling process of PMMA in sc-CO 2 at 333 K and 25 MPa. The snapshots have been made by using the visualization program VMD [64]. We have calculated the density of CO 2 and PMMA in x, y, and z directions for each nanosecond (Figure 7). Figure 7a shows dependences ρ x (r), ρ y (r), and ρ z (r) after 2, 3, 4, 5, and 10 ns. Figure 7b,c show the density profiles of PMMA versus time during polymer swelling in sc-CO 2 . As we can see, in the first four nanoseconds after the beginning of the simulation, more dramatic changes in the PMMA structure are observed. The CO 2 molecules gradually diffuse into the polymer and the difference in the solvent density at the edges of the box (i.e., outside the polymer) and closer to the center of the box (i.e., in the area where the polymer is located) decreases. The density profiles of PMMA become wider and lower (Figure 7b,c), meaning the polymer adsorbs the solvent molecules and swells.
The ability of polymers to swell is characterized by the degree of swelling (α), which is defined as the amount of solvent absorbed by the polymer, per unit mass or volume of the polymer: where m 0 and V 0 are the mass and volume of the original polymer, respectively; m and V are, respectively, the mass and volume of the swollen polymer, and m ads is the mass of the solvent adsorbed by the polymer. We try to estimate the degree of PMMA swelling by using the first method, namely, focusing on the mass. The calculation was carried out as follows: the mass of carbon dioxide molecules, located at a distance of 0.5 nm from any PMMA atom was calculated at the initial time (m 0 (CO 2 )) and every ns (m t (CO 2 )). Then we found the change in the mass of carbon dioxide near PMMA, compared to the initial value: m ads (CO 2 ) = m t (CO 2 )-m 0 (CO 2 ), i.e., the mass of the solvent adsorbed by the polymer, and then used Equation (4).  Figure 7a shows dependences ρx(r), ρy(r), and ρz(r) after 2, 3, 4, 5, and 10 ns. Figure 7b and 7c show the density profiles of PMMA versus time during polymer swelling in sc-CO2. As we can see, in the first four nanoseconds after the beginning of the simulation, more dramatic changes in the PMMA structure are observed. The CO2 molecules gradually diffuse into the polymer and the difference in the solvent density at the edges of the box (i.e., outside the polymer) and closer to the center of the box (i.e., in the area where the polymer is located) decreases. The density profiles of PMMA become wider and lower (Figure 7b,c), meaning the polymer adsorbs the solvent molecules and swells.
The ability of polymers to swell is characterized by the degree of swelling (α), which is defined as the amount of solvent absorbed by the polymer, per unit mass or volume of the polymer: where m0 and V0 are the mass and volume of the original polymer, respectively; m and V are, respectively, the mass and volume of the swollen polymer, and mads is the mass of the solvent adsorbed by the polymer. We try to estimate the degree of PMMA swelling by using the first method, Limited swelling occurs over a long time and is determined by the rate of diffusion of the solvent molecules into the polymer. In the simplest case, the swelling process proceeds as a first-order reaction, therefore the swelling rate is equal to: where k is the swelling rate constant, α max is the maximum degree of swelling. Thus, the swelling kinetics Equation has the following form: Using Equation (7), the degree of PMMA swelling obtained from MD by means of Equation (4) was described as a function of time (Figure 8). The standard deviation was about 0.99. Table 5 shows the coefficients of the swelling equation for all of the systems.

(7)
Using Equation (7), the degree of PMMA swelling obtained from MD by means of Equation (4) was described as a function of time (Figure 8). The standard deviation was about 0.99. Table 5   1.23 ± 0.01 0.47 ± 0.02 Certainly, the maximum degree of PMMA swelling in sc-CO2obtained from MD simulations is significantly different from the experimental results. For example, in reference [65][66][67] it has been shown that the equilibrium degree of swelling of PMMA in sc-CO2 can be as high as 20 wt %, and the authors of reference [26] obtained a volume degree of PMMA swelling of 32 ± 6% at 323K and 12.5 MPa and 37 ± 5% at 311 K and 12.5 MPa. In the recent work of R. Li et al. [24] the swelling ratio, which is defined as the ratio of the volume change under isobaric conditions during the swelling with the initial volume of the polymer sample at ambient temperatures, was found to be 0.8654 at 353 K and 10 MPa and the highest value of the swelling ratio 1.1910 was achieved at 363 K and 12 MPa. Such discrepancies between the MD and the experiment are the result of differences in the sample molecular weight and size: in reference [26] the molecular weight of PMMA was MW = 4098 kDa, in reference [24] the average molecular weight was 500 kDa, while in the MD simulations it was a few  Table 5. The rate (α max -the maximum degree of swelling, k -the swelling rate constant) of swelling from Equation (7) for all thermodynamic points. The numbers after ± denote the standard deviation. Certainly, the maximum degree of PMMA swelling in sc-CO 2 obtained from MD simulations is significantly different from the experimental results. For example, in reference [64][65][66] it has been shown that the equilibrium degree of swelling of PMMA in sc-CO 2 can be as high as 20 wt %, and the authors of reference [26] obtained a volume degree of PMMA swelling of 32 ± 6% at 323K and 12.5 MPa and 37 ± 5% at 311 K and 12.5 MPa. In the recent work of R. Li et al. [24] the swelling ratio, which is defined as the ratio of the volume change under isobaric conditions during the swelling with the initial volume of the polymer sample at ambient temperatures, was found to be 0.8654 at 353 K and 10 MPa and the highest value of the swelling ratio 1.1910 was achieved at 363 K and 12 MPa. Such discrepancies between the MD and the experiment are the result of differences in the sample molecular weight and size: in reference [26] the molecular weight of PMMA was M W = 4098 kDa, in reference [24] the average molecular weight was 500 kDa, while in the MD simulations it was a few orders of magnitude lower than in the experiments, M W = 10 kDa. Based on the analysis of the data presented in Table 5, it can be concluded that temperature has a positive effect on the PMMA swelling rate, but has a negative effect on the maximum swelling degree of the polymer, which is probably more highly influenced by the density of the solvent rather than by the temperature or pressure. Based on experimental study of swelling and impregnation process of PMMA in supercritical CO 2 , the authors of reference [67] have concluded that volume expansion of PMMA increases with increase in pressure and decreases with increases in temperature. Moreover, the effects of pressure and temperature on the extent of volume increase were directly related to the increase in solvent density.
The addition of a PMMA sample to sc-CO 2 decreases the diffusion of the CO 2 molecules compared to the pure fluid, which indicates that carbon dioxide actively interacts with the PMMA functional groups (Table 6). These data clearly show that the diffusion coefficients of the CO 2 molecules which are located inside the PMMA sample near the ester groups are 1.5-2 times lower than the system averaged diffusion coefficients of the CO 2 molecules, and even up to 8 times lower than in bulk supercritical carbon dioxide. Among other things, the latter fact demonstrates the polymer thickening ability in relation to the solvent. One can see that while in the pure supercritical carbon dioxide, the CO 2 diffusion coefficients are strongly dependent on the pressure (or, in other words, on the solvent density) and vary from ≈23 × 10 −5 cm 2 /s to ≈111 × 10 −5 cm 2 /s, in the mixture with the polymer, the diffusion coefficients values of the solvent molecules vary around one order (10-14) × 10 −5 cm 2 /s. CO 2 diffusion inside PMMA does not almost depend on temperature and pressure and D CO2 values vary from ≈6 × 10 −5 cm 2 /s to ≈7 × 10 −5 cm 2 /s. The latter may indicate that internal structure of PMMA samples (for example, size of cavities where CO 2 molecules are located) is similar. Through analysis of the movement of optical boundaries and the kinetics of swelling, the diffusion coefficients of CO 2 in the PMMA were calculated by Gallyamov M. O. and colleagues [26]. For instance, they found that at 311 K and 15 MPa, the diffusion coefficients of CO 2 obtained with using optical technique and from volume swelling kinetics are equal at (0.20 ± 0.03) × 10 −5 cm 2 /s and (0.07 ± 0.03) × 10 −5 cm 2 /s, respectively. Taking into account the differences in the temperature and the sample size in the experiment and our simulations, we can conclude that the diffusion coefficients from MD simulations are consistent with the experimental data. Gallyamov M. O. and colleagues [26] also noted that as the temperature increased from 311 K to 338 K, the diffusion coefficient in PMMA increased by 20%-40%. In our case when the temperature increased from 333 K to 353 K, the diffusion coefficient in PMMA-CO 2 system also increased by 20%-27% depending on the pressure. The radial distribution functions (RDFs) between the PMMA and CO 2 atoms were obtained by averaging over the last 2 ns of the trajectory. Figure  As clearly shown at Figure 9a, in the range of r < 0.5 nm, peaks on the RDFs carbon(PMMA)oxygen(CO2)/carbon(CO2) are well distinguished. These peaks determine the arrangement of solvent molecules around the polymer. The most pronounced and highest peaks on the RDFs C5-CCO2, C5-OCO2, O2-CCO2, and H6,7,8-OCO2 compared to each other allow us to conclude that the ester groups of the polymer are more solvated by carbon dioxide than the methyl groups (including C3, H3, H4, H5 atoms), as well as the carbon and hydrogen atoms of the chain (C1, C2, H1, H2). In the literature [69], it has been found that there are specific interactions between CO2 and PMMA which are most probably of a Lewis acid-base nature. However, the strength of such a specific interaction is very weak. Actually the first peak on the O2-CCO2 RDF from 0.3 to 0.5 nm could be attributed to electron donor-acceptor (EDA) interactions because the maximum is located around 0.35 nm and is within the geometric criterion (0.26 ≤ RC--O ≤ 0.43 nm) for EDA interactions used by Xu W. et al. and Saharay M. et al. [70,71]. For the O1-CCO2 RDF, the first peak is wider and shifted to 0.5 nm and the small shoulder located in the range 0.3-0.43 nm indicates low probability of O1 atom participation in the EDA interactions with the CO2 molecules.
For a more accurate description of the intermolecular interaction between CO2 and PMMA, we carried out an ab initio Car-Parrinello molecular dynamics simulation of a PMMA fragment consisting of 3 monomer units surrounded by 58 CO2 molecules. Figure 10 illustrates the mutual As clearly shown at Figure 9a, in the range of r < 0.5 nm, peaks on the RDFs carbon(PMMA)-oxygen(CO 2 )/carbon(CO 2 ) are well distinguished. These peaks determine the arrangement of solvent molecules around the polymer. The most pronounced and highest peaks on the RDFs C5-C CO2 , C5-O CO2 , O2-C CO2 , and H6,7,8-O CO2 compared to each other allow us to conclude that the ester groups of the polymer are more solvated by carbon dioxide than the methyl groups (including C3, H3, H4, H5 atoms), as well as the carbon and hydrogen atoms of the chain (C1, C2, H1, H2). In the literature [68], it has been found that there are specific interactions between CO 2 and PMMA which are most probably of a Lewis acid-base nature. However, the strength of such a specific interaction is very weak. Actually the first peak on the O2-C CO2 RDF from 0.3 to 0.5 nm could be attributed to electron donor-acceptor (EDA) interactions because the maximum is located around 0.35 nm and is within the geometric criterion (0.26 ≤ R C-O ≤ 0.43 nm) for EDA interactions used by Xu W. et al. and Saharay M. et al. [69,70]. For the O1-C CO2 RDF, the first peak is wider and shifted to 0.5 nm and the small shoulder located in the range 0.3-0.43 nm indicates low probability of O1 atom participation in the EDA interactions with the CO 2 molecules.
For a more accurate description of the intermolecular interaction between CO 2 and PMMA, we carried out an ab initio Car-Parrinello molecular dynamics simulation of a PMMA fragment consisting of 3 monomer units surrounded by 58 CO 2 molecules. Figure 10 illustrates the mutual arrangement of two CO 2 molecules interacting with the polymer fragment at the end of the simulation. Visual inspection of the CPMD trajectory shows that such a triple complex existed for ≈120 fs. It is interesting to note that the main distances between the atoms of PMMA and CO 2 are in agreement with the positions of the peaks on the respective RDFs obtained by classical MD. In Figure 10, one can see that the carbon dioxide molecules are oriented in such a way that they can interact with the hydrogen atoms of the PMMA methyl groups. Previously, P. Raveendran and S. L. Wallen [71] investigated the role of cooperative C-H· · · O hydrogen bonds (HBs) as a stabilization factor in addition to the EDA interactions between CO 2 and carbonyl group by using ab initio calculations at the second-order Møller-Plesset (MP2) [72,73] level. Despite the fact that the C-H· · · O hydrogen bonds (HBs) are very weak, they may have an important stabilizing effect. Assuming that the interaction between the PMMA O atoms and the CO 2 C atoms is the only EDA type, it can be expected that two C=O bond lengths of CO 2 should be identical. However, the asymmetry of the C = O bonds in the CO 2 molecules (1.15 Å and 1.19 Å, 1.17 Å and 1.20 Å for two CO 2 molecules, respectively) and, in addition, the shortening of the C-H bonds participating in the HBs (1.09 Å) compared with "free" C-H (1.14-1.16Å) allow us to suppose that there is interaction like an HB. Moreover, the distance between C of the methyl groups and O of CO 2 is less than the upper limit (0.4 nm) for C-H· · · O hydrogen bonding [74]. On the classical RDFs C5-O CO2 , and H6,7,8-O CO2 , the peaks are located at 0.35 and 0.28 nm, respectively, which is also within the geometric criterion of weak C-H· · · O HB. arrangement of two CO2 molecules interacting with the polymer fragment at the end of the simulation. Visual inspection of the CPMD trajectory shows that such a triple complex existed for ≈120 fs. It is interesting to note that the main distances between the atoms of PMMA and CO2 are in agreement with the positions of the peaks on the respective RDFs obtained by classical MD. In Figure  10, one can see that the carbon dioxide molecules are oriented in such a way that they can interact with the hydrogen atoms of the PMMA methyl groups. Previously, P. Raveendran and S. L. Wallen [72] investigated the role of cooperative C-H···O hydrogen bonds (HBs) as a stabilization factor in addition to the EDA interactions between CO2 and carbonyl group by using ab initio calculations at the second-order Møller-Plesset (MP2) [73,74] level. Despite the fact that the C-H ···O hydrogen bonds (HBs) are very weak, they may have an important stabilizing effect. Assuming that the interaction between the PMMA O atoms and the CO2C atoms is the only EDA type, it can be expected that two C=O bond lengths of CO2 should be identical. However, the asymmetry of the С = О bonds in the СО2 molecules (1.15 Å and 1.19 Å, 1.17 Å and 1.20 Å for two CO2 molecules, respectively) and, in addition, the shortening of the С-Н bonds participating in the HBs (1.09 Å) compared with "free" C-H (1.14-1.16Å) allow us to suppose that there is interaction like an HB. Moreover, the distance between C of the methyl groups and O of CO2 is less than the upper limit (0.4 nm) for C−H···O hydrogen bonding [75]. On the classical RDFs C5-OCO2, and H6,7,8-OCO2, the peaks are located at 0.35 and 0.28 nm, respectively, which is also within the geometric criterion of weak C-H···O HB. Figure 10. The complex of a PMMA fragment with two CO2 molecules according to the CPMD. All the distances are given in nm. The averaging of the distances was made over the last 120 fs. (The red spheres denote the oxygen atoms, the gray ones represent the hydrogen atoms, the cyan ones represent the carbon atoms.).
For PMMA, the research of J. R. Fried and W. Li [76] indicated that there are weak dipole-dipole interactions between the ester moieties in PMMA. Also, they demonstrated that carbon dioxide is capable of overcoming the internal interactions of the carbonyl groups. Our findings concerning specific interactions between CO2 and PMMA confirm that they are responsible for the observed plasticization of the polymer.

Conclusions
The simulations performed on bulk PMMA at atmospheric pressure have been used to validate the model and the force field employed in the classical molecular dynamics simulations. The simulations were able to satisfactorily reproduce PMMA's experimental density, glass transition temperature and coefficient of thermal expansion. An equilibrated PMMA sample was then placed in a sc-CO2 medium in different thermodynamic states. In comparison with bulk PMMA, we Figure 10. The complex of a PMMA fragment with two CO 2 molecules according to the CPMD. All the distances are given in nm. The averaging of the distances was made over the last 120 fs. (The red spheres denote the oxygen atoms, the gray ones represent the hydrogen atoms, the cyan ones represent the carbon atoms).
For PMMA, the research of J. R. Fried and W. Li [75] indicated that there are weak dipole-dipole interactions between the ester moieties in PMMA. Also, they demonstrated that carbon dioxide is capable of overcoming the internal interactions of the carbonyl groups. Our findings concerning specific interactions between CO 2 and PMMA confirm that they are responsible for the observed plasticization of the polymer.

Conclusions
The simulations performed on bulk PMMA at atmospheric pressure have been used to validate the model and the force field employed in the classical molecular dynamics simulations. The simulations were able to satisfactorily reproduce PMMA's experimental density, glass transition temperature and coefficient of thermal expansion. An equilibrated PMMA sample was then placed in a sc-CO 2 medium in different thermodynamic states. In comparison with bulk PMMA, we observed an increase in the end-to-end distance and radius of gyration, following PMMA swelling in sc-CO 2 . Moreover, in the sc-CO 2 medium, the chains mobility and the distance between them were significantly increased. The kinetics of polymer swelling has also been studied. The results have shown that the PMMA swelling degree is overestimated compared to the experimental values due to the limitations inherent in MD simulations, which are related to the size of the polymer sample. Nevertheless, the present calculations demonstrate that the combined use of the classical and ab initio MD provides a good way to better understand the intermolecular interactions in the polymer/fluid system. In particular, the PMMA intermolecular interactions with the solvent molecules are not only of the electron donor-acceptor type but also weak C-H· · · O hydrogen bonds. The results of MD indicate that sc-CO 2 could be a desirable swelling agent in the impregnation of PMMA with additives and, therefore, we will make the molecular mechanism of the diffusion of small organic molecules into the polymer matrix in supercritical media the subject of further publications.