Microscopic Structure and Solubility Predictions of ... - ACS Publications

Subscriber access provided by Fudan University

Article

Microscopic Structure and Solubility Predictions of Multi-Functional Solids in Supercritical Carbon Dioxide: A Molecular Simulation Study Javad Noroozi, and Andrew S Paluch J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b12390 • Publication Date (Web): 23 Jan 2017 Downloaded from http://pubs.acs.org on January 23, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Microscopic Structure and Solubility Predictions of Multi-Functional Solids in Supercritical Carbon Dioxide: A Molecular Simulation Study Javad Noroozi† and Andrew S. Paluch∗,‡ Department of Chemical and Petroleum Engineering, Sharif University of Technology, Tehran, Iran, and Department of Chemical, Paper and Biomedical Engineering, Miami University, Oxford, Ohio 45056, USA E-mail: [email protected] Phone: (513) 529-0784. Fax: (513) 529-0761

∗

To whom correspondence should be addressed Sharif University ‡ Miami University †

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Molecular dynamics simulations were employed to both estimate the solubility of the nonelectrolyte solids acetanilide, acetaminophen, phenacetin, methylparaben and lidocaine in supercritical carbon dioxide, and to understand the underlying molecular-level driving forces. The solubility calculations involve the estimation of the solute’s limiting activity coefficient which may be computed using conventional staged free energy calculations. For the case of lidocaine wherein the infinite dilution approximation is not appropriate, we demonstrate how the activity coefficient at finite concentrations may be estimated without additional effort using the dilute solution approximation, and how this may be used to further understand the solvation process. Combining with experimental pure solid properties, namely, the normal melting point and enthalpy of fusion, solubilities were estimated. The results are in good quantitative agreement with available experimental data, suggesting that molecular simulations may be a powerful tool for understanding supercritical processes and the design of carbon dioxide-philic molecular systems. Structural analyses were performed to shed light on the microscopic details of the solvation of different functional groups by carbon dioxide and the observed solubility trends.

Introduction Tunable solvency and the ability to offer both liquid- and vapor-like properties make supercritical fluids (SCFs) a versatile solvent in a variety of applications including bioactive compound synthesis, chemical reactions, and material processing 1–4 . Carbon dioxide (CO2 ), in particular, is the most widely used supercritical solvent because of its environmentally friendly nature, abundance, low toxicity and mild critical parameters (Tc ≈304 K and pc ≈7.4 MPa). Due to these key attributes, numerous biomolecular separation and formulation processes have been developed to produce ultrafine particles which enhances the bioavailability of the solvated compounds 5,6 . Knowledge of the phase behavior of the system of interest (or specifically here the equilibrium solubility) is of central importance in all of these processes as it controls the morphology and 2 ACS Paragon Plus Environment

Page 2 of 50

Page 3 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


size of the particles. As a result, there is an impressive amount of studies reporting solubility in SCFs 4 . In addition to the need to install high-pressure equipment, the experimental determination of solubilities in SCFs is an expensive and laborious task, as one needs to measure very low concentrations 7 ( up to 10−7 ) and usually there are large discrepancies in reported solubilities by different researchers and laboratories 8,9 . Currently, most supercritical CO2 (sc-CO2 ) solubility predictions rely on empirical and semiempirical density-based correlations 10,11 which require extensive experimental solubility data to train and often are of limited transferability. Moreover, such approaches are unable to provide molecular-level insight into the underlying mechanism of the solvation phenomena in sc-CO2 , which could help drive rational design processes. As a consequence, the development of predictive approaches relying on less experimental input would be of great value for both minimizing experimental attempts and optimizing the solvency of supercritical fluids prior to conducting extensive experimental work. In this context, molecular simulation has proven to be very effective in reproducing the experimentally observed solubility trends and providing a molecular level explanation for the behavior of the system 12–16 . Microscopic information derived from simulations can help one to understand the interaction of CO2 with various functional groups, identifying the general type of the molecular fragments solvated by sc-CO2 , thereby facilitating the design of CO2 -philic materials. However, the central challenge here is the development of a molecular potential model which is sufficiently detailed to capture the role of different competing inter-molecular forces, at the same time, requiring a reasonable amount of computational resources. Once the soundness of the force field parameters are established the model can be readily extended to mixtures without additional efforts, for instance, for the screening of a best suited co-solvent for the extraction of a target solute. Compared to conventional organic liquids, poor solvency of sc-CO2 is one of the main limitations of utilizing CO2 as an industrial solvent for polar compounds, and usually co-solvents/entrainers are employed to enhance the solubility 17–19 . Recently, Forolov et al. computed the co-solvent induced solubility enhancement of a series of mono-functional organic compounds in sc-CO2 using



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

molecular dynamics (MD) simulations to explain how co-solvents modulate solubility 20,21 . Such studies provide guidance for co-solvent selection and are crucial for improving supercritical fluid processes. While molecular simulation is widely used to investigate thermodynamic and transport properties of binary sc-CO2 mixtures 22,23 , its application to solid-liquid equilibrium is less routine 24–30 . The main reasons that hinder the application of molecular simulation to phase equilibria calculations involving crystalline solids include the lack of robust computational algorithms, limitations of the accuracy of conventional force fields for crystalline thermodynamic calculations, and crystallographic data for a solute of interest may not be readily available. To address the issue of molecular model accuracy for solid phase simulations, Albo and Müller 31 probed the dissolution of naphthalene in sc-CO2 using direct contact molecular dynamics simulations. Initially, they found that the adopted model for naphthalene was able to well capture the solution phase behavior in sc-CO2 , as well as the pure component vapor-liquid equilibrium properties. On the other hand, the model performed poorly when used for solid phase property calculations. Albo and Müller 31 therefore showed that the excellent performance of the molecular model for modeling pure component vapor-liquid equilibrium and the solution phase behavior of the solute is not enough to ensure accurate solid-liquid equilibrium calculations, and care must therefore be taken when applying the model parameters for solid phase property estimation. Similar observation have been made by other researchers (see for example ref. 32). As a result, previous studies to predict solid solubility in sc-CO2 generally avoid direct simulation of the solid phase and instead couple molecular simulations in the solution phase (i.e, the calculation of the solute’s Henry’s law constant) with the experimental vapor pressure of the solid solute to estimate solubility 12,13,20,21,33 . There is, however, a great difficulty in the experimental determination of vapor pressure of low-volatile compounds as it deals with very low pressures, and to obtain directly measurable vapor pressures, high temperatures are often required which may result in the thermal decomposition of the compound of interest. To obviate the need for sublimation pressure, one may instead estimate the (Lewis-Randall normalized) activity coefficient of the


Page 4 of 50

Page 5 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


solute using molecular simulation and combine this with a limited number of pure solute fusion properties, that is the melting point and fusion enthalpy, to make solubility predictions 34–36 . Compared to sublimation pressure, fusion properties are far more simple to measure experimentally and are available for a wide range of compounds, or may be predicted with reasonable accuracy using group-contribution methods 37 . The challenge here therefore remains to accurately account for both solute-solvent and solutesolute interactions via a proper force field parameterization. To the best of our knowledge, there have been no previous studies on the prediction of activity coefficients of complex molecular solids in sc-CO2 , and for the first time, the current study attempts to investigate the accuracy-level of a standard all-atom force field for this end. To accomplish this task, the solubility of five solid solutes, acetanilide, acetaminophen, phenacetin, methylparaben and lidocaine (see fig. 1) in sc-CO2 was computed. The studied solutes contain various functional groups, and some of the solutes share similar functional groups allowing one to estimate the influence of adding or removing a functional groups on solubilization. Lastly, structural analyses through radial, spatial and combined angular-radial distribution functions were performed to investigate the solvation structure of the functional groups by sc-CO2 and to provide a molecular level explanations for the observed solubilities.

Methodology For a non-electrolyte solid solute (component 2) dissolved in a pure solvent (component 1), the equilibrium solubility is given by the classical equations of phase equilibrium thermodynamics 38,39

ln x2 = − ln γ2 (T, p, x2 ) + ln

