Water Partition Coefficients of Alcohol Ethoxylate

Article pubs.acs.org/IECR

Predicting Octanol/Water Partition Coefficients of Alcohol Ethoxylate Surfactants Using a Combination of Molecular Dynamics and the Conductor-like Screening Model for Realistic Solvents Peyman Yamin, Rolf Isele-Holder, and Kai Leonhard* Chair of Technical Thermodynamics, RWTH Aachen University, 52062 Aachen, Germany S Supporting Information *

ABSTRACT: We present a hybrid modeling strategy for the prediction of octanol/water partition coefficients for alcohol ethoxylate surfactants of varying chain lengths. The strategy makes use of molecular dynamics (MD) simulations for the generation of molecular conformations in the presence of a solvent. A clustering of the conformations from the MD trajectories was carried out based on principal component analysis of their dihedral angles. Representative conformations thus selected were then used in the conductor-like screening model for realistic solvents (COSMO-RS). Each conformation has been assigned a weight using an equation derived on the basis of its probability of occurrence in the MD trajectory. Experimental partition coefficients were reproduced within conformation-independent accuracy of COSMO-RS despite the size and flexibility of the ethoxy chain otherwise posing a challenge on their solvation modeling.

■

INTRODUCTION As the entire chemical industry is gradually following the pharmaceutical branch in the application of molecular design methods for the discovery of new functional molecules, computational chemical strategies need to be devised for the wider class of molecules encountered. A large number of molecules dealt with in the chemical industries are small enough to be tractable using available computational chemical methods but have the characteristic of being partially made of flexible polar chains. This high flexibility is often a fundamental problem encountered in theoretical and computational study of their chemistry. Many surfactant classes include such mediumsize and flexible molecules. From a practical point of view, accurately predicting thermodynamic properties of surfactant systems plays an important role in modeling their micelle formation,1 estimating their environmental impact2 and toxicology.3,4 For instance, the extent of aggregation of micellar systems is highly sensitive to the activity coefficient of their monomeric surfactants. For a micelle with an aggregation number of 110, this extent almost triples with a 1% increase in the monomer’s activity coefficient (1.01110 = 2.99). Consequently, when designing a new surfactant with a desired critical micelle concentration in mind, an accurate knowledge of their monomeric activity coefficients is required. Another property of broad practical importance is the octanol/water partition coefficient of surfactant molecules. This property is conventionally expressed as the ratio of their infinite dilution activity coefficients in water and octanol and plays a crucial role in models dealing with their biological and environmental impact and toxicity. As an example, the chemical hazard assessment and risk management (CHARM) model5,6 can be © XXXX American Chemical Society

mentioned. CHARM is used for the assessment of surfactants’ environmental and toxicological behavior wherein various partitioning properties including their octanol/water partition coefficient play a significant role. In this work, we limit ourselves to the study of alcohol ethoxylate (AE) surfactants. AE surfactants are a major class of surfactants with a high annual production rate and diverse applications.7 Figure 1 shows the characteristic chain of carbon and etheric oxygen atoms for the AE surfactant C8EO5. The flexibility of this chain leads to complicated possibilities of interand intramolecular interactions such as associative bonding in the presence of solvent and renders the calculation of their liquid state conformations difficult. We present a hybrid approach for the computational thermodynamic modeling of AE surfactants particularly addressing the problem caused by their structural flexibility. We make use of the quantum mechanical conductor-like screening model for realistic solvents (COSMO-RS).8,9 A commonly encountered challenge in the application of solvation models in general is the acquirement of molecular conformations, often the only input required for their application. One remedy is a systematic conformational search. This involves discrete increments in the geometrical degrees of freedom of the molecule and would, in theory, scan the whole conformational space available to it. In practice, however, such an approach is only applicable to either very small molecules or Received: December 27, 2015 Revised: March 9, 2016 Accepted: March 28, 2016

A

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

Article

Industrial & Engineering Chemistry Research

Figure 1. AE surfactant C8EO5. AE surfactants are used in this work as model systems for the study of the influence of conformations on COSMORS calculations and are referred to using a general formula CxEOy, where x and y are the number of carbon atoms and ethoxy groups in the alkane and polyethylene glycol chains, respectively. All AE surfactants studied are linear and have one alcohol group at the beginning of their hydrophilic region.

■

