Generalized Energy-Conserving Dissipative Particle Dynamics with Mass Transfer. Part 2: Applications and Demonstrations.

We present the second part of a two-part paper series intended to address a gap in computational capability for coarse-grain particle modeling and simulation, namely, the simulation of phenomena in which diffusion via mass transfer is a contributing mechanism. In part 1, we presented a formulation of a dissipative particle dynamics method to simulate interparticle mass transfer, termed generalized energy-conserving dissipative particle dynamics with mass transfer (GenDPDE-M). In the GenDPDE-M method, the mass of each mesoparticle remains constant following the interparticle mass exchange. In part 2 of this series, further verification and demonstrations of the GenDPDE-M method are presented for mesoparticles with embedded binary mixtures using the ideal gas (IG) and van der Waals (vdW) equation-of-state (EoS). The targeted readership of part 2 is toward practitioners, where applications and practical considerations for implementing the GenDPDE-M method are presented and discussed, including a numerical discretisztion algorithm for the equations-of-motion. The GenDPDE-M method is verified by reproducing the particle distributions predicted by Monte Carlo simulations for the IG and vdW fluids, along with several demonstrations under both equilibrium and non-equilibrium conditions. GenDPDE-M can be generally applied to multi-component mixtures and to other fundamental EoS, such as the Lennard-Jones or Exponential-6 models, as well as to more advanced EoS models such as Statistical Associating Fluid Theory.

Generalized Energy-Conserving Dissipative Particle Dynamics with Mass Transfer. Part 1: Theoretical Foundation and Algorithm.

An extension of the generalized energy-conserving dissipative particle dynamics method (GenDPDE) that allows mass transfer between mesoparticles via a diffusion process is presented. By considering the concept of the mesoparticles as property carriers, the complexity and flexibility of the GenDPDE framework were enhanced to allow for interparticle mass transfer under isoenergetic conditions, notated here as GenDPDE-M. In the formulation, diffusion is described via the theory of mesoscale irreversible processes based on linear relationships between the fluxes and thermodynamic forces, where their fluctuations are described by Langevin-like equations. The mass exchange between mesoparticles is such that the mass of the mesoparticle remains unchanged after the transfer process and requires additional considerations regarding the coupling with other system properties such as the particle internal energy. The proof-of-concept work presented in this article is the first part of a two-part article series. In Part 1, the development of the GenDPDE-M theoretical framework and the derivation of the algorithm are presented in detail. Part 2 of this article series is targeted for practitioners, where applications, demonstrations, and practical considerations for implementing the GenDPDE-M method are presented and discussed.

Generalized Energy-Conserving Dissipative Particle Dynamics with Reactions.

We present an extension of the generalized energy-conserving dissipative particle dynamics method (J. Bonet Avalos, et al., Phys Chem Chem Phys, 2019, 21, 24891-24911) to include chemical reactivity, denoted GenDPDE-RX. GenDPDE-RX provides a means of simulating chemical reactivity at the micro- and mesoscales, while exploiting the attributes of density- and temperature-dependent many-body force fields, which include improved transferability and scalability compared to two-body pairwise models. The GenDPDE-RX formulation considers intra-particle reactivity via a coarse-grain reactor construct. Extent-of-reaction variables assigned to each coarse-grain particle monitor the temporal evolution of the prescribed reaction mechanisms and kinetics assumed to occur within the particle. Descriptions of the algorithm, equations of motion, and numerical discretization are presented, followed by verification of the GenDPDE-RX method through comparison with reaction kinetics theoretical model predictions. Demonstrations of the GenDPDE-RX method are performed using constant-volume adiabatic heating simulations of three different reaction models, including both reversible and irreversible reactions, as well as multistep reaction mechanisms. The selection of the demonstrations is intended to illustrate the flexibility and generality of the method but is inspired by real material systems that span from fluids to solids. Many-body force fields using analytical forms of the ideal gas, Lennard-Jones, and exponential-6 equations of state are used for demonstration, although application to other forms of equation of states is possible. Finally, the flexibility of the GenDPDE-RX framework is addressed with a brief discussion of other possible adaptations and extensions of the method.

Universal Scaling for the Exit Dynamics of Block Copolymers from Micelles at Short and Long Time Scales.

