Estimating Octanol–Water Partition Coefficients for ... - ACS Publications

The use of the conductor-like screening with segment activity coefficient (COSMO-SAC) model to predict octanol–water partition coefficients (log KO/...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/IECR

Estimating Octanol−Water Partition Coefficients for Selected Nanoscale Building Blocks Using the COSMO-SAC Segment Contribution Method Patrick S. Redmill* Department of Chemical Engineering, Vanderbilt University, Nashville, Tennessee 37235, United States ABSTRACT: The use of the conductor-like screening with segment activity coefficient (COSMO-SAC) model to predict octanol−water partition coefficients (log KO/W) for selected nanoscale building blocks (NBBs) has been evaluated. COSMO calculations have been performed using charge density information obtained from density functional theory calculations of the compounds that, in turn, provide input data that allow the calculation of the activity coefficient in the COSMO-SAC model. The NBBs of interest in this study are C60, the hydrogen functionalized polyhedral oligomeric silsesquioxane Si8O12H8 (H-POSS), and their functionalized variants C60(OH)32, Si8O12F8 (F-POSS), and Si8O12(OH)8 (OH-POSS). It is found that COSMO-SAC, while being unable to accurately predict the solubility of very hydrophobic compounds in water, can quickly and efficiently provide qualitative estimates of log KO/W for the selected solutes. The COSMO-SAC estimation of log KO/W for C60 is 6.79, which is consistent with the experimentally reported value and reflects the general consensus that C60 is strongly hydrophobic. No partitioning data are available for H-POSS or the functionalized NBBs studied. In this work, COSMO-SAC calculations on these particles indicate shifts toward water partitioning for C60(OH)32, F-POSS, and OH-POSS relative to their respective base particles (C60 and H-POSS). These shifts in water solubility are consistent with the calculation of the Gibbs free energy of hydration for each particle using molecular dynamics simulations and thermodynamic integration.

1. INTRODUCTION With the burgeoning interest in nanoscience,1 there is a growing need for solubility related parameters for nanoparticles and nanoscale building blocks in common solvents, such as water and organic compounds. Such parameters can describe the partitioning behavior of the solute in an equilibrium aqueous/organic environment, a characteristic that plays an important role in fields such as toxicology and pharmacology. In particular, accurate descriptions of the infinite dilution activity coefficient (γU∞) of these particles in water and octanol solvents allow for the determination of KO/W, which is commonly considered as an indicator of the tendency for a particle to absorb into sedimentary or tissue systems. For a given solute, labeled U, KO/W is defined by assuming that U is present in two immiscible liquid phases, one octanol-rich, the other water-rich, in equilibrium with each other. If the concentration of U in the octanol-rich phase is CUoctanol and the concentration of U in the water-rich phase is CUwater, then: KO/W =

octanol CU

water = CU

octanol-rich phase can be approximated using the following equation:3 ρmix w 1 ρ̂ V = ; =∑ i ∑i xiMi ρmix ρi (2) i

where ρmix is the mixture mass density, xi is the mole fraction of solvent i, Mi is the molecular weight of solvent i, ρi is the pure component mass density of solvent i, and wi is the mass fraction of solvent i. Many empirical methods, such as the Universal Functional Activity Coefficient (UNIFAC) method,4 quantitative structure−activity relationship (QSAR) approaches,5 and neural network (e.g., ALOGPs 6−8) approaches9−11 have been developed for the estimation of KO/W either directly or indirectly through the calculation of activity coefficients. However, each of these methods relies on the availability of phase equilibria data for the components of interest, or at least phase equilibria data of components that share common functional groups with the components of interest. Given the structural and chemical novelty of many nanoparticles, these data are often unavailable. Alternatively, some molecularsimulation-based techniques, such as thermodynamic integration or Widom particle insertion,12 have been developed for the calculation of the Gibbs free energy of solvation (ΔGsolv) of a solute. Typically, ΔGsolv is used to calculate the Henry’s law

octanol ρ̂octanol XU water ρ̂water xU

(1)

xUV

where is the mole fraction of solute U in the solvent, denoted generally as V, and ρ̂V is the molar density of solvent V. The “wet” KO/W, where both the water and the octanol phases are mutually saturated, is the generally accepted form of the octanol−water partition coefficient. In an equilibrium water−octanol system, the water-rich phase is effectively devoid of octanol; however, the octanol-rich phase contains about 23 mol % water.2 The molar density of the multicomponent © 2012 American Chemical Society

Received: Revised: Accepted: Published: 4556

September 14, 2011 February 29, 2012 March 8, 2012 March 8, 2012 dx.doi.org/10.1021/ie202107t | Ind. Eng. Chem. Res. 2012, 51, 4556−4566

Industrial & Engineering Chemistry Research

Article

COSMO-based solvation methods use the screening charges on the surface of the solvent-accessible area of a given molecule in a perfect conductor, which can be obtained using modest quantum mechanical calculations. In turn, charge densities on the surface of the molecules can be calculated, and the interaction of these charge densities between the solute and the solvent can be quantified, which allows for the prediction of activity coefficients. The ability to use an activity coefficient model to calculate various solvation parameters, such as γU∞ and log KO/W, would greatly assist in quickly and inexpensively ascertaining the potential environmental impact of a large array of NBBs and their increasing number of functionalized variants. The use of a segment contribution activity coefficient model, such as COSMO-SAC, is particularly suited for the calculation of log KO/W for novel NBBs, as no system specific phase equilibria data are needed for reliable results. Furthermore, the use of a segment contribution activity coefficient model circumvents the limitation inherent in a molecular dynamics/thermodynamic integration approach, which requires a sparingly soluble solute in a given solvent. In this study, the applicability of COSMO-SAC toward the prediction of log KO/W for NBBs is examined by comparing the partitioning shifts upon particle functionalization, as calculated by COSMO-SAC, to shifts in the Gibbs free energy of solvation in water (ΔGhyd) as calculated by thermodynamic integration. The ability to perform log KO/W calculations on NBB systems using activity coefficient models is beneficial because it has been shown that ΔGUsolv,octanol for moderate to large NBBs is very low.30 As a result, many NBBs, such as those in this study, are likely to fail the sparingly soluble requirement of eqs 3−5, and the application of an activity coefficient model will be necessary to estimate log KO/W. The comparison between the COSMOSAC log KO/W and ΔGhyd is significant, as the low density of polar groups in the octanol phase will result in small changes in both cavity formation free energy and enhanced solute−solvent interactions upon NBB functionalization. Therefore, changes in water solubility at infinite dilution will directly drive changes in log KO/W. The COSMO-SAC calculation of log KO/W for C60 gives results that reflect the well-documented extreme hydrophobicity of C60. Additionally, COSMO-SAC calculations for H-POSS, F-POSS, OH-POSS, and C60(OH)32 predict hydrophilic shifts of varying degrees upon the addition of functional groups to the base particles. The COSMO-SAC results are corroborated by ΔGhyd values of C60(OH)32, H-POSS, OHPOSS, and F-POSS that were calculated using thermodynamic integration in this work.

coefficient, which is a common solubility descriptor, using the relationship defined by Ben-Naim.13 ⎡ ΔGsolv ⎤ ⎥ HU,V = ρ̂RT exp⎢ ⎢⎣ RT ⎥⎦

(3)

where HU,V is the Henry’s law coefficient of solute U in solvent V, R is the gas constant, and T is the system temperature. The fugacity of U dissolved in V, f U,V, can in turn be calculated using the following relationship. fU,V = lim HU,Vx UV x UV → 0

(4)

xUV

where is the mole fraction of solute U in solvent V. Assuming that the octanol and water phases are in equilibria, the fugacity of solute U in each phase can be equated, and the ratio of mole fractions seen in eq 1 can be calculated. octanol xU

HU,water water = H xU U,octanol

(5)

Such techniques can often yield very accurate solubilities with relatively little phase equilibria data required. However, these approaches may be limited due to the size of the particle, lack of relevant force-fields, and computational expense. Additionally, the application of the Henry’s law coefficient for the calculation of the fugacity is limited to conditions of a sparingly soluble solute, where solute−solute interactions become negligible. In cases where ΔGsolv indicates significant solute relative solubility, alternative methodologies that take into account the solute−solute interactions must be used to accurately predict partition coefficients. Traditional activity coefficient models, such as UNIFAC, have terms that account for solute−solute interactions and are often required when calculating the partitioning of solutes of significant solubility;14 the parameters in such activity coefficient models are either fitted to experimental data or calculated from molecular attributes found in quantum mechanical computations. Although somewhat less effective than models that are fitted to experimental data,15 COSMObased activity coefficient models have shown remarkable accuracy for the prediction of γU∞16,17 and log KO/W18−22 for a variety of solutes, which exhibit significant solubility in the given solvent. From a practical standpoint, COSMO-based solvation methods, such as COSMO-RS23,24 and COSMOSAC,25 seem to be an excellent choice, as the calculations take only moderate computational effort, and no specific phase equilibria data are needed to fit energetic parameters for functional groups. COSMO-based solvation methods are not without their limitations. One such limitation is the tendency of COSMObased methods to underpredict γU∞,water for hydrophobic particles.26−28 Therefore, the predictions of γU∞,water for extremely hydrophobic solutes, such as C60, should be considered qualitative, and the resulting log KO/W values should be regarded as a lower bound. Furthermore, COSMO-based methods typically only consider a single configuration for a given molecule. It has been shown that COSMO-SAC produces appreciably different values of γ for conformational isomers of some molecules.29 Therefore, when applying COSMO-SAC to molecules that have multiple stable configurations, a somewhat larger deviation from experimental results should be anticipated.

2. METHODOLOGY 2.1. The COSMO-SAC Model. Using liquid phase equilibria expressions, the KO/W of a component can be estimated if γU∞ is known for both phases, as shown below. Starting with the standard equation for liquid−liquid equilibria and focusing on the solute U, if the mole fraction of U in the octanol-rich liquid phase is xUoctanol and the mole fraction of the solute in the water-rich phase is xUwater, then equilibrium requires that octanolf 0 = γwaterx waterf 0 γoctanol xU U U U U U

(6)

f U0 ,

where the standard state fugacity of the solute, is the same for both phases, and γUoctanol and γUwater are the activity coefficients of the solute in the octanol-rich and water-rich phases, 4557

dx.doi.org/10.1021/ie202107t | Ind. Eng. Chem. Res. 2012, 51, 4556−4566

Industrial & Engineering Chemistry Research

Article

calculation divides the solvent-accessible surface area into ni segments, and the probability density is then expressed as:

respectively. KO/W can be calculated by rearranging eq 6 and multiplying by bulk solvent properties: KO/W =

octanol CU

n (σ) A (σ) pi (σ) = i = i ni Ai

water CU x octanol ρ̂ = Uwater octanol ρ̂water xU water γ ρ̂octanol = U ρ̂water γoctanol U

where ni is the number of segments whose charge density is σ and Ai(σ) is the area that corresponds to a charge density of σ. The segment activity coefficient, Γs(σ), is defined by the following expression. ⎧ ⎫ ⎡ −ΔW (σm, σn) ⎤⎪ ⎪ lnΓs(σm) = −ln⎨∑ ps (σn)Γs(σn)exp⎢ ⎥⎬ ⎣ ⎦⎪ kT ⎪σ ⎭ ⎩ n

(7)

The water-rich phase can be assumed to contain no octanol; however, it is generally accepted that n-octanol in equilibrium with water contains about 23 mol % water.2 If one assumes that the solute is at low dilution in each phase, eq 7 can be approximated by KO/W =

(15)

Equation 15 requires the calculation of an aggregate charge density probability, ps(σ), for the mixture, which is defined as

,water γ∞ ρ̂octanol U ,octanol ρ̂ γ∞ water U

ps (σ) =

(8)

where γUoctanol and γU∞,water are the activity coefficients of the solute at infinite dilution in the octanol-rich and water-rich phases, respectively. The COSMO-SAC model expresses the activity coefficient of a chemical species i in a mixture using a linear combination of two contributions, the combinatorial contribution (γi C) and the residual contribution (γiR):

li =

∑j qjxj

∑i xiA i

(16)

0, σacc − σHB]min[0, σdon + σHB]

