Definition and Computation of Intermolecular Contact in Liquids Using

Article pubs.acs.org/JPCA

Definition and Computation of Intermolecular Contact in Liquids Using Additively Weighted Voronoi Tessellation Rolf E. Isele-Holder,*,†,‡ Brooks D. Rabideau,† and Ahmed E. Ismail*,†,‡ †

Aachener Verfahrenstechnik: Molecular Sciences and Transformations, Faculty of Mechanical Engineering, and ‡AICES Graduate School, RWTH Aachen University, Schinkelstrasse 2, 52062 Aachen, Germany S Supporting Information *

ABSTRACT: We present a definition of intermolecular surface contact by applying weighted Voronoi tessellations to configurations of various organic liquids and water obtained from molecular dynamics simulations. This definition of surface contact is used to link the COSMO-RS model and molecular dynamics simulations. We demonstrate that additively weighted tessellation is the superior tessellation type to define intermolecular surface contact. Furthermore, we fit a set of weights for the elements C, H, O, N, F, and S for this tessellation type to obtain optimal agreement between the models. We use these radii to successfully predict contact statistics for compounds that were excluded from the fit and mixtures. The observed agreement between contact statistics from COSMORS and molecular dynamics simulations confirms the capability of the presented method to describe intermolecular contact. Furthermore, we observe that increasing polarity of the surfaces of the examined molecules leads to weaker agreement in the contact statistics. This is especially pronounced for pure water.

■

INTRODUCTION There is a strong need for models predicting thermodynamic properties of substances based on their molecular structure. A particularly powerful predictive model for organic compounds, with a strong physical background and a wide application range, is the conductor-like screening model for real solvents (COSMO-RS).1−3 COSMO-RS treats a liquid phase as an ensemble of densely packed molecules. As a simplification, it is assumed that the “surface” of a molecule can be divided into pieces, each of which has a direct contact partner. It is further assumed that all interactions in the liquid can be modeled as interactions between pairs of contacting surface pieces. To define interactions between these surface pieces, one must have well-defined properties describing the interaction potentials. This required information is provided by the COSMO4 model, a quantum-mechanical continuum solvation model. To avoid expensive sampling of phase space, another approximation is made within the COSMO-RS model: all surface segments are allowed to pair with one another on a completely random basis.5 By applying this approximation, the information on the relationships between neighboring segments in the original molecule is lost, and an inexpensive selfconsistent solution of the partition function for calculating the chemical potential of compounds in an arbitrary mixture can be determined, from which various thermodynamic properties of interest can be derived. Aside from the possibility of quickly performing thermodynamic property calculations, the free-segment approximation enables a simplified representation of molecules. Because interactions are determined by the surface charges and surface © 2012 American Chemical Society

segments can be paired completely randomly, the threedimensional structure of the surface (see Figure 10 below) of a molecule is not needed. Instead, so-called σ profiles, which denote the amount of surface area corresponding to a given surface charge σ, provide the information required to perform thermodynamic property calculations. Some examples of σ profiles are depicted in Figure 1. The strength of COSMO-RS is not the accuracy of its results, but rather its applicability to a wide range of compounds and its computational efficiency. However, COSMO-RS is also known to provide highly inaccurate predictions for some compounds, such as diethylamine and triethylamine,2 which we examine in the present study. The recently published COSMOtherm release6 made some major improvements for these compounds. However, the description of, for example, octanol−water partition coefficients of these compounds is still far below the mean expected accuracy of COSMO-RS (see Figure 3 below). The assumptions made within the COSMO-RS model have already been subjected to various tests for model molecules with simple shapes. These tests have usually been performed with regard to the model’s ability to determine thermodynamic properties.7−9 An alternative quantity to examine is the contact statistics of the paired surfaces produced by COSMO-RS. As all interactions are determined by contacting surface segments, the thermodynamic quantities are fully determined by the contact statistics. In fact, the contact statistics (or pair distribution functions) in COSMO-RS are the equivalent to radial Received: March 6, 2012 Published: April 5, 2012 4657

dx.doi.org/10.1021/jp3021886 | J. Phys. Chem. A 2012, 116, 4657−4666

The Journal of Physical Chemistry A

Article

As in many other areas,11 Voronoi tessellation has already been successfully applied to molecular systems. Kim and coworkers12,13 and Bajaj et al.14 used additively weighted and radical tessellation, respectively, as an auxiliary tool for the efficient construction of solvent-excluded surfaces. The application of unweighted tessellation for examining the structures of MD-simulated material has been studied, for example, by Ruocco et al.,15 Montoro and Abascal,16 Brostow,17 and Jedlovszky.18 Studies on the presence of voids were performed by Medvedev and co-workers using standard19 and later additively weighted20,21 tessellation. Voids were also studied by Debenedetti and co-workers22−25 using standard tessellation. Additively weighted tessellation for the examination of voids was also performed by Zhang et al.26 Fern et al.27,28 used unweighted tessellation as a tool for distinguishing whether a particle belongs to a condensed or vapor phase in two-phase MD simulations. David29 used standard tessellation to introduce surfaces between amino acids. Kim et al.30 used additively weighted Voronoi tessellation for creating surfaces in proteins. In this work, we use additively weighted Voronoi tessellation to construct intermolecular contact surfaces in MD simulations. These surfaces are used to determine contact statistics in MD, which we compare to contact statistics obtained from COSMORS for various compounds in the pure liquid state and in mixtures with water and octanol. In the next section, we provide an overview of the examined compounds, the calculation of contact statistics from COSMO-RS, the details of the MD simulations, and more details on the Voronoi tessellations used. An examination of the influence of the weights needed for Voronoi tessellation and an optimized set of weights are presented in the third section. This section also includes the results for the similarities between COSMO-RS and MD statistics. In the fourth section, we discuss the observed results before briefly summarizing our findings.

Figure 1. σ profiles of some of the studied compounds in the pure state.

distribution functions in molecular dynamics (MD) simulations, because both types of distribution functions specify the likelihood of all interactions in the system. In this study, we use contact statistics from COSMO-RS and compare them to surface contact statistics from MD simulations. Because surface contacts are not normally required in MD, we use the concept of Voronoi tessellation10 to introduce molecular surfaces in the simulation results. In Voronoi tessellation, each point in space is assigned to the object to which it is closest. For spherical particles, such as atoms, different types of Voronoi tessellation exist; these differ in the definition of “closest”. The most commonly used tessellations are unweighted, or standard, tessellation; radical tessellation; additively weighted tessellation; and multiplicatively weighted tessellation. Figure 2 shows the different