METHODS Thermodynamic Modeling Using COSMO-RS. COSMO-RS 8,9 is an analytical model using statistical thermodynamics and is built on top of the COSMO solvation model. COSMO16 is a solvation model where the solute molecule is placed in a cavity surrounded by an ideal conductor continuum resembling a liquid solvent of very high polarity such as water. As a result, a boundary condition is created where the electrostatic potential on the surface of the cavity and beyond in the conductor vanishes, making possible remarkably simple calculations of the surface charge distribution and its interaction energy.16 COSMO calculations were done using TURBOMOLE (version 6.0.2)17 at a density-functional theory (DFT)18 level with the Becke-Perdew (BP)19,20 functional, and the triple-ζ valence plus polarization (TZVP)21,22 basis set was used to calculate the screening charge distribution on the surface of the cavity in which the molecule is placed. The distribution of the surface charge hence obtained is known as the σ-profile of the molecule. In COSMO-RS, solute and solvent molecules are both present and treated at the same level by means of their surface segment interactions. As a result, the analytical formulation of COSMO-RS renders possible fast calculation of the electrostatic component of solvation energy. For systems where localized interactions with the solvent such as hydrogen bonding play an important role, additional empirical terms are then added to take them into account. In COSMO-RS, then, the ensemble of interacting molecules in a liquid mixture is replaced by the ensemble of freely interacting surface segments. The calculation of molecular interactions in COSMO-RS is hence reduced to the calculation of the pairwise interactions of these surface segments with each other. Electrostatic and other contributions to the free energy of solvation such as association and dispersion are then taken into account to calculate the chemical potential of all the molecules in the mixture. This is done, however, based on the freesegment approximation. According to this approximation, all surface segments interact with each other without any geometrical constraints. One implication of this assumption on the calculation of solvation entropy relevant to EA surfactants will be discussed later. COSMO-RS calculations for the prediction of thermodynamic properties were all carried out using COSMOtherm (release C.2.1, COSMOlogic GmbH). In COSMO-RS calculations σ-profiles from preceding COSMO calculations are used to calculate the chemical potential of a solute in a given liquid mixture. To calculate partition coefficients, we have used the automatic calculation of activity coefficients in COSMOtherm rather than a direct calculation of chemical potentials and the subsequent calculation of activity coefficients from them. Thereby, COSMOtherm always calculates the reference state as a pure COSMO-RS liquid of the solute. Even though a separate calculation of chemical potentials by COSMOtherm would be more accurate, in many cases where the methodology can be of interest the pure solute has a large molecular size that makes the MD simulations of its pure form impractically expensive or is not a liquid at all.

special cases where only a small portion of it is subject to the search process. As the molecular size increases, however, so does the dimensionality of the conformational free energy surface. Consequently, a systematic search of molecular conformations becomes impractical for many substances of practical interest. Other examples of conformational search strategies include stochastic and heuristic search methods. Stochastic search methodologies rely on a random change in the molecular degrees of freedom, whereas a heuristic search is based on a set of rules for the generation of hypothetical conformations. Generally, conformations generated using any of these methods should then be evaluated using physicochemical or geometrical criteria resulting in a set of valid conformations. A review of common conformational search methods in use can be found in the literature.10 Addressing this conformational search problem in the context of solvation modeling using COSMO-RS, Buggert et al.11 have taken initial steps using explicit-solvent molecular dynamics (MD) simulations for the search of conformations. They have used short MD simulations (0.5 ns) from which conformations were extracted every 10 ps. The chosen conformations were not geometry-optimized afterward and were used directly for COSMO-RS calculations. Mokrushina et al.12 then probed other possibilities for combining MD simulations and COSMO-RS modeling using longer simulations. Conformations were extracted on the basis of their onedimensional free energy profiles expressed in terms of molecular radius of gyration, a mechanical property. Also in that work a direct use of conformations without further geometry optimization was employed. The present work can be seen as a continuation in which conformational free energy surfaces are calculated in terms of molecular dihedral angles, a purely geometrical property. Additionally, geometry optimizations were done for all conformations to the same level of quantum mechanical theory used in COSMO-RS. We further make use of MD simulations for conformation weighting prior to their use in COSMO-RS. In line with these previous works, we are also concerned with the modeling of the surfactants’ octanol/water partition coefficient, KOW defined as i K iOW =

ciO ciW

=

γi∞ ,W v W γi∞ ,O vO

(1)