f2S (T, p) f20 (T, p)

(1)

where x2 and γ2 are the solute mole fraction and Lewis-Randall normalized activity coefficient, T is the temperature, p is the pressure, and f2S and f20 are the fugacity of the pure solid solute and pure subcooled liquid solute. 5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 50

The term f2S /f20 is a pure component property of the solute. Assuming the solid and subcooled liquid phases are incompressible, it can be related to the value at a reference pressure (p0 ), which we will take here to be atmospheric pressure of 0.1 MPa, using the standard Poynting correction 38 . This leads to [ (( )] ) v2S − v20 (p − p0 ) f2S (T, p) f2S (T, p0 ) f S (T, p0 ) ln 0 = ln 0 exp ≈ ln 20 f2 (T, p) f2 (T, p0 ) RT f2 (T, p0 )

(2)

where v2S − v20 is the difference in the pure solid and subcooled liquid solute molar volume, and R is the molar gas constant. Assuming the change in molar volume is small, we obtain the final expression. The term f2S /f20 at atmospheric pressure (p0 ) may readily be estimated using standard fusion properties. We begin by using the definition of fugacity 38,39

ln

f2S (T, p0 ) 1 S 1 0 1 S 1 0 1 = µ2 (T, p0 )− µ2 (T, p0 ) = g2 (T, p0 )− g2 (T, p0 ) = − ∆g2fus (T, p0 ) 0 f2 (T, p0 ) RT RT RT RT RT (3)

where µS2 , g2S and µ02 , g20 are the chemical potential and molar Gibbs free energy of component 2 in the pure solid and (subcooled) liquid states, respectively, and ∆g2fus is the molar Gibbs free energy change of transferring the pure solute from its crystalline form to a (subcooled) liquid state at the same temperature and pressure. The temperature derivative of ∆g2fus / (RT ) may be related to the 38,39 enthalpy change of fusion, ∆hfus 2 , via the Gibbs-Helmholtz relation

(

∂∆g2fus (T, p0 ) / (RT ) ∂T

) =− p

∆hfus 2 (T, p0 ) RT 2

(4)

Integrating eq. (4) from the normal melting point, Tm , to T , and assuming that the difference in isobaric heat capacity between the (subcooled) liquid and solid (∆Cpfus = Cp0 − CpS ) is constant over the temperature range of interest yields


Page 7 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


f S (T, p0 ) ∆hfus 2 (T, p0 ) ln 20 = f2 (T, p0 ) RTm

( ) ( ) ( ) ∆Cpfus (T, p0 ) ∆Cpfus (T, p0 ) Tm Tm Tm − ln − 1− 1− T R T R T (5)

As noted in refs. 38,39, in eq. (5), the first term is often dominant, while the second and third terms have opposite signs and tend to cancel each other. Therefore, if the difference between the heat capacities of the solute in the solid and liquid state is small, the following expression is commonly adopted for phase equilibrium calculations, and is used in the present study 40 f S (T, p0 ) ∆hfus 2 (T, p0 ) ln 20 = f2 (T, p0 ) RTm

( ) Tm 1− T

(6)

We may readily calculate f2S /f20 using limited experimental fusion data for the pure solute at atmospheric pressure, or it may be predicted with reasonable accuracy using the group-contribution method of Hukkerikar et al. 37 . Focus here will therefore be on computing the challenging, composition dependent activity coefficient using molecular simulation.

Activity Coefficient Evaluation of eq. (1) requires that we calculate the activity coefficient of the solute (γ2 ), which is defined as 38

ln γ2 (T, p, x2 ) = ln

f2 (T, p, x2 ) − ln x2 f20 (T, p)

(7)

where f2 is the fugacity of the solute in solution at the same conditions as γ2 , and f20 is the fugacity of the pure (subcooled) liquid solute at the same T and p. Using conventional free energy calculations in the isothermal-isobaric (N pT ) ensemble, one may compute the dimensionless residual 41 chemical potential of the solute in solution, βµres 2 , defined as

