Developing a Predictive Form of MOSCED for Nonelectrolyte Solids

Article pubs.acs.org/IECR

Developing a Predictive Form of MOSCED for Nonelectrolyte Solids Using Molecular Simulation: Application to Acetanilide, Acetaminophen, and Phenacetin Ryan T. Ley, Georgia B. Fuerst, Bryce N. Redeker, and Andrew S. Paluch* Department of Chemical, Paper and Biomedical Engineering, Miami University, Oxford, Ohio 45056, United States S Supporting Information *

ABSTRACT: The Modified Separation of Cohesive Energy Density Model (MOSCED) is an efficient, analytic method to predict infinite dilution activity coefficients over a range of temperatures. Its predictability makes MOSCED an attractive engineering design tool. However, its use is limited. When trying to model a novel compound, reference data must first be available to regress the necessary MOSCED parameters. Here, we propose the use of molecular simulation to generate the reference dataset. In this fashion, MOSCED can be made a truly predictive engineering design tool. This combines the predictive strength of molecular simulation with the efficiency of MOSCED to create a powerful new tool. Here, we use molecular simulation to generate MOSCED parameters for the nonelectrolyte solid solutes acetanilide, acetaminophen, and phenacetin. Adopting the melting point temperature and enthalpy of fusion of these compounds from available experimental data, we are able to predict equilibrium solubilities. Predictions using the new predictive MOSCED are in good agreement with available experimental solubility data for acetaminophen in nonaqueous solvents.

■

INTRODUCTION Knowledge of the phase behavior of nonelectrolyte solids in solution is critical in the design of separation processes and to model and understand a wide range of physical processes. The equilibrium solubility and relative solubility (or partition coefficient) are important to model both the transport of environmental contaminants and the design of extraction processes for their removal from wastewater streams. Similarly, it is crucial for the design of drug delivery processes and to understand bioavailability. Moreover, solubility of pharmaceuticals in a wide range of solvents is crucial for the design of synthesis, separation, and purification processes to manufacture the drug for general use. Countless other physical and industrial processes exist. For the design and optimization of separation processes, it is important that predictive methods for equilibrium solubility and relative solubility calculations exist. Historically, these methods are often limited to empirical and semiempirical correlations and group contribution methods, with typical focus on aqueous systems.1 UNIFAC has long been a workhorse for chemical engineering design.2−5 However, challenges exist to use UNIFAC to model novel, complex nonelectrolyte solids (such as pharmaceutical compounds). Parameters are missing for many complex solid solutes, and predictions with existing parameters have shown appreciable error.5−7 Efforts are currently being made to fill these voids.5,8 Predictive, group contribution cubic equations of state are also available and are capable of predicting a range of phase equilibriums, including solid−liquid equilibrium.4,9 As compared to UNIFAC (or use of an excess Gibbs free-energy model), an equation of state is capable of modeling the pressure © XXXX American Chemical Society

dependence of the excess Gibbs free energy, and is therefore capable of modeling solubility in supercritical solvents. However, in addition to requiring the availability of a limited number of pure component properties (which include the critical temperature, pressure, and accentric factor), these methods are dependent on a predictive excess Gibbs freeenergy model, such as UNIFAC or one of its modified forms, to employ predictive excess Gibbs free-energy mixing rules. This leads to the same shortcomings as mentioned above for UNIFAC when modeling complex nonelectrolyte solids. Good solubility predictions have been made with the theoretically based NRTL-SAC,10 Modified Separation of Cohesive Energy Density Model (MOSCED),11,12 and PCSAFT equation of state.13−15 Normally, however, solubility data are required for the solute of interest to determine the necessary model parameters prior to making predictions in the solvent of interest. Thus, while the methods are theoretically based and can yield accurate predictions, the methods themselves are not truly predictive. That is, if we seek to screen a massive number of solvents to design a separation process for a novel compound for which experimental data are not available, these methods would fail. Therefore, while the methods are rapid, efficient, can be highly accurate, and are analytic in nature (which can be useful for process design and optimization), they suffer in that they are not truly predictive. Note that, in principal, all of the necessary PC-SAFT Received: December 16, 2015 Revised: March 1, 2016 Accepted: April 12, 2016

A

DOI: 10.1021/acs.iecr.5b04807 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

state are capable of modeling solid−liquid equilibrium, as is done here. In the present study, molecular simulation is used to parametrize MOSCED for acetanilide, acetaminophen, and phenacetin. This is accomplished by computing the solute infinite dilution activity coefficient in 11 or 12 solvents, which may then be used to regress the missing MOSCED parameters. For all three solutes, sigma profiles and FORTRAN programs are available to perform COSMO-SAC calculations for comparison. 27,28,37,38 For the case of acetaminophen, MOSCED parameters additionally exist, which were determined by fitting directly to experimental solubility data.11 This allows us to benchmark our results. We find that the proposed method is able to make reasonable predictions for acetaminophen in nonaqueous solvents. However, we find that the predictions are sensitive to small changes in the solvent force field and the solvent set used to generate the reference data used to parametrize MOSCED. While promising, further studies are required to better understand these challenges.

parameters could be computed using electronic structure methods.13 However, as common with equation of state methods, the necessary mixing rules are difficult to predict a priori.13,16 On the other hand, several recent studies have demonstrated the use of molecular simulation to predict solubility devoid of experimental data.17−23 These methods can predict solubility in the absence of any experimental data as long as the crystal structure is known or may be predicted. Predictions could also be made with use of a limited amount of data of the pure solid solute. Alternatively, relative and excess solubility predictions may be made in the absence of any experimental data and without modeling the solid phase.24−26 Nonetheless, a few hurdles still exist. Primarily, it requires access to highperformance computing resources to perform the calculations, and during solvent screening, when a large number of pure solvents and even larger number of solvent mixtures are potential candidates, the computational expense might be prohibitive. They also lack an analytic form that would be advantageous for process development and optimization. Therefore, we find that, while molecular simulation is attractive, because of its purely predictive nature, for many applications, it may be computationally expensive and impractical. We seek a rapid, efficient, and accurate predictive method to predict the equilibrium solubility and relative solubility of nonelectrolyte solids in a wide range of solvents. In light of the previous discussion, we therefore propose to combine molecular simulation with MOSCED. Molecular simulations may be performed to compute the infinite dilution activity coefficient for a solute of interest in a limited range of reference solvents. This reference data may then be used to directly parametrize MOSCED. By combining these two methods, we are effectively developing a predictive form of MOSCED that is highly useful for design processes. How the methods are combined and the benefit of doing so will be explored further in the section entitled “Methodology”. Note that MOSCED and NRTL-SAC have many similarities, and the suggested approach could just as well be applied to parametrize NRTL-SAC. The proposed method is similar in spirit to COSMO-RS and COSMO-SAC.27−31 In these methods, electronic structure calculations are performed to compute a molecule’s sigma profile. The sigma profiles may then be used within a thermodynamic model to predict composition dependent activity coefficients. Sigma profiles are pure component properties, allowing sigma profile databases to be readily generated. Here, we instead use classical molecular simulation along with MOSCED and an appropriate excess Gibbs freeenergy model. Molecular simulations are used to regress a set of solute MOSCED parameters, analogous to a sigma profile, as in COSMO-RS and COSMO-SAC. Similarly, the proposed method is comparable to the multiscale Gibbs−Helmholtz constrained equation of state approach of Lucia and co-workers.32−36 In both the present study and the work of Lucia and co-workers, molecular simulation is used to make a thermodynamic model predictive. While the present study uses molecular simulation free-energy calculations to parametrize MOSCED (and, ultimately, an excess Gibbs free-energy model), Lucia and co-workers use simulations to compute the internal energy of departure (or equivalently the residual internal energy) which can be used directly to parametrize a cubic equation of state, applicable to both pure component systems and mixtures. Cubic equation of

■

METHODOLOGY The equilibrium solubility of a nonelectrolyte solid solute in pure and mixed solvents may be described by the classical equations of solid−liquid equilibrium. For the case of a solid solute (component 2) in a pure solvent (component 1), we have16 ⎛ f S (T , p ) ⎞ ⎟⎟ ln x 2 = −ln γ2(T , p , x 2) + ln⎜⎜ 2 0 ⎝ f2 (T , p) ⎠

(1)

where x2 and γ2 are the solute mole fraction and the Lewis− Randall normalized activity coefficient at equilibrium, respectively; T is the temperature, p is the pressure, and f 2S and f 20 are the fugacity of the pure solid solute and pure subcooled liquid solute, respectively. Equation 1 is rigorously correct and assumes just that the solid phase in equilibrium with the solution phase is pure. The latter term is a property only of the pure solute, which may readily be estimated using limited properties of the solute at the melting point16,39−43 or may be predicted using molecular simulation as long as the solute solid crystal structure is known.19,44−46 The present study focuses on the calculation of the challenging composition-dependent parameter γ2. We seek an efficient predictive method useful for solvent selection during design processes. This is accomplished here by combining the efficient and analytic MOSCED with the predictive power of molecular simulation. In the present study, the proposed methodology is applied to the nonelectrolyte solids acetanilide, acetaminophen, and phenacetin (see Figure 1). MOSCED. The MOSCED model is an extension of regular solution theory applicable to polar and associating compounds.11,47 Within MOSCED, a compound is modeled using the following six pure component parameters: the liquid molar volume at 20 °C (v), the dispersion parameter (λ), the induction parameter (q), the polar parameter (τ), and the acidity and basicity parameters (α and β, respectively). The liquid molar volume is typically known, and q is commonly fixed; for all of the aromatic centered solutes, we will study q = 0.9.11 The remaining four parameters may be viewed as adjustable and are fit to referenced infinite dilution activity coefficient data. The parameters all have a physical basis; α and B


Article