where c and γ denote the molar concentration and infinite dilution activity coefficient of surfactant i, and v is the corresponding molar volume of the solvent. Similar attempts have been recently reported in the literature, where a combination of MD simulations and continuum solvation methods including COSMO-RS was used to model membrane−water partition coefficients13 and solvation free energies.14 Furthermore, and for the sake of comparison, results obtained using the heuristic search algorithm COSMOconf15 were also introduced and discussed. B


Article

Industrial & Engineering Chemistry Research Thermodynamic Treatment of Multiple Conformations. As mentioned before, molecular conformations are the only input required for COSMO and the subsequent COSMORS calculations. A thermodynamically correct use of multiple conformations of the same compound in COSMOtherm requires their individual weighting according to their free energy in the solution by a Boltzmann distribution and COSMOtherm can automatically carry out such a calculation of the conformational weights based on their COSMO-RS free energy in the mixture. Although these calculations are based on a simplified model of the liquid phase, they are very fast. Calculation of a molecule’s free energy by means of MD, on the other hand, requires multiple simulations, takes much longer and consumes cosiderably larger computational resources. Since in the methodology presented here the extraction of conformations of each molecule is done from a single MD simulation’s trajectory, the idea is to use the same trajectory to also calculate conformational weights. Hence, in the following we derive an equation showing how to use the weights from MD to weight the chemical potentials of the conformations calculated using COSMO-RS to obtain the chemical potential of the whole surfactant molecule. To derive this equation, we write the total free energy G for one phase of the mixture with components i each having conformations j as, G=

a conformation, while nk is a molecule, or as we have called it throughout the text, a component of the mixture. Two cases can then be distinguished in a mixture: when component i is the same as the component k, and when it is not. In the latter case, the differentiation is simply zero because no change in the number of moles of component i occurs. However, when i = k, the number of moles of component i is equal to the fraction of n nk that counts as a conformation ni,j of i, or i ,j . Hence, it is nk

shown that the chemical potential of a component, μk, can be written as a sum over the chemical potentials of its conformations weighted by their individual probability of occurrence, according to the following equation,

μk =

j

(2)

j

Further, the chemical potential of a component k can be denoted by μk and written as ⎛ ∂G ⎞ μk = ⎜ ⎟ ⎝ ∂nk ⎠T , P , n

(3)

k̅ ≠k

Substituting G from eq 2 into eq 3 and carrying out the partial derivative we then have, ⎛ ∂ni , j

∑ ∑ ⎜⎜

μk =

i

j

⎝ ∂nk

μi , j +

⎞ ni , j⎟⎟ ∂nk ⎠

∂μi , j

(4)

The total number of moles for component k is in turn given by the sum over all conformations j belonging to that component so that we can write,

nk =

∑ nk ,j (5)

j

From Gibbs−Duhem equation, however, we know that,

∑ ∑ ni ,j i

j

∂μi , j ∂nk

(8)

Conformational Search Using COSMOconf. Along with the methodology developed in this work, COSMOconf (version 2.1), was also used for finding relevant molecular conformation. COSMOconf applies a multistep search-andexclude strategy. We have used linear initial structures similar to the one shown in Figure 1. COSMOconf then continues to generate its starting conformations using a database of precalculated σ-profiles based on their similarity to the given initial structure by means of the COSMOfrag algorithm.23 The structures then go through multiple steps of subsequent geometry optimization and clustering. At each step a fraction of these clusters is eliminated and the rest enters the next step where the geometry optimization is done at a higher level. This process repeats until only a few conformations remain which are then taken to be suitable as the input for COSMO-RS calculations. Initial geometry optimization steps are carried out in the ideal gas and ideal conductor phases at an AM124 level. These are then followed by clustering based on structural and energetic similarities for conformations of the ideal gas phase and σ-match similarity (SMS)25 for those in the ideal conductor phase. SMS is equal to one for two identical σ-profiles and asymptotically approaches zero the more dissimilar they get. Toward the end of the algorithm, final energy calculations and geometry optimizations are again carried out in the gas and COSMO phases and the resulting conformations are then proposed for use in COSMO-RS. Generation and clustering of conformations in COSMOconf was done using the default values for all its adjustable parameters but SMSbond, the SMS value above which two structures are taken to represent a single conformation, which was set to 0.975. A sensitivity analysis of the final partition coefficients has been done with respect to the other adjustable COSMOconf parameters. It has been observed that although the number of conformations obtained could be increased, except for SMSbond, varying other parameters did not show any significant influence on the quality of the final partition coefficient predictions. MD Simulations for the Generation of Conformations. MD simulations were carried out using GROMACS26,27 program package and the GROMOS (G53A6) force field.28,29 TIP4P-Ew30 was used to model the water molecule whose rigidity throughout the simulations was guaranteed by means of the SHAKE31 algorithm. Detailed comparison with experimentally available physical properties for various water models can be found elsewhere.30 A time-step of 1 fs was used for all production simulations. Each MD simulation was equilibrated for 5 ns prior to a 100 ns final production run. Every 10 ps,