The correlation function for the exit of poloxamer copolymers from equilibrated micelles is found to show up to four regimes depending on the chain flexibility: an initial fast reorganization, a logarithmic intermediate regime, followed by an exponential intermediate regime, and a final exponential decay. The logarithmic intermediate regime has been observed experimentally and attributed to the polydispersity of the polymer samples. However, we present dynamic single-chain mean-field theory simulations with chains of variable flexibility which show the same logarithmic relaxation but with strictly monodisperse systems. In agreement with our previous studies, we propose that this logarithmic response arises from a degeneracy of energy states of the hydrophobic block in the micelle core. For this to occur, a sufficiently large number of degenerate conformational states are required, which depend on the polymer flexibility and therefore should not be present for rigid polymers. Experimental results for monodisperse polymeric samples claiming the absence of such a logarithmic response may also lack a sufficient number of hydrophobic blocks for the required number of configurational states for this type of response to be seen. The insight gained from analyzing the simulation results allows us to propose a modified Eyring equation capable of reproducing the observed dynamic behavior. On scaling experimental results from different sources and systems according to this equation, we find a unique master curve showing a universal nature of the intermediate regimes: the logarithmic regime together with the secondary exponential decay. The terminal exponential regime at long times proposed by the standard Halperin and Alexander model is beyond the range of the data analyzed in this article. The universality observed suggests an entropic origin of the short-time dynamic response of this class of systems rather than the polydispersity.

Shear-viscosity-independent bulk-viscosity term in smoothed particle hydrodynamics.

An angular momentum conservative pure bulk viscosity term for smoothed particle hydrodynamics (SPH) is proposed in the present paper. This formulation permits independent modeling of shear and bulk viscosities, which is of paramount importance for fluids with large bulk viscosity in situations where sound waves or large Mach numbers are expected. With this aim a dissipative term proportional to the rate of change of the volume is considered at the particle level. The equations of motion are derived from the minimization of a Lagrangian combined with an appropriate dissipation function that depends on this rate of change of particle volume, in analogy with the corresponding entropy production contribution in fluids. Due to the Galilean invariance of the formulation, the new term is shown to exactly conserve linear momentum. Moreover, its invariance under solid-body rotations also ensures the conservation of angular momentum. Two verification cases are proposed: the one-dimensional propagation of a sound pulse and a two-dimensional case, modeling the time decay of an accelerating-decelerating pipe flow. The SPH solutions are compared to exact ones, showing that the newly proposed term behaves indeed as a viscosity associated only with the local expansion-compression of the fluid. In view of these considerations, we conclude that the method presented in this paper allows for setting up a bulk viscosity independently of the shear one and as large as any particular problem may require. At the same time, together with the prescribed momentum conservation to reproduce the Navier-Stokes equation, the new term also keeps the angular momentum conservation required to properly model free interfaces or overall rotations of the bulk fluid.

Generalised dissipative particle dynamics with energy conservation: density- and temperature-dependent potentials.

We present a generalised, energy-conserving dissipative particle dynamics (DPDE) method appropriate for the non-isothermal simulation of particle interaction force fields that are both density- and temperature-dependent. A detailed derivation is formulated in a bottom-up manner by considering the thermodynamics of small systems with the appropriate consideration of the fluctuations. Connected to the local volume is a local density and corresponding local pressure, which is determined from an equation-of-state based force field that depends also on a particle temperature. Compared to the original DPDE method, the formulation of the generalised DPDE method requires a change in the independent variable from the particle internal energy to the particle entropy. As part of the re-formulation, the terms dressed particle entropy and the corresponding dressed particle temperature are introduced, which depict the many-body contributions in the local volume. The generalised DPDE method has similarities to the energy form of the smoothed dissipative particle dynamics method, yet fundamental differences exist, which are described in the manuscript. The basic dynamic equations are presented along with practical considerations for implementing the generalised DPDE method, including a numerical integration scheme based on the Shardlow-like splitting algorithm. Demonstrations and validation tests are performed using analytical equation-of-states for the van der Waals and Lennard-Jones fluids. Particle probability distributions are analysed, where excellent agreement with theoretical estimates is demonstrated. As further validation of the generalised DPDE method, both equilibrium and non-equilibrium simulation scenarios are considered, including adiabatic flash heating response and vapour-liquid phase separation.

Simulation Analysis of the Kinetic Exchange of Copolymer Surfactants in Micelles.