■

METHODS

Examined Compounds. The contact statistics of various compounds, listed in Table 1, were determined in this study. The chosen substances cover various functional groups frequently encountered in chemical engineering, such as water, alkanes, ethers, and esters. Furthermore, COSMO-RS is known to provide proper results for all of the compounds except diethylamine and triethylamine, as can be seen in Figure 3, which shows octanol−water partition coefficients obtained from experiment and from COSMO-RS calculations. Experimental partition coefficients for compounds containing S or F were taken from ref 31, whereas all other values were taken from ref 32. Additionally, all of the examined compounds are rather small and therefore easy to study in MD simulations. The simulated and experimental densities presented in Table 1 demonstrate the feasibility of simulating the compounds in MD. Details of the COSMO-RS calculations and MD simulations are described below. Contact Statistics from COSMO-RS. The calculation of contact probabilities in arbitrary mixtures is already included in COSMOtherm, the implementation of COSMO-RS. We used COSMOtherm, version 2.1,55 and the TZVP parametrization, which is recommended for obtaining accurate results. The probability, pab, that a segment a of molecule A is in contact with a segment b of molecule B is

Figure 2. Different tessellation types applied to the same set of spheres in two-dimensional space: (upper left) unweighted tessellation, (upper right) radical tessellation, (lower left) additively weighted tessellation, and (lower right) multiplicatively weighted tessellation.

tessellation types in two-dimensional space applied to the same set of spheres. The gray disks represent the atoms, which are the objects around which the Voronoi cells are constructed; the black lines are the bounds of the cells. Depending on the chosen tessellation type, the cells differ in size and shape. In some cases, the connectivity between neighboring atoms can be different as a result of the choice of tessellation. 4658



Article

The activity coefficients, γa and γb, of segments a and b in eq 1 are defined as

Table 1. Overview of the Examined Compounds and Their Experimental and Simulated Densities compound

formula H2O

water

experimental density (g/L)

C6H14 C6H5CH3 C2H5OH (C2H5)2O C4H8O CH3COCH3 C2H5OCOCH3 C4H4O CH3CN C4H9NH2 C5H5N OCHNH2 CH3NO2 C2H2OCHN (C2H5)2S CH3S2CH3 C2H5SH C6F14 CF3CH2OH (C2H5)2NH (C2H5)3N C8H17OH

ln γi(σ ) =

997.433

SPC/E TIP3P TIP3P-C TIP3P-Ew hexane toluene ethanol diethyl ether tetrahydrofuran (THF) acetone ethyl acetate furan acetonitrile n-butylamine pyridine formamide nitromethane oxazole diethyl sulfide dimethyl disulfide ethanethiol perfluorohexane trifluoroethanol diethylamine triethylamine octanol

simulated densitya (g/L)

655.234 862.235 785.136 707.837 883.338

998.4 985.7 1015.7 999.7 653.7 866.1 797.2 713.7 856.0

784.639 894.040 931.641 776.742 732.343 979.044 1129.145 1133.446 1080.047 831.148 1055.848 833.249 1676.150 1382.451 699.052 723.153 821.554

799.8 913.1 946.4 751.7 755.3 971.3 1122.9 1101.1 1100.3 808.5 1020.1 859.8 1672.5 1383.0 715.3 715.3 823.9

∑∑ i ∈ SA j ∈ S N

Figure 3. Octanol−water partition coefficients obtained from experiment and COSMO-RS calculations. The triangles represent partition coefficients that were calculated with the newly released BPTZVPD-FINE HB2012 parametrization.

A A total A total

⎛ E ⎞ exp⎜ − ab ⎟ ⎝ kT ⎠

pij = 1

(3)

where SA are all surface segments of molecule A and SN are all surface segments of all molecules in the mixture. Note that the contact probability pab is not necessarily equal to pba. This becomes obvious when studying an infinitely diluted solute in a pure solvent: the probability that a solvent molecule is in contact with the solute is zero, whereas the probability that the solute molecule is in contact with the solvent is unity. COSMOtherm provides segment−segment contact statistics. Because we do not expect molecular dynamics simulations to return the same set of surfaces as a COSMOtherm calculation, we summed over the contact probabilities for individual segments corresponding to a given atom to obtain atomspecific contact probabilities. Furthermore, it is not always possible to distinguish among all atoms of a given type in a COSMO-RS calculation. For example, the six hydrogen atoms at the ends of an n-hexane molecule are chemically identical and therefore their contact probabilities have to be considered together. This can be done, however, by simply summing the segment-specific contact probabilities. If conformers are present, they need to be accounted for separately, with all segments for each conformer treated individually. When combining the different conformers, the contact probabilities have to be weighted by the surface area and the occurrence probability of the conformer to which they belong. Details of the MD Simulations. All MD simulations were performed with LAMMPS.56 All molecules except water were modeled using the OPLS-AA force field.57−63 The contact statistics of different popular water models were examined: SPC/E,64 TIP3P,65 TIP3P as implemented in CHARMM66 (hereafter TIP3P-C), and the TIP3P-Ew67 and TIP4P-Ew68 extensions for use with long-range solvers. In all water simulations, the bond distance and bond angle of water were kept rigid using the SHAKE69 algorithm. Simulations were executed for all of the compounds from Table 1 except octanol in the pure state and for all of the compounds from Table 1 except the fluorine and sulfur compounds at infinite dilution in mutually saturated water and octanol. The TIP3P-Ew67 model was used in the mixtures. The numbers of molecules (nmolecules) used in the simulations for pure compounds are listed in Table 2. For simulations at infinite dilution, the water-rich phase contained 4123 water molecules and one solute molecule, whereras the octanol-rich phase contained 470 octanol molecules, 173 water molecules, and one solute molecule. The initial size of all simulation cells was set to 50 Å in each box dimension. All initial configurations of the simulations were created using Packmol.70 For all simulations, to avoid crashing, a soft potential was used to remove overlaps between nearby atoms. Afterward, the Lennard-Jones (LJ) and Coulomb interactions with a cutoff of 10 Å were switched on, including long-range LJ corrections71

