Convergence of KirkwoodâBuff Integrals of Ideal and Nonideal

Subscriber access provided by READING UNIV

Article

Convergence of Kirkwood-Buff Integrals of Ideal and NonIdeal Aqueous Solutions Using Molecular Dynamics Simulations Jasmin Milzetti, Divya Nayar, and Nico F. A. van der Vegt J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b11831 • Publication Date (Web): 17 Jan 2018 Downloaded from http://pubs.acs.org on January 18, 2018

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

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

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

The Journal of Physical Chemistry

Convergence of Kirkwood-Buff Integrals of Ideal and Non-ideal Aqueous Solutions Using Molecular Dynamics Simulations Jasmin Milzetti, Divya Nayar, and Nico F. A. van der Vegt∗ Eduard-Zintl-Institut für Anorganische und Physikalische Chemie, Center of Smart Interfaces, Technische Universität Darmstadt, Alarich-Weiss-Straße 10, 64287, Darmstadt, Germany E-mail: [email protected]

Abstract

The computation of Kirkwood-Buff integrals (KBIs) using molecular simulations of closed systems is challenging due to finite system-size effects. One of the problems involves the incorrect asymptotic behavior of the radial distribution function. Corrections to rectify such effects have been proposed in the literature. This study reports a systematic comparison of the proposed corrections (as given by Ganguly et al. J. Chem. Theory Comput. 2013, 9 1347-1355 and Krüger et al. J. Phys. Chem. Lett. 2013, 4, 4-7) to assess the asymptotic behavior of the RDFs, the KBIs as well as the estimation of thermodynamic quantities for ideal urea–water and non-ideal modified-urea–water mixtures using molecular dynamics simulations. The results show that applying the KBI correction suggested by Krüger et al. on the RDF corrected with the Ganguly et ∗ To

whom correspondence should be addressed

1

ACS Paragon Plus Environment

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

Page 2 of 37

al. correction (denoted as B-KBI) yields improved KBI convergence for the ideal and non-ideal aqueous mixtures. Different averaging regions in the running KBIs (correlated or long-range) are assessed and averaging over the correlated region for large system sizes is found to be robust towards the change in the degree of solvent non-ideality and concentration, providing good estimates of thermodynamic quantities. The study provides new insights into improving the KBI convergence, suitability of different averaging regions in KBIs to estimate thermodynamic properties, as well as the applicability of correction methods to achieve KBI convergence for non-ideal aqueous binary mixtures.

Introduction The rigorous Kirkwood-Buff (KB) theory of solutions relates the microscopic structure of a multicomponent system to its macroscopic thermodynamic quantities. 1 KB theory is gaining increasing attention in the molecular simulation community as it can be used to extract thermodynamic quantities such as activity coefficients from radial distribution functions (RDF). The applications range from studies of binary mixtures 2–12 to studies of solvation and cosolvent effects on the conformational equilibria of peptides and macromolecules. 13–22 One of the challenges in using KB theory in molecular simulations is the computation of its central quantity, the Kirkwood-Buff integral (KBI), referred to as Gij . In the thermodynamic limit, V → ∞ (V is the volume), of a macroscopically large system with constant chemical potentials of solution components i and j, the KBI is defined as 1 Gij = V µVT

where gij

Z Z h V

V

µVT

gij

i (r ) − 1 dr1 dr2

(1)

is the grand-canonical RDF. Using the coordinate transformation r2 → r =

2


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


r1 − r2 , Eq. 1 reduces to Gij =

Z ∞h µVT

gij

0

i (r ) − 1 4πr2 dr

(2)

KBIs are frequently computed using Eq. 2 in molecular simulations of finite-sized, µVT

closed systems. Correspondingly, one then replaces gij

N pT

with gij

(or gijNVT ) and typically

considers a running integral up to distances R smaller than the half box dimension, i.e. GijR ( R)

=

Z Rh N pT 0

gij

i (r ) − 1 4πr2 dr

(3)

We shall refer to the running Kirkwood-Buff integral in Eq. 3 as R-KBI. Application of Eq. 3 may potentially introduce finite size errors, the magnitudes of which however depend on the chosen system size. The KBI (Eq. 1) characterizes cross fluctuations in the number of particles in an open µVT

system. Because fluctuations transform between ensembles, 23,24 replacing gij

N pT

by gij

is

formally not allowed and, in practice, leads to non-convergence of the R-KBI when the box size is small. When, however, a sufficiently large box is chosen, and the correlation lengths of the components’ density fluctuations are small compared to the linear dimension of the box, this problem disappears. The R-KBI then typically approaches a plateau value at sufficiently large R. Under these conditions, the volume of the simulation box surrounding the subvolume (4/3)πR3 in which the integral is calculated can be considered as a particle bath, mimicking an osmotic boundary condition for the subvolume considered. Based on this interpretation of the R-KBI in a sufficiently large box, Ganguly et al. 9 proposed a tail N pT

correction for gij

which ensures that the limiting behavior of the RDF corresponds to µVT

the limiting behavior of gij

(gij → 1 for large r).

It is interesting to ask if a converged R-KBI, obtained by means of this correction, approximates the KBI obtained from Eq. 1, calculated in a finite subvolume V = (4/3)πR3 embedded in an infinitely extended open system. Krüger et al. 25 have pointed out that for a finite subvolume V, the double integral in Eq. 1 cannot be reduced to a single one (Eq. 2)

3



by the transformation r2 → r = r1 − r2 . These authors provided the correct single-integral expression for the finite-volume KBI, in which 4πr2 in Eq. 3 is replaced by a geometrical function w(r, R) which we shall discuss later on. The corresponding finite, or small system, KBI, referred to as GijV , could be demonstrated to scale with the linear size R of V as 1/R. 25 Hence, the thermodynamically limiting KBI (Eq. 1) is obtained by linear extrapolation of GijV according to 1/R → 0. Different methods and corrections have been suggested in the literature to overcome the inaccuracies associated with the integral truncation or finite-size effects. 7,9,25,26 These corrections to KBIs and RDFs are explained in the theory section. The aim of this study is to compare different correction methods to estimate R-KBIs as well as their limiting values and to achieve improved convergence of the KBIs as a function of radius, herein denoted by G(R). We perform this study for ideal and non-ideal aqueous binary mixtures, which to the best of our knowledge have not been systematically investigated previously. Non-ideal solutions face the challenge of adequate sampling due to local inhomogeneities in the system and therefore G(R) suffers poor convergence. 7 Such systems require large simulation boxes and long simulation times, which has become computationally affordable in recent years; however, there is still a need to examine these systems for obtaining accurate KB integrals. We first examine the influence of different correction terms and simulation parameters, such as run length and box size, on the convergence of G(R) for nearly ideal behaving aqueous urea mixtures. The influence of non-ideality is studied by systematically scaling down the partial charges on all atoms of the original urea model. 3 Such a scaling weakens water-urea and urea-urea electrostatic interactions, leading to clustering of urea in the solution. This scaling approach provides a possibility to selectively tune the solution structure without any major changes in other properties such as the excluded volume of the particles. Our results provide an understanding of how to compute reliable KBIs and associated thermodynamic quantitities for ideal and non-ideal binary solutions over a concentration range.

4


Page 4 of 37

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