( ) ig (T, p, N , N ) = βµ (T, p, N , N ) − βµ T, ⟨V ⟩ , N βµres 1 2 2 1 2 2 2 2 T,p,N1 ,N2 −1 7 ACS Paragon Plus Environment

(8)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 50

where N1 and N2 are the number of molecules of component 1 and 2, respectively, in the molecular simulation, β −1 = kB T where kB is Boltzmann’s (or the molecular gas) constant, µ2 is the chemig ical potential of component 2 in solution at the same conditions as µres 2 , and µ2 is the chemical

potential of component 2 in a hypothetical, non-interacting, ideal gas state at the same T , but at a fixed density. The relevant volume is given by the ensemble average volume of the system in the absence of the solute molecule that is being coupled/decoupled to the system, ⟨V ⟩T,p,N1 ,N2 −1 . By definition, the change in dimensionless chemical potential of component 2 between two states at the same T is equivalent to the log ratio of fugacity at the same conditions, where the fugacity of an ideal gas is equal to pressure, which may be computed using the ideal gas equation of state

βµres 2 (T, p, N1 , N2 ) = ln = ln

f2ig

( (

f2 (T, p, N1 , N2 ) T, ⟨V ⟩T,p,N1 ,N2 −1 , N2 f2 (T, p, N1 , N2 )

pig T, ⟨V ⟩T,p,N1 ,N2 −1 , N2 2

)

)

f2 (T, p, N1 , N2 ) = ln N k T 2 B /⟨V ⟩ T,p,N1 ,N2 −1

(9)

⟨V ⟩T,p,N1 ,N2 −1 N2 kB T ⟨V ⟩T,p,N1 ,N2 −1 = ln f2 (T, p, N1 , N2 ) + ln x2 (N1 + N2 ) kB T

= ln f2 (T, p, N1 , N2 ) + ln

which leads to the expression

ln f2 (T, p, N1 , N2 ) = βµres 2 (T, p, N1 , N2 ) − ln

⟨V ⟩T,p,N1 ,N2 −1 + ln x2 (N1 + N2 ) kB T

(10)

An equivalent expression can also be written for the case of performing free energy calculations for component 2 in a solution of pure component 2 (N1 = 0)

ln f20 (T, p, N2 ) = βµres,pure (T, p, N2 ) − ln 2

⟨V ⟩T,p,N2 −1 N2 kB T

and if we further assume that the number of molecules is large, then 8 ACS Paragon Plus Environment

(11)

Page 9 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


⟨V ⟩T,p,N2 −1 ⟨V ⟩T,p,N2 ≈ = v2 (T, p) N2 /Navo N2 /Navo

(12)

where v2 is the molar volume of pure (subcooled) liquid component 2, and Navo is Avogadro’s number. For substances that are solid at the temperature of interest, extreme care must to be taken when performing simulations at subcooled conditions as direct calculations may provide erroneous results. As our group has suggested previously for these cases 42 , evaluation of eq. (11) requires performing a series of simulations at temperatures above the normal melting point and then extrapolating to the subcooled liquid conditions. In our previous work, the fugacity of the pure (subcooled) liquid was computed at atmospheric pressure (p0 = 0.1 MPa). Assuming the liquid is incompressible, the fugacity at the pressure of interest may be computed using the standard Poynting correction 38

ln f20 (T, p) = ln f20 (T, p0 ) +

v2 (p − p0 ) RT

(13)

Infinite Dilution Approximation For many solutes of interest, it is reasonable to assume that the concentration is sufficiently small so as to be considered infinitely dilute. The infinite dilution limit (x2 → 0) is achieved in a molecular simulation when N2 = 1 and N1 is large. For this case we can make the simplifying assumption ⟨V ⟩T,p,N1 ⟨V ⟩T,p,N1 ≈ = v1 (T, p) (N1 + 1) /Navo N1 /Navo

(14)

where v1 is the molar volume of the pure solvent. This leads to the final expression for the infinite dilution (or limiting) activity coefficient of the solute (component 2) of

ln γ2∞ (T, p) = βµres,∞ (T, p, N1 , N2 = 1) − ln 2 =