Statistical uncertainties in simulated densities are less than 0.3 g/L for all species.

x BA a Abγaγb

(2)

where μS(σ) is the σ potential of the given mixture. The contact probabilities calculated with eq 1 are weighted by the area of the patches. Furthermore, the contact probabilities are normalized such that

a

pab =

μ S (σ ) RT

(1)

where xB is the molar fraction of species B; Aa and Ab are the surface areas of segments a and b, respectively; AAtotal is the total area of molecule A; Atotal is the surface area of all molecules in the mixture; k is the Boltzmann constant; T is the temperature; and Eab is the interaction energy between segments a and b. 4659



Article

Table 2. Examined Molecules, Numbers of Molecules Used in the Simulations of the Pure Liquids, Sets of Radii Found in the Scans for the Radical, Additively Weighted Tessellations and Corresponding rms Deviations, and Radii from the Overall Fit and Corresponding rms Deviations molecule

nmolecules ΔradrCH

ΔradrOH

ΔradrNH

σrms

ΔaddrCH

ΔaddrOH

ΔaddrNH

σrms

σrms

acont

(Å)

(Å)

(Å)

(rad, scan)

(Å)

(Å)

(Å)

(add, scan)

(add, opt)

(Å)

− − − − − 1.75 1.7 1.7 1.7 1.7a 1.7 1.7 1.7 1.7 1.7a 1.7 1.7a 1.7a 1.7a − − − − − − − −

1.15 − − − − − − 1.2 1.5 1.3 1.4 1.4 1.2 − − − 1.3 0.9 1.3 − − − − − − − −

− − − − − − − − − − − − − 1.4 1.0 1.4 1.2 1.0 1.4 − − − − − − − −

0.111 − − − − 0.034 0.052 0.046 0.039 0.028 0.047 0.036 0.027 0.048 0.035 0.027 0.061 0.049 0.035 − − − − − − − −

− − − − − 0.7 0.7 0.7 0.7 0.7a 0.7 0.7 0.7 0.7 0.7a 0.7 0.7a 0.7a 0.7a − − − − − − − − ΔaddrCH (Å) 0.71

0.35 0.35 0.35 0.35 0.35 − − 0.5 0.6 0.5 0.5 0.5 0.4 − − − 0.5 0.3 0.4 − − − − − − − − ΔaddrOH(Å) 0.47

− − − − − − − − − − − − − 0.4 0.4 0.5 0.4 0.4 0.5 − − − − − − − − ΔaddrNH(Å) 0.48

0.102 0.106 0.105 0.100 0.092 0.020 0.020 0.032 0.022 0.020 0.037 0.028 0.025 0.041 0.028 0.022 0.051 0.046 0.038 − − − − − − − − ΔaddrSH(Å) 0.93

− − − − − 0.021 0.013 0.033 0.023 0.018 0.038 0.031 0.026 0.044 0.031 0.020 0.052 0.053 0.038 0.013 0.014 0.021 0.0001 0.022 0.022 0.023 − ΔaddrFH(Å) 0.44

− − − − − 14.32 12.64 8.02 12.08 10.52 9.54 11.58 9.40 7.74 11.72 10.38 6.36 7.91 8.63 12.84 11.25 9.35 19.66 9.53 11.82 14.98 −

0.70

0.42

0.53

0.86

0.42

water SPC/E TIP3P TIP3P-C TIP3P-Ew TIP4P-Ew hexane toluene ethanol diethyl ether THF acetone ethyl acetate furan acetonitrile n-butylamine pyridine formamide nitromethane oxazole diethyl sulfide dimethyl disulfide ethanethiol perfluorohexane trifluoroethanol diethylamine triethylamine octanol

4123 4123 4123 4123 4123 575 704 1291 718 922 1020 848 1032 1426 755 931 1889 1397 1177 700 848 1019 370 1042 723 537 477

globally optimized radii COSMO-RS radii a

Value fixed before the scan.

and the PPPM72 solver for long-range electrostatics. The time step was set to 1 fs. After several short NVE simulations with limited movement of the atoms, Nosé−Hoover barostats and thermostats73−75 were switched on. The damping factors were set to τp = 1 ps for the pressure and τT = 100 fs for the temperature. The equations of motion were solved using the velocity Verlet algorithm.76 For systems of pure furan, pyridine, and oxazole, the initial temperature of the thermostat was set to 25 °C. The initial pressure was set to 51 atm. During a simulation of 1 ns, the pressure was slowly reduced to 1 atm. The use of the increased pressure was required to force the simulations into the liquid state. Afterward, the system was equilibrated for another 1 ns at a pressure of 1 atm, and then the simulation was run for 5 ns. During these runs, the density of the system was measured every 1000 time steps. Finally, the barostat was switched off, and simulations were run for another nanosecond. Snapshots of this part of the simulation were stored every 20 ps. These snapshots were used for determining the contact probabilities for the MD simulations. For all systems containing octanol, the initial temperature of the thermostat was set to 170 °C, and the initial pressure was set to 1 atm. After a simulation of 500000 time steps using

these conditions, the temperature was slowly reduced until a temperature of 25 °C was reached after another 1000000 time steps. The system was then equilibrated for another 500000 time steps. Simulations with pure octanol were run for a further 5 ns to determine the densities. Simulations with mixtures including octanol were switched to the NVT ensmble. The simulations were run for 2 ns, and snapshots stored every 0.5 ps were used to determine the contact statistics in the mixtures with octanol. For all other systems except simulations with TIP4P-Ew, the initial values of the thermostat and barostat were set to 25 °C and 1 atm, respectively. The simulations were equilibrated for 1 ns. In simulations with pure compounds, the simulations were run for additional 5 ns to determine the densities of the fluids. Afterward, simulations were performed in the NVT ensemble. Simulations of pure compounds were run for 1 ns. Snapshots for determining the contact statistics were stored every 20 ps. Simulations with solutes in water at infinite dilution were run in the NVT ensemble for another 2 ns. Snapshots for determining contact statistics were stored every 0.5 ps in this case. As the initial configuration for TIP4P-Ew, the final shapshot of the simulation for determining the density of the TIP3P-Ew model was chosen. The simulation was switched to the NVT 4660