∑ ∑ nijμij i

∑ pj μk ,j

=0 (6)

and hence the second term in the right-hand side of eq 4 vanishes. On the other hand, regarding the first term we have the following condition, ⎧ ni , j for i = k ∂ni , j ⎪ = ⎨ nk ∑∑ ∂nk ⎪ i j ⎩ 0 for i ≠ k (7) In the context of the MD simulation trajectories corresponding to an ensemble with a constant number of total particles, ni , j pj = n corresponds to the probability of occurrence of the k

conformation j. Please note that the point of reference for ni,j is C


Article


measured partition coefficients are reported.45 This difference in the temperature is assumed to have a negligible effect on the ratio of the molecular volumes. Conformations Obtained Using COSMOconf. For each AE, i, Figure 2 shows a comparison between COSMO-RS predictions using COSMOconf conformations and the experimental values.

snapshots of the trajctory were taken for further analyses. All simulations used for production were carried out as NPT ensembles. The control of temperature and pressure was achieved using the Nosé−Hoover32,33 thermostat and the Parrinello−Rhaman34,35 barostat, respectively. Time constants of the thermostats and barostats were 1.5 and 2.5 ps, respectively. All simulations for the extraction of conformations were run at 23 °C and 1 bar. The smooth particle-mesh Ewald method was used for the electrostatic interactions.36,37 For AE molecules, the model generated by PRODRG38 with default force field parameters was used except for the proposed atomic charges. These were in turn obtained by applying the CHELPG39 algorithm on the molecular geometries optimized at the DFT level with the B3LYP functional and the TZVP basis set, except for C10EO8 where due to the convergence problems the optimization was carried out with the same basis set but instead at the Hartree−Fock (HF) level. Since GROMOS is a united atom force field, unpolar hydrogen charges were added to the corresponding carbon atoms. The starting conformations were linear and the geometry optimization was carried out in the vacuum. These calculations were carried out using the Gaussian 09 program.40 After each simulation the structures were again saturated with hydrogen atoms and DFT (B3LYP, BP-TZVP) level geometry optimizations were carried out with all the atoms taken into account. Force field parameters used for the octanol model correspond to those generated by PRODRG for all but the atomic charges and the harmonic potential parameters used for the description of the bond stretching. Respectively, these were taken from the charges proposed for alcohols in the protein databank implemented in GROMACS28 and parameters suggested in the literature41 for the octanol molecule. Principal component analysis of the dihedral angles, that is, the off-plane angles formed between each four atoms of a conformation, was then made according to a method proposed elsewhere.42 Hereby, all molecular dihedral angles, ϕn, were taken into account. To address the periodicity of the dihedral angles,42 an appropriate conversion of the system of coordinates was carried out and each dihedral angle was transformed into two new coordinates, q2n − 1 = cos(ϕn)

Figure 2. Octanol/water partition coefficients predicted using conformations obtained by COSMOconf. Experimental values are obtained from the literature.45

Experimental partition coefficients used45 are reported to have a maximum error of 0.09 log units, barely visible in the scale shown in Figure 2. The height of the error bars on each side for the calculated partition coefficients corresponds to the absolute value of the maximum error reported46 as slightly less than 0.8 log units. Since the molecules studied for obtaining this error estimate were both small and rigid, these error bars can be taken to correspond to the conformation-independent uncertainty in the calculation of partition coefficients. It can be seen that the use of conformations obtained using COSMOconf led to overestimations in the partition coefficient values and their deviation from the experimental measurements are larger than the conformation-independent uncertainty46 expected. It can be hence argued that this is due to the quality of the chosen conformations as opposed to that of the inherent COSMO-RS error. However, it is worth noting that these calculations could nevertheless reproduce the increase in the octanol/water partition coefficients as a function of their size. Extraction of Conformations from MD Trajectories. Trajectories obtained from simulating each AE separately in water and in water saturated octanol were used for principal component analysis of the converted set of coordinates. The first two greatest principal components were used to cluster available conformations. Here, the basic assumption is that each conformation is equivalent to the others belonging to the same cluster. If too rigorous a definition of a conformation, for example, a distinct set of three-dimensional coordinates was assumed, nearly all conformations in the trajectory would turn out to be unique. Hence, it is obvious that using a weaker definition is practically essential. We have therefore divided each principal component coordinate in equal segments and