must be available from which the parameters may be regressed. Early in the design process, these data likely are unavailable. Moreover, it is challenging to obtain the necessary data for a solid solute. The infinite dilution activity may not be directly measured. Rather, one would measure the equilibrium solubility. Subsequently, eq 1 may be applied in reverse with reference data for the pure solid, along with an excess Gibbs free-energy model (to capture the composition dependence of the activity coefficient) to estimate the infinite dilution activity coefficient. From these data, MOSCED may be parametrized. This approach was recently adopted to parametrize MOSCED for a wide range of nonelectrolyte solid solutes.11,48 However, it is limited for design applications, because it is not truly predictive. In this study, we will circumvent this challenge by using molecular simulation to directly compute infinite dilution activity coefficients to generate the reference data necessary to parametrize MOSCED. This offers the advantage of making MOSCED predictive in nature, while minimizing the number of time-intensive molecular simulations that otherwise must be performed. Once parametrized, MOSCED allows for the prediction of infinite dilution activity coefficients. To predict the equilibrium solubility using eq 1, the composition dependence of the activity coefficient of the solute is required. With knowledge of ∞ γ∞ 1 and γ2 , parameters may be obtained for an excess Gibbs free-energy model containing two binary parameters per binary pair, such as Wilson’s equation and the UNIQUAC excess Gibbs free-energy models.16 For Wilson’s equation applied to a binary mixture, we have16,40

Figure 1. Chemical structures of acetanilide (CAS No. 103-84-4), acetaminophen (CAS No. 103-90-2), and phenacetin (CAS No. 6244-2).

β are only assigned nonzero values if deemed physically reasonable.11 The infinite dilution activity coefficient of component 2 in component 1 (γ∞ 2 , in a binary mixture) may be computed using the following series of equations: ln γ2∞ =

⎡ q 2q 2(τ1(T ) − τ2(T ))2 v2 ⎢ (λ1 − λ 2)2 + 1 2 RT ⎢⎣ ψ1 +

(α1(T ) − α2(T ))(β1(T ) − β2(T )) ⎤ ⎥ + d12 ⎥⎦ ξ1

⎛ v ⎞aa ⎛ v ⎞aa d12 = ln⎜ 2 ⎟ + 1 − ⎜ 2 ⎟ ⎝ v1 ⎠ ⎝ v1 ⎠

⎡ ⎤ Λ12 Λ 21 ln γ1 = − ln(x1 + Λ12x 2) + x 2⎢ − ⎥ Λ 21x1 + x 2 ⎦ ⎣ x1 + Λ12x 2

aa = 0.953 − 0.002314[(τ2(T ))2 + α2(T )β2(T )] αi(T )

⎡ ⎤ Λ12 Λ 21 ln γ2 = − ln(x 2 + Λ 21x1) − x1⎢ − ⎥ Λ 21x1 + x 2 ⎦ ⎣ x1 + Λ12x 2

⎛ 293 K ⎞0.8 (T ) ⎛ 293 K ⎞0.8 (T ) ⎛ 293 K ⎞0.4 ⎟ ⎟ ⎟ , βi = βi ⎜ , τi = τi⎜ = αi⎜ ⎝ T ⎠ ⎝ T ⎠ ⎝ T ⎠

(3)

(where i = 1 or 2)

where Λ12 and Λ21 are adjustable parameters. References 11 and 48 compare the use of Wilson’s equation and UNIQUAC for modeling solid−liquid equilibrium with MOSCED and suggest the use of Wilson’s equation, motivating our decision here to use only Wilson’s equation. Wilson’s equation and UNIQUAC may both be generalized to multicomponent mixtures, allowing for solubility predictions to additionally be made in multicomponent solvents. The parameters necessary to model multicomponent mixtures may be obtained for each binary pair that constitutes the mixture.16,40 The ability to predict MOSCED parameters in the absence of experimental data would be of great value. Next, we show how molecular simulation may be used to obtain the necessary data, allowing one to take advantage of the efficiency of MOSCED while using the predictive nature of molecular simulation. Molecular Simulation. Compared to MOSCED, molecular simulation is a predictive tool. As shown previously, the infinite dilution activity coefficient of the solute (component 2) may be computed using conventional molecular simulation free-energy calculations as49,50

ψ1 = POL + 0.002629α1(T )β1(T ) 2

ξ1 = 0.68(POL − 1) + [3.4 − 2.4 exp(−0.002687(α1β1)1.5 )](293 K/ T ) POL = q14[1.15 − 1.15 exp(−0.002337(τ1(T ))3 )] + 1

(2)

where R is the molar gas constant and T is the absolute temperature in Kelvin, and where an equivalent expression may be written for component 1 by switching the subscript indices. The central feature of these equations is that six pure component parameters are needed for each component in the binary mixture. Of these six parameters, four may be viewed as adjustable (λ, τ, α, and β). In addition, while v, q, and λ are assumed to be independent of temperature, MOSCED accounts for the temperature dependence of α, β, and τ, allowing γ1∞ and γ2∞ to be estimated as a function of temperature. However, before this may be done, parameters are needed for the compounds of interest. This requires having reference data to regress the necessary parameters. In the recent 2005 MOSCED revision, parameters were obtained for 133 solvents (including water) using an extensive database of experimental infinite dilution activity coefficients.11 While MOSCED is parametrized for a wide range of 133 unique solvents, the difficulty of applying MOSCED for solvent selection during design processes results from the absence of parameters for novel solid solutes of interest. Reference data

ln γ2∞(T , p , 0) = βμ2res, ∞(T , p , N1 , N2 = 1) ⎛ RT ⎞ + ln⎜ ⎟ − ln f2 0 (T , p) ⎝ v1(T , p) ⎠ C

(4)


Article

Industrial & Engineering Chemistry Research where βμres,∞ is the dimensionless residual chemical potential of 2 the solute in solution at infinite dilution (β−1 = kBT, where kB is the Boltzmann constant); v1 is the molar volume of the pure solvent; and N1 and N2 are the number of solvent and solute molecules in the molecular simulation, respectively. The parameter f 20 is the pure liquid fugacity of the solute, which may be computed as51 ⎛ RT ⎞ ln f2 0 (T , p) = βμ2res (T , p) + ln⎜ ⎟ ⎝ v2(T , p) ⎠

TraPPE-UA.58−61 The use of UA models is highly advantageous, because it results in a significant reduction of the CPU requirement. Benzene and phenol were modeled using TraPPE-EH,54−56 with intramolecular parameters adopted from the General AMBER Force Field (GAFF),62 which is consistent with the model for the studied solutes. All of these TraPPE models were developed with the goal of reproducing reference pure component vapor−liquid coexistence data. While TraPPE does not contain a water model, we adopted the TIP4P model,63 which has been shown previously to work well with TraPPE.64−66 Chloroform was modeled using the TraPPE-related force field of Potoff and co-workers,67,68 which was optimized to reproduce pure and binary vapor−liquid equilibrium data with a TraPPE-UA-modified acetone model. The same work demonstrated the model worked well with the TraPPE-UA methanol model to predict binary vapor−liquid equilibrium data. We additionally considered the case of using the acetone model of Potoff and co-workers67,68 and its extension to butanone, in place of TraPPE-UA.68 The second set of solvent parameters were adopted based on the Kirkwood−Buff-derived force field (KBFF) of Smith and co-workers.69−72 While KBFF does not explicitly contain parameters for linear and branched alkanes, consistent with KBFF nonbonded parameters were taken from the (improved) united atom GROMOS96 force field.73 Methanol and acetone were united atom models taken directly from refs 69 and 70, respectively. These models were extended to ethanol, 1propanol, 2-propanol, and butanone by adopting missing nonbonded parameters from the (improved) united atom GROMOS96 force field. Partial charges were the same as methanol for ethanol and 1-propanol, and the same as acetone for butanone. For 2-propanol, partial charges were taken from ref 74. The molecular geometry, angle bending force constants, and the dihedral potentials for the studied alcohols were adopted from TraPPE-UA.58,60,61,75 For acetone and butanone, the geometry and intramolecular parameters for the carbonyl group were taken from ref 70, and missing intramolecular parameters for butanone were adopted from TraPPE-UA.59 Benzene and phenol were explicit hydrogen models taken from ref 72 with intramolecular parameters adopted from GAFF.62 KBFF methanol and acetone were parametrized to reproduce the composition-dependent Kirkwood−Buff integrals (and, hence, excess Gibbs free energy) for aqueous methanol and acetone solutions under ambient conditions with the SPC/E water model.69,70,76 KBFF benzene and phenol were parametrized to reproduce the composition-dependent Kirkwood− Buff integrals for mixtures with methanol.72 The use of KBFF to model the solvents in the present study was motivated by the common use of solvent mixtures to solvate nonelectrolyte solids and the importance of modeling solvent mixtures accurately.25 Consistent with KBFF, water was modeled with the SPC/E water model.76 Solvent Set Summary. Before moving on, it is important to clarify the solvent sets that will be considered when generating the reference data necessary to parametrize MOSCED. (1) TraPPE I consists of the following 11 solvents: n-hexane, 2,5-dimethylhexane, methanol, ethanol, 1-propanol, 2-propanol, acetone, butanone, benzene, phenol, and water. Water is modeled with TIP4P, benzene and phenol are modeled with TraPPE-EH, and all other solvents are modeled with TraPPEUA. (2) TraPPE II consists of the 11 solvents included in TraPPE I plus chloroform. Water is modeled with TIP4P, benzene and

(5)

where βμres 2 and v2 are the dimensionless residual chemical potential and molar volume of the pure liquid solute. For a solute that is liquid under the conditions of interest, these quantities may be computed using a single molecular simulation. However, challenges exist for solutes that are solid under the conditions of interest for which these correspond to properties of a pure subcooled liquid. Our recent work has demonstrated that, for these cases, direct calculations of the subcooled liquid should be avoided; instead, a series of simulations should be conducted at elevated temperatures and extrapolated to the conditions of interest.51 The three solutes considered here (acetanilide, acetaminophen, and phenacetin) are all solid under ambient conditions and were studied extensively in ref 51, wherein f 20 was computed for each solute. These calculations were not repeated in the present study; instead, values of f 20 were taken directly from ref 51.

■

COMPUTATIONAL DETAILS Force Fields. All of the interactions for the molecular simulations were modeled using a conventional “Class I” potential energy function with all the nonbonded intermolecular interactions (Unb) treated using a combined LennardJones (LJ) plus fixed-point-charge model of the form52,53 ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σij σij 1 ⎛⎜ qiqj ⎞⎟ Unb(rij) = 4εij⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥ + r 4πε0 ⎜⎝ rij ⎟⎠ ⎝ rij ⎠ ⎦ ⎣⎝ ij ⎠

(6)

where rij, εij, and σij are the site separation distance between atoms i and j, well depth of the LJ interaction, and distance at which the LJ interaction is zero, respectively. The parameters qi and qj are the partial charge values of atoms i and j, respectively. Molecular models for acetanilide, acetaminophen, and phenacetin (see Figure 1) were taken directly from our previous work.51 The models are based on the Explicit Hydrogen Transferable Potentials for Phase Equilibria (TraPPE-EH) force field54−57 with additional parameters from the TraPPE United Atom (TraPPE-UA) force field.58,59 GROMACS force field files (type itp and xyz) for the studied solutes are provided in the Supporting Information of this manuscript. The solvents studied were n-hexane, 2,5-dimethylhexane, methanol, ethanol, 1-propanol, 2-propanol, acetone, butanone, benzene, phenol, chloroform, and water. The solvents were selected based on the availability of force field parameters and the desire to use solvents containing a range of chemical functionality. The solvents were all modeled with two different force field families. First, parameters were adopted from the TraPPE force field. N-Hexane, 2,5-dimethylhexane, methanol, ethanol, 1-propanol, 2-propanol, acetone and butanone were modeled using D


Article


error) was computed using the GROMACS tool “g_energy”82 by breaking the remaining 12.5 ns into five blocks of 2.5 ns each. Solute in Solution. The MD simulations for each solute in each solvent involved five steps: initial structure generation, minimization, NpT equilibration and production, and freeenergy calculations. First, initial structures were generated. The (final) equilibrated solvent boxes described above were used to solvate a single solute molecule in a cubic box with edge length of 45 Å using the GROMACS tool “genbox”.82 Steps 2−4 were exactly as described above for the pure solvents. In step 5, the infinite dilution residual chemical potential (μres,∞ ) of the solute in solution was computed using a 2 multistage free-energy perturbation approach92−96 with the multistate Bennett’s acceptance ratio method (MBAR)97−100 in an almost-identical manner, as in our recent publication for which this work builds upon.51 Intermolecular LJ and electrostatic interactions are regulated by the stage (m)dependent coupling parameter λmLJ and λmelec, respectively, elec LJ which vary from 0 ≤ λLJ m ≤ 1 and 0 ≤ λm ≤ 1. When λm = LJ λelec = 1, the solute is fully coupled to the system. When λ m m = elec λm = 0, the solute is decoupled from the system. Solute− solvent intermolecular LJ interactions were decoupled using a “soft-core” potential of the form101−103

phenol are modeled with TraPPE-EH, n-hexane, 2,5dimethylhexane, methanol, ethanol, 1-propanol, and 2-propanol are modeled with TraPPE-UA, and chloroform, acetone, and butanone are modeled using the TraPPE-related models of Potoff and co-workers.67,68 (3) TraPPE III consists of the 11 solvents included in TraPPE I, modeled exactly the same as in TraPPE I, plus chloroform (i.e., TraPPE III is TraPPE I plus chloroform). Chloroform is modeled using the TraPPE-related model of Potoff and co-workers.67,68 (4) KBFF consists of the 11 solvents included in TraPPE I. Here, however, water is modeled with SPC/E and all of the other solvents are modeled using the KBFF force field. Molecular Simulations. The present study requires the calculation of configurational properties, allowing one to use either Monte Carlo or molecular dynamics (MD) simulations. Here, we chose to use MD to sample the configurational phase space. All of the MD simulations were performed with GROMACS 4.6.377−79 at 298.15 K and 1 bar. The MD simulations may be split into two categories: simulations involving pure solvent and simulations for the solute in solution. Pure Solvents. The MD simulations for each pure solvent involved four steps: initial structure generation, minimization, and NpT equilibration and production. First, initial structures were generated using Packmol.80,81 The number of molecules was chosen to obtain an equilibrated cubic box with an edge length of ∼45 Å. Second, up to 3000 steps of steepest descent minimization were performed on the initial structure to remove any bad contacts that may have been created during packing. For steps 3 and 4 for 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.82−84 The third step consisted of 100 ps of NpT equilibration with the Berendsen barostat,82,85 followed in step 4 with a 15 ns NpT production simulation with the Parrinello−Rahman barostat.82,86,87 The time constant for the stochastic (or Langevin) thermostat was 1.0 and 2.0 ps in steps 3 and 4, respectively, the time constant for the barostat was 5.0 ps, and all bond lengths were constrained using P-LINCS.82,88,89 Note that, for the case of phenol modeled with TraPPE-EH, the time constant for the thermostat was increased to 8.0 ps, because of the long correlation times encountered when performing free-energy calculations. The Verlet cutoff scheme in GROMACS was employed with the LJ interactions truncated at 14.0 Å. Longrange analytic dispersion corrections were applied to the energy and pressure to account for the truncation of these interactions.52,53,82,84 For interactions between unlike LJ sites, Lorentz−Berthelot combining rules were employed for simulations with TraPPE solvents and geometric combining rules for simulations with KBFF solvents.52 Electrostatic interactions were evaluated using the smooth particle-meshEwald (SPME) method with tin-foil boundary conditions82,84,90,91 with real space interactions truncated at 14.0 Å. A SPME B-spline order of 6, a Fourier spacing of 1 Å, and a relative tolerance between long- and short-range energies of 10−6 were used. The molar volume of the solvent necessary for eq 4 was computed during the 15 ns NpT production simulation (step 4). During these runs, the first 2.5 ns was discarded as additional equilibration, and the average and uncertainty (or

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

⎧ ⎫ σij12 σij 6 ⎪ ⎪ ⎬ ×⎨ − LJ 6 6 2 LJ 6 6 ⎪ ⎪ [(1 − λm )αLJσij + rij ] ⎭ ⎩ [(1 − λm )αLJσij + rij ]

(7)

where αLJ is a constant, taken in this study to be 1/2. The intermolecular electrostatic interactions are decoupled linearly as ⎛ 1 ⎞⎛ qiqj ⎞ ⎟⎟ Uelec(rij ; m) = λmelec⎜ ⎟⎜⎜ ⎝ 4πε0 ⎠⎝ rij ⎠

(8)

An independent MD simulation is performed in each stage m. Periodically throughout the simulation, the change in the Hamiltonian with the current configuration between stage m and every other stage is computed. This information is saved to disk for post-simulation analysis with MBAR100,104−106 to compute μres,∞ . 2 In the present study, the solute was taken from noninteracting (m = 0) to fully interacting (m = 9) by first bringing the intermolecular LJ interaction to full strength over 5 stages (1 ≤ m ≤ 5), and then adding in intermolecular electrostatic interactions in the final 4 stages (6 ≤ m ≤ 9) for a total of 10 stages (including the ideal gas reference state m = 0). In stage m = 0, the intermolecular LJ and electrostatic interactions elec between the solute and the system are turned off (λLJ 0 = λ0 = 0). Over the next 5 stages, the intermolecular electrostatic interactions remain off and the intermolecular LJ interactions are strengthened linearly as λLJ m = {0.2, 0.4, 0.6, 0.8, 1.0}. Next, with the LJ interactions fully restored (λmLJ = 1), the intermolecular electrostatic interactions were strengthened following a square root relation as λelec m = {0.50, 0.71, 0.87, 1.00}.91 During this process, all bonded and intramolecular LJ and electrostatic interactions are unchanged. The final configuration from the corresponding 15 ns NpT production simulation was taken as the starting configuration in each stage. The NpT ensemble MD simulation in each stage was 9.0 ns in length, with the first 0.5 ns being discarded from E


Article

Industrial & Engineering Chemistry Research analysis as equilibration. During these simulations, the differences in the Hamiltonians was computed every 0.20 ps. The parameter μres,∞ and the associated error was computed 2 using the GROMACS analysis script distributed with the Python implementation of MBAR (PyMBAR).100,104−106 Simulation parameters for this step were the same as step 4, except the efficient Verlet cutoff scheme is not available for use with the free-energy calculations in GROMACS 4.6.3. Instead, the LJ interactions were switched off smoothly between 13.5 Å and 14.0 Å, and the real space part of the electrostatic interactions were switched off smoothly between 13.8 Å and 14.0 Å. As already mentioned, when performing free-energy calculations with TraPPE-EH phenol, we encountered long correlation times. To overcome this problem, we made three changes to these simulations. First, as in step 4, the time constant for the thermostat was increased from 2.0 ps to 8.0 ps. Second, the simulation length was increased from 9.0 ns to 17.5 ns, where the first 0.5 ns was discarded from analysis as equilibration. And third, the number of LJ stages were doubled, so that the LJ interactions were strengthened as λLJ m = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0}. Lastly, computing the solute infinite dilution activity coefficient (normalized with respect to a Lewis−Randall standard state, eq 4) requires knowledge of the pure solute subcooled liquid fugacity. This is accomplished by performing free-energy calculations for the solute in itself at temperatures above the normal melting point. These results may then be extrapolated to the temperature of interest. The procedure is discussed in detail in our recent work from which the pure component fugacity at 298.15 K and (subcooled liquid) molar volume at 293.15 K (20 °C) for acetanilide, acetaminophen, and phenacetin is adopted.51

