Efficient simulation of chemical potentials and phase equilibria in associating fluids: monomer/dimer insertion versus gradual particle insertion in primitive water models

Widom's test particle method for simulation of chemical potentials fails for associating (network forming) fluids. Therefore we study more sophisticated insertion methods, such as gradual particle insertion and the recent unbonded particle insertion of Tripathi and Chapman. We model strongly associating fluids by short-ranged primitive models (PMs) due to Nezbeda and co-workers. Recently, we showed that gradual insertion, in principle, is applicable to PM water. Here, we study systematically subcritical chemical potential isotherms, determine vapour–liquid coexistence densities and estimate the critical point density and temperature. For comparison we implement two variants of the Tripathi–Chapman algorithm. In the first case we use as unbonded test particles only particles which remain single molecules (monomers), when inserted and in the second case such test particles which remain as single molecules or form a two-particle cluster with an existing free molecule (monomer/dimer), when inserted. While monomer insertion improves Widom's method slightly, monomer/dimer insertion extends substantially the range of application to the subcritical temperatures studied by gradual insertion. But, for low temperatures monomer/dimer insertion requires extremely long Markov chains and significantly more computer time than gradual insertion. The Tripathi–Chapman algorithm—a direct extension of Widom's method—is much simpler than gradual insertion and therefore should be preferred at supercritical and critical conditions.     

Chemical potential calculations in dense liquids using metadynamics

The calculation of chemical potential has traditionally been a challenge in atomistic simulations. One of the most used approaches is Widom's insertion method in which the chemical potential is calculated by periodically attempting to insert an extra particle in the system. In dense systems this method fails since the insertion probability is very low. In this paper we show that in a homogeneous fluid the insertion probability can be increased using metadynamics. We test our method on a supercooled high density binary Lennard-Jones fluid. We find that we can obtain efficiently converged results even when Widom's method fails.     

A new study of associating inhomogeneous fluids with classical density functional theory

ABSTRACT

A new density functional for the study of associating inhomogeneous fluids based on Wertheim's first-order thermodynamic perturbation theory is presented and compared to the most currently used associating density functionals. This functional is developed using the weighted density approximation in the range of association of hard spheres. We implement this functional within the framework of classical density functional theory together with modified fundamental measure theory to account for volume exclusion of hard spheres. This approach is tested against molecular simulations from literature of pure associating hard spheres and mixtures of non-associationg and associating hard spheres with different number of bonding sites close to a hard uniform wall. Furthermore, we compare and review our results with the performance of associating functionals from literature, one based on fundamental measure theory and the inhomogeneous version of Wertheim's perturbation theory. Results obtained with classical DFT and the three functionals show excellent agreement with molecular simulations in systems with one hard wall. For the cases of small pores where only one or two layers of fluid are allowed discrepancies between results with classical DFT and molecular simulations were found.     

Associating fluids with four bonding sites against a hard wall: density functional theory

We present two new perturbation density functional theories to investigate non-uniform fluids of associating molecules. Each fluid molecule is modelled as a spherical hard core with four highly anisotropic square well sites placed in tetrahedral symmetry on the hard core surface. In one theory we apply the weighting from Tarazona's hard sphere density functional theory to Wertheim's bulk first-order perturbation theory. The other theory uses the inhomogeneous form of Wertheim's theory as a perturbation to Tarazona's hard-sphere density functional theory. Each theory approaches Tarazona's theory in the limit of zero association. We compare results from theory and simulation for density profiles, fraction of monomers, and adsorption of an associating fluid against a hard, smooth wall over a range of temperatures and densities. The non-uniform fluid theory which uses Tarazona's weighting of Wertheim's theory in the bulk is in good agreement with computer simulation results.     

Ground state structures of a lattice model: honeycomb and tetrahedral lattices

The ground state structures of a lattice model with nearest-neighbor and next-nearest-neighbor interactions were examined by the method of “partitioning” the hamiltonian. This model is Widom's lattice model of microemulsions. Previous analysis had examined only cubic lattices, however the analysis is more general as this paper illustrates.     

The O-H distance in ice

The thermodynamic properties of normal and para-hydrogen are computed from multiple time-step path integral hybrid Monte Carlo (PIHMC) simulations. Four different isotropic pair potentials are evaluated by comparing simulation results with experimental data. The Silvera–Goldman potential is found to be the most accurate of the potentials tested for computing the density and internal energy of fluid hydrogen. Using the Silvera–Goldman potential, simulation and experimental data are compared on isobars ranging from 0.1 to 100 MPa and for temperatures from 18 to 300 K. The Gibbs free energy is calculated from the PIHMC simulations by an adaptation of Widom's particle insertion technique to a path integral fluid. A new method is developed for computing phase equilibria for quantum fluids directly by combining PIHMC with the Gibbs ensemble technique. This Gibbs–PIHMC method is used to calculate the vapour–liquid phase diagram of hydrogen from simulations. Agreement with experimental data is good.     

TPT2 and SAFTD equations of state for mixtures of hard chain copolymers

We present the second-order thermodynamic perturbation theory (TPT2) and the dimer statistical associating fluid theory (SAFTD) equations of state for mixtures consisting of hetero-nuclear hard chain molecules based on extensions of Wertheim's theory for associating fluids. The second-order perturbation theory, TPT2, is based on the hard sphere mixture reference fluid. SAFTD is an extension of TPT1 (= SAFT) and is based on the non-spherical (hard disphere mixture) reference fluid. The TPT2 equation of state requires only the contact values of the hard sphere mixture site-site correlation functions, while the SAFTD equation of state requires the contact values of site-site correlation functions of both hard sphere and hard disphere mixtures. We test several approximations for site-site correlation functions of hard disphere mixtures and use these in the SAFTD equation of state to predict the compressibility factor of copolymers. Since simulation data are available only for a few pure copolymer systems, theoretical predictions are compared with molecular simulation results for the compressibility factor of pure hard chain copolymer systems. Our comparisons show a very good performance of TPT2, which is found to be more accurate than TPT1 (= SAFT). Using a modified Percus-Yevick site-site correlation function SAFTD is found to represent a significant improvement over SAFT and is slightly more accurate than TPT2. Comparison of SAFTD with generalized Flory dimer (GFD) theory shows that both are equivalent at intermediate to high densities for the compressibility factor of copolymer systems investigated here.     

Cross-association of multi-component systems

The design and optimization of equipment in chemical industry (e.g. heat exchanger) and also process simulations require the knowledge of physical properties of mixtures, for instance the involved phase equilibria, enthalpies and the heat capacities. Most experimental data on these properties exists for pure compounds (e.g. water) and for binary mixtures (e.g. water ethanol). The database is however, very limited for mixtures of more than two species. Physically sound equations of state, like the Perturbed-Chain Statistical Associating Theory (PC-SAFT) have been used successfully to provide information about these thermodynamic properties for a wide variety of substances including systems of associating and non-associating or systems of associating and cross-associating species. One of the main challenges using this Wertheim-type equation of state is the mathematically implicit form of the underlying nonlinear system of equations, if association occurs. This article provides in depth information about our recently developed fast and stable algorithm to solve this system of equations numerically for multi-component systems, as well as a new method to find good initial values for the numerical algorithm. Furthermore, the numerical results are compared to experimental data on several properties of interest and found to be in good to excellent agreement.     

Existence of Weak Solutions for a Diffuse Interface Model for Viscous,Incompressible Fluids with General Densities

We study a diffuse interface model for the flow of two viscous incompressible Newtonian fluids in a bounded domain. The fluids are assumed to be macroscopically immiscible, but a partial mixing in a small interfacial region is assumed in the model. Moreover, diffusion of both components is taken into account. In contrast to previous works, we study the general case that the fluids have different densities. This leads to an inhomogeneous Navier-Stokes system coupled to a Cahn-Hilliard system, where the density of the mixture depends on the concentration, the velocity field is no longer divergence free, and the pressure enters the equation for the chemical potential. We prove existence of weak solutions for the non-stationary system in two and three space dimensions.     

Surface phase transitions of multiple-site associating fluids

Surface phase transitions of Lennard–Jones (LJ) based two- and four-site associating fluids have been studied for various associating strengths using grand-canonical transition matrix Monte Carlo simulations. Our results suggest that, in the case of a smooth surface, represented by a LJ 9-3-type potential, multiple-site associating fluids display a prewetting transition within a certain temperature range. However, the range of the prewetting transition decreases with increasing associating strength and increasing number of sites on the fluid molecules. With the addition of associating sites on the surface, a quasi-2D vapor–liquid transition may appear, which is observed at a higher surface site density for weaker associating fluids. The prewetting transition at lower associating strength is found to shift towards the quasi-2D vapor–liquid transition with increasing surface site density. However, for highly associating fluids, the prewetting transition is still intact, but shifts slightly towards the lower temperature range. Adsorption isotherms, chemical potentials and density profiles are used to characterize surface phase transitions.     

Coupled flow-polymer dynamics via statistical field theory: Modeling and computation

Field-theoretic models, which replace interactions between polymers with interactions between polymers and one or more conjugate fields, offer a systematic framework for coarse-graining of complex fluids systems. While this approach has been used successfully to investigate a wide range of polymer formulations at equilibrium, field-theoretic models often fail to accurately capture the non-equilibrium behavior of polymers, especially in the early stages of phase separation. Here the “two-fluid” approach serves as a useful alternative, treating the motions of fluid components separately in order to incorporate asymmetries between polymer molecules. In this work we focus on the connection of these two theories, drawing upon the strengths of each of the approaches in order to couple polymer microstructure with the dynamics of the flow in a systematic way. For illustrative purposes we work with an inhomogeneous melt of elastic dumbbell polymers, though our methodology will apply more generally to a wide variety of inhomogeneous systems. First we derive the model, incorporating thermodynamic forces into a two-fluid model for the flow through the introduction of conjugate chemical potential and elastic strain fields for the polymer density and stress. The resulting equations are composed of a system of fourth order PDEs coupled with a non-linear, non-local optimization problem to determine the conjugate fields. The coupled system is severely stiff and with a high degree of computational complexity. Next, we overcome the formidable numerical challenges posed by the model by designing a robust semi-implicit method based on linear asymptotic behavior of the leading order terms at small scales, by exploiting the exponential structure of global (integral) operators, and by parallelizing the non-linear optimization problem. The semi-implicit method effectively removes the fourth order stability constraint associated with explicit methods and we observe only a first order time-step restriction. The algorithm for solving the non-linear optimization problem, which takes advantage of the form of the operators being optimized, reduces the overall simulation time by several orders of magnitude. We illustrate the methodology with several examples of phase separation in an initially quiescent flow.     

Statistical approach to the kinetics of nonuniform fluids

We present a new formalism in Fourier space for the study of spatially nonuniform fluids in nonequilibrium states which generalizes previous work on uniform fluids. Starting from the Liouville equation we obtain a hierarchy of equations for the reduced distribution functions which gives their rate of change at any given order of the system mean density as a sum of a finite number of terms. Using a finite-ranged repulsive interaction potential we derive, as a first application of the formalism, the Boltzmann integrodifferential equation for an infinite system which is initially in a “weakly” inhomogeneous state. This is accomplished introducing an initial statistical assumption, namely initial molecular chaos; this condition is seen to hold during the time evolution described by the resulting kinetic equation.     

Modelling the effect of methanol,glycol inhibitors and electrolytes on the equilibrium stability of hydrates with the SAFT-VR approach

In this work we integrate the statistical associating fluid theory for fluids interacting through potentials of variable range (SAFT-VR) into a traditional van der Waals and Platteeuw framework for modelling clathrate hydrates. We incorporate a new water–guest cell potential for the hydrate phase that can be related to the potential adopted in the familiar SAFT-VR equation of state for modelling fluids. We show how the ability of this equation of state to treat a wide range of complex fluids increases the scope of hydrate modelling to incorporate, in a single framework, the presence of various inhibitors (alcohols, glycols) or brines – or, indeed, any fluid for which a model is available (for use within SAFT-VR) or can be conveniently obtained. Agreement with experimental results is good throughout and, in many cases, excellent.     

Auxiliary field Monte-Carlo for quantum many-body ground states

We develop an algorithm for determining the exact ground state properties of quantum many-body systems which is equally applicable to bosons and fermions. The Schroedinger eigenvalue equation for the ground state energy is recast as a many-dimensional integral using the Hubbard-Stratonovitch representation of the imaginary-time many-body evolution operator. The integral is then evaluated stochastically. We test the algorithm for an exactly soluble boson system with an attractive potential and then extend it to fermions and repulsive potentials. Importance sampling is crucial to the success of the method, particularly for more complex systems. Computational efficiency is improved by performing the calculations in Fourier space.     

A two-dimensional model associating fluid with spherically symmetric intracore square-well shell. Integral equations and Monte Carlo simulation study

The structure and thermodynamics of a monolayer of an associating fluid in the framework of the primitive model of Cummings and Stell is studied by using a two-dimensional approximation. The model permits formation of dimer species for small values of the bonding length parameter, the formation of chains, if the bonding length is slightly larger, and also the vulcanization of species for bonding length values close to the diameter of particles. The structure and thermodynamics of the model that are of interest for statistical mechanics of surface chemical reactions, are studied by computer simulations in the canonical, grand canonical and isobaric ensembles and from the two-dimensional Ornstein—Zernike-like or Wertheim's Ornstein—Zernike integral equation. We have shown that the theory is satisfactory for the case of dimerization if the fluid density is low. For higher densities one must apply a correction for the cavity distribution functions to describe the fraction of unbonded species more adequately. For the case of chain formation the theory just resembles trends following from the simulation data. For the model with vulcanization of species we have obtained the fractions of differently bonded particles and have shown that chemical ordering of species is manifested in the antiphase oscillations of the pair distribution functions of species.     

Systems of oscillators with statistical energy exchange in collisions

Classical transition probabilities are derived for the exchange of energy between sets of weakly coupled harmonic oscillators, subject to the assumption that energy is ‘statistically redistributed’ in all collisions. In the case of single oscillators the result is particularly simple and an explicit solution of the relaxation equation for such a system is derived. A ‘mean first passage time’ for the dissociation of truncated harmonic oscillators is also obtained using Widom's theory and the limiting case of equilibrium reaction is discussed.

The same approach is sketched briefly for the case of quantized systems in which statistical redistribution of energy is allowed.

Within the general scope of the model, these transition probabilities represent the most efficient conceivable coupling between system and heat bath.     

The classical fluid of associating hard rods in an external field

Two models of one-dimensional fluids of associating hard rods in an arbitrary external field are investigated. In the first model particles can only form dimers, while in the second model, which has been solved previously by Percus, aggregates of any size coexist. In both cases the grand canonical potential and the external potential are found exactly as functionals of the density. It is shown that Wertheim's thermodynamic perturbation theory of polymerization provides a straightforward route to the exact solution by expanding the functional space to include more density parameters. This suggests that Wertheim's theory should be used also for studying the structure (and not only the thermo-dynamics) of real associating fluids.     

Predicting the solubility of xenon in n-hexane and n-perfluorohexane: a simulation and theoretical study

The solubility of xenon in n-hexane and n-perfluorohexane has been studied using both molecular simulation and a version of the SAFT approach (SAFT-VR). The calculations were performed close to the saturation line of each solvent, between 200 K and 450 K, which exceeds the smaller temperature range where experimental data are available in the literature. Molecular dynamics simulations, associated with Widom's test particle insertion method, were used to calculate the residual chemical potential of xenon in n-hexane and n-perfluorohexane and the corresponding Henry's law coefficients. The simulation results overestimate the solubility of xenon in both solvents when simple geometric combining rules are used, but are in good agreement if a binary interaction parameter is included. With the SAFT-VR approach we are able to reproduce the experimental solubility for xenon in n-hexane, using simple Lorentz-Berthelot rules to describe the unlike interaction. In the case of n-perfluorohexane as a solvent, a binary interaction parameter was introduced, taken from previous work on (xe + C₂F₆) mixtures. Overall, good agreement is obtained between the simulation, theoretical and experimental data.     

Monte-Carlo simulation of an inhomogeneous reaction-diffusion system in the biophysics of receptor cells

The numerical simulation of Markov processes is usually performed by means of the minimal process method or Gillespie algorithm. In reaction-diffusion systems including extremely inhomogeneous situations, the direct application of this algorithm meets with severe difficulties which eventually cause extremely large computing times. We present a modification of the minimal process method which make it applicable to such situations even on small size computers within moderate computing times. Our modifications include the use of logarithmic classes for transition probabilities in order to increase the acceptance rate of von Neumann rejection methods, a non-local storage of the lattice state, hash tables and dynamical storage management in order to save memory capacity. Our actual example for demonstrating our modified algorithm is signal transduction in biological receptor cells where transmitter molecules are released on a two-dimensional structure (cell membrane) after a quantum reception event, e.g. a photon capture in photoreceptor cells, and interact with ionic transport channels. We find satisfying agreement of our simulation results with the experimental data for the ventral nerve photoreceptor of Limulus.