Article pubs.acs.org/IECR
Discrete Modeling of Lattice Systems: The Concept of Shannon Entropy Applied to Strongly Interacting Systems Thomas Wallek,*,† Martin Pfleger,† and Andreas Pfennig‡ †
Institute of Chemical Engineering and Environmental Technology, NAWI Graz, Graz University of Technology, Inffeldgasse 25/C/I, 8010 Graz, Austria ‡ Department of Chemical Engineering, Université de LiègeSart-Tilman, 4000 Liège, Belgium S Supporting Information *
ABSTRACT: Discrete modeling is a novel approach that uses the concept of Shannon entropy to develop thermodynamic models that can describe fluid-phase behavior. While previous papers have focused on reviewing its theoretical background and application to the ideal-gas model as one limiting case for fluid phases, this paper addresses its application to lattice models for strongly interacting condensed phase systems, which constitute the other limiting case for fluids. The discrete modeling approach is based on the discrete energy classes of a lattice system of finite size, represented by a distribution of discrete local compositions. In this way, the model uses the same level of discretization as classical statistical thermodynamics in terms of its partition functions, yet avoids (1) a priori averaging of local compositions by utilizing a distribution, and (2) confinement to systems of infinite size. The subsequent formulation of the probabilities of discrete energy classes serves as the basis for introducing the concept of Shannon information, equivalent to thermodynamic entropy, and for deriving the equilibrium distribution of probabilities by constrained maximation of entropy. The results of the discrete model are compared to those derived from Monte Carlo simulations and by applying the Guggenheim model of chemical theory. We point out that this applicability of discrete modeling to systems of finite size suggests new possibilities for model development.
1. INTRODUCTION Chemical and process engineering is currently facing a change in resources, whereby increasing fractions of biogenic components are appearing in many lines of products. As a consequence, thermodynamic models such as equations of state and gE approaches are being challenged by the need to describe more complex and/or oxygenated molecules, which are characterized by multiple, spatially distributed interaction sites. For engineering calculations, mixtures of such complex components also pose a challenge to frequently used models such as UNIQUAC or similar state-of-the-art approaches commonly used in molecular thermodynamics.1 Such models are intended to account for strongly interacting systems, in particular liquid mixtures with miscibility gaps, and are mostly based on such elements as lattice theories, statistical thermodynamics, and the local composition concept. However, the fact that local compositions represent averaged system values limits the applicability of such approaches in accounting for detailed molecular information, which is needed to cope with the wide spectrum of complex molecules. Hence, the basic elements of current models must be thoroughly scrutinized and, where necessary, extended. From the authors’ points of view, two steps of first priority, among many others, need to be taken: first a consideration of three-dimensional, orientation-dependent molecular interactions, and second a refinement of the discretization level beyond averaged local compositions, if possible, toward the detailed (discrete) states of individual molecules. The first demand was recently incorporated into the COSMOsim3D model,2 using local σ-profiles in a spatial grid, and the MOQUAC model,3 which considers orientation-dependent © XXXX American Chemical Society
interactions in three dimensions, yet is still based on averaged local compositions. To address the second need, this paper extends the local composition concept from averaged local compositions to a distribution of discrete local compositions, which represent the energy classes inherent in the system. This is intended to be the first step of refinement taken toward discrete states of molecules and offers researchers the opportunity to draw direct comparisons with classical statistical thermodynamics. The method of discrete modeling was introduced in two previous papers as a novel approach to incorporate a more detailed picture of molecular conditions in thermodynamic models.4,5 As a proof of concept, the thermal and caloric equations of state, the heat capacity, and the Maxwell− Boltzmann distribution of energies for ideal gas were properly derived on the basis of the discretized states of individual molecules. In this paper, the discretized states of a lattice system as a whole are used as the basis for the application of discrete modeling to spherical molecules in binary condensed-phase mixtures, with the long-term goal of gE-model development. This comprises the following steps: • First is an analysis of the lattice system, considering the discrete states and the combination of states of equal internal energy in energy classes, and a formulation of the corresponding probabilities of those classes. Received: November 22, 2015 Revised: February 1, 2016 Accepted: February 2, 2016
A
DOI: 10.1021/acs.iecr.5b04430 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research • Next is the use of these probabilities as a joint set of variables to express the Shannon information on the system and all relevant constraints, such as internal energy due to molecular interactions. • The last step is the use of the Shannon information in an equivalent way to thermodynamic entropy (in this paper referred to as “Shannon entropy”) by determining the equilibrium distribution of the probabilities of the energy classes by constrained maximation of Shannon entropy. The latter represents a concrete application of the maximum entropy principle as described in a general form by previous seminal works.6−8 The discrete model is applied to a one-dimensional and simple cubic lattice, and the results are compared to those of the well-known Guggenheim model of chemical theory as well as Monte Carlo simulations. In the Conclusion, we present arguments that critically assess the new possibilities offered by discrete modeling.
molecules of types 1 and 2. For a specific contact number, N12,j, the internal energy of the lattice is then given by9 Uj =
(2)
where ε11, ε22, and ε12 represent the interaction energies between the molecules of different types, and the abbreviation ω12 * usually represents the interchange energy. It is evident from eq 2 that, for a system of given composition N1 and N2, the contact number N12,j fully specifies the internal energy Uj of the lattice and can, therefore, be used as a criterion for assessing discrete classes of internal energy. This is essentially the same level of discretization used in partition functions of classical statistical thermodynamics. 2.3. Degeneracy of States. Because various different states of the lattice may show one and the same contact number, N12,j, with respectively internal energy, Uj, it is convenient to combine all these states in an energy class j. The number of system states resulting in one and the same value for N12,j is called degeneracy, gj. Hence, a certain class j is characterized by a pair of variables for the contact number and the corresponding degeneracy: internal energy class j:
2. DISCRETE STATES OF A LATTICE SYSTEM The lattice framework used in this paper is a one-dimensional or simple cubic lattice with periodic boundary conditions and fixed sites, each of which is occupied by one spherically symmetrical molecule. The model is confined to a binary system comprising N1 molecules of type 1 and N2 molecules of type 2. Only z nearest (touching) neighbors contribute to the energetic interactions, with z representing the coordination number. 2.1. Discrete States. The discrete state of a lattice system is defined by detailed knowledge about the occupation of each single lattice site by a certain molecule type, with different states representing different arrangements. For a binary system, the number of feasible states is given by ⎛ N1 + N2 ⎞ ⎛ N1 + N2 ⎞ (N1 + N2)! ⎟=⎜ ⎟= m=⎜ N1! N2! ⎝ N1 ⎠ ⎝ N2 ⎠
1 1 1 zN1ε11 + zN2ε22 + N12, j (2ε12 − ε11 − ε22) 2 2 2 * ω12
{N12, j , gj}
∀ (1 ≤ j ≤ n)
(3)
where n is the number of energy classes, and the sum of all degeneracies equals the number of possible system states, m, given by (1): n
∑ {gj} = m
(4)
j=1
As stated in eq 3, the complete combinatorial description of energy classes of a given system (N1, N2, z) includes knowledge of all n feasible energy classes N12,j as well as their corresponding degeneracies, gj:
(1)
which equals the “thermodynamic probability” when using Boltzmann’s statistics. As illustrated in Figure 1, the number of N possible system states reaches a maximum at x1 = N +1N = 0.5, 1
N12, j = N12, j(j , N1 , N2 , z)
2
rapidly decreasing with other compositions. 2.2. Discrete Energy Classes. To characterize the internal energy of a binary lattice system, it is convenient to use the contact number, N12, which designates all contacts between
∀ (1 ≤ j ≤ n)
gj = gj(N12, j(j , N1 , N2 , z), N1 , N2 , z)
(5a) (5b)
The availability of both specifications, (5a) and (5b), is an indispensable prerequisite for discrete modeling, either in tabular form or preferably as functions. Of course, degeneracy functions (5b) are also used in conventional models that are based on classical statistical thermodynamics. But in that case they are continuous functions, and not connected with discrete N12,j using function (5a).10 2.4. Probabilities of States and Probabilities of Energy Classes. Let pi be the probability of the system existing in either of m possible states, with m given by (1) and m
∑ {pi } = 1 i=1
as a natural constraint. Considering the energy classes, one can group these probabilities using two indices: j for the energy class, and k for one specific system arrangement that represents energy class j, among many others. The system probabilities pj,k then underlie the natural constraint
Figure 1. Ratio of the number m(x1) of system states at global N composition x1 = N +1N and the number m(0.5) of system states at 1
(6)
2
global composition x1 = 0.5 for a system of N1 + N2 = 105 molecules. Apart from x1 = 0.5, the function shows a rapid decline. B
DOI: 10.1021/acs.iecr.5b04430 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research n
1=
gj
∑ ∑ {pj ,k }
3.1. The Random System. The random system is defined as having a zero interchange energy and/or infinite temperature. From the standpoint of Shannon entropy, it is characterized by the absence of any constraints to the Shannon entropy, except for the normalization condition of the probabilities, described in eq 10. As outlined in the Introduction, discrete modeling uses the constrained maximation of Shannon entropy6−8 to determine the equilibrium distribution of the probabilities of energy classes. For a random system, the following system of equations has to be solved for p12,j:
(7)
j=1 k=1
It is justifiable to assume equal probabilities for all states within one and the same energy class, which is in line with the principle of equal a priori probabilities of classical statistical thermodynamics11 and suggests that the notation p12, j pj ,1 = pj ,2 = ... = pj , g ≡ j gj (8) gj
⇒ ∑ {pj , k } = gj k=1
p12, j gj
n ⎧ ⎛ p ⎞⎫ ⎪ 12, j ⎟⎪ ! ⎬ S(p12, j ) = −kB ∑ ⎨p12, j ln⎜⎜ ⎟ = max j=1 ⎪ ⎝ g j ⎠⎪ ⎩ ⎭
= p12, j (9)
target function
would be appropriate, defining p12,j as the system’s probability of existing in energy class j. Combining eqs 7−9 yields the natural constraint for p12,j ⎧ ⎫ ⎪ p12, j ⎪ ⎬= 1 = ∑ ∑ {pj , k } = ∑ ∑ ⎨ ⎪ g ⎪ j=1 k=1 j=1 k=1 ⎩ j ⎭ n
gj
n
gj
n
f1 (p12, j ) =
∑ {p12,j } = 1
constraint (14b)
j=1
n
∑ {p12,j } j=1
(14a)
Constrained Maximation of Entropy. The starting point for this constrained maximation is the Lagrange function, 3(p12, j ),
(10)
which links the target function and the constraint by means of the Lagrangian multiplier, λ1:
In a next step, these probabilities will be used as a joint set of variables to express the Shannon entropy and all relevant constraints for both random systems and nonrandom systems that consider molecular interactions.
n
3(p12, j ) = S(p12, j ) − λ1(1 −
∑ {p12,j }) j=1
3. SHANNON ENTROPY As in eqs 6−10, we start from the Shannon entropy that is based on the probabilities of system states, pi, which can be subsumed in energy classes j and represented by a specific system arrangement k:
The Lagrangian multiplier is a constant that must be determined during the course of maximation, which is achieved by setting the first partial derivatives to zero: ∂3(p12, j )
m
S /kB = −∑ {pi ln(pi )} i=1
n
∂p12, j (11)
j=1 k=1
(12)
The assumption of equal probabilities for all states within one and the same energy class according to (8) yields an expression for Shannon entropy that is based on probabilities of the energy classes
(16)
Inserting (17) in (14a) and applying (4) and (1) allows the determination of the Shannon entropy for the random system: ∞ ⎫ ⎧ ⎪ ∞ ⎛ p12, j ⎞⎪ ⎟⎬ ⎜ ⎨ S /kB = −∑ p12, j ln⎜ ⎟ ⎪ ⎝ g j ⎠⎪ j=1 ⎩ ⎭ n
∞
gj ⎧ ⎫ n ⎪ p12, j ⎛ p12, j ⎞⎪ ⎟⎬ S / kB = − ∑ ∑ ⎨ ln⎜⎜ ⎟ ⎪ g ⎝ g j ⎠⎪ j=1 k=1 ⎩ j ⎭
⎧ g ⎛ g ⎞⎫ ⎪ j j ⎟⎪ ⎬ = −∑ ⎨ ln⎜⎜ ⎟⎪ ⎪ m mg ⎝ j ⎠⎭ j=1 ⎩ n
⎫ n ⎧ p ⎪ 12, j ⎛ p12, j ⎞⎪ ⎟⎬ = −∑ ⎨gj ln⎜⎜ ⎟ ⎪ gj ⎝ g j ⎠⎪ j=1 ⎩ ⎭ ⎛ p ⎞⎫ n ⎧ ⎪ 12, j ⎟⎪ ⎬ = −∑ ⎨p12, j ln⎜⎜ ⎟ ⎪ ⎝ g j ⎠⎪ j=1 ⎩ ⎭
∀j
=0
This calculation is presented in detail in the Supporting Information. The results are the equilibrium probabilities for the random system gj ∞ p12, = j (17) m
gj
S /kB = −∑ ∑ {pj , k ln(pj , k )}
(15)
=
ln(m) m
n
∑ {gj} j=1
⎛ (N + N2)! ⎞ = ln(m) = ln⎜ 1 ⎟ ⎝ N1! N2! ⎠
(13)
which will be used as a general entropy function throughout this paper. A more detailed argument supporting eq 13 was published in a previous paper.5 In the next step, the constraints to this function are determined for both random systems and nonrandom systems that consider intermolecular forces.
(18)
Equation 18 is equivalent to Boltzmann’s entropy approach, S/kB = ln(W) . By applying Stirling’s approximating formula in a truncated form, ln(x!) = x ln(x) − x, it is possible to determine the entropy of ideal mixing: C
DOI: 10.1021/acs.iecr.5b04430 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research S iM /kB = −N1 ln(x1) − N2 ln(x 2)
n
(19)
UE ≡ U −
On one hand, this result confirms the number of discrete system states determined using eq 1, because it yields the correct thermodynamic equivalent. On the other hand, it must be emphasized that a homogeneous function of the number of molecules is only achieved after application of Stirling’s approximation. [The complete series in Stirling’s approximating formula is 1 1 1 ln(x! ) ≈ x ln(x) − x + ln(2πx) + − + ... 2 12x 360x 3 B 2k 1 + + ... 2k − 1 (2k − 1)2k x
⎛1 ⎞ 1 * 1 ⎜ zN ε + zN ε ⎟ = ω12 ∑ {N12,jp12,j } ⎝ 2 1 11 2 2 22⎠ 2 j=1 (22)
To determine the equilibrium distribution, the following system of equations must be solved for p12,j with given values U, N1, N2 and the interaction energies and their respective interchange energy: n ⎧ ⎛ p ⎞⎫ ⎪ 12, j ⎟⎪ ! ⎬ S(p12, j ) = −kB ∑ ⎨p12, j ln⎜⎜ ⎟ = max j=1 ⎪ ⎝ g j ⎠⎪ ⎩ ⎭
target function
for large x.] Figure 2 illustrates the coincidence of eqs 18 and 19 with increasing system size. The essential conclusion is that a lattice
(23a)
n
f1 (p12, j ) =
∑ {p12,j } = 1
first constraint (23b)
j=1
n
U (p12, j ) = =U
1 1 1 * zN1ε11 + zN2ε22 + ω12 ∑ {N12,jp12,j } 2 2 2 j=1 (23c)
second constraint
Constrained Maximation of Entropy. As an analogue of eq 15, the Lagrange function has been formulated to link the target function and the constraints by means of the Lagrangian multipliers λ1 and λ2: n
3(p12, j ) = S(p12, j ) − λ1(1 −
j=1
Figure 2. Ratio of eqs 18 and 19 as a function of system size, log(N1 + N2), for various global compositions x1. For systems comprising 105 molecules or more, the two equations coincide sufficiently.
− λ 2(U − U (p12, j ))
edge length of 3 105 ≈ 50, which is the typical system size chosen in Monte Carlo simulations. 3.2. System with Energy Constraint. Again we use eq 13 as an entropy function along with the natural constraint, the normalization condition (10). As a second constraint, the internal energy U of the lattice system must be considered by applying 1 1 1 * U = zN1ε11 + zN2ε22 + ω12 N12 ̅ (20) 2 2 2 an analogue to eq 2, where N̅ 12 designates the mean contact number of the system in equilibrium. Given that N̅ 12 is a homogeneous function of N1 and N2, eq 20 represents a homogeneous constraint. As shown in a previous paper,5 this ensures the homogeneity of the Shannon entropy when the maximum entropy principle is applied. To link the constraint (20) with the distribution variables, p12,j, the mean value of the contact number is expressed in terms of the probabilities
p12, j =
gjq N12,j n
∑k = 1 {gk q N12,k }
∀j (25)
analogous to the random system (17), with the abbreviation ⎛ λ ω* ⎞ q = exp⎜ 2 12 ⎟ ⎝ kB 2 ⎠
(26)
and an equation for the determination of q for a given value of U n
n
1 1 1 * zN1ε11 + zN2ε22 + ω12 ∑ {N12,jp12,j } 2 2 2 j=1
(24)
where U denotes a concrete value for the corresponding function U(p12,j). As in the random system, cf. eq 16, the Lagrangian multipliers are determined during the course of maximation which is achieved by setting the first partial derivatives to zero. The goal is to determine the distribution variables p12,j, the two Lagrangian multipliers, λ1 and λ2, and the entropy, S. For determination of these (j + 3) variables, j equations (eq 16) and 3 equations (eqs 23a, 23b, and 23c) are available. This calculation is presented in detail in the Supporting Information. The results are the equilibrium probabilities
of finite size, at least one that comprises at least 105 sites, can be regarded as representative, and can be used to draw conclusions about real thermodynamic systems, ensuring the homogeneity of Shannon entropy with a systematic error that is negligible. Considering a cubic lattice, this corresponds to a minimum
U=
∑ {p12,j })
U=
N12, j 1 1 1 * ∑ j = 1 {N12, jgjq } zN1ε11 + zN2ε22 + ω12 n 2 2 2 ∑k = 1 {gk q N12,k }
(21)
(27)
in which the last summand represents the excess internal energy
The true meaning of λ2 in eq 26 is revealed when the definition of thermodynamic temperature is recalled D
DOI: 10.1021/acs.iecr.5b04430 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research ⎛ ∂S(U , V , N ) ⎞ 1 ⎜ ⎟ ≡ ⎝ ⎠V , N T ∂U
section 4.2 to a cubic lattice will be based on two different degeneracy functions: the first one is based on combinatorial considerations as derived by Guggenheim;9,10 the second one is obtained from Monte Carlo simulations. 4.1. Application to a One-Dimensional Lattice. Degeneracy Function. For the one-dimensional lattice, a rigorous solution based on combinatorial considerations was used as a degeneracy function:12
(28)
and the Lagrange function (24) is inserted into (28), resulting in the identity λ2 ≡ −
1 T
(29)
which was already pointed out by other authors in a more general context.6,8 Finally, by inserting (29) into (26), the abbreviation q assumes the meaning of a Boltzmann factor
N12, j = jz
gj(N1 , N2 , N12, j) =
⎛ 1 ω* ⎞ 12 q ≡ exp⎜ − ⎟ 2 k ⎝ BT ⎠
(30)
gj q
p12, j =
n
∑k = 1 {gk exp( −UkE/kBT )}
( 12 N12,j)!⎤⎦ ⎡⎣ 12 (2N2 − N12,j)⎤⎦! 2
An outline of the derivation of eq 33b is part of the Supporting Information. The contact pairs can also be expressed in terms of local compositions,9 which are formulable both as discrete values and as averaged system values x12, j ≡
N12, j zN2
x12 ̅ ≡
N12 ̅ zN2
(34)
whereby, in random systems, the mean local composition equals the global composition ∞ x12 ̅ = x1
(35)
Using these local compositions, Table 1 illustrates the number of discrete energy classes, (33a), for different systems, showing the increasing refinement of discretization with the increasing system size.
(31)
where the term UEj designates the excess internal energy of a system of energy class j; cf. eq 22. In this notation, the equilibrium probabilities (25) can be reformulated as gj exp( −UjE/kBT )
2N1N2 N1! N2!
(33b)
⎡ ⎛ ω * ⎞⎤ N12,j = gj⎢exp⎜ − 12 ⎟⎥ ⎢⎣ ⎝ 2kBT ⎠⎥⎦ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ * 1 ω12 ⎜ = gj exp − N ⎟ ⎜ 2 kBT 12, j ⎟ ⎜ ⎟ ⎜ U E/k T ⎟ ⎝ ⎠ B j
(33a)
(N1 + N2)N12, j
⎡ 1 (2N − N )⎤ ! ⎡ ⎣2 1 12, j ⎦ ⎣
As in eq 18, an expression for the Shannon entropy can be derived by inserting (30) and (25) into (23a). Interpretation from the Classical Viewpoint. Using eq 30, the determining factor gjqN12,j of the equilibrium probabilities expressed through eq 25 can be rewritten as follows: N12, j
(1 ≤ j ≤ N1 , N1 ≤ N2)
Table 1. Refinement of Discretization with Increasing System Size, in a One-Dimensional Lattice (z = 2) of Given Global Composition, x1 = 0.3a
∀j (32)
range of x12,j
The denominator of (32) is equivalent to the classical statistical-mechanical partition function of a system with discrete energy classes, and the entire expression represents the interpretation of probability as the number of desired outcomes divided by the total number of all possible outcomes of a random experiment.
a
4. APPLICATION OF THE MODEL From the standpoint of constrained maximation, the given values of U, N1, N2, and the interaction energies respectively interchange energy are needed to calculate the equilibrium probabilities and the corresponding Shannon entropy of the system. To be precise, the specification of U allows (27) to be solved for q respectively temperature, which is used for the subsequent calculation of the equilibrium distribution (25) and entropy (23a). However, in terms of practical application, it is more common to specify the temperature of a system rather than its internal energy. In this case, eq 27 can be used to calculate U at a given temperature, respectively q, and all other values are determined as stated above. In section 4.1, the model will first be applied to a onedimensional lattice because, in this case, a rigorously derived degeneracy function is available. The subsequent application in
The degeneracies, (33b), are illustrated in Figure 3, using relative frequencies of degeneracy. Although this is initially suggested by the graph, it should be noted that this function only approximates a normal distribution, and significant deviations from a normal distribution apart from its maximum are seen. According to eq 35, the maximum occurs at a discrete x12,j equal to global composition, x1. Apart from the global composition, the number of possible states is reduced significantly. The bigger the system, the more rapidly this function declines apart from the global composition. From the classical viewpoint, this fact is taken into account when the sum over all degeneracies is replaced by its maximum term for infinite systems.9,10 Resulting Equilibrium Distribution. The degeneracy functions (33b) were inserted into the results of discrete modeling, cf. section 3.2, to calculate the equilibrium probability E
N1 + N2
no. system states
min
max
no. discrete x12,j
103 104 105
5.43 × 10263 7.76 × 102650 8.70 × 1026526
0.001 428 57 0.000 142 857 0.000 014 285 7
0.428 571 0.428 571 0.428 571
300 3000 30000
The number of discrete x12,j equals the number of energy classes.
DOI: 10.1021/acs.iecr.5b04430 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research
needed. Therefore, particular attention must be paid to correctly describing the degeneracy function over the entire range of the discrete local compositions, not only around the peak of this function. Resulting Equilibrium Distribution in the Context of Boltzmann’s Law. From the previous equations, it is evident that the specification of ω12 * /kBT determines the internal energies of each energy class, (Uj/kBT), according to (2) and the equilibrium probabilities, p12,j, according to (25), for a system of given composition N1, N2, and interaction energies. To compare two equilibrium states of the lattice, we can choose two arbitrary dimensionless interchange energies, such as (ω*12/kBT)(1) and (ω*12/kBT)(2), for which the according (Uj/ (2) kBT)(1), (Uj/kBT)(2) and p(1) 12,j, p12,j can be calculated as shown above. Plotting the latter equilibrium probabilities against the discrete local compositions, cf. Figure 6, the logarithmic ratio of
Figure 3. For a one-dimensional lattice system of given composition, the relative frequencies of degeneracy show a maximum at a discrete x12,j equal to global composition. With increasing system size, the relevant contributions to this function increasingly concentrate around the peak.
distributions of the energy classes, (25), for a system of given N1 and N2 but different values of ω12 * /kBT. These are illustrated in Figure 4, showing the maximum probability at the global composition for the random system and the shift of the maximum away from the global composition due to molecular interactions.
Figure 6. For given (ω*12/kBT)(1) and (ω*12/kBT)(2), the logarithmic (2) ratio of the two resulting probability distributions p(1) 12,j and p12,j is a linear function of x12,j.
the two distributions is revealed to be a linear function of x12,j, with the limiting value
Figure 4. Equilibrium probability distributions for a one-dimensional lattice at various dimensionless interchange energies: repulsion (+0.7), random mixing (0.0), and attraction (−0.7).
⎧ ⎛ p(2) ⎞⎫ ⎪ 12, j ⎪ lim ⎨ln⎜ (1) ⎟⎬ = const ⎜ ⎟⎪ x12, j → 0⎪ ⎩ ⎝ p12, j ⎠⎭
As a consequence of this shift, the contributions to Shannon entropy, (23a), also shift accordingly as illustrated in Figure 5. Figure 5 also illustrates that, in the case of strong molecular interactions, the decaying outer parts of the degeneracy function, cf. Figure 3, become more relevant, in the sense that degeneracies somewhat apart from the peak of the degeneracy function, which is found at global composition, are
(36)
in which the limit x12,j → 0 can be reached by enlarging the system toward infinity, as illustrated in Table 1. Combining (36) with the internal energies, (Uj/kBT)(1) and (Uj/kBT)(2), the following correlation between the two equilibrium probability distributions can be established: ⎧ ⎛ p(2) ⎞⎫ ⎛ p(2) ⎞ ⎪ 12, j 12, j ⎪ ln⎜ (1) ⎟ − lim ⎨ln⎜ (1) ⎟⎬ ⎜ p ⎟ x12,j → 0⎪ ⎜ p ⎟⎪ ⎝ 12, j ⎠ ⎩ ⎝ 12, j ⎠⎭ const ⎡⎛ U ⎞(2) ⎛ U ⎞(1)⎤ j j = − ⎢⎜ ⎟ −⎜ ⎟ ⎥ ⎢⎝ kBT ⎠ k T ⎝ ⎠ ⎥⎦ B ⎣
∀j (37)
Since the probabilities of energy classes can be interpreted as the relative frequencies of energetic populations, cf. discussion of eq 32, eq 37 can be considered to be an analogue of Boltzmann’s law in classical statistics
Figure 5. Contributions to Shannon entropy at strong interaction levels, along with the degeneracy function of the system (N1 + N2 = 105, x1 = 0.3). With increasing deviations from random mixing, the decaying parts of the degeneracy function become more relevant.
⎛ U − U2 ⎞ n1 = exp⎜ − 1 ⎟ n2 kBT ⎠ ⎝ F
(38) DOI: 10.1021/acs.iecr.5b04430 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research where n1 and n2 represent the mean occupation numbers at energies U1 and U2.11 Sensitivity of the Result. The shape of the resulting equilibrium distributions, cf. Figure 4, indicates that the sensitivity of Shannon entropy, eq 23a, should be investigated, considering a limited number of energy classes, as opposed to all, around the mean local composition, designating the peak of the distribution. This is shown in Table 2. Table 2. Consequences of Considering Only Energy Classes in a Limited Range around x1̅ 2 When Summing Up the Terms of the Shannon Entropy (23a), ΔSrel = (Srange − Sall states)/Sall states·100a N1 + N2 range around x1̅ 2 [±%] 9 103 10 11 105 0.9 1 1.1 107 0.09 0.1 0.11
% of all states 0.93 1.38 2.83 6.29 9.80 1.42 1.66 1.13 7.22
× × × × × ×
10−446 10−441 10−435 10−48777 10−48723 10−48670
% of all energy classes
ΔSrel [%]
10.67 11.67 13.00 1.07 1.19 1.31 0.11 0.12 0.13
−1.93 −1.05 −0.43 −1.94 −0.93 −0.42 −1.93 −0.93 −0.43
Figure 7. Absolute values of the relative deviations between entropies predicted by the discrete model, SSh, and those given by the Guggenheim model, SGu, for a one-dimensional lattice of 103 respectively 105 sites at a global composition of x1 = 0.3. Minimum * /kBT = 0. deviations occur for the random system, ω12
model, and those given by the Guggenheim model. The corresponding deviations of calculated energies are at most 0.13% for the 103 lattice and at most 0.0013% for the 105 lattice, both maxima occurring at strongest repulsion, ω*12/kBT = 0.7. Hence, the models practically coincide for systems of ≥105 sites, which is in conformity with the expectations; cf. Figure 2 and Table 3. 4.2. Application to a Cubic Lattice. To determine all possible discrete energy classes within a simple cubic system, represented by discrete N12,j or x12,j, a rigorous algorithm was developed.12 As a rule of thumb, neglecting irregularities that mainly occur in the lower region, energy classes occur at an increment of 2 within
a
The absolute numbers of states and energy classes are given in Table 1.
For example, in a system with 105 sites, considering only states in the region x1̅ 2 ± 1% (i.e., only 357 of 3 × 104 possible energy classes (1.19%)), a relative deviation in entropy of N2
⎡ 1 (zN − N )⎤ ! ( 12 N12)!⎤⎦ ⎣2 2 12 ⎦ 2
N22 (40)
where N11 and N22 represent the number of contacts between molecules of the same type and N12 = N21 between molecules of different types. Since this function is then used within the framework of classical statistics, in the original paper,10 first Stirling’s formula was applied to replace the factorial terms by logarithms, and second the overall sum of degeneracies was replaced by its maximum term. Neither of these simplifications, which were performed for systems of infinite size, are used here because the degeneracy function will be applied to a system of finite size. Instead, normalization of (40) to fulfill eq 4 was
·100 [%]
−0.61 −0.08 −0.01 −0.001 −0.0001 G
DOI: 10.1021/acs.iecr.5b04430 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research conducted as the next step. This is achieved by determining a function h(N1,N2) in such a way that n
gj = h(N1 , N2) g (̃ N1 , N2 , N12, j)
and
!
∑ {gj} = m j=1
(41)
which results in h(N1 , N2) =
m n ∑k = 1 {gk̃ }
(42)
Combining (41) and (42) yields the discrete degeneracy function ∼ gj(N1 , N2 , N12, j) gj = m n ∑k = 1 {gk̃ } (43)
Figure 8. Relative deviations between Shannon entropy calculated by using the discrete model, SSh, and entropy from Monte Carlo simulations, SMC, using two different degeneracy functions, based on a cubic lattice with 105 sites and a global composition of x1 = 0.3.
Using a Degeneracy Function Obtained from Monte Carlo Simulations. To explore an alternative way to determine degeneracies apart from combinatorial analysis, Monte Carlo simulations were initially used to derive degeneracies, using the classical Metropolis algorithm as a basis.14−16 Starting from two directions, on one hand from a closepacked arrangement of one component, and on the other hand from a checkered, widespread arrangement, 108 random steps were performed in each case to canvas the range between N12,min and N12,max. In this way, approximately 2/3 of all possible energy classes could be determined. Due to the fragmentary nature of this method, a normal distibution was used as an approximating function for degeneracy. As an example, for an edge length of 50 and a global composition of x1 = 0.3, the following degeneracy function was derived and used for the subsequent model comparisons: dx12 =
dN12 2 = zN2 6N2
5. CONCLUSION Compared to the classical statistical-mechanical approach, the discrete model is based on systems of finite size for which degeneracy functions can be determined by Monte Carlo simulation and need not be estimated from combinatorial considerations, which are difficult to apply to more complex systems. The inclusion of Shannon entropy as an essential part of discrete modeling suggests that a different viewpoint must be taken when considering nonideal behavior in thermodynamic models: rather than making a strict distinction between “combinatorial” and “residual” contributions, the system should rather be understood as an aggregation of possible states, that are subject to constraints due to energetic interactions. While the discrete modeling of the ideal gas was based on the discrete states of noninteracting molecules, which represent statistically independent subsystems, the sites of a lattice represent concentrations of strongly interacting molecules that are not statistically independent. To account for this situation, discrete energy classes of the system were used as a basis in the present paper, which represent the most simple, rather coarsegrained, level of discretization, which is also the basis of classical statistical thermodynamics. Whereas classical thermodynamics a priori restricts itself to infinite systems, the key benefit of discrete modeling is its applicability to systems of finite size. The lower limit of system size is determined by the demand for the homogeneity of the Shannon entropy; for the case investigated in this paper (spherical molecules), a lattice edge length of 50 was shown to be sufficient, representing a system size that is typical for Monte Carlo simulations. Monte Carlo simulations can effectively be used to determine both the degeneracy functions and the possible energy classes of the system, when extending the model to nonspherical molecules for which degeneracies based on combinatorial considerations are difficult to derive. These functions are the main prerequisite for the application of discrete modeling and
(44a)
μ = x1 = 0.3
(44b)
σ = 0.00052
(44c)
g (x12, j) =
discussion about Figure 5, this hints at the inaccuracy of the Monte Carlo degeneracy in the sense that degeneracies somewhat apart from the peak of this function, which is found at global composition, are not accurate enough. Hence, there is potential to improve the aforementioned method to determine the Monte Carlo degeneracy, (44d), to get correct degeneracies over the entire range of the discrete local compositions, not only around the peak of this function.
⎛ 1 ⎛ x12, j − μ ⎞2 ⎞ (N1 + N2)! dx12 ⎟ ⎟ exp⎜ − ⎜ ⎠⎠ σ N1! N2! σ 2π ⎝ 2⎝ (44d)
Comparison with Monte Carlo Simulations. Monte Carlo simulations for a lattice with an edge length of 50 were the basis for the model comparison.16 The Monte Carlo entropy, derived from the internal energy using the Gibbs−Helmholtz equations, was compared to Shannon entropy, first using the Guggenheim degeneracy function, (43), and second using the Monte Carlo function, (44d). The results of this comparison are illustrated in Figure 8. Using the Guggenheim degeneracy function within the discrete model yields exactly the same results as when the Guggenheim model itself is used. This is an expected consequence, due to the equivalence of the maximum entropy principle, which is the basis of the discrete model, and statistical mechanics, which is the basis of the Guggenheim approach, as shown in a more general context in a previous work.6 Using the Monte Carlo degeneracy function is slightly advantageous as compared to the use of the Guggenheim model at increasing interaction levels, yet displays somewhat higher deviations at low interaction levels. With reference to the H
DOI: 10.1021/acs.iecr.5b04430 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research
3(p12, j ) = Lagrange function, used for constrained
the crucial point for its success, because they contain the entire combinatorial picture of the system, implicitly accounting for three-dimensional, orientation-dependent interactions between the species considered. As a first step toward application in chemical engineering and to approach more realistic systems, the model can immediately be extended to nonspherical molecules. This can be achieved by determining degeneracy functions for nonspherical molecules by Monte Carlo simulations of a finite system, and then applying the same procedure as shown in this paper for spherically symmetrical molecules. The parameters of the discrete model are the parameters of the degeneracy function. Future work will focus on two key topics: first, the refinement of Monte Carlo techniques with the aim of determining degeneracy functions for nonspherical molecules, in order to cover the range of possible energy classes as thoroughly as possible while spending reasonable amounts of simulation time; second, a further refinement of the discretization level from discrete energy classes of a system toward discrete states of individual molecules. One obvious approach could be to link the parameters of degeneracy functions to molecular properties of the components, initially to the van der Waals constants for the molecular structures as in UNIQUAC, and next, to spatial charge distributions as supplied by quantum-chemical calculations. In addition, the use of occupation numbers for energy classes of individual molecules will be investigated as a basis for discrete modeling.
■
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.5b04430. Constrained maximation of Shannon entropy in detail, for both random systems and nonrandom systems considering an energy constraint; outline of the combinatorial analysis of the one-dimensional lattice (PDF)
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Tel.: +43 (0)316 873 7966. Fax: +43 (0)316 873 7469. Notes
The authors declare no competing financial interest.
■ ■
■
ACKNOWLEDGMENTS The authors gratefully acknowledge support from NAWI Graz. Sara Crockett is thanked for English language editing.
maximation of entropy λ1, λ1 = Lagrangian multipliers, used for constrained maximation of entropy m = number of feasible states of a lattice system of given N1 and N2 n = number of feasible energy classes of a lattice system of given N1 and N2 n1, n2 = mean occupation numbers of systems representing energies U1 and U2 N = total number of molecules in the system N1, N2 = number of molecules of types 1 and 2 in the system N12 = total number of contacts between molecules of types 1 and 2 in the system N̅ 12 = mean contact number of the system in equilibrium N12,j = number of contacts between molecules of types 1 and 2 in energy class j ω*12 = interchange energy, a combination of the interaction energies pi = probability of the system residing in either of m possible states pj,k = probability of the system representing energy class j and the specific arrangement k p12,j = probability of the system residing in energy class j, which is characterized by the discrete contact number N12,j p∞ 12,j = probability of the random system residing in energy class j q = Boltzmann factor S = entropy of the lattice system SiM = entropy of ideal mixing T = temperature U = internal energy of the lattice system Uj = internal energy of the lattice system residing in energy class j UE = excess internal energy of the lattice system UEj = excess internal energy of the lattice system residing in energy class j V = volume of the system x1, x2 = molar fractions of molecules of types 1 and 2 in the lattice system, representing global compositions x12,j = discrete local composition of the lattice system residing in energy class j x1̅ 2 = mean local composition of the lattice system x∞ 1̅ 2 = mean local composition of the random system z = coordination number of the lattice system
REFERENCES
(1) Prausnitz, J. M.; Lichtenthaler, R. N.; de Azevedo, E. G. Molecular Thermodynamics of Fluid-Phase Equilibria, 3rd ed.; Prentice Hall, Inc.: 1999. (2) Thormann, M.; Klamt, A.; Wichmann, K. COSMOsim3D: 3DSimilarity and Alignment Based on COSMO Polarization Charge Densities. J. Chem. Inf. Model. 2012, 52, 2149−2156. (3) Bronneberg, R.; Pfennig, A. MOQUAC, a new expression for the excess Gibbs energy based on molecular orientations. Fluid Phase Equilib. 2013, 338, 63−77. (4) Pfleger, M.; Wallek, T.; Pfennig, A. Constraints of Compound Systems: Prerequisites for Thermodynamic Modeling Based on Shannon Entropy. Entropy 2014, 16, 2990−3008. (5) Pfleger, M.; Wallek, T.; Pfennig, A. Discrete Modeling: Thermodynamics Based on Shannon Entropy and Discrete States of Molecules. Ind. Eng. Chem. Res. 2015, 54, 4643−4654. (6) Jaynes, E. T. Information Theory and Statistical Mechanics. Phys. Rev. 1957, 106, 620−630.
LIST OF SYMBOLS dx12, μ, σ = parameters of a normal distibution used as approximating degeneracy function ε11, ε22, ε12 = interaction energies between the molecules f1(p12,j) = constraint function, used for maximization of entropy gj = degeneracy of energy class j of the lattice system g̃, h = elements of the Guggenheim degeneracy function i = index of the state, i.e. a specific arrangement of the lattice system j = index of the energy class of the lattice system k = index kB = Boltzmann constant I
DOI: 10.1021/acs.iecr.5b04430 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
Article
Industrial & Engineering Chemistry Research (7) Jaynes, E. T. Information Theory and Statistical Mechanics. II. Phys. Rev. 1957, 108, 171−189. (8) Ingarden, R. S. Information Theory and Variational Principles in Statistical Theories. Bull. Acad. Pol. Sci., Ser. Sci., Math., Astron. Phys. 1963, 11, 541−547. (9) Guggenheim, E. A. Mixtures; Clarendon Press: 1952. (10) Guggenheim, E. A.; McGlashan, M. L. Statistical Mechanics of Regular Mixtures. Proc. R. Soc. London, Ser. A 1951, 206, 335−353. (11) Maczek, A. Statistical Thermodynamics, 3rd ed.; Oxford University Press: 1998. (12) Mayer, C. Kombinatorische Diskretisierung von Gittersystemen. Thesis, Graz University of Technology, 2015 (in German). (13) Pielen, G. Detaillierte molekulare Simulationen und Parameterstudie für ein ternäres Gemisch zur Weiterentwicklung des GEQUACModells. Ph.D. Thesis, RWTH Aachen, 2005 (in German). (14) Zapf, F. Monte-Carlo Verfahren zur Diskretisierung von Gittersystemen. Thesis, Graz University of Technology, 2015 (in German). (15) Maltezos, G. Simulation lokaler Zusammensetzungen. Thesis, RWTH Aachen, 1987 (in German). (16) König, L. Auswertung von Monte Carlo-Simulationen zur Validierung thermodynamischer Modelle. Thesis, Graz University of Technology, 2013 (in German).
J
DOI: 10.1021/acs.iecr.5b04430 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX
