Local density dependent potentials for an underlying van der Waals equation of state: A simulation and density functional theory analysis.

There is an ever increasing use of local density dependent potentials in the mesoscale modeling of complex fluids. Questions remain, though, about the dependence of the thermodynamic and structural properties of such systems on the cutoff distance used to calculate these local densities. These questions are particularly acute when it comes to the stability and structure of the vapor/liquid interface. In this article, we consider local density dependent potentials derived from an underlying van der Waals equation of state. We use simulation and density functional theory to examine how the bulk thermodynamic and interfacial properties vary with the cutoff distance, rc, used to calculate the local densities. We show quantitatively how the simulation results for bulk thermodynamic properties and vapor-liquid equilibrium approach the van der Waals limit as rc increases and demonstrate a scaling law for the radial distribution function in the large rc limit. We show that the vapor-liquid interface is stable with a well-defined surface tension and that the interfacial density profile is oscillatory, except for temperatures close to critical. Finally, we show that in the large rc limit, the interfacial tension is proportional to rc and, therefore, unlike the bulk thermodynamic properties, does not approach a constant value as rc increases. We believe that these results give new insights into the properties of local density dependent potentials, in particular their unusual interfacial behavior, which is relevant for modeling complex fluids in soft matter.

Microstructural and thermodynamic characterization of wormlike micelles formed by polydisperse ionic surfactant solutions.

For industrial applications of self-assembled wormlike micelles, measurement and characterization of a micellar material's microstructure and rheology are paramount for the development and deployment of new high-performing and cost-effective formulations. Within this workflow, there are significant bottlenecks associated with experimental delays and a lack of transferability of results from one chemistry to another. In this work, we outline a process to predict microscopic and thermodynamic characteristics of wormlike micelles directly from rheological data by combining a more robust and efficient fitting algorithm with a recently published constitutive model called the Toy Shuffling model [J. D. Peterson and M. E. Cates, J. Rheol. 64, 1465-1496 (2020) and J. D. Peterson and M. E. Cates, J. Rheol. 65, 633-662 (2021)]. To support this work, linear rheology measurements were taken for 143 samples comprising a common base formulation of commercial sodium lauryl ether sulfate, cocamidopropyl betaine, and salt (NaCl). The steady state zero shear viscosity evident in linear rheology was measured in duplicate via direct steady and oscillatory shear experiments. Fitting the collected data to the model, we found trends in the microstructural and thermodynamic characteristics that agree with molecular dynamics simulations. These trends validate our new perspective on the parameters that inform the study of the relationship between chemical formulation and rheology. This work, when implemented at scale, can potentially be used to inform and test strategies for predicting self-assembled micellar structures based on chemical formulation.

Hierarchical Aggregation in a Complex Fluid─The Role of Isomeric Interconversion.

There is an ever-increasing body of evidence that metallic complexes involving amphiliphic ligands do not form normal solutions in organic solvents. Instead, they form complex fluids with intricate structures. For example, the metallic complexes may aggregate into clusters, and these clusters themselves may aggregate into superclusters. To gain a deeper insight into the mechanisms at play, we have used an improved force field to conduct extensive molecular dynamics simulations of a system composed of zirconium nitrate, water, nitric acid, tri-n-butyl phosphate, and n-octane. The important new finding is that a dynamic equilibrium between the cis and trans isomers of the metal complex is likely to play a key role in the aggregation behavior. The isolated cis and trans isomers have similar energies, but simulation indicates that the clusters consist predominantly of cis isomers. With increasing metal concentration, we hypothesize that more clustering occurs and the chemical equilibrium shifts toward the cis isomer. It is possible that such isomeric effects play a role in the liquid-liquid extraction of other species and the inclusion of such effects in flow sheet modeling may lead to a better description of the process.

A Telescoping View of Solute Architectures in a Complex Fluid System.