The organization of the paper, therefore, is as follows. First, different correction terms for RDFs and G(R)s and the calculation of relevant thermodynamic quantities are described, with an overview of prior molecular dynamics simulations studies reported. Second, the computational details are discussed, followed by the results for aqueous urea solutions and subsequently for modified urea solutions. At the end, main conclusions are summarized.

Theory KBIs from closed systems As described in the previous section, KBIs for closed systems are defined by Eq. 3 that can be derived from Eq. 1. It is often attempted to countervail the inaccuracies arising from N pT

the ensemble change in RDFs by corrections applied to the closed-system gij

(r ). The

asymptotic behavior of closed system RDFs at large distances, due to finite-size effects, can be explained as follows. If in an isobaric-isothermal ensemble there is an excess of particle type i in a volume VΩ around a central particle j, there will be a complementary depletion of i in the remaining volume Vtot − VΩ , with Vtot being the total system volume, because the total number of particles is fixed. In an analogous situation in the grand-canonical ensemble, this depletion is not apparent. Stemming from this explanation, Ganguly and van der Vegt have derived a correction to rectify the asymptotic behavior of the RDF in finite-sized systems, 9 Ganguly

gij

N pT

Nj f (r ) Nj f (r ) − ∆Nij (r ) − δij

(r ) = gij

(r )

f (r ) = 1 −

4 πr3 3 Vtot

(4)

with (5)

As stated in the paper, ∆Nij (r) is “the excess number of particles of type j within a sphere of radius r around particle type i” 9 and δij is the Kronecker delta. ∆Nij (r) is typically 5



Page 6 of 37

calculated from the R-KBI. The correction effectively only acts on the “tail” i.e. at large distances of the RDF. If the volume VΩ is small compared to the total volume, the effect explained above is negligible. This can be understood from the notion that at small radii r the simulation box effectively is a bath for the test volume VΩ . G(R) calculated from this corrected RDF in this work is denoted as G-KBI. An alternative RDF correction 27 based on the principle that deviations of the RDF from unity will decrease with

1 Ntot

Krüger

gij

is shown in Eq. 6,

(r ) = gijA (r ) −

1 1 c(r ) = gijB (r ) − c (r ) NA NB

(6)

with NA NB B A c (r ) = gij (r ) − gij (r ) NA − NB

(7)

This correction requires the knowledge of the RDF for two systems A and B with total number of molecules NA and NB , respectively, at the same thermodynamic state point. The system sizes NA and NB should be sufficiently large and different from each other in order to minimize the numerical errors associated with subtraction of the two nearly identical numbers. 27,28 In the literature, more RDF corrections can be found which will not be examined in this work. For instance, the software package GROMACS used for simulations in this study 29 provides a RDF correction, which asymptotically damps the RDF to unity using an exponential function and starting at a predefined radius. Perera et al. used a similar switch-function but of a tanh shape and chose an effective radius as starting point for the correction to operate. 7 Ben-Naim described the closed system RDFs to depend on “closure correlations” or the non-local correlations that become effective beyond the correlation distance. 30 This means that placing a particle at a fixed position in solution affects the average number of particles in the correlation volume. This changes the density of particles outside the correlation volume, affecting the asymptotic behavior of the RDF, which needs

6


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


to be corrected for. This argument is very similar to the reasoning by Ganguly et al. as both examine the influence of the fixed number of particles within the correlation volume, on the asymptotic behavior of the RDF. Based on the ideas of Ben-Naim, Cortes-Huerto et al. alternatively suggested a constant correction term for the RDF. 26 Applying the coordinate transformation r2 → r = r1 − r2 to Eq. 1 in order to derive Eq. 3 for systems with finite volume V leads to inaccuracies because the integration domain of r is not independent of r1 . To account for this, Krüger et al. first developed an expression for KBIs in finite volumes, termed as the "exact" finite KBI, which is defined by equation 6 in the original paper and denoted by GijV ( R), 25 GijV ( R)

=

Z R 0

w(r, R)

N pT gij (r ) − 1

dr

(8)

w(r, R) = 4πr2 (1 − 3r/2R + r3 /2R3 ) This finite KBI reduces to Eq. 2 (in this work) for R → ∞. The analytical weighting function denoted by w(r,R) accounts for the geometry of the finite volume. Here, this function has been derived for a spherical subvolume in three dimensions, however, other shapes of the subvolume are also possible. Sufficiently far away from the critical point, the finite KBI in Eq. 8 varies asymptotically with the inverse size 1/R. Therefore, the calculation of GijV ( R) allows to obtain the KBI (Eq. 1) by considering the limiting value of GijV ( R) for 1/R → 0. Alternatively, one can use the extrapolated expression of the exact KBI as given by equation 7 in the orginal paper of Krüger et al., 25 that we use in this work. This extrapolated KBI will be referred to in our work as Krüger-KBI or K-KBI: GijK ( R)

=

Z R 0

N pT ω (r, R) gij (r ) − 1 dr

(9)

with ω (r, R) = 4πr

2

1−

7

r 3 R


(10)


Page 8 of 37

Note the change in notation compared to equation 7 in the original paper. 25 Here, ω (r, R) in terms of the notation in the original paper is defined as, r 3 ∂w 2 ω (r, R) = (w − x ) = 4πr 1 − ∂x R

(11)

with R being the large distance and w(r,x) depending on x = r/R. The function ω (r, R) depends on the upper integration limit R and always damps the oscillations of the RDF close to this limit, i.e. ω (r → R) = 0. This function in three-dimensions is derived for spherical subvolumes, but other shapes of the subvolume are also possible. Unlike the R-KBI, the K-KBI does not depend on the position of the central particle and the average density of a finite-sized volume (4/3πR3 ) will be equal to the average macroscopic density ("density effect"). 25 As Krüger et al. have shown for an analytically described test system, applying ω (r, R) leads to the suppression of fluctuations in the RDF that would otherwise be amplified in G(R). Furthermore, it leads to a shift of the phase with which G(R) is fluctuating. As will be illustrated in the results section, the improved convergence of the K-KBI over the R-KBI is not only associated with the density effect. The weighting function also attenuates the influence of the tail region of the RDF, which is subject to greater statistical deviations from unity due to limited sampling, on the final value of the KBI. To treat both effects impairing the convergence of G(R), we define a new correction. Since the Ganguly RDF tail correction improves the convergence of G(R) (discussed later in Figure 2), we apply the Krüger correction for integration (Eq. 9) to the Ganguly corrected RDFs (Eq. 4). We denote KBIs with both corrections applied as B-KBIs. Other methods not relying on correction terms are also available in the literature for the calculation of the KBI. A method that is mostly used in engineering applications is the fitting of the structured part of the RDF 31 or the KBI 32 to obtain improved results over RKBIs. An equivalent approach to the computation of Eq. 8 (with the ω (r, R) corresponding

8


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


to the same shape as that of the subvolume) is to directly compute the particle number fluctuations in inserted sub-boxes as suggested by Schnell et al. 33,34 These methods were implemented and tested by Ganguly and van der Vegt 9 to urea-water and methanol-water systems. A set of sub-boxes of different sizes is inserted in each frame of the trajectory and the finite KBI is calculated from the particle number fluctuations in the sub-box. The finite KBI is then plotted as a function of the inverse radius or length of the test box. The KBI for a infinite sized system can be obtained as the extrapolated intercept of a linear fit. A drawback of this approach is that the range of linear scaling with 1/R is limited, 25 as observed in applications of, e.g., simple Lennard-Jones particle mixtures. 35 Thus, this method is not further discussed in this work.