(9)

and

q2n = sin(ϕn)

(10)

As discussed later, the free energy surface of each AE molecule was generated using the first two principal components as independent variables of the free energy. For each system, the number of conformations used for the COSMO-RS calculations is equal to the number of conformation clusters found in each free energy surface (FES).

■

RESULTS AND DISCUSSION Infinite dilution activity coefficients are calculated in pure water (W) and water-saturated octanol phase (O) with a water mole W fraction of 0.268.43 Approximate molar volume ratios v O were v taken from literature44 to be equal to 0.151. Both equilibrium concentrations and molar volume fractions are reported for 298.15 K while the temperature in all calculations is set equal to 296.15 K, the temperature for which the experimentally D


Article

Industrial & Engineering Chemistry Research used each two-dimensional interval to define a conformation. Having thus clustered conformations of the MD trajectories, probability ratio method47 was used to calculate the free energy of each cluster using the following equation: ΔG = −kBT (ln P − ln Pmax )

(11)

where ΔG is the free energy difference relative to the cluster with the minimum free energy. The conformation with the minimum free energy, that is, the one with the maximum probability of occurrence Pmax can then be assigned an arbitrarily chosen free energy value (Gmin = 0 in Figure 3). P

Figure 3. Free energy surface for C14EO8 in octanol at three different resolutions 25 × 25, 100 × 100, and 500 × 500, in order from left to right and top to bottom. A magnified region of the highest resolution graph is also shown at the bottom right. The position of the conformations obtained using COSMOconf is shown in the free energy surface at the bottom left.

Figure 4. Three conformations of C6EO5 corresponding to the lowest energy conformation of each method obtained using COSMOconf (1), MD simulations in water (2) and MD simulations in octanol (3). Aliphatic hydrogens are not shown. Carbon, oxygen, and the polar hydrogen atoms are shown in red, gray, and white, accordingly.

is then the probability with which a certain conformation, for which the relative free energy is to be calculated, appears in the trajectory of the simulated NPT ensemble. kB and T denote the Boltzmann constant and temperature, respectively. Figure 3 shows the free energy surface obtained for the simulation of C14EO8 in octanol at three different resolutions along with a magnification for a minimum’s region. Positions of the obtained COSMOconf conformations are shown in the pricipal component coordinates in the bottom left graph. It can be seen that they are almost in the same region but also that this region is far from the minima obtained from the MD simulations. It can be seen that they are at the center of the system of the two greatest principal components and hence do not span the conformational space of the molecule adequately. This was observed for the other AE molecules studied and with the other COSMOconf parameters used. For the case of C6EO5, Figure 4 shows three conformations obtained using COSMOconf and the MD simulations in octanol and water. Their corresponding σ-profiles are depicted in Figure 5. Figure 4 in which these conformations are depicted further shows that only a few dihedral angles are twisted as a result of which the conformations are more or less linear and appear as if they were rigid. For the systems studied then COSMOconf conformations are not a subset of the molecular dynamics conformations and vice versa.

Figure 5. σ-profiles of three C6EO5 conformations shown in Figure 4 and the octanol molecule. In each case, only the σ-profile of the lowest energy conformation is depicted.

The complete set of C6EO5 conformations is given in the Supporting Information (SI). Because of the tight positioning of a relatively large number of oxygen atoms in the molecules E


Article

