Semi-grand Canonical Monte Carlo Simulations Research Articles

Miscibility and exchange chemical potential of ring polymers in symmetric ring-ring blends.

Generally, differences of polymer topologies may affect polymer miscibility even with the same repeated units. In this study, the topological effect of ring polymers on miscibility was investigated by comparing symmetric ring-ring and linear-linear polymer blends. To elucidate the topological effect of ring polymers on mixing free energy, the exchange chemical potential of binary blends was numerically evaluated as a function of composition ϕ by performing semi-grand canonical Monte Carlo and molecular dynamics simulations of a bead-spring model. For ring-ring blends, an effective miscibility parameter was evaluated by comparing the exchange chemical potential with that of the Flory-Huggins model for linear-linear polymer blends. It was confirmed that in the mixed states satisfying χN > 0, ring-ring blends are more miscible and stable than the linear-linear blends with the same molecular weight. Furthermore, we investigated finite molecular weight dependence on the miscibility parameter, which reflected the statistical probability of interchain interactions in the blends. The simulation results revealed that the molecular weight dependence on the miscibility parameter was smaller in ring-ring blends. The effect of the ring polymers on miscibility was verified to be consistent with the change in the interchain radial distribution function. In ring-ring blends, it was indicated that the topology affected miscibility by reducing the effect of the direct interaction between the components of the blends.

Polyelectrolyte-nanoparticle mutual charge regulation and its influence on their complexation

The complexation between a charged nanoparticle and an oppositely charged linear polyelectrolyte is studied using coarse-grained Semi-Grand Canonical Monte Carlo simulations. This is done for systems where one or both of the macromolecules are titratable, i.e. they have a pH-dependent charge. The presence of charge fluctuation is found to have a significant impact on the ionization of both macromolecules, the conformational properties of the adsorbed polyelectrolyte and the charge of the resulting complex. When both macromolecules are titratable, we find that the impact is larger when the absolute difference in the pKa values of the nanoparticle and polyelectrolyte increases. The physicochemical concepts here discussed are found to be in good agreement with previous experimental observations, providing a fundamental basis to rationalize these experiments.

Semi-grand canonical Monte Carlo simulation of the acrolein induced surface segregation and aggregation of AgPd with machine learning surrogate models.

The single atom alloy of AgPd has been found to be a promising catalyst for the selective hydrogenation of acrolein. It is also known that the formation of Pd islands on the surface will greatly reduce the selectivity of the reaction. As a result, the surface segregation and aggregation of Pd on the AgPd surface under reaction conditions of selective hydrogenation of acrolein are of great interest. In this work, we lay out a workflow that can predict the surface segregation and aggregation of Pd on a FCC(111) AgPd surface with and without the presence of acrolein. We use machine learning surrogate models to predict the AgPd bulk energy, AgPd slab energy, and acrolein adsorption energy on AgPd slabs. Then, we use the semi-grand canonical Monte Carlo simulation to predict the surface segregation and aggregation under different bulk Pd concentrations. Under vacuum conditions, our method predicts that only trace amount of Pd will exist on the surface at Pd bulk concentrations less than 20%. However, with the presence of acrolein, Pd will start to aggregate as dimers on the surface at Pd bulk concentrations as low as 6.5%.

Influence of macromolecular crowding on the charge regulation of intrinsically disordered proteins.

In this work we study the coupling between ionization and conformational properties of two IDPs, histatin-5 and β-amyloid 42, in the presence of neutral and charged crowders. The latter is modeled to resemble bovine serum albumin (BSA). With this aim, semi-grand canonical Monte Carlo simulations are performed, so that the IDP charge is a dynamic property, undergoing protonation/deprotonation processes. Both ionization properties (global and specific amino acid charge and binding capacitance) and radius of gyration are analyzed in a large range of pH values and salt concentrations. Without crowder agents, the titration curve of histatin-5, a polycation, is salt-dependent while that of β-amyloid 42, a polyampholyte, is almost unaffected. The salt concentration is found to be particularly relevant at pH values where the protein binding capacitance (directly linked with charge fluctuation) is larger. Upon addition of neutral crowders, charge regulation is observed in histatin-5, while for β-amyloid 42 this effect is very small. The main mechanism for charge regulation is found to be the effective increase in the ionic strength due to the excluded volume. In the presence of charged crowders, a significant increase in the charge of both IDPs is observed in almost all the pH range. In this case, the IDP charge is altered not only by the increase in the effective ionic strength but also by its direct electrostatic interaction with the charged crowders.