where α′, cHB, and σHB are constants defined by the developers of COSMO-RS,23,34,35 and σHB and σdon are the larger and smaller values of σm and σn, respectively. If one uses the resulting charge densities from COSMO calculations and the constants in eq 17, eq 15 can be solved iteratively, and, in turn, the residual portion of the activity coefficient can be calculated using eq 13. 2.2. Thermodynamic Integration Technique. One major hurdle in establishing COSMO-based activity coefficient models as an effective tool for calculating KO/W parameters for NBB systems is the general lack of solubility data available for comparison. Some work correlating the solubility of C60 using QSPR36,37 has been published, along with some experimental partitioning data of the C60 and C70 fullerenes in water38 and C60 in organic solvents.39 However, the availability of solubility parameters for novel NBBs, such as POSS or functionalized fullerenes, is very limited. To verify that predictive techniques can adequately model systems of NBBs at infinite dilution, an alternate means of ascertaining the relative solubility is desirable. In this work, ΔGhyd values for the hydoxylated fullerene C60(OH)32 and three POSS variants, H-POSS, FPOSS, and OH-POSS, have been calculated in a water solvent using molecular dynamics simulations and thermodynamic integration (TI).12 Figure 1 shows the NBB variants studied. The simulations were performed in the isothermal−isobaric (NPT) ensemble, in which the number of particles (N), the system pressure (P), and the system temperature (T) are constant. To avoid periodic effects in the molecular dynamics simulations, each system is of sufficient size, consisting of one solute particle constrained to the center of the simulation cell and 3334 water molecules. Because many experimental KO/W studies are carried out at ambient conditions, the pressure and temperature of the simulated systems were held at 1 bar and

∑ xjl j j

rixi ∑j rjxj

z (ri − qi) − (ri − 1) 2

(10)

(11)

(12)

In these equations, xi is the mole fraction of species i, z is a constant commonly taken as 10, and ri and qi are the normalized molecular volume and surface area for molecule i. The residual contribution is due to the charging free energy, which arises due to the restoration of solute charges in the solvent. COSMO-SAC uses the form of Lin and Sandler to estimate the residual effects: lnγiR = ni ∑ p(σm)[lnΓs(σm) − lnΓi(σm)] σm

∑i xiA i pi (σ)

ΔW (σm, σn) = (α′/2)(σm + σn)2 + cHB max[

(9)

Φi Φi z θi lnγC i = ln x + 2 ln Φ + li − x i i i

; Φi =

∑i xini

=

where xi is the mole fraction of molecule i in the system. The exchange energy between segments m and n, ΔW(σm,σn), is calculated from:

The combinatorial contribution is due to the so-called cavity formation free energy, which is the energy that would be required to insert the uncharged solute into the solvent. The COSMO-SAC method uses the Staverman−Guggenheim equation31,32 to account for combinatorial effects.

qixi

∑i xinipi (σ)

(17)

R lnγi = lnγC i + lnγi

θi =

(14)

(13)

where p(σ) is the probability of finding charge density σ(q/A2) on the surface of molecule i. The p(σ) profile can be generated from a COSMO33 calculation performed on a molecule using electronic structure information that can be obtained from a quantum mechanical point energy calculation. The COSMO 4558

dx.doi.org/10.1021/ie202107t | Ind. Eng. Chem. Res. 2012, 51, 4556−4566

Industrial & Engineering Chemistry Research

Article

simulations. When cross terms were unavailable, Lorentz− Berthelot mixing rules were employed: σi + σj σij = ; εij = εiεj (20) 2 where σi and εi are the vdW radius and well depth, respectively. It has been shown that allowing conformational changes in the solute can lead to significant noise in the thermodynamic integration calculations;48 therefore, all solutes in this study have been fixed rigidly in the center of the simulation cell. The lowest, in vacuo, energy configuration for each solute was found using a DFT level optimization with a minimum 6-31G basis set and a B3LYP exchange/correlation functional. With respect to the C60(OH)32 particle, which is believed to exhibit several hydroxyl adsorption patterns on its surface, the lowest energy C60(OH)32 adsorption pattern, calculated by Rodriguez-Zavala et al.,49 was used for the initial configuration for the DFT optimization. Most of the solutes in this study have no corresponding electrostatic force-field parameters. The DFT optimized structures used in the thermodynamic integration analysis are also used to obtain partial atomic charges by performing a Mulliken population analysis on a DFT/6-311G** point energy calculation. There is a considerable amount of debate regarding the validity of using the Mulliken population analysis50 for the calculation of partial atomic charges instead of other methods,51 which typically base partial charges on the fitting of electrostatic surfaces, such as ChelpG52 and Merz− Kollman53 schemes. However, these methods often fail for large molecules, as they often predict unrealistic charges for atoms that lie far from the electrostatic potential surface (ESP).54 Given the unrealistic charges computed in this work for the molecules of interest using ESP techniques, it seems likely that such methods are inadequate due to the inability to fit separate inner and outer surfaces of cage-like particles; hence, the particles in this work are calculated using the Mulliken population analysis. The integration used in eq 19 was performed in two steps. First, the vdW interactions between the solute and the solvent were grown in, with no solute partial charges, using a modified 12-6 potential:55

Figure 1. Nanoparticles compared in TI and COSMO studies. (a) HPOSS, (b) F-POSS, (c) OH-POSS, (d) C60(OH)32.

298 K, respectively. In a NPT system, the ΔGhyd can be expressed as follows: ⎛ ∂G ⎞ ⎜ ⎟ = ⎝ ∂λ ⎠ N , P , T

∂U ∂λ λ

(18)

so that ΔG hyd =

1 ⎛ ∂G ⎞

∫0





⎝ ∂λ ⎠ N , P , T

dλ =

1

∫0

∂U dλ ∂λ λ

(19)

where U is the potential between the solute and the solvent and λ is an arbitrarily defined variable such that λ = 0 represents a pure solvent state and λ = 1 represents a fully solvated NBB in solution. All intermediate values of λ indicate a “ghost” state for a solute particle, in which the solute−solvent interactions are softened. Therefore, ΔGhyd can be calculated given an explicit form of the solvent−solute potential, which should decay to 0 as λ → 0 and reduce to the standard potential form at λ = 0, and the resulting MD trajectories. When possible, with regard to the form of the solvent−solute potential, a CHARMM compatible force-field was used to model the systems studied.40 The CHARMM force-field was chosen due to the fact that it has been successfully applied to a wide array of compounds. Additionally, CHARMM is often considered the standard force-field in terms of phase equilibria calculations at ambient temperature and pressure,41 which are the conditions of this work. In all simulations, the TIP3P forcefield42 was used to model water, as it is the only water potential optimized for CHARMM. For systems with a POSS solute, a CHARMM compatible force-field for silica/water systems43 was used for the inner cage. For the functionalized groups in POSS, the CHARMM fluorocarbon potential was used to model van der Waals interactions (vdW) of the fluorine group,44 and the CHARMM22 force-field for lipids40 was used to model vdW interactions of the functional groups for HPOSS and OH-POSS. Because there is no CHARMM parametrization for graphitic materials, the potential of Jiang and Sandler,45 which utilizes the Bojan/Steel force-field for graphene atoms46 and the TraPPE force-field47 for hydroxyl groups, was used for the application to the C60(OH)32

⎛ ⎜ ⎜ ⎜ 1 Uij(r , λ1) = λ1n4εij⎜ n ⎜⎡ ⎛ r ⎞6⎤ ⎜ ⎢αLJ(1 − λ12) + ⎜ ⎟ ⎥ ⎜⎢ ⎝ σij ⎠ ⎥⎦ ⎝⎣ ⎞ ⎟ ⎟ ⎟ 1 − n⎟ ⎛r ⎞ ⎟ αLJ(1 − λ12) + ⎜ ⎟ ⎟ ⎝ σij ⎠ ⎟ ⎠

(21)

where r is the distance between particles i and j, λ1 is the van der Waals softening parameter defined before, αLJ is a positive constant, chosen to be 0.25, and n is an integer chosen to be 2. This particular modification was chosen, as it gives finite energies at low distances for soft potentials, minimizing 4559

dx.doi.org/10.1021/ie202107t | Ind. Eng. Chem. Res. 2012, 51, 4556−4566

Industrial & Engineering Chemistry Research

Article

Table 1. Infinite Dilution Activity Coefficients Calculated by COSMO-SAC for Various Solvents with Resulting log KO/Wa γU∞,water

γU∞,octanol

Alkanes propane 7.86 × 102 pentane 6.88 × 103 octane 1.41 × 105 undecane 2.33 × 106 Alcohols methanol 2.50 propanol 1.41 × 101 pentanol 9.09 × 101 heptanol 5.58 × 102 Aromatics benzene 3.05 × 102 toluene 8.65 × 102 o-xylene 1.98 × 103 m-xylene 2.43 × 103 p-xylene 2.58 × 103 Aromatic Alcohols phenol 4.67 o-cresol 6.92 × 101 m-cresol 1.45 × 101 b benzenetetrol 3.74 × 10−4 benzenehexolb 2.30 × 10−6 Fluorocarbons fluorobenzene 4.06 × 102 perfluorobenzene 9.96 × 103 perfluoroethane 2.75 × 103 vinyl fluoride 4.95 × 101 Polyaromatic Hydrocarbons kekuleneb 5.38 × 104 coroneneb 2.14 × 104 chryseneb 5.00 × 103 b anthracene 2.95 × 103

log KO/W

log KO/W exp.

% error from exp.

ALOGPS

AC_logP

1.05 1.27 1.43 1.45

2.03 2.89 4.15 5.36

2.36 3.39 5.18 −

14.0 14.8 19.9 −

2.19 3.41 4.73 5.96

1.84 2.77 4.16 5.55

0.724 0.542 0.461 0.400

−0.306 0.571 1.45 2.30

−0.770 0.250 1.51 2.62

60.2 128 3.94 12.2

−1.38 0.210 1.47 2.53

−0.01 0.890 1.82 2.74

1.18 1.14 1.14 1.12 1.11

1.57 2.04 2.40 2.49 2.52

2.13 2.73 3.12 3.20 3.15

26.4 25.4 23.2 22.1 19.9

2.03 2.56 3.16 3.15 3.15

1.98 2.30 2.61 2.61 2.61

1.23 1.53 1.65 0.184 0.327

1.46 1.95 1.96 − −

15.8 21.7 16.0 − −

1.39 1.89 1.93 0.060 0.18

1.68 2.00 2.00 0.790 0.19

1.06 1.09 1.13 0.859

1.74 3.12 2.54 0.916

2.27 2.55 2.00 −

23.4 22.2 27.1 −

2.18 2.33 2.46 1.02

2.04 2.34 2.34 1.32

5.85 × 10−3 1.12 1.19 1.12

6.11 3.45 2.78 2.58

− 7.64 5.81 4.45

− 55.0 52.2 42.1

10.41 7.26 5.71 4.56

− 7.59 5.53 4.34

3.94 × 0.295 4.69 × 3.50 × 1.55 ×

10−2 10−2 10−5 10−7

a

Provided for comparison are experimental log KO/W values, with the associated COSMO-SAC percent error, and log KO/W values from dedicated log KO/W algorithms. bCOSMO-SAC activity coefficients and log KO/W calculated by σ profiles produced with Gaussian 03.

point is sampled over a 1 ns MD run following an initial 1 ns equilibration period, which were performed using the DL_POLY56 program.

numerical instabilities. The thermodynamic analyses tend to be very sensitive to changes in λ1 at soft potentials; therefore, for λ1 < 0.50, an increment of 0.05 is used, whereas for λ1 ≥ 0.50, an increment of 0.10 is used. The solvent−solute electrostatic contribution is calculated, with full vdW interactions, using a simple variation of the point−point Coulombic interaction: qiqj Uijelec(r , λ2) = λ22 (22) r

3. RESULTS AND DISCUSSION 3.1. Application of the COSMO-SAC Method to the Computation of log KO/W for Selected Nanoscale Building Blocks. Table 1 compares the log KO/W calculated by COSMO-SAC for various compounds to values from experiment 57 and those from the dedicated log K O/W algorithms, ALOPs and AC_logP.7 When available, σ-profile data obtained from the VT-2005 σ-profile database26 were used to calculate γU∞,water and γU∞,octanol. For particles for which the σprofile was not available in the database, σ-profiles were generated using the COSMO capability in the Gaussian 03 software.58 The solute structures were relaxed using the NWCHEM code59 at the DFT level of theory with a 6-31G basis set and the B3LYP exchange/correlation functional. The COSMO calculation was performed after a point energy calculation using the B3PW91 exchange/correlation functional with a 6-31G(d, p) basis set, which is considered to be comparable60 to the basis set used by Mullins et al.26 The COSMO-SAC calculation was carried out using a modified version of source code provided by Oldland.61 Table 1 shows that, although COSMO-SAC is generally not as accurate as

where λ2 is the electronic softening parameter. Because ⟨(dUelec)/(dλ)⟩λ is generally less sensitive relative to λ2, λ2 is sampled in increments of 0.10 over 0.10 ≤ λ2 ≤ 1.0. If one applies eqs 19, 21, and 22, the ΔGhyd can be solved as follows: ΔG hyd =

dU vdW dλ1

1

∫0 +

1

∫0

dU elec dλ2

dλ1 λ1, λ 2 = 0

dλ2 λ1= 1, λ 2

(23)

where λ1 is the vdW interaction softening parameter, λ2 is the electrostatic softening parameter, and ⟨...⟩ denotes the time average of the summation of all solute−solvent interactions. For both the vdW and the electrostatic analyses, each data 4560

dx.doi.org/10.1021/ie202107t | Ind. Eng. Chem. Res. 2012, 51, 4556−4566

Industrial & Engineering Chemistry Research

Article

dedicated log KO/W algorithms, it performs favorably for the calculation of log KO/W over a wide range of solutes. For most components, COSMO-SAC predicts KO/W within an order of magnitude of experimental or predicted values. This is a very powerful observation as, unlike the log KO/W prediction algorithms, COSMO-SAC requires no phase equilibria data specific to the system of interest, making it more widely applicable. It can be seen in Table 1 that the values of log KO/W for polyaromatic hydrocarbons (PAHs) as predicted by COSMOSAC are considerably lower than that of experiment and dedicated log KO/W algorithms. This underestimation is likely due to the well-documented inability of COSMO-based methods to quantitatively predict γU∞,water for very hydrophobic solutes. The COSMO-SAC calculation of log KO/W for PAHs does qualitatively agree with values from experiment and dedicated log KO/W algorithms, as COSMO-SAC predicts the correct order of relative solubility for this family of compounds. Another limitation is that the values of γU∞ for aromatic alcohols calculated by COSMO-SAC appear artificially low. Fortunately, this seems to not affect the calculation of log KO/W, as the calculated values are comparable to those obtained from experiment and the dedicated algorithms. The activity coefficients of the aromatic alcohols are uncharacteristically small for both solvents, and any error in the calculation of the activity coefficients seems to cancel when the ratio of activity coefficients is taken in eq 7 to obtain log KO/W. Therefore, the error associated with the calculation of γU∞ for aromatic alcohols is likely due to a limitation common to both systems, such as the inability to accurately calculate the energy of the pure solute phase. Using the same methodology as applied on the compounds in Table 1, γU∞ values of C60 in water and octanol-rich environments are calculated. The σ-profile for C60 is not available in literature, and, therefore, the σ-profile was produced from a COSMO calculation, using atomic radii optimized for COSMO-RS,23 after a point energy DFT calculation. When performing a COSMO calculation on cage-like particles, such as POSS or C60, it is imperative to ensure that the inner cage areas are not being sampled. Given the size of the POSS cage, COSMO atomic radii, and COSMO probe radius, it is not possible to sample the interior of the POSS family of molecules. However, the same is not true for C60 and C60 derived compounds. To ensure that the Gaussian 03 COSMO calculation only considers the solvent-accessible area of C60, the COSMO area of C60 was compared to a Connolly surface area calculation62 with an excluded interior area and atomic radii consistent with that of COSMO standards. The COSMO calculated area was found to be 406.42 Å2, which is comparable to the Connolly area of 405.04. This comparison suggests that Gaussian 03 indeed excludes the interior area of cage-like molecules in COSMO calculations. The DFT calculation was performed using the 6-31G(d,p) basis set with the B3PW91 exchange-correlation functional. The resulting σ-profile can be found in Figure 2. The COSMO-SAC calculation yields activity coefficients of γC∞,water = 8.26 × 108 and γC∞,octanol = 19.0, 60 60 resulting in log KO/W = 6.79. This value is significantly lower than theoretical estimation of log KO/W = 17.6,38 yet higher than the experimental value of log KO/W = 6.67 found by Jafvert et al.63 The comparison of the COSMO-SAC log KO/W of C60 to the experimental value should be examined cautiously, as the experimental methods of that study may bias the water

Figure 2. Area weighted σ-profile for C60.

solubility measurements toward a more hydrophilic result. In particular, the C60 samples are dissolved in toluene, which is subsequently evaporated off before dissolution in water. However, it is unlikely that trace amounts of toluene were removed from the surface of C60 during evaporation. In turn, the toluene on the surface of C60 may behave as a weak surfactant, enhancing water solubility. When comparing the COSMO-SAC log KO/W of C60 against the log KO/W estimate from theoretical calculations, the tendency of COSMO-SAC to underestimate γU∞,water of hydrophobic solutes, shown in Table 1 and in multiple studies,26−28 must be considered. The underestimation of γU∞,water for C60 will produce an artificially low log KO/W, and, as a result, the estimation of log KO/W for C60 using COSMO-SAC should be considered a lower bound. Nevertheless, even given the possible underestimation of , COSMO-SAC predicts C60 to be exceptionally γC∞,water 60 hydrophobic, which is a conclusion reached in many studies.64−66 In addition to C60, COSMO-SAC calculations were also performed on the C60(OH)32, H-POSS, F-POSS, and OHPOSS nanoparticles. Figure 3 presents the resulting σ-profiles and Table 2 the results of these calculations. It is shown that COSMO-SAC predicts the hydrophilic shifts associated with the addition of hydroxyl groups to the C60 surface found in other studies.67−69 Like C60, H-POSS is shown to be very hydrophobic and is likely to preferentially partition into the organic phase. The addition of fluorine atoms to the POSS cage results in a modest hydrophilic shift; however, this particle is also predicted to partition heavily into the octanol phase of an octanol−water system. In contrast to this, the fully hydroxylated POSS molecule is predicted to have an extreme hydrophilic shift toward water solubility relative to the HPOSS molecule and is shown to be roughly neutral with regard to octanol−water partitioning. The C60(OH)32 molecule has been shown to have several stable isomers,49 and OH-POSS is likely to have low-energy rotations of the hydroxyl groups about the Si−O bond. It has been shown29 that COSMO-SAC may produce moderately different results for different configurations of the same molecule. When possible, a sensitivity analysis that investigates the influence of multiple configurations on COSMO-SAC estimations of γ should be performed. Given the size of these 4561

dx.doi.org/10.1021/ie202107t | Ind. Eng. Chem. Res. 2012, 51, 4556−4566

Industrial & Engineering Chemistry Research

Article

Figure 3. Area weighted profiles for (a) C60(OH)32, (b), H-POSS, (c) F-POSS, and (d) OH-POSS.

3.2. Solubility Response Associated with Adding Hydrophilic Groups Using Thermodynamic Integration and COSMO-SAC. The Gibbs free energies of hydration for C60(OH)32, H-POSS, F-POSS, and OH-POSS were calculated using the methodology described in section 2.2. As a test of the presented molecular-simulation-based methodology for computing the free energy of hydration, ΔGhyd is calculated for benzene in water using thermodynamic integration with the CHARMM vdW potential and electrostatics derived from a DFT/6-311G** Mulliken population analysis. Because the solubility of benzene in water is well-known, the accuracy of the calculation should indicate the efficacy of using Mulliken partial charges in a TI study. The resulting partial charges are shown in Table 3, where it is seen that the Mulliken partial charges are

Table 2. Infinite Dilution Activity Coefficients Calculated for Various C60 and POSS Functionalized Variants with Resulting log KO/W C60 C60(OH)32 H-POSS F-POSS OH-POSS

γU∞,water

γU∞,octanol

log KO/W

8.26 × 108 3.04 × 103 1.17 × 106 9.89 × 104 2.93 × 10−8

1.90 × 101 1.10 × 101 3.26 1.52 4.14 × 10−8

6.79 1.60 4.71 3.97 −0.995

solutes and the incomplete knowledge of their various isomer structures, the DFT geometry optimizations needed for a thorough sensitivity analysis are prohibitively computationally expensive and beyond the scope of this study. As such, it should be noted that a COSMO-SAC analysis on higher-energy configurations of C60(OH)32, and to a lesser degree OH-POSS, may produce different results for γU∞ and log KO/W. The more basic cage-like molecules studied in this work, such as C60 and H-POSS, are known to be quite rigid70−74 and are unlikely to exhibit other meaningful conformers. The OH-POSS study, like that of the hydroxylated aromatics in Table 1, indicates γU∞ values that are likely artificially low. Given the geometric and chemical similarity to the hydroxylated aromatic solutes, which exhibit the same behavior, this may be due to a poor estimation of the segment activity coefficient in the pure solute phase. However, as shown before, the inability to predict the energy of the pure solute reference state should not affect log KO/W calculations.

Table 3. Mulliken Partial Charges for Benzene CHARMM atom name

Mulliken (qi)

CHARMM22

GROMOS 43A1

GROMOS 53A6

CA HP

−0.15 0.15

−0.115 0.115

−0.10 0.10

−0.14 0.14

comparable to those of both the CHARMM22 potential and the GROMOS potential,75 the latter of which has been shown to produce accurate values of ΔGhyd of benzene in water.76 A 250 ps thermodynamic integration of the benzene molecule with Mulliken partial charges produces the vdW and electrostatic TI curves shown in Figure 4. Integrating Figure 4a over λ1 yields the first term in eq 23 resulting in hyd ΔGvdW = +1.29 kcal/mol, where integrating Figure 4b over λ2 4562

dx.doi.org/10.1021/ie202107t | Ind. Eng. Chem. Res. 2012, 51, 4556−4566

Industrial & Engineering Chemistry Research

Article

Figure 4. Thermodynamic integration curves for CHARMM benzene with Mulliken charges: (a) step 1, (b) step 2. hyd gives ΔGelec = −2.03 kcal/mol. Therefore, ΔGhyd = −0.74 kcal/ mol, which is in better agreement with the experimental value of ΔGhyd = −0.86 kcal/mol, than either GROMOS 43A1 (ΔGhyd = −1.15 kcal/mol) or GROMOS 53A6 (ΔGhyd = −1.60 kcal/mol).76 Applying the same technique to C60(OH)32 results in the partial charges shown in Table 4, shown with the partial charges

Table 5. Mulliken Partial Charges for H-POSS, F-POSS, and OH-POSS

Table 4. Comparison of Mulliken C60(OH)32 Partial Charges with TraPPE and CHARMM Partial Charges atom

Mulliken (qi)

TraPPE qi

CHARMM22a qi

CAd CB O H

0.09000 0.26125 −0.73000 0.39000

0.000 0.265 −0.700 0.434

0.00 0.23 −0.66 0.43

atom

H-POSS

F-POSS

OH-POSS

Si OP HP FP O H

1.1900 −0.6400 −0.2300 − − −

0.21005 −0.06670 − −0.11000 − −

0.57505 −0.27170 − − −0.44500 0.27750

functional groups of F-POSS and OH-POSS. The hydration energies from the resulting thermodynamic integration calculation are shown in Table 6. The COSMO-SAC Table 6. Solubility Parameters for POSS Molecules in H2O Calculated Using Thermodynamic Integration

a

The CAd and CB qi terms for the explicit-atom CHARMM22 force field are taken as the aggregate of the carbon charge and the charges of any bonded hydrogen atoms.

for similar atoms in TraPPE and CHARMM22 alcohol molecules, which are provided in an effort to compare the parametrization to an established potential. In Table 4, CB represents a carbon atom bonded to an oxygen atom, and CAd is a carbon atom bonded to a CB, but not to an oxygen atom. Despite the often-cited limitation of the Mulliken technique,77 and the fundamental differences between the hydroxyl group on a fullerenol and the hydroxyl group on a normal alcohol, the partial charges are quite comparable. The thermodynamic integration of C60(OH)32 in water yields hyd hyd ΔGvdW = +8.11 kcal/mol and ΔGelec = −53.00 kcal/mol, hyd = −44.89 kcal/mol. This result is in resulting in ΔG accordance with experimental studies suggesting that fullerenols with more than 12 hydroxyl groups are water-soluble.78,79 As expected, the thermodynamic integration of C60(OH)32 in water predicts a substantial hydrophilic shift relative to the base C60 molecule. The COSMO-SAC analysis of C60(OH)32 yields a log KO/W = 1.60, which implies a sharp increase in water solubility relative to C60. This response is at least a qualitative indicator that the COSMO-SAC method can distinguish between strongly hydrophobic and hydrophilic fullerenes. A similar thermodynamic integration analysis has been done for H-POSS, F-POSS, and OH-POSS. Table 5 summarizes the partial charges found from the Mulliken population analysis, with FP, OP, and HP representing the atoms that constitute the

solute

solv ΔGvdW (kcal/mol)

solv ΔGelec (kcal/mol)

ΔGsolv (kcal/mol)

H-POSS F-POSS OH-POSS

−6.78 −9.02 −7.51

−2.04 −1.96 −12.43

−8.82 −10.98 −19.94

calculation predicts a slightly lower log KO/W for F-POSS relative to the H-POSS particle, which is in agreement with the relatively subtle hydrophilic shift from H-POSS to F-POSS shown in the thermodynamic integration analysis. Furthermore, the COSMO-SAC calculation indicates a significantly higher degree of partitioning into the aqueous phase for OH-POSS relative to H-POSS and F-POSS. This calculation is also supported by the thermodynamic integration analysis, which shows that OH-POSS has a significantly lower ΔGhyd than does F-POSS or H-POSS.

4. CONCLUSION The ability of the COSMO-SAC method to accurately estimate the octanol/water partition coefficient for several variants of the C60 and POSS NBBs was analyzed. It was found that while the COSMO-SAC method may underestimate log KO/W for the very hydrophobic C60 particle, the predicted results still lie between the experimental and theoretical values reported in the literature. Additionally, thermodynamic integration studies were performed on C60(OH)32, H-POSS, F-POSS, and OHPOSS in water using parametrizations from DFT studies. The shift in relative solubility from the base NBB unit to the 4563

dx.doi.org/10.1021/ie202107t | Ind. Eng. Chem. Res. 2012, 51, 4556−4566

Industrial & Engineering Chemistry Research

Article

(12) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Application; Academic Press: San Diego, CA, 2002; Vol. 1. (13) Ben-Naim, A. Solvation Thermodynamics; Plenum Press: New York, 1987. (14) Smith, F. L.; Harvey, A. H. Avoid common pitfalls when using Henry’s Law. Chem. Eng. Prog. 2005, 33−39. (15) Grensemann, H.; Gmehling, J. Performance of a conductor-like screening model for real solvents model in comparison to classical group contribution methods. Ind. Eng. Chem. Res. 2005, 44, 1610− 1624. (16) Putnam, R.; Taylor, R.; Klamt, A.; Eckert, F.; Schiller, M. Prediction of infinite dilution activity coefficients using COSMO-RS. Ind. Eng. Chem. Res. 2003, 42, 3635−3641. (17) 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. (18) Hsieh, C. M.; Lin, S. T. Prediction of 1-octanol-water partition coefficient and infinite dilution activity coefficient in water from the PR plus COSMOSAC model. Fluid Phase Equilib. 2009, 285, 8−14. (19) Buggert, M.; Cadena, C.; Mokrushina, L.; Smirnova, I.; Maginn, E. J.; Arlt, W. COSMO-RS calculations of partition coefficients: Different tools for conformational search. Chem. Eng. Technol. 2009, 32, 977−986. (20) Buggert, M.; Mokrushina, L.; Smirnova, I.; Schomacker, R.; Arlt, W. Prediction of equilibrium partitioning of nonpolar organic solutes in water-surfactant systems by UNIFAC and COSMO-RS models. Chem. Eng. Technol. 2006, 29, 567−573. (21) Schroeder, B.; Santos, L.; Rocha, M. A. A.; Oliveira, M. B.; Marrucho, I. M.; Coutinho, J. A. P. Prediction of environmental parameters of polycyclic aromatic hydrocarbons with COSMO-RS. Chemosphere 2010, 79, 821−829. (22) Klamt, A.; Eckert, F.; Diedenhofen, M. Prediction or partition coefficients and activity coefficients of two branched compounds using COSMOtherm. Fluid Phase Equilib. 2009, 285, 15−18. (23) Klamt, A.; Jonas, V.; Burger, T.; Lohrenz, J. C. W. Refinement and parametrization of COSMO-RS. J. Phys. Chem. A 1998, 102, 5074−5085. (24) Eckert, F.; Klamt, A. Fast solvent screening via quantum chemistry: COSMO-RS approach. AIChE J. 2002, 48, 369−385. (25) Wang, S.; Sandler, S. I.; Chen, C. C. Refinement of COSMOSAC and the applications. Ind. Eng. Chem. Res. 2007, 46, 7275−7288. (26) 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. (27) 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. (28) Athes, V.; Paricaud, P.; Ellaite, M.; Souchon, I.; Furst, W. Vapour-liquid equilibria of aroma compounds in hydroalcoholic solutions: Measurements with a recirculation method and modelling with the NRTL and COSMO-SAC approaches. Fluid Phase Equilib. 2008, 265, 139−154. (29) 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. (30) Redmill, P. S.; Capps, S. L.; Cummings, P. T.; McCabe, C. A molecular dynamics study of the Gibbs free energy of solvation of fullerene particles in octanol and water. Carbon 2009, 47, 2865−2874. (31) Stavermann, A. J. The entropy of high polymer solutions. Recl. Trav. Chim. Pays-Bas 1950, 69, 163. (32) Guggenheim, E. A. Mixtures: The Theory of the Equilibrium Properties of Some Simple Classes of Mixtures, Solutions and Alloys; Clarendon Press: Oxford, UK, 1952. (33) Klamt, A.; Schuurmann, G. COSMO − A new approach to dielectric screening in solvents with explicit expressions for the

functionalized particle was tested with both the thermodynamic integration method and the COSMO-SAC method. The COSMO-SAC method is shown to be capable of predicting both subtle shifts in hydrophilicity, such as that seen upon the addition of fluorine atoms to the POSS cube, and extreme hydrophilic shifts, such as that seen from the hydroxylation of C60. These results demonstrate that COSMO-SAC provides a fast, simple, and inexpensive means of obtaining log KO/W for nanoparticle systems. This methodology is not meant to replace theoretical or experimental means of obtaining log KO/W values for NBB systems; rather, it is to be used as a tool to quickly estimate octanol/water partitioning for novel, functionalized NBBs. Given the rate at which new NBB variants are being synthesized, this should prove to be a powerful tool in terms of tuning the solubility of new particles.



AUTHOR INFORMATION

Corresponding Author

*Tel.: (615) 319-6477. Fax: (979) 690-3250. E-mail: psr@htri. net. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS I thank P. T. Cummings, C. McCabe, and H. Docherty for valuable comments provided in the preparation of this paper. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.



REFERENCES

(1) Roco, M. C. International perspective on government nanotechnology funding in 2005. J. Nanopart. Res. 2005, 7, 707−712. (2) Leo, A.; Hansch, C.; Elkins, D. Partiton coefficients and their uses. Chem. Rev. 1971, 71, 525−616. (3) Felder, R. M.; Rousseau, R. W. Elementary Principles of Chemical Processes; John Wiley & Sons: New York, 2000. (4) Fredenslund, A.; Jones, R. L.; Prausnitz, J. M. Group-contribution estimation of activity-coefficients in nonideal liquid-mixtures. AIChE J. 1975, 21, 1086−1099. (5) Hansch, C.; Leo, A. J. Substituent Constants for Correlation Analysis in Chemistry and Biology; John Wiley and Sons: New York, 1979. (6) Tetko, I. V.; Tanchuk, V. Y. Application of associative neural networks for prediction of lipophilicity in ALOGPS 2.1 program. J. Chem. Inf. Comput. Sci. 2002, 42, 1136−1145. (7) Tetko, I. V.; Tanchuk, V. Y.; Villa, A. E. P. Prediction of noctanol/water partition coefficients from PHYSPROP database using artificial neural networks and E-state indices. J. Chem. Inf. Comput. Sci. 2001, 41, 1407−1421. (8) Tetko, I. V.; Tanchuk, V. Y.; Kasheva, T. N.; Villa, A. E. P. Estimation of aqueous solubility of chemical compounds using E-state indices. J. Chem. Inf. Comput. Sci. 2001, 41, 1488−1493. (9) Breindl, A.; Beck, B.; Clark, T.; Glen, R. C. Prediction of the noctanol/water partition coefficient, logP, using a combination of semiempirical MO-calculations and a neural network. J. Mol. Model. 1997, 3, 142−155. (10) Taskinen, J.; Yliruusi, J. Prediction of physicochemical properties based on neural network modelling. Adv. Drug Delivery Rev. 2003, 55, 1163−1183. (11) Huuskonen, J. J.; Livingstone, D. J.; Tetko, I. V. Neural network modeling for estimation of partition coefficient based on atom-type electrotopological state indices. J. Chem. Inf. Comput. Sci. 2000, 40, 947−955. 4564

dx.doi.org/10.1021/ie202107t | Ind. Eng. Chem. Res. 2012, 51, 4556−4566

Industrial & Engineering Chemistry Research

Article

screening energy and its gradient. J. Chem. Soc., Perkin Trans. 2 1993, 799−805. (34) Klamt, A. Conductor-like screening model for real solvents − a new approach to the quantitative calculation of solvation phenomena. J. Phys. Chem. 1995, 99, 2224−2235. (35) Klamt, A.; Eckert, F. COSMO-RS: a novel and efficient method for the a priori prediction of thermophysical data of liquids. Fluid Phase Equilib. 2000, 172, 43−72. (36) Marcus, Y.; Smith, A. L.; Korobov, M. V.; Mirakyan, A. L.; Avramenko, N. V.; Stukalin, E. B. Solubility of C-60 fullerene. J. Phys. Chem. B 2001, 105, 2499−2506. (37) Hansen, C. M.; Smith, A. L. Using Hansen solubility parameters to correlate solubility of C-60 fullerene in organic solvents and in polymers. Carbon 2004, 42, 1591−1597. (38) Heymann, D. Solubility of fullerenes C60 and C70 in water. Lunar Planet. Inst. 1996, 27, 543−544. (39) Heymann, D. Solubility of C60 in alcohols and alkanes. Carbon 1995, 34, 627−631. (40) Merz, K. M.; Roux, B. Biological Membranes: A Molecular Perspective from Computation and Experiment; Birkhauser: Boston, 1996. (41) Martin, M. G. Comparison of the AMBER, CHARMM, COMPASS, GROMOS, OPLS, TraPPE and UFF force fields for prediction of vapor-liquid coexistence curves and liquid densities. Fluid Phase Equilib. 2006, 248, 50−55. (42) 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−935. (43) Cruz-Chu, E. R.; Aksimentiev, A.; Schulten, K. Water-silica force field for simulating nanodevices. J. Phys. Chem. B 2006, 110, 21497− 21508. (44) Chen, I. J.; Yin, D. X.; MacKerell, A. D. Combined ab initio/ empirical approach for optimization of Lennard-Jones parameters for polar-neutral compounds. J. Comput. Chem. 2002, 23, 199−213. (45) Jiang, J. W.; Sandler, S. I. Adsorption and phase transitions on nanoporous carbonaceous materials: insights from molecular simulations. Fluid Phase Equilib. 2005, 228, 189−195. (46) Bojan, M. J.; Steele, W. A. Interactions of diatomic-molecules with graphite. Langmuir 1987, 3, 1123−1127. (47) 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. (48) Wescott, J. T.; Fisher, L. R.; Hanna, S. Use of thermodynamic integration to calculate the hydration free energies of n-alkanes. J. Chem. Phys. 2002, 116, 2361−2369. (49) Rodriguez-Zavala, J. G.; Guirado-Lopez, R. A. Stability of highly OH-covered C-60 fullerenes: Role of coadsorbed O impurities and of the charge state of the cage in the formation of carbon-opened structures. J. Phys. Chem. A 2006, 110, 9459−9468. (50) Mulliken, R. S. Electronic population analysis on LCAO-MO molecular wave functions. J. Chem. Phys. 1955, 23, 1833, 1841−2338, 2343. (51) Wiberg, K. B.; Rablen, P. R. Comparison of atomic charges derived via different procedures. J. Comput. Chem. 1993, 14, 1504− 1518. (52) Breneman, C. M.; Wiberg, K. B. Determining atom-centered monopoles from molecular electrostatic potentials − The need for high sampling density in formamide conformational-analysis. J. Comput. Chem. 1990, 11, 361−373. (53) Besler, B. H.; Merz, K. M.; Kollman, P. A. Atomic charges derived from semiempirical methods. J. Comput. Chem. 1990, 11, 431− 439. (54) Martin, F.; Zipse, H. Charge distribution in the water molecule− A comparison of methods. J. Comput. Chem. 2004, 26, 97−105. (55) Beutler, T. C.; Mark, A. E.; Vanschaik, R. C.; Gerber, P. R.; Vangunsteren, W. F. Avoiding singularities and numerical instabilities in free-energy calculations based on molecular simulations. Chem. Phys. Lett. 1994, 222, 529−539.

(56) Smith, W.; Forrester, T. R. DL_POLY, 2.14; Daresbury Laboratory: Warrington, England, 2003. (57) Sangster, J. LOGKOW: A databank of evaluated octanol-water partition coefficients (LogP); http://www.vcclab.org/lab/alogps/start. html. (58) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scusseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vrevren, T.; Kudin, K. A.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Menucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakakima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Strattman, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pompelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenburg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Mallick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanyakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision C.01; Gaussian, Inc.: Wallingford, CT, 2004. (59) Kendall, R. A.; Apra, E.; Bernholdt, D. E.; Bylaska, E. J.; Dupuis, M.; Fann, G. I.; Harrison, R. J.; Ju, J. L.; Nichols, J. A.; Nieplocha, J.; Straatsma, T. P.; Windus, T. L.; Wong, A. T. High performance computational chemistry: An overview of NWChem a distributed parallel application. Comput. Phys. Commun. 2000, 128, 260−283. (60) Otto, A. H.; Steiger, T.; Schrader, S. HAlCl4 in the gas phase is stronger than HTaF6. Chem. Commun. 1998, 391−392. (61) Oldland, R. J. Predicting Phase Equilibria Using COSMO-based Thermodynamic Models and the VT-2004 Sigma-Profile. M.S. Thesis, Virginia Polytechnic Institute and State University, Blacksburg, VA, 2004. (62) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38. (63) Jafvert, C. T.; Kulkarni, P. P. Buckminsterfullerene’s (C-60) octanol-water partition coefficient (K-ow) and aqueous solubility. Environ. Sci. Technol. 2008, 42, 5945−5950. (64) Ruoff, R. S.; Tse, D. S.; Malhotra, R.; Lorents, D. C. Solubility of C-60 in a variety of solvents. J. Phys. Chem. 1993, 97, 3379−3383. (65) Walther, J. H.; Jaffe, R. L.; Kotsalis, E. M.; Werder, T.; Halicioglu, T.; Koumoutsakos, P. Hydrophobic hydration of C-60 and carbon nanotubes in water. Carbon 2004, 42, 1185−1194. (66) Terashima, M.; Nagao, S. Solubilization of 60 fullerene in water by aquatic humic substances. Chem. Lett. 2007, 36, 302−303. (67) Nakamura, E.; Isobe, H. Functionalized fullerenes in water. The first 10 years of their chemistry, biology, and nanoscience. Acc. Chem. Res. 2003, 36, 807−815. (68) Sayes, C. M.; Fortner, J. D.; Guo, W.; Lyon, D.; Boyd, A. M.; Ausman, K. D.; Tao, Y. J.; Sitharaman, B.; Wilson, L. J.; Hughes, J. B.; West, J. L.; Colvin, V. L. The differential cytotoxicity of water-soluble fullerenes. Nano Lett. 2004, 4, 1881−1887. (69) Da Ros, T.; Prato, M. Medicinal chemistry with fullerenes and fullerene derivatives. Chem. Commun. 1999, 663−669. (70) Chan, E. R.; Striolo, A.; McCabe, C.; Cummings, P. T.; Glotzer, S. C. Coarse-grained force field for simulating polymer-tethered silsesquioxane self-assembly in solution. J. Chem. Phys. 2007, 127. (71) Ionescu, T. C.; Qi, F.; McCabe, C.; Striolo, A.; Kieffer, J.; Cummings, P. T. Evaluation of force fields for molecular simulation of polyhedral oligomeric silsesquioxanes. J. Phys. Chem. B 2006, 110, 2502−2510. (72) Kroto, H. W. The stability of the fullerenes C-24, C-28, C-32, C36, C-50, C-60 and C-70. Nature 1987, 329, 529−531. (73) Li, H. C.; Lee, C. Y.; McCabe, C.; Striolo, A.; Neurock, M. Ab initio analysis of the structural properties of alkyl-substituted polyhedral oligomeric silsesquioxanes. J. Phys. Chem. A 2007, 111, 3577−3584. 4565

dx.doi.org/10.1021/ie202107t | Ind. Eng. Chem. Res. 2012, 51, 4556−4566

Industrial & Engineering Chemistry Research

Article

(74) Lu, X.; Chen, Z. F. Curved Pi-conjugation, aromaticity, and the related chemistry of small fullerenes (