Short- and long-range correlations between solutes in solvents can influence the macroscopic chemistry and physical properties of solutions in ways that are not fully understood. The class of liquids known as complex (structured) fluids-containing multiscale aggregates resulting from weak self-assembly-are especially important in energy-relevant systems employed for a variety of chemical- and biological-based purification, separation, and catalytic processes. In these, solute (mass) transfer across liquid-liquid (water, oil) phase boundaries is the core function. Oftentimes the operational success of phase transfer chemistry is dependent upon the bulk fluid structures for which a common functional motif and an archetype aggregate is the micelle. In particular, there is an emerging consensus that mass transfer and bulk organic phase behaviors-notably the critical phenomenon of phase splitting-are impacted by the effects of micellar-like aggregates in water-in-oil microemulsions. In this study, we elucidate the microscopic structures and mesoscopic architectures of metal-, water-, and acid-loaded organic phases using a combination of X-ray and neutron experimentation as well as density functional theory and molecular dynamics simulations. The key conclusion is that the transfer of metal ions between an aqueous phase and an organic one involves the formation of small mononuclear clusters typical of metal-ligand coordination chemistry, at one extreme, in the organic phase, and their aggregation to multinuclear primary clusters that self-assemble to form even larger superclusters typical of supramolecular chemistry, at the other. Our metrical results add an orthogonal perspective to the energetics-based view of phase splitting in chemical separations known as the micellar model-founded upon the interpretation of small-angle neutron scattering data-with respect to a more general phase-space (gas-liquid) model of soft matter self-assembly and particle growth. The structure hierarchy observed in the aggregation of our quinary (zirconium nitrate-nitric acid-water-tri-n-butyl phosphate-n-octane) system is relevant to understanding solution phase transitions, in general, and the function of engineered fluids with metalloamphiphiles, in particular, for mass transfer applications, such as demixing in separation and synthesis in catalysis science.

A Novel Microemulsion Phase Transition: Toward the Elucidation of Third-Phase Formation in Spent Nuclear Fuel Reprocessing.

We present evidence that the transition between organic and third phases, which can be observed in the plutonium uranium reduction extraction (PUREX) process at high metal loading, is an unusual transition between two isotropic bicontinuous microemulsion phases. As this system contains so many components, however, we have been seeking first to investigate the properties of a simpler system, namely, the related metal-free, quaternary water/n-dodecane/nitric acid/tributyl phosphate (TBP) system. This quaternary system has been shown to exhibit, under appropriate conditions, three coexisting phases: a light organic phase, an aqueous phase, and the so-called third phase. In the current work, we focused on the coexistence of the light organic phase with the third phase. Using Gibbs ensemble Monte Carlo (GEMC) simulations, we found coexistence of a phase rich in nitric acid and dilute in n-dodecane (the third phase) with a phase more dilute in nitric acid but rich in n-dodecane (the light organic phase). The compositions and densities of these two coexisting phases determined using the simulations were in good agreement with those determined experimentally. Because such systems are generally dense and the molecules involved are not simple, the particle exchange rate in their GEMC simulations can be rather low. To test whether a system having a composition between those of the observed third and organic phases is indeed unstable with respect to phase separation, we used the Bennett acceptance ratio method to calculate the Gibbs energies of the homogeneous phase and the weighted average of the two coexisting phases, where the compositions of these phases were taken both from experimental results and from the results of the GEMC simulations. Both demixed states were determined to have statistically significant lower Gibbs energies than the uniform, mixed phase, providing confirmation that the GEMC simulations correctly predicted the phase separation. Snapshots from the simulations and a cluster analysis of the organic and third phases revealed structures akin to bicontinuous microemulsion phases, with the polar species residing within a mesh and with the surface of the mesh formed by amphiphilic TBP molecules. The nonpolar n-dodecane molecules were observed in these snapshots to be outside this mesh. The only large-scale structural differences observed between the two phases were the dimensions of the mesh. Evidence for the correctness of these structures was provided by the results of small-angle X-ray scattering (SAXS) studies, where the profiles obtained for both the organic and third phases agreed well with those calculated from simulations. Finally, we looked at the microscopic structures of the two phases. In the organic phase, the basic motif was observed to be one nitric acid molecule hydrogen-bonded to a TBP molecule. In the third phase, the most common structure was that of the hydrogen-bonded TBP-HNO3-HNO3 chain. A cluster analysis provided evidence for TBP forming an extended, connected network in both phases. Studies of the effects of metal ions on these systems will be presented elsewhere.

Macroscopic chiral symmetry breaking in monolayers of achiral nonconvex platelets.

The fabrication of chiral structures using achiral building blocks is a fundamental problem that remains a challenge in materials science. In this work we present a molecular dynamics simulation study of nonconvex polygonal platelets, interacting via soft-repulsive interactions, that are confined in two-dimensional space. These particle models are designed to promote, even at moderate densities, a natural offset displacement between the edges of neighbouring particles. In particular we demonstrate that nonconvex platelets exhibit macroscopic chiral symmetry breaking when the symmetry of the particles equals (or is multiple of) the number of nearest neighbours in the condensed crystalline phase, corresponding to the situation of platelets with 4-, 6-, and 12-fold symmetries.