Article

ensemble and run for 2 ns. In the last nanosecond of the simulation, snapshots for calculating contact statistics were stored every 20000 time steps. Contact Statistics from MD Simulations through Voronoi Tessellation. The differences of the four Voronoi tessellations mentioned in the Introduction are their definitions of the term closest. In standard or unweighted tessellation, the distance between a point in space and an object is defined as dstd = d =

(xc − x)2 + (yc − y)2 + (zc − z)2

σrms =

(4)

(5)

Δrad rij =

In additively weighted tessellation, distance is defined as78,79 dadd = d − radd

d rmult

i

j≥i

(8)

rrad, i 2 − rrad, j 2

(9)

whereas in additively weighted tessellation, the relevant quantity is

(6)

Δadd rij = radd, i − radd, j

The distance measure in multiplicatively weighted tessellation is80 dmult =

∑ ∑ (pij ,CRS − pij ,MD )2

between the COSMO-RS contact probabilities, pij,CRS, and the MD contact probabilities, pij,MD, was used, where n is the total number of contact types. Alternative approaches, such as using the relative deviation of the contact probabilities, were not used because small deviations in very small contact probabilities would lead to huge calculated error measures; matches between surfaces with small contact probabilities would dominate the fitting procedure in this case. As can be seen from eqs 5 and 6, the results of the tessellations depend on the radii assigned to the atoms. Strictly speaking, the results of the tessellation depend not on the absolute value of the radii but rather on their relation to each other. In radical tessellation, the relevant quantity is

where x, y, and z are the coordinates of the point in threedimensional Cartesian space and xc, yc, and zc are the coordinates of the center of the object. The distance in radical tessellation77 is defined as

drad 2 = d 2 − rrad 2

1 n

(10)

If, for a given tessellation, the atomic radii are adjusted in a way that conserves the corresponding Δr values for all pairs of atoms, the resulting tessellations, and thus the resulting contact statistics, will remain unchanged. As a consequence, the number of independent variables to tune is the number of unique atom types minus 1. In this study, the radii were assumed to be element-specific. In the calculations, the radius of the hydrogen atoms was set to an arbitrarily chosen small value. All other atomic radii were chosen and are reported relative to the hydrogen radius. We performed a series of calculations both to develop an understanding of the effect of atomic radius values and to determine an optimal set of radii. The step width of the scans was chosen as 0.05 Å for compounds containing two elements and 0.1 Å for compounds with three or more different elements. Because the fluctuations in the contact statistics for different snapshots were relatively small, as shown in the Supporting Information, only one snapshot was evaluated in the scans for each set of radii. The molecule-specific optimal radii and the corresponding rms deviations are reported in Table 2. Because of the overall better agreement of the contact statistics in the scans and the more reasonable values for the radii, further studies were performed using only additively weighted tessellation. An optimal set of radii for the CHON elements was determined such that the sum of the rms deviations for all compounds in Table 1 except water, octanol, diethylamine, and triethylamine was minimized. From the results obtained from the scan, a suitable starting point for an iteration was chosen. Contact statistics and rms deviations at the starting point and at six points in its surroundings were calculated for all compounds. For the points in the surroundings, one of the three independent variables was increased or decreased by 0.01 Å. A step with a width of 0.01 Å was performed in all directions with a negative gradient. Ten snapshots were used for calculating the contact statistics in the optimization process. The process was stopped when the rms deviation of the current point became smaller than the rms deviations of all of its neighbors. The points found in this way are not necessarily global minima. However, because of the shapes of the rms

(7)

In these equations, rrad, radd, and rmult are weights that are assigned to the objects. In radical and additively weighted tessellation, the radii of objects are commonly used as their weights. Even though atoms in molecular dynamics simulations do not have well-defined radii, this convention is also used in the rest of this article. Examining the results of the four different tessellations shown in Figure 2, we decided to focus on additively weighted and radical tessellations. A weakness of unweighted tessellations is that the surfaces of the Voronoi cells can pass through atoms (compare Figure 2), which is counterintuitive to the concept of an atomic radius. Multiplicatively weighted tessellation was excluded because the larger atoms dominated the smaller atoms too much in areas of low atom density. To create additively weighted tessellations, we used the algorithm and implementation provided by Medvedev et al.78 Radical tessellations were created using Voro++.81,82 Using the tessellations, we calculated the surface areas and derived the contact statistics from them. Thus, the statistics are dependent on the configurations obtained from the MD simulations and the radii assigned to the atoms. Kim et al.30 suggested using Bondi radii83 as atomic radii in additively weighted tessellation; however, as shown below, such an assumption does not represent an optimal solution in the present study. Instead, we fit the radii to the results for pure liquids in a way that contact statistics from MD simulation best reproduced COSMO-RS contact statistics. Using the optimized radii, contact statistics were computed for pure diethylamine and triethylamine, as well as for various compounds at infinite dilution in water and water-saturated octanol.

■

RESULTS Optimization of Radii. To quantify the similarity between surface contact statistics from COSMO-RS and MD simulations, the root-mean-square (rms) deviation 4661



Article

deviations as a function of the radii, which resemble broad valleys, we are confident that the calculated radii were sufficiently close to the global optimum. The radii for S and F were determined in a subsequent fit with already fixed values for the radii of C, H, O, and N. The optimal radii and the corresponding rms deviations are given in Table 2. In addition, the differences in the values of the COSMO-RS radii are included in this table. There is a good agreement between the reported values, with the largest deviation of 0.07 Å occurring for S. The contact probabilities obtained with the optimized radii are reported in the Supporting Information. Contact Statistics of Pure Compounds. Results of the scans using additively weighted tessellation for hexane are depicted in Figure 4. The strong influence of the chosen radius

Figure 5. COSMO-RS contact statistics and MD contact statistics of pure water as a function of ΔaddrOH. The different colors represent the different water models. From dark to bright: SPC/E, TIP3P, TIP3P-C, TIP3P-Ew, TIP4P-Ew.