Industrial & Engineering Chemistry Research studied, the electrostatics of the polyether chain is highly sensitive to the choice of conformations. As the area under the σ-profiles correlates with the molecular surface area, it can be seen from Figure 5 that the conformation obtained using COSMOconf has a larger surface due to its less folded structure. The σ-profiles corresponding to the conformations from the MD simulations and from COSMOconf in Figure 5 are qualitatively different in two major points indicated by the dotted circles in the graph. The left circle is a sign of a larger unpolar chain. The right circle indicates a peak at the midrange of the positive screening charge corresponding to negative atomic charges. This peak is additional to the other more positive one but is comparable to it in size. It hence indicates a larger separation of the oxygen atoms in the COSMOconf conformations compared to those extracted from the MD simulations, another sign of them being less folded. In the case of MD conformations and due to the twists in their dihedral angles the oxygen atoms are much more closely spaced and thus lead to a single concentrated negatively charged area. Difference in Conformational Weight and Free Energy in COSMO-RS and MD. For compounds with multiple conformations, each conformation should be assigned a weight with which it contributes to the overall chemical potential. This weight is calculated in COSMO-RS based on each conformation’s free energy in mixture. Comparing each conformation’s free energy obtained from the free energy surfaces and COMSO-RS, a discrepancy can be seen between the two methods. If there was no such discrepancy, the conformations obtained using MD simulations could be used in COSMO-RS without any further consideration. Their weights would have been thus calculated by COSMO-RS in a way corresponding to the weights obtained using the MD simulations and could be used to calculate the final chemical potential and activity coefficients with no further concern. There is, however, at least one known issue with COSMO-RS that makes a second thought about its calculation of conformational weights for the given systems worthwhile. Because of the free segment approximation mentioned earlier, in COSMO-RS, the probability with which a hydrogen bound water molecule can undergo a second bond with a neighboring oxygen atom is equal to the probability with which a free water molecule would undergo such a bond with that atom. According to the difference between the entropy loss for the second coupled hydrogen bond and the first single one, however, the probability of occurrence of a coupled hydrogen bond, given that the geometrical constraints of the conformational structure allow it, is higher than that of the first, single hydrogen bond. An unbound free water molecule undergoing a hydrogen bond loses all its translational and two out of three rotational degrees of freedom. The only degree of freedom left to such a molecule after the first hydrogen bond is the rotation around the newly formed hydrogen bond. If now the same molecule undergoes a second hydrogen bond, however, the loss in entropy is only due to the loss in one degree of freedom while the energy effect remains comparable to that of the first bond. This leads to the second hydrogen bonds being entropically more favorable than the initial ones. A schematic representation of such a coupled hydrogen bond is shown in Figure 6. Figure 7 shows the frequency of occurrence of single and double hydrogen bonds throughout 100 ns of the MD trajectory of C14EO8 in water captured by 10k snapshots. The GROMACS default criteria for the existence of a single hydrogen bond were taken corresponding to a donor−acceptor

Figure 6. Schematic representation of a double hydrogen bond between C6EO5 and a water molecule. The formation of the second bond is energetically similar to the first one but has a lower entropic cost.

Figure 7. Occurrence of single (black) and double (gray) hydrogen bonding constellations for C14EO8 in water. The single and double hydrogen bonds were counted independently and for each case (number and type) only the occurrence of that case is counted and shown. Hence, the height of the gray column at 0 indicates the number of snapshots with no incident of a double hydrogen bond while single hydrogen bonds might have occurred. Similarly, the height of the black column at 1 indicates the number of snapshots where only a single hydrogen bond was formed while double hydrogen bonds might have occurred as well.

distance of 0.35 nm and a hydrogen-donor−acceptor angle of 30°. The criterion for a double hydrogen bond is then that two single hydrogen bonds occur simultaneously and on the same water molecule in the MD trajectory snapshot. Figure 7 hence shows cases where a different number of single and double hydrogen bonds are formed including the number of times where no single or no double bonds could be identified. It can be seen that almost half the snapshots contain at least one double hydrogen bonding association. Further by looking at the frequencies of occurrence shown in Figure 7 corresponding to the formation of a single double and a single single hydrogen bond, it can be seen that the probability of having a single double association is more than double the probability of having two single ones. Since in COSMO-RS, however, all hydrogen bonds are treated as single ones, their probability will be always calculated as equal. From Figure 7, the number of times where some double hydrogen bonding occurs in the trajectory can be calculated as the sum of all gray columns’ heights except for that of the zeroth one. This in turn is equal to 10k minus the height of the zeroth gray column and is slightly more than 5k. Consequently, more than half of the F


Article


conformations are then weighted according to their frequency of occurrence in the MD trajectory and used for COSMO-RS modeling. All experimental points for a series of alcohol ethoxylates could be reproduced within the COSMO-RS error range, implying a prediction free of the error introduced by the choice of conformations. Despite the surprisingly good experimental agreement between the results of this hybrid modeling scheme for the studied ethoxy alcohols, a generalization of the application of this scheme should be done with care. Upon application of this modeling scheme on other chemical systems, the choice of MD force field should be done according to the chemical system used, so that the reliability of the obtained conformation distribution can be guaranteed beforehand.