Table 1. Tabulated Values of the Solvent Molar Volume (v1) and the Dimensionless Infinite Dilution Residual Chemical Potential of Acetaminophen in Each Solvent (βμres,∞ ) 2 Computed Using Molecular Simulation with Four Different Solvent Sets

■

RESULTS AND DISCUSSION The overarching goal of the present study is to investigate the use of molecular simulation to develop a predictive form of MOSCED useful for nonelectrolyte solids in solution. However, before it may be applied to a particular solute, reference data must be available to regress the necessary MOSCED parameters. Here, the use of molecular simulation to generate the necessary reference data is proposed. We consider four sets of reference data to regress the four necessary MOSCED parameters (λ, τ, α, and β), using the solvent sets summarized in the section entitled “Solvent Set Summary”. Parameters are regressed for a solute following a four-step procedure: (1) In each solvent, a free-energy calculation is performed to compute βμ2res,∞. The results of all of the free-energy calculations may be found in Table 1 for acetaminophen and in Table 2 for acetanilide and phenacetin. (2) For each solute, f 20 must be computed by performing free-energy calculations for the pure liquid solute at a series of temperatures above the normal melting point (see eq 5), and then extrapolated to the temperature of interest as described in ref 51. The solute force field in steps 1 and 2 should be consistent (the same). In this step, the liquid molar volume of the solute at 20 °C (MOSCED parameter v) may also be computed. Here, f 20 and v are taken directly from ref 51. (3) In each solvent, γ∞ 2 may be computed using eq 4. (4) With this set of reference data, MOSCED parameters can be regressed. We obtain v in step 2, and for the three aromatic centered solutes studied in this work, the solute induction

solvent

solvent set

v1 (cm3/mol)a

a −βμres,∞ 2

2,5-dimethylhexane

TraPPE I, II, III KBFF

130.81913 136.87318

15.5808 15.3508

n-hexane


130.29612 129.50814

13.9405 14.0305

water


18.13601 18.04101

27.9307 27.5908

methanol


40.91403 41.37904

33.0906 33.6606

ethanol


58.91406 60.47407

32.6108 33.0408

2-propanol


76.89716 80.48318

31.4010 32.6424

1-propanol


75.67910 77.73012

31.9710 32.7112

acetone

TraPPE I, III TraPPE II KBFF

74.71415 74.08312 71.31707

27.5105 31.2505 32.2806

butanone

TraPPE I, III TraPPE II KBFF

91.43006 90.93810 89.10506

29.2706 33.7307 31.7406

phenol


86.22724 88.43207

31.9019 32.9715

benzene


89.38114 88.78212

20.1507 23.3506

chloroform

TraPPE II, III

80.54414

25.9306

a

The subscripted values correspond to the uncertainty in the last two decimal places.

parameter (q) was set to 0.9, as suggested in ref 11. This left the dispersion parameter (λ), the polar parameter (τ), and the acidity and basicity parameters (α and β) unknown. These four parameters were found by minimizing the objective function (OBJ), which is defined as the squared difference between the natural logarithm of the infinite dilution activity coefficient of the solute computed using molecular simulation (“sim”) and the natural logarithm of the infinite dilution activity coefficient of the solute computed using MOSCED, OBJ = (ln γ2∞ ,sim − ln γ2∞ ,MOSCED)2

(9)

Logarithmic terms are used because it is the logarithm of the infinite dilution activity coefficient that is directly related to βμres,∞ . The optimization was performed using the differential 2 evolution method,107 as implemented in GNU Octave.108 All of the MOSCED parameters regressed and used in the present study may be found in Table 3. F


Article


is not the focus of the present work; further details may be found in the Supporting Information. Acetaminophen. We begin with acetaminophen, which has the most experimental reference data for comparison. Figure 2

Table 2. Tabulated Values of the Dimensionless Infinite Dilution Residual Chemical Potential of Acetanilide and Phenacetin in Each Solvent (βμres,∞ ) Computed Using 2 Molecular Simulation with Two Different Solvent Setsa b −βμres,∞ 2

Acetanilide

Phenacetin

solvent

TraPPE III

KBFF

TraPPE III

KBFF

2,5-dimethylhexane n-hexane water methanol ethanol 2-propanol 1-propanol acetone butanone phenol benzene chloroform

13.5307 12.2405 19.4107 23.7705 23.5107 22.7109 22.8308 21.1305 21.9106 24.3627 16.6507 22.1606

13.7308 12.4605 19.0708 24.0406 23.4807 23.5214 23.2309 23.4005 23.0105 24.3713 19.0506

15.9410 14.7306 20.2610 26.1106 25.9508 24.8610 25.3209 24.0706 24.7707 27.6936 19.6109 25.6307

16.0110 15.0206 19.8311 26.1607 25.5708 24.9514 25.7410 26.0906 25.5307 27.1617 21.8608

a

Note that the solvent molar volumes are the same as those given in Table 1. bThe subscripted values correspond to the uncertainty in the last two decimal places.

The focus of this study is on calculating the difficult, composition-dependent activity coefficient of a solute in solution. With regressed MOSCED parameters, infinite dilution activities of the solute in any of the 133 parametrized solvents at any temperature (for which the solvent is liquid) may be predicted. This alone would be sufficient to predict relative solubilities or partition coefficients.16 However, we may additionally use MOSCED to predict the (hypothetical) infinite dilution activity coefficient of the solvent in the (subcooled) liquid solute. We can then parametrize an excess Gibbs freeenergy model (here, Wilson’s equation (eq 3)) to predict composition-dependent activity coefficients. With knowledge of ln( f1S/f10), we can predict the equilibrium solubility using eq 1. Computing the equilibrium solubility is beneficial to the present study, because it allows us to quantitatively compare the results of the present study to available experimental data. To accomplish this, ln( f1S/f10) was calculated using the expression proposed by Nordström and Rasmuson,42 with the solute enthalpy of fusion and the melting point temperature taken from available experimental data. The expression of Nordström and Rasmuson42 permits the calculation of ln( f1S/ f10) at 10, 15, 20, 25, and 30 °C, thereby limiting solubility predictions to these temperatures. Although this is important, it

Figure 2. Parity plot of the infinite dilution activity coefficient predicted using MOSCED with our simulation regressed parameters versus the simulation computed infinite dilution activity coefficients used to regress the MOSCED parameters for acetaminophen. Errors in the simulation results are roughly the size of the symbols or smaller. The subscripts correspond to the uncertainty in the last digit for AAD and the uncertainty in AAD%, where the uncertainty is computed as one standard deviation about the mean.

is a parity plot of ln γ∞,MOSCED versus ln γ∞,sim for the four sets 2 2 of reference data. The deviation between ln γ∞,MOSCED and ln 2 γ∞,sim is quantified using the percent average absolute deviation 2 (AAD%), AAD% =

1 N

N

∑

ln γ2∞ ,MOSCED − ln γ2∞ ,sim ln γ2∞ ,sim

i=1

× 100 (10)

and the average absolute deviation (AAD),

Table 3. MOSCED Parameters of the Present Study Regressed Using Molecular Simulation Generated Reference Data solute

solvent set

v (cm3/mol)

λ ((J/cm3)1/2)

τ ((J/cm3)1/2)

q

α ((J/cm3)1/2)

β ((J/cm3)1/2)

acetaminophen

TraPPE I TraPPE II TraPPE III KBFF

121.1 121.1 121.1 121.1

15.646 15.742 15.368 16.706

1.898 5.536 4.403 2.762

0.9 0.9 0.9 0.9

30.158 33.625 21.222 39.988

9.888 7.501 12.028 7.065

acetanilide

TraPPE III KBFF

121.6 121.6

18.053 17.610

5.786 4.495

0.9 0.9

11.675 19.138

13.073 8.743

phenacetin

TraPPE III KBFF

161.8 161.8

17.141 16.682

5.627 4.823

0.9 0.9

8.671 13.625

12.057 8.477

G


Article

Industrial & Engineering Chemistry Research AAD =

1 N

N

∑ |ln γ2∞ ,MOSCED − ln γ2∞ ,sim| i=1

(11)

where the summation is over all N solvents used to regress the MOSCED parameters. In all cases, the simulation results appear to be reasonably correlated with MOSCED with R2 values of ≥0.89. The AAD and AAD% values are noticeably larger. The AAD is found to be between ∼1 and ∼2 logarithmic units. Comparing relative values of the dimensionless residual chemical potential at infinite dilution of solutes in various nonaqueous solvents (which is related to the relative solubility in molar concentrations) using conventional force fields, we obtained average errors of ∼1 logarithmic unit in a recent study.26 Here, we expect errors both due to the force fields (as in the previous study), in addition to errors in the ability of MOSCED to reproduce experimental data, which is reflected in the larger AAD value. We also note that the regressed MOSCED parameters are sensitive to the training set, which will be further discussed momentarily. Also, it is important to notice that the reported AAD% values are rather large. This is due to the nature of the calculations, which causes the representation of the data by AAD and AAD% to be quite different. For the case of TraPPE III, for = acetaminophen in n-hexane, the simulations give ln γ∞,sim 2 15.4 ± 0.2, while ln γ∞,MOSCED = 12.89. This results in AAD and 2 AAD% values of 2.51 and 16.3, respectively (for this single point). On the other hand, for acetaminophen in 2-butanone, = 0.4 ± 0.2, while ln γ∞,MOSCED = the simulations give ln γ∞,sim 2 2 1.84. This results in an AAD of 1.44, which is ∼1 logarithmic unit smaller than with n-hexane. However, this results in an extremely large AAD% value of 360, which is 22 times larger than that observed with n-hexane. For this reason, we define our objective function as described in eq 9, to minimize AAD, which is consistent with the recent MOSCED revision.11 Otherwise, ln γ∞,sim values that are close to zero would be given 2 greater weight. Note that the plot corresponds to logarithmic values and spans approximately 16 log units. We are able to correlate the results over a wide range of solvent chemical functionality using just four adjustable parameters. In Figure 3, we compare the solubility of acetaminophen predicted using MOSCED with our simulation predicted parameters using the TraPPE I, II, and III sets of reference data against the experimental solubility. Comparing TraPPE I, II, and III will allow us to understand the sensitivity of the predictions on the solvent force fields used to generate the reference data. The experimental data presented in Figure 3, as well as Figures 6, 7, 10, and 11, are all of the nonaqueous, single-component solubility data at 10, 15, 20, 25, and 30 °C available in DECHEMA’s “Solubility and Related Properties of Large Complex Chemicals” series (Part 1109 and Part 2110) for which MOSCED solvent parameters exist. For acetaminophen, this corresponds to 77 experimental solubilities over all 5 temperatures. The deviation between the MOSCED predictions and experiment are quantified using the percent average absolute deviation, AAD% =

