Determination of the Residual Entropy of Simple Mixtures by Monte

Article

Determination theby the Langmuir of is published American Chemical Society. residual entropy of N.W., 1155 Sixteenth Street Washington, DC 20036

Subscriber access provided by University of Florida Published by American | Smathers Libraries Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to

simple mixtures by Monte Carlo simulation

Langmuir is published by the

Eric Detmar,American Simin Chemical Yazdi Society. 1155 Sixteenth Street N.W., Nezhad, and Ulrich K. Deiters Washington, DC 20036


Langmuir, Just Accepted Manuscript • DOI: 10.1021/acs.langmuir.7b02000 • Publication Date (Web): 21 Jul 2017 Langmuir is published by the

DownloadedAmerican from http:// Chemical Society. 1155 Sixteenth Street N.W., pubs.acs.org on July 23, 2017 Washington, DC 20036


Just Accepted

“Just Accepted” manuscripts have by been Langmuir is published the peer-re American Chemical Society. online prior to technical editing, formatting for pu 1155 Sixteenth Street N.W., Society provides “Just Accepted” as a free s Washington, DC 20036 Subscriber access provided by University of Florida dissemination of scientific material as soon as p Published by American | Smathers Libraries Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to

appear in full in PDF format accompanied by an H fully peer reviewed, but should not be considered readers and citable Langmuir by the isDigital Object Identifie published by the to authors. Therefore, the Chemical “Just Accepted” Web American Society. Sixteenth StreetisN.W., in the journal. After1155 a manuscript technically e Washington, DC 20036 Subscriber access provided by University of Florida Accepted” Web sitePublished and published as an ASAP by American | Smathers Libraries Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to

changes to the manuscript text and/or graphics and ethical guidelines that apply to the journa or consequences arising from the use of inform Langmuir is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Subscriber access provided by University of Florida Published by American | Smathers Libraries Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to

Page 1 of 16

–Smr

1 2 3 4 5 6 7 8

Langmuir

(Ne + Kr) ACS Paragon Plus Environment

λ–1

Langmuir

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

Determination of the residual entropy of simple mixtures by Monte Carlo simulation

Eric Detmar, Simin Yazdi Nezhad, and Ulrich K. Deiters∗

Institute of Physical Chemistry, University of Cologne, Luxemburger Str. 116, 50939 Koln, ¨ Germany

Abstract Extensive Monte Carlo simulations were performed for (neon + krypton) mixtures for temperatures between 200 and 600 K and pressures up to 1 GPa, using Lennard-Jones potentials to describe the intermolecular interactions. The residual entropies were obtained via Widom’s insertion method as well as by an integration technique. At high pressures, the residual entropy is, to a very good approximation, a linear function of λ−1 a , the reciprocal value of the average Monte Carlo displacement parameter that gives the acceptance ratio a for translational moves. The slope of this linear function varies linearly with the mole fraction and is related to the effective collision diameters of the molecules. As the displacement parameter is available during Monte Carlo simulations of fluids, its linearity with the residual entropy can be used to compute this properties with negligible computational effort at high densities, when particle insertion methods become unreliable.

Keywords: residual entropy; Monte Carlo simulation; high pressure; fluid mixtures

1 Introduction Of the many thermophysical properties of fluids which can be obtained by computer simulation, the group of entropy-dependent properties—entropy, Gibbs and Helmholtz energies, chemical potentials—is of particular importance, as its members control phase equilibria as well as chemical reaction equilibria; furthermore they are relevant for transport phenomena. Their calculation by computer simulation, however, requires a special computational effort, as they cannot be expressed as canonical ensemble averages like the internal energy or the virial. It is therefore not surprising that several calculation methods have been developed over the previous decades. Some of them will be discussed below; for an overview, the reader is referred to modern textbooks on computer simulation [1]. In this work we revisit the work of Yoon, Scheraga and their coworkers from more than a quarter of a century ago, who published a method for the direct calculation of the entropy ∗

Corresponding author; e-mail: [email protected]

1

ACS Paragon Plus Environment

Page 2 of 16

Page 3 of 16

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

Langmuir