from MD and COSMO-RS was observed. In contrast, increasing polarity and less weakly charged surfaces led to poorer agreement. This can be clearly seen by comparing the σ profiles of ethanol, formamide, and water in Figure 1 with the calculated rms deviations in Table 2. A closer examination of the contact statistics revealed that, for all compounds except n-butylamine, energetically favorable surface contacts were favored too much in COSMO-RS compared to MD, and energetically unfavorable contacts were avoided too much, as observed for ethanol, for example (see Figure 6). The energetically favorable contact between the

Figure 4. COSMO-RS contact statistics and MD contact statistics of pure hexane as a function of ΔaddrCH.

on the tessellation becomes obvious from the figure. For small values of ΔaddrCH, almost all intermolecular contact in MD simulations is between hydrogen atoms. This result is far from the contact statistics obtained using COSMO-RS. The limit ΔaddrCH = 0 Å corresponds to contact statistics obtained from the unweighted tessellation. Obviously, this tessellation type cannot be used for establishing the link between COSMO-RS and MD. The difference in the radii of the hydrogen and carbon atoms is 0.5 Å for the Bondi radii.83 For this value, the agreement between COSMO-RS and MD is better than for the unweighted tessellation; however, the differences are still large. For a difference in the radii of 0.7 Å, a reasonable agreement between COSMO-RS and MD was achieved. Results for the scans of pure water are depicted in Figure 5. A strong influence of the chosen radii can be detected for this compound (and also for all others). There is a slight difference among the different water models. However, the differences between the MD and COSMO-RS contact statistics are large for any chosen value of ΔaddrOH. The minimum rms deviation was obtained for ΔaddrOH = 0.35 Å and was almost equal to or higher than 0.1 for any of the examined water models, which is by far the highest minimum rms deviation for any of the pure compounds studied. These results show that a good agreement in contact statistics could be achieved for the nonpolar compound hexane, whereas the agreement was less satisfactory for the very polar water molecule; we found this to be true for other molecules as well. There is a correlation between the shape of the σ profiles and the rms deviation obtained from the calculations. For compounds that have only weakly charged surfaces or that have strongly charged surfaces with only one polarity, such as THF (see Figure 1), good agreement between the contact statistics

Figure 6. Contact statistics of pure ethanol from COSMO-RS and MD simulations.

oxygen and hydrogen atoms is preferred in COSMO-RS, whereas the unfavorable O−O contact is avoided. As the entire surface of the oxygen atoms must be in contact with some other atom, the aforementioned error in the contact probabilities induces an error in the C−O contact, which is energetically neither strongly advantageous nor disadvantageous. A similar behavior was found for each of the examined pure liquids except n-butlylamine. For n-butylamine, the opposite was observed. As can be seen from Figure 7, the energetically 4662



Article

One of the main approximations underlying the COSMO-RS model is that surface segments can be paired completely randomly. To examine this approximation, we determined conditional probabilities of surface contact of neighboring atoms. The conditional probability is defined as the probability that, given an arbitrary “fixed pair” of atoms already in contact, another pair of atoms of the same molecules (later called the conditional pair) in the vicinity of the first pair is in contact. We calculated these conditional probabilities up to the fifth nearest neighbor. Results for the calculations are depicted in Figure 9 for ethanol. Obviously, when the conditional pair is the “zeroth”

Figure 7. Contact statistics of pure n-butylamine from COSMO-RS and MD simulations.

favorable N−H contact occurred more often in MD than in COSMO-RS. A possible reason for this circumstance is presented in the Discussion section. Using the optimized radii, contact statistics of diethylamine and triethylamine were calculated from MD and compared to results from COSMO-RS. As can be seen in Table 2, the rms deviation between the MD and COSMO-RS contact probabilities is low, indicating good agreement. Molecular Contact Areas and Neighborhood Relations. One of the parameters introduced in COSMO-RS is the effective pairwise contact area of molecules, which is set to aeff = 6.15 Å2 in fully published COSMO-RS versions.3 This parameter is of critical importance for the results obtained from COSMO-RS calculations. We calculated the pairwise contact areas between molecules from MD simulations, denoted acont. The arithmetic means of these contact areas are listed in Table 2. The molecular contact areas calculated from MD are larger than those assumed in COSMO-RS. Except for formamide, the difference in the contact areas is greater than 1 Å2. In addition to the mean molecular contact areas, the distribution of the contact areas was determined. A sample distribution is shown in Figure 8 for n-butylamine. The figure has a high peak for very small contact areas and a second peak for intermediate contact areas. About 25% of the molecular contact areas are larger than 20 Å2. The shape of the depicted histogram is typical of all of the compounds examined here.

Figure 9. Conditional contact probabilities of ethanol.

nearest neighbor relative to the fixed pair for both molecules, the conditional pair is the fixed pair and therefore the conditional probability has to be unity. For more distant nearest neighbors, the conditional probability decays to small values. In the vicinity of the original pair, the conditional probabilities are significantly higher than zero. The range of values obtained from our analysis for the first three nearest neighbors are listed in Table 3. Because of their substantial Table 3. Ranges of Values for Conditional Contact Probabilities of All Compounds except Perfluorohexane nearest neighbor zeroth first second third

zeroth

first

second

third

1.00

0.49−0.60 0.28−0.45

0.23−0.34 0.18−0.28 0.13−0.26

0.09−0.24 0.11−0.20 0.11−0.20 0.09−0.32

differences from all other values, the results for perfluorohexane are reported separately in Table 4. Contact Statistics in Mixtures. The optimized radii were used to study compounds at infinite dilution in water and water-saturated octanol, which are the mixtures used in Table 4. Conditional Contact Probabilities of Perfluorohexane

Figure 8. Histogram of pairwise molecular contact areas of pure nbutylamine. 4663

nearest neighbor

zeroth

first

second

third

zeroth first second third

1.00

0.23 0.06

0.18 0.05 0.06

0.15 0.04 0.06 0.07



Article

models and the good agreement between simulated and experimental densities in Table 1 indicate that the configurations obtained from the simulations are reasonable. We have no proof that our introduced definition of surface contact in MD is correct; however, its ability to approach COSMO-RS contact probabilities suggests that it is an appropriate method for defining molecular surfaces and that our calculated radii are reasonable. The possible errors on the COSMO-RS side are related to the main approximations explained in the Introduction: the free-segment approximation and the surface-contact approximation. The close connection between independent pairing of surface patches and the preference for energetically favorable pairs, which was observed for all compounds in the pure state except n-butylamine, is obvious. For example, the COSMO-RS model allows for the OH hydrogen atom of one ethanol molecule to form a hydrogen bond with another ethanol molecule, without the two molecules’ oxygens or OH hydrogens interacting with each other. As a consequence, the energetically favorable O−H contact does not “notice” that a close approach of O−H also requires the approaches of O−O and H−H, which are unfavorable. As a consequence, the O−H contact is preferred too much in the COSMO-RS model. The results obtained for n-butylamine seem to contradict this explanation. However, examining the σ surface of n-butylamine depicted in Figure 10, one can see that the two negatively

COSMO-RS to calculate partition coefficients. The contact probabilities obtained from COSMO-RS and MD simulations are given in the Supporting Information. For the majority of the contact probabilities of all compounds, the agreement between COSMO-RS and MD is quite reasonable. Strong deviations can be found, for example, for the contact probabilities of the hydrogen atoms in ethanol, whose contact with the octanol oxygen is preferred strongly in MD compared to COSMO-RS; for the contact probabilities of the nitromethane oxygen in pure water; and for some of the formamide contact probabilities in either mixture.

■

DISCUSSION For almost all of the compounds studied, the agreement found in the COSMO-RS and MD statistics is good both in the pure state and in mixtures. This finding is not self-evident, because the number of independent contact probabilities is higher than the number of variables “tuned” in the tessellations. This is especially true for the uniform set of radii optimized for all compounds. The stability of the molecule-specific chosen radii and the corresponding rms deviations reported in Table 2 demonstrates that additively weighted tessellation is superior to radical tessellation in terms of matching contact probabilities from MD to those from COSMO-RS. Except for oxazole, the rms deviations obtained from additively weighted tessellation were smaller than those obtained from radical tessellation. For oxazole, radical tessellation worked only slightly better. Because of its inability to distinguish between different elements, unweighted tessellation was inappropriate for this study. As can be seen from Figures 4 and 5, the radii chosen for performing the tessellations have a strong influence on the resulting contact statistics. It is encouraging to see that a fixed set of radii was able to reproduce contact statistics from COSMO-RS quite accurately, especially for the contact statistics that were fully predictivethat is, those for pure diethylamine and triethylamine and for all mixtures. For the compounds used to fit the radii, the rms deviations obtained from the optimized set were usually only slightly higher than the values observed for the molecule-specific optimum set of radii obtained from the scanand in some cases were even smalleralthough the fluctuations in the molecule-specific optimal values for ΔaddrOH are large. Furthermore, the strong agreement between contact statistics of mixtures including water, whose COSMO-RS contact probabilities in the pure liquid could not be approached properly with results from MD for any value of ΔaddrOH, supports the use of the optimized set of radii. Two trends can be recognized form the results presented in the previous section. The first is that the stronger the surface charges in COSMO-RS, and thus the stronger the interactions, the worse the agreement with the contact statistics from MD. The second important trend is that, with one exception (nbutylamine), energetically favorable contact was preferred and energetically unfavorable contact was avoided too much in COSMO-RS compared to MD. Inaccuracies in the contact statistics derived from either model could be responsible for the observed trends. On the MD side, the contact statistics depend on the configurations provided from the simulation itself, the choice of weighted Voronoi tessellation to introduce surface contact, and the choice of weights assigned to the atoms. However, both the use of the popular OPLS-AA force field and the different water

Figure 10. σ surface of n-butylamine. Red, strongly positive surface charge densities; blue, strongly negative surface charge densities; green, neutral or nearly neutral surfaces.

charged surfaces of the amine hydrogens are closely spaced. Because of their proximity, these hydrogen atoms have an increased drive to approach nitrogen atoms from other molecules, which leads to a high N−H contact probability. The information that the amine hydrogens are closely spaced is lost in COSMO-RS, meaning that the increased drive cannot be detected in this model. The weakness of the free-segment approximation is also apparent from the results in Figure 8. We found the average contact area between molecules to be significantly higher in MD than what is assumed in COSMO-RS. Furthermore, the contact areas between molecules exhibit strong fluctuations. The conditional probabilities reported in the same section 4664


The Journal of Physical Chemistry A underline the strong dependence of neighboring atoms and therefore contradicts the free-segment approximation. As can be seen in Figure 3, the use of diethylamine and triethylamine in COSMO-RS is related to inaccuracies in calculated thermodynamic properties. Consequently, studying their contact probabilities is particularly interesting. For pure compounds, we found that the contact probabilities predicted from MD simulations using the optimized radii were in very good agreement with those from COSMO-RS. This agreement matches the discovered relation between the shape of the σ profiles and the agreement in the contact statistics, because neither diethylamine nor triethylamine has a strongly negatively charged surface. On the other hand, this result could not be foreseen, because the description of the partition coefficients is bad in COSMO-RS for these compounds. The inability of COSMO-RS to describe secondary and tertiary amines accurately might be related to the mixtures in which they are studied, which indicates that their contact statistics in the mixtures for which COSMO-RS is known to fail, such as mutually saturated water and octanol, should be different from those obtained from MD. In practice, however, the contact statistics from MD and COSMO-RS were very similar in the mixtures. Hence, the bad results for calculated partition coefficients in COSMO-RS are not related to incorrect pairing of surface segments.

■

REFERENCES

