Research Group of Molecular and Mesoscopic Modelling

Micro- and mesoscale (below 1 μm) phenomena play a dominant role in the resulting macroscopic properties for a wide variety of material types, e.g., soft matter, biological matter, complex fluids, composite materials or additively-manufactured materials.  Predictive modelling and simulation of materials response requires multiscale modelling at length and time scales ranging from atomistic to continuum.  Modelling and simulation at this scale is far from amenable to atomistic scale approaches, while continuum scale simulations lack the fidelity to properly include material microstructure and rely on phenomenological models to effectively recover and/or estimate the physical processes occurring at the micro and mesoscales.  Thus, a key gap in the theoretical foundation and computational methods exists at the microscale/mesoscale across various scientific fields, where studies have been performed over an extensive scope of materials and applications, including, but not limited to, the life sciences (proteins, colloidal suspensions, bio-membranes, and micelles), industrial applications (surfactants, asphaltenes, and viscoelastic fluids), national defence applications (energetic material composites and liquid propellants), and novel materials (self-assembled block copolymers/nanoparticles). We are developing novel particle-based mesoscale modelling tools for predicting the spatial-temporal energy release, heat transfer, mass transfer and chemical reactivity of energetic materials.  We aim to establish a rigorous, well-founded, and general computational approach for systematic study directed toward identifying and characterising fundamental aspects of the dynamic response of energetic materials.

The physics of fluids in nanoconfinement such as micropores or at interfaces with solid surfaces such as graphene differs significantly from the physics in the bulk. The significant and often counterintuitive changes in fluid behaviour in nanospace are caused by dominant effects of interactions between the fluid molecules and surface atoms. Results of strong fluid-solid interactions are inhomogeneities in fluid density, appearance of specific fluid phases, slow-down of fluid diffusion and transport, e.g. aqueous solutions in hydrophilic mineral micropores, or speed-up of fluid diffusion and transport, e.g. water in carbon nanotubes. We use molecular dynamics and Monte Carlo simulations to understand, describe and quantify at the molecular level effects of nanospace and fluid-solid interactions on the phase, structural, diffusion and transport behaviour of nanoconfined fluids. A molecular-level understanding of nanoconfined fluid behaviour plays a critical role in rational design of porous materials or in development of graphene-based technologies.

Our understanding of natural hazards such as earthquakes, landslides or soil liquefaction is tightly connected to research of granular media. From a mechanical perspective, granular materials stand between solids and fluids, as they are able to sustain considerable shear stresses under some conditions, while they may flow as fluids under other conditions. Questions of friction and the mechanical stability of granular systems are therefore the key to understanding natural hazards. Ongoing observations have shown that sites of hydrothermal activity are most prone to occurrence of the hazards. The presence of fluids, specifically the pore-fluid pressure, appears to be an important agent in triggering instable deformation. However, the mechanisms by which pore pressure rises and triggers an instability is not fully understood. We therefore develop physical models of coupled deformation of rock grains and pore fluid to understand their macroscopic behaviour.

There is a lot of unambiguous experimental evidence that modification of solid surfaces may dramatically alter the structure and phase behaviour of the ambient fluid. It is already well known that geometrical or chemical modification of solid surfaces can induce completely new types of phase transitions and critical phenomena, related to but distinct from wetting occurring at planar (and ideally smooth) surfaces. Using methods of statistical mechanics (Density Functional Theory, Effective Hamiltonian Models, Scaling,…) we endeavour to capture the essential physics of these phenomena and understand the subtle interplay between the surface geometry, molecular interactions, interfacial fluctuations and measureable characteristics of the system (such as adsorption). Our main focus is on surfaces whose heterogeneity is characterized by a nanoscopically small length-scale, a research area which is still in its infancy and revealing a number of counter-intuitive phenomena that cannot be captured by standard tools. Beside its theoretical importance, the research has also tight connections with modern nanotechnologies dealing with fabrication of functional materials or for a development of devices controlling flows of extremely tiny volumes of liquid, such as “lab-on-a-chip”.

Graphene-based applications, such as supercapacitors or capacitive deionization, take place in an aqueous environment, and they benefit from molecular-level insights into the behavior of aqueous electrolyte solutions in single-digit graphene nanopores with a size comparable to a few molecular diameters. Under single-digit graphene nanoconfinement (smallest dimension < 2 nm), water and ions behave drastically different than in the bulk. Most aqueous electrolytes in graphene-based applications as well as in nature contain a mix of electrolytes. We study several prototypical aqueous mixed alkali-chloride electrolytes containing an equimolar fraction of Li/Na, Li/K, or Na/K cations confined between neutral and positively or negatively charged parallel graphene sheets. The strong hydration shell of small Li⁺ vs a larger Na⁺ or large K⁺ with weaker or weak hydration shells affects the interplay between the ions’ propensity to hydrate or dehydrate under the graphene nanoconfinement and the strength of the ion–graphene interactions mediated by confinement-induced layered water. We perform molecular dynamics simulations of the confined mixed-cation electrolytes using the effectively polarizable force field for electrolyte–graphene systems and focus on a relation between the electrochemical adsorption and structural properties of the water molecules and ions and their diffusion behavior. The simulations show that the one-layer nanoslits have the biggest impact on the ions’ adsorption and the water and ions’ diffusion. The positively charged one-layer nanoslits only allow for Cl⁻ adsorption and strengthen the intermolecular bonding, which along with the ultrathin confinement substantially reduces the water and Cl⁻ diffusion. In contrast, the negatively charged one-layer nanoslits only allow for the adsorption of weakly hydrated Na⁺ or K⁺ and substantially break up the non-covalent bond network, which leads to the enhancement of the water and Na⁺ or K⁺ diffusion up to or even above the bulk diffusion. In wider nanoslits, cations adsorb closer to the graphene surfaces than Cl⁻’s with preferential adsorption of a weakly hydrated cation over a strongly hydrated cation. The positive graphene charge has an intuitive effect on the adsorption of weakly hydrated Na⁺’s or K⁺’s and Cl⁻’s and a counterintuitive effect on the adsorption of strongly hydrated Li⁺’s. On the other hand, the negative surface charge has an intuitive effect on the adsorption of both types of cations and only mild intuitive or counterintuitive effects on the Cl⁻ adsorption. The diffusion of water molecules and ions confined in the wider nanoslits is reduced with respect to the bulk diffusion, more for the positive graphene charge, which strengthened the intermolecular bonding, and less for the negative surface charge, which weakened the non-covalent bond network.

• E. Rezlerová, F. Moučka, M. Předota, M. Lísal: Structure and self-diffusivity of mixed-cation electrolytes between neutral and charged graphene sheets, J. Chem. Phys. 160, 094701, 2024. DOI

This article addresses the debate about the correct application of Green–Kubo expressions for transport coefficients from dissipative particle dynamics simulations. We demonstrate that the Green–Kubo expressions are valid provided that (i) the dynamic model conserves the physical property, whose transport is studied, and (ii) the fluctuations satisfy detailed balance. As a result, the traditional expressions used in molecular dynamics can also be applied to dissipative particle dynamics simulations. However, taking the calculation of the shear viscosity as a paradigmatic example, a random contribution, whose strength scales as 1/dt^1/2, with dt the timestep, can cause difficulties if the stress tensor is not separated into the different contributions. We compare our expression to that of Ernst and Brito (M. H. Ernst and R. Brito, Europhys. Lett. 73, 183-189, 2006), which arises from a diametrically different perspective. We demonstrate that the two expressions are completely equivalent and find exactly the same result both analytically and numerically. We show that the differences are not due to the lack of time-reversibility but instead from a preaveraging of the random contributions. Despite the overall validity of Green–Kubo expressions, we find that the Einstein–Helfand relations (D. C. Malaspina et al. Phys. Chem. Chem. Phys. 25, 12025-12040, 2023) do not suffer from the need to decompose the stress tensor and can readily be used with a high degree of accuracy. Consequently, Einstein–Helfand relations should be seen as the preferred method to calculate transport coefficients from dissipative particle dynamics simulations.

• D. C. Malaspina, M. Lísal, J. P. Larentzos, J. K. Brennan, A. D. Mackie, J. Bonet Avalos, Green-Kubo expressions for transport coefficients from dissipative particle dynamics simulations revisited. Phys. Chem. Chem. Phys. 26, 1328-1339, 2024. DOI