[βµres,∞ 2

(T, p, N1 , N2 = 1) −

v1 (T, p) − ln f20 (T, p) RT

βµres,pure 2

v2 (T, p) (T, p, N2 )] + ln v1 (T, p)


(15)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 50

Dilute Solution Approximation Unfortunately, as the solute concentration increases, the infinite dilution approximation quickly breaks down. For example, the simple two-suffix Margules equation takes the form 38 : ln γ2 = x21 ln γ2∞ . For the relatively dilute concentration of x2 = 0.1, ln γ2 = 0.81 ln γ2∞ or γ2 = (γ2∞ )0.81 , which would result in a large change in the predicted solubility via eq. (1). The infinite dilution approximation assumes x2 → 0 and that solute-solute interactions are negligible. If we relax our assumptions to only solute-solute interactions are negligible, we can adopt the dilute solution approximation of Paluch and Maginn 43,44 . We can write 38

ln γ2 (T, p, x2 ) = ln γ2∞ (T, p) + [ln γ2 (T, p, x2 ) − ln γ2∞ (T, p)]

(16)

= ln γ2∞ (T, p) + ln γ2∗ (T, p, x2 ) where γ2∗ is the unsymmetric activity coefficient of component 2 normalized with respect to an ideal dilute solution or Henry’s law standard state. Following the work of Prausnitz and co-workers in the development of UNIQUAC 38 , ln γ2 may be interpreted as a sum of a residual and combinatorial term. The residual term is the result of intermolecular interactions whereas the combinatorial term results from the size dissimilarity of the solute and solvent. If we assume solute-solute interactions are negligible, the residual term (or the intermolecular interactions of the system) are fully captured by the infinite dilution limit. The composition dependent ln γ2∗ can be fully captured using an analytic combinatorial activity coefficient (or excess Gibbs free energy) model (see ref. 44 for three examples). Two well known models are the athermal Flory-Huggins equation (which is adopted by Wilson’s equation) and the extension of Staverman and Guggenheim (which is adopted by UNIQUAC) 38 . Adopting here the athermal Flory-Huggins equation

ln γ2∗ (T, p, x2 ) = − ln [x2 r + (1 − x2 )] + (r − 1) ϕ2

(17)

where r is a measure of the size dissimilarity of the solute and solvent and is computed as the ratio


Page 11 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of the (subcooled) liquid molar volume of the solute relative to the liquid molar volume of the solvent, r = v2 /v1 , and where ϕ2 = x2 r/ [x2 r + (1 − x2 )] is the solute volume fraction. The dilute solution approximation results in a simple, analytic expression (eq. (17)) that only uses parameters v1 and v2 which we are already computing via molecular simulation to calculate ln γ2∞ . As we will show for the case of lidocaine, this simple addition drastically improves our solubility predictions over the infinite dilution approximation. When x2 → 0 or r → 1, γ2∗ → 1 and we recover the infinite dilution approximation.

Computational Details Force Fields All nonbonded intermolecular interactions were modeled using a standard Lennard-Jones (LJ) plus Coulombic potential with fixed point charges 45 [( Unb (rij ) = 4εij

σij rij

)12

( −

σij rij

)6 ] +

1 qi qj 4πε0 rij

(18)

where εij and σij are the LJ energy and size parameters, respectively, qi and qj are the atomic partial charge of atom i and j, respectively, and rij is the separation distance between atoms i and j. Unlike intermolecular interactions were calculated using Lorentz-Berthelot combining rules 45 . Molecular models for acetanilide, acetaminophen, phenacetin, methylparaben and lidocaine (see fig. 1) were taken directly from our previous work 42,46 . The models were constructed with the LJ parameters derived from the closest atoms in the Explicit Hydrogen Transferable Potentials for Phase Equilibria (TraPPE-EH) force field 47–50 and all intramolecular parameters were taken from the general AMBER force field (GAFF) 51 . Methyl (CH3 and CH2 ) groups, if present in the molecular structure, were treated as pseudo-atoms with LJ parameters from the TraPPE United Atom (TraPPE-UA) force field 52,53 . To be consisted with the solute models, CO2 was modeled using the TraPPE model of Potoff and Siepmann 54 developed to predict pure component phase



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

equilibrium properties and its binary mixture with propane. The model has previously been shown to capture transport and thermodynamics properties of sc-CO2 with excellent agreement with experiment 22,55 . All C=O bonds were kept fixed using P-LINCS 56–58 , however, due to algorithmic difficulties in maintaining the linear shape of CO2 molecules 56 , we used a flexible O=C=O angle with bending force constant taken from ref. 59; angle bending has been shown previously to have a negligible effect on predicted phase equilibrium of CO2 54 . All solute bonds containing H atoms were kept fixed using P-LINCS. All of the GROMACS force field files used in the present study are provided in the Supporting Information.

Molecular Dynamics All of the MD simulations were performed using GROMACS 4.6.3 60–62 with standard periodic boundary condition in all directions. Simulations were performed for pure sc-CO2 and for sc-CO2 with a single solvated solute molecule (corresponding to the infinite dilution limit) at a temperature of 308.15 K and pressure of 8, 10, 12, 14, 20, 30, 40, and 50 MPa 63 . Initial configurations were generated using Packmol 64,65 with the number of CO2 molecules chosen to yield an equilibrated cubic box with edge length of approximately 4 nm; an initial estimate of the required number of molecules was obtained from reference sc-CO2 density data 66 . Next, 3000 steps of steepest descent minimization were performed to remove any bad contacts that might have been formed during packing, followed by a short 200 ps N pT simulation with the Berendsen barostat 56,67 . During the dynamics, the equations of motion were integrated with the GROMACS stochastic dynamics integrator, corresponding to stochastic or velocity Langevin dynamics integrated with the leapfrog algorithm 56,68,69 with a time step of 2 fs and friction coefficient of 5 ps−1 . The resulting structure was used to perform an additional production N pT simulation using the Parrinello-Rahman barostat 56,70,71 with a characteristic oscillation time of 5.0 ps and the compressibility was set to 7.4×10−4 bar−1 . The simulations of pure sc-CO2 were 15 ns in length and were used to compute the molar volume; during these runs, the first 2.5 ns was discarded as additional equilibration, and 12 ACS Paragon Plus Environment

Page 12 of 50

Page 13 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the average and uncertainty (or error) was computed using the GROMACS tool “g_energy” by breaking the remaining 12.5 ns into five blocks of 2.5 ns each. The simulations with a single solvated solute molecule were 22.5 ns in length with the first 2.5 ns discarded as equilibration. These trajectories were used for subsequent structural analysis. All of the structural analyses performed in the present study were performed with the TRAVIS package 72,73 . Free energy calculations were subsequently performed for the solute, which will be discussed momentarily. Short-range LJ interactions were truncated at 1.4 nm and a mean-field long-range correction was applied to the energy and pressure to account for the truncation of these interactions 45,56,69 . A Verlet neighbor list was used with a target energy drift of 0.005 kJ/mol/ns per atom 56 . Long-range electrostatic interactions were evaluated with the smooth particle-mesh-Ewald (SPME) method 56,69,74 with tin-foil boundary conditions with real space interactions truncated at 1.4 nm. A SPME Bspline order of 4, a grid spacing of 0.12 nm, and a relative accuracy of 10−4 was used.

Free Energy Calculations The dimensionless residual chemical potential of the solutes infinitely dilute in sc-CO2 , βµres,∞ , 2 was computed using a multi-stage free energy perturbation method 14,75–78 with the multi-state Bennett’s acceptance ratio method (MBAR) 79–82 . A “soft-core” potential was used to decouple the solute-sc-CO2 intermolecular LJ interactions. Stage (m) dependent decoupling parameters, λLJ m and λelec m , controlled the LJ and electrostatic intermolecular interactions, respectively. The decouelec pling parameters varied from 0 to 1. When λLJ m = λm = 1, the solute is fully coupled to the elec system. When λLJ m = λm = 0, the solute is decoupled from the system. The “soft-core” potential