Comparative Molecular Dynamics Study on Tri-n-butyl Phosphate in Organic and Aqueous Environments and Its Relevance to Nuclear Extraction Processes.

A refined model for tri-n-butyl phosphate (TBP), which uses a new set of partial charges generated from our ab initio density functional theory calculations, has been proposed in this study. Molecular dynamics simulations are conducted to determine the thermodynamic properties, transport properties, and the microscopic structures of liquid TBP, TBP/water mixtures, and TBP/n-alkane mixtures. These results are compared with those obtained from four other TBP models, previously described in the literature. We conclude that our refined TBP model appears to be the only TBP model from this set that, with reasonable accuracy, can simultaneously predict the properties of TBP in bulk TBP, in organic diluents, and in aqueous solution. The other models only work well for two of the three systems mentioned above. This new TBP model is thus appropriate for the simulation of liquid-liquid extraction systems in the nuclear extraction process, where one needs to simultaneously model TBP in both aqueous and organic phases. It is also promising for the investigation of the microscopic structure of the organic phase in these processes and for the characterization of third-phase formation, where TBP again interacts simultaneously with both polar and nonpolar molecules. Because the proposed TBP model uses OPLS-2005 Lennard-Jones parameters, it may be used with confidence to model mixtures of TBP with other species whose parameters are given by the OPLS-2005 force field.

Note: A dual temperature closed loop batch reactor for determining the partitioning of trace gases within CO2-water systems.

An experimental approach is presented which can be used to determine partitioning of trace gases within CO2-water systems. The key advantages of this system are (1) The system can be isolated with no external exchange, making it ideal for experiments with conservative tracers. (2) Both phases can be sampled concurrently to give an accurate composition at each phase at any given time. (3) Use of a lower temperature flow loop outside of the reactor removes contamination and facilitates sampling. (4) Rapid equilibration at given pressure/temperature conditions is significantly aided by stirring and circulating the water phase using a magnetic stirrer and high-pressure liquid chromatography pump, respectively.

Thermodynamics of the self-assembly of non-ionic chromonic molecules using atomistic simulations. The case of TP6EO2M in aqueous solution.

Atomistic molecular dynamic simulations have been performed for the non-ionic chromonic liquid crystal 2,3,6,7,10,11-hexa-(1,4,7-trioxa-octyl)-triphenylene (TP6EO2M) in aqueous solution. TP6EO2M molecules consist of a central poly-aromatic core (a triphenylene ring) functionalized by six hydrophilic ethyleneoxy (EO) chains, and have a strong tendency to aggregate face-to-face into stacks even in very dilute solution. We have studied self-assembly of the molecules in the low concentration range corresponding to an isotropic solution of aggregates, using two force fields GAFF and OPLS. Our results reveal that the GAFF force field, even though it was successfully used previously for modelling of ionic chromonics, overestimates the attraction of TP6EO2M molecules in water. This results in an aggregation free energy which is too high, a reduced hydration of EO chains and, therefore, molecular self-assembly into compact disordered clusters instead of stacks. In contrast, use of the OPLS force field, leads to self-assembly into ordered stacks in agreement with earlier experimental studies of triphenylene-based chromonics. The free energy of association follows a "quasi-isodesmic" pattern, where the binding free energy of two molecules to form a dimer is of the order of 2.5 RT larger than the corresponding energy of addition of a molecule into a stack. The obtained value for the binding free energy, ΔG=-12 RT, is found to be in line with the published values for typical ionic chromonics (-7 to -12 RT), and agrees reasonably well with the experimental results for this system. The calculated interlayer distance between the molecules in a stack is 0.37 nm, which is at the top of the range found for typical chromonics (0.33-0.37 nm). We suggest that the relatively large layer spacing can be attributed to the repulsion between EO side chains.

Self-assembly and mesophase formation in a non-ionic chromonic liquid crystal system: insights from dissipative particle dynamics simulations.