snapshots contain some degree of double hydrogen bonding leading to erroneous calculation of their entropy and free energy. Obviously, geometrical constraints and the number of oxygen atoms limit the maximum number of double hydrogen bonds leading to a low number and a vanishingly small one for the probability of having three and four simultaneous double associations for C14EO8. Figure 8 shows the resulting prediction of the partition coefficients and their comparison with experimental values,

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.5b04955. Conformations obtained using the application of COSMOconf algorithm; those corresponding to the FES minima and their COSMOfiles are given in the three archive files COSMOconf-conformations.zip, MD-conformations.zip, and COSMOfiles-optimized-MD-conformations.zip, respectively. Details of the MD system size and the values of the octanol/water partition coefficient obtained using the methods discussed are given in Table S1 and S2. A series of tables (S3) is further provided containing the number of members for each FES cluster for all AE surfactants considered along with the calculated conformation weight for all the representative conformations of all clusters of the FES minima (ZIP)

Figure 8. Octanol/water partition coefficients calculated with conformations obtained using the MD-based strategy presented in this work using an automatic (gray) and manual (black) weighting procedure by COSMOtherm and eq 8, respectively.

using both COSMO-RS and MD-based weights calculated according to the formulation derived in the previous section. Regarding MD simulations, we know that if they are done for a long enough time to allow for reasonable sampling of the phase space an exact statistics for the given force field is obtained. Whether the computationally exact statistics thus obtained also corresponds to the physical system, however, on the applicability of the chosen force field on the system considered. In this work, we use the statistics obtained from the application of the general purpose force field GROMOS (G53A6) to obtain the probability of occurrence of various conformations. Although we cannot be sure of the quality of the selected force field for reproducing the correct statistics for alcohol ethoxylates, the final results are very promising. Nevertheless, it should be noted that a generalization of the application of this force field for other chemical classes should be done with care. The error bars shown in Figure 8 correspond to the expected COSMO-RS accuracy, provided the right choice of conformations and their weights is made. It can be seen that all predictions using the MD-based weighting procedure reproduce the experimental data within the mentioned conformation independent error range. The automatic weighting, however, does so only for two out of five surfactants.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Fax: +49 241 8092255. Tel: +49 241 8098174. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was performed as part of the Cluster of Excellence “Tailor-Made Fuels from Biomass”, which is funded by the Excellence Initiative by the German federal and state governments to promote science and research at German universities.

■

REFERENCES

(1) Clint, J. H. Surfactant Aggregation; Springer Science & Business Media: 2012. (2) Britton, L. N. J. Surfactants Deterg. 1998, 1, 109−117. (3) Maibach, H. I. Toxicology of Skin; CRC Press: 2001. (4) Adams, W. J.; Chapman, G. A. Aquatic Toxicology and Hazard Assessment: 10th Vol.; ASTM International: 100 Barr Harbor Drive, PO Box C700, West Conshohocken, PA, 1988. (5) Foekema, E. M.; Bos, A.; Verstappen, P.; Karman, C. C. Field validation of the CHARM-model; TNO report R98/317; The Netherlands Organization: 1998. (6) Schobben, H. P. M.; Vik, E. A.; Hutjes, G. G.; Karman, C. C.; Øfjord, G. D. Produced Water 2; Springer US: Boston, MA, 1996; pp 295−301.

■

CONCLUSIONS Prediction of octanol/water partition coefficients of selected alcohol ethoxylate surfactants was carried out using a combination of MD simulations and COSMO-RS. MD simulations were used to obtain the conformational free energy surface from which conformations are extracted that correspond to the minima of this free energy surface. These G


Article