A semi-GCMC simulation study of electrolytic capacitors with adsorbed titrating peptides.

We use semi-grand canonical Monte Carlo simulations to study an electrolytic capacitor with an adsorbed peptide on the electrode surfaces. Only homogeneous peptides are considered, consisting of only a single residue type. We find that the classical double-hump camel-shaped differential capacitance in such systems is augmented by the addition of a third peak, due to the capacitance contribution of the peptide, essentially superimposed on the salt contribution. This mechanistic picture is justified using a simple mean-field analysis. We find that the position of this third peak can be tuned to various surface potential values by adjusting the ambient pH of the electrolyte solution. We investigate the effect of changing the residue type and the concentration of the adsorbed peptide and of the supporting electrolyte. Varying the residue species and pH allows one to modify the capacitance profile as a function of surface potential, facilitating the design of varying discharging patterns for the capacitor.

Interface kinetics of rapid solidification of binary alloys by atomistic simulations: Application to Ti-Ni alloys

Using molecular dynamics and Monte Carlo simulations we investigate the non-equilibrium interfacial kinetics during rapid solidification of Ti-Ni alloys. According to the existing theories, the kinetic coefficient is related, via an analytical expression, to the equilibrium and non-equilibrium solute concentration profiles across the crystal-melt interface, the interface temperature and velocity, and the drag coefficient. The kinetic coefficient was obtained by deriving these properties from specifically designed atomistic simulations and then fitting using the analytical expression. The results show that the kinetic coefficient is only weakly anisotropic and increases with increasing temperature. The velocity-dependent partition coefficient, as described by two solute trapping models, the continuous growth and the local non-equilibrium models, were fitted to the molecular dynamics simulation results. In addition, molecular dynamics and semi-grand canonical Monte Carlo simulations suggest that complete solute trapping might occur only under certain conditions. This can be explained by the fact that the maximum achievable effective free energy driving force for solidification, which defines the complete solute trapping, is limited by the chemical potential and free energy profiles for the Ti-Ni alloy. The investigations, using molecular dynamics simulations, of the dependence of crystal-melt interface width on solidification velocity show, at high velocities, similar trend to that predicted by the hyperbolic phase field model which is suitable for studies of rapid solidification of alloys.

Genetic algorithm-guided deep learning of grain boundary diagrams: Addressing the challenge of five degrees of freedom

Grain boundaries (GBs) often control the processing and properties of polycrystalline materials. Here, potentially transformative research is represented by constructing GB property diagrams as functions of temperature and bulk composition, also called “complexion diagrams,” as a general materials science tool on par with phase diagrams. However, a GB has five macroscopic (crystallographic) degrees of freedom (DOFs). It is essentially a “mission impossible” to construct property diagrams for GBs as a function of five DOFs by either experiments or modeling. Herein, we combine isobaric semi-grand canonical ensemble hybrid Monte Carlo and molecular dynamics (hybrid MC/MD) simulations with a genetic algorithm (GA) and deep neural network (DNN) models to tackle this grand challenge. The DNN prediction is ∼108 faster than atomistic simulations, thereby enabling the construction of the property diagrams for millions of distinctly different GBs of five DOFs. Notably, excellent prediction accuracies have been achieved for not only symmetric-tilt and twist GBs, but also asymmetric-tilt and mixed tilt-twist GBs; the latter are more complex and much less understood, but they are ubiquitous and often limit the performance properties of real polycrystals as the weak links. The data-driven prediction of GB properties as function of temperature, bulk composition, and five crystallographic DOFs (i.e., in a 7D space) opens a new paradigm.

Semi-grand canonical Monte Carlo simulation for derivation of thermodynamic properties of binary alloy

Semi-grand canonical Monte Carlo (SGCMC) simulations are performed to derive thermodynamic properties of binary alloy from atomistic-based simulations. Particularly, solidus and liquidus compositions are directly derived for Fe-Cr alloy described by two different EAM potentials. Although the SGCMC simulation can derive relationship between the free energy and composition at any temperature straightforwardly, partial phase diagram obtained from SGCMC simulations strongly depends on the choice of interatomic potential.

Structure of lateral heterogeneities in a coarse-grained model for multicomponent membranes.

