CET Volume 86 DOI: 10.3303/CET2186214 Paper Received: 10 October 2020; Revised: 25 January 2021; Accepted: 20 April 2021 Please cite this article as: Mayer C., Wallek T., 2021, Discrete Excess Gibbs-energy Modeling Approach Based on Clusters of Molecules, Chemical Engineering Transactions, 86, 1279-1284 DOI:10.3303/CET2186214 CHEMICAL ENGINEERING TRANSACTIONS VOL. 86, 2021 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Sauro Pierucci, Jiří Jaromír Klemeš Copyright © 2021, AIDIC Servizi S.r.l. ISBN 978-88-95608-84-6; ISSN 2283-9216 Discrete Excess Gibbs-energy Modeling Approach Based on Clusters of Molecules Christoph Mayer*, Thomas Wallek Institute of Chemical Engineering & Environmental Technology, Graz University of Technology, Inffeldgasse 25/C, 8010 Graz, Austria cmayer@tugraz.at This paper presents a discrete modeling approach that is based on clusters of molecules. The aim of the approach is to provide the excess Gibbs function (gE) resp. activity coefficients of condensed phase mixtures, in particular, liquids. The modeling approach uses the probabilities of the discrete states of molecular clusters as its main variables. The clusters consist of four molecules. The compounds are modeled as dice-like molecules, giving them the option of having one to six different energetic interaction sites. A molecular sampling algorithm links the abstract dice-like representations to real molecules. The model results are compared to experimental data for mixtures of butanal + n-heptane, acetone + n-heptane, and acetone + methanol. The comparison shows that the deviations of this approach are similar to those of the well-known UNIFAC model, which indicates the approach as a promising alternative for the development of gE models. 1. Introduction Equations of state and activity coefficient models, which are thermodynamic approaches to calculate fluid phase equilibria, are being challenged to describe mixtures with pronounced interactions. These systems are hard to describe using traditional approaches, due to their significant divergences from random mixing behavior. However, due to the upcoming change towards more bio-based raw materials, the modeling of these systems is gaining in importance. Previous work introduced ‘discrete modeling’ as an attempt to include more fine-grained molecular information into thermodynamics on a fundamental level. The use of Shannon’s measure of information as an equivalent to the thermodynamical entropy (Pfleger et al. 2015; Wallek et al. 2016) is the defining characteristic of this approach. The basis of the approach is a former paper using discrete Markov chains to model the thermodynamics of planar (2D) lattices for solid solutions (Vinograd and Perchuk, 1996). This modeling concept was further developed towards a three-dimensional system (Wallek et al., 2018), where the compounds were modeled with a geometry similar to dice, providing one to six distinguishable interaction sites per molecule (Mayer and Wallek, 2020). The basis for thermodynamic modeling of this approach is a cluster comprising several molecules that is intended as a representative sub-system of the lattice and consists of four molecules and three pairwise contacts. In order to link this modeling approach to real molecules, the energetic interactions between molecules are determined using a sampling algorithm similar to the PAC-MAC approach (Sweere and Fraaije, 2015). 2. Model for dice-like molecules The modeling approach is based on a simple cubic lattice which is populated with the dice-like molecule representations. The core idea behind this approach based on Vinograd and Perchuk (1996) is that the process of inserting molecules into this lattice is sequential in nature, inserting the molecules one by one. Each newly inserted molecule is surrounded by a partially complete proximity of formerly inserted molecules. Figure 1 shows an example of such an insertion step. The insertion process itself is influenced by the energetic interactions between the newly inserted molecule and its immediate neighbors. Therefore, conditional probabilities can be assigned to each insertion step based on the selection of component type and 1279 orientation of the newly inserted molecule with respect to the specific neighborhood it is placed into. The aim of this modeling approach is to determine the averaged probabilities for the construction of a lattice system that has been thermodynamically equilibrated. The probabilities are determined though a minimization of the free energy expressed as a function of these probabilities. The following paragraphs provide an overview of this modeling approach on a conceptional level. A detailed description including the full system of equations can be found in Mayer and Wallek (2020). Figure 1: Sequential assembly of a lattice by the insertion of a new molecule at the position A with neighbors B, C, and D already being present from previous insertions. 2.1 Cluster as a modeling basis To achieve the introduced modeling concept, a small cluster is chosen from the entire lattice system and used as a representative modeling basis. The cluster is chosen in a way that incorporates the immediate neighborhood that the new molecule is placed into. An example of such a cluster is given in Figure 2. The cluster given in this figure can be interpreted such that dice A represents the newly inserted molecule and dice B, C, and D represent the next neighbors which are present at the time of insertion. With this choice of cluster, the probabilities of the sequential lattice insertion steps can also be modeled based on this small cluster. Therefore, the thermodynamic property of interest, i.e., the free energy, can be expressed by making use of the probabilities from the small cluster. Figure 2: Example of a cluster comprising four dice-like molecules where the dice color signals the component type. 2.2 Entropy and internal energy The entropy of the model is based on Shannon entropy applied to the conditional insertion probabilities. It is important that the entropy of the sequential insertion step is used as the system’s entropy. If the entropy based only on the clusters were to be used, the described system would be an ensemble of independent clusters instead of the desired lattice-system. The system’s internal energy is assumed to be additive, and interactions further apart than nearest neighbors are assumed to have a negligible contribution to the entire internal energy. Considering these assumptions, 1280 the internal energy of this system equals the average cluster energy, that is the sum of the energies of all possible clusters weighted by their respective probability of occurrence. Figure 3: From the cluster given in Figure 2, the contacting pairs of sites are highlighted. The contacts are between molecules: (a) A-B, (b) A-C, (c) A-D. The cluster interaction energy in turn is calculated from the sum of the interactions between the three contacting pairs. These pairs are highlighted in Figure 3. The model requires a pairwise interaction energy for every possible pair of sites, which are 78 in total, as the coefficients defining the component types used in the binary mixture. 2.3 Free energy The Helmholtz function is used as objective function for the optimization that determines the probability distribution. Its minimum describes the desired thermodynamic equilibrium. Eq (1) shows the relationship between the free energy, a, the internal energy, u, and the entropy, s. Through this it can be immediately expressed as a function of the model variables, which are the cluster probabilities. = − (1) Furthermore, the primary target is the determination of the excess Gibbs-energy. It can be directly determined from the Helmholtz function. Because the model considers all molecules to be placed at fixed positions on a regular lattice, the excess volume can only be zero. Therefore, the excess free energy and the excess Gibbs- energy can be considered identical under the chosen model assumptions. 2.4 Constraints Aside from the target function, there are three types of constraints that take part in the optimization. The first type can be summarized as mathematical constraints. They are based on the law of total probability and provide the model with a link between the global composition, which is a desired input of the model, and the optimization variables. The second type of constraints reflects some of the fundamental properties of the lattice system. A key part of the model is that the average probabilities that lead to the construction of a lattice system are independent of the direction that the lattice system is viewed from. This means that every local directional preference in a single lattice has counterparts when looking at an ensemble of lattice systems that have a similar preference in other directions. This property of the lattice system means that the lattice has to be considered isotropic. It is this isotropic property that informs this type of constraints, which is important for reducing the number of independent variables that the numerical optimization has to handle. Preliminary investigations into the modeled system have yielded the fact that additional constraints are necessary to fully reach the desired results. The approach chosen for these additional constraints is based on the considerations regarding the additivity property of the cluster energy. Figure 3 introduced the fact that the cluster energy is equal to the sum of the three pairwise interactions present in the chosen cluster shape. It seems natural to also consider this connection between the cluster and the three pairs in terms of probabilities. The third type of constraints is based on approximating the entire cluster by sequentially constructing it out of the three pairs in a fashion that is inspired by the sequential construction of the entire lattice-system. This type of constraints completes the optimization system that, for a specific mixture of components, takes the global composition as molar fractions and the desired temperature as input and returns the Helmholtz free energy as well as the probability distribution in equilibrium as output. 1281 3. Link between dice-like and real molecules While the consideration of abstract dice-like molecules seems like a natural intermediate step during model development between simple uniform, spherical molecules and more realistic molecule models, it is not entirely intuitive how real molecules can be linked with these abstract representations. The following sections introduce one methodology that shows how this can be achieved. For the representation of real molecules, the ‘optimized potential for liquid simulations’ in its ‘all-atoms’ version (OPLS-AA, cf. Jorgensen et al., 1996) is used. One orientation per component is chosen as the default orientation. From this point on, other molecule orientations are limited to 90 degree rotations of the default orientation. This is analogous to the possible 24 orientations of dice. Therefore, every component has a discrete set of 24 possible orientations in which they can exist in the cluster. Figure 4 illustrates three of the twenty-four feasible orientations for acetone. Figure 4: Three example orientations of acetone when modeled with dice-like geometry. 3.1 Sampling procedure for pairwise interaction energies The interaction energies are input into the model as a matrix of pairwise interactions for every combination of dice faces. Because of this, the interaction energies can be determined by individual pairs of molecules. The sampling procedure for pairs of molecules is as follows: A specific orientation is selected for each molecule. Then, one molecule is fixed in the origin and the other one starts sufficiently far apart that overlapping cannot occur. The second molecule is then gradually moved closer to the first one, until a certain distance is reached. In the case of the example systems presented here, the second molecule is moved until the van der Waals- surfaces touch. The interaction energy is then evaluated and stored with a reference to the orientations that the molecules are in. After that, a new combination of orientations is selected and the procedure starts over. The whole process is performed until all pairings of orientations have been sampled. For a binary mixture this amounts to 1728 samples. The model itself differentiates only between the interaction sites involved in the contact and not the specific orientations of other sites on these molecules. Therefore, all samples that have the same sites involved in the contact are combined to classes and the mean interaction energies of these 78 classes are used as model input. 4. Results The results of the proposed modeling approach are illustrated using three example systems. They comprise butanal + n-heptane, acetone + n-heptane, and acetone + methanol. The model results are then compared with measurement data, results from the UNIFAC model, originally proposed by Fredenslund et al. (1975), as frequently used approach, as well as Monte Carlo simulations. The latter are based on the same dice-like molecule geometry as the model. The molecules are placed on a lattice and random rotations and location exchanges are performed using a standard Metropolis algorithm. The simulation receives the same set of pairwise interaction energies as the model and is used as a verification of model assumptions. Figure 5 gives the excess Gibbs function for a mixture of butanal + n-heptane at ambient temperature (298.2 K). Here, the resulting curve is drawn over the molar fraction of butanal. Comparisons are made to experimental data which are taken from Krenzer (1985) in Redlich-Kister polynomial form. It can be seen that the side of the curves that is rich in butanal shows similar deviations from the experimental data for both UNIFAC and the model presented in this work. A higher n-heptane composition leads to noticeably larger deviations. An explanation for this is the fact that n-heptane has an elongated form and is therefore not that well described using a dice-like geometry. Beside the section between a composition of 0.5 and 0.6, the model returns data which is slightly larger than that for the Monte-Carlo simulations. 1282 Figure 5: Excess Gibbs-energy for the mixture butanal + n-heptane is plotted over the molar fraction of butanal. Measurements (Krenzer 1985), UNIFAC (Fredenslund et al. 1975), and Monte Carlo data are also shown for comparison. The mixture acetone + n-heptane, shown in Figure 6 at ambient temperature (298.15 K), shows a similar behavior as the previous mixture. The deviations between this model and the measured data from literature (Krenzer 1985) are again of a similar magnitude than those of the UNIFAC model. The model calculates slightly larger results than the Monte-Carlo simulation. Figure 6: Excess Gibbs-energy for the system acetone + n-heptane is plotted over acetone composition. Measurements (Krenzer 1985), UNIFAC (Fredenslund et al. 1975), and Monte Carlo data are also shown for comparison. Figure 7 gives the excess Gibbs function of the system acetone + methanol plotted over the molar fraction of acetone at a temperature of 323.15 K. From Gmehling & Onken (2005), experimental data are taken again by using a Redlich-Kister polynomial expansion. For an acetone composition below 0.5, the model results are slightly beneath the experimental data and slightly higher at larger concentrations. Overall, the model is in satisfying agreement with measurement data. In comparison to results from Monte Carlo simulations, a slightly positive deviation can be observed. For the UNIFAC model, however, a more significant discrepancy from measurement data becomes evident for this two-component mixture. 1283 Figure 7: Excess Gibbs-energy for a mixture of acetone + methanol, plotted over acetone composition. Measurement data (Gmehling & Onken 2005), results from the UNIFAC model (Fredenslund et al., 1975), and data from Monte Carlo simulations are also shown for comparison. 5. Conclusions In this work, a model based on dice-like abstractions is coupled with a molecular-sampling procedure to allow real molecules to be projected into this abstract form. The model assumes that several restrictions are in place to account for the selected shape and properties. The first assumption is that all molecules have the same size and shape. This means that the model yields the best results when applied to components that are of similar size and are fairly spherical in shape. The next assumption is that each molecule is limited to the 24 possible orientations that can be observed from dice. Furthermore, the molecules are forced on fixed positions in a lattice. First results of the example systems of butanal + n-heptane, acetone + n-heptane, and acetone + methanol have shown that the abstract dice-like molecule representations with all the mentioned restrictions provide promising results regarding the underlying modeling approach. The applicability of this modeling approach to different mixtures consisting of a diverse variety of components is topic of ongoing research investigations. Acknowledgments This research was funded by the Austrian Science Fund (FWF), grant number P 32609-N. References Fredenslund A., Jones R. L., Prausnitz J. M., 1975, Group‐contribution estimation of activity coefficients in nonideal liquid mixtures, AIChE J., vol. 21, no. 6, 1086–1099. Gmehling J., Onken U., 2005, Vapor-Liquid Equilibrium Data Collection, Dechema Chemistry Data Series, Vol. I, Part 2g. Jorgensen W. L., Maxwell D. S., Tirado-Rives J., 1996, Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids, J. Am. Chem. Soc. 118, 11225–11236. Krenzer L., 1985, [in German] Untersuchungen zur Beeinflussung der Exzessenthalpie, der freien Exzessenthalpie und der Exzessentropie von binären Mischungen mit polaren Komponenten durch die Moleküleigenschaften, PhD thesis, Technische Hochschule Darmstadt. Mayer C., Wallek T., 2020, Cluster-Based Thermodynamics of Interacting Dice in a Lattice, Entropy, 22(10), 1111. Pfleger, M., Wallek, T., and Pfennig, A., 2015, Discrete Modeling: Thermodynamics Based on Shannon Entropy and Discrete States of Molecules., Ind. Eng. Chem. Res., 54, 4643–4654. Sweere, A. J. M.; Fraaije, J. G. E. M., 2015, Force-Field Based Quasichemical Method for Rapid Evaluation of Binary Phase Diagrams., J. Phys. Chem. B, 119, 14200−14209. Vinograd, V.L., Perchuk, L.L., 1996, Informational Models for the Configurational Entropy of Regular Solid Solutions: Flat Lattices., J. Phys. Chem., 100, 15972–15985. Wallek, T., Pfleger, M., and Pfennig, A., 2016, Discrete Modeling of Lattice Systems: The Concept of Shannon Entropy Applied to Strongly Interacting Systems., Ind. Eng. Chem. Res., 55, 2483–2492. Wallek, T., Mayer, C., Pfennig, A., 2018, Discrete Modeling Approach as a Basis of Excess Gibbs-Energy Models for Chemical Engineering Applications, Ind. Eng. Chem. Res., 57, 1294–1306. 1284