1 N

N

∑ i=1

exp |x 2,pred i − x 2, i |

x 2exp

Figure 3. Parity plot of the equilibrium solubility (nonaqueous) of acetaminophen predicted versus experimental data. Predictions are made using MOSCED with simulation regressed parameters using the three TraPPE (I, II, and III)-based reference sets. The symbols are colored orange, blue, green, red, and black, corresponding to data obtained at temperatures of 30, 25, 20, 15, and 10 °C, respectively.

AAD =

1 N

N exp ∑ |x2,pred i − x 2, i | i=1

(13)

where the summation is over all N available data, xpred 2,i is the predicted solubility, and xexp 2,i is the experimental solubility. First, consider the low solubility cases of xpred < 1 × 10−3. 2 This corresponds to the solvents toluene (10, 15, 20, 25, and 30 °C), benzene (25 °C), cyclohexane (25 °C), dichloromethane (30 °C), chloroform (or trichloromethane; 25 and 30 °C), and tetrachloromethane (30 °C). In Figure 4, we compare the AAD value computed using TraPPE III versus TraPPE I and TraPPE II over this composition range. First, consider the case of TraPPE III versus TraPPE I. Recall, TraPPE I and TraPPE III are identical, except the TraPPE III solvent set additionally includes chloroform modeled with the force field of Kamath et al.67 Examining Figure 4, points below and above the y = x line indicate that the corresponding AAD value with TraPPE III is less than and greater than that obtained with TraPPE I, respectively, and points on the y = x line indicate that both perform equally well. Generally, the two perform equally well, except for three points, which are noticeably below the y = x line. These correspond to chloroform at 25 and 30 °C, and dichloromethane at 30 °C. Therefore, this finding suggests that, by including chloroform in the reference data set used to regress MOSCED parameters, resulting predictions in similar solvents is improved. Both TraPPE III and TraPPE I use the same force fields to model the alkane and aromatic solvents, and both appear to be in close agreement for these species. Figure 4 also compares TraPPE III versus TraPPE II. TraPPE III and TraPPE II are identical, except that TraPPE II models acetone and butanone using the modified TraPPE-UA parameters of Kamath et al.,67 which results in a more polar

× 100 (12)

and the average absolute deviation H


Article


Figure 4. Parity plot of the AAD computed using the TraPPE III solvent set versus the TraPPE I or TraPPE II solvent sets for pred 0.92. Also, AAD was found to be less than 1 logarithmic unit in all cases. On the other hand, AAD% is large with a large uncertainty, just as we had found with acetaminophen. A major contribution to the large AAD% value is from ln γ∞,sim values close to zero. For example, the 2 largest contribution to AAD% is for phenacetin in ethanol, using KBFF. For this case, the simulations give ln γ∞,sim = 2 −0.002 ± 0.914, while ln γ∞,MOSCED = 0.173. This results in 2 AAD and AAD% values of 0.175 and 8750, respectively (for this single point), reaffirming our choice of objective function (eq 9). With our simulation parametrized MOSCED models, we proceed to make solubility predictions. In Figures 10 and 11, we compare the solubility of acetanilide and phenacetin, respectively, predicted using our simulation-trained MOSCED model and COSMO-SAC versus experimental solubility data. For both acetanilide and phenacetin, the experimental systems are the same, and we obtain similar results. In both cases, the largest two experimental solubilities are in 1-octanol at 30 °C and 1,4-dioxane at 25 °C. The other six data points correspond to alkanes at 25 °C. For acetanilide, quantitatively for TraPPE III, KBFF, and COSMO-SAC, respectively, we find that AAD% = {61.5, 79.7, and 2415.9} and AAD = {0.006, 0.027, and 0.039}. Quantitatively for phenacetin, for TraPPE III, KBFF, J


Article


Figure 9. Parity plot of the infinite dilution activity coefficient predicted using MOSCED with our simulation regressed parameters versus the simulation computed infinite dilution activity coefficients used to regress the MOSCED parameters for phenacetin. Errors in the simulation results are roughly the size of the symbols or smaller. The subscripts correspond to the uncertainty in the last digit for AAD and the uncertainty in AAD%, where the uncertainty is computed as one standard deviation about the mean.

Figure 8. Parity plot of the infinite dilution activity coefficient predicted using MOSCED with our simulation regressed parameters versus the simulation computed infinite dilution activity coefficients used to regress the MOSCED parameters for acetanilide. Errors in the simulation results are roughly the size of the symbols or smaller. The subscripts correspond to the uncertainty in the last digit for AAD and the uncertainty in AAD%, where the uncertainty is computed as one standard deviation about the mean.

and COSMO-SAC, respectively, we find that AAD% = {78.1, 71.9, and 654.8} and AAD = {0.017, 0.008, and 0.007}. Overall, we find that our simulation-parametrized MOSCED models perform well, with top performance again obtained using MOSCED with the TraPPE III reference data. Water. Until this point in the study, aqueous solubility predictions have been excluded from our comparison, because we have found that the procedure to predict aqueous solubilities must be modified slightly. Figure 12 shows the temperature-dependent solubility of acetaminophen in water predicted using MOSCED 2005 (MOSCED-trained using experimental data). First, consider solubility predictions as performed in the nonaqueous solvents. MOSCED is used to predict infinite dilution activity coefficients of the solute in the solvent and the solvent in the (subcooled) liquid solute, which may then be used to parametrize Wilson’s equation (eq 3). The equilibrium solubility may then be predicted using eq 1. However, as shown in Figure 12, by adopting this approach, we predict that the solubility decreases as the temperature increases. This trend is contrary to that which is expected. This is unacceptable; if the model is not in excellent quantitative agreement with experiment, we expect that it at least be able to predict correct physical trends. The aqueous solubility is typically small. Therefore, we assume that the concentration is sufficiently small to be considered infinitely dilute, and instead use the infinite dilution activity coefficient (which is composition-independent) from Wilson’s equation to predict the solubility via eq 1. As shown in Figure 12, in doing this, we are able to capture the correct physical trend that the solubility increases as the temperature increases. We therefore conclude that, when making predictions

Figure 10. Parity plot of the equilibrium solubility (nonaqueous) of acetanilide predicted versus experimental data. Predictions are made using MOSCED with simulation regressed parameters using the TraPPE III (△) and KBFF (▽) reference sets, and using COSMOSAC (◇). The symbols are colored orange and blue, corresponding to data at 30 and 25 °C, respectively.

K


Article


using MOSCED with Wilson’s equation in water, infinite dilution activity coefficients should be used. By using infinite dilution activity coefficients from Wilson’s equation, we still maintain the ability to model multicomponent solvents while improving the accuracy of the predictions. MOSCED 2005 results were shown, as we expected, to have the greatest agreement with experiment. Analogous findings were obtained with MOSCED using parameters from the TraPPE III and KBFF reference datasets. Having found that infinite dilution activity coefficients should be used to make solubility predictions in water with MOSCED, we next compare predicted solubilities of acetaminophen versus experiment. Solubilities are predicted using MOSCED 2005, MOSCED parametrized using the TraPPE III and KBFF reference datasets, and COSMO-SAC. The solubility data, again, are obtained at temperatures of 10, 15, 20, 25, and 30 °C from Parts 1 and 2 of DECHEMA’s “Solubility and Related Properties of Large Complex Chemicals” work.109,110,113 Figure 13 is a parity plot of the predicted solubility versus available experimental data. We find that all four methods predict the same trend, with the predictions all offset from each other. The accuracy of agreement, from most to least, is Figure 11. Parity plot of the equilibrium solubility (nonaqueous) of phenacetin predicted versus experimental data. Predictions are made using MOSCED with simulation regressed parameters using the TraPPE III (△) and KBFF (▽) reference sets, and using COSMOSAC (◇). The symbols are colored orange and blue, corresponding to data at 30 and 25 °C, respectively.

MOSCED 2005 > KBFF > TraPPE III > COSMO‐SAC

Figure 12. Equilibrium solubility predictions for acetaminophen in water using MOSCED with experimentally regressed parameters (MOSCED 2005) with the composition-dependent activity coefficient modeled using Wilson’s equation (○) and using (compositionindependent) infinite dilution activity coefficients (△).

Figure 13. Parity plot of the equilibrium aqueous solubility of acetaminophen predicted versus experimental data. Predictions are made using MOSCED with experimentally regressed parameters (MOSCED 2005 (○)), simulation regressed parameters using the TraPPE III (△) and KBFF (▽) reference sets, and using COSMOSAC (◇). The symbols are colored orange, blue, green, red and black corresponding to data obtained at temperatures of 30, 25, 20, 15, and 10 °C, respectively. L


Article


current suggestions are that all of the solvent force fields should be consistent with each other and should be compatible with the solute force field, and the solute force field should be the same for both the pure liquid solute and solute in solution calculations. We stress this last point. Our first attempt adopted values of the pure component fugacity from experimental vapor pressure data. However, changes in the pure liquid fugacity cause a shift in the calculated infinite dilution activity coefficients used to train MOSCED (see eq 4). This led to drastically different MOSCED parameters and poor predictions. The infinite dilution activity coefficient, as shown in eqs 4 and 5, directly compares the intermolecular solute−solvent and solute−solute interactions. Therefore, the solute should be modeled consistently. We found that we were able to make reasonable solubility predictions for acetaminophen in nonaqueous solvents using simulation-regressed solute MOSCED parameters. The models were able to well correlate solubility data over a wide range, with quantitative predictions superior to those of COSMOSAC using the VT Sigma Sigma Profile Databases. Comparing to 77 nonaqueous experimental solubilities at 5 temperatures for acetaminophen, we find that MOSCED parametrized using the TraPPE III reference set gives AAD and AAD% values of 0.0028 and 76.9, respectively. This dataset spans 5 orders of magnitude. The deviation is not far off from MOSCED 2005 parametrized using experimental solubility data, and it is comparable to errors exhibited for the experimental data itself. Lastly, we found that, when making solubility predictions in water, the infinite dilution (composition-independent) solute activity coefficient should be used. However, we reiterate that we find that the predictions are sensitive to small changes in the solvent force field and the set of reference solvents used to generate the reference data to parametrize MOSCED. While promising, further studies are required to better understand these challenges. Overall, using molecular simulation to develop a predictive form of MOSCED is promising. Here, we have focused on its use to predict the solubility of nonelectrolyte solids. However, a similar procedure could be used to parametrize MOSCED for use to model other types of phase equilibria, such as vapor− liquid and liquid−liquid. The ability to extend predictions to multicomponent systems by using MOSCED to parametrize an excess Gibbs free-energy model is attractive. Work is underway in our research group to investigate these possibilities and to additionally parametrize MOSCED using electronic structure methods, instead of using molecular simulation. In addition, while we have focused on regressing parameters for MOSCED, the method could just as well be applied to NRTL-SAC.

Quantitatively, for MOSCED 2005, TraPPE III, KBFF, and COSMO-SAC, respectively, we find that AAD% = {14.5, 261.9, 53.2, and 444.3} and AAD = {0.0002, 0.0040, 0.0008, and 0.0068}.112 Given the predictive nature of the MOSCED calculations trained with TraPPE III and KBFF, the results are very good. Again, note that errors in the experimental data themselves also exist. As shown in Figure 12 with the blue symbols, experimental data at 25 °C has a range of 0.0009, and comparing the two extrema, we would have an AAD% value of 90%. Importantly, in all cases, we are able to predict the same trend.

■

SUMMARY AND CONCLUSIONS The present study focuses on developing an efficient predictive method to calculate the equilibrium solubility of nonelectrolyte solids in a range of solvents at various temperatures useful for solvent selection during design processes. To accomplish this, we began with MOSCED. MOSCED is an efficient analytic method to compute infinite dilution activity coefficients over a range of temperatures that has been parametrized for 133 solvents. However, MOSCED is not predictive; rather, it requires reference infinite dilution activity coefficient data to obtain parameters for a compound of interest that has not previously been parametrized. This prevents MOSCED from being used with novel solutes. On the other hand, molecular simulation is a predictive tool that may be used to directly predict activity coefficients. Its use alone for design purposes would be computationally exhaustive. Here, we propose performing a limited number of molecular simulations to compute the infinite dilution activity coefficient in a range of solvents to obtain the reference data necessary to parametrize MOSCED. MOSCED may then be used to make predictions in other solvents and at other temperatures. In this fashion, we are able to leverage the strength of molecular simulation (pure predictions) with the efficiency of MOSCED. Infinite dilution activity coefficients alone are not sufficient when predicting the equilibrium solubility. Using MOSCED, we are able to predict both the infinite dilution activity coefficient of the solute in the solvent and for the solvent in the (subcooled) liquid solute. The subcooled liquid solute corresponds to a hypothetical state, and would require great care if one were to attempt to perform a molecular simulation to compute this activity coefficient directly. On the other hand, it can readily be computed using MOSCED, which allows one to parametrize an excess Gibbs free-energy model such as Wilson’s equation or UNIQUAC. With a parametrized excess Gibbs free-energy model, we readily have an analytic expression to predict the composition-dependent activity coefficient. Molecular simulations were performed to generate reference data to regress MOSCED parameters for the nonelectrolyte solid solutes acetanilide, acetaminophen, and phenacetin. Three unique reference sets were constructed using TraPPE-based solvents, and one reference set was constructed using the Kirkwood−Buff force field. Once parametrized, MOSCED was used to parametrize Wilson’s equation to capture the composition dependence of the activity coefficient. Computing ln( f1S/f10) using experimental data, equilibrium solubilities were then predicted. We first found that, despite their similarities, the solubility predictions for acetaminophen were sensitive to the TraPPE reference set used to regress MOSCED parameters. As more systems are studied in the future, we expect to gain a better understanding of how best to construct the reference set. Our

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.5b04807. Proof of the equivalence of the solvation free energy computed using molecular simulation free-energy calculations and the solute infinite dilution residual chemical potential and detailed discussion of the calculation of ln( f1S/f10) of the solutes (PDF) GROMACS force field files (type itp and type xyz) for the studied solutes for simulations with the TraPPE solvent (TXT) M


Article


■

GROMACS force field files (type itp and type xyz) for the studied solutes for simulations with the KBFF solvent (TXT)

Equilibria Properties of Complex Fluid Mixtures. Ind. Eng. Chem. Res. 2002, 41, 953−962. (15) Spyriouni, T.; Krokidis, X.; Economou, I. G. Thermodynamics of pharmaceuticals: Prediction of solubility in pure and mixed solvents with PC-SAFT. Fluid Phase Equilib. 2011, 302, 331−337. (16) Prausnitz, J. M.; Lichtenthaler, R. N.; de Azevedo, E. G. Molecular Thermodynamics of Fluid-Phase Equilibria, 2nd Edition; Prentice−Hall: Englewood Cliffs, NJ, 1986. (17) Ferrario, M.; Ciccotti, G.; Spohr, E.; Cartailler, T.; Turq, P. Solubility of KF in water by molecular dynamics using the Kirkwood integration method. J. Chem. Phys. 2002, 117, 4947−4953. (18) Sanz, E.; Vega, C. Solubility of KF and NaCl in water by molecular simulation. J. Chem. Phys. 2007, 126, 014507. (19) Paluch, A. S.; Jayaraman, S.; Shah, J. K.; Maginn, E. J. A method for computing the solubility limit of solids: Application to sodium chloride in water and alcohols. J. Chem. Phys. 2010, 133, 124504. (20) Moucka, F.; Lisal, M.; Skvor, J.; Jirsak, J.; Nezbeda, I.; Smith, W. R. Molecular Simulation of Aqueous Electrolyte Solubility. 2. Osmotic Ensemble Monte Carlo Methodology for Free Energy and Solubility Calculations and Application to NaCl. J. Phys. Chem. B 2011, 115, 7849−7861. (21) Schnieders, M. J.; Baltrusaitis, J.; Shi, Y.; Chattree, G.; Zheng, L. Q.; Yang, W.; Ren, P. Y. The Structure, Thermodynamics, and Solubility of Organic Crystals from Simulation with a Polarizable Force Field. J. Chem. Theory Comput. 2012, 8, 1721−1736. (22) Aragones, J. L.; Sanz, E.; Vega, C. Solubility of NaCl in water by molecular simulation revisited. J. Chem. Phys. 2012, 136, 244508. (23) Mester, Z.; Panagiotopoulos, A. Z. Mean ionic activity coefficients in aqueous NaCl solutions from molecular dynamics simulations. J. Chem. Phys. 2015, 142, 044507. (24) Paluch, A. S.; Parameswaran, S.; Liu, S.; Kolavennu, A.; Mobley, D. L. Predicting the excess solubility of acetanilide, acetaminophen, phenacetin, benzocaine, and caffeine in binary water/ethanol mixtures via molecular simulation. J. Chem. Phys. 2015, 142, 044508. (25) Paluch, A. S.; Maginn, E. J. Efficient Estimation of the Equilibrium Solution-Phase Fugacity of Soluble Nonelectrolyte Solids in Binary Solvents by Molecular Simulation. Ind. Eng. Chem. Res. 2013, 52, 13743−13760. (26) Liu, S.; Cao, S.; Hoang, K.; Young, K. L.; Paluch, A. S.; Mobley, D. L. Using MD Simulations To Calculate How Solvents Modulate Solubility. J. Chem. Theory Comput. 2016, 12, 1930−1941. (27) Mullins, E.; Oldland, R.; Liu, Y. A.; Wang, S.; Sandler, S. I.; Chen, C. C.; Zwolak, M.; Seavey, K. C. Sigma-profile database for using COSMO-based thermodynamic methods. Ind. Eng. Chem. Res. 2006, 45, 4389−4415. (28) Mullins, E.; Liu, Y. A.; Ghaderi, A.; Fast, S. D. Sigma profile database for predicting solid solubility in pure and mixed solvent mixtures for organic pharmacological compounds with COSMO-based thermodynamic methods. Ind. Eng. Chem. Res. 2008, 47, 1707−1725. (29) Lin, S. T.; Sandler, S. I. A Priori Phase Equilibrium Prediction from a Segment Contribution Solvation Model. Ind. Eng. Chem. Res. 2002, 41, 899−913. (30) Klamt, A.; Eckert, F.; Arlt, W. COSMO-RS: An Alternative to Simulation for Calculating Thermodynamic Properties of Liquid Mixtures. Annu. Rev. Chem. Biomol. Eng. 2010, 1, 101−122. (31) Klamt, A.; Eckert, F.; Hornig, M.; Beck, M. E.; Bürger, T. Prediction of Aqueous Solubility of Drugs and Pesticides with COSMO-RS. J. Comput. Chem. 2002, 23, 275−281. (32) Lucia, A. Multi-scale methods and complex processes: A survey and look ahead. Comput. Chem. Eng. 2010, 34, 1467−1475. (33) Lucia, A.; Bonk, B. M.; Waterman, R. R.; Roy, A. A multi-scale framework for multi-equilibrium flash. Comput. Chem. Eng. 2012, 36, 79−98. (34) Lucia, A.; Henley, H. Thermodynamic consistency of the multiscale Gibbs-Helmholtz constrained equation of state. Chem. Eng. Res. Des. 2013, 91, 1748. (35) Henley, H.; Thomas, E.; Lucia, A. Density and phase equilibrium for ice and structure I hydrates using the Gibbs−

AUTHOR INFORMATION

Corresponding Author

*Tel.: (513) 529-0784. Fax: (513) 529-0761. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS A.S.P. gratefully acknowledges funding from the Miami University Committee on Faculty Research. A.S.P. and R.T.L. are grateful for summer financial support from the Miami University College of Engineering and Computing. Computing support was provided by the Ohio Supercomputer Center and Miami University’s Research Computing Support group.

■

REFERENCES

(1) Jouyban, A. Handbook of Solubility Data for Pharmaceuticals; CRC Press/Taylor and Francis Group: Boca Raton, FL, 2010. (2) Poling, B. E.; Prausnitz, J. M.; O’Connell, J. P. The Properties of Gases and Liquids, 5th Edition; McGraw−Hill: New York, 2001. (3) Gmehling, J. Present status and potential of group contribution methods for process development. J. Chem. Thermodyn. 2009, 41, 731−747. (4) Gmehling, J.; Constantinescu, D.; Schmid, B. Group Contribution Methods for Phase Equilibrium Calculations. Annu. Rev. Chem. Biomol. Eng. 2015, 6, 267−292. (5) O’Connell, J. P.; Gani, R.; Mathias, P. M.; Maurer, G.; Olson, J. D.; Crafts, P. A. Thermodynamic Property Modeling for Chemical Process and Product Engineering: Some Perspectives. Ind. Eng. Chem. Res. 2009, 48, 4619−4637. (6) Gracin, S.; Brinck, T.; Rasmuson, A. C. Prediction of Solubility of Solid Organic Compounds in Solvents by UNIFAC. Ind. Eng. Chem. Res. 2002, 41, 5114−5124. (7) Hojjati, H.; Rohani, S. Measurement and Prediction of Solubility of Paracetamol in Water−Isopropanol Solution. Part 2. Prediction. Org. Process Res. Dev. 2006, 10, 1110−1118. (8) Diedrichs, A.; Gmehling, J. Solubility Calculation of Active Pharmaceutical Ingedients in Alkanes, Alcohols, Water and Their Mixtures Using Various Activity Coefficient Models. Ind. Eng. Chem. Res. 2011, 50, 1757−1769. (9) Schmid, B.; Gmehling, J. The universal group contribution equation of state VTPR present status and potential for process development. Fluid Phase Equilib. 2011, 302, 213−219. (10) Chen, C.; Crafts, P. A. Correlation and Prediction of Drug Molecule Solubility in Mixed Solvent Systems with the Nonrandom Two-Liquid Segment Activity Coefficient (NRTL-SAC) Model. Ind. Eng. Chem. Res. 2006, 45, 4816−4824. (11) Lazzaroni, M. J.; Bush, D.; Eckert, C. A.; Frank, T. C.; Gupta, S.; Olson, J. D. Revision of MOSCED Parameters and Extension to Solid Solubility Calculations. Ind. Eng. Chem. Res. 2005, 44, 4075−4083. (12) Draucker, L. C.; Janakat, M.; Lazzaroni, M. J.; Bush, D.; Eckert, C. A.; Frank, T. C.; Olson, J. D. Experimental determination and model prediction of solid solubility of multifunctional compounds in pure and mixed nonelectrolyte solvents. Ind. Eng. Chem. Res. 2007, 46, 2198−2204. (13) Cassens, J.; Ruether, F.; Leonhard, K.; Sadowski, G. Solubility Calculation of Pharmaceutical CompoundsA Priori Parameter Estimation Using Quantum-Chemistry. Fluid Phase Equilib. 2010, 299, 161−170. (14) Economou, I. G. Statistical Associating Fluid Theory: A Successful Model for the Calculation of Thermodynamic and Phase N


Article

Industrial & Engineering Chemistry Research Helmholtz constrained equation of state. Chem. Eng. Res. Des. 2014, 92, 2977−2991. (36) Kelly, R. B.; Lucia, A. On the linear approximation of mixture internal energies of departure. Comput. Chem. Eng. 2016, 85, 72−75. (37) Liu, Y. A. VT Sigma Profile Database; available via the Internet at: http://www.design.che.vt.edu/VT-Databases.html (accessed May 1, 2015). (38) Mullins, P. E. Application of COSMO-SAC to Solid Solubility in Pure and Mixed Solvent Mixtures for Organic Pharmacological Compounds; M.Sc. Thesis, Virginia Polytechnic Institute and State University: Blacksburg, VA, 2007. (39) Hildebrand, J. H.; Prausnitz, J. M.; Scott, R. L. Regular and Related Solutions; Van Nostrand Reinhold Company: New York, 1970. (40) Walas, S. M. Phase Equilibria in Chemical Engineering; Butterworth Publishers: Stoneham, MA, 1985. (41) Nordström, F. L.; Rasmuson, A. C. Determination of the activity of a molecular solute in saturated solution. J. Chem. Thermodyn. 2008, 40, 1684−1692. (42) Nordström, F. L.; Rasmuson, A. C. Prediction of solubility curves and melting properties of organic and pharmaceutical compounds. Eur. J. Pharm. Sci. 2009, 36, 330−344. (43) Yang, H.; Thati, J.; Rasmuson, A. C. Thermodynamics of molecular solids in organic solvents. J. Chem. Thermodyn. 2012, 48, 150−159. (44) Eike, D. M.; Maginn, E. J. Atomistic simulation of solid-liquid coexistence for molecular systems: Application to triazole and benzene. J. Chem. Phys. 2006, 124, 164503. (45) Jayaraman, S.; Maginn, E. J. Computing the melting point and thermodynamic stability of the orthorhombic and monoclinic crystalline polymorphs of the ionic liquid 1-n-butyl-3-methylimidazolium chloride. J. Chem. Phys. 2007, 127, 214504. (46) Recall that ln( f S2/f 02) = [1/(RT)](μS2 − μS2). (47) Thomas, E. R.; Eckert, C. A. Prediction of limiting activity coefficients by a modified separation of cohesive energy density model and UNIFAC. Ind. Eng. Chem. Process Des. Dev. 1984, 23, 194−209. (48) Frank, T. C.; Anderson, J. J.; Olson, J. D.; Eckert, C. A. Application of MOSCED and UNIFAC to screen hydrophobic solvents for extraction of hydrogen-bonding organics from aqueous solution. Ind. Eng. Chem. Res. 2007, 46, 4621−4625. (49) Paluch, A. S.; Maginn, E. J. Predicting the Solubility of Solid Phenanthrene: A Combined Molecular Simulation and Group Contribution Approach. AIChE J. 2013, 59, 2647−2661. (50) The relation between the solvation free energy computed using molecular simulation and the residual chemical potential is provided in the section entitled “Relating Solvation Free Energy to Residual Chemical Potential” in the Supporting Information. (51) Fuerst, G. B.; Ley, R. T.; Paluch, A. S. Calculating the Fugacity of Pure, Low Volatile Liquids via Molecular Simulation with Application to Acetanilide, Acetaminophen, and Phenacetin. Ind. Eng. Chem. Res. 2015, 54, 9027−9037. (52) Leach, A. R. Molecular Modelling: Principles and Applications, 2nd Edition; Pearson Education Limited: Harlow, England, 2001. (53) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications, 2nd Edition; Academic Press: San Diego, CA, 2002. (54) Rai, N.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 9. Explicit Hydrogen Description of Benzene and FiveMembered and Six-Membered Heterocyclic Aromatic Compounds. J. Phys. Chem. B 2007, 111, 10790−10799. (55) Rai, N.; Bhatt, D.; Siepmann, J. I.; Fried, L. E. Monte Carlo simulations of 1,3,5-triamino-2,4,6-trinitrobenzene (TATB): Pressure and temperature effects for the solid phase and vapor−liquid phase equilibria. J. Chem. Phys. 2008, 129, 194510. (56) Rai, N.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 10. Explicit-Hydrogen Description of Substituted Benzenes and Polycyclic Aromatic Compounds. J. Phys. Chem. B 2013, 117, 273−288. (57) Wick, C. D.; Stubbs, J. M.; Rai, N.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 7. Primary, Secondary, and Tertiary