Results are presented from a dissipative particle dynamics (DPD) simulation of a model non-ionic chromonic system, TP6EO2M, composed of a poly(ethylene glycol) functionalised aromatic (triphenylene) core. The simulations demonstrate self-assembly of chromonic molecules to form single molecule stacks in solution at low concentrations, the formation of a nematic mesophase at higher concentrations and a columnar phase in the more concentrated regime. The simulation model used allows very large system sizes, of many thousands of particles, to be studied. This provides, for the first time, an opportunity to study chromonic phase behaviour by simulation without severe restrictions imposed by system size. In the low concentration limit, the simulations demonstrate approximate isodesmic association from which a binding energy can be obtained, allowing the simulations to be tuned to reproduce the behaviour of the real experimental system.

Screening properties of Gaussian electrolyte models, with application to dissipative particle dynamics.

We investigate the screening properties of Gaussian charge models of electrolyte solutions by analysing the asymptotic behaviour of the pair correlation functions. We use a combination of Monte Carlo simulations with the hyper-netted chain integral equation closure, and the random phase approximation, to establish the conditions under which a screening length is well defined and the extent to which it matches the expected Debye length. For practical applications, for example, in dissipative particle dynamics, we are able to summarise our results in succinct rules-of-thumb which can be used for mesoscale modeling of electrolyte solutions. We thereby establish a solid foundation for future work, such as the systematic incorporation of specific ion effects.

Phase behaviour and the random phase approximation for ultrasoft restricted primitive models.

Phase separation of the ultrasoft restricted primitive model (URPM) with gaussian charges is re-investigated in the random phase approximation (RPA)--the "Level A" approximation discussed by Nikoubashman, Hansen, and Kahl [J. Chem. Phys. 137, 094905 (2012)]. We find that the RPA predicts a region of low temperature vapour-liquid coexistence, with a critical density much lower than that observed in either simulations or more refined approximations (we also remark that the RPA critical point for a related model with Bessel charges can be solved analytically). This observation suggests that the hierarchy of approximations introduced by Nikoubashman et al. should be analogous to those introduced by Fisher and Levin for the restricted primitive model [Phys. Rev. Lett. 71, 3826 (1993)], which makes the inability of these approximations to capture the observed URPM phase behaviour even more worthy of investigation.

High-level virial theory of hard spheroids.

We present a method for calculating high-order virial expansions of the isotropic-nematic phase transition which we apply here to hard spheroids. Studying a range of aspect ratios, for both oblate and prolate particles, we obtain equations of state, coexistence densities, and nematic order parameters, using expansions truncated at up to eighth virial level. For particles of large aspect ratios our results show rapid convergence, with truncation at sixth order sufficient to give excellent agreement with simulation data. For more spherical particles the convergence is less rapid, with results for up to eighth-order theory approaching but still not reaching simulation data. Our results indicate that high-order viral expansions are better suited to predicting equations of state than coexistence densities. We also test the validity of using the Onsager trial function to approximate the orientational distribution function, finding only small errors when making this approximation.

Thermodynamic stability of the cubatic phase of hard cut spheres evaluated by expanded ensemble simulations.

The system of hard cut spheres (disk-shaped particles formed by symmetrically truncating the end caps of a sphere) exhibits an intriguing "cubatic" phase with cubic orientational symmetry. However, it is unclear whether this phase is metastable with respect to the columnar phase. We attempt to provide an answer to this question by carrying out free energy calculations by the expanded ensemble Monte Carlo method. We conclude that there may be a very small region of cubatic stability in the vicinity of the isotropic-cubatic phase transition, but that that transition would need to be determined more accurately to obtain a definitive answer. We also comment on the efficacy of the expanded ensemble method for these kinds of calculations.

Theory and computer simulation for the cubatic phase of cut spheres.

The phase behavior of a system of hard-cut spheres has been studied using a high-order virial theory and by Monte Carlo simulation. The cut-sphere particles are disks of thickness L formed by symmetrically truncating the end caps of a sphere of diameter D . The virial theory predicts a stable nematic phase for aspect ratio LD=0.1 and a stable cubatic phase for LD=0.15-0.3 . The virial series converges rapidly on the equation of state of the isotropic and nematic phases, while for the cubatic phase the convergence is slower, but still gives good agreement with the simulation at high order. It is found that a high-order expansion (up to B8 ) is required to predict a stable cubatic phase for LD> or =0.15 , indicating the importance of many-body interactions in stabilizing this phase. Previous simulation work on this system has focused on aspect ratios LD=0.1 , 0.2, and 0.3. We expand this to include also LD=0.15 and 0.25, and we introduce a fourth-rank tensor to measure cubatic ordering. We have applied a multiparticle move which dramatically speeds the attainment of equilibrium in the nematic phase and therefore is of great benefit in the study of the isotropic-nematic phase transition. In agreement with the theory, our simulations confirm the stability of the nematic phase for LD=0.1 and the stability of the cubatic phase over the nematic for LD=0.15-0.3 . There is, however, some doubt about the stability of the cubatic phase with respect to the columnar. We have shown that the cubatic phase found on compression at LD=0.1 is definitely metastable, but the results for LD=0.2 were less conclusive.

