Article pubs.acs.org/JPCB
Cite This: J. Phys. Chem. B 2019, 123, 5316−5323
Alcohol Clustering Mechanisms in Supercritical Carbon Dioxide Using Pulsed-Field Gradient, Diffusion NMR and Network Analysis: Feedback on Stepwise Self-Association Models Trent R. Graham,† Daniel J. Pope,‡ Yasaman Ghadar,‡ Sue Clark,‡ Aurora Clark,*,†,‡ and Steven R. Saunders*,†,‡ The Gene and Linda Violand School of Chemical Engineering and Bioengineering and ‡Department of Chemistry, Washington State University, Pullman, Washington 99164, United States
Downloaded via BUFFALO STATE on July 25, 2019 at 18:59:26 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.
†
S Supporting Information *
ABSTRACT: Co-solvent clustering in complex fluids is fundamental to solution phase processes, influencing speciation, reactivity, and transport. Herein, methanol (MeOH) clustering in supercritical carbon dioxide is explored with pulsed-field gradient, diffusionordered nuclear magnetic resonance spectroscopy (DOSY-NMR), and molecular dynamics (MD) simulations. Refinements on the application of self-association models to DOSY-NMR experiments on clustering species are presented. Network analysis of MD simulations reveals an elevated stability of cyclic tetrameric clusters across MeOH concentrations, which is consistent with experimental DOSY-NMR molecular cluster distributions calculated with selfassociation models that include both cooperative cluster assembly and entropic penalties for the formation of large clusters. Simulations also detail the emergence of cluster-assembly and cluster-disassembly reactions that deviate from stepwise monomer addition or removal. This combination of experiment, simulation, and novel analyses facilitates refinement of models that describe co-solvent aggregation with far-reaching impact on the prediction of solution phase properties of complex fluids.
1. INTRODUCTION The addition of co-solvents to supercritical carbon dioxide (sCO2) improves the efficiency of solution phase processes such as metal extraction.1−5 In this case, simple alcohols (e.g., MeOH) form molecular clusters6−12 that enhance the solubility of chelating ligands and metal−ligand complexes.13,14 The behavior of co-solvent clusters and their size distribution are fundamentally important in sCO2 chemistry. Nuclear magnetic resonance spectroscopy (NMR) has proven useful in the description of the behavior and size distribution of cosolvent clusters.6,7,12,15,16 Established strategies for evaluating alcohol cluster distributions include the use of the ensemble hydroxyl chemical shift, which is measured as a mass average due to fast proton exchange between clusters compared to the timescale of the NMR measurement.12 The hydroxyl chemical shift of alcohol clusters is sensitive to hydrogen bonding and exhibits a cluster size-dependent NMR shielding tensor, which can be estimated using theory-based NMR chemical shift calculations.6,16,17 In addition to techniques involving NMR spin relaxation,18 cluster distributions are also investigated via ensemble diffusometry with pulsed-field gradient, diffusion-ordered NMR spectroscopy (DOSY-NMR).12,15,19 The acquired diffusion coefficients are often analyzed to describe clustering through the definition of an average aggregation number, yet many studies do not determine the distribution of cluster sizes.20−22 In addition, accounting for the size dependence of molecular diffusion is a challenge because solvent clustering © 2019 American Chemical Society
impacts the diffusion of all chemical species within the sCO2 solution.15,23,24 Typical analyses to evaluate the distribution of clusters utilize the apparent hydrodynamic radius to approximate aggregation numbers volumetrically with two assumptions: (1) the size-dependent transition from stick to slip boundary conditions at molecular scales is independent of cluster size, and (2) the apparent radius calculated from the ensemble diffusion coefficient via Stokes−Einstein relationship is equivalent to the mass-average radius.15 In this work, these two assumptions are investigated in detail. DOSY-NMR and molecular dynamics simulations, in tandem, enumerate MeOH cluster distributions in sCO2 at 19.7 MPa and 323 K. Using the cluster distributions observed in complementary molecular dynamics (MD) simulation, the stepwise self-association models applied to DOSY-NMR data are varied in complexity to include cooperativity in the formation of small clusters and entropic penalties for the formation of large clusters. The presence of relatively stable, cyclic, and tetrameric MeOH clusters is observed in MD simulation, and further analysis reveals alternative cluster formation and dissolution pathways in percolating mixtures of MeOH and sCO2 that deviate from stepwise reactions of monomers. Received: March 11, 2019 Revised: May 9, 2019 Published: June 3, 2019 5316
DOI: 10.1021/acs.jpcb.9b02305 J. Phys. Chem. B 2019, 123, 5316−5323
Article
The Journal of Physical Chemistry B
2. EXPERIMENTAL AND THEORETICAL METHODS 2.1. Materials. Molecular sieves (3 Å, 1−2 mm beads) were obtained from Alfa Aesar and flame-dried before use. NMR-quality tetramethylsilane (TMS, Si(CH3)4, 99.9+%) was obtained from ACROS Organics and stored at 5 °C. Methanol (HPLC grade, 0.2 μm filtered) was purchased from Fisher Chemicals and stored under molecular sieves in a desiccator. Deuterated methanol (D, 99.8%) was purchased from Cambridge Isotope Laboratories, Inc. and stored under molecular sieves in a desiccator. Carbon dioxide (99.9+%) was purchased from NorLab. 2.2. Sample Preparation. Samples were prepared in a large-volume, ceramic, high-pressure Daedalus Innovations NMR cell equipped with an integrated valve. Mixtures of protonated and deuterated methanol and a trace amount of 5 °C tetramethylsilane were added to the high-pressure NMR tube. The high-pressure NMR cell was immediately assembled to manufacturer’s specification, connected to a manifold containing a temperature and pressure probe in close proximity to the NMR cell, and purged with CO2 gas. The manifold and NMR cell were then heated in a heated water bath to 50 ± 0.1 °C. The cell was then pressurized to 19.7 ± 0.02 MPa and allowed to equilibrate for 6 h. The standard deviation on the sample temperature and pressure describes the deviation between the measured temperature and pressure proximal to the NMR cell after closing the NMR cell valve for the measured samples. 2.3. NMR Acquisition Parameters. NMR spectra were acquired on a Varian-Inova nuclear magnetic resonance spectrometer at 11.7 T at a calibrated temperature of 50 °C.25 The spectrometer was equipped with a Nalorac HF dualtuned probe with Z-axis gradients. NMR spectra were acquired 6 h after insertion to minimize convection. For onedimensional 1H NMR experiments, 32 transients were collected. The delay between experiments (d1) was set to 25 s. Concentrations of MeOH were calculated using the qNMR toolset in vNMRJ (v4.2), which utilizes a probe-specific external calibration to relate the integrated signal intensity of the methyl protons of MeOH to absolute concentration. Further information regarding the qNMR toolset and implementation of the program in vNMRJ is available.26 For mixtures of MeOH and deuterated MeOD, the concentration of MeOH was determined gravimetrically. For DOSY-NMR experiments, a bipolar pulse pair with stimulated echo and convection compensation experiment (Dbppste_cc) was conducted at 19.7 MPa and 50 °C. Sixteen transients were recorded to produce each spectrum. The delay between experiments (d1) was set to 25 s. The diffusionencoding gradient pulse length was optimized for each sample and varied between 0.3 and 1 ms. The diffusion gradient amplitude varied linearly in 15 steps between 2.96 and 59.91 G·cm−1. For each sample, the delay between diffusionencoding gradient pulses was set to 50 ms. Diffusivities were calculated using a nonuniform gradient (NUG)-corrected spline model (SPLMOD) from the attenuation of the integrated CH3 resonance as a function of gradient strength in vNMRJ (v4.2). Typical 1H NMR and 1H DOSY-NMR spectra are included in the Supporting Information. 2.4. Simulation Protocol. The analysis of the cluster distributions determined from DOSY-NMR experiments was complimented by classical molecular dynamics (MD) simulations. Simulation boxes with different mole fractions of
MeOH and sCO2 were constructed using the Packmol program27 that span 6.4, 12.5, 18.8, 25, and 65 mol % MeOH. The number of CO2 and MeOH molecules within each system is shown in Table S1. Each simulation box was initially constructed of dimension Lx = Ly = Lz = 46 Å that obeys periodic boundary conditions in all three dimensions and then was brought to mechanical and thermal equilibrium at 323 K using alternating NPT and NVT runs followed by a production run using the NVT ensemble of Nose and Hoover.28−30 For each system, multiple runs of 1 ns were performed to confirm that equilibration was reached. The production run was performed for 2 ns, with 1.5 ns used for cluster formation analysis. The MD simulations employed EPM231 and OPLS32 benchmarked force fields for MeOH and sCO2, respectively.9 The equilibrated clusters have average densities and pressures in good accord with the experimental conditions, as presented in Table S2. Clustering was defined by the presence of a hydrogen bond between MeOH hydroxyl groups using a geometric definition of 2.55 Å ascertained from the first minima of the radial distribution between oxygen and hydroxyl hydrogens of MeOH as shown in the Supporting Information and further justified by the consistent angle distributions shown in Figure S5. The obtained cluster distributions were also tested against other force fields employed within prior literature studies (see the Supporting Information).8,9,31,32 The hydrogen bond network created by aggregated MeOH molecules was analyzed using the ChemNetworks software program33 with results shown in the Supporting Information. All simulations were performed using the LAMMPS software (version 16 February 2016)34 with the Verlet algorithm35 and a 1 fs time step. The long-range electrostatic interactions were modeled by Ewald summation in K-space with a cutoff of 11 Å and an accuracy of 1 × 10−6 for determining the root-mean-squared (RMS) error when the forces are being evaluated. This reference value was chosen as a representative of the magnitude of electrostatic forces in atomic systems and is the default in LAMMPS. Thus, an accuracy value of 1 × 10−6 means that the RMS error will be a factor of 1,000,000 smaller than the reference force. The accuracy setting is used in conjunction with the pairwise cutoff of 11 Å to determine the number of K-space vectors. For further information regarding the estimation of the RMS force errors in real space for Ewald summation, we refer the reader to ref 36 and the LAMMPS documentation for the default parameters employed in this work. Additional details regarding the analysis of MD data to determine cluster distributions, the distribution of hydrogen bond angles, comparison of diffusion coefficients from the experiment and simulation, and the generation of transition matrices for quantifying association and dissociation pathways of MeOH clusters are provided in the Supporting Information.
3. RESULTS AND DISCUSSION While the Stokes−Einstein relationship was developed for systems in which the diffusing object is larger than the solvent, systematic studies of tracer diffusion have generated empirical models to correct the Stokes−Einstein relationship at the nanoscale.23 DOSY-NMR experiments of clustering compounds often determine cluster distributions through hydrodynamic radii calculations with the modified Stokes−Einstein relationship15 5317
DOI: 10.1021/acs.jpcb.9b02305 J. Phys. Chem. B 2019, 123, 5316−5323
Article
The Journal of Physical Chemistry B Da =
kBT caπηR a
applied to eq 6, resulting in an expression for the apparent hydrodynamic radius ij n α yz R a = jjjjca∑ i zzzz j i = 1 ciR i z k {
(1)
−1
where Da is the mass-average diffusion coefficient of the clusters, kB is the Boltzmann constant, T is the temperature, Ra is the hydrodynamic radius of a sphere exhibiting the massaverage diffusion of the clusters, ca is an empirical friction coefficient of a sphere that exhibits the mass-average diffusion of the clusters,23 and η is the viscosity of the mixture. With some notable exceptions,37−40 indefinite self-association models are applied to DOSY-NMR experiments with the limiting assumption that the apparent hydrodynamic radius is equivalent to the mass-average radius.15 In the Supporting Information, we detail the application of self-assembly models to determine cluster distributions with this assumption. As an alternative, the apparent hydrodynamic radius can instead be considered the radius of a sphere exhibiting the cluster massaverage diffusion. We clarify that the apparent hydrodynamic radius is not equivalent to the mass-average hydrodynamic radius, with an example case where a cluster distribution is described by an isodesmic, stepwise self-association model presented in eq 2.
where ci is an empirical friction coefficient of a cluster of size i,23 and Ri is the hydrodynamic radius of a cluster of i molecules. In contrast, the mass-average radius of the MeOH clusters, Rm, is given by n
Rm =
n i=1
(3)
n i=1
(4)
where CT is the total analytical concentration of MeOH. The kinetics of cluster assembly and disassembly are rapid on the NMR timescale, which leads to the acquisition of mass-average diffusion coefficients. The cluster distribution in terms of the mass fraction of MeOH clusters composed of i molecules (αi) is given as αi =
iCi iK i − 1(C1)i = n CT ∑ j = 1 (jK j − 1(C1) j )
(5)
The mass-average diffusion Da is related to the cluster distribution via n
Da =
∑ αiDi i=1
i=1
n
∑ j = 1 (jK j − 1(C1) j )
(8)
R a3 R13
,where Nv is the aggregation
number calculated with the limiting assumption) is explored for a wide range of cluster distributions corresponding with various solvent and solute sizes. Importantly, the error becomes more pronounced for cluster distributions in which the monomer is the same size or larger than the solvent and for systems exhibiting extensive clustering, as is the case for the current system consisting of sCO2/MeOH. To address such issues, this manuscript goes beyond prior studies utilizing DOSY-NMR to investigate clustering. Specifically, the hydrodynamic radii and aggregation numbers are directly calculated from the ensemble diffusion coefficient with consideration of (1) the size-dependent transition from stick to slip boundary conditions at molecular scales for discrete cluster sizes,23 (2) the apparent radius calculated via the Stokes−Einstein relationship from the ensemble diffusion coefficient being not equivalent to the mass-average radius, and (3) the reported indistinguishability43 between self-assembly models being significantly reduced by comparison with complimentary theoretical simulations. Based on the information obtained from MD simulation, including hydrogen bonding network analysis that results in a detailed depiction of cluster distributions and reaction mechanisms (vide infra), the DOSY-NMR data were analyzed using four stepwise self-association models that were independently applied to investigate cluster distributions and the effects of nonideal self-association. Nonideal, stepwise selfassociation models allow for the equilibrium constant to vary depending on the size of the cluster. These self-association models likewise only describe clusters that are generated from the addition or removal of monomers, as shown in eq 9. While it is true that equilibrium coefficients can be construed such that eq 9 can describe the thermodynamics and stability of clusters including many-body, cluster−cluster reaction pathways (coagulation and fragmentation), here, the equation is limited to describe only the binary stepwise addition and subtraction of monomers. Simple mathematical relationships are used to reduce the dimensionality of the equilibrium coefficients in self-assembly models.42 The equilibrium coefficients were determined from the DOSY-NMR exper-
where Ci is the molar concentration of MeOH clusters that are composed of i molecules, and C1 is the molar concentration of MeOH monomers that are not in a cluster. The clusters must sum to a finite analytical concentration
∑ iCi = ∑ iK i− 1(C1)i
∑
i=1
iK i − 1(C1)i R i
methodology (i.e., Nv =
where K is the equilibrium coefficient, Ai is a cluster of size i, and i is an integer between 1 and a maximum cluster size of n. In general, these stepwise self-association models are steadystate solutions of discrete Smoluchowski equations simplified to restrict cluster reaction pathways to those involving monomers.41 In this example, the assumption of a constant equilibrium coefficient for the addition of a monomer to any sized cluster is given by the isodesmic, stepwise self-association model.42 The concentration of each cluster with this assumption is given as
CT =
n
∑ αiR i =
We suggest that the apparent radius of the MeOH clusters given in eq 7 corresponds to the radius of a sphere that would exhibit the mass-average diffusivity of the distribution of MeOH clusters. In contrast, the radius found in eq 8 is the mass-average radius of the clusters. In the Supporting Information, the corresponding error in calculating the aggregation number based on the previous volumetric ratio
(2)
Ci = KCi − 1C1 = K i − 1(C1)i
−1
(7)
K
Ai − 1 + A1 ↔ Ai
ij n yz iK i − 1(C1)i j zz z = jjjca∑ n j−1 j z jj ciR i ∑ j = 1 (jK (C1) ) zz i = 1 k {
(6)
where Di is the diffusion coefficient of a cluster composed of i monomers. The Stokes−Einstein relationship (eq 1) is then 5318
DOI: 10.1021/acs.jpcb.9b02305 J. Phys. Chem. B 2019, 123, 5316−5323
Article
The Journal of Physical Chemistry B imental data with an iterative nonlinear least-squares regression.44 Experimental results and additional details regarding the use of internal tracer molecules to evaluate the viscosity of the mixture and implementation of the Stokes− Einstein relationship to DOSY-NMR data to generate cluster distributions are presented in the Supporting Information. Ki
Ai − 1 + A1 ↔ Ai
(9)
The four self-association models systematically vary the equilibrium coefficient governing the stepwise self-association of clusters to approximate complex phenomena such as cooperative self-assembly and regions of stable cluster sizes. In the first approach, the isodesmic K model describes ideal self-association and assumes that there are no differences in cluster stability.42 K i = KISO
(10)
In the case of the increasing K model, the clusters are assumed to exhibit cooperativity in self-association.42 Ki =
(i − 1)KINC i
(11)
Figure 1. (A) Apparent hydrodynamic radius of MeOH in sCO2 at 19.7 MPa and 323 K where the line is the regressed isodesmic equilibrium coefficient. Semitransparent lines are included to show effects of a 15% deviation in the regressed equilibrium coefficient, and error bars on the data points correspond to uncertainty arising from the standard error of diffusion obtained by fitting the Stejskal−Tanner relationship46,47 in vNMRJ (v. 4.2). (B) Cluster size-dependent equilibrium coefficients (Ki) for the stepwise self-association models.
The third model is based on a decreasing K, where cluster stability decreases at large cluster sizes so as to approximate the effect of entropic penalty for large cluster formation.45 Ki =
iKDEC i−1
(12)
Finally, in the piecewise K model, cooperative self-association at small cluster sizes and entropic penalties at larger cluster sizes produce a local maximum in cluster stability at a cluster size of 4. l (i − 1)K o 2K , A o o , 2≤i