The exchange of surfactants in micelles involves several processes that are difficult to characterize experimentally. Microscopic simulations have the potential to reveal some of the key activities that occur when a surfactant spontaneously exits a micelle. In this work, we present a quantitative analysis of the kinetic exchange process over a large range of time. This study is based on a dynamic version of single-chain mean-field theory using a coarse-grained model for poly(ethylene oxide)-poly(propylene oxide)-poly(ethylene oxide) triblock copolymer systems. The kinetics described in our simulations involves three different regimes. After a fast initial rearrangement of the labeled chains, the system undergoes a logarithmic relaxation, which has been experimentally observed. Contrary to what has been reported in previous analyses, our simulations indicate that this regime is caused by the intrinsic physical behavior of the system and is not due only to the polydispersity of the samples. Finally, the terminal regime is characterized by an exponential decay. The exit rates predicted by our simulations are in good agreement with the values reported experimentally. In addition, we address the sequence of microscopic conformational changes undergone by the surfactants when leaving the micellar aggregates. We found a subtle variation in the radius of gyration of the hydrophobic block, which challenges the image of either a complete collapse or a full stretching commonly accepted in the current theoretical and experimental literature.

Hydrophobic effect as a driving force for host-guest chemistry of a multi-receptor Keplerate-type capsule.

The effectiveness of the interactions between various alkylammonium cations and the well-defined spherical Keplerate-type {Mo132} capsule has been tracked by (1)H DOSY NMR methodology, revealing a strong dependence on the self-diffusion coefficient of the cationic guests balanced between the solvated and the plugging situations. Analysis of the data is fully consistent with a two-site exchange regime involving the 20 independent {Mo9O9} receptors of the capsule. Furthermore, quantitative analysis allowed us to determine the stability constants associated with the plugging process of the pores. Surprisingly, the affinity of the capsule for a series of cationic guests increases continuously with its apolar character, as shown by the significant change of the stability constant from 370 to 6500 for NH4(+) and NEt4(+), respectively. Such observations, supported by the thermodynamic parameters, evidence that the major factor dictating selectivity in the trapping process is the so-called "hydrophobic effect". Computational studies, using molecular dynamics simulations, have been carried out in conjunction with the experiments. Analysis of the radial distribution functions g(r) reveals that NH4(+) and NMe4(+) ions behave differently in the vicinity of the capsule. The NH4(+) ions do not exhibit well-defined distributions when in close vicinity. In contrast, the NMe4(+) ions displayed sharp distributions related to different scenarios, such as firmly trapped or labile guest facing the {Mo9O9} pores. Together, these experimental and theoretical insights should aid in the exploitation of these giant polyoxometalates in solution for various applications.

Molecular dynamics algorithm enforcing energy conservation for microcanonical simulations.

A reversible algorithm [enforced energy conservation (EEC)] that enforces total energy conservation for microcanonical simulations is presented. The key point is the introduction of the discrete-gradient method to define the forces from the conservative potentials, instead of the direct use of the force field at the actual position of the particle. We have studied the performance and accuracy of the EEC in two cases, namely Lennard-Jones fluid and a simple electrolyte model. Truncated potentials that usually induce inaccuracies in energy conservation are used. In particular, the reaction field approach is used in the latter. The EEC is able to preserve energy conservation for a long time, and, in addition, it performs better than the Verlet algorithm for these kinds of simulations.

Model shape transitions of micelles: spheres to cylinders and disks.

We present a microscopic analysis of shape transitions of micelles of model linear nonionic surfactants. In particular, symmetric H(4)T(4) and asymmetric H(3)T(6) surfactants have been chosen for the study. In a previous work, it has been observed that symmetric surfactants have a strong tendency to prefer spherical micelles over a wide range of chemical potentials, while asymmetric surfactants undergo shape transitions between a spherical micelle at low concentration to other forms, mainly finite cylindrical micelles. This study combines the application of a two-dimensional single-chain mean-field theory (SCMFT) with Monte Carlo (MC) simulations of exactly the same systems. On the one hand, the characteristics of the SCMFT make this method suitable for free energy calculations, especially for small surfactants, due to the incorporation of relevant microscopic details in the model. On the other hand, MC simulations permit us to obtain a complete picture of the statistical mechanical problem, for the purpose of validation of the mean-field calculations. Our results reveal that the spherical shape for the symmetric surfactant is stable over a large range of surfactant concentrations. However, the asymmetric surfactant undergoes a complex shape transition that we have followed by calculating the standard chemical potential as a function of the aggregation number. The results indicate that the system forms prolate spheroids prior to developing short capped cylinders that gradually grow in length, with some oscillations in the energy of formation. The most important result of our work is the evidence of a bifurcation where, together with the elongated objects, the system can develop oblate aggregates and finally a torus shape similar to a red blood cell.