The structure of aqueous solutions of alkali metal halides is studied under ambient thermodynamic conditions and concentrations from infinite dilution to supersaturation using molecular dynamics simulations of phase-transferable polarizable models. The results of the solution densities, radial distribution functions, 3D spatial distribution functions, the properties of hydrogen and other noncovalent bonds in solutions, hydration numbers, coordination numbers, numbers of contact cation-anion pairs, and other statistics of the number of ions hydrated simultaneously by a shared water molecule are systematically presented. In particular, the results show and quantify how the strengths of the hydration bonds of different ions vary and how the hydration numbers decrease with increasing concentration in parallel with an increase in the number of contact cation-anion pairs. In most cases, they completely compensate for the loss of water-ion bonds by an increase in cation-anion bonds. The exceptions are solutions based on the Li+ cation, which retain a solid hydration shell even at high concentrations. This behavior is conceptualized based on three imaginary driving forces: the first dominating at low concentrations and causing full hydration of the ions, the second representing a lack of water necessary for full hydration of the ions and increasing with increasing concentration, and the third attracting counterions to the water-unoccupied sites of the hydration shells and also increasing with concentration. This concept can be used not only to understand the structural behavior of homogeneous electrolytes in thermodynamic equilibrium but also to study phenomena that involve preferential adsorption of ions on electrodes, nanochannels, or porous materials. The data obtained for the number and strength of hydration bonds and ion pairs can also be used in further studies to elucidate the diffusion behavior, viscosity, and conductivity of aqueous electrolyte solutions.

• J. Dočkal, P. Mimrová, M. Lísal, F. Moučka, Structure of aqueous alkali metal halide electrolyte solutions from molecular simulations of phase-transferable polarizable models. J. Mol. Liq. 394, 123797, 2024. DOI
  

Understanding the microscopic behavior of aqueous electrolyte solutions in graphene-based ultrathin nanochannels is important in nanofluidic applications such as water purification, fuel cells, and molecular sensing. Under extreme confinement (< 2 nm), the properties of water and ions differ drastically from those in the bulk phase. We studied the structural and diffusion behavior of prototypical aqueous solutions of electrolytes (LiCl, NaCl, and KCl) confined in both neutral and positively-, and negatively-charged graphene nanochannels. We performed molecular dynamics simulations of the solutions in the nanochannels with either one, two- or three-layer water structures using the effectively polarisable force field for graphene.

We analyzed the structure and intermolecular bond network of the confined solutions along with their relation to the self-diffusivity of water and ions. The simulations show that Na and K cations can more easily rearrange their solvation shells under the graphene nanoconfinement and adsorb on the graphene surfaces or dissolve in the confinement-induced layered water than the Li cation. The negative surface charge together with the presence of ions orient water molecules with hydrogens towards the graphene surfaces, which in turn weakens the intermolecular bond network. The one-layer nanochannels have the biggest effect on the water structure and intermolecular bonding as well as on the adsorption of ions with only co-ions entering these nanochannels. The self-diffusivity of confined water is strongly reduced with respect to the bulk water and decreases with diminishing nanochannel heights except for the negatively-charged one-layer nanochannel. The self-diffusivity of ions also decreases with the reducing the nanochannel heights except for the self-diffusivity of cations in the negatively-charged one-layer nanochannel, evidencing cooperative diffusion of confined water and ions. Due to the significant break-up of the intermolecular bond network in the negatively-charged one-layer nanochannel, self-diffusion coefficients of water and cations exceed those for the two- and three-layer nanochannels and become comparable to the bulk values.

• E. Rezlerová, F. Moučka, M. Předota, M. Lísal, Structure and self-diffusivity of alkali-halide electrolytes in neutral and charged graphene nanochannels, Phys. Chem. Chem. Phys. 25, 21579-21594, 2023. DOI

Dissipative particle dynamics (DPD) is a promising candidate technique for modeling the rheological properties of soft matter systems. However, several methodological issues inhibit its exploitation as a computational rheology tool. In this work, we focus on the development of an automated method to compute transport properties from equilibrium simulation with particular attention to the assessment of the Green-Kubo approach reliability and computational feasibility for a large set of simple DPD systems with increasing Schmidt number. Furthermore, we investigate the time step size dependency of dynamic properties and the role of different time integration schemes. In particular, we assess the performance of the Shardlow-splitting algorithm against the most popular modified velocity-Verlet algorithm. We consider, for the first time, the application of the Shardlow-splitting algorithm to the transverse DPD thermostat in different friction regimes relying on systematic numerical experiments. In addition, we make use of these findings to perform a multi-parametric study aiming to investigate the Schmidt number relationship with the effective friction coefficient.

• N. Lauriello, G. Boccardo, D. Marchisio, M. Lísal, A. Buffo, Development of an automated reliable method to compute transport properties from DPD equilibrium simulations: Application to simple fluids, Comput. Phys. Commun. 291, 108843, 2023. DOI

Fatty alcohols (CnFAs) combined with cationic surfactants are common ingredients of lamellar-phase personal care liquids. We employ mesoscopic modelling to study how surfactant-CnFA monolayers originating from the corresponding bilayers of the personal care liquids interact with the human hair surface above and below the fluid-gel transition temperature as well as in- and out-of-equilibrium. For the mono-layer model, we consider the single-tail cationic surfactant cetyltrimethylammonium chloride (CTAC) and an excess of CnFAs with their alkyl tail length equal to or longer than the CTAC alkyl tail length.