of fluids by means of Monte Carlo simulations [2–4]. The method involves the so-called “radial free-space distribution function”. Unfortunately, this function cannot always be computed easily. But we could recently [5] show that, in the limit of high densities, this theory predicts an approximately linear relationship between the residual entropy and the reciprocal value of λa , the mean displacement parameter of a Monte Carlo simulation at the acceptance ratio a, c1 r + c0 . (1) S ≈ −N kB λa

N is the number of particles in the simulation ensemble. For the definition of λa it is necessary to recall the principle of a typical Monte Carlo simulation of a fluid: During such a simulation, a huge number of configurations is produced by randomly displacing the molecules according to !

~rinew = ~riold + λ~urnd

(2)

or a similar recipe, where ~riold and ~rinew denote the location of molecule i before and after the displacement, respectively; ~urnd is a unit vector with random direction, and λ a properly chosen displacement parameter or step width. Usually λ is adjusted by means of a feed-back loop in such a way that the acceptance ratio for displacements matches a pre-set value a—typically 30–50%. Consequently, the mean value of this displacement parameter, λa ≡ hλi, becomes a measure of the “lacunarity” of the fluid and, ultimately, of its residual entropy. An appealing aspect of this approach is that λ must be calculated during the simulation anyway, and that λa can therefore be obtained practically without additional computational costs. Moreover, the linearity exists at high densities, where insertion schemes for the calculation of chemical potential, e.g., Widom’s method [6], are difficult to apply. The linearity could be demonstrated for hard spheres as well as particles with Mie-(n, 6) pair potentials (8 ≤ n ≤ 16), which includes the Lennard-Jones potential. Furthermore, it could be demonstrated for spherical particles as well as for particles consisting of two or three spherical sites. In this work we investigate the relation between the displacement parameter and the residual entropy for mixtures. In order to obtain a clear picture, we shall consider a simple example only, namely the binary mixture (neon + krypton). In the literature, systems of this kind are sometimes called simple mixtures, because their molecules are spherical, interact with dispersion potentials only, and do not exhibit complications like association or chemical reactivity. But “simple” neither means “easy” nor “ideal”: Neon and krypton do not mix well at low temperatures, and the system exhibits a so-called gas–gas equilibrium of the 2nd kind [7]. One of the reasons for this behavior is the size ratio of the atoms: The volume of a krypton atom (with respect to its Lennard-Jones sizes parameter) is more than twice that of a neon atom. The calculation of phase envelopes at high pressures from equations of state requires fractal mixing rules [8, 9]. Therefore the (neon + krypton) system is an interesting test case for the relation between the Monte Carlo displacement parameter and the residual entropy.

2 Monte Carlo simulations of the (neon + krypton) system 2.1 Simulation details The Monte Carlo computer simulations reported here were carried out with the program mc++ 2


Langmuir

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], using an isobaric (N pT ) ensemble. The ensemble contained N = 512 molecules in a cubic box. Periodic boundary conditions were assumed and the minimum image convention was applied; for the latter, a fast algorithm was used [11]. Test calculations with a large particle number (1000) did not give significantly different results. The simulation runs typically consisted of an equilibration phase of 1–2 × 104 cycles and a production phase of 2.5 × 105 cycles; a cycle consists of N attempts to displace a randomly chosen molecule of the ensemble, followed by one attempt to change the box size. The acceptance or rejection of moves was decided with the usual Metropolis criterion. The control parameter governing the volume fluctuations was adjusted to obtain an acceptance ratio of 50%. The displacement parameter was also adjusted to obtain an acceptance ratio of 50%. This, however, was not possible at low densities. In this case, the maximum displacement was confined to half the box size, which resulted in an acceptance ratio a > 50%. Such a run could not be used to test the linearity, but it could be still be used to determine the chemical potential by integration, as explained below. During the production phase, mc++ splits into several threads running in parallel. This permits the comparison of the results of the threads as well as ensuring their statistical independence. For most of the simulation runs of this work three threads were used. The internal energies of these threads usually differed by 0.1% or less.