had the form 83–85

{ sc (rij ; m) = 4λLJ ULJ m εij

σij6

σij12

] ] −[ [ LJ ) α σ 6 + r 6 6 6 2 (1 − λ LJ ) α σ + r (1 − λLJ m ij ij LJ m ij ij

} (19)

where αLJ is a constant, which had a value of 1/2. The advantage of using a “soft-core” potential to



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 50

decouple the LJ interactions is that while it yields the correct limiting value of the potential (when λLJ m = 0 and 1), it additionally allows nearly decoupled molecules to overlap with a finite energy (and hence finite probability). The electrostatic term in the intermolecular potential was decoupled linearly as

Uelec (rij ; m) = λelec m

1 qi qj 4πε0 rij

(20)

At each stage m, an independent 6.5 ns MD simulation was performed, where the first 0.5 ns was discarded from analysis as equilibration. The initial structure for these simulations was taken as the final structure from the production N pT simulation of a single solute in solution. The change in the Hamiltonian with the current configuration between stage m and all other stages was computed every 0.20 ps. This information was saved for subsequent post-simulation analysis with MBAR 82 to determine βµres,∞ . This analysis was performed using the Python implementation of 2 MBAR (PyMBAR) and the GROMACS analysis script distributed with it 86–88 . A total of 10 different stages were used for the free energy calculations where m = 0 corresponds to a non-interacting (ideal gas) state and m = 9 is a fully interacting system. The LJ interactions were increased linearly following λLJ m ={0.2, 0.4, 0.6, 0.8, 1.0} from m = 1 to 5. Electrostatic interactions were increased in a square root fashion following λelec ={0.50, 0.71, m 0.87, 1.00} from m = 6 to 9. The simulation parameters for the free energy calculations were the same as the production N pT simulation except the LJ interactions were switched off smoothly between 1.35 and 1.4 nm, and the real space part of the electrostatic interactions were switched off smoothly between 1.38 and 1.4 nm. Long range corrections were applied as described previously. Additionally, the efficient Verlet neighbor list was not available for free energy calculations in GROMACS 4.6.3. The system sizes (N1 ), the calculated pure sc-CO2 molar volumes, and the results of all of the free energy calculations are provided in Tables S1 and S2 of the Supporting Information.


Page 15 of 50

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Results and Discussion Structural Analysis To investigate the solvation mechanism of the studied solutes in sc-CO2 , partial (site-site) radial distribution functions (RDFs) for the C and O atoms of CO2 with the various sites of interest in acetaminophen, lidocaine and methylparaben are shown in the top, middle, and bottom pane of fig. 2, respectively. Acetanilide, acetaminophen and phenacetin share similar functional groups in the N-phenylethanamide scaffold, so for the sake of briefness, only the RDFs for the sites of acetaminophen are presented 89 . All of the studied solutes have a carbonyl group in their chemical structure. As shown by the red curves, for the studied solutes the =C(CO2 )· · · O=(solute carbonyl) RDFs are all very similar and show a strong first peak at 0.3 nm, indicating a highly favorable interaction. This is possibly a Lewis acid-base interactions for which the =C· · · O= separation distance was previously proposed to be around 0.29 nm 90 . For the case of acetaminophen, a second weak peak positioned at 0.6 nm is also visible suggesting that CO2 molecules are more structured and form a second solvation shell around the carbonyl group. To obtain further details of the spatial orientation of CO2 molecules around the carbonyl group, the combined angular-radial distribution function (ARDF) for the acetaminophen-CO2 system is shown in panel a) of fig. 3. Within this figure, the separation distance (r) is between the carbon atom of CO2 (=C) and the carbonyl oxygen of acetaminophen (O=), while the angle is defined by the acetaminophen carbonyl C=O bond vector and the vector connecting the oxygen atom in the carbonyl group of acetaminophen to the carbon atom of CO2 (=C· · · O=). From fig. 3, it is clear that CO2 molecules tend to align perpendicular to the C=O bond vector ( i.e., 135

Microscopic Structure and Solubility Predictions of ... - ACS Publications

Recommend Documents