Industrial & Engineering Chemistry Research (7) Cowan-Ellsberry, C.; Belanger, S.; Dorn, P.; Dyer, S.; McAvoy, D.; Sanderson, H.; Versteeg, D.; Ferrer, D.; Stanton, K. Crit. Rev. Environ. Sci. Technol. 2014, 44, 1893−1993. (8) Klamt, A.; Krooshof, G. J. P.; Taylor, R. AIChE J. 2002, 48, 2332−2349. (9) Klamt, A.; Eckert, F.; Arlt, W. Annu. Rev. Chem. Biomol. Eng. 2010, 1, 101−122. (10) Hatfield, M. P. D.; Lovas, S. Curr. Pharm. Des. 2014, 20, 3303− 3313. (11) Buggert, M.; Cadena, C.; Mokrushina, L.; Smirnova, I.; Maginn, E. J.; Arlt, W. Chem. Eng. Technol. 2009, 32, 977−986. (12) Mokrushina, L.; Yamin, P.; Sponsel, E.; Arlt, W. Fluid Phase Equilib. 2012, 334, 37−42. (13) Jakobtorweihen, S.; Ingram, T.; Smirnova, I. J. Comput. Chem. 2013, 34, 1332−1340. (14) Kolár,̌ M.; Fanfrlík, J.; Lepšík, M.; Forti, F.; Luque, F. J.; Hobza, P. J. Phys. Chem. B 2013, 117, 5950−5962. (15) Klamt, A.; Eckert, F.; Diedenhofen, M. J. Phys. Chem. B 2009, 113, 4508−4510. (16) Klamt, A.; Schuurmann, G. J. Chem. Soc., Perkin Trans. 2 1993, 2, 799−805. (17) Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Chem. Phys. Lett. 1989, 162, 165−169. (18) Koch, W.; Holthausen, M. C. A Chemist’S Guide to Density Functional Theory; Wiley-VCH Verlag GmbH: Weinheim, FRG: 2001. (19) Becke, A. D. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098− 3100. (20) Perdew, J. P. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 33, 8822−8824. (21) Schäfer, A.; Horn, H.; Ahlrichs, R. J. Chem. Phys. 1992, 97, 2571−2577. (22) Schäfer, A.; Huber, C.; Ahlrichs, R. J. Chem. Phys. 1994, 100, 5829−5835. (23) Hornig, M.; Klamt, A. J. Chem. Inf. Model. 2005, 45, 1169−1177. (24) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. J. Am. Chem. Soc. 1993, 115, 5348−5348. (25) Thormann, M.; Klamt, A.; Hornig, M.; Almstetter, M. J. Chem. Inf. Model. 2006, 46, 1040−1053. (26) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS: A message-passing parallel molecular dynamics implementation. Comput. Phys. Commun. 1995, 91, 43−56. (27) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. J. Chem. Theory Comput. 2008, 4, 435−447. (28) Oostenbrink, C.; Villa, A.; Mark, A. E.; van Gunsteren, W. F. J. Comput. Chem. 2004, 25, 1656−1676. (29) 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, 11205−1218. (30) Horn, H. W.; Swope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Head-Gordon, T. J. Chem. Phys. 2004, 120, 9665− 9678. (31) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. Comput. Phys. 1977, 23, 327−341. (32) Hoover, W. G. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (33) Nosé, S. Mol. Phys. 1984, 52, 255−268. (34) Nosé, S.; Klein, M. L. Mol. Phys. 1983, 50, 1055−1076. (35) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182−7190. (36) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577−8593. (37) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089−10092. (38) Schüttelkopf, A. W.; van Aalten, D. M. F. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2004, 60, 1355−1363. (39) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361−373.

(40) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision D.01; Gaussian, Inc.: Wallingford, CT, 2009 (41) Siwko, M. E.; de Vries, A. H.; Mark, A. E.; Kozubek, A.; Marrink, S. J. Biophys. J. 2009, 96, 3140−3153. (42) Mu, Y.; Nguyen, P. H.; Stock, G. Proteins: Struct., Funct., Genet. 2005, 58, 45−52. (43) Arce, A.; Blanco, A.; Souza, P.; Vidal, I. J. Chem. Eng. Data 1994, 39, 378−380. (44) Berti, P.; Cabani, S.; Mollica, V. Fluid Phase Equilib. 1987, 32, 195−203. (45) Wiggins, H.; Karcher, A.; Wilson, J. M.; Robb, I. Determining Octanol/Water Partitioning Coefficients for Industrial Surfactants; IPEC Conference, Albuquerque, NM, November 10−13, 2008; Vol. 13. (46) Klamt, A.; Jonas, V.; Bürger, T.; Lohrenz, J. C. W. Refinement and Parametrization of COSMO-RS. J. Phys. Chem. A 1998, 102, 5074−5085. (47) Knight, J. L.; Brooks, C. L. J. Comput. Chem. 2009, 30, 1692− 1700.

H


Water Partition Coefficients of Alcohol Ethoxylate

Recommend Documents