Dynamics of encapsulated water inside Mo132 cavities.

The structure and dynamics of water confined inside a polyoxomolybdate molecular cluster [{(Mo)Mo(5)O(21)(H(2)O)(6)}(12){Mo(2)O(4)(SO(4))}(30)](72-) metal oxide nanocapsule have been studied by means of molecular dynamics simulations under ambient conditions. Our results are compared to experimental data and theoretical analyses done in reverse micelles, for several properties. We observe that the characteristic three-dimensional hydrogen bond network present in bulk water is distorted inside the cavity where water organizes instead in concentric layered structures. Hydrogen bonding, tetrahedral order, and orientational distribution analyses indicate that these layers are formed by water molecules hydrogen bonded with three other molecules of the same structure. The remaining hydrogen bond donor/acceptor site bridges different layers as well as the whole structure with the hydrophilic inner side of the cavity. The most stable configuration of the layers is thus that of a buckyball with 12 pentagons and a variable number of hexagons. The geometrical constraints make it so that the bridges between the layers display a significant degree of frustration. The main modes of motion at short times are correlated fluctuations of the entire system with a characteristic frequency. Switches of water molecules between layers are rare events, due to the stability of the layers. At long times, the system shows a power law decay (pink noise) in properties like the fluctuations in the number of molecules in the structures and the total dipole moment. Such behavior has been attributed to the complex relaxation of the hydrogen bond network, and the exponents found are close to those encountered in bulk water for the relaxation of the potential energy. Our results reveal the importance of the competition between the confinement and the long-range structure induced in this system by the hydrogen bond network.

Aggregation of amphiphilic polymers in the presence of adhesive small colloidal particles.

The interaction of amphiphilic polymers with small colloids, capable to reversibly stick onto the chains, is studied. Adhesive small colloids in solution are able to dynamically bind two polymer segments. This association leads to topological changes in the polymer network configurations, such as looping and cross-linking, although the reversible adhesion permits the colloid to slide along the chain backbone. Previous analyses only consider static topologies in the chain network. We show that the sliding degree of freedom ensures the dominance of small loops, over other structures, giving rise to a new perspective in the analysis of the problem. The results are applied to the analysis of the equilibrium between colloidal particles and star polymers, as well as to block copolymer micelles. The results are relevant for the reversible adsorption of silica particles onto hydrophilic polymers, used in the process of formation of mesoporous materials of the type SBA or MCM, cross-linked cyclodextrin molecules threading on the polymers and forming the structures known as polyrotaxanes. Adhesion of colloids on the corona of the latter induce micellization and growth of larger micelles as the number of colloids increase, in agreement with experimental data.

Gated and differently functionalized (new) porous capsules direct encapsulates' structures: higher and lower density water.

By the deliberate choice of the internal ligands of the porous nanocapsules [{(Mo)Mo(5)}(12){Mo(2)(ligand)}(30)](n-), the respective cavities' shells can be differently sized/functionalized. This allows one to trap the same large number of water molecules, that is, 100 in a capsule cavity with formate ligands having a larger space available, as well as in a cavity containing sulfates and hypophosphites, that is, with less space. Whereas the 100 molecules fill the space completely in the second case in which they are organized in three shells, a four-shell system with underoccupation and broken hydrogen bonds is observed in the other case. This is an unprecedented result in terms of the structurally well defined special forms of "higher and lower density" water molecule assemblies. Precisely, by replacing the larger ligands in the mentioned nanocapsule type by formates, voids in the capsule cavity of (HC(NH(2))(2))(22)[{(HC(NH(2))(2))(20)+(H(2)O)(100)} subset{(Mo)Mo(5)O(21)(H(2)O)(6)}(12){Mo(2)O(4)(HCO(2))}(30)]ca. 200 H(2)O are generated that get filled with water molecules concomitant with an expansion of the three to four shell {H(2)O}(100) cluster. The water shells in both capsules containing different ligands are organized in the form of dodecahedra (partly with underoccupation) and a strongly distorted rhombicosidodecahedron spanned by a {H(2)O}(60)={(H(2)O)(5)}(12) aggregate. The well-defined water shells only emerge if cations cannot enter into the capsules, which is achieved by closing the pores with plugs/guests such as formamidinium cations. The work is based on the syntheses of two new compounds, related single-crystal X-ray diffraction studies, and molecular dynamics simulations, which show remarkably that water molecule shell structuring occurs in the capsules due to the confined conditions even in the case of open pores and at room temperature if cation uptake is prevented.