Structure and stability of isotropic states of hard platelet fluids.

We study the thermodynamics and the pair structure of hard, infinitely thin, circular platelets in the isotropic phase. Monte Carlo simulation results indicate a rich spatial structure of the spherical expansion components of the direct correlation function, including nonmonotonical variation of some of the components with density. Integral equation theory is shown to reproduce the main features observed in simulations. The hypernetted chain closure, as well as its extended versions that include the bridge function up to second and third order in density, perform better than both the Percus-Yevick closure and Verlet bridge function approximation. Using a recent fundamental measure density functional theory, an analytic expression for the direct correlation function is obtained as the sum of the Mayer bond and a term proportional to the density and the intersection length of two platelets. This is shown to give a reasonable estimate of the structure found in simulations, but to fail to capture the nonmonotonic variation with density. We also carry out a density functional stability analysis of the isotropic phase with respect to nematic ordering and show that the limiting density is consistent with that where the Kerr coefficient vanishes. As a reference system, we compare to simulation results for hard oblate spheroids with small, but nonzero elongations, demonstrating that the case of vanishingly thin platelets is approached smoothly.

Structure of molecular liquids: hard rod-disk mixtures.

The structure of hard rod-disk mixtures is studied using Monte Carlo simulations and integral equation theory, for a range of densities in the isotropic phase. By extension of methods used in single component fluids, the pair correlation functions of the molecules are calculated and comparisons between simulation and integral equation theory, using a number of different closure relations, are made. Comparison is also made for thermodynamic data and phase behavior.

Structure of molecular liquids: closure relations for hard spheroids.

We present the results of Monte Carlo simulations of hard spheroids of revolution of different elongations. Both prolate and oblate shapes are examined. A systematic study of the bridge function b(1,2), and direct comparison with the indirect correlation function gamma(1,2)=h(1,2)-c(1,2) at densities spanning the isotropic fluid range, allow us to evaluate the accuracy of various proposed closure relations for integral equations.

Structure of molecular liquids: cavity and bridge functions of the hard spheroid fluid.

We present methodologies for calculating the direct correlation function c(1,2), the cavity function y(1,2), and the bridge function b(1,2), for molecular liquids, from Monte Carlo simulations. As an example we present results for the isotropic hard spheroid fluid with elongation e = 3. The simulation data are compared with the results from integral equation theory. In particular, we solve the Percus-Yevick and hypernetted chain equations. In addition, we calculate the first two terms in the virial expansion of the bridge function and incorporate this into the closure. At low densities, the bridge functions calculated by theory and from simulation are in good agreement, lending support to the correctness of our numerical procedures. At higher densities, the hypernetted chain results are brought into closer agreement with simulation by incorporating the approximate bridge function, but significant discrepancies remain.

Exploring reaction pathways with transition path and umbrella sampling: application to methyl maltoside.

The transition path sampling (TPS) method is a powerful approach to study chemical reactions or transitional properties on complex potential energy landscapes. One of the main advantages of the method over potential of mean force methods is that reaction rates can be directly accessed without knowledge of the exact reaction coordinate. We have investigated the complementary nature of these two differing approaches, comparing transition path sampling with the weighted histogram analysis method to study a conformational change in a small model system. In this case study, the transition paths for a transition between two rotational conformers of a model disaccharide molecule, methyl beta-D-maltoside, were compared with a free energy surface constrained by the two commonly used glycosidic (phi,psi) torsional angles. The TPS method revealed a reaction channel that was not apparent from the potential of mean force method, and the suitability of phi and psi as reaction coordinates to describe the isomerization in vacuo was confirmed by examination of the transition path ensemble. Using both transition state theory and transition path sampling methods, the transition rate was estimated. We have estimated a characteristic time between transitions of approximately 160 ns for this rare isomerization event between the two conformations of the carbohydrate. We conclude that transition path sampling can extract subtle information about the dynamics not apparent from the potential of mean force method. However, in calculating the reaction rate, the transition path sampling method required 27.5 times the computational effort than was needed by the potential of mean force method.