Polyacrylamide in glycerol solutions from an atomistic perspective of the energetics, structure, and dynamics

Polyacrylamide (PAM) is a biocompatible and water-soluble polymer that can be tailored to meet a broad range of industrial applications, most of them based on its well above room temperature, glass transition temperature of 400 K.^1,2 The polymer is synthesized either as a simple linear chain or as a cross-linked structure. PAM increases the viscosity of water and belongs to the super water-absorbent polymer (SAP) family. When hydrated, PAM forms a soft gel used in gel electrophoresis for protein separation. Indeed, PAM is hydrophilic and can form aqueous solutions of very high concentrations.³ Because of their gel-like properties, these aqueous solutions are employed as flocculants in the removal of suspended particles from sewage and industrial effluents such as paper mill wastewater. Through the highly reactive amide NH₂ groups, the polymer can be chemically modified to produce cationic or anionic polymers, which are particularly useful in mineral-processing and metallurgical operations for the separation of metals from residues.⁴ PAM increases the viscosity of other fluids and may be the cause of unexpected turbulent behavior in the flow of otherwise viscoelastic fluids at low Reynolds numbers.⁵ 

In this article, we will focus on the behavior of the structural fate of PAM oligomers as solutes in pure glycerol and water solvents and in a 90:10 mixture of these two liquids. The approach is an analysis at the all-atom modeling of the components using Molecular Dynamics (MD) for inspection of the system thermodynamic and energetic properties. The past couple of decades have seen a dramatic increase in computational power and high performance computer algorithms among which MD simulations have emerged as valuable tools for studying macromolecules and large molecular systems at the atomic scale. From MD simulations, it is possible to follow in time the interfacial dynamics of complex 3D molecular structures both localized around particular macromolecules and interacting with molecular liquids, which are not yet possible to be observed experimentally.⁶ For example, using a coarse grained force field (FF) and MD,⁷ the effect on the PAM oligomers’ structure by surfactant molecules in aqueous solutions was investigated. The study revealed that short oligomers in water curled into compact clusters in the absence of surfactant molecules contrasting with the stretched-out conformations at the hydrophilic interface created by the surfactant molecules. Another investigation⁸ demonstrated that the addition of PAM oligomers at the interface of water and a foam system increased the foam stability. This study used MD and the all-atom COMPASS FF. More recently, de Oliveira and co-workers⁹ conducted extensive all-atom MD simulations of N-propylacrylamide solvated in water for inspecting the impact of copolymerization with acrylamide on the lower critical transition temperature (LCST). In a combined experimental–computational work,¹⁰ a PAM oligomer with 40 monomers solvated in alcohols, water, and their mixture was simulated with all-atom MD with the objective of tuning the upper critical transition temperature (UCST) upon the alcohol–water relative concentration.

The atomic scale behavior of a large number of liquids has been the subject of several decades of discovery, which spans from the time that MD became a clearly useful method for studying the dynamics and structure of systems in the fluid phases¹¹ to the current research of more complex liquids. For example, in Refs. 9 and 10, several alcohols were MD simulated, both pure and mixed with water. In another study,¹² a new FF was developed for ethyl acetate and its mixture with water. Propane-1,2,3-triol or glycerol is a sugar alcohol with three hydroxyl groups, which was MD stimulated in aqueous solutions.¹³ These authors provide an extensive literature review of MD simulations of glycerol. Chelli et al.¹⁴ carried out the first glycerol, all-atom, MD simulation using the Generalized Amber Force Field (GAFF).^15,16 The latter was re-parameterized with different atomic charges and Lennard-Jones parameters,¹⁷ yielding a glycerol self-diffusion coefficient that agreed better with the experimental results. The CHARMM FF was proposed to model glycerol;^18,19 however, the published simulations give MD values for the diffusion coefficient about half of the experimental values. More recently, the AMBER force field of Ref. 17 with modifications for the glycerol bending and torsion angle constants²⁰ gave rise to reasonable MD-simulated densities of glycerol, although diffusivities were low when compared to experiments. A detailed MD analysis of several glycerol–water mixtures²⁰ pointed to the formation of local pockets of water and of glycerol induced by the hydrogen bonding network occurring in both liquids. Jahn et al.²¹ conducted a comparative study of glycerol MD calculated densities and thermodynamic properties including five FFs, the GAFF-AMBER,^14,17 the CHARMM,¹⁸ and three versions of the OPLS.²² This study indicated that the use of the AMBER FF resulted in reasonable thermodynamic properties. Over the years, the AMBER package¹⁵ has had periodical improvements to their FFs, which, currently, enables users for custom-inclusion of two types of atomic charges to the GAFF.

Our study is motivated, in part, by an experiment⁵ that investigates the flow of a viscoelastic fluid, (90:10) glycerol–water containing diluted PAM, in microchannels at low Reynolds numbers. In that study, it was observed that, downstream, the flow became unstable, displaying the features of elastic turbulence. In that experimental setup, as the polymers flow initially around cylindrical obstacles, the PAM is thought to elongate when flowing along the cylinders, while the turbulent flow structures might be associated with sudden coiling–stretching events due to fluctuations in the stream-wise velocity gradients. A recent fluid dynamics continuum simulation²³ corroborated that vortices can be produced in viscoelastic fluids if there is a spatiotemporal localization of energy in the neighborhood of flowing polymers modeled by FENE-P dumbbells.²⁴ 

In this work, we present extensive all-atom MD simulations of PAM oligomers of various molecular weights solvated in pure glycerol, pure water, and mixed glycerol–water (90:10) by weight. The goal is to investigate the structural fate of the solvated oligomers in the three chosen solvents. Our simulation protocol is explained in Sec. II. The first modeling scenario considered (Sec. III A) involves the structural fate of 40 PAM oligomers containing an even number of monomers (between 10 and 50) in the three implicit solvents. In the second modeling scenario (Secs. III B and III C), the structural changes of a selected set of PAM oligomers containing 10, 20, and 30 monomers are analyzed when solvated in explicit solvents, modeled at the all-atom level. To achieve the latter, a new parameterization of the GAFF was carried out for pure glycerol and the mixed glycerol–water solution. The results concerning the structural properties of the three selected PAM oligomers in the explicit solvents are summarized in Sec. III C, including the radius of gyration, hydrodynamic radius, end-to-end distance, Z orientation order parameter, and interaction energies between the oligomers and the solvents. Summarizing observations are cast in the discussion of Sec. IV. We predict that PAM oligomers of all sizes under ambient conditions collapse into a compact globule in implicit solvents (modeled as a continuum) but remain as randomly coiled oligomers of variable elongations in the three all-atom explicit solvents considered. We additionally predict that the conformations of PAM oligomers considered here are basically trapped in either of the glycerol-based solvents and their structure is maintained without visible changes over tens of nanoseconds. This paper is concluded in Sec. V.

Three solvents are considered, pure glycerol, pure water, and their mixture at 90:10 concentration in molecular weight (M_w). The solutes in each of these solvents are polyacrylamide oligomers of various lengths, n-PAM, formed by n acrylamide monomers [—CH₂CH(CONH₂)—]. The polymer backbone has two carbon atoms, one of them bonded to the amide group (NH₂)—C=O. All generated polymers are syndiotactic with the amide groups arranged alternatively on both sides of the chain backbone. Two scenarios are considered. First, the focus is on the atomic charges for the GAFF modeling of n-PAM oligomers, followed by a survey of the oligomers’ energetic and structural behavior in implicit solvents under ambient conditions. Second, the focus is on the effect that the three considered solvents modeled within the all-atom explicit solvent approach have a subset of the n-PAM oligomers: 10-PAM (712.8 u), 20-PAM (1423.6 u), and 30-PAM (2134.4 u). Additionally. the atomic charges for GAFF are custom-generated for the glycerol and glycerol–water solvents. The MD AMBER package¹⁵ is used in all simulations. The simulation strategy is an adaptation of our more general MD modeling process for macromolecular systems.²⁵ In what follows, we adopt the terminology oligomer for the n-PAM because these polymeric chains have M_w < 10 ku.²⁶ 

The fractional atomic charges for n-PAM containing even numbers of monomers between 10 and 50 (Mw 710.8 u–3555.9 u) are calculated within the restrained electrostatic potential (RESP)²⁷ approach. These atomic charges are generated for each n-PAM as a full molecule at the B3LYP 6-31G^* density functional theory (DFT) level, including the polarizable continuum model (PCM),²⁸ based on the Merz–Singh–Kollman population analysis,^29,30 and using Gaussian09.³¹ The atomic charges are later ported into the AMBER Tools16¹⁵ to generate their corresponding RESP values. Additionally, for a selected set of n-PAM oligomers, an alternative set of atomic charges was generated at the BCC (additive bond charge corrections) approach.³² The latter are based on semiempirical AM1 quantum approaches and are calculated within the AMBERTools16 routines.

Concerning the all-atom MD simulations of pure glycerol and the mixed 90:10 glycerol–water solvents, we used our modified GAFF with the two types of atomic charges. For water, the SPC-E model³³ is adopted.

Employing the newly generated FF, a battery of 50 ns NVE-MD simulations were performed for 20 n-PAM lengths (n = 10–50) in implicit solvents (generalized Born solvation model³⁴) with ε = 42.48 for pure glycerol and ε = 48.66 for 90:10 glycerol–water.³⁵ These systems were started from the fully elongated oligomers, first thermalized at T = 298 K during 1 ns followed by 50 ns NVE MD equilibration, for finally collecting data during production runs of 50 ns. Five consecutive 50 ns production simulations were carried out. At the end of each run, the RESP atomic charges were recalculated, as described in Sec. II A. This round robin strategy was designed for adopting charges associated with systems that had dynamically modified initial conformations.

Along the MD trajectories, the energetics and several structural properties of the n-PAM are examined such as the radius of gyration $Rg2=∑i=1Nri−rcm2/N$⁠, where r_i is the position of atom i referred to the oligomer center of mass and N is the number of atoms in the oligomer. A useful property related to the polymer flexibility is the end-to-end distance $Ree2=∑i=1Nbb∑j=1Nbbli⋅lj$⁠, where l_i are vectors along the backbone bonds and N_bb is the number of backbone bonds.³⁶ Also useful is the orientation order parameter Z that describes the distribution of orientation angles ϕ_i between vectors connecting contiguous monomers with respect to the polymer director vector joining the end monomers, $Z=32(n−1)∑i=1n−1cos2ϕi−12$⁠.^37,38 A fully stretched-out polymer chain has Z = 1, while Z = 0 indicates random orientations with respect to the director vector. The oligomers’ principal moments of inertia I_A, I_B, and I_C obtained from diagonalizing the molecular inertia tensor I are calculated. Moments of inertia are useful for providing the rotating shape identity of molecular systems.

Glycerol is in the liquid phase between 291 K and 563 K and is often used mixed with water in a large variety of relative concentrations. Here, we validate that the adopted modeling of the pure glycerol solvent and its 90:10 mixture with water is acceptable. These systems were equilibrated with NPT-MD at T = 298 K and 101.32 Pa via the Berendsen thermostat and barostat³⁹ along 40 ns long trajectories using a 2 fs time step, a 20 Å cutoff, and periodic boundary conditions (PBCs). Ewald sums are used in all calculations for the long-range electrostatics within the particle mesh implementation (PME). This equilibration stage is followed by a 20 ns NVE-MD equilibration run, followed by a 10 ns production run for each system at the equilibrium densities and temperatures of 298 ± 2 K. As described in Secs. III and IV, the NVE simulations are used for the calculation of the pure solvents’ radial distribution functions and their self-diffusion coefficients from