2.2 The interaction potential For the intermolecular interactions a slightly modified Lennard-Jones pair potential was used,  6  σ σij 12 hs 4ǫ − rijij rij ≥ σij ij rij , (3) uij (rij ) =  hs ∞ rij < σij

with rij denoting the center–center distance between two molecules i and j, ǫij the well depth, hs is the minimal distance of approach; the particles are quasi and σij the collision diameter. σij assumed to have rigid cores. This is advantageous in Monta Carlo simulations, as it prevents numerical overflow conditions when particles are moved to close distances as well as speeds hs = 0.8σ was used. At that distance, the energy of up the simulations. For this work, σij ij the (original) Lennard-Jones potential exceeds 42ǫij . Hence the effect of the rigid core on the thermodynamic properties of the mixtures studied here is negligible. The values of the pair potential parameters are given in Table 1; they were taken from a previous publication [12]. It should be noted that additivity of the collision diameters was assumed (σ12 = 21 (σ11 + σ22 )), whereas the well depth of the cross interaction amounts to √ merely 85% of the estimate from Berthelot’s geometric-mean rule, ǫij ≈ ǫii ǫjj . This relatively large deviation is necessary to represent the phase equilibria in the range 160–220 K. The more sophisticated combining rules of Fender and Halsey [13] and Kong [14] give values between 57 and 75 K. Calculations with equations of state require similar deviations from the geometric mean [9]. A spherical potential cut-off at 4σ was used, and the usual “tail corrections” were applied. We emphasize that the objective of this work is to demonstrate a new way to calculate the residual entropy of mixtures—not to give an accurate thermodynamic representation of the (neon + krypton) system over the full range of pressure and temperature, let alone to obtain 3


Page 4 of 16

Page 5 of 16

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

Langmuir

Table 1: Pair potential parameters of the (neon + krypton) system [12].

interaction Ne–Ne Ne–Kr Kr–Kr

ǫij /kB K 34.2 64.6 167.0

σij nm 0.2860 0.3265 0.3670

such a representation from ab initio interaction potentials. The latter would require not only accurate pair potentials for neon [15,16] and for krypton [17–20], but also the neon–krypton cross pair potential [19], corrections for three-body effects, and quantum corrections. The interested reader is referred to the recent publication of Vlasiuk et al. [21] and the literature cited therein. 2.3 Calculation of chemical potentials and residual entropies (a) Calculation with Widom’s test particle method This method can be applied during a simulation run: A particle of species i is inserted into the simulation box at a random location; the (i) required energy for this step, ∆UN , is recorded, and the particle removed. This procedure is repeated frequently. According to Widom [6], the residual chemical potential of the species i is then obtained from !+ * (i) ∆UN r (4) . µi = −kB T exp − kB T For the simulations reported here, 1–2 insertion attemps for each species were made at the end of each cycle. The residual molar entropy is then obtained from 1 1 r r Sm = (Hm − Grm ) = T T

r Hm −

Nc X

xi µri

i=1

!

.

(5)

r , can of course be obtained from the simulation run directly. The residual molar enthalpy, Hm

This method works well at low and medium densities, but becomes inefficient at high densities, because then successful particle insertions become unlikely. Another property that is of interest in this context is the residual partial molar entropy, r 1 ∂S r = (Hir − µri ) . Si ≡ ∂ni nj6=i ,p,T T

(6)

Here ni denotes the amount of substance of species i. For binary mixtures, the partial molar enthalpy, Hir , can be obtained from r ∂H r dHm r r , = Hm + (1 − xi ) Hi ≡ (7) ∂ni nj6=i ,p,T dxi where the enthalpy slope can be read off a diagram like Fig. 1a. For practical purposes it is r data and then to differentiate these advisable to perform a cubic-spline interpolation of the Hm functions analytically. 4


Langmuir

2

3

(Ne + Kr)

(Ne + Kr) 2.5

0

2 HE/(kJmol –1)