Dewetting of a stratified two-component liquid film on a solid substrate.

Dissipative particle dynamics simulations for the dewetting of a stratified thin film, composed of two immiscible fluids of different viscosity, reveal a nontrivial dewetting kinetics depending on the location of the more viscous liquid with respect to the solid substrate. Our simulations show that when the layer of higher viscosity is in contact with the substrate, the energy loss through viscous dissipation is concentrated near the contact line. The kinetics of dewetting is then inversely proportional to the viscosity of the lower layer eta_{A} (with respect to eta_{B} for the upper layer). However, when the liquid of higher viscosity is on the upper layer, a new dynamics appears with a dewetting velocity V approximately eta_{B};{-0.44}eta_{A};{-0.56}. The difference between these two scenarios lies in the different routes through which the interfacial energy is converted into heat by viscous dissipation in the rim.

Histogram reweighting method for dynamic properties.

The histogram reweighting technique, widely used to analyze Monte Carlo data, is shown to be applicable to dynamic properties obtained from molecular dynamics simulations. The theory presented here is based on the fact that the correlation functions in systems in thermodynamic equilibrium are averages over initial conditions of trajectory functions, the latter depending on the volume of the system, the total number of particles, and the classical Hamiltonian. Thus, the well-known histogram reweighting method can be almost straightforwardly applied to reconstruct the probability distribution of initial states at different thermodynamic conditions, without extra computational effort. Correlation functions and transport coefficients are obtained with this method from few simulation data sets.

Computing the Soret coefficient in aqueous mixtures using boundary driven nonequilibrium molecular dynamics.

We have computed the Soret coefficient in aqueous mixtures using a boundary driven nonequilibrium molecular dynamics algorithm and standard molecular force fields. The choice of this specific approach is justified by the nature of the mixtures studied here. Four aqueous solutions, including methanol, ethanol, acetone, and dimethyl-sulfoxide (DMSO) have been studied at ambient conditions for different compositions. The experimental behavior of water-alcohol mixtures was reproduced, including the change of sign of the Soret coefficient with composition, in excellent agreement with existing experimental data. The methodology has been applied to obtain pure predictions for water-acetone and water-DMSO where no experimental data are accessible. A change of sign is also observed in the same range of composition as in water-alcohol mixtures. It is suggested that the nature and strength of the molecular interactions, rather than the mass or shape ratio of the components, dominates the behavior of the Soret coefficient versus composition for the aqueous associating mixtures studied here.

Prediction of the critical micelle concentration in a lattice model for amphiphiles using a single-chain mean-field theory.

A single-chain mean-field theory is used to predict the properties of binary surfactant solutions including the critical micelle concentration (cmc). In particular, the cmc of two symmetric nonionic amphiphiles is calculated as a function of temperature in order to analyze the validity of the ideal mixing assumption, often employed in the mass action model. On comparing against literature Monte Carlo results for the same lattice model we find that although it is applicable at low temperatures and hence cmcs at low amphiphile concentrations, at higher temperatures it becomes necessary to correct for the nonideal mixing of the free chain-free chain bulk interaction. We find that a simplistic model taking into account only the repulsive interaction is sufficient to restore the excellent quantitative agreement found between a single-chain mean-field theory calculations and literature molecular simulation results at the low temperature limit.

Polyoxometalates in solution: Molecular dynamics simulations on the alpha-PW12O40(3-) Keggin anion in aqueous media.

The Keggin anion, PW12O40(3-), is one of the most representative polyoxometalates (POMs). In recent years increasing theoretical work focused on this family of compounds has explained or even predicted some of their properties using quantum mechanics methods. In this report we applied for the first time molecular dynamics (MD) to the title compound to analyze its interactions with water. We used three atomic charge definitions (Mulliken, ChelpG, and formal charges) to carry out MD simulations. The results show that the terminal oxygens of the cluster are invariably most effectively solvated by water because of their prominent position within the framework. On the other hand, bridging oxygens, which are confined in more internal positions, concentrate a smaller portion of the whole solvation. We investigated the hydrogen bonds existing between water and the cluster, confirming that the terminal positions form more contacts with H2O than any other site of the cluster. In the end, the lifetime of such contacts is longer with bridging oxygens, presumably due to their higher atomic charge.