Thermodynamic properties of binary mixtures The following expressions have been derived for binary mixtures. 30 The notation is adopted for the system at hand in this work: The subscripts “U” and “W” refer to urea or modified urea and water respectively. All KBIs Gij are assumed to be given in cm3 mol−1 , so that the molar concentration c in mol cm−3 has to be used. The isothermal compressibility κ T can be calculated as: 30 κ T = ( R T ) −1

2 ) 1 + cU GUU + cW GWW + cU cW ( GWW GUU − GUW cW + cU + cW cU ( GWW + GUU − 2GUW )

(12)

where R is the gas constant and T is the temperature. The isothermal compressibility is accessible through volume fluctuations, too: 24

κT =

1

hV i N pT k B T

hδV 2 i N pT

hδV 2 i N pT = hV 2 i N pT − hV i2N pT

(13) (14)

where h· · · i N pT denotes the isobaric-isothermal ensemble average. The derivative aUU 9



Page 10 of 37

of the urea activity aU on a molar concentration scale in a closed binary system, with yU being the corresponding activity coefficient can be calculated as shown in Eqs. 15 and 16. aUU = yU ( cU ) = −

∂ ln aU ∂ ln c˜U Z cU 0

dc˜U

= 1+ p,T

∂ ln yU ∂ ln c˜U

= p,T

c˜U ( GUU − GUW ) 1 + c˜U ( GUU − GUW )

1 1 + c˜U ( GUU − GUW )

(15) (16)

The KBIs can also be interpreted rather directly. They are a measure of the molecular attraction. Most often, the excess number ∆Nij is calculated, which indicates the “change in the average number of i molecules in the sphere of radius R, caused by placing a j molecule at the center of the sphere. This interpretation applies to any R.”. 36 The excess coordination number is related to the KBI and the particle number density for R → ∞, according to; ∆Nij = Gij c j

(17)

Overview of previous KBI calculations Several parameters influence the convergence of the G(R), some of which have already been discussed in the literature. 9 The most obvious are the box size and the simulation time. The larger the box, the less important the finite-size effect becomes. On the other hand, for larger boxes equilibration and sampling of the configuration space requires longer trajectories and more computational effort. The KBIs for urea-water system have been reported in the literature before. In the original paper introducing the urea force field, which will also be used in this study, Weerasinghe and Smith used 4 nm boxes and short simulation times of 2.5 ns. 3 They averaged uncorrected R-KBIs (Eq. 3) in the window from 0.95 nm to 1.30 nm, which is a region spanning one period in the fluctuations of the RDF. They used the same method of averaging also in the study of other mixtures such as methanol-water (6 nm, average between 0.95 nm and 1.20 nm) 5 and acetone-water (most mixtures: 2-3 ns, average 0.85 - 1.25 nm). 4 Slightly larger boxes (4.4 nm to 4.9 nm) and 10


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


longer runs (20 ns) were used by Chiba et al. 12 De Oliveira et al. used large boxes of ≈ 8 nm box length and 25 ns trajectories. 37 Both groups used the same method of averaging as suggested by Weerasinghe and Smith. Horinek et al. used 5 ns trajectories from ≈ 4 nm box length systems for the calculation of activity coefficients. 15 Correction terms and subbox methods have been mostly used in the development or comparison of new methods. Ganguly and van der Vegt 9 compared R-KBIs and G-KBIs for small (4 nm) and large (≈ 7 nm) systems and long trajectories of 100 ns with the small-box-method introduced by Schnell et al. 33 using aqueous urea and methanol-water mixtures as testsystems. For other than methodological studies, the usage of elaborate corrections and sub-box methods 38 is rarely seen independent of the system. Ploetz et al. used the subbox method in a study analysing the influence of the polarisability of the force field on different macroscopic properties for aqueous methanol. 39 To the best of the authors’ knowledge, the extrapolated KBI-correction (Equation 9) introduced by Krüger et al. has not been used before. There have been, however, previous studies making use of the exact finite KBI correction introduced by the same authors. 10,40,41 In this work, the focus will be on comparing the correction methods by Krüger et al. and Ganguly et al.

Computational Details System setup and simulations MD simulations of aqueous urea mixtures were carried out using the SPC/E water force field 42 and the KB urea force field. 3 The urea force field reproduces the ideality of aqueous urea solutions reasonably, 3,15,43 as it has been designed to reproduce KBIs in mixtures of concentrations between 2 M and 8 M. Urea bond lengths were constrained with the LINCS algorithm; 44 for water molecules SETTLE was used. 45 All simulations were performed using Gromacs 4.6.7. 29 The leap-frog integrator was used with a timestep of 2 fs. Periodic boundary conditions were applied in all three 11



directions. In the original force field a twin range cut-off between 0.8 and 1.5 nm was used. 3 For technical reasons, the cut-off was modified and a plain cut-off at a distance of 1.2 nm was applied for both the electrostatic (real-space) interactions and van der Waals interactions. The influence of the cut-off on the KBI has been examined (see below and SI). The electrostatics were treated using the particle mesh Ewald method. 46 Long range corrections were applied to energy and pressure. In the development of the force field Weerasinghe and Smith varied the atomic partial charges to reproduce the nearly ideal behavior of the urea-water mixture. In this study, the absolute values of all partial charges of all atoms C, O, N and H in urea were scaled down to modify the balance of interactions (urea-urea and urea-water) towards non-ideal behavior in steps from full charge (100 %) to 70 %, keeping the van der Waals interaction parameters the same as in urea. We denote these systems as modified urea or ’m-urea’ systems and specify the partial charge in XX % of the original charge as U-qXX. A range of concentrations of aqueous solutions, 0.5 to 9 M were studied. The number of particles used in the simulations is given in Table 1. Some systems could not be simulated because the solutions became unstable. A full overview over all simulations performed and parameters is given in the SI (Tables S1, S2, S3). In general, cubic boxes with a side length of about 7 nm (“large”) were used. In addition, some simulations were performed with “small” boxes of 4 nm box length. The following protocol was used for equilibration for one system of each concentration: At first, an energy minimization using a steepest descent method was performed. This was followed by equilibration at 300 K and 1 atm, under NVT conditions using the Berendsen thermostat 47 (τT = 2.0 ps) for 1 ns and under NPT conditions using the Berendsen thermostat and barostat (τp = 1.0 ps) for 100 ps. Further equilibration for 1 ns was done with ParrinelloRahman barostat 48 (τp = 1.0 ps) and Nosé-Hoover thermostat 49 (τT = 1.0 ps), which were also used for production runs. The urea force field parameters were then changed to the targeted partial charge. The first 10 ns of the production runs were discarded for

12


Page 12 of 37

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


Table 1: Overview of concentrations of the systems for which the simulations were performed. ctarget denotes the target concentration and c is the concentration averaged over all NPT simulations of all systems of aqueous urea and modified urea. NU and NW are the number of urea and water molecules, respectively. The upper part of the table lists simulations in boxes of 7 nm side length; the lower part describes details for the 4 nm systems. 7 nm ctarget [M] 0.5 1 1.5 2 3 4 5 6 7 9 4 nm ctarget [M] 1 2 5 6

c [M]

NU

NW

0.486 0.983 1.481 1.985 3.013 4.081 5.159 6.208 7.229 8.715

103 207 310 413 620 826 1033 1239 1446 1652

11475 11131 10787 10442 9754 9065 8377 7800 7373 6311

c [M]

NU

NW

0.981 1.990 5.181 6.308

39 77 193 231

2077 1948 1563 1435

equilibration due to the change in partial charge. 50 ns runs were used for further analysis. The frames were written every 2000 steps (4 ps, 12500 frames total). For the calculation of the water-water RDF every 5th frame (every 20 ps) was used.

Computation of the RDFs and KBIs All the pair correlations presented in this work have been calculated using the carbon atom as reference in urea and the oxygen atom in water. As has been tested, this procedure leads to nearly the same RDFs as taking the center of mass (COM) because the displacement

13



between the central atoms and COMs are small. For the KBI it is irrelevant which site is taken. 50 Most commonly, KBIs are calculated by numerically integrating RDFs. These RDFs are generated from evaluating pair-pair distances. In the “direct approach”, which has also been described by Chiba et al., 12 the number of atoms in concentrical spheres is counted right away. This reduces numerical errors. However, no noteworthy differences between the methods could be observed in this work. The uncorrected G(R) (R-KBI) is calculated by the direct method. The corrected G(R)s (G-KBI, K-KBI, B-KBI) cannot be calculated in this way because the methods require modifications to the RDF and thus the integration method is chosen. Despite using different integration methods for different KBIs, no inconsistencies are expected as no difference between the methods was apparent for the R-KBI. The bin size is 0.005 nm, for which the RDF was found to be insensitive to minor changes. For computation of average KBI values, two common approaches are used. First, the structured regime or the "correlated regime" (CR) of the R-KBI is considered, typically between two maxima in the RDF around 1.0 nm. 3,5 Such a region has been chosen to lie between 0.97 to 1.31 nm in our work after examining the RDFs of the systems. In the second approach, called the "long range" (LR) method, the thermodynamic limiting value is taken from the (converged) region with zero slope at large distance in the tail of a corrected KBI (unlike the commonly used R-KBI). We chose the B-KBI for this purpose as it exhibits superior convergence as will be discussed later. If no such point can be identified, an average is taken over a window of certain distance in the tail of the B-KBI. For this purpose, different averaging windows have been tested and 2.0-2.5 nm has been selected for large box systems. This region shows plateau values such that final values are insensitive to small variations (see Figure 3, which will be discussed later). For the small systems, the average between 1.2-1.6 nm is taken. Note that there are no oscillations in the B-KBI at these radii although the R-KBI shows oscillations up to radii of 1.9 nm. In our study, these

14


Page 14 of 37

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


two approaches are examined for the accuracy in estimation of thermodynamic quantities. Taking into account the different box sizes, the four methods are denoted as 4nm LR, 4nm CR, 7nm LR and 7nm CR.

Choice of simulation parameters Figure 1 shows checks performed during the choice of simulation parameters for a 5 M aqueous urea solution. Improvement in the convergence of R-KBIs and B-KBIs can be observed for simulation lengths up to 50 ns, but no significant further improvement is observed for the 70 ns simulation. The lower panel shows that larger boxes lead to more converged tail regions in the B-KBI. Furthermore, differences in the correlated region (up to 1.3 nm) in uncorrected KBIs (R-KBI) can be seen. Figure S1 (in the supporting information) depicts the influence of cut-offs on R-KBIs and B-KBIs. A value of 1.0 nm underestimates the KBI, while cut-offs of 1.2 nm and 1.5 nm show consistent results with better convergence than with value of 1.0 nm. The value of 1.2 nm is, therefore, chosen due to higher computational affordability and convergence of gUU (r) to unity at this distance, as also reported by Weerasinghe et al. 3

Results and Discussion Effect of corrections on the convergence of RDFs and KBIs The effect of RDF corrections by Ganguly et al. (Eq. 4) and Krüger et al. (Eq. 9) on the urea-urea RDF (gUU (r )) tail convergence and corresponding GUU ( R) convergence for 5 M aqueous urea solution (7 nm box) are compared in Figure 2 (corresponding results for 2 M urea in Figure S4). Both the corrected RDFs converge to a value of one, although the Krüger correction introduces a substantial fluctuation in values. The larger noise in Krüger corrected RDFs may be attributed to the method of computation. As described above, it

15



Figure 1: Influence of choice of simulation parameters. The upper and lower panels show the dependence of GUU on simulation run length and the box length respectively for 5 M urea solution. involves the subtraction of two large numbers. Moreover, as the RDFs for two system sizes are subtracted, the deviations from unity due to limited sampling will add up. The Ganguly correction performs well in terms of retaining the shape of gUU (r ) and improving the convergence of both the gUU (r ) and GUU ( R). The effect of the Ganguly correction is more prominent at low urea concentrations (2 M), where insufficient sampling causes poor convergence of the uncorrected RDF and the corresponding GUU ( R) (see Figure S4). The good performance of the Ganguly RDF correction is supported by the decreasing magnitude of difference between corrected and uncorrected RDFs, ∆gUU , as well as the root mean squared deviation of the corrected RDF from the expected value of unity (σ(gUU )) with increasing urea concentration (see Figure S5). Both ∆gUU and σ(gUU ) are close to zero beyond 4 M concentration. The performance of different correction methods (R-KBI vs. G-KBI, K-KBI and B-KBI) in improving Gij (R) convergence is assessed next. The uncorrected urea-urea R-KBI shows poor convergence, especially in the case of the 2 M system due to poorer sampling than in the 5 M system (see Figure 3). The effect of the Krüger correction is evident in the 16


Page 16 of 37

Page 17 of 37

1.010

(a)

gUU(r)

1.005

● 1.000

0.995 no correction Krüger (rdf) Ganguly

0.990 1.5

2.0

2.5

3.0

distance (nm) 50

GUU (cm3 mol−1)

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


0

−50

−100

(b) ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ●● ● ● ●● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●

no correction Krüger (rdf) Ganguly

−150 0.0

0.5

1.0

1.5

2.0

2.5

3.0

distance (nm)

Figure 2: Influence of different RDF corrections on (a) the RDF, gUU (r) and (b) the corresponding KBI, GUU for 5 M aqueous urea (7 nm system size). Note that the Krüger corrected GUU ( R) shown here has been computed using Krüger corrected RDF (Eq. 6) which was used in Eq. 3 to obtain the R-KBI. displacement of peaks, faster convergence and more converged tail regions in Gij ( R) as compared to R-KBI, as has been indicated in the theory section. G-KBI and B-KBI show remarkable improved convergence, especially for the 2 M system, with the B-KBI converging faster than the G-KBI. This is also evident from the root mean square deviations of the GUU from the mean value computed between 2.0 and 2.5 nm, as shown in Figure 4. The deviation is close to zero for concentrations above 2 M and shows higher deviations at lower concentrations, due to larger fluctuations in the tail of GUU (R) as discussed before. The deviation is smallest in the case of the B-KBI at all concentrations indicating reliable performance of this method to obtain improved convergence in the KBIs for nearly 17



Figure 3: Effect of correction terms on the KBI convergence of (a) 2 M and (b) 5 M aqueous urea solutions (7 nm system size). ideal aqueous urea solutions. The effect of the Krüger correction in inducing faster KBI convergence can be seen in the overlap of the tail of B-KBI (or K-KBI) with the correlated region of R-KBI between 1.3 and 1.5 nm distance. Beyond 1.5 nm, B-KBI (and K-KBI) shows improved convergence. This is due to the ω (r, R) function in the extrapolated KBI that dampens the effect of the oscillations in the RDF at large distances caused by the long-range correlations due to weak urea-water interactions as discussed later. Average KBI values for aqueous urea solutions of different concentrations obtained from the CR and LR methods are shown in Figure 5. GWW agrees well with experimental values at all urea concentrations. GUW and GUU at lower concentrations, however, exhibit deviations from experimental values. 4nm CR shows systematically higher GUW values 18


Page 18 of 37

Page 19 of 37

40 ●

σ(GUU) / cm3mol−1

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


R−KBI G−KBI K−KBI B−KBI

●

30

●

●

20 ●

10 ● ● ●

●

●

●

0 0

2

4 6 [urea] (mol l−1)

8

Figure 4: Root mean squared deviation (rmsd) in the tail (2.0 nm to 2.5 nm) of the GUU relative to the mean value in this region for different KBIs of aqueous urea systems. and better agreement with experiment than the three other methods. However, as the force-field was parametrized using the 4 nm CR method and the three other methods all lead to the same values, it can be concluded that 4nm CR probably overestimates GUW rather than 4nm LR, 7nm CR and 7nm LR underestimating it. The 7nm CR method produces GUU values closest to the experimental data at all urea concentrations and 7nm LR overestimates the values at low concentrations. It is intriguing that the two methods 7nm CR and 7nm LR give similar GUU values at higher urea concentrations, although these methods are based on different approximations. This is inherent in the nature of corrections applied to the KBI, that are discussed later in detail. The methods for small system size perform well except that 4nm CR underestimates the GUU values at low urea concentrations. Deviations of GUU from experimental data at low concentrations does not strongly affect aUU , as it will be weighted with the concentration (see Eq.15 and discussion below). At concentrations above 4 M, the values of GUU obtained from all the methods agree well with experiment.

19



Figure 5: Comparison of different methods for the calculation of KBI values for aqueous urea solutions as a function of concentration. The CR and LR data corresponds to that obtained from R-KBI and B-KBI respectively. The dimensions on all y-axis are cm3 mol −1 . Black lines and dots are experimental data as compiled by Chiba et al. 12

Estimating thermodynamic quantities We next assess the performance of the correction methods to estimate the thermodynamic quantities. The urea activity derivatives (aUU ) computed using different methods, as shown in Figure 6 (a), are in agreement with each other except for the 4nm CR method, which gives higher aUU values at all urea concentrations. As Weerasinghe et al. used this method for parameterization, it is unexpected that the associated aUU values show the largest offset from the experiment. This may be due to the differences in barostat/thermostat used, as Weerasinghe et al. used weak coupling schemes. 3 The qualitative trend in aUU computed from all the methods is in agreement with experiment. As has been reported before, 3,37 there are two sets of experimental data 51,52 and the qualitative trend of the 20


Page 20 of 37

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


data by Stokes is more accurately depicted. Quantitatively, the deviation of aUU obtained from different methods is within 1% and 24% from the experiment. Considering the data presented in Figures 5 and 6 The CR methods perform slightly better than the LR methods since aUU effectively depends on the short-range molecular interactions which are well-captured by the CR methods. Among the 7nm CR and 4nm CR, the former performs better in terms of estimating GUU and GUW especially at low concentrations and therefore will be used to estimate the Gij for non-ideal m-urea systems in the further discussion. Figure 6 (b) compares the isothermal compressibility (κ T ) values as obtained from different KBI methods with those computed from volume fluctuations. The values obtained from volume fluctuations are in close agreement with the experimental values at low urea concentrations and deviate at higher urea concentrations. The KBI methods do not agree well among each other and difficulties are observed in reproducing the values obtained from volume fluctuations. The 7nm LR method, however, outperforms the CR methods in estimating the κ T with values closer to those obtained from volume fluctuations and experiment. This is expected since the κ T depends on the macroscopic density fluctuations in the system, which are determined by long-range correlations. Difficulties in calculation of compressibilities from KBIs have also been observed previously by Weerasinghe and Smith 3 and Chiba et al. 12

Effect of solvent non-ideality on KBI estimation The effect of modifying the balance of urea-water and urea-urea interactions on the KBIs is examined in Figure 7. Reduced urea-water interactions enhance urea clustering, impeding GUU R-KBI convergence as shown for the 5 M U-q85 system in Figure 7 (a) (compare urea 5 M, Figure 3). Severe aggregation as for U-q75 systems or lower concentrations enhances the problems in G(R) convergence, as shown in Figure S6. However, B-KBIs still show improved convergence over R-KBIs. Figure 7 (b) shows Gij values obtained using LR (B-KBI) and CR (R-KBI) methods for both system sizes for U-q85 systems. At low m-urea 21



1.6

7nm CR 7nm LR 4nm CR 4nm LR exp. 1 exp. 2

aUU

1.4

(a)

1.2 ideal

1.0 0.8 0

10

2

7nm CR 7nm LR

4

-1

6

[urea] (mol l ) 4nm CR 4nm LR

8

V-fluct exp.

(b)

8

κT (10-10Pa-1)

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

Page 22 of 37

6 4 2 0 0

2

4

6

[urea] (mol l-1)

8

Figure 6: Comparison of thermodynamic quantities computed from different methods. (a) shows the urea activity derivatives calculated using different KBI methods. The CR and LR data corresponds to that obtained from R-KBI and B-KBI respectively. The experimental data is taken from (1) Stokes 51 and (2) Miyawaki et al., 52 as reprinted in [3]. (b) shows isothermal compressibilities calculated from volume fluctuation and KBI method. Experimental data taken from reference [53]. concentrations they do not agree well with each other, with large error bars at low m-urea concentrations due to poor sampling in urea-agglomerated systems. Lack of experimental data on the KBIs for such systems makes it difficult to predict the accuracy of the LR and CR methods in estimating these thermodynamic quantities. Since the fluctuation in the

22


Page 23 of 37

GUU values is least for 7nm CR method at low m-urea concentrations, it can be used as a reliable method for estimating thermodynamic quantities for such non-ideal systems. At higher concentrations, the error bars in GUU reduce for LR and CR methods and the values obtained for two system sizes and both LR and CR, agree with each other, with slight underestimation of values for the small system sizes. (a) GUU (cm3 mol−1)

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


200 150 100 50 0 −50 −100 −150

R−KBI G−KBI K−KBI B−KBI

0.0

0.5

1.0

1.5

2.0

2.5

3.0

distance (nm)

Figure 7: (a) Effect of correction methods on the KBI convergence for U-85 at 5 M. (b) KBI values obtained from different methods for U-q85 systems. The CR and LR data corresponds to that obtained from R-KBI and B-KBI respectively. Figure 8 summarizes the average Gij values for the m-urea systems. The trends in GUU 23



and GUW with increasing solvent non-ideality (decreasing atom partial charges), agree well with the expected behavior of the non-ideal solutions examined here. The weakened attractive urea-water interactions result in aggregation of urea and depletion of water from urea hydration shells (as also reflected in decreasing ∆NUW excess numbers, see SI Figure S9 (b)), leading to increasing GUU (and ∆NUU , Figure S9 (a)) and decreasing GUW values with decreasing partial charges. This causes water to engage stronger with itself, resulting in increasing GWW values with increasing solvent non-ideality (and increasing ∆NWW , Figure S9 (c)). GUU diverges for U-q75 and U-q70 beyond 2 M concentration, where the structural micro-heterogeneities increase due to urea clustering in the system. GUU decreases and GWW increases with concentration for all (m-)urea systems. The GUW as a function of (m-)urea concentration shows a non-monotonic trend for urea and U-q95, with a minimum around 4 M concentration. For other m-urea systems, GUW monotonically decreases. The urea activity derivatives (aUU ), shown in Figure 9 computed with the 7nm CR method, show a non-monotonic dependence on urea concentration with a minimum around 3 M for high partial charges. For nearly ideal urea solutions, aUU is close to one till 5 M concentration, beyond which deviations are observed. With increasing solvent non-ideality, aUU decreases and the non-monotonic concentration dependence disappears for partial charges scaling lower than 85%. The aUU reaches a value close to zero for U-q70, denoting instability of such systems. The performance of different methods in estimating aUU for U-q85 systems is shown in Figure S7 (a) of the SI. The values obtained for the 4 nm and 7 nm system sizes agree well at low concentrations, with 4nm CR showing an offset. At higher concentrations, the offset between the methods for the two system sizes increases. Figure S7 (b) compares the robustness of R-KBI and B-KBI derived values of aUU . The isothermal compressibilities for m-urea systems as computed with the 7nm LR method are compared with those computed using volume fluctuations in Figure S8 (a). Although in many cases the volume fluctuation values can be reproduced within error bars

24


Page 24 of 37

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


Figure 8: The average Gij for UU, UW and WW for urea and m-urea systems at different concentrations as computed using 7nm CR (R-KBI) method. by the KBI method, in general the sensitivity of this method is insufficient to accurately calculate compressibilities. Figure S8 (b) shows that especially at high concentrations the LR methods perform better than CR methods. Similar results could be found for urea (Figure 6). The applicabiliy of KBIs to quantify the degree of aggregation in non-ideal aqueous mixtures has been shown in a previous study. 43 In support, our results highlight

25



the sensitivity of different correction methods as well as averaging methods in providing converged and accurate KBIs for such systems.

Figure 9: (a) Variation of urea activity derivatives for urea and m-urea, as computed using 7nm CR, with different partial charge scaling and urea concentration. (b) Snapshots of some systems of 5 M concentration. They are marked with rectangles in (a) except for U-q60 which represents an unstable system.

Discussion and Conclusions The application of Kirkwood-Buff (KB) theory in molecular simulations has rapidly increased in recent years, yet the accurate calculation of KB integrals from MD trajectories, especially in the case of isothermal-isobaric conditions, remains a debatable subject. In the literature, corrections have been proposed which lead to an improved convergence of G(R). In this work, different RDF corrections as well as the boundary correction introduced by Krüger et al. 25 have been tested for binary aqueous mixtures. The RDF correction by Ganguly and van der Vegt 9 has been found to be most effective. As Krüger et al. also suggested in their original paper, it has been found that it is necessary to apply RDF corrections prior to the use of the “Krüger correction” (Equation 9). Consequently, our results 26


Page 26 of 37

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


illustrate that the B-KBI (both corrections) performed better than the sole application of the Krüger correction (K-KBI) or the RDF correction by Ganguly and van der Vegt (G-KBI). By combining two different system sizes (4 nm and 7 nm box length) with two approaches to compute KBIs (CR and LR), four different methods were compared. Additionally, applicability of the corrections in estimating KBIs and thermodynamic properties of non-ideal systems has been examined. The methods have been tested for a series of ideal to non-ideal solutions of varying concentrations. To this end, we uniformly scaled down the partial charges on the urea model developed by Weerasinghe and Smith, known to form nearly ideal aqueous solutions. 3 This has provided new insights into existing understanding of the G(R) convergence for such systems. 43 We analyzed the performance of each KBI method by examining the quantities Gij , aUU and κ T . The performance of methods for accurately computing Gij and aUU could be identified, computation of κ T , however, was found to be highly sensitive to variations in KBIs. Difficulties have been reported in previous studies as well. 3,12 Activity derivatives aUU proved meaningful to compare different methods. Direct comparison of results obtained by different methods used in this work with the experimental data may not be the best approach since the original force field has been parametrized using 4nm CR. However, trends could be extracted from comparison among methods for all urea and m-urea systems. Our results show that the computational effort can be minimized by using small boxes and B-KBI method i.e. 4nm LR to estimate accurate thermodynamic quantities for nearly ideal aqueous mixtures at all concentrations. However, for the systems where convergence of G(R) is challenging due to solvent non-ideality and/or low concentrations such difficulties can be circumvented by using big boxes and using the correlated regime in R-KBIs (7nm CR). It can safely be assumed that the smaller boxes are not large enough to accommodate and accurately depict micro-heterogeneities, which are the molecular basis for non-ideal behavior. This trend was expected and has been reported before

27



such as by Perera et al. for aqueous methanol and methanol-acetone mixtures and by Weerasinghe and Smith for aqueous acetone. 4,7 Inadequate sampling at low concentrations and high computational effort associated with elongating simulation run lengths restrict the performance of correction terms for both ideal and non-ideal systems. This has been exemplified in Figure S6. Therefore, LR methods are found to be unsuitable to compute properties at low concentrations. We, therefore, conclude that the 7nm CR method yields consistent results across concentrations and degrees of ideality for binary aqueous mixtures studied in this work. Interestingly, at higher urea concentrations, the 7nm CR (R-KBI) and 7nm LR (B-KBI) methods yield similar results despite that the former method does not involve any corrections to the KBI and the averaging is done in the correlated regime of the KBI. On closer examination, this contradiction is resolvable. The analytical function ω (r, R) in the B-KBI dampens the influence of the RDF at large distances. It suppresses the long-range oscillations in the resulting KBI, which converge very slowly on a typical MD time scale and are caused by long-range correlations due to the weak urea-water fluctuating hydrogen bonding interactions. At smaller distances, the RDF correction does not have significant influence on the KBI convergence, as can be observed by comparing the R-KBI and G-KBI or equivalently K-KBI and B-KBI in Figure 3. This is reasonable because for probe sub-volumes with smaller radii, the surrounding bulk can effectively act as a particle bath emulating the grand-canonical ensemble. For instance, the sub-volume used to compute R-KBIs with r = 1.3 nm corresponds to 2.6 % of the total volume in the 7 nm box, effectively reducing the finite-size effects, consistent with the observations in a previous study. 26 As a consequence, it is reasonable to average over one oscillation in the correlated regime of the KBI, resulting in values similar to those obtained by the B-KBI LR methods. Interestingly, the KBI values obtained using our 7nm CR and 7nm LR methods coincide with those obtained by extrapolation of the finite KBI (Eq. 8 with Ganguly corrected RDF) with 1/R → 0, especially for aqueous urea systems at high concentrations (see Figure 10). The conclusions on the reliability of the methods are summarized in Table

28


Page 28 of 37

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


2 below.

Figure 10: Urea-urea (top panel), urea-water (middle panel) and water-water (bottom panel) exact, finite KBI (GV , Eq. 8) as a function of inverse distance. The ’100%’ and ’85%’ labels denote the aqueous urea and aqueous U-q85 m-urea systems. The point symbols indicate the average KBI values obtained using different methods used in this work.

29



Table 2: Description of the performance of different methods for the calculation of the KBI. The [c] denotes the concentration range in molarity for aqueous binary mixtures studied in this work. The nearly ideal mixtures are aqueous urea and U-q95 up to 6 M concentration. The non-ideal mixtures are aqueous urea above 6 M and U-qxx m-urea systems with xx below 95%. CR

LR Corrected KBI (B-KBI) with Ganguly Uncorrected KBI (R-KBI) with the (Eq. 4) and Krüger (Eq. 9) correcaverage taken in the structured (cortions. The average is taken in the related) region of 0.97 - 1.31 nm converged region at large radii.

4 nm

Results are consistent with other methods and meet expected trends only for nearly ideal mixtures at higher concentrations ([c] ≥ 4 M, urea or U-q95).

Can be used for nearly ideal mixtures independent of concentration and can substantially save associated computational effort.

7 nm

Robust to both changes in concentration and degree of non-ideality of the system. Yields similar results as 7nm LR for urea systems at higher concentrations ([c] ≥ 4 M). Allows using shorter trajectories than 7nm LR method, saving the computational effort. The region for averaging has to be chosen for every system individually to span one oscillation of the RDF.

Improved KBI convergence and accurate determination of thermodynamic quantities if the sampling is sufficient, i.e. at high concentrations ([c] ≥ 4 M). Longer trajectories are necessary for non-ideal mixtures. Poor convergence (a slope of > 10cm3 mol−1 nm−1 of the B-KBI in the range 2.0 nm to 3.0 nm) for low concentrations ([c] < 4 M).

30


Page 30 of 37

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


Supporting Information: The Supporting Information file describes system details, computation of the error bars, additional results on influence of cut-off and simulation length on KBI estimation, results on the RDFs of urea and m-urea systems, comparison of RDF tail correction methods for 2 M urea system, performance of Ganguly RDF correction on all the systems studied, comparison of KBI convergence computed using different methods for m-urea systems, thermodynamic quantities for m-urea systems (urea activity coefficient derivatives and isothermal compressibility), comparison of the excess numbers for nearly ideal and non-ideal solutions.

Acknowledgement The authors thank David Rosenberger for insightful discussions. Computations for this research were conducted on the Lichtenberg high performance computer of the TU Darmstadt. Part of the research described herein was supported by the German Research Foundation (DFG) within the Collaborative Research Center "Interaction between Transport and Wetting Processes" (SFB 1194).

31



References (1) Kirkwood, J. G.; Buff, F. P. The statistical mechanical theory of solutions. I. J. Chem. Phys. 1951, 19, 774–777. (2) Chitra, R.; Smith, P. E. Molecular Association in Solution: A Kirkwood-Buff Analysis of Sodium Chloride, Ammonium Sulfate, Guanidinium Chloride, Urea, and 2,2,2Trifluoroethanol in Water. J. Phys. Chem. B 2002, 106, 1491–1500. (3) Weerasinghe, S.; Smith, P. E. A Kirkwood - Buff Derived Force Field for Mixtures of Urea and Water. J. Phys. Chem. B 2003, 107, 3891–3898. (4) Weerasinghe, S.; Smith, P. E. Kirkwood-Buff Derived Force Field for Mixtures of Acetone and Water. J. Chem. Phys. 2003, 118, 10663. (5) Weerasinghe, S.; Smith, P. E. A Kirkwood-Buff Derived Force Field for Methanol and Aqueous Methanol Solutions. J. Phys. Chem. B 2005, 109, 15080–15086. (6) Lee, M.-E.; van der Vegt, N. F. A. A New Force Field for Atomistic Simulations of Aqueous Tertiary Butanol Solutions. J. Chem. Phys. 2005, 122, 114509. (7) Perera, A.; Zorani´c, L.; Sokoli´c, F.; Mazighi, R. A comparative Molecular Dynamics study of water–methanol and acetone–methanol mixtures. J. Mol. Liqu. 2011, 159, 52–59. (8) Ganguly, P.; Mukherji, D.; Junghans, C.; van der Vegt, N. F. A. Kirkwood-Buff CoarseGrained Force Fields for Aqueous Solutions. J. Chem. Theory Comput. 2012, 8, 1802– 1807. (9) Ganguly, P.; van der Vegt, N. F. A. Convergence of Sampling Kirkwood - Buff Integrals of Aqueous Solutions with Molecular Dynamics Simulations. J. Chem. Theory Comput. 2013, 9, 1347–1355.

32


Page 32 of 37

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


(10) Schnell, S. K.; Skorpa, R.; Bedeaux, D.; Kjelstrup, S.; Vlugt, T. J. H.; Simon, J.-M. Partial Molar Enthalpies and Reaction Enthalpies from Equilibrium Molecular Dynamics Simulation. J. Chem. Phys. 2014, 141, 144501. (11) Ploetz, E. A.; Smith, P. E. Particle and Energy Pair and Triplet Correlations in Liquids and Liquid Mixtures from Experiment and Simulation. J. Phys. Chem. B 2015, 119, 7761–7777. (12) Chiba, S.; Furuta, T.; Shimizu, S. Kirkwood-Buff Integrals for Aqueous Urea Solutions Based upon the Quantum Chemical Electrostatic Potential and Interaction Energies. J. Phys. Chem. B 2016, 120, 7714–7723. (13) Aburi, M.; Smith, P. E. A Combined Simulation and Kirkwood-Buff Approach to Quantify Cosolvent Effects on the Conformational Preferences of Peptides in Solution. J. Phys. Chem. B 2004, 108, 7382–7388. (14) Pierce, V.; Kang, M.; Aburi, M.; Weerasinghe, S.; Smith, P. E. Recent Applications of Kirkwood-Buff Theory to Biological Systems. Cell Biochem. Biophys. 2008, 50, 1–22. (15) Horinek, D.; Netz, R. R. Can Simulations Quantitatively Predict Peptide Transfer Free Energies to Urea Solutions? Thermodynamic Concepts and Force Field Limitations. J. Phys. Chem. A 2011, 115, 6125 – 6136. (16) Ben-Naim, A. Theoretical Aspects of Pressure and Solute Denaturation of Proteins: A Kirkwood-Buff- Theory Approach. J. Chem. Phys. 2012, 137, 235102. (17) Mukherji, D.; van der Vegt, N. F.; Kremer, K. Preferential Solvation of Triglycine in Aqueous Urea: An Open Boundary Simulation Approach. J. Chem. Theory Comput. 2012, 8, 3536–3541. (18) Ben-Naim, A. Theoretical Aspects of Self-Assembly of Proteins: A Kirkwood-Bufftheory Approach. J. Chem. Phys. 2013, 138, 224906. 33



(19) Moeser, B.; Horinek, D. Unified Description of Urea Denaturation: Backbone and Side Chains Contribute Equally in the Transfer Model. J. Phys. Chem. B 2014, 118, 107–114. (20) Kanduc, M.; Chudoba, R.; Palczynski, K.; Kim, W. K.; Roa, R.; Dzubiella, J. Selective Solute Adsorption and Partitioning Around Single PNIPAM Chains. Phys. Chem. Chem. Phys. 2017, 19, 5906–5916. (21) Nayar, D.; Folberth, A.; van der Vegt, N. F. A. Molecular Origin of Urea Driven Hydrophobic Polymer Collapse and Unfolding Depending on Side Chain Chemistry. 2017, 19, 18156–18161. (22) Ganguly, P.; Boserman, P.; van der Vegt, N. F. A.; Shea, J.-E. Trimethylamine N-oxide Counteracts Urea Denaturation by Inhibiting Protein-Urea Preferential Interaction. J. Am. Chem. Soc. 2018, 140, 483–492. (23) Lebowitz, J. L.; Percus, J. K. Long-Range Correlations in a Closed System with Applications to Nonuniform Fluids. Phys. Rev. 1961, 122, 1675–1961. (24) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford Science Publications, 2008. (25) Krüger, P.; Schnell, S. K.; Bedeaux, D.; Kjelstrup, S.; Vlugt, T. J. H.; Simon, J.-M. Kirkwood-Buff Integrals for Finite Volumes. J. Phys. Chem. Lett. 2013, 4, 4–7. (26) Cortes-Huerto, R.; Kremer, K.; Potestio, R. Communication: Kirkwood-Buff Integrals in the Thermodynamic Limit from Small-Sized Molecular Dynamics Simulations. J. Chem. Phys. 2016, 141103. (27) Salacuse, J. J.; Denton, A. R.; Egelstaff, P. A. Finite-Size Effects in Molecular Dynamics Simulations: Static Structure Factor and Compressibility. I. Theoretical Method. Phys. Rev. E 1996, 53, 2382–2389.

34


Page 34 of 37

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


(28) Dawass, N.; Krüger, P.; Schnell, S. K.; Bedeaux, D.; Kjelstrup, S.; Simon, J. M.; Vlugt, T. J. H. Finite-size Effects of Kirkwood–Buff Integrals from Molecular Simulations. Molecular Simulation 2017, doi:10.1080/08927022.2017.1416114. (29) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435–447. (30) Ben-Naim, A. Molecular Theory of Solutions; Oxford University Press, 2006. (31) Christensen, S.; Peters, G.; Hansen, F.; O’Connell, J.; Abildskov, J. Generation of Thermodynamic Data for Organic Liquid Mixtures from Molecular Simulations. Mol. Simulat. 2007, 33, 449–457. (32) Wedberg, R.; Peters, G. H.; Abildskov, J. Total Correlation Function Integrals and Isothermal Compressibilities from Molecular Simulations. Fluid Phase Equilib. 2008, 273, 1–10. (33) Schnell, S. K.; Liu, X.; Simon, J.-M.; Bardow, A.; Bedeaux, D.; Vlugt, T. J. H.; Kjelstrup, S. Calculating Thermodynamic Properties from Fluctuations at Small Scales. J. Phys. Chem. B 2011, 115, 10911–10918. (34) Schnell, S. K.; Vlugt, T. J. H.; Simon, J. H.; Bedeaux, D.; Kjelstrup, S. Thermodynamics of a small system in a µT reservoir. Chem. Phys. Lett. 2011, 504, 199–201. (35) Rosenberger, D.; Hanke, M.; van der Vegt, N. F. A. Comparison of Iterative Inverse Coarse-Graining Methods. Eur. Phys. J. Special Topics 2016, 225, 1323–1345. (36) Ben-Naim, A. Comment on “The Kirkwood-Buff Theory of Solutions and the Local Composition of Liquid Mixtures”. J. Phys. Chem. B 2008, 112, 5874–5875. (37) de Oliveira, T. E.; Netz, P. A.; Kremer, K.; Junghans, C.; Mukherji, D. C–IBI: Targeting

35



Cumulative Coordination within an Iterative Protocol to Derive Coarse-Grained Models of (Multi-Component) Complex Fluids. J. Chem. Phys. 2016, 144, 174106. (38) Reif, M. M.; Winger, M.; Oostenbrink, C. Testing of the GROMOS Force-Field Parameter Set 54A8: Structural Properties of Electrolyte Solutions, Lipid Bilayers, and Proteins. J. Chem. Theory Comput. 2013, 9, 1247–1264. (39) Ploetz, E. A.; Geerke, D. P.; Smith, P. E. To Polarize or Not to Polarize? Chargeon-Spring versus KBFF Models for Water and Methanol Bulk and Vapor - Liquid Interfacial Mixtures. J. Chem. Theory Comput. 2016, 12, 2373–2387. (40) Schnell, S. K.; Englebienne, P.; Simon, J.-M.; Krüger, P.; Balaji, S. P.; Kjelstrup, S.; Bedeaux, d.; Bardow, A.; Vlugt, T. J. H. How to Apply the Kirkwood–Buff Theory to Individual Species in Salt Solutions. Chem. Phys. Lett. 2013, 582, 154–157. (41) Liu, X.; Schnell, S. K.; Simon, J.-M.; Krüger, P.; Bedeaux, D.; Kjelstrup, S.; Bardow, A.; Vlugt, T. J. H. Diffusion Coefficients from Molecular Dynamics Simulations in Binary and Ternary Mixtures. Int. J. Thermophys. 2013, 34, 1169–1196. (42) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269–6271. (43) Weerasinghe, S.; Smith, P. E. Cavity Formation and Preferential Interactions in Urea Solutions: Dependence on Urea Aggregation. J. Chem. Phys. 2003, 118, 5901. (44) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463–1472. (45) Miyamoto, S.; Kollman, P. A. Settle: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13, 952–962. (46) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577–8592. 36


Page 36 of 37

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


(47) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684–3690. (48) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182–7190. (49) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255–268. (50) Hansen, J.-P.; McDonals, I. R. Theory of Simple Liquids; University of Cambridge, 2007. (51) Stokes, R. H. Thermodynamics of Aqueous Urea Solutions. Aust. J. Chem. 1967, 20, 10. (52) Miyawaki, O.; Saito, A.; Nakamura, K. Activity and Activity Coefficient of Water in Aqueous Solutions and Their Relationships with Solution Structure Parameters. Biosci. Biotechnol. Biochem. 1996, 61, 466–469. (53) Endo, H. The Adiabatic Compressibility of Nonelectroyle Aqueous Solutions in Relation to the Structures of Water and Solutions. Bull. Chem. Soc. Jpn. 1973, 46, 1106– 1111.

TOC Graphic

37


Convergence of KirkwoodâBuff Integrals of Ideal and Nonideal

Recommend Documents

Convergence of KirkwoodâBuff Integrals of Ideal and Nonideal