The hair surface mimics keratin surface proteins covered by a film of lipid chains covalently bonded to the proteins. Our modelling shows the formation of a dense adsorbed layer due to the interactions of the CTAC and CnFA alkyl tails with the hydrophobic hair surface. The adsorption and the behaviour of the adsorbed layer is different under fluid and gel conditions. The differences are related to the structure of the adsorbed layer as characterised by density profiles across the adsorbed layer and the orientational order parameters of the chains within the adsorbed layer. Under steady-state shearing (an approximation of real, non-equilibrium conditions), increasing the shear rate above a threshold leads to continuous or abrupt desorption of the CTAC and CnFA chains under fluid or gel conditions, respectively; the desorbed chains can then form various self-assembled structures in the bulk solution. The underlying mechanism of CTAC and CnFA desorption from the adsorbed layers is closely related to the corresponding adsorption mechanism.

• K. Šindelka, A. Kowalski, M. Cooke, C. Mendoza, M. Lísal, Interactions of cationic surfactant-fatty alcohol monolayers with natural human hair surface: Insights from dissipative particle dynamics, J. Mol. Liq. 375, 121385, 2023. DOI

Cores of geological faults are granular layers often filled with fluid, which tend to dilate on increase in slip rate, for example, during an earthquake nucleation. Dilation of the fault causes a drop in pore fluid pressure which in turn increases strength of the fault. This effect, known as dilatancy hardening, may suppress generation of earthquakes. However, the magnitude of fluid depressurization during dynamic slip is difficult to measure. We systematically study pore pressure drop following the onset of slip in fluid-saturated granular layers, using theory and numerical simulations. We find two regimes with distinct evolutions of pore pressure and dilatation. For slow slip rate, high permeability or high stress acting upon the layer (drained conditions), the pore pressure drop is quickly restored and its ability to suppress earthquake nucleation is weak. Faster slip rate, lower permeability or stress (undrained conditions) lead to an enduring pore pressure drop, which significantly increases strength of the layer. The growing magnitude of the pore pressure drop with slip rate promotes the stability of shear rupture. Our results can be used to better constrain drainage conditions associated with changes in slip rate, the magnitude of the generated pore pressure and the corresponding fault strengthening.

• Parez, S.,  Kozakovic, M., and  Havlica, J., Pore pressure drop during dynamic rupture and conditions for dilatancy hardening. Journal of Geophysical Research: Solid Earth 128, e2023JB026396, 2023.  https://doi.org/10.1029/2023JB026396

This article deals with numerical simulation via Discrete Element Method (DEM) in a flighted rotary drum with a varying number of flights from 1 to 25 at a constant load of particles. The granular kinematics is studied in the three regimes: under-, design- and over-loading. It was shown that the behavior of granular material during the unloading process from the flights governs the regime transitions. The cross-sectional distribution analysis of the granular material in the different parts of the drum established the existence of two dense media (bed and flights) and one dilute medium (in the airborne portion). Since one of the criteria that affect the efficiency of heat transport between the gas and the particle media is the number of particles in the airborne phase, three recognition methods are proposed to detect the airborne portion. Based on geometry, velocity magnitude, and minimal separation distance, the comparison between them provides a basis for discussion about the relevance of each method in detecting airborne particles.

• M. Kozakovic, J. Havlica, L. Le Guen, S. Parez, F. Huchet, Recognition of the granular airborne portion in a flighted rotary drum, Powder Technology 425, 118565, 2023. DOI

We consider fluid adsorption on a planar wall patterned by a single stripe of width L of a wet material embedded within a partially wet surface. If L were infinite we suppose that the wall would exhibit a first-order wetting transition at temperature T_w, and an associated line of pre-wetting transitions extending off bulk-coexistence. For finite widths L, we a) determine the finite-size scaling shift of the wetting temperature T_w(L), b) derive an equation, analogous to the Kelvin equation for capillary condensation, for the shift of the pre-wetting line but involving boundary tensions, and a boundary contact angle and c) highlight the scaling and conformally invariant properties of condensed liquid drops lying above the stripe at bulk coexistence. We point out analogies between these predictions and those for capillary condensation and criticality in a parallel plate geometry and test them against numerical results obtained from a microscopic non-local classical density functional theory.

• M. Pospíšil, A. O. Parry, A. Malijevský, Phase transitions and droplet shapes above a wetting stripe, J. Mol. Liq. 381, 121834, 2023. DOI

It has been known that classical molecular simulations of concentrated aqueous electrolyte solutions are limited by the ability of microscopic models (force fields) to predict the correct values of salt solubility, which also requires their ability to faithfully model crystalline salts. Previous simulation studies have often focused on the solubility of anhydrous crystalline salts, but virtually never on crystalline hydrates, except for hydrohalite, NaCl^.2H₂O, despite there are at least 23 experimentally known different hydrates that can precipitate from alkali-halide solutions. This work attempts to fill this gap in hydrate simulation studies by systematically investigating the ability of the best force fields selected to qualitatively capture the stability of the individual phases of various alkali-halide hydrates and to quantitatively predict their lattice parameters. First, we show that the nonpolarizable force fields studied often fail to model hydrates containing the Li⁺ cations, whereas the polarizable force fields recently refined by us are able to model all the hydrates except for LiCl^.H₂O. Second, we further refine our FFs for Li⁺ to yield stable LiCl^.H₂O. Third, our simulations clarify the positions of the Li⁺ cations in the phases of LiBr^.H₂O and LiI^.H₂O, whose distributions were previously described only as stochastic. As a byproduct, a simple and reliable simulation methodology suitable also for complex polarizable models and nonorthorhombic crystal lattices is proposed and tested, based on simulations of finite crystals floating in a vacuum.

• P. Matysová, M. Lísal, F. Moučka, Molecular simulations of alkali metal halide hydrates, J. Mol. Liq. 384, 122197, 2023. DOI

We demonstrate that contrary to general belief, the standard Einstein–Helfand (EH) formulas are valid for the evaluation of transport coefficients of systems containing dissipative and random forces provided that for these mesoscopic systems: (i) the corresponding conservation laws are satisfied, and (ii) the transition probabilities satisfy detailed balance. Dissipative particle dynamics (DPD) and energy-conserving DPD methods (DPDE), for instance, are archetypical of such mesoscopic approaches satisfying these properties. To verify this statement, we have derived a mesoscopic heat flux form for the DPDE method, suitable for the calculation of the thermal conductivity from an EH expression. We have compared EH measurements against non-equilibrium simulation values for different scenarios, including many-body potentials, and have found excellent agreement in all cases. The expressions are valid notably for systems with density- and temperature-dependent potentials, such as the recently developed generalized DPDE method (GenDPDE) [Avalos et al., Phys. Chem. Chem. Phys. 21, 24891, 2019]. We thus demonstrate that traditional EH formulas in equilibrium simulations can be widely used to obtain transport coefficients, provided that the appropriate expression for the associated flux is used.

• C. Malaspina, M. Lísal, J. P. Larentzos, J. K. Brennan, A. D. Mackie, J. Bonet Avalos, Transport coefficients from Einstein-Helfand relations using standard and energy-conserving dissipative particle dynamics methods, Phys. Chem. Chem. Phys. 25, 12025-12040, 2023. DOI

Organic-shale formations are unconventional gas reservoirs with broad pore size distributions. Shale consists of two distinct components: organic matter and clay minerals. The size of pores in the organic matter is mostly concentrated at less than six nanometres and these micropores, and small mesopores provide the majority of adsorption surface area and gas storage volume. In these nanometre-sized pores, the geofluid behaviour becomes significantly different from the bulk behaviour due to the strong solid-fluid interactions and other confinement effects. Understanding the properties of fluids such as methane, ethane, propane and carbon dioxide in narrow shale pores is critical for identifying ways to deploy shale gas technology with reduced environmental impact. Specifically, methane is a proxy of the shale gas and ethane, and propane are minor shale-gas components. Further, carbon dioxide is used for the enhanced shale-gas recovery. We employ molecular-level simulations to explore adsorption, diffusion, and transport of methane, ethane, propane and carbon dioxide in realistic models of organic-shale materials at a representative shale-reservoir temperature and pressures. We first use Hybrid Reverse Monte Carlo with experimental pair distribution functions to build dual-porosity kerogen models corresponding to an immature marine kerogen from the Eagle Ford Play and a mature marine kerogen from the clay-rich Marcellus Play. We then employ Grand Canonical Monte Carlo simulations to study the fluid adsorption in the porous kerogen structures. We complete the adsorption studies by simulating the self-diffusivity, collective diffusivity, and transport diffusivity of the adsorbed fluid molecules in the shale kerogens using equilibrium and non-equilibrium molecular dynamics.

• E. Rezlerová, S. K. Jain, M. Lísal, Adsorption, diffusion, and transport of C1 to C3 alkanes and carbon dioxide in dual-porosity kerogens: Insights from molecular simulations, Energy Fuels 37, 492-508, 2023. DOI