where r_i is the position of the ith molecule center of mass at time t and N is the number of molecules in the solvent. Each NVE run is split into m time series, each series starting from a reference position $ri(t0k)$⁠, and their average is taken as indicated in Eq. (1). The last term is the correction due to the PBC,^40,$DPBC=2.837297kBT6πηL$⁠, with k_B being Boltzmann’s constant, T being the temperature, L being the computational box length, and η being the solvent viscosity. The values for the latter are taken from experiments at 298 K: η_glycerol = 945 mpa s,^41,η_water = 0.8937 mpa s,⁴² and η_90:10 = 0.977 mpa s.⁴³ 

The next step is the preparation of systems with one n-PAM, n = 10, 20, 30, solvated into each of the three different solvents. Each of the nine systems is brought to equilibrium with NPT-MD at 298 K and 101.325 Pa along 100 ns. The equilibrium densities of the solvated systems are within the standard deviation of the solvent equilibrated densities for the three PAM sizes considered. These simulations are followed by 25 ns NVE-MD production runs at temperatures around 298 K. It is from these NVE simulations that the solvated n-PAM structural properties (as defined in Sec. II B) are calculated, including their diffusion coefficients in the solution. For the latter, the center of mass of the full n-PAM is followed in time within the solution and calculated from Eq. (1).

The n-PAM properties presented below result from a 298 K thermalization Langevin-MD⁴⁴ lasting 1 ns–5 ns depending on the size, followed by a production NVE-MD of 50 ns for every n-PAM size in either the implicit solvent corresponding to glycerol or 90:10 glycerol–water. Figure 1 shows the averages of the potential energy per monomer, the radius of gyration R_g, and the moments of inertia along the principal axes of n-PAM as a function of the number of monomers n. The Z order parameter was in all cases very close to zero. Within the analyzed size range, n-PAM potential energy per monomer is not a smooth plateauing property for either solvent modeled implicitly, as shown in Fig. 1(a). Additionally, the n-PAM energetics is basically the same in either of the two solvents, where oligomers with n = 10, 22, 28, 32, 38, 40, 48 stand out as the least stable oligomers, while those with n = 14, 18, 26, 42, 44, 46, 50 are relatively the most stable sizes.

Concerning the radius of gyration at T = 298 K, all n-PAM sizes have collapsed with radii ranging between 6 Å and 11 Å for n between 10 and 50 when solvated in glycerol, as shown in Fig. 1(b). The effect of water in the 90:10 mixed solvent does not alter significantly this behavior, also illustrated in Fig. 1(b). Analysis of the moments of inertia indicate that n-PAM oligomers are prolate tops with a tendency to be slightly more spherical with increasing size, as shown in Fig. 1(c).

In the 1980s, it was determined that when a system has fractal characteristics, its mass does not fill densely the 3D space, and consequently, its volume $Rdf$ is characterized by a power law of the system size R less than 3. The exponent d_f is referred to as the fractal dimension⁴⁵ such that the mass M of the system is expressed by $M∝Rdf$⁠. For solvated polymers, this type of scaling exists,³⁶ where the size of the system is represented by its radius of gyration R_g and its mass by the molecular weight M_w as $Rg=κMwα$⁠, where κ is a proportionality constant. The latter implies that for α = 0.5, the fractal dimension is 2, and the space occupied by the polymer in 3D is equivalent to a surface. A “random coiled” polymer has an exponent α = 0.5. Exponents α > 0.5 identify polymers that fill the 3D space more sparsely than the random coiled polymer, while polymers with α < 0.5 identify polymers that fill the space more compactly. Figure 2 is a log–log plot of the n-PAM radii of gyration as a function of their number of monomers, clearly evidencing two size regimes, n ≤ 20 and n > 20. The slopes of the straight lines in Fig. 2 are 0.4 and 0.33 corresponding to the exponent α of oligomers with n ≤ 20 and n > 20, respectively. This result indicates that in the implicit solvent, the n-PAM with n ≤ 20 are less densely packed than the larger oligomers that fill densely the 3D space. Figure 3 illustrates the molecular structure of several n-PAM sizes. These cluster-like structures are not spherical; in all cases, the ratio R_g/R_hyd fluctuates between 0.43 and 0.57, while a value of 0.77 is accepted to characterize spherical proteins. Therefore, these compact structures are prolate globules.