Amines, Nitroalkanes and Nitrobenzene, Nitriles, Amides, Pyridine, and Pyrimidine. J. Phys. Chem. B 2005, 109, 18974−18982. (58) Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 1. United-Atom Description of n-Alkanes. J. Phys. Chem. B 1998, 102, 2569−2577. (59) Stubbs, J. M.; Potoff, J. J.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 6. United-Atom Description for Ethers, Glycols, Ketones, and Aldehydes. J. Phys. Chem. B 2004, 108, 17596−17605. (60) Martin, M. G.; Siepmann, J. I. Novel Configurational-Bias Monte Carlo Method for Branced Molecules. Transferable Potentials for Phase Equilibria. 2. United-Atom Description of Branched Alkanes. J. Phys. Chem. B 1999, 103, 4508−4517. (61) Chen, B.; Potoff, J. J.; Siepmann, J. I. Monte Carlo Calculations for Alcohols and Their Mixtures with Alkanes. Transferable Potentials for Phase Equilibria. 5. United-Atom Description of Primary, Secondary, and Tertiary Alcohols. J. Phys. Chem. B 2001, 105, 3093−3104. (62) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (63) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926. (64) Chen, B.; Siepmann, J. I. Microscopic structure and solvation in dry and wet octanol. J. Phys. Chem. B 2006, 110, 3555−3563. (65) Rafferty, J. L.; Sun, L.; Siepmann, J. I.; Schure, M. R. Investigation of the driving forces for retention in reversed-phase liquid chromatography: Monte Carlo simulations of solute partitioning between n-hexadecane and various aqueous-organic mixtures. Fluid Phase Equilib. 2010, 290, 25−35. (66) Bai, P.; Siepmann, J. I. Gibbs ensemble Monte Carlo simulations for the liquid-liquid phase equilibria of dipropylene glycol dimethyl ether and water: A preliminary report. Fluid Phase Equilib. 2011, 310, 11−18. (67) Kamath, G.; Georgiev, G.; Potoff, J. J. Molecular Modeling of Phase Behavior and Microstructure of Acetone−Chloroform− Methanol Binary Mixtures. J. Phys. Chem. B 2005, 109, 19463−19473. (68) The acetone model of Potoff and co-workers differs from TraPPE-UA in that the partial charge values and the carbonyl carbon (CO) LJ ε value are different. For extension to butanone, we took TraPPE-UA butanone and changed the partial charge distribution and carbonyl carbon LJ ε to agree with the modifications of Potoff and coworkers. Also, from personal communication with Jeff Potoff, the chloroform carbon LJ ε was changed from ε/kB = 68.94 to ε/kB = 53.0. (69) Weerasinghe, S.; Smith, P. E. A Kirkwood−Buff Derived Force Field for Methanol and Aqueous Methanol Solutions. J. Phys. Chem. B 2005, 109, 15080−15086. (70) Weerasinghe, S.; Smith, P. E. Kirkwood−Buff derived force field for mixtures of acetone and water. J. Chem. Phys. 2003, 118, 10663. (71) Ploetz, A. E.; Bentenitis, N.; Smith, P. E. Developing force fields from the microscopic structure of solutions. Fluid Phase Equilib. 2010, 290, 43−47. (72) Ploetz, A. E.; Smith, P. E. A Kirkwood−Buff force field for the aromatic amino acids. Phys. Chem. Chem. Phys. 2011, 13, 18154− 18167. (73) Schuler, L. D.; Daura, X.; van Gunsteren, W. F. An Improved GROMOS96 Force Field for Aliphatic Hydrocarbons in the Condensed Phase. J. Comput. Chem. 2001, 22, 1205−1218. (74) Chen, F. Molecular dynamics simulations of solution mixtures and solution/vapor interfaces; Ph.D. Thesis, Kansas State University: Manhattan, KS, 2003. (75) The KBFF methanol model adopts intramolecular parameters from OPLS-UA, which were also adopted by TraPPE-UA.61 Therefore, the intramolecular parameters for the other alcohols should be consistent. (76) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. O


Article

Industrial & Engineering Chemistry Research (77) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (78) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: A high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 2013, 29, 845−854. (79) GROMACS: Fast, Flexible, Free; available via the Internet at: http://www.gromacs.org/ (accessed Aug. 1, 2013). (80) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. Packmol: A package for building initial configurations for molecular dynamics simulations. J. Comput. Chem. 2009, 30, 2157−2164. (81) Packmol: Packing Optimization for Molecular Dynamics Simulations; available via the Internet at: http://www.ime.unicamp. br/~martinez/packmol/ (accessed Aug. 1, 2013). (82) van der Spoel, D.; Lindahl, E.; Hess, B. (The GROMACS Development Team). GROMACS User Manual, version 4.6.3; available via the Internet at: ftp://ftp.gromacs.org/pub/manual/manual-4.6.3. pdf, 2013. (83) van Gunsteren, W. F.; Berendsen, H. J. C. A leap-frog algorithm for stochastic dynamics. Mol. Simul. 1988, 1, 173−185. (84) Berendsen, H. J. C. Simulating the Physical World: Hierarchial Modeling from Quantum Mechanics to Fluid Dynamics; Cambridge University Press: New York, 2007. (85) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684−3690. (86) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182−7190. (87) Nosé, S.; Klein, M. L. Constant pressure molecular dynamics for molecular systems. Mol. Phys. 1983, 50, 1055−1076. (88) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for molecular simulations. J. Comput. Chem. 1997, 18, 1463−1472. (89) Hess, B. P-LINCS: A Parallel Linear Constraint Solver for molecular simulation. J. Chem. Theory Comput. 2008, 4, 116−122. (90) Deserno, M.; Holm, C. How to mesh up Ewald sums. I. A theoretical and numerical comparison of various particle mesh routines. J. Chem. Phys. 1998, 109, 7678−7693. (91) The alkanes n-hexane and 2,5-dimethylhexane were modeled with united atom force fields and contain no partial charges. For these pure solvent systems, use of SPME for long-range electrostatic interactions is not necessary. When a solute is added, the solute contains partial charges so use of SPME is necessary. When performing free-energy calculations, since there are no solvent−solute electrostatic interactions, stages to decouple the electrostatic interactions are not necessary. (92) Shing, K. S.; Chung, S. T. Computer Simulation Methods for the Calculation of Solubility in Supercritical Extraction Systems. J. Phys. Chem. 1987, 91, 1674−1681. (93) Kofke, D. A.; Cummings, P. T. Quantitative comparison and optimization of methods for evaluating the chemical potential by molecular simulation. Mol. Phys. 1997, 92, 973−996. (94) Shirts, M. R.; Pitera, J. W.; Swope, W. C.; Pande, V. S. Extremely precise free energy calculations of amino acid side chain analogs: Comparison of common molecular mechanics force fields for proteins. J. Chem. Phys. 2003, 119, 5740−5761. (95) Kofke, D. A.; Cummings, P. T. Precision and accuracy of staged free-energy perturbation methods for computing the chemical potential by molecular simulation. Fluid Phase Equilib. 1998, 150− 151, 41−49. (96) Chipot, C., Pohorille, A., Eds. Free Energy Calculations: Theory and Applications in Chemistry and Biology; Springer Series in Chemical Physics, Vol. 86; Springer: New York, 2007. (97) Bennett, C. H. Efficient Estimation of Free Energy Differences from Monte Carlo Data. J. Comput. Phys. 1976, 22, 245−268.

(98) Shirts, M. R.; Bair, E.; Hooker, G.; Pande, V. S. Equilibrium Free Energies from Nonequilibrium Measurements Using MaximumLikelihood Methods. Phys. Rev. Lett. 2003, 91, 140601. (99) Lu, N.; Singh, J. K.; Kofke, D. A. Appropriate methods to combine forward and reverse free-energy perturbation averages. J. Chem. Phys. 2003, 118, 2977−2984. (100) Shirts, M. R.; Chodera, J. D. Statistically optimal analysis of samples from multiple equilibrium states. J. Chem. Phys. 2008, 129, 124105. (101) Beutler, T. C.; Mark, A. E.; van Schaik, R. C.; Gerber, P. R.; van Gunsteren, W. F. Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations. Chem. Phys. Lett. 1994, 222, 529−539. (102) Shirts, M. R.; Pande, V. S. Solvation free energies of amino acid side chain analogs for common molecular mechanics water models. J. Chem. Phys. 2005, 122, 134508. (103) Steinbrecher, T.; Mobley, D. L.; Case, D. A. Nonlinear scaling schemes for Lennard-Jones interactions in free energy calculations. J. Chem. Phys. 2007, 127, 214108. (104) PyMBAR: Python implementation of the multistate Bennett acceptance ratio (MBAR); available via the Internet at: https://github. com/choderalab/pymbar (accessed May 1, 2014). (105) Chodera, J. D.; Swope, W. C.; Pitera, J. W.; Seok, C.; Dill, K. A. Use of the Weighted Histogram Analysis Method for the Analysis of Simulated and Parallel Tempering Simulations. J. Chem. Theory Comput. 2007, 3, 26−41. (106) Klimovich, P. V.; Shirts, M. R.; Mobley, D. L. Guidelines for analysis of free energy calculations. J. Comput.-Aided Mol. Des. 2015, 29, 397−411. (107) Storn, R.; Price, K. Differential EvolutionA Simple and Efficient Heuristic for Global Optimization over Continuous Spaces. J. Global. Optim. 1997, 11, 341−359. (108) Eaton, J. W.; Bateman, D.; Hauberg, S. GNU Octave Version 3.0.1 Manual: A High-Level Interactive Language for Numerical Computations; CreateSpace Independent Publishing Platform, 2009 (ISBN No. 1441413006). (109) Marrero, J., Abildskov, J., Eds. Solubility and Related Properties of Large Complex Chemicals, Part 1: Organic Solutes ranging from C4 to C40; DECHEMA: Frankfurt am Main, Germany, 2003. (110) Abildskov, J., Ed. Solubility and Related Properties of Large Complex Chemicals, Part 2: Organic Solutes ranging from C2 to C41; DECHEMA: Frankfurt am Main, Germany, 2005. (111) Mu, T.; Rarey, J.; Gmehling, J. Performance of COSMO-RS with Sigma Profiles from Different Model Chemistries. Ind. Eng. Chem. Res. 2007, 46, 6612−6629. (112) Note that Mullins et al.28 considered three different conformations (a, b, and c) of acetaminophen when generating its sigma profile. The results presented here are for conformation c, which we found to perform best. For conformations a, b, and c, we find nonaqueous solubility predictions of AAE% = 1425.9, 104.2, and 89.3, respectively, and AAE = 0.1590, 0.1561, and 0.1354, respectively. For aqueous predictions, we have AAE% = 567.2, 708.6, and 444.3, respectively, and AAE = 0.0087, 0.0109, and 0.0068, respectively. (113) The experimental data are all of the aqueous solubility data from DECHEMA,109,110,113 with one exception. At 30 °C, they report three values of the mole fraction solubility, one of which is 2.5981 × 10−5. This is 2 orders of magnitude smaller than all of the other data reported at 30 °C or any other temperature. Therefore, we suspect an error with this datum and have excluded it.

P


Developing a Predictive Form of MOSCED for Nonelectrolyte Solids

Recommend Documents