(1) Klamt, A. J. Phys. Chem. 1995, 99, 2224−2235. (2) Klamt, A.; Jonas, V.; Bürger, T.; Lohrenz, J. C. W. J. Phys. Chem. A 1998, 102, 5074−5085. (3) Klamt, A.; Eckert, F. Fluid Phase Equilib. 2000, 172, 43−72. (4) Klamt, A.; Schüürmann, G. J. Chem. Soc., Perkin Trans. 2 1993, 2, 799−805. (5) Klamt, A.; Eckert, F.; Wolfgang, A. Annu. Rev. Chem. Biomol. Eng. 2010, 1, 101−122. (6) Klamt, A.; Reinisch, J.; Eckert, F.; Hellweg, A.; Diedenhofen, M. Phys. Chem. Chem. Phys. 2012, 14, 955−963. (7) Klamt, A.; Krooshof, G. J. P.; Taylor, R. AIChE J. 2002, 2332− 2349. (8) Klamt, A.; Leonhard, K. Fluid Phase Equilib. 2007, 261, 162−167. (9) Neiman, M.; Cheng, H.; Parekh, V.; Peterson, B.; Klier, K. Phys. Chem. Chem. Phys. 2004, 6, 3474−3483. (10) Voronoi, G. J. Reine Angew. Math. 1907, 133, 97−178. (11) Okabe, A.; Boots, B.; Sugihara, K.; Chiu, D. S. N.; Chiu, S. N. Spatial Tessellations: Concepts and Application of Voronoi Diagrams, 2nd ed.; Wiley: Chichester, U.K., 2000. (12) Kim, D.-S.; Cho, Y.; Sugihara, K.; Ryu, J.; Kim, D. Comput.Aided Des. 2010, 42, 911−929. (13) Ryu, J.; Kim, D.; Cho, Y.; Park, R.; Kim, D.-S. Comput.-Aided Des. 2005, 2, 1−10. (14) Bajaj, C. L.; Pascucci, V.; Shamir, A.; Holt, R. J.; Netravali, A. N. Discrete Appl. Math. 2003, 127, 23−51. (15) Ruocco, G.; Sampoli, M.; Vallauri, R. J. Chem. Phys. 1992, 96, 6167−6176. (16) Montore, J. C. G.; Abascal, J. L. F. J. Phys. Chem. 1993, 97, 4211−4215. (17) Brostow, W. Phys. Rev. B 1998, 57, 13448−13458. (18) Jediovsyky, P. J. Chem. Phys. 2000, 113, 9113−9121. (19) Voloshin, V. P.; Beaufils, S.; Medvedev, N. N. J. Mol. Liq. 2002, 96−97, 101−112. (20) Alinchenko, M. G.; Anikeenko, A. V.; Medvedev, N. N.; Voloshin, V. P.; Mezei, M.; Jedlovzky, P. J. Phys. Chem. B 2004, 108, 19056−19067. (21) Anikeenko, A. V.; Alinchenko, M. G.; Voloshin, V. P.; Medvedev, N. N.; Gavrilova, M. L.; Jedlovszky, P. Comput.-Aided Des. 2004, 3045, 217−226. (22) Corti, D. S.; Debenedetti, P. G.; Sastry, S.; Stillinger, F. H. Phys. Rev. E 1997, 55, 5522−5534. (23) Sastry, S.; Corti, D. S.; Debenedetti, P. G.; Stillinger, F. H. Phys. Rev. E 1997, 56, 5524−5532. (24) Sastry, S.; Debenedetti, P. G.; Stillinger, F. H. Phys. Rev. E 1997, 56, 5533−5543. (25) Debenedetti, P. G.; Truskett, T. M. Fluid Phase Equilib. 1999, 158−160, 549−556. (26) Zhang, C.; Duan, Y.; Li, M. Geochim. Cosmochim. Acta 2010, 74, 4140−4149. (27) Fern, F. T.; Keffer, D. J.; Steele, W. V. J. Phys. Chem. B 2007, 111, 3469−3475. (28) Fern, F. T.; Keffer, D. J.; Steele, W. V. J. Phys. Chem. B 2007, 111, 13278−13286.

CONCLUSIONS Additively weighted and radical Voronoi tessellation were used to introduce intermolecular surfaces in MD-simulated liquids to compare contact statistics from MD and COSMO-RS. Our analysis showed that additively weighted tessellation is more suitable for this purpose than radical tessellation and much better than unweighted tessellation. We obtained a set of weights for the elements C, H, O, N, S, and F with which reasonable agreement in the contact statistics was observed for all compounds except pure water. Furthermore, we observed that increasing polarity in the surface charges in COSMO-RS leads to less agreement in the contact statistics derived from the different simulation methods. This becomes especially obvious for the study of pure water. Additionally, in almost all cases, energetically favorable surface contact is preferred and energetically unfavorable contact is avoided too much in COSMO-RS. This systematic deviation is likely to be a result of the approximations made within the COSMO-RS model. ASSOCIATED CONTENT

S Supporting Information *

Tables containing element-specific contact statistics from COSMO-RS and MD simulations with the optimized set of radii of all studied liquids. Statistical uncertainty for the results from MD simulations. This material is available free of charge via the Internet at http://pubs.acs.org.

■

ACKNOWLEDGMENTS

The authors thank Kai Leonhard, Peyman Yamin, and Andreas Klamt for helpful discussions and comments on this research. This work was performed as part of the Cluster of Excellence “Tailor-Made Fuels from Biomass” funded by the Excellence Initiative by the German federal and state governments to promote science and research at German universities. R.E.I.-H. was supported in part through a fellowship from the Aachen Institute for Advanced Study in Computational Engineering Science (AICES); financial support from the Deutsche Forschungsgemeinschaft (German Research Association) through Grant GSC 111 is gratefully acknowledged.

■

■

■

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] (R.E.I.-H.), [email protected]. edu (A.E.I.). Notes

The authors declare no competing financial interest. 4665



Article