The polymer structure fate in implicit solvents is somehow expected, since the chains also collapse from elongated to compact structures in vacuo. Evidently, the implicit solvent environment does not modify significantly the collapsed structure of the isolated n-PAM chains at ambient temperature. We also tested the AMBER ff14SB with both RESP and BCC atomic charges. This FF is better suited for proteins and biomolecules, and the AMBER manual¹⁵ states that the generalized Born solvation model as implemented in the software package is optimal for the various existing ffXXSBs. For PAM oligomers, however, repeating the full set of simulations with ff14SB and an implicit solvent yielded basically the same energetics and structural behavior as described above for the GAFF with both types of atomic charges.

An additional evaluation of the newly parameterized GAFF is done by correlating the potential energy of the 20 MD-generated geometries with their corresponding DFT energies, which gives a reasonable linear regression with correlation coefficients above 70%.

Based on the analysis in this section, the n-PAM conformations of the last MD simulation and their calculated RESP atomic charges are adopted as the starting initial values for oligomers of size n = 10, 20, 30 in extensive explicit solvent simulations described in Sec. III C. The BCC atomic charges for these three selected PAM oligomers are calculated for the same molecular geometry as the newly determined RESP counterparts.

Pure solvent simulations are performed for quantifying the appropriate behavior of the solvents at 298 K. The simulations of pure glycerol contain 2000 glycerol molecules, while for the mixed 90:10 glycerol:water system, the simulations involved 1800 glycerol molecules and 1018 water molecules. Meanwhile, the pure water system contains 9500 water molecules. From the NPT-MD simulations, the solvent systems attain equilibrium densities of 1258 ± 2 kg/m³ for glycerol and 1227 ± 2 kg/m³ for mixed glycerol–water with RESP atomic charges. With BCC atomic charges, the attained equilibrium densities are 1265 ± 2 kg/m³ for glycerol and 1241 ± 2 kg/m³ for glycerol–water. The calculated densities are in excellent agreement with the experiments at 298 K, 1257.91 kg/m³ for glycerol and 1253 kg/m³ for the mixed glycerol–water solvent.²⁰ The SPC-E water equilibrium density is 998 ± 2 kg/m³, consistent with the previous simulations.⁴⁶ 

At these equilibrium densities and 298 K, the radial distribution function, rdf or g(r), of each solvent is calculated, showing a very good agreement with the experimental results. Figure 4(a) shows the rdf of glycerol between atom pairs in different molecules. The inset in this figure provides the label given to the glycerol molecule atoms entering in the considered pairs between molecules. Our calculated peak positions for the six atom pairs depicted are 1.89 Å, 1.89 Å, 2.85 Å, 2.85 Å, 2.85 Å, 2.85 Å (RESP case) and 1.83 Å, 1.83 Å, 2.79 Å, 2.79 Å, 2.79 Å, 2.79 Å (BCC case), which are in excellent agreement with the experimental values.⁴⁷ As expected, the BCC-based calculation (dashed lines) yields peaks at slightly shorter distances than the RESP case due to the system higher density. For example, the O—H peak is at 1.89 Å for the RESP case and at 1.83 Å for the BCC case. Otherwise, the glycerol liquid structure is unchanged if the glycerol atomic charges are RESP or BCC modeled. Figure 4(b) depicts the rdf of the O—H pairs between glycerol–water molecules, showing that the hydrogen bonds formed with glycerol oxygens are slightly longer than the bonds formed with water oxygens. This figure includes the rdf from distances between the glycerol molecules’ centers of mass and the water oxygens (red), which shows that the two liquids form a hydrogen-bonded network (first peak) with the second peak associated with a first coordination shell of glycerol–water molecules at 2.8 Å and a second coordination shell at 3.5 Å. The SPC/E rdf is known to reproduce well the structure of liquid water,^33,48 and our results are consistent with the previous simulations.