We study the lateral domain structure in a coarse-grained molecular model for multicomponent lipid bilayers by semi-grandcanonical Monte Carlo simulations. The membranes are filled with liquid ordered (lo) domains surrounded by a liquid disordered (ld) matrix. Depending on the membrane composition and temperature, we identify different morphological regimes: one regime (I) where the lo domains are small and relatively compact, and two regimes (II, II') where they are larger and often interconnected. In the latter two regimes, the ld matrix forms a network of disordered trenches separating the lo domains, with a relatively high content of interdigitated line defects. Since such defects are also a structural element of the modulated ripple phase in one component membranes, we argue that the regimes II, II' may be amorphous equivalents of the ripple phase in multicomponent membranes. We also analyze the local structure and provide evidence that the domains in regime I are stabilized by a monolayer curvature mechanism postulated in earlier work [S. Meinhardt et al., PNAS, 2013, 110, 4476].

Atomistic and mesoscale simulation of sodium and potassium adsorption in cement paste.

An atomistic and mesoscopic assessment of the effect of alkali uptake in cement paste is performed. Semi-grand canonical Monte Carlo simulations indicate that Na and K not only adsorb at the pore surface of calcium silicate hydrates (C-S-H) but also adsorb in the C-S-H hydrated interlayer up to concentrations of the order of 0.05 and 0.1 mol/kg, respectively. Sorption of alkali is favored as the Ca/Si ratio of C-S-H is reduced. Long timescale simulations using the Activation Relaxation Technique indicate that characteristic diffusion times of Na and K in the C-S-H interlayer are of the order of a few hours. At the level of individual grains, Na and K adsorption leads to a reduction of roughly 5% of the elastic moduli and to volume expansion of about 0.25%. Simulations using the so-called primitive model indicate that adsorption of alkali ions at the pore surface can reduce the binding between C-S-H grains by up to 6%. Using a mesoscopic model of cement paste, the combination of individual grain swelling and changes in inter-granular cohesion was estimated to lead to overall expansive pressures of up to 4 MPa-and typically of less than 1 MPa-for typical alkali concentrations observed at the proximity of gel veins caused by the alkali-silica reaction.

Unusual temperature dependence of the solid-liquid interfacial free energy in the Cu-Zr system

The solid-liquid portion of the equilibrium phase diagram at the Cu rich side of a Finnis-Sinclair description of Cu-Zr was determined by a combination of molecular dynamics and semi-grand canonical Monte Carlo simulations and the model system was found to be in very good agreement with experiment. Molecular dynamics simulations and the capillary fluctuation method were used to compute the solid-liquid interfacial free energy, γ, and its anisotropy as a function of undercooling. Unlike the case of pure metals and in contrast to most previous atomistic studies in alloy systems, γ was found to increase with decreasing temperature and the increase is quite pronounced. To verify the unusual results, the excess energy and number of Zr atoms were computed for the [1 0 0] interface and from the excess thermodynamic quantities γ was obtained through a thermodynamic integration. The results of the two approaches for computing γ are in agreement to within the statistical uncertainties. The implications of the temperature behavior of γ is discussed in terms of crystal nucleation from the melt and glass formation.

Coupling of Charge Regulation and Conformational Equilibria in Linear Weak Polyelectrolytes: Treatment of Long-Range Interactions via Effective Short-Ranged and pH-Dependent Interaction Parameters.

The classical Rotational Isomeric State (RIS) model, originally proposed by Flory, has been used to rationalize a wide range of physicochemical properties of neutral polymers. However, many weak polyelectrolytes of interest are able to regulate their charge depending on the conformational state of the bonds. Recently, it has been shown that the RIS model can be coupled with the Site Binding (SB) model, for which the ionizable sites can adopt two states: protonated or deprotonated. The resulting combined scheme, the SBRIS model, allows for analyzing ionization and conformational equilibria on the same foot. In the present work, this approach is extended to include pH-dependent electrostatic Long-Range (LR) interactions, ubiquitous in weak polyelectrolytes at moderate and low ionic strengths. With this aim, the original LR interactions are taken into account by defining effective Short-Range (SR) and pH-dependent parameters, such as effective microscopic protonation constants and rotational bond energies. The new parameters are systematically calculated using variational methods. The machinery of statistical mechanics for SR interactions, including the powerful and fast transfer matrix methods, can then be applied. The resulting technique, which we will refer to as the Local Effective Interaction Parameters (LEIP) method, is illustrated with a minimal model of a flexible linear polyelectrolyte containing only one type of rotating bond. LEIP reproduces very well the pH dependence of the degree of protonation and bond probabilities obtained by semi-grand canonical Monte Carlo simulations, where LR interactions are explicitly taken into account. The reduction in the computational time in several orders of magnitude suggests that the LEIP technique could be useful in a range of areas involving linear weak polyelectrolytes, allowing direct fitting of the relevant physical parameters to the experimental quantities.

Quantum–Continuum Modeling of Pd-Au(111) Surface Alloys Under Electrochemical Conditions

Bimetallic surface alloys often exhibit enhanced electrocatalytic performance due to ligand, ensemble, and bifunctional effects that arise from the presence of different transition metals along the catalyst surface. Segregation, dealloying, and poisoning events are common modes of catalyst deactivation that are difficult to suppress, limiting the durability of bimetallic catalysts in harsh environments. Efforts to estimate the stability of these novel catalyst formulations under operational conditions via first principles techniques is challenging due to the critical influence of solvation, surface electrification, and finite temperature effects. In this talk, we present a quantum-continuum model of the electrode-solution interface [Weitzner & Dabo, npj Comp. Mater. 3, 1 (2017)] that directly accounts for the influence of the electrochemical environment and pair this model with a recently developed voltage-sensitive cluster expansion technique to perform semi-grand canonical Monte Carlo simulations of the interface [Weitzner & Dabo, Phys. Rev. B (in press)]. We apply this method to study the voltage-dependent stability of Pd-Au(111) surface alloys under realistic electrochemical conditions and discuss the role of environmental effects on the durability of the catalyst.

Soluble Oligomeric Nucleants: Simulations of Chain Length, Binding Strength, and Volume Fraction Effects.

Recent theories and simulations suggest that molecular additives can bind to the surfaces of nuclei, lower the surface energy, and accelerate nucleation. Experiments have shown that oligomeric and polymeric additives can also modify nucleation rates of proteins, ice, and minerals; however, general design principles for oligomeric or polymeric promoters do not yet exist. Here we investigate oligomeric additives for which each segment of the oligomer can bind to surfaces of nuclei. We use semigrand canonical Monte Carlo simulations in a Potts lattice gas model to study the effects of oligomer chain length, volume fraction, and binding strength. We find that increasing each of those parameters lowers the nucleation barrier. At extremely low oligomer concentrations, the nucleation kinetics can be modeled as though each oligomer is a heterogeneous nucleation site in solution.

Structure and dynamics of binary liquid mixtures near their continuous demixing transitions

The dynamic and static critical behavior of a family of binary Lennard-Jones liquid mixtures, close to their continuous demixing points (belonging to the so-called model H' dynamic universality class), are studied computationally by combining semi-grand canonical Monte Carlo simulations and large-scale molecular dynamics (MD) simulations, accelerated by graphic processing units (GPU). The symmetric binary liquid mixtures considered cover a variety of densities, a wide range of compressibilities, and various interactions between the unlike particles. The static quantities studied here encompass the bulk phase diagram (including both the binodal and the λ-line), the correlation length, and the concentration susceptibility, of the finite-sized systems above the bulk critical temperature Tc, the compressibility and the pressure at Tc. Concerning the collective transport properties, we focus on the Onsager coefficient and the shear viscosity. The critical power-law singularities of these quantities are analyzed in the mixed phase (above Tc) and non-universal critical amplitudes are extracted. Two universal amplitude ratios are calculated. The first one involves static amplitudes only and agrees well with the expectations for the three-dimensional Ising universality class. The second ratio includes also dynamic critical amplitudes and is related to the Einstein-Kawasaki relation for the interdiffusion constant. Precise estimates of this amplitude ratio are difficult to obtain from MD simulations, but within the error bars our results are compatible with theoretical predictions and experimental values for model H'. Evidence is reported for an inverse proportionality of the pressure and the isothermal compressibility at the demixing transition, upon varying either the number density or the repulsion strength between unlike particles.

Grain-boundary layering transitions and phonon engineering

We employ semi-grand canonical Monte Carlo simulation to investigate layering transitions at grain boundaries in a prototypical binary alloy. We demonstrate the existence of such transitions among various interfacial states and examine the role of elastic fields in dictating state equilibria. The results of these studies are summarized in the form of diagrams that highlight interfacial state coexistence in this system. Finally, we examine the impact of layering transitions on the phononic properties of the system, as given by the specific heat and, by extension, the thermal conductivity. Thus, it is suggested that by inducing interfacial layering transitions via changes in temperature or pressure, one can thereby engineer thermodynamic and transport properties in materials.

Calculation of open and closed system elastic coefficients for multicomponent solids

Thermodynamic equilibrium in multicomponent solids subject to mechanical stresses is a complex nonlinear problem whose exact solution requires extensive computations. A few decades ago, Larch\'e and Cahn proposed a linearized solution of the mechanochemical equilibrium problem by introducing the concept of open system elastic coefficients [Acta Metall. 21, 1051 (1973)]. Using the Ni-Al solid solution as a model system, we demonstrate that open system elastic coefficients can be readily computed by semigrand canonical Monte Carlo simulations in conjunction with the shape fluctuation approach. Such coefficients can be derived from a single simulation run, together with other thermodynamic properties needed for prediction of compositional fields in solid solutions containing defects. The proposed calculation approach enables streamlined solutions of mechanochemical equilibrium problems in complex alloys. Second order corrections to the linear theory are extended to multicomponent systems.

Hierarchical self-assembly: Self-organized nanostructures in a nematically ordered matrix of self-assembled polymeric chains.

We report many different nanostructures which are formed when model nanoparticles of different sizes (diameter σn) are allowed to aggregate in a background matrix of semiflexible self-assembled polymeric wormlike micellar chains. The different nanostructures are formed by the dynamical arrest of phase-separating mixtures of micellar monomers and nanoparticles. The different morphologies obtained are the result of an interplay of the available free volume, the elastic energy of deformation of polymers, the density (chemical potential) of the nanoparticles in the polymer matrix, and, of course, the ratio of the size of self-assembling nanoparticles and self-avoidance diameter of polymeric chains. We have used a hybrid semi-grand-canonical Monte Carlo simulation scheme to obtain the (nonequilibrium) phase diagram of the self-assembled nanostructures. We observe rodlike structures of nanoparticles which get self-assembled in the gaps between the nematically ordered chains, as well as percolating gel-like network of conjoined nanotubes. We also find a totally unexpected interlocked crystalline phase of nanoparticles and monomers, in which each crystal plane of nanoparticles is separated by planes of perfectly organized polymer chains. We identified the condition which leads to such interlocked crystal structure. We suggest experimental possibilities of how the results presented in this paper could be used to obtain different nanostructures in the laboratory.

Semigrand Canonical and Kinetic Monte Carlo simulations of binary B2-ordered nano-films with triple defects

Equilibrium atomic configurations and the kinetics of “order–order” and surface segregation processes in B2-ordering stoichiometric A-50 at.%B binary thin films are investigated by means of Semigrand Canonical Monte Carlo (SGCMC) and Kinetic Monte Carlo (KMC) simulations. The (100)-oriented films are modeled with an Ising-type Hamiltonian with previously evaluated pair-interaction energy parameters yielding the effect of “triple-defect disordering”. The SGCMC simulations provide equilibrium vacancy concentration and atomic configuration in the films with B-atom termination of both free surfaces achieved at high temperatures by the generation of an antiphase boundary. Despite strong vacancy surface segregation, the thermodynamic activation energy for their formation inside the films is the same as in the bulk material. KMC simulations implemented with the SGCMC-determined equilibrium vacancy concentration reveal very slow relaxation of the films towards equilibrium configuration. The B-termination of the (100) free surfaces is produced by A-atom diffusion inwards into the films mediated by vacancies segregating on surfaces.

Phase separation of binary nonadditive hard sphere fluid mixture confined in random porous media

I analyze the fluid-fluid phase separation of nonadditive hard sphere fluid mixture absorbed in random porous media. An equation of state is derived by using the perturbation theory to this complex system with quenched disorders. The results of this theory are in good agreement with those obtained from semi-grand canonical ensemble Monte Carlo simulations. The contact value of the fluid-fluid radial distribution functions of the reference which is the key point of the perturbation process is derived as well, the comparison against Monte Carlo simulations shows that it has an excellent accuracy.