(29) David, C. W. Biopolymers 1988, 27, 339−344. (30) Kim, C.-M.; Won, C.-I.; Cho, Y.; Kim, D.; Lee, S.; Bhak, D.-S.; an Kim, J. Comput.-Aided Des. 2006, 38, 1192−1204. (31) Howard, P. H.; Meylan, W. Handbook of Physical Properties of Organic Chemicals; CRC Press: Boca Raton, FL, 1996. (32) Sangster, J. J. Phys. Chem. Ref. Data 1989, 18, 1111−1227. (33) Pátek, J.; Hrubý, J.; Klomfar, J.; Součková, M.; Harvey, A. H. J. Phys. Chem. Ref. Data 2009, 38, 21−29. (34) Nayak, J. N.; Aralaguppi, M. I.; Aminabhavi, T. M. J. Chem. Eng. Data 2003, 48, 1152−1156. (35) Harris, K. R. J. Chem. Eng. Data 2000, 45, 893−897. (36) Djojoputro, H.; Ismadji, S. J. Chem. Eng. Data 2005, 50, 2003− 2007. (37) Nath, J. J. Chem. Thermodyn. 1996, 28, 167−170. (38) Toti, U. S.; Kariduraganavar, M. Y.; Aralaguppi, M. I.; Aminabhavi, T. M. J. Chem. Eng. Data 2000, 45, 920−925. (39) Estrada-Baltazar, A.; León-Rodriguez, A. d.; Hall, K. R.; RamosEstrada, M.; Iglesias-Silva, G. A. J. Chem. Eng. Data 2003, 48, 1425− 1431. (40) Gardas, R. L.; Johnson, I.; Vaz, D. M. D.; Fonseca, I. M. A.; Ferreira, A. G. M. J. Chem. Eng. Data 2007, 52, 737−751. (41) Naorem, H.; Suri, S. K. Can. J. Chem. 1989, 67, 1672−1675. (42) Roy, M. N.; Sarkar, B. K.; Chanda, R. J. Chem. Eng. Data 2007, 51, 1630−1637. (43) Domíniguez, M.; Rodriguez, S.; López, M. C.; Royo, F. M.; Urieta, J. S. J. Chem. Eng. Data 1996, 41, 37−42. (44) Nayak, J. N.; Aralaguppi, M. I.; Toti, U. S.; Aminabhavi, T. M. J. Chem. Eng. Data 2003, 48, 1483−1488. (45) Cases, A. M.; Marigliano, A. C. G.; Bonatti, C. M.; Sólimo, H. N. J. Chem. Eng. Data 2001, 46, 712−715. (46) Garcia-Miaja, G.; Troncoso, J.; Romani, L. J. Chem. Eng. Data 2007, 52, 2261−2265. (47) McCormick, D. G.; Hamilton, W. S. J. Chem. Thermodyn. 1978, 10, 275−278. (48) Marongiu, B.; Marras, G.; Porcedda, S. J. Chem. Eng. Data 1993, 38, 638−643. (49) Haines, W. E.; Helm, R. V.; Bailey, C. W.; Ball, J. S. J. Phys. Chem. 1954, 58, 270−278. (50) Dias, A. M. A.; Gonçalves, C. M. B.; Caço, A. I.; Santos, L. M. N. B. F.; Piñeiro, M. M.; Vega, L. F.; Coutinho, J. A. P.; Marrucho, I. M. J. Chem. Eng. Data 2005, 50, 1328−1333. (51) Atik, Z. J. Chem. Eng. Data 2007, 52, 1128−1130. (52) Resa, J. M.; Gonzáles, C.; Landaluce, S. O., de; Lanz, J. J. Chem. Eng. Data 2000, 45, 867−871. (53) Nayak, J. N.; Aralaguppi, M. I.; Aminabhavi, T. M. J. Chem. Eng. Data 2003, 48, 1152−1156. (54) Bhattacharjee, A.; Roy, M. N. J. Chem. Eng. Data 2010, 55, 5914−5920. (55) Eckert, F. COSMOtherm Users Manual; COSMOlogic GmbH & Co. KG: Leverkusen, Germany, 2009. (56) Plimpton, S. J. Comput. Chem. 1995, 117, 1−19. (57) Jorgensen, W. J.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225−11236. (58) Damm, W.; Frontera, A.; Tirado-Rives, J.; Jorgensen, W. L. J. Comput. Chem. 1997, 18, 1955−1970. (59) Rizzo, R. C.; Jorgensen, W. L. J. Am. Chem. Soc. 1999, 121, 4827−4836. (60) Price, M. L. P.; Ostrovsky, D.; Jorgensen, W. J. J. Comput. Chem. 2001, 22, 1340−1352. (61) Jorgensen, W. L.; McDonald, N. A. J. Mol. Struct. (THEOCHEM) 1998, 424, 145−155. (62) McDonald, N. A.; Jorgensen, W. L. J. Phys. Chem. B 1998, 102, 8040−8059. (63) Watkins, E. K.; Jorgensen, W. L. J. Phys. Chem. A 2001, 105, 4118−4125. (64) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269−6271. (65) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Phys. Chem. 1983, 79, 926−935.

(66) Durell, S. R.; Brooks, B. R.; Ben-Naim, A. J. Phys. Chem. 1994, 98, 2198−2202. (67) Price, J. D.; Brooks, C. L., III. J. Chem. Phys. 2004, 121, 10096− 10103. (68) Horn, H. W.; Swope, J. C.; Pitera, J. W.; Madera, J. D.; Dick, T. J.; Hura, G. L.; Head-Gordon, T. J. Chem. Phys. 2004, 120, 9665− 9678. (69) Ryckaert, J. P.; Berendsen, H. J. C.; Ciscotti, G. J. Comput. Phys. 1977, 23, 327−341. (70) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. J. Comput. Chem. 2009, 30, 2157−2164. (71) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K., 1987. (72) Hockney, R.; Eastwood, J. Computer Simulations Using Particles; McGraw-Hill Inc.: New York, 1981. (73) Martyna, G. J.; Tobias, D. J.; Klein, M. L. J. Chem. Phys. 1994, 101, 4177−4189. (74) Shinoda, W.; Shiga, M.; Mikami, M. Phys. Rev. B 2004, 69, 134103. (75) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182−7190. (76) Tuckerman, M. E.; Alejandre, J.; López-rendón, R.; Jochim, A. L.; Martyna, G. J. J. Phys. A 2006, 39, 5629−5651. (77) Gervois, A. Lect. Notes Comput. Sci. 2002, 2331, 95−104. (78) Medvedev, N. N.; Voloshin, V. P.; Luchnikov, V. A.; Gavrilova, M. L. J. Comput. Chem. 2006, 27, 1676−1692. (79) Kim, D.-S.; Cho, Y.; Kim, D. Comput.-Aided Des. 2005, 37, 1412−1424. (80) Dong, P. Comput. Geosci. 2008, 34, 411−421. (81) Rycroft, C. H. Voro++: A Three-Dimensional Voronoi Cell Library in C++; Lawrence Berkeley National Laboratory: Berkeley, CA, 2009; available at http://math.lbl.gov/voro++/. (82) Rycroft, C. H. Multiscale Modeling in Granular Flow. Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, MA, 2007. (83) Bondi, A. J. Phys. Chem. 1964, 68, 441−451.

4666


Definition and Computation of Intermolecular Contact in Liquids Using

Recommend Documents