The self-diffusion coefficients of glycerol and water are estimated from Eq. (1), considering 40 different time origins, each of 0.5 ns of NVE time evolution. The PBC corrected self-diffusion coefficient from the simulations at 298 ± 1 K are (3.35 ± 0.03) × 10⁻⁷ cm²/s for the RESP case and (1.96 ± 0.02) × 10⁻⁷ cm²/s for the BCC case. As expected, the diffusivity is lower in the BCC case because the liquid glycerol density is higher. These diffusivities compare well with the experimental value of 1.7 × 10⁻⁷ cm²/s obtained from the Taylor dispersion method.⁴⁹ However, this experimental value is larger than the diffusion coefficient obtained from the nuclear magnetic resonance (NMR) pulsed magnetic field gradient⁵⁰ or modulated gradient spin echo method.⁵¹ The corrected self-diffusion coefficient of water from the simulation at 300 ± 1 K is 2.860 × 10⁻⁵ cm²/s, which compares well with the experiment at 298 K of 2.299 × 10⁻⁵ cm²/s.⁵² 

The density of glycerol as a function of temperature is investigated as well. Figure 5 shows this behavior for the two types of charges as compared with experiments. The agreement of the two FFs is reasonable around 300 K. At higher temperatures, both FFs yield a lower density than the experimental value. We assess that the FF is modeling well the glycerol solutions in the temperature range adequate for the goals of this work.

All simulations are started from the initial n-PAM geometries and FF atomic charges, as described in Sec. III A, with the initial R_g values of 5.1 Å, 6.8 Å, and 9.4 Å for n = 10, 20, and 30, respectively. These solutes in the glycerol solutions have dilution compositions by mass of 0.38%, 0.77%, and 1.16%, while in water, the dilution is 0.42%, 0.84%, and 1.26% for n-PAM of increasing n.

Once the PAM solutes in glycerol are equilibrated at 298 K, the system acquires slightly different equilibrium densities of 1294 ± 3 kg/m³ (RESP) and 1268 ± 2 kg/m³ (BCC), which are slightly higher than the glycerol density. Similarly, the density of n-PAM solvated in the other two solvents is also slightly higher than the solvent densities. An energetics evaluation of the oligomers in the various solvents yields results, which are provided in Table I. The oligomer potential energy averages per monomer over the last 10 ns of MD-NVE runs showing for the RESP case an increasing stabilization with increasing oligomer size are given in Table I. However, the BCC case yields E_PAM/monomer basically equal for all the solvated oligomers in the three solvents, which is indicative that the electrostatic terms of the FF with BBC atomic charges are not efficient in tracking the effect of the solvent on the oligomers. The reported error corresponds to the standard deviation. The interaction energy E_int between each n-PAM and the solvents represents the balance between the total potential energy of the system and the sum of the individually separated potential energies of the solvent and the oligomer: E_int = E_full-system − (E_solvent + E_PAM). As observed in Table I, the trend of the interaction energy in each of the three solvents is for the system to become less cohesive as the size of the oligomer increases. These energies present large fluctuations in the mixed glycerol–water system attributed to the mobility of the water molecules in the neighborhood of the oligomers that results in frequent changes in the solution surrounding the solutes. Comparison between the interaction energies of n-PAM in the different solvents depends on the system sizes, with a decreasing PAM stabilizing propensity as the number of water molecules increase in the solvent.

Analysis of the various MD trajectories indicate that the structural fingerprints of the solvated oligomers are acquired during the NPT MD equilibration runs. On the other hand, along the subsequent 40 ns NVE runs, those fingerprints are basically locked-in and do not change substantially with time. The calculated polymer structural properties include gyration and hydrodynamic radii, end-to-end distance, and orientational order parameter. These properties are very similar between the different solvents, although there are evident fluctuations. Indeed, analysis of these properties’ time evolution reveals that the random coiled conformations remain as such for longer periods of time in the glycerol and glycerol–water systems than in water. These properties are not always distributed normally; therefore, calculation of averages has little meaning. Figure 6 illustrates the R_g distribution of the three oligomers modeled with RESP (full line) and BCC (dashed line) when solvated in the three solvents considered. The RESP case results in narrow Rg distributions, portraying oligomers that are unable to change their conformation rapidly, which contrasts with the broad and multipeaked distribution of Rg in the BBC case. The latter is another evidence that the BCC atomic charges are not as well suited for representing interactions between linear molecules and molecular solvents such as glycerol.

The n-PAM oligomers are massive, compared to the solvent molecules. Despite this fact, we additionally evaluated the oligomers’ diffusivity in each system at 298 ± 1 K. Therefore, the motion of the oligomer center of mass is followed in time. Using Eq. (1), our estimates of n-PAM diffusion coefficients in glycerol are 0.41, 0.30, and 0.20 ± 0.01 × 10⁻⁸ cm²/s in the RESP case and 2.97, 1.83, and 1.59 ± 0.01 × 10⁻⁸ cm²/s in the BCC case for 10-, 20-, and 30-PAM, respectively. In the mixed glycerol–water solvent, the diffusivity increases approximately by a factor of 4 yielding diffusion coefficients of 1.72, 0.98, and 0.93 ± 0.01 × 10⁻⁸ cm²/s (RESP case) and 8.58, 6.20, and 5.41 ± 0.01 × 10⁻⁸ cm²/s (BCC case) for the three oligomers with n = 10, 20, and 30, respectively. These estimates indicate that the oligomers diffuse about ten times less than the glycerol molecules where they are solvated. Finally, within water as a solvent, the oligomers’ diffusion coefficients increase by about 200 times, and the difference between the two atomic charges is strongly reduced with values 460, 306, and 237 ± 1 × 10⁻⁸ cm²/s (RESP) and 432, 286, and 267 ± 1 × 10⁻⁸ cm²/s (BCC) for 10-, 20-, and 30-PAM, respectively.

PAM of high molecular weights is part of the family of thermo-responsive polymers, and we analyzed a battery of simulations for evaluating the glycerol solvated oligomers’ sensitivity induced by temperature changes. In one of them, the runs are initiated from fully stretched-out chains, equilibrated first at a high temperature of 500 K and finally cooled down in steps of 10 K until reaching 298 K. This cooling process was carried out with MD NPT. At each cooling step of the ladder, the system was run for 20 ns trajectories, while the high temperature and final temperature runs were 40 ns long. The in silico cooling rate is 0.8 K/ns. Another alternative explored consisted in starting the simulation at 298 K from the extended oligomers, and let the system run for hundreds of nanoseconds. In yet another alternative setup, the solvated extended oligomers were heated up in 5 K steps from 298 K to about 350 K and cooled back down. What was learned from these attempts is the enormous variability in conformations with different R_g that these n-PAM oligomers can take when they become trapped in the solvent. For example, in the RESP case and pure glycerol, 30-PAM ended these processes with R_g ≈ 9 Å–10 Å, which contrasts with other attempts in which the final R_g ≈ 12 Å–14 Å. Not only these outcomes for the 30-PAM radii are significantly different between them but also the distributions are peaked at different values than that shown in Fig. 6. A commonality across simulations, however, is that the RESP case gives rise to vey compact R_g distributions when compared to the BCC equivalents. In the mixed glycerol–water system, a similar situation is encountered. We conjecture that the size of these oligomers is short for having a well-defined lower critical transition temperature (LCST) in which the polymer changes abruptly from the coiled conformation to collapsed globule. Therefore, when cooling between 350 K and 300 K, the chain can be in any intermediate configuration, and by lowering the temperature, the glycerol solvents trap the instantaneous geometry, and the fate of the chain is to be locked-in an instantaneous structure. When solvated in water, these oligomers are flexible, and they do not acquire a locked-in R_g at 298 K.

The end-to-end distance R_e−e is another interesting structural property commonly associated with flexibility. In this respect and in the glycerol solvents, we observe a marked difference between the RESP and BCC cases. Compounding the observations, we conclude that the glycerol molecules have the ability to keep the oligomer ends fairly localized for extended periods of 50 ns or more while simultaneously promoting global torsional motions of the macromolecule that changes its overall size. In all the trial simulations and analyzed temperatures, the orientational order parameter Z fluctuates around zero. Therefore, the oligomer structure has a random distribution of angles with respect to the oligomer chain orientation vector, irrespective of being a random coil or a collapsed globule.

Attempting a scaling as described in Sec. III A is elusive due to the ability of these oligomers to lock-in more than one value of the R_g. To compensate for that lack, we conducted an equivalent battery of simulations for two other oligomer sizes, 16-PAM and 24-PAM. Considering averages between the various simulations and based on only five n-PAM sizes, we predict the α exponent in $Rg∝Mwα$ to be approximately 0.7, 0.55, and 0.5 for oligomers in glycerol, mixed glycerol–water, and water, respectively. These values contrast with the smaller α exponents obtained for the implicit solvent (Sec. III A), implying that as more water is present in the glycerol solvents, the more compact the oligomer structure becomes. In turn, these exponent values validate that glycerol is a good solvent for PAM and, as water is added, the solvent approaches a θ solvent. We also predict that, from the simulations carried out at 298 K, the n-PAM are not as compact as the structures obtained from the implicit solvent. Considering the scaling laws of polymer conformations for 30-PAM, from $Rgcoil/Rgglobule∼n/ℓν−1/3$⁠, the revised Flory exponent ν = 0.588,⁵⁴ and the persistence length ℓ ∼ 2.5, n = 30 ∼ 12ℓ, which suggests $Rgcoil/Rgglobule=1.83$⁠. From our simulations, we have $Rgcoil/Rgglobule=1.8$⁠. Therefore, even though adventurous, we predict that PAM of any length in the implicit solvent defines the compact globular size, while the simulation-based R_g of PAM in all-atom explicit solvents is identified with the randomly coiled chains of variable length.

Additionally, efforts were carried out for identifying solvent molecules near the n-PAM. From several radial distribution functions, on the average, we assert that glycerol molecules in the neighborhood of 30-PAM form a broad coordination shell of about 5.3 Å width and a minimum 2.5 Å approach of any glycerol molecule to the oligomer monomers. As water is added to the glycerol solvent, a small amount of water clumps, forming a thin coordination shell at 2.7 Å.

In this work, we have presented a dynamical modeling of the structure and energetics of PAM oligomers in glycerol solvents and in water with the goal of elucidating the sensitivity of the PAM structure to viscous non-Newtonian liquids such as glycerol and its mixture with water at high concentrations. In order to achieve this goal, we have modified the atomic charges of the GAFF and introduced a set of RESP atomic charges that effectively permit a strong temporal localization of the polymer chain in these glycerol solvents. We have also modified the established GAFF atomic charges of the glycerol molecules by the RESP values, obtaining a very good representation of the glycerol liquid at ambient temperature. It is the combination of these two modeling strategies that enables the PAM oligomers to acquire locked-in, swollen, and elongated structures in the glycerol solvent, while they are more flexible and prone to remain as less elongated random coils in the glycerol–water solvent. In water, the oligomers are very flexible, changing frequently their random coil structure without fully collapsing into a compact globule. Both the glycerol–water and the pure water solvents behave as θ solvents for the 10-, 20-, and 30-PAM oligomers considered here. In contrast, the BCC atomic charges’ case fails to provide a clear and distinct behavior of the solvated oligomers in the different solvents. Furthermore, we predict that the compact globular conformations of PAM oligomers are acquired in the implicit solvent MD, and these conformations set the limit of maximum compactness for the PAM oligomers. This is the first simulation of its type for PAM oligomers solvated in glycerol solvents. Our protocol and simulation strategy is portable and will be useful when applied to other modeling environments.

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

We acknowledge partial support from the Commonwealth of Virginia 4-VA grant and from the George Mason University (MDR Seed Grant No. 215001). R.A.H. acknowledges the partial support of the National Science Foundation (Grant No. 1904953). Simulations were performed on ARGO, the computing cluster provided by the Office of Research Computing, George Mason University.