400 K Hr/(kJmol –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 6 of 16

–2 200 K –4

1.5 200 K 1

–6

–8

0.5

0

0.2

0.4

0.6

0.8

1

0

400 K

0

0.2

0.4

x1

0.6

0.8

1

x1

(a) Residual molar enthalpy.

(b) Excess enthalpy.

Fig. 1: Molar enthalpy of (neon + krypton) mixtures as function of the mole fraction. Symbols: Monte Carlo simulation, this work, ▽: 200 K, △: 400 K; open symbols: 10 MPa, filled symbols: 50 MPa. The curves are spline interpolations which merely serve to guide the eye.

For completeness’ sake we also present some excess enthalpy curves (Fig. 1b), which were obtained from Fig. 1a by subtracting the pure-fluid contributions. Evidently, the excess enthalpies are positive and quite large, which corroborates the fact that this system is far from being ideal. (b) Calculation by integration over density As already explained in our previous article, the integration method of Deiters and Hoheisel [12,22] permits the calculation of the residual Gibbs energy and then the entropy by integrating the pressure over density, Z ρ Z −1 r dρ + Z − 1 , (8) Gm = RT ρ 0 where Z = p/(ρRT ) denotes the compression factor, p the pressure, and ρ = 1/Vm the molar density. At low densities the integrand converges against the 2nd virial coefficient and hence remains finite, so that the integration can be carried out without problems. The integration was carried out by computing a smoothing cubic spline function [23] using the statistical uncertainties of the integrand data and then evaluating the integral values analytically. For this work, the virial coefficients were obtained from the exact formula B2 =

Nc Nc X X

xi xj B2,ij

(9)

i=1 j=1

with B2,ij = −2π

Z

0

∞

uij (rij ) 2 − 1 rij drij . exp − kB T 5


(10)

Page 7 of 16

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

Langmuir

The integration in this equation was performed numerically with Romberg’s method. The integration method is inefficient in the sense that, in order to determine a single Grm value, one needs a series of Monte Carlo simulations at various densities. But it gives reliable results for all densities, particularly at densities where particle insertion schemes encounter difficulties. The chemical potentials of a binary mixture can be obtained from the well-known thermodynamic relation r dGr ∂G = Grm + (1 − xi ) m . µri ≡ (11) ∂ni nj6=i ,p,T dxi The mole fraction derivative of the Gibbs energy can be obtained by means of spline interpolations in a similar way as the derivative of the enthalpy in the previous section.

3 Results and discussion 3.1 Residual entropy Figs. 2 and 3 show the residual entropy of (neon + krypton) mixtures as a function of λa over the full range of mole fractions, where λa is the average value of the displacement parameter for the mixture components, λa =

Nc X

xi λa,i .

(12)

i=1

The diagrams contain states only at which the translation acceptance ratio could be maintained at a = 50%. The simulations were performed for various temperatures between 200 and 600 K and for pressures up to 1000 MPa. In all cases the linear behavior at high pressures (= larger values of λ−1 a ) is evident. Moreover, all entropy data fall onto the same line, irrespective of temperature and pressure. For moderate densities, the entropy values calculated with the particle insertion scheme have statistical uncertainties (r.m.s. errors) of 0.2–1.5% (estimated from the results of different simulation threads or simulation runs). But, as expected, the uncertainties get large (typically 10–50%) when λa falls below 0.1 nm, whereas the integration method continues to work well. The slopes of the entropy vs. λ−1 a are presented in Table 2 as well as in Fig. 4. It appears that c1 varies linearly with the mole fraction. In our previous work [5] we could use several reference equations of state or fundamental equations (Am (ρ, T )) to obtain the residual entropies according to r r ∂ Arm Um Arm Ar Sm = − =− (13) − m . RT RT RT ∂T RT RT To our knowledge, however, there exists no high-precision fundamental equation of state for Lennard-Jones mixtures. We therefore use the equation of state for pure Lennard-Jones fluids

6


Langmuir

50

50

(0.00 Ne + 1.00 Kr)

(1.00 Ne + 0.00 Kr) 40

–Sr/(J mol–1K –1)

40


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 8 of 16

30

20

10

0

30

20

10

0

10

20

30

40

50

0

0

λ–1/nm–1

10

20

30

40

50

λ–1/nm–1

(a) Pure krypton.

(b) Pure neon.

Fig. 2: The residual molar entropy of neon or krypton as a function of the reciprocal value of the mean displacement parameter λa , a = 50%. +: MC simulation results, insertion method; : integration method. The color indicates the temperature (red . . . blue = 100 . . . 600 K); —: linear fitting function. Table 2: Slope c1 of the linearity relation Eq. (1) as a function of the neon mole fraction.

x1 0.00 0.25 0.50 0.75 1.00

c1 nm 0.10860 0.10070 0.09387 0.08628 0.07885

by Johnson et al. [24, 25] and extend it to mixtures with the van der Waals 1-fluid model: ǫσ 3 =

Nc Nc X X

3 xi xj ǫij σij

i=1 j=1

3

σ =

Nc Nc X X

(14) 3 xi xj σij

.

i=1 j=1

The predictions from this equation of state are compared with the simulation results in Fig. 5. It turns out that the equation of state describes the entropy of the (neon + krypton) system 7


Page 9 of 16

50

50

(0.25 Ne + 0.75 Kr)

(0.50 Ne + 0.50 Kr) 40



40

30

20

10

0

30

20

10

0

10

20

30

40

50

0

0

10

λ–1/nm–1

20

30

40

λ–1/nm–1

(a) Mole fraction x1 = 0.25 .

(b) Mole fraction x1 = 0.50 .

50

(0.75 Ne + 0.25 Kr) 40


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

Langmuir

30

20

10

0

0

10

20

30

40

50

λ–1/nm–1

(c) Mole fraction x1 = 0.75 . Fig. 3: The residual molar entropy of (neon + krypton) mixtures as a function of the reciprocal value of the mean displacement parameter λa . +: MC simulation results, insertion method; : integration method. The color indicates the temperature (yellow . . . blue = 200 . . . 600 K); —: linear fitting function.

remarkably well. There are some slight systematic deviations particularly at high densities, but

8


50

Langmuir

0.1

0.08

c1/nm

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 10 of 16

0.04

0.02

0

0

0.2

0.4

0.6

0.8

1

x1

Fig. 4: Slope c1 of the linearity relation Eq. (1) as a function of the neon mole fraction.

for an application to mixtures—without any adjustable parameters—the agreement is excellent. In order to explain the data in Table 2, it is helpful to recall the equation of van der Waals, p=

a RT − 2 , Vm − b Vm

(15)

from which an expression for the residual entropy can be derived, r Sm = R ln(1 − ρ˜) with

ρ˜ =

b . Vm

(16)

According to this equation, the residual entropy does not depend on temperature, but only on the reduced density ρ˜. Of course, the equation of van der Waals is too simple, but it is known to capture essential features: The most relevant contribution to the residual entropy comes from the repulsive part of the equation of state; furthermore, at the same reduced density, neon and krypton should have approximately the same residual entropy. The c1 values of these two species differ because their λa values are different—naturally, the displacement should scale with the the molecular diameter. As a first approximation, one might therefore expect the ratio of the c1 values of krypton and neon to be equal to the ratio of their σ parameters. This, however, is not quite correct, since in the temperature range studied here, 200–600 K, krypton is at a reduced temperature of 1.20–3.59, whereas neon is at 5.85–17.54, so that it is necessary to account for the temperature dependence of the effective collision diameter. A 9


Page 11 of 16

40

(0.50 Ne + 0.50 Kr)

35 30 –Sr/(J mol–1K –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

Langmuir

25 20 15 10 5 0

0

10

20

30

40

50

λ–1/nm–1

Fig. 5: The residual molar entropy of an equimolar (neon + krypton) mixture as a function of the reciprocal value of the mean displacement parameter λa . : MC simulation results, integration method. The color indicates the temperature (yellow . . . blue = 200 . . . 600 K). —: equation of state + 1-fluid theory.

correlation function for this property for a truncated Lennard-Jones potential [26, Eq. (22)], based on perturbation theory, yields a collision diameter ratio range of 1.350–1.379, while the c1 ratio from Table 2 is 1.377—a fair agreement in view of the many simplifying assumptions underlying this comparison. 3.2 Residual partial entropy Fig. 6 shows the dependence of the residual molar entropy on the mole fraction. Again, large deviations from ideality are evident. From the curves in this diagram, the partial molar entropy, Sir , can be obtained in similar way as the partial molar enthalpy from Eq. (7). In view of the linearity revealed in Figs. 2 and 3 one might be tempted to expect a linear relation between Sir and λ−1 a,i , the reciprocal mean displacement for the species i. This, however, is not correct. Differentiating Eq. (1) while observing Eq. (12) yields λa,1 − λa,2 r ′ ′ 1 S1 = −N kB (c0 + x2 c0 ) + (c1 + x2 c1 ) − c1 x2 λa λ2a (17) λa,1 − λa,2 1 , S2r = −N kB (c0 − x1 c′0 ) + (c1 − x1 c′1 ) + c1 x1 λa λ2a where the prime ′ denotes differentiation with respect to x1 . This equation is evidently quadratic i −1 in λ−1 a . Fig. 7 indeed shows that the S vs. λa functions are curved, but that Eq. (17) provides excellent fits over a wide pressure range. 10


Langmuir

25

(Ne + Kr) 20

–Sr/(kJmol –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 12 of 16

200 K

15

10 –S2r 5 400 K –S1r 0

0

0.2

0.4

0.6

0.8

1

x1

Fig. 6: Residual molar entropy of a (neon + krypton) mixtures as functions of the mole fraction. Symbols: Monte Carlo simulation, this work, ▽: 200 K, △: 400 K; open symbols: 10 MPa, filled symbols: 50 MPa. The curves are spline interpolations which merely serve to guide the eye. The arrows show the construction of the partial molar entropies for an equimolar mixture at 400 K, 50 MPa.

3.3 Other considerations Our results show that, if the intermolecular interaction in a mixture can be described by LennardJones pair potentials, the residual molar entropy at high densities is a linear function of the reciprocal mean displacement parameter at fixed acceptance ratio. This naturally leads to the question whether such a linear relation exists for other pair potentials, too. It was previously observed [5] for pure fluids, however, that the linear relation is not very sensitive to the shape of the pair potential function, as long as it is short-ranged. The linearity exists not only for the Lennard-Jones potential, but also for a range of Mie-(m, 6) potentials and the hard-sphere potential; furthermore, it exists for spherical particles as well as for multi-center particles. It therefore appears likely that the linearity will also be found in mixtures of non-Lennard-Jones fluids. If, as in our example, neon is present in a mixture, one might wonder about the influence of thermodynamic quantum effects. At the high densities considered here, the dominant contribution is that of the translational quantum effect, which amounts to an increase of the internal energy that occurs when the free path lengths of the particles fall below their thermal de Broglie wavelengths. This gives rise to a repulsive contribution to the equation of state. At high densities (where the hard-sphere fluid is a reasonably good model for the structure of the fluid) and at not too low temperatures (as in this work), the effect can be approximately accounted for by using an effective, temperature-dependent molecular diameter [27, Eq. (6.15)]. 11


Page 13 of 16

50

(0.50 Ne + 0.50 Kr) 40

30 –Si r/(J mol–1K –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

Langmuir

Kr 20 Ne 10

0

–10

0

10

20

30

40

50

λ–1/nm–1

Fig. 7: The residual partial molar entropies of an equimolar (neon + krypton) mixture as a function of the reciprocal value of the mean displacement parameter λa . : MC simulation results, integration method. The color indicates the temperature (yellow . . . blue = 200 . . . 600 K). —: Eq. (17).

Consequently, the thermodynamic functions, if expressed in reduced units, approximately remain the same. As the temperature dependence of the effective diameter is rather weak (an additive correction ∝ T −1/2 ), one might expect that the linearity rule is not much disturbed. This, however, only applies to the thermodynamic conditions studied here; simulations at lower temperatures would require perturbation theory corrections or path integral techniques, both being beyond the scope of this work.

4 Conclusion We conclude that, for short-ranged pair potentials, the mean displacement parameter of Monte Carlo simulations is closely related to the residual molar entropy not only for pure fluids, but also for mixtures. There is a linear, practically temperature-independent dependence of the residual molar entropy on the reciprocal average displacement parameter (linear average of the diplacement parameters of all mixture components). The slope of this linear function is a linear function of the mole fraction. The Monte Carlo simulations simulations reported here were performed with an translational acceptance ratio of 50%. From our previous work it is known, however, that the choice r –λ−1 of a different acceptance ratio does not afftect the linearity, but merely displaces the Sm a graph. For the partial molar entropies the dependence on the reciprocal average displacement pa12


Langmuir

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 14 of 16

rameter is quadratic. The relations represent a limiting behavior at high densities. They may be useful to calculate entropies and entropie-dependent thermodynamic properties at such densities, when other computational techniques become inefficient.

Acknowledgments The authors thank the computer center of the University of Cologne (RRZK) for the permission to use the CHEOPS supercomputer (funded by Deutsche Forschungsgemeinschaft) as well as for technical support. S. Y. N. gratefully acknowledges a scholarship of the German Academic Exchange Service (DAAD).

Symbols Arm

residual molar Helmholtz energy

a

MC displacement acceptance ratio

avdW

attarction parameter of the van der Waals equation of state

bvdW

covolume, size parameter of the van der Waals equation of state

B2

2nd virial coefficient

Grm r Hm Hir

residual molar Gibbs energy residual partial molar enthalpy of species i

kB

Boltzmann’s constant

N

number of molecules, simulation ensemble size

Nc

number of components of a mixture

ni

amount of substance of species i

p

pressure

R

universal gas constant

r

distance

~r

location vector of a molecule

r Sm

residual molar entropy

Sir

residual partial molar entropy of species i

T

temperature

r Um

residual molar internal energy

residual molar enthalpy

(i)

∆UN

energy needed to insert a particle of species i into an N -particle simulation ensemble

u(r)

pair potential

~urnd

unit vector with random orientation 13


Page 15 of 16

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

Langmuir

Vm

volume

xi

mole fraction of species i

Z

compression factor, Z = pVm /(RT )

ǫ

energy parameter of the Lennard-Jones pair potential

λ

MC displacement parameter

λa

mean displacement parameter at constant acceptance ratio a, averaged over all species

λa,i

mean displacement parameter for species i at constant acceptance ratio a

µ

chemical potential

ρ

molar density, ρ = Vm−1

ρ˜

reduced density, ρ˜ = bvdW ρ

σ

size parameter of the Lennard-Jones pair potential

σ hs

inner hard core of the pair potential

References [1] Frenkel, D.; Smit, B. Understanding Molecular Simulation, 2 edition. Academic Press, London, 2002. [2] Yoon, B. J.; Scheraga, H. A. Calculation of the entropy of a fluid by a Monte-Carlo simulation based on free volume. J. Molec. Struct. 1989, 199, 33–54. [3] Yoon, B. J.; Hong, S. D.; Jhon, M. Sh.; Scheraga, H. A. Calculation of the entropy and the chemical potential of fluids and solids from the radial free-space distribution function. Chem. Phys. Lett. 1991, 181, 73–77. [4] Yoon, B. J. Calculation of entropies and chemical potentials of hard-sphere solutes solvated in hard-sphere solids using the radial free-space distribution function. Bull. Korean Chem. Soc. 1999, 20, 1209–1212. [5] Yazdi Nezhad, S.; Deiters, U. K. Estimation of the entropy of fluids with Monte Carlo computer simulation. Mol. Phys. 2017, 115, 1074–1085. [6] Widom, B. Some topics in theory of fluids. J. Chem. Phys. 1963, 39, 2808–2812. [7] Trappeniers, N. J.; Schouten, J. A. Vapour–liquid and gas–gas equilibria in simple systems. III. The system neon–krypton. Physica 1974, 73, 546–555. [8] Deiters, U. K. Density-dependent mixing rules for the calculation of fluid phase equilibria at high pressures. Fluid Phase Equilib. 1987, 33, 267–293. [9] Deiters, U. K. Correlation and prediction of high-pressure phase equilibria and related thermodynamic properties of simple fluid mixtures. In Brunner, G., editor, Supercritical Fluids as Solvents and Reaction Media, chapter 1.8, pages 185–209. Elsevier, Amsterdam, 2004. [10] Deiters, U. K. mc++ project homepage: http://www.uni-koeln.de/deiters/mc++/index.html, accessed on 2016-02-17. 14


Langmuir

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 16 of 16

[11] Deiters, U. K. Efficient coding of the minimum image convention. Z. Phys. Chem. 2013, 227, 345–352. [12] Hoheisel, C.; Deiters, U. High pressure molecular dynamics of the partially miscible fluid mixture neon/krypton. Mol. Phys. 1979, 37, 95–109. [13] Fender, B. E. F.; Halsey, jr., G. D. Second virial coefficients of argon, krypton, and argon– krypton mixtures at low temperatures. J. Chem. Phys. 1962, 36, 1881–1888. [14] Kong, C. L. Combining rules for intermolecular potential parameters. II. Rules for the Lennard-Jones (12–6) potential and the Morse potential. J. Chem. Phys. 1973, 59, 2464–2467. [15] Leonhard, K.; Deiters, U. K. Monte Carlo simulations of neon and argon using ab initio potentials. Mol. Phys. 2000, 98, 1603–1616. [16] Hellmann, R.; Bich, E.; Vogel, E. Ab initio potential energy curve for the neon atom pair and thermophysical properties of the dilute neon gas. I. Neon–neon interatomic potential and rovibrational spectra. Mol. Phys. 2008, 106, 133–140. [17] Nasrabad, A. E.; Deiters, U. K. Prediction of thermodynamic properties of krypton by Monte Carlo simulation using ab initio interaction potentials. J. Chem. Phys. 2003, 119, 947–952. [18] Slav´ıcˇ ek, P.; Kalus, R.; Paˇska, P.; Odvárková, I.; Hobza, P.; Malijevsky, ´ A. State-of-the-art correlated ab initio potential energy curves for heavy rare gas dimers: Ar2 , Kr2 and Xe2 . J. Chem. Phys. 2003, 119, 2102–2119. [19] Haley, T. P.; Cybulski, S. M. Ground state potential energy curves for He–Kr, Ne–Kr, Ar– Kr, and Kr2 : Coupled-cluster calculations and comparison with experiment. J. Chem. Phys. 2003, 119, 5487–5496. [20] Jäger, B.; Hellmann, R.; Bich, E.; Vogel, E. State-of-the-art ab initio potential energy curve for the krypton atom pair and thermophysical properties of dilute krypton gas. J. Chem. Phys. 2016, 144, 114304:1–21. [21] Vlasiuk, M.; Frascoli, F.; Sadus, R. J. Molecular simulation of the thermodynamic, structural, and vapor–liquid equilibrium properties of neon. J. Chem. Phys. 2016, 145, 104501:1– 16. [22] Ram´ırez-González, P. V.; Quinones-Cisneros, ˜ S. E.; Deiters, U. K. Thermophysical properties of the Lennard-Jones chain fluid. Mol. Phys. 2015, 113, 28–35. [23] Reinsch, C. M. Smoothing by spline functions. Numer. Math. 1967, 10, 177–183. [24] Johnson, J. K.; Zollweg, J. A.; Gubbins, K. E. The Lennard-Jones equation of state revisited. Mol. Phys. 1993, 78, 591–618. [25] Johnson, J. K.; Muller, ¨ E. A.; Gubbins, K. E. Equation of state for Lennard-Jones chains. J. Phys. Chem. 1994, 98, 6413–6419. [26] Boshkova, O. L.; Deiters, U. K. Soft repulsion and the behavior of equations of state at high pressures. Int. J. Thermophys. 2010, 31, 227–252. [27] Singh, Y.; Sinha, S. K. Semiclassical statistical mechanics of fluids. Phys. Rep. 1981, 79, 213–329.

15


Determination of the Residual Entropy of Simple Mixtures by Monte

Recommend Documents