Discrete Modeling Approach as a Basis of Excess ... - ACS Publications

Subscriber access provided by READING UNIV

Article

A Discrete Modeling Approach as a Basis of Excess GibbsEnergy Models for Chemical Engineering Applications Thomas Wallek, Christoph Mayer, and Andreas Pfennig Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b04415 • Publication Date (Web): 26 Dec 2017 Downloaded from http://pubs.acs.org on January 3, 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.

Industrial & Engineering Chemistry Research 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 41 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

Industrial & Engineering Chemistry Research

A Discrete Modeling Approach as a Basis of Excess Gibbs-Energy Models for Chemical Engineering Applications Thomas Wallek,

∗,†

†

Christoph Mayer,

and Andreas Pfennig

‡

Institute of Chemical Engineering and Environmental Technology, NAWI Graz, Graz University of Technology, Ineldgasse 25/C/I, 8010 Graz, Austria, and Department of Chemical Engineering, Université de Liège - Sart-Tilman, 4000 Liège, Belgium E-mail: [email protected]

Phone: +43 (0)316 873 7966. Fax: +43 (0)316 873 7469

Abstract Discrete modeling is a concept to establish thermodynamics on Shannon entropy expressed by variables that characterize discrete states of individual molecules in terms of their interacting neighbors in a mixture. To apply this method to condensed-phase lattice uids, this paper further develops an approach proposed by Vinograd which features discrete Markov-chains for the sequential lattice construction and rigorous use of Shannon information as thermodynamic entropy, providing an in-depth discussion of the modeling concept evolved. The development comprises (1) improved accuracy compared to Monte-Carlo data and (2) an extension from a two-dimensional to a threedimensional simple lattice. The resulting model outperforms the quasichemical approximation proposed by Guggenheim, a frequently-used reference model for the simple case ∗ † ‡

To whom correspondence should be addressed Graz University of Technology, Austria Université de Liège - Sart-Tilman, Belgium

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 41

of spherical molecules with uniform energetic surface properties. To illustrate its potential as a starting point for developing g E -models in chemical engineering applications, the proposed modeling methodology is extended, using the example of a simple approach for dice-like lattice molecules with multiple interaction sites on their surfaces, to adress more realistic substances. A comparison with Monte-Carlo simulations shows the model's capability to distinguish between isomeric congurations, which is a promising basis for future g E -model development in view of activity coecients for liquid mixtures.

1. Introduction

Discrete modeling was introduced as a novel approach that uses the concept of Shannon entropy 1 to develop thermodynamic models that can describe uid-phase behavior. 24 Focusing on strongly-interacting condensed-phase systems, a preceding paper addressed the application of discrete modeling to lattice systems of nite size, which were represented by a distribution of discrete local compositions, synonymous with discrete energy classes. This rst step used the same, coarse-grained level of discretization as is used in classical statistical thermodynamics in terms of its partition functions, yet it avoids

a-priori averaging of

local compositions and limitation to systems of innite size. 4 One challenge associated with this approach is the need for appropriate degeneracy functions which subsume system states (i.e., congurations having identical energy), into one respective energy class. Generally, such functions are dicult to determine, and, particularly for systems containing nonspherical molecules even though it has been shown that Monte-Carlo simulations can be used for their determination. 4 To avoid this complication, a further renement of the discretization level from discrete energy classes of a system towards discrete states of individual molecules is investigated in this paper on the basis of simple lattices in two and three dimensions. Strongly interacting molecules positioned on the sites of a lattice are not statistically independent, and thus the molecules themselves cannot be seen as independent subsystems in the context of information theory. 2 To account for this, various approximations have been 2


Page 3 of 41 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


proposed in the previous literature. One early, simplifying viewpoint is the assumption of statistically-independent contact pairs, which is the combinatorial basis of the well-known quasichemical approach taken by Guggenheim. 57 This model produces exact results for a one-dimensional lattice and approximate ones for more dimensions due to the approximate character of the degeneracy function used. It is important to note that many actual excess Gibbs-energy models like GEQUAC, 8 MOQUAC, 9 COSMO-RS, 1014 and the pair-sampling method PAC-MAC 1518 can be reduced to the Guggenheim approach for mixtures of spherical molecules of equal size. However, the latter shows considerable deviations from Monte-Carlo data especially in the case of strong eective repulsion near the miscibility gap. To improve accuracy, cluster variation methods (CVM) were developed, 1921 mainly intended to compute phase diagrams for alloys 22 and focus on two-dimensional and threedimensional, face-centered, cubic systems. Such approaches show considerable improvements over pair approximations and are well-developed for moderate cluster sizes. 23 However, with increasing cluster size, CVM models tend to suer disadvantages due to their immense number of variables and considerable complexity, which require special algorithms for constrained minimization of the Helmholtz free energy over the course of equilibration. 24,25 Most likely due to the latter reasons, CVM methods were not adopted for the development of activity coecent models in chemical engineering applications. To overcome such limitations, Vinograd et al. proposed a dierent approach that is based on the hypothetical, sequential construction of a lattice using the Markov-chain theory. 26,27 This model was intended to describe two-dimensional crystal structures and can be seen as an approximation of rigorous CVM methods. In this approach, the entropy of the system is calculated as the Shannon information of a cluster of sites that occurs over the course of the sequential construction process of a lattice, which fully agrees with the intention of the discrete modeling concept mentioned above. Consequently, it is the goal of this paper to adopt and further develop the Vinograd 3



approach, which was developed for solid solutions in two-dimensional lattices, for threedimensional solutions, in view of providing a basis for rening g E -models in chemical engineering applications. For this purpose, rst, a renement of the Vinograd method is developed in two dimensions. Second, the Vinograd method and the rened approach are transferred to three dimensions, resulting in three-dimensional models for spherical molecules with uniform energetic surface properties. Although it is of little practical signicance, this simple case represents an important reference model for all state-of-the art g E -models and activity coecients, respectively. Third, the simplest possible extension from spherical molecules to dice-like molecules with multiple interaction sites on their surfaces is derived, to briey illustrate how the proposed modeling methodology can be extended towards more realistic substances. The results of the models developed are compared to those derived from Monte-Carlo simulations. Conclusions are drawn regarding the next steps towards developing a g E -approach that is comparable to state-of-the-art concepts and applicable to real experimental data.

2. Lattice Modeling Using Discrete Markov Chains

Lattice modeling, as introduced by Vinograd, 26 is based on the sequential construction of a lattice in terms of discrete Markov-chains, using conditional probabilities to account for the addition of a new site under the condition that one or more nearest-neighbors of the new site are already present due to previous insertions. In the following, the basic statistical laws needed to characterize this construction process and formulate internal energy and entropy of the lattice system are outlined.

Site probabilities. Using non-italicized capital letters to specify the sites themselves and italicized lowercase letters to specify the types of sites, corresponding to molecule types, the

4


Page 4 of 41

Page 5 of 41 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


probability that site A is type a is denoted

pa = Pr[A = a] .

(1)

Cluster probabilities. The combined probability of a group of sites, known as cluster, is expressed by linking the types together in a chain with centered dots. An example of a three-site cluster is

pa·b·c = Pr[A = a ∧ B = b ∧ C = c] ,

(2)

which is illustrated in Figure 1 for a one-dimensional lattice. Both sides of (2) represent the probability that site A is type a and site B is type b and site C is type c. A

B

C

Figure 1: A three-site cluster, ABC, in a one-dimensional lattice.

Conditional probabilities. In order to describe the nearest-neighbor interactions, conditional probabilities are required. The probability that site A is type a, given that site B is type b, is dened as 28

pa|b = Pr[A = a|B = b] =

Pr[A = a ∧ B = b] pa·b = . Pr[B = b] pb

(3)

Various combinations of cluster probabilities and conditional probabilities occur over the course of modeling. One example is that site A is dependent on sites B and C:

pa|b·c = Pr[A = a|B = b ∧ C = c] =

Pr[A = a ∧ B = b ∧ C = c] pa·b·c = . Pr[B = b ∧ C = c] pb·c

5


(4)


Page 6 of 41

The insertion of a further conditional probability can be written in a similar way as the denition of the conditional probability: (5)

pa·b|c = pa|b·c · pb|c .

Eq (5) shows that the probability that the a · b cluster is dependent on c can be calculated from the probability that a is dependent on the b · c cluster times the probability that b is dependent on c.

The law of total probability. Applied to cluster probabilities, this law states that the sum over all possible states of a site yields the probability that the cluster lacks that site. Using the example of a simple two-site cluster, this means that the sum over all states of b of the a · b cluster results in the probability of a:

pa =

X

pa|b · pb =

X

(6)

pa·b .

b

b

Again, a dependency on c can be considered, yielding 29

pa|c =

X

pa|b·c · pb|c =

X

(7)

pa·b|c .

b

b

Discrete Markov-chains. Taking the one-dimensional case as an example, we assume that a chain lattice such as the one illustrated in Figure 2 contains a suciently large number of molecules such that the boundary eects are negligible. The addition of site A to the end of the chain represents a step in the construction of the lattice. Since only nearestneighbor interactions are considered, the insertion process of site A is only dependent on the previously-placed site B, so that the insertion procedure represents a discrete Markov-chain. B

A

Figure 2: Addition of site A to the end of the one-dimensional lattice. 6


Page 7 of 41 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 probability that site A is type a, dependent on site B being of type b, pa|b , forms the transition matrix of the discrete Markov process, using the conditional probabilities in terms of transition probabilities from one to the next state. Considering a two-component system that comprises molecules of types 1 and 2, this transition matrix is

1

2

1 p1|1 p2|1 .

(8)

2 p1|2 p2|2 The application of the law of total probability to the rows of the transition matrix results in the expressions

p1|1 + p2|1 = 1 and

(9)

p1|2 + p2|2 = 1. This Markov-chain is ergodic except for one special case: p1|1 = 1, p2|2 = 0. This means that the probabilities that a certain site is type 1 or 2, p1 or p2 , take on constant values after an innite number of steps. These values equal those of the global composition of the mixture:

p 1 = x1

and

(10)

p 2 = x2 . The law of total probability can also be applied to a column of the transition matrix:

p1 p1|1 + p2 p1|2 = p1 .

(11)

Replacing p1|1 with expression (9) yields

p1 p2|1 = p2 p1|2

(12)

which can be interpreted as a condition of symmetry. Using cluster probabilities, eq (12)

7



Page 8 of 41

can be rewritten as

p1·2 = p2·1 .

(13)

Eqs (9), (11) and (12) will be used later during the course of the minimization of the Helmholtz free energy, representing general mathematical rather than physical constraints, in the following form:

p1|1 = 1 −

p1|2 p2 p1

p2|2 = 1 − p1|2 p2|1 =

(14)

p1|2 p2 . p1

With the three types of probabilities dened above, explicit expressions for the internal energy and entropy of the lattice system as functions of these probabilities can be formulated. The latter are the basis for modeling the Markov-chain process for the lattice. As shown by Vinograd, 26 the values of these probabilities for a given system in equilibrium are obtained by applying the extremum principle (i.e., constrained maximization of entropy or, alternatively, minimization of the internal energy or the Helmholtz free energy, respectively). As an illustration, this is exemplied in the Supporting Information for the one-dimensional lattice, yielding explicit expressions for the conditional probabilities in equilibrium and obtaining results identical to those produced by Guggenheim's quasichemical approach. 30 In the following, major aspects of the Vinograd model for two-dimensional lattices are outlined insofar as necessary to clarify the further development of this approach proposed in this paper.

3. Simple Square Lattice

Simulating a Markov-chain process, the two-dimensional square lattice is built up by the sequential construction procedure illustrated in Figure 3, assuming that the lattice is large enough that boundary eects are negligible. In this context, a decisive aspect of the Markov-

8


Page 9 of 41 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


chain is the fact that the actual insertion of a site only depends on its nearest-neighbors.

D

C

→

B

D

C

B

A

Figure 3: Sequential construction of the square lattice by addition of site A to the lattice. Each row of the lattice is lled with molecules from left to right. After one row has been completed, the next row below it is begun. This one is again lled from left to right. Because of the exclusive consideration of nearest-neighbor interactions and neglect of next-to-nearest neighbors, the addition of site A is only dependent on sites B and C. It is worth mentioning that the insertion of site B is independent of the insertion of site C: All site additions only depend on the sites directly above and to the left of them, and, consequently, sites B and C are both dependent on site D but not on one another; in fact, the additions of sites B and C are separated by an innite number of steps. With this construction in mind, the Shannon entropy of the two-dimensional square lattice can be modeled as

S = −N kB

X

pb·c

X

pa|b·c ln(pa|b·c ) .

(15)

a

b,c

In eq (15), pb·c represents the probability of the neighborhood into which site A is placed, and pa|b·c is the probability that a is inserted into the b · c neighborhood. After this insertion, the resulting a · b · c cluster represents a statistically-independent subsystem 2 of the lattice. Consequently, the total entropy of the system is then given by summations over all a, b and

c, respectively. The internal energy for nearest-neighbor interactions is calculated by adding up the 9



Page 10 of 41

interaction energies of all molecule pairs in the lattice. In terms of probabilities, this is equal to the sum of all interaction energies weighted by each of the respective probabilities of occurrence

U =N

X

pb·c

X

pa|b·c (εab + εac ) ,

(16)

a

b,c

where εij represents the interaction energies between the molecules of types i and j , which ∗ , as follows: are commonly subsumed under the term ìnterchange energy', 7 ω12

ωij∗ ≡ (εij + εji − εii − εjj ) .

(17)

In eq (17), εii and εjj designate the interaction energies between molecules of the same type, and εij = εji between molecules of dierent types. Hence, for example, a two-component system is dened by specifying ε11 , ε22 , and ε12 = ε21 . Eqs (15) and (16), which are based on conditional probabilities, can be rewritten equivalently in terms of cluster probabilities by applying the stochastic rules summarized in the previous section. The neighborhood probability multiplied by the probability of the insertion of site A yields the probability of the ABC cluster which contains both contact pairs contributing to the internal energy, A-B and A-C:

pa·b·c = pb·c pa|b·c .

(18)

Using (18), the internal energy (16) can be reformulated in terms of cluster probabilities:

U =N

X

pa·b·c (εab + εac ) .

(19)

a,b,c

The same can be achieved for the entropy function (15). In addition to eq (18) this also requires the calculation of the neighborhood probabilities from the cluster probabilities. This relationship can be found by applying the law of total probability, where the summation of

10


Page 11 of 41 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 cluster probabilities over all possible states of a yields the neighborhood probabilities:

pb·c =

X

(20)

pa·b·c .

a

Applying eqs (18) and (20) to (15) results in the entropy as a function of the probabilities of the three-site a · b · c cluster:

S = −N kB

X

pa·b·c ln

a,b,c

p P a·b·c a pa·b·c

(21)

.

The cluster probabilities of the system in equilibrium are found by constrained minimization of the Helmholtz free energy, A = U − T S , which is formulated using eqs (19) and (21):

X A = pa·b·c kB N T a,b,c

εac εab + + ln kB T kB T

p P a·b·c a pa·b·c

(22)

Eq (22) constitutes the target function for the minimization. To account for the necessary constraints, Vinograd proposes two variants, as outlined in the following discussion.

Vinograd's rst approximation. The rst approximation proposed by Vinograd treats the cluster formed by the sites B, C, and D as a one-dimensional chain, which is illustrated by the rst row in Figure 4. Site D is the starting point. Subsequently, sites B and C are each placed in a manner dependent on site D. They are, as explained before, independent from one another. Therefore, the starting point of this chain, as well as its construction direction, are arbitrary:

pb·c·d = pd pc|d pb|d = pb pd|b pc|d = pc pd|c pb|d .

(23)

The probability of the b · c · d cluster can be used to calculate the probability of the b · c cluster via summation over all possible paths that lead to this cluster. In the case of a two-component system, this implies that the summation occurs over two possible states of 11



D

D

C

→

D

C

→

D

C

→

D

C

B

A

D

C

B

A

→ B

C

Page 12 of 41

D

C

→

→ A

Figure 4: The sequential construction process of the square lattice. The rst row illustrates the construction in Vinograd's rst approximation; the rst and the second row together illustrate the two relevant directions for the construction in Vinograd's second approximation.

d. The resulting cluster represents the neighborhood in which the molecule a is to be placed:

pb·c =

X

pb·c·d =

X

d

pd pc|d pb|d .

(24)

d

Together with eq (18), which accounts for the insertion of molecule a into this neighborhood, this yields

pa|b·c = pa·b·c

X

(25)

pb·c·d .

d

In addition, making use of the fact that the sum over all options for the addition of molecule

a must equal unity, the following relation holds: X a

pa|b·c = 1 =

X

pa·b·c

X

a

pb·c·d .

(26)

d

From eq (26), Vinograd deduces in his rst approximation that the probability of the a · b · c cluster can be expressed by

pa·b·c = pa pb|a pc|a .

(27)

In a two-component system, where the indices a, b, and c each take on the states 1 and 2, the relation (27) comprises 23 = 8 equations. Together with the general relations, eq (14), these

12


Page 13 of 41 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


consitute the set of constraints for the minimization of the Helmholtz free energy, eq (22), which are necessary to determine the cluster probabilities for a given system in equilibrium. One drawback of the rst approximation, as Vinograd mentions, is that the sequential algorithm of lattice construction causes an asymmetry in next-to-nearest neighbor pairs. This means that the direction of lattice growth has an inuence on the neighborhood probability and, therefore, on the whole a · b · c · d cluster. This is illustrated in more detail in the Supporting Information.

Vinograd's second approximation. Vinograd proposes a second approximation to remedy the asymmetry of the rst approximation. The basic principle here is to keep the structure of the rst approximation while making the cluster probabilities independent of the lattice growth direction. The proposed method to achieve this goal is to formulate the probability of the square cluster, pa·b·c·d , as the arithmetic mean over all possible lattice growth directions:

pa·b·c·d =

1 pb·c·d pa|b·c + pa·c·d pb|a·d + pa·b·d pc|a·d + pa·b·c pd|b·c . 4

(28)

Due to the structure of the underlying equations, the dierent construction directions show a diagonal symmetry, so that it is only necessary to consider the two relevant directions

pa·b·c·d =

1 pb·c·d pa|b·c + pa·c·d pb|a·d . 2

(29)

These are illustrated in Figure 4. The neighborhood probabilities of the rst direction are constructed by starting with site D. Then site C is inserted in a manner which is dependent on site D. Next, site B is inserted depending on site D but not site C because they are not nearest-neighbors and, considering the whole lattice, site C is placed an innite number of steps before site B:

pb·c·d = pd pc|d pb|d . 13


(30)


Page 14 of 41

The following insertion of site A is dependent on sites B and C. The corresponding probability, pa|b·c , can be calculated starting from a cluster probability:

pa·b·c = pa pb|a pc|a pa pb|a pc|a pa·b·c =P pa|b·c = P . a pa·b·c a pa pb|a pc|a

(31)

Inserting eqs (30) and (31) into eq (29) yields the four-site cluster probabilities expressed by conditional pair probabilities as the main constraint of Vinograd's second approximation:

pa·b·c·d

1 = 2

pa pb|a pc|a pb pa|b pd|b pd pb|d pc|d P + pc pa|c pd|c P . a pa pb|a pc|a b pb pa|b pd|b

(32)

In order to reformulate the expression for entropy, eq (15), in terms of these four-site cluster probabilities, an additional fourth site must be considered in the outer sum:

S = −N kB

X

pb·c·d

X

pa|b·c·d ln(pa|b·c·d ) .

(33)

a

b,c,d

Using the denition of conditional probabilities, the conditional insertion probability of site

a is given by pa|b·c·d =

pa·b·c·d , pb·c·d

(34)

where the neighborhood probabilities, pb·c·d , are obtained from the four-site cluster probabilities via the law of total probability:

pb·c·d =

X

pa·b·c·d .

a

14


(35)

Page 15 of 41 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


Substituting the conditional probabilities with eq (34) and the neighborhood probabilities with eq (35), the entropy with regard to the four-site cluster probabilities results in

S = −N kB

X

pa·b·c·d ln

b,c,d

p P a·b·c·d a pa·b·c·d

.

(36)

Starting from eq (19), the internal energy can similarly be rewritten as a function of the four-site cluster probabilities by replacing the probabilities and expanding the summation to incorporate the molecule d:

U =N

X

pa·b·c·d (εab + εac ).

(37)

a,b,c,d

It is worth mentioning that, although the new four-site cluster contains two additional contact pairs, b · d and c · d, their interaction energies are not added to the equation. The reason is that eq (37) already sums over all contact pairs in the lattice, and the addition of the new interactions would merely double the resulting value of the internal energy. Combining eqs (36) and (37) to the Helmholtz free energy yields the target function for the minimization in terms of the four-site cluster probabilities

X εac pa·b·c·d A εab . = pa·b·c·d + + ln P p kB N T k T k T a·b·c·d B B a a,b,c,d

(38)

The constraints for this minimization comprise the general relations, eq (14), and the averaged four-site cluster probabilities, eq (32), contributing 24 = 16 equations because each of the indices a, b, c, and d can take on states 1 or 2 in a two-component system.

Renement of the two-dimensional Vinograd approach. In the following section, a rened version of Vinograd's second approximation is developed. First, this modication ensures the symmetry of the cluster probabilities, analogous to Vinograd's second approximation, by constructing the lattice in terms of four-site clusters. Second, it utilizes the 15



Page 16 of 41

same formulae for internal energy and entropy as the rst approximation, eqs (19) and (21), in terms of three-site cluster probabilities. For this purpose, the downscaling from four-site cluster probabilities to three-site cluster probabilities is achieved by applying the law of total probability to the four-site cluster probabilities from the second approximation, as shown in Figure 5. D

C

C



→

d

B

B

A

A

Figure 5: Calculation of the a-neighborhood in terms of a three-site cluster, starting from four-site cluster probabilities. Thus, the probability of the a · b · c cluster is calculated via a summation of pa·b·c·d over all possible states of d:

pa·b·c =

X

pa·b·c·d .

(39)

d

The Helmholtz free energy, the target function for constrained minimization, is that of Vinograd's rst approximation, eq (22). One constraint is obtained by inserting the four-site cluster probabilities from the second approximation, eq (32), in eq (39) to express the cluster probabilities in terms of conditional pair probabilities:

pa·b·c

pa pb|a pc|a pb pa|b pd|b 1X = pd pb|d pc|d P + pc pa|c pd|c P . 2 d a pa pb|a pc|a b pb pa|b pd|b

(40)

In a two-component system, the formula (40) comprises 23 = 8 equations, where each of the indices a, b, c, and d can take on the states 1 or 2. The remaining constraints are again the three general relations, eq (14), nally resulting in the following system of 12 equations for

16


Page 17 of 41 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 rened two-dimensional approach:

X A εac pa·b·c εab ! = + + ln P pa·b·c = min. kB N T kB T kB T a pa·b·c a,b,c pa pb|a pc|a pb pa|b pd|b 1X pd pb|d pc|d P + pc pa|c pd|c P pa·b·c = 2 d a pa pb|a pc|a b pb pa|b pd|b p1|1 = 1 −

(41)

p1|2 p2 p1

p2|2 = 1 − p1|2 p2|1 =

p1|2 p2 . p1

For given dimensionless interaction energies, εij /kB T , and composition, p1 = x1 and p2 = x2 , the 13 variables of the system (41) are the Helmholtz free energy, A, the eight three-site cluster probabilities, pa·b·c , and the four conditional pair probabilities, p1|1 , p2|2 , p1|2 , and p2|1 . Although insertion of all constraints into the target function would reduce the system to one equation for one unknown (i.e., p1|2 ), an explicit solution of this rather complex formula was not found, which is why (41) was solved numerically by means of an interior-point method. 31 The proof of the thermodynamic consistency of the rened approach using a GibbsHelmholtz equation, together with an analysis of limiting cases of model behavior, is outlined in the Supporting Information. Numerical results for the rened approach, the quasichemical model, and the two Vinograd approximations are illustrated in Figure 6, showing the relative deviations of the calculated internal energy and entropy from the Monte-Carlo data. The Monte-Carlo data for the internal energy were directly obtained from simulations of a simple square lattice with an edge length of 150, using the classical Metropolis algorithm as a basis. 3234 The Monte-Carlo entropy was derived from the internal energy using a GibbsHelmholtz integration. 34 ∗ All models perform better in the region of eective attraction (ω12 /kB T < 0) than in ∗ the region of eective repulsion (ω12 /kB T > 0), and the largest deviations occur when ap∗ proaching the miscibility gap which is around ω12 /kB T = +0.7. It is evident that the rened

17



approach proposed in this work represents a considerable improvement over the existing ∗ /kB T ), models, in particular in cases of high eective repulsion energies (i.e., positive ω12

approaching the miscibility gap.

18


Page 18 of 41

Page 19 of 41

1 .0

G u g V in o V in o th is

0 .8 0 .6

g e g r g r w o

n h e im a d 's f ir s t a p p r o x im a t io n a d 's s e c o n d a p p r o x im a t io n rk

∆U [ % ]

0 .4 0 .2 0 .0 -0 .2 -0 .4 -0 .7

-0 .5

-0 .3

-0 .1

0

0 .1

0 .3

0 .5

0 .7

0 .5

0 .7

ω1 2 * / k B T 0 .3 5


0 .3 0 0 .2 5

g e g r g r w o


0 .2 0

∆S [ % ]

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 .1 5 0 .1 0 0 .0 5 0 .0 0 -0 .0 5 -0 .7

-0 .5

-0 .3

-0 .1

0

0 .1

0 .3

ω1 2 * / k B T Figure 6: Relative deviations for internal energy, ∆U , and entropy, ∆S , between calculations (calc) by four dierent models and Monte-Carlo (MC) data, ∆x = (xcalc − xMC )/xMC , for a simple square lattice system, coordination number z = 4, with 1502 sites at a global ∗ composition of x1 = 0.3. ω12 /kB T is the dimensionless interchange energy. 19



Page 20 of 41

4. Simple Cubic Lattice

Just like the square lattice, the three-dimensional lattice is also constructed in a sequential manner. The system is built up layer by layer, whereby each layer is constructed like the square lattice, and the number of nearest-neighbors for every new site increases from two to three. Figure 7 depicts the insertion of the new site A into the lattice. The three construction directions are from back to front, top to bottom, and left to right. The addition of site A is dependent on its nearest-neighbors B, C, and D. The sites B, C, and D have no direct dependencies on one another. However, each pair of them has one common neighbor besides A. These neighbors are the sites E, F, and G, which are all dependent on site H.

H G

E C

F

H

→

G

B

E C

F

D

D

B A

Figure 7: The sequential construction of a simple cubic cluster by adding site A to the lattice. With this construction in mind, the Shannon entropy of the simple cubic lattice can be written as

S = −N kB

X b,c,d

pb·c·d

X

pa|b·c·d ln(pa|b·c·d ) ,

(42)

a

where pb·c·d represents the probability of the neighborhood into which site A is placed, and

pa|b·c·d is the probability that a is inserted into the b · c · d neighborhood. After this insertion, the resulting a · b · c · d cluster represents a statistically-independent subsystem 2 of the lattice. Consequently, the total entropy of the system is calculated by summations over all respective

a, b, c, and d. 20


Page 21 of 41 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


By applying the denition of conditional probabilities to the probability of this insertion,

pa|b·c·d , it can be expressed as the a · b · c · d cluster probability divided by the neighborhood probability:

pa|b·c·d =

pa·b·c·d . pb·c·d

(43)

The neighborhood probabilities can be calculated from the cluster probabilities via a summation over all possible states of molecule a:

pb·c·d =

X

(44)

pa·b·c·d .

a

Inserting eqs (43) and (44), eq (42) can be rewritten in terms of cluster probabilities instead of conditional and neighborhood probabilities:

S = −N kB

X

pa·b·c·d ln

a,b,c,d


.

(45)

As in the two-dimensional case, the internal energy is modeled as the sum of all contact pairs weighted by each respective probability of occurrence

U =N

X

pb·c·d

X

pa|b·c·d (εab + εac + εad ) ,

(46)

a

b,c,d

which can also be expressed in terms of cluster probabilities by the use of eq (43):

U =N

X

pa·b·c·d (εab + εac + εad ) .

(47)

a,b,c,d

All three-dimensional models are solved by numerical minimization of the Helmholtz free energy while considering the dependencies between the pair probabilities, eq (14), as general constraints. To approximate the cluster probabilities with pair probabilities, the threedimensional analogs of the two Vinograd approximations and the rened approach will be developed in the following section. 21



Page 22 of 41

The three-dimensional analog of Vinograd's rst approximation. The starting point from which the rst approximation to three dimensions is extended is eq (27) for the square lattice, which denes the cluster probability, pa·b·c , of the new site a and its two nearest-neighbors, b and c. For the three-dimensional case, this can be expanded to include all three neighbors, b, c, and d, as shown in Figure 7:

pa·b·c·d = pa pb|a pc|a pd|a .

(48)

In a two-component system, eq (48) connects the 24 = 16 cluster probabilities, pa·b·c·d , with the four conditional pair probabilities, pa|b , and represents the simplest possible approximation for the three-dimensional lattice, however, it suers from the disadvantage of an asymmetry that is similar to that of the rst approximation for the square lattice. The eect of this asymmetry is even more signicant for the cubic lattice than for the square lattice. The system in equilibrium is calculated by minimizing the Helmholtz free energy. Using the expressions for the internal energy, (47), and the entropy, (45), yields the target function for the minimization

X εac εad pa·b·c·d εab A = pa·b·c·d + + + ln P . kB N T k k k p BT BT BT a·b·c·d a a,b,c,d

(49)

The constraints of this optimization are the 24 = 16 eqs (48) and the 3 correlations between the conditional pair probabilities, eq (14). The resulting system of equations consists of 20 equations and 21 variables, which are the Helmholtz free energy, A, the 4 conditional pair probabilities, pi|j , and the 16 cluster probabilities, pa·b·c·d . It can be solved numerically as described above.

The three-dimensional analog of Vinograd's second approximation. The eight-site cluster shown in Figure 7, comprising sites A to H, is considered to transform Vinograd's 22


Page 23 of 41 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


second approximation into three dimensions. This cubic cluster will be approximated by use of pair probabilities, starting from a symmetric square cluster which has already been developed in eqs (28)-(32), and accounts for the two relevant construction directions illustrated in Figure 4. In the context of the three-dimensional cluster, the square cluster is marked with a superscript

psq a·b·c·d =

1 pd pc|d pb|d pa|b·c + pc pd|c pa|c pb|a·d . 2

(50)

The next step in the construction process is the addition of two sites to the square cluster, whereby the two newly inserted sites have two nearest-neighbors as illustrated in Figure 8. E

H

C

G

→ F D

→

B

F

A

E C

D

B A

F D

B A

Figure 8: The construction of the cubic cluster starting from a square cluster. To account for the rst step of Figure 8, the according conditional probability of adding two new sites, c · e, depending on two given sites of the square cluster, a · b, can be calculated by applying the law of total probability to this cluster

pc·e|a·b

psq psq c·e·a·b = P c·e·a·b = . sq pa·b c,e pc·e·a·b

(51)

In the second step shown in Figure 8, the resulting L-shaped cluster can be used to calculate the probability of g · h depending on c, e, d, and f

pa·b|c·e·d·f

psq a·b·d·f pc·e|a·b . =P sq a,b pa·b·d·f pc·e|a·b

(52)

This procedure has to be averaged over the dierent construction directions. Again, because of the structure of the equations, not all construction directions must be considered. Figure 9 shows the six relevant directions which suce to generate a symmetric cluster, resulting 23



Page 24 of 41

in the nal expression for the eight-site cluster probability

1 pc·e|a·b pg·h|c·e·d·f pa·b·c·d·e·f ·g·h = (psq 6 a·b·d·f +psq a·c·b·e pd·g|a·c pf ·h|d·g·b·e +psq a·d·b·f pc·g|a·d pe·h|b·f ·c·g +psq c·e·a·b

(53)

pg·h|c·e pd·f |g·h·a·b

+psq b·e·a·c pf ·h|b·e pd·g|f ·h·a·c +psq c·g·a·d pe·h|c·g pb·f |e·h·a·d ) .

One might think that, as an alternative, it would be possible to construct the cubic cluster in a similar way to that of the square cluster by starting in one corner, sequentially adding sites until the opposite corner is reached, and averaging the eight construction directions starting at each corner. This method would, however, still return an asymmetric cluster. The entropy equation is formulated using the conditional probability of inserting a into the b, c, d, e, f , g , and h neighborhood

S = −kB N

X

pb·c·d·e·f ·g·h

X

pa|b·c·d·e·f ·g·h ln pa|b·c·d·e·f ·g·h .

(54)

a

b,c,d,e,f,g,h

This conditional probability can be expressed as a ratio of the cluster probability to the neighborhood probability

pa|b·c·d·e·f ·g·h =

pa·b·c·d·e·f ·g·h . pb·c·d·e·f ·g·h

(55)

Then again, the neighborhood probability can be expressed by summation of the cluster probabilities over all possibilities of a

pb·c·d·e·f ·g·h =

X

pa·b·c·d·e·f ·g·h .

(56)

a

Finally, by combining eqs (54), (55), and (56), the entropy equation can be rewritten in 24


Page 25 of 41 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


terms of the eight-site cluster probabilities

X

S = −kB N

pa·b·c·d·e·f ·g·h ln

a,b,c,d,e,f,g,h

p P a·b·c·d·e·f ·g·h a pa·b·c·d·e·f ·g·h

.

(57)

The internal energy is calculated as a weighted sum of the interaction energies of the three new pairs that are formed during the insertion of molecule a

U =N

X

pb·c·d·e·f ·g·h

X

pa|b·c·d·e·f ·g·h (εab + εac + εad ) .

(58)

a

b,c,d,e,f,g,h

As in the case of entropy, the internal energy can be expressed using the eight-site cluster probabilities by combining eqs (55) and (58)

U =N

X

pa·b·c·d·e·f ·g·h (εab + εac + εad ) .

(59)

a,b,c,d,e,f,g,h

The equilibrium probabilities are calculated by minimizing the Helmholtz free energy, using the expressions for the entropy, (57), and the internal energy, (59)

X εab A pa·b·c·d·e·f ·g·h εac εad pa·b·c·d·e·f ·g·h = + + + ln P . kB N T k k k BT BT BT a pa·b·c·d·e·f ·g·h a,b,c,d,e,f,g,h

(60)

Just as in eq (48) in Vinograd's rst approximation, the cluster probabilities shall next be expressed in terms of conditional pair probabilties. For this purpose, the cubic cluster

25



Page 26 of 41

probabilties will rst be expressed by the square cluster probabilties

 pa·b·c·d·e·f ·g·h

1 =  P 6

sq sq sq psq a·d·b·f pc·g·a·d pb·f ·e·h pe·h·c·g P P psq psq sq sq b·f ·e·h e·h·c·g P c,g pc·g·a·d b,f pb·f ·e·h e,h psq b,f

+

+

+

+

b·f ·e·h

sq sq sq psq b·f ·a·d pc·g·a·d pe·h·b·f pe·h·c·g P psq psq P P sq sq b·f ·a·d e·h·b·f P sq b,f e,h pe·h·c·g e,h pe·h·b·f e,h pe·h·b·f sq sq sq psq c·e·a·b pa·b·d·f pc·e·g·h pg·h·d·f P P P psq psq sq sq c·e·g·h g·h·d·f P sq c,e pc·e·a·b c,e pc·e·g·h g,h c,e pc·e·g·h sq sq sq psq a·c·b·e pd·g·a·c pf ·h·b·e pd·g·f ·h P P P psq psq sq sq f ·h·b·e d·g·f ·h P sq d,g pd·g·a·c d,g pd·g·f ·h f,h d,g pd·g·f ·h sq sq sq sq pb·e·a·c pd·g·a·c pf ·h·b·e pf ·h·d·g P psq psq P P sq sq d·g·a·c f ·h·d·g P p p sq d,g f,h f ·h·d·g f,h f ·h·b·e f,h pf ·h·d·g

+ P

sq sq sq psq c·e·a·b pd·f ·a·b pg·h·c·e pg·h·d·f sq g,h pg·h·c·e

P

sq g,h pg·h·d·f

P

sq psq d·f ·a·b pg·h·d·f P sq d,f g,h pg·h·d·f

(61)

   .

Second, the use of eq (50) replaces the square cluster probabilities by the desired pair probabilities

psq a·b·c·d

1 = 2

pa pb|a pc|a pb pa|b pd|b + pc pa|c pd|c P . pd pb|d pc|d P a pa pb|a pc|a b pb pa|b pd|b

(62)

In a two-component system, the 28 = 256 eqs (61) together with the 24 = 16 eqs (62) and the 3 correlations between the conditional pair probabilities, eq (14), constitute the constraints of the target function, eq (60), to be minimized. Consequently, the resulting system of equations consists of 276 equations and 277 variables which are the Helmholtz free energy,

A, the 256 cluster probabilities, pa·b·c·d·e·f ·g·h , the 16 square cluster probabilities, psq a·b·c·d , and the 4 conditional pair probabilities, pi|j . It can be solved numerically as described previously.

26


Page 27 of 41 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


E

H

C

G

→ F D

F D

B

F

A

E

D E

C

G

B A

H

C

G

→

E C

→

B

B

A

C

→

B A

E

D

F

A

D

B A

H G

C

G

→ F D

F D

E

B

G

D E

E C

→ B

A

F

A E

H

C

D E

A

G

→

B

H

C

E C

→

B

F

A

B

F

A H

C

A

G

→

B

H

C

B

G

F

A H

C

C

→

B A

E

G

D E

A H

C

→

B

G

E C

→ F

D

A

D

A

D

B A

Figure 9: The six relevant construction directions to generate a symmetric cubic cluster.

27



Page 28 of 41

The three-dimensional analog of the rened approach. As in the two-dimensional case, the renement in three dimensions comprises two steps: First, it uses the eight-site cubic cluster previously introduced to ensure the symmetry of the cluster probabilities. Second, it utilizes the same formulae for entropy and internal energy as the rst approximation, eqs (45) and (47), to emulate the insertion of a into its direct neighborhood. Downscaling from the cubic cluster to a smaller cluster is achieved by applying the law of total probability via summing up the probabilities of all the possible states of e, f , g , and h

pa·b·c·d =

X

(63)

pa·b·c·d·e·f ·g·h ,

e,f,g,h

which is illustrated in Figure 10. Using the expressions for the internal energy and the H G

E

C

C



→ F

e,f,g,h

D

B

B D

A

A

Figure 10: Calculation of the a neighborhood from the cubic cluster entropy given by eqs (47) and (45) yields the Helmholtz free energy as the target function for constrained minimization, which is analogous to eq (38) of the two-dimensional case:

X A εab εac εad pa·b·c·d = pa·b·c·d + + + ln P . kB N T kB T kB T kB T a pa·b·c·d a,b,c,d

(64)

The constraints of this optimization are the equations that connect the cluster probabilities,

pa·b·c·d , with the conditional pair probabilities, pa|b . This connection is again derived in two

28


Page 29 of 41 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


steps. First, the cluster probabilities are expressed by the square cluster probabilities

 pa·b·c·d

sq sq sq psq a·d·b·f pc·g·a·d pb·f ·e·h pe·h·c·g P P psq psq sq sq b·f ·e·h e·h·c·g P c,g pc·g·a·d b,f pb·f ·e·h e,h psq

1 X  =  P 6 e,f,g,h +

+

+

+

b,f

b·f ·e·h

sq sq sq psq b·f ·a·d pc·g·a·d pe·h·b·f pe·h·c·g P psq psq P P sq sq b·f ·a·d e·h·b·f P sq b,f e,h pe·h·c·g e,h pe·h·b·f e,h pe·h·b·f sq sq sq psq c·e·a·b pa·b·d·f pc·e·g·h pg·h·d·f P P P psq psq sq sq c·e·g·h g·h·d·f P sq c,e pc·e·a·b c,e pc·e·g·h g,h c,e pc·e·g·h sq sq sq psq a·c·b·e pd·g·a·c pf ·h·b·e pd·g·f ·h P P P psq psq sq sq f ·h·b·e d·g·f ·h P sq d,g pd·g·a·c d,g pd·g·f ·h f,h d,g pd·g·f ·h sq sq sq sq pb·e·a·c pd·g·a·c pf ·h·b·e pf ·h·d·g P psq psq P P sq sq d·g·a·c f ·h·d·g P p p sq d,g f,h f ·h·d·g f,h f ·h·b·e f,h pf ·h·d·g

+ P

sq sq sq psq c·e·a·b pd·f ·a·b pg·h·c·e pg·h·d·f sq g,h pg·h·c·e

P

sq g,h pg·h·d·f

P

sq psq d·f ·a·b pg·h·d·f P sq d,f g,h pg·h·d·f

(65)

   .

Second, eq (62) is used to replace the square cluster probabilities by the conditional pair probabilities

psq a·b·c·d

1 = 2

pa pb|a pc|a pb pa|b pd|b + pc pa|c pd|c P . pd pb|d pc|d P a pa pb|a pc|a b pb pa|b pd|b

(66)

The correlations between the conditional pair probabilities, eq (14), nally complete the set of constraints. In a two-component system, the resulting system of equations consists of the target function, eq (64), the 24 = 16 eqs (65), the 24 = 16 eqs (66), and the 3 eqs (14), equaling a total of 36 equations. The 37 variables comprise the Helmholtz free energy, A, the 16 cluster probabilities, pa·b·c·d , the 16 square cluster probabilities, psq a·b·c·d , and the 4 conditional pair probabilities, pi|j . This system can be solved numerically as described previously. The proofs of thermodynamic consistency for all three-dimensional models developed in this paper, together with the results of the analyses of limiting cases of model behavior, 29



are outlined in the Supporting Information. Numerical results for the three-dimensional approaches and Guggenheim's quasichemical model are illustrated in Figure 11, and show the relative deviations of calculated internal energy and entropy from the Monte-Carlo data. The Monte-Carlo data for the internal energy were directly obtained from simulations of a simple cubic lattice with an edge length of 50, using the classical Metropolis algorithm as a basis. 3234 The Monte-Carlo entropy was derived from the internal energy using a GibbsHelmholtz integration. 34 As in the two-dimensional case, all models perform better in the region of eective attrac∗ ∗ /kB T > 0), and the largest /kB T < 0) than in the region of eective repulsion (ω12 tion (ω12 ∗ /kB T = +0.7. deviations occur when approaching the miscibility gap which is around ω12

Again, it is evident that the rened approach proposed in this work represents a considerable improvement over those used in existing models, in particular at high eective repulsion ∗ energies (i.e., positive ω12 /kB T ), approaching the miscibility gap.

30


Page 30 of 41

Page 31 of 41

3 .0


2 .5 2 .0

g e g r g r w o


∆U [ % ]

1 .5 1 .0 0 .5 0 .0 -0 .5 -1 .0 -0 .7

-0 .5

-0 .3

-0 .1

0

0 .1

0 .3

0 .5

0 .7

0 .5

0 .7

ω1 2 * / k B T 1 .4


1 .2 1 .0

g e g r g r w o


0 .8

∆S [ % ]

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 .6 0 .4 0 .2 0 .0 -0 .7

-0 .5

-0 .3

-0 .1

0

0 .1

0 .3

ω1 2 * / k B T Figure 11: Relative deviations for internal energy, ∆U , and entropy, ∆S , between calculations (calc) by four dierent models and Monte-Carlo (MC) data, ∆x = (xcalc − xMC )/xMC , for a simple cubic lattice system, coordination number z = 6, with 503 sites at a global ∗ composition of x1 = 0.3. ω12 /kB T is the dimensionless interchange energy. 31



Extending the modeling concept towards excess Gibbs-energy models for chemical engineering applications. The case discussed hitherto, i.e., spherical molecules with uniform energetic properties, is of great signicance as a reference point for excess Gibbsenergy models, and in this context the rened approach could be used as a more precise alternative to Guggenheim's quasichemical approximation. However, as uniform spheres represent a gross simplication of real molecules, in the following, it is briey illustrated how the proposed modeling methodology can be extended towards more realistic substances with multiple interaction sites on their surfaces, targeting better g E -models for chemical engineering applications, like activity coecients for liquid mixtures. For that purpose, the molecules are simplied as six-sided dice, where each side of a die can represent an individual energetic property. We investigate two examplary cases of pure component systems, namely molecules with two dierent charge concentrations, as illustrated in Figure 12: One àngled' conguration, cf. Figure 12(a), and one `stretched' conguration, cf. Figure 12(b). Regarding the intermolecular forces, repulsion between equal sites, ε11 /R =

ε22 /R = +1250 K, attraction between dierent sites, ε12 /R = ε21 /R = −1250 K and zero interactions with neutral areas, εi0 /R = ε0i /R = 0, are assumed, mimicing a polar molecule. The reason for choosing this specic example, which was suggested by other authors, 8 is to basically investigate the model's ability to distinct between isomeric congurations, which is a particular challenge for the state-of-the-art g E models discussed in the introduction. 30,35 As a basis for extending from spheres to dice, the three-dimensional analog of Vinograd's rst approximation is chosen because it reects the simplest possible approach and can vividly be adopted for dice-like molecules by taking the following four steps: 1. The sequential lattice construction in Vinograd's rst approximation, illustrated in Figure 7, is based on the cluster probabilities pa·b·c·d , cf. eq (48). In the case of a purecomponent system of six-sided dice there are 4 · 6 = 24 possibilities to insert a die into

32


Page 32 of 41

Page 33 of 41 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


1

1 0

0

0 2

0

0

IV B

0 0

0

I

C

D

A

II III VI

2

(a)

V

(b)

(c)

(d)

Figure 12: Surface properties of molecules modeled as six-sided dice, with two charge concentrations (labeled 1, 2) and neutral areas (labeled 0) on their surfaces. Case (a) represents an àngled' conguration, and (b) a `stretched' conguration of an isomeric molecule. Figures (c) and (d) illustrate the nomenclature used in eq (69): face I of the die on site A interacts with face VI of the die on site C, face IV of the die on site A interacts with face III of the die on site B, and face V of the die on site A interacts with face II of the die on site D. the lattice, which expands the indices of the cluster probabilities according to

pa·b·c·d

∀ {a, b, c, d ∈ N|1 ≤ a, b, c, d ≤ 24} .

(67)

2. Because each of the 24 possible orientations of a die has an equal probability of occurrence, the three correlations between the conditional pair probabilities, eq (14), change to

pa = 1/24 X X pa = pa·b = pa|b pb b

X

b

(68)

pa = 1 .

a

3. To calculate the internal energy as a function of the cluster probabilities, it must be considered that, once a molecule has been inserted in a particular orientation, the remaining contact sites are xed. This results in the following modication of eq (47)

U =N

X

pa·b·c·d (εaIV bIII + εaI cVI + εaV dII ) ,

(69)

a,b,c,d

where, for example, the index aIV bIII denes that face IV of the die a interacts with 33



Page 34 of 41

face III of the die b, as illustrated in Figure 12(c)-(d). At this point, the model intrinsically introduces a kind of geometric information which is the reason why the model is expected to distinct between the two isomeric congurations illustrated in Figure 12(a)-(b). 4. Finally, the formula for the entropy, eq (45), has to be corrected by the surplus of information which is introduced by the number of distinguishable possibilities of inserting a die, using the symbol ζ

S/N kB = ln(ζ) +

X

pa·b·c·d ln

a,b,c,d


,

(70)

with ζ = 24 distinguishable orientations of the angled isomer, cf. Figure 12(a), and

ζ = 6 distinguishable orientations of the stretched isomer, cf. Figure 12(b). Using these four modications to formulate the Helmholtz free energy, the model can be solved numerically by constrained minimization as outlined before. Figure 13 compares the internal energy predicted by the resulting model to that of MonteCarlo simulations for dice-like molecules which are discribed elsewhere. 8,30,35 It is evident that the dierence between angled and stretched isomers, which becomes apparent at low temperatures, is clearly reected by the model. Although this rst approach, based on Vinograd's rst approximation, is not yet comparable to established models like GEQUAC or COSMORS, the ability to distinct between isomers is a promising basis for further developing the method towards an activity coecient model for liquids.

5. Conclusions and Outlook

This paper illustrates the transfer of an approach from the eld of theoretical physics describing solid soultions to the eld of g E -models aiming at liquid mixtures in chemical engineering applications: Namely, the further development of the Vinograd approach, which was devel34


Page 35 of 41

-0 .4 -0 .6 -0 .8

U /k B T [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


-1 .0 -1 .2 -1 .4

a n s tr a n s tr

-1 .6

g le e tc g le e tc

d , h e d , h e

M C - s im u la tio n d , M C - s im u la tio n c a lc u la te d d , c a lc u la te d

-1 .8 4 5 0

5 0 0

5 5 0

6 0 0

6 5 0

7 0 0

7 5 0

8 0 0

8 5 0

9 0 0

9 5 0

1 0 0 0

T [K ] Figure 13: Results of the three-dimensional analog of Vinograd's rst approximation, adopted for six-sided dice, compared to Monte-Carlo (MC) simulation data, comparing a system of angled isomers with a system of stretched isomers, as illustrated in Figure 12(a)(b): The dimensionless internal energy as a function of temperature. At low temperatures, the internal energy of a pure-component system of angled isomers diers from that of a pure-component system of stretched isomers, which is also predicted by the model. oped for solid solutions in two-dimensional lattices, for three-dimensional solutions in view of describing liquid mixtures. First, the resulting rened approach for spherical molecules with uniform energetic surface properties represents a considerable improvement over the quasichemical approximation, in particular at high eective repulsion energies approaching the miscibility gap, which, for example, is an important criterion for liquid-liquid separations. Consequently, the rened approach could be used as a more precise alternative to Guggenheim's quasichemical approximation which is an important reference model for all state-of-the art g E -models and activity coecients. 35



Second, a rst extension from uniform spheres to more realistic molecules, in particular, substances with multiple interaction sites on their surfaces, illustrates the intrinsic ability of the approach to distinct between isomers. This exciting feature is a consequence of the sequential lattice construction using conditional probabilities: Once a molecule has been inserted in a particular orientation, the remaining contact sites are xed, which introduces a kind of geometric information into the model. As a next step for the further development of the modeling methodology, the latter aspect suggests to adopt the proposed rened approach for dice-like molecules instead of Vinograd's rst approximation: Compared to the latter, the rened approach introduces even more geometric information by using an eight-site cluster, where three interacting surfaces are xed for each site. This aspect oers a promising outlook, particularly as compared to state-of-the-art g E -models, which merely rely on the number densities of surface fragments without revealing the spatial information that is provided by a three-dimensional cluster.

Acknowledgement The authors thank R. Tichy for inspiring discussions, S. Crockett for English language editing, and gratefully acknowledge support from NAWI Graz.

Supporting Information Available The following le is available free of charge.

• SuppInfo.pdf

The Markov-chain process for a one-dimensional lattice Asymmetry in Vinograd's rst approximation Thermodynamic consistency and model behavior in limiting cases This material is available free of charge via the Internet at http://pubs.acs.org/. 36


Page 36 of 41

Page 37 of 41 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


References (1) Shannon, C. E. A Mathematical Theory of Communication. The

Bell System Technical

Journal 1948, 27, 379423, 623656. (2) Peger, M.; Wallek, T.; Pfennig, A. Constraints of compound systems: Prerequisites for thermodynamic modeling based on Shannon entropy. Entropy 2014, 16, 29903008. (3) Peger, M.; Wallek, T.; Pfennig, A. Discrete modeling: Thermodynamics based on shannon entropy and discrete states of molecules.

Ind. Eng. Chem. Res. 2015, 54,

46434654. (4) Wallek, T.; Peger, M.; Pfennig, A. Discrete modeling of lattice systems: The concept of Shannon entropy applied to strongly interacting systems. Ind. Eng. Chem. Res. 2016,

55, 24832492. (5) Guggenheim, E. A. Statistical thermodynamics of mixtures with zero energies of mixing. Proc. Royal Soc. of London A: Math. Phys. Eng. Sci. 1944; pp 203212. (6) Guggenheim, E. A.; McGlashan, M. L. Statistical Mechanics of Regular Mixtures. Proc.

Royal Soc. London A: Math., Phys. and Eng. Sci. 1951, 206, 335353. (7) Guggenheim, E. A.

Mixtures ; Oxford at the Clarendon Press, 1952.

(8) Egner, K.; Gaube, J.; Pfennig, A. GEQUAC, an excess Gibbs energy model for simultaneous description of associating and non-associating liquid mixtures.

Ber. Bunsenges.

Phys. Chem. 1997, 101, 209218. (9) Bronneberg, R.; Pfennig, A. MOQUAC, a new expression for the excess Gibbs energy based on molecular orientations.

Fluid Phase Equil. 2013, 338, 6777.

(10) Klamt, A. Conductor-like screening model for real solvents: a new approach to the quantitative calculation of solvation phenomena. 37

J. Phys. Chem. 1995, 99, 22242235.



Page 38 of 41

(11) Klamt, A.; Jonas, V.; Bürger, T.; Lohrenz, J. C. Renement and parametrization of COSMO-RS.

J. Phys. Chem. A 1998, 102, 50745085.

(12) Klamt, A.; Eckert, F. COSMO-RS: a novel and ecient method for the a priori prediction of thermophysical data of liquids. (13) Klamt, A.

Fluid Phase Equil. 2000, 172, 4372.

Encycl. Comp. Chem.; John Wiley & Sons, Ltd, 2002.

(14) Klamt, A.; Krooshof, G. J. P.; Taylor, R. COSMOSPACE: Alternative to conventional activity-coecient models.

AIChE J. 2002, 48, 23322349.

(15) Sweere, A. J. M.; Fraaije, J. G. E. M. Force-eld based quasi-chemical method for rapid evaluation of binary phase diagrams.

J. Phys. Chem. B 2015, 119, 1420014209.

(16) Sweere, A. J. M.; Gracia, R. S.; Fraaije, J. G. E. M. Extensive accuracy test of the forceeld-based quasichemical method PAC-MAC.

J. Chem. Eng. Data 2016, 61, 3989

3997. (17) Sweere, A. J.; Fraaije, J. G. Prediction of polymer-solvent miscibility properties using the force eld based quasi-chemical method PAC-MAC.

Polymer 2016, 107, 147153.

(18) Sweere, A. J. M.; Fraaije, J. G. E. M. Accuracy test of the OPLS-AA force eld for calculating free energies of mixing and comparison with PAC-MAC.

J. Chem. Theory

Comp. 2017, 13, 19111923. (19) Kikuchi, R. A theory of cooperative phenomena.

Phys. Rev. 1951, 81, 9881003.

(20) Kikuchi, R.; Brush, S. G. Improvement of the cluster-variation method.

The Journal

of Chemical Physics 1967, 47, 195203. (21) Sanchez, J. M.; Ducastelle, F.; Gratias, D. Generalized cluster description of multicomponent systems.

Physica A: Stat. Mech. Applic. 1984, 128, 334350.

38


Page 39 of 41 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


(22) Pelizzola, A. Cluster variation method in statistical physics and probabilistic graphical models.

J. Phys. A: Math. Gen. 2005, 38, R309R339.

(23) Kikuchi, R. Natural iteration method and boundary free energy. J.

Chem. Phys. 1976,

65, 45454553. (24) Gratias, D.; Sanchez, J. M.; De Fontaine, D. Application of group theory to the calculation of the congurational entropy in the cluster variation method.

Physica A: Stat.

Mech. Applic. 1982, 113, 315337. (25) Pelizzola, A. Generalized belief propagation for the magnetization of the simple cubic Ising model.

Nucl. Phys. B 2014, 880, 7686.

(26) Vinograd, V. L.; Perchuk, L. L. Informational models for the congurational entropy of regular solid solutions: Flat lattices.

J. Phys. Chem. 1996, 100, 1597215985.

(27) Vinograd, V. L.; Perchuk, L. L. Markov's chains and the congurational thermodynamics of a solid solution.

Vestnik Moskovskogo Universiteta (ser. Geol) 1992, 3, 4458; in

Russian. (28) Kemeny, J. G.; Snell, J. L. (29) Ferschl, F.

Finite Markov Chains ; Springer-Verlag, 1976.

Markovketten ; Lecture notes in operations research and mathematical sys-

tems, Vol. 35. Springer-Verlag, 1970. (30) Pielen, G. Detaillierte molekulare Simulationen und Parameterstudie für ein ternäres Gemisch zur Weiterentwicklung des GEQUAC-Modells. Ph.D. thesis, RWTH Aachen, 2005; in German. (31) Mehrotra, S. On the implementation of a primal-dual interior point method.

SIAM

Journal on optimization 1992, 2, 575601. (32) Zapf, F. Monte-Carlo Verfahren zur Diskretisierung von Gittersystemen. Thesis, Graz University of Technology, 2015; in German. 39



(33) Maltezos, G. Simulation lokaler Zusammensetzungen. Thesis, RWTH Aachen, 1987; in German. (34) König, L. Auswertung von Monte Carlo-Simulationen zur Validierung thermodynamischer Modelle. Thesis, Graz University of Technology, 2013; in German. (35) Ehlker, G. H. Entwicklung der Gruppenbeitragsmethode GEQUAC zur thermodynamischen Beschreibung ausgeprägt nichtidealer Gemische. Ph.D. thesis, RWTH Aachen, 2001; in German.

40


Page 40 of 41

Page 41 of 41

Graphical TOC Entry

M o n te - C a r lo d a ta

3 .0

G u g g e n h e im 's q u a s ic h e m ic a l a p p r o x im a t io n th is w o r k

2 .5 2 .0 1 .5 1 .0

∆U [ % ] d e v i a t i o n f r o m

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 .5 0 .0 -0 .5 -1 .0 -0 .7

-0 .5

-0 .3

-0 .1

ω1 2 * / k B T

in te r c h a n g e e n e r g y

0

41

0 .1

0 .3

0 .5


0 .7

Discrete Modeling Approach as a Basis of Excess ... - ACS Publications

Recommend Documents