Computing induced charges in inhomogeneous dielectric media: application in a Monte Carlo simulation of complex ionic systems

GCMe: Efficient Implementation of the Gaussian Core Model with Smeared Electrostatic Interactions for Molecular Dynamics Simulations of Soft Matter Systems

In recent years, molecular dynamics (MD) simulations have emerged as an essential tool for understanding the structure, dynamics, and phase behavior of charged soft matter systems. To explore phenomena across greater length and time scales in MD simulations, molecules are often coarse-grained for better computational performance. However, commonly used force fields represent particles as hard-core interaction centers with point charges, which often overemphasizes the packing effect and short-range electrostatics, especially in systems with bulky deformable organic molecules and systems with strong coarse-graining. This underscores the need for an efficient soft-core model to physically capture the effective interactions between coarse-grained particles. To this end, we implement a soft-core model uniting the Gaussian core model with smeared electrostatic interactions that is phenomenologically equivalent to recent theoretical models. We first parametrize it generically using water as the model solvent. Then, we benchmark its performance in the OpenMM toolkit for different boundary conditions to highlight a computational speedup of up to 34 × compared to commonly used force fields and existing implementations. Finally, we demonstrate its utility by investigating how boundary polarizability affects the adsorption behavior of a polyelectrolyte solution on perfectly conducting and nonmetal boundaries.

Electrostatic interactions between rough dielectric particles

From colloid suspension to particle aggregation in protoplanetary formation, electrostatic attraction and repulsion between particles is a key mechanism behind the aggregation and clustering of particles. Although most studies have focused on canonical spherical particles, it remains unclear how nonspherical and rough dielectric particles interact and whether the complicated interplay between roughness and charge distribution affects their force couplings. Here a boundary-element method model was leveraged to study electrostatic interactions between charged dielectric particles with modular, axisymmetric surface features. Charge accumulation at convex surface asperities decreases the strength of electrostatic interactions between particles, and the sensitivity of the electrostatic force to the particle surface roughness and orientation is especially apparent at small particle separations. Surface interactions between the particle near-contact regions were isolated to determine the degree that near-contact interactions dictate the relationship between the net electrostatic force and the particle roughness and orientation. A correction factor ΔF is introduced to recover higher order dielectric effects from a low order analytical model. Finally, implications of surface charge asymmetries produced for different particle orientations and surface roughnesses on the long-standing problem of triboelectrification are discussed.

Ideal conductor/dielectric model (ICDM): A generalized technique to correct for finite-size effects in molecular simulations of hindered ion transport

Molecular simulations serve as indispensable tools for investigating the kinetics and elucidating the mechanism of hindered ion transport across nanoporous membranes. In particular, recent advancements in advanced sampling techniques have made it possible to access translocation timescales spanning several orders of magnitude. In our prior study [Shoemaker et al., J. Chem. Theory Comput. 18, 7142 (2022)], we identified significant finite size artifacts in simulations of pressure-driven hindered ion transport through nanoporous graphitic membranes. We introduced the ideal conductor model, which effectively corrects for such artifacts by assuming the feed to be an ideal conductor. In the present work, we introduce the ideal conductor dielectric model (Icdm), a generalization of our earlier model, which accounts for the dielectric properties of both the membrane and the filtrate. Using the Icdm model substantially enhances the agreement among corrected free energy profiles obtained from systems of varying sizes, with notable improvements observed in regions proximate to the pore exit. Moreover, the model has the capability to consider secondary ion passage events, including the transport of a co-ion subsequent to the traversal of a counter-ion, a feature that is absent in our original model. We also investigate the sensitivity of the new model to various implementation details. The Icdm model offers a universally applicable framework for addressing finite size artifacts in molecular simulations of ion transport. It stands as a significant advancement in our quest to use molecular simulations to comprehensively understand and manipulate ion transport processes through nanoporous membranes.

Charge equilibration model with shielded long-range Coulomb for reactive molecular dynamics simulations

Atomic description of electrochemical systems requires reactive interaction potential to explicitly describe the chemistry between atoms and molecules and the evolving charge distribution and polarization effects. Calculating Coulomb electrostatic interactions and polarization effects requires a better estimate of the partial charge distribution in molecular systems. However, models such as reactive force fields and charge equilibration (QEq) include Coulomb interactions up to a short-distance cutoff for better computational speeds. Ignoring long-distance electrostatic interaction affects the ability to describe electrochemistry in large systems. We studied the long-range Coulomb effects among charged particles and extended the QEq method to include long-range effects. By this extension, we anticipate a proper account of Coulomb interactions in reactive molecular dynamics simulations. We validate the approach by computing charges on a series of metal-organic frameworks and some simple systems. Results are compared to regular QEq and quantum mechanics calculations. The study shows slightly overestimated charge values in the regular QEq approach. Moreover, our method was combined with Ewald summation to compute forces and evaluate the long-range effects of simple capacitor configurations. There were noticeable differences between the calculated charges with/without long-range Coulomb interactions. The difference, which may have originated from the long-range influence on the capacitor ions, makes the Ewald method a better descriptor of Coulomb electrostatics for charged electrodes. The approach explored in this study enabled the atomic description of electrochemical systems with realistic electrolyte thickness while accounting for the electrostatic effects of charged electrodes throughout the dielectric layer in devices like batteries and emerging solid-state memory.

Probing Accuracy-Speedup Tradeoff in Machine Learning Surrogates for Molecular Dynamics Simulations

The performance promise of machine learning surrogates of molecular dynamics simulations of soft materials is significant but generally comes at the cost of acquiring large training datasets to learn the complex relationships between input soft material attributes and output properties. Under the constraint of limited high-performance computing resources, optimizing the size of the training datasets becomes paramount. Using an artificial neural network based surrogate for molecular dynamics simulations of confined electrolytes, we explore the tradeoff between surrogate accuracy and computational gains. Accuracy is assessed by computing the root-mean-square errors between the surrogate predictions and the ground truth results obtained via molecular dynamics simulations. The computational performance is judged by evaluating the speedup which incorporates the training dataset creation time. Improvement in accuracy occurs with a loss of speedup, which scales as the inverse of the training dataset size. The link between surrogate generalizability and the accuracy-speedup tradeoff is assessed by examining the errors incurred in surrogate predictions on unseen, interpolated input variables and developing a net speedup metric to capture the associated gains.

A molecular perspective on induced charges on a metallic surface

Understanding the response of the surface of metallic solids to external electric field sources is crucial to characterize electrode-electrolyte interfaces. Continuum electrostatics offer a simple description of the induced charge density at the electrode surface. However, such a simple description does not take into account features related to the atomic structure of the solid and to the molecular nature of the solvent and of the dissolved ions. In order to illustrate such effects and assess the ability of continuum electrostatics to describe the induced charge distribution, we investigate the behavior of a gold electrode interacting with sodium or chloride ions fixed at various positions, in a vacuum or in water, using all-atom constant-potential classical molecular dynamics simulations. Our analysis highlights important similarities between the two approaches, especially under vacuum conditions and when the ion is sufficiently far from the surface, as well as some limitations of the continuum description, namely, neglecting the charges induced by the adsorbed solvent molecules and the screening effect of the solvent when the ion is close to the surface. While the detailed features of the charge distribution are system-specific, we expect some of our generic conclusions on the induced charge density to hold for other ions, solvents, and electrode surfaces. Beyond this particular case, the present study also illustrates the relevance of such molecular simulations to serve as a reference for the design of improved implicit solvent models of electrode-electrolyte interfaces.

Molecular Simulation of Electrode-Solution Interfaces

Many key industrial processes, from electricity production, conversion, and storage to electrocatalysis or electrochemistry in general, rely on physical mechanisms occurring at the interface between a metallic electrode and an electrolyte solution, summarized by the concept of an electric double layer, with the accumulation/depletion of electrons on the metal side and of ions on the liquid side. While electrostatic interactions play an essential role in the structure, thermodynamics, dynamics, and reactivity of electrode-electrolyte interfaces, these properties also crucially depend on the nature of the ions and solvent, as well as that of the metal itself. Such interfaces pose many challenges for modeling because they are a place where quantum chemistry meets statistical physics. In the present review, we explore the recent advances in the description and understanding of electrode-electrolyte interfaces with classical molecular simulations, with a focus on planar interfaces and solvent-based liquids, from pure solvent to water-in-salt electrolytes.

Dynamics of ionic flows via Poisson-Nernst-Planck systems with local hard-sphere potentials: Competition between cations

We study a quasi-one-dimensional steady-state Poisson-Nernst-Planck type model for ionic flows through a membrane channel with three ion species, two positively charged with the same valence and one negatively charged. Bikerman's local hard-sphere potential is included in the model to account for ion sizes. The problem is treated as a boundary value problem of a singularly perturbed differential system. Under the framework of a geometric singular perturbation theory, together with specific structures of this concrete model, the existence and uniqueness of solutions to the boundary value problem for small ion sizes is established. Furthermore, treating the ion sizes as small parameters, we derive an approximation of individual fluxes, from which one can further study the qualitative properties of ionic flows and extract concrete information directly related to biological measurements. Of particular interest is the competition between two cations due to the nonlinear interplay between finite ion sizes, diffusion coefficients and boundary conditions, which is closely related to selectivity phenomena of open ion channels with given protein structures. Furthermore, we are able to characterize the distinct effects of the nonlinear interplays between these physical parameters. Numerical simulations are performed to identify some critical potentials which play critical roles in examining properties of ionic flows in our analysis.

Tunable interactions between particles in conically rotating electric fields

Tunable interactions between colloidal particles in external conically rotating electric fields are calculated, while the (vertical) axis of the field rotation is normal to the (horizontal) particle motion plane. The comparison of different approaches, including the methods of noninteracting, self-consistent dipoles, and the boundary element method, indicates that the last method is the most suitable for tunable interaction analysis. Thorough analysis, performed for interactions in pairs and clusters of colloidal particles, indicate that two- and three-body interactions make the main contributions in the interaction energy, while the effect of high-order terms is negligible. The tunable interactions are determined by the dielectric properties of the particles and solvent and can be changed in a wide range, providing a rich variety for the experimental "design" of different interactions, including repulsion, attraction, combination of short-range repulsion with long-range attraction, barrier-type interactions with short-range attraction and long-range repulsion, and double-scale repulsive (core-shoulder) interactions. These conclusions can be generalized for magnetically induced tunable interactions. The results indicate that tunable interactions can be widely applied in self-assembly and particle-resolved studies of generic phenomena in fluids and crystals, and, therefore, are of broad interest in the fields of chemical physics, physical chemistry, materials science, and soft matter.

Boundary Layer Effects on Ionic Flows Via Classical Poisson-Nernst-Planck Systems

A quasi-one-dimensional steady-state Poisson-Nernst-Planck model of two oppositely charged ion species through a membrane channel is analyzed. The model problem is treated as a boundary value problem of a singularly perturbed differential system. Our analysis is based on the geometric singular perturbation theory but, most importantly, on specific structures of this concrete model. The existence and (local ) uniqueness of solutions to the boundary value problem is established. In particular, an approximation of both the individual flux and the I-V (current-voltage) relation are derived explicitly from the zeroth order approximation (in ") solutions, from which the boundary layer effects on ionic flows are studied in great details.

Simulations of Coulomb systems confined by polarizable surfaces using periodic Green functions

We present an efficient approach for simulating Coulomb systems confined by planar polarizable surfaces. The method is based on the solution of the Poisson equation using periodic Green functions. It is shown that the electrostatic energy arising from the surface polarization can be decoupled from the energy due to the direct Coulomb interaction between the ions. This allows us to combine an efficient Ewald summation method, or any other fast method for summing over the replicas, with the polarization contribution calculated using Green function techniques. We apply the method to calculate density profiles of ions confined between the charged dielectric and metal surfaces.

Controlling ion transport through nanopores: modeling transistor behavior

We present a modeling study of a nanopore-based transistor computed by a mean-field continuum theory (Poisson–Nernst–Planck, PNP) and a hybrid method including particle simulation (Local Equilibrium Monte Carlo, LEMC) that is able to take ionic correlations into account including the finite size of ions.

Multiscale modeling of a rectifying bipolar nanopore: Comparing Poisson-Nernst-Planck to Monte Carlo

In the framework of a multiscale modeling approach, we present a systematic study of a bipolar rectifying nanopore using a continuum and a particle simulation method. The common ground in the two methods is the application of the Nernst-Planck (NP) equation to compute ion transport in the framework of the implicit-water electrolytemodel. The difference is that the Poisson-Boltzmann theory is used in the Poisson-Nernst-Planck (PNP) approach, while the Local Equilibrium Monte Carlo (LEMC) method is used in the particle simulation approach (NP+LEMC) to relate the concentration profile to the electrochemical potential profile. Since we consider a bipolar pore which is short and narrow, we perform simulations using two-dimensional PNP. In addition, results of a non-linear version of PNP that takes crowding of ions into account are shown. We observe that the mean field approximation applied in PNP is appropriate to reproduce the basic behavior of the bipolar nanopore (e.g., rectification) for varying parameters of the system (voltage, surface charge,electrolyte concentration, and pore radius). We present current data that characterize the nanopore's behavior as a device, as well as concentration, electrical potential, and electrochemical potential profiles.

Effective short-range Coulomb correction to model the aggregation behavior of ionic surfactants

We present a short-range correction to the Coulomb potential to investigate the aggregation of amphiphilic molecules in aqueous solutions. The proposed modification allows to quantitatively reproduce the distribution of counterions above the critical micelle concentration (CMC) or, equivalently, the degree of ionization, α, of the micellar clusters. In particular, our theoretical framework has been applied to unveil the behavior of the cationic surfactant C24H49N2O2 (+) CH3SO4 (-), which offers a wide range of applications in the thriving and growing personal care market. A reliable and unambiguous estimation of α is essential to correctly understand many crucial features of the micellar solutions, such as their viscoelastic behavior and transport properties, in order to provide sound formulations for the above mentioned personal care solutions. We have validated our theory by performing extensive lattice Monte Carlo simulations, which show an excellent agreement with experimental observations. More specifically, our coarse-grained model is able to reproduce and predict the complex morphology of the micelles observed at equilibrium. Additionally, our simulation results disclose the existence of a transition from a monodisperse to a bidisperse size distribution of aggregates, unveiling the intriguing existence of a second CMC.

Hybrid Monte Carlo and continuum modeling of electrolytes with concentration-induced dielectric variations

The distribution of ions near a charged surface is an important quantity in many biological and material processes, and has been therefore investigated intensively. However, few theoretical and simulation approaches have included the influence of concentration-induced variations in the local dielectric permittivity of an underlying electrolyte solution. Such local variations have long been observed and known to affect the properties of ionic solution in the bulk and around the charged surface. We propose a hybrid computational model that combines Monte Carlo simulations with continuum electrostatic modeling to investigate such properties. A key component in our hybrid model is a semianalytical formula for the ion-ion interaction energy in a dielectrically inhomogeneous environment. This formula is obtained by solving for the Green's function Poisson's equation with ionic-concentration-dependent dielectric permittivity using a harmonic interpolation method and spherical harmonic series. We also construct a self-consistent continuum model of electrostatics to describe the effect of ionic-concentration-dependent dielectric permittivity and the resulting self-energy contribution. With extensive numerical simulations, we verify the convergence of our hybrid simulation scheme, show the qualitatively different structures of ionic distribution due to the concentration-induced dielectric variations, and compare our simulation results with the self-consistent continuum model. In particular, we study the differences between weakly and strongly charged surfaces and multivalencies of counterions. Our hybrid simulations conform particularly the depletion of ionic concentrations near a charged surface and also capture the charge inversion. We discuss several issues and possible further improvement of our approach for simulations of large charged systems.

Cation Selectivity in Biological Cation Channels Using Experimental Structural Information and Statistical Mechanical Simulation

Cation selective channels constitute the gate for ion currents through the cell membrane. Here we present an improved statistical mechanical model based on atomistic structural information, cation hydration state and without tuned parameters that reproduces the selectivity of biological Na+ and Ca2+ ion channels. The importance of the inclusion of step-wise cation hydration in these results confirms the essential role partial dehydration plays in the bacterial Na+ channels. The model, proven reliable against experimental data, could be straightforwardly used for designing Na+ and Ca2+ selective nanopores.

Dielectric response of the water hydration layer around spherical solutes

We calculate the local dielectric function ɛ(r) inside the hydration layer around a spherical solute (i) from molecular dynamics simulations including explicit solutes and (ii) theoretically using the nonlocal dielectric function of bulk water which includes the radial electric field, but not the explicit solute. The observed agreement between the two concepts shows that while ɛ(r) is strongly different from bulk, this difference is not due to restructuring of the hydrogen bond network but is mostly a consequence of the field geometry. The dielectric response differs for anions and cations, yielding a natural explanation for the well-known charge asymmetry of ionic solvation in agreement with experimental data.

Electrolytes between dielectric charged surfaces: Simulations and theory

We present a simulation method to study electrolyte solutions in a dielectric slab geometry using a modified 3D Ewald summation. The method is fast and easy to implement, allowing us to rapidly resum an infinite series of image charges. In the weak coupling limit, we also develop a mean-field theory which allows us to predict the ionic distribution between the dielectric charged plates. The agreement between both approaches, theoretical and simulational, is very good, validating both methods. Examples of ionic density profiles in the strong electrostatic coupling limit are also presented. Finally, we explore the confinement of charge asymmetric electrolytes between neutral surfaces.

Iterative spectral method for solving electrostatic or magnetostatic problems in complex and evolving heterostructures

Electrostatic or magnetostatic problems involving complex heterogeneity are nontrivial for modeling and simulation. Most existing numerical methods focus on sharp interface models and the computational cost increases with increasing complexity of the geometry. Here we develop an iterative spectral method, the bound charge successive approximation algorithm, to solve electrostatic or magnetostatic heterogeneity problems in the context of diffuse-interface modeling. As tests and verifications, this algorithm is applied to calculation of the depolarization factor of an ellipsoid, and simulation of random dielectric mixtures and the dielectophoretic motion of multiple particles. The algorithm shows excellent efficiency and the computational cost mainly depends on the permittivity or permeability contrast in the whole system, regardless of the complexity of the geometry. In particular, for evolving heterostructures the solution of bound charge in one time step can be used as input for the next, which could further significantly shorten the iteration (approximation) process, making it practical to simulate the long-range electrostatic or magnetostatic interaction in complex and evolving heterostructures.

Dielectric effects in the self-assembly of binary colloidal aggregates

Electrostatic interactions play an important role in numerous self-assembly phenomena, including colloidal aggregation. Although colloids typically have a dielectric constant that differs from the surrounding solvent, the effective interactions that arise from inhomogeneous polarization charge distributions are generally neglected in theoretical and computational studies. We introduce an efficient technique to resolve polarization charges in dynamical dielectric geometries, and demonstrate that dielectric effects qualitatively alter the predicted self-assembled structures, with surprising colloidal strings arising from many-body effects.

Singular perturbation solutions of steady-state Poisson-Nernst-Planck systems

We study the Poisson-Nernst-Planck (PNP) system with an arbitrary number of ion species with arbitrary valences in the absence of fixed charges. Assuming point charges and that the Debye length is small relative to the domain size, we derive an asymptotic formula for the steady-state solution by matching outer and boundary layer solutions. The case of two ionic species has been extensively studied, the uniqueness of the solution has been proved, and an explicit expression for the solution has been obtained. However, the case of three or more ions has received significantly less attention. Previous work has indicated that the solution may be nonunique and that even obtaining numerical solutions is a difficult task since one must solve complicated systems of nonlinear equations. By adopting a methodology that preserves the symmetries of the PNP system, we show that determining the outer solution effectively reduces to solving a single scalar transcendental equation. Due to the simple form of the transcendental equation, it can be solved numerically in a straightforward manner. Our methodology thus provides a standard procedure for solving the PNP system and we illustrate this by solving some practical examples. Despite the fact that for three ions, previous studies have indicated that multiple solutions may exist, we show that all except for one of these solutions are unphysical and thereby prove the existence and uniqueness for the three-ion case.

Generating true minima in constrained variational formulations via modified Lagrange multipliers

Variational principles are important in the investigation of large classes of physical systems. They can be used both as analytical methods as well as starting points for the formulation of powerful computational techniques such as dynamical optimization methods. Systems with charged objects in dielectric media and systems with magnetically active particles are important examples. In these examples and other important cases, the variational principles describing the system are required to obey a number of constraints. These constraints are implemented within the variational formulation by means of Lagrange multipliers. Such constrained variational formulations are in general not unique. For the application of efficient simulation methods, one must find specific formulations that satisfy a number of important conditions. An often required condition is that the functional be positive-definite, in other words, its extrema be actual minima. In this article, we present a general approach to attack the problem of finding, among equivalent variational functionals, those that generate true minima. The method is based on the modification of the Lagrange multiplier which allows us to generate large families of effective variational formulations associated with a single original constrained variational principle. We demonstrate its application to different examples and, in particular, to the important cases of Poisson and Poisson-Boltzmann equations. We show how to obtain variational formulations for these systems with extrema that are always minima.

Ion permeation in the NanC porin from Escherichia coli: free energy calculations along pathways identified by coarse-grain simulations

Using the X-ray structure of a recently discovered bacterial protein, the N-acetylneuraminic acid-inducible channel (NanC), we investigate computationally K(+) and Cl(-) ions' permeation. We identify ion permeation pathways that are likely to be populated using coarse-grain Monte Carlo simulations. Next, we use these pathways as reaction coordinates for umbrella sampling-based free energy simulations. We find distinct tubelike pathways connecting specific binding sites for K(+) and, more pronounced, for Cl(-) ions. Both ions permeate the porin preserving almost all of their first hydration shell. The calculated free energy barriers are G(#) ≈ 4 kJ/mol and G(#) ≈ 8 kJ/mol for Cl(-) and K(+), respectively. Within the approximations associated with these values, discussed in detail in this work, we suggest that the porin is slightly selective for Cl(-) versus K(+). Our suggestion is consistent with the experimentally observed weak Cl(-) over K(+) selectivity. A rationale for the latter is suggested by a comparison with previous calculations on strongly anion selective porins.

Electrostatic interaction in the presence of dielectric interfaces and polarization-induced like-charge attraction

Electrostatic polarization is important in many nano- and micro-scale physical systems such as colloidal suspensions, biopolymers, and nanomaterials assembly. The calculation of polarization potential requires an efficient algorithm for solving 3D Poisson's equation. We have developed a useful image charge method to rapid evaluation of the Green's function of the Poisson's equation in the presence of spherical dielectric discontinuities. This paper presents an extensive study of this method by giving a convergence analysis and developing a coarse-graining algorithm. The use of the coarse graining could reduce the number of image charges to around a dozen, by 1-2 orders of magnitude. We use the algorithm to investigate the interaction force between likely charged spheres in different dielectric environments. We find the size and charge asymmetry leads to an attraction between like charges, in agreement with existing results. Furthermore, we study three-body interactions and find, in the presence of an external interface, that the interaction force depends on the curvature of the interface and performs a nonmonotonic electrostatic force.

Simulation of charged systems in heterogeneous dielectric media via a true energy functional

For charged systems in heterogeneous dielectric media, a key obstacle for molecular dynamics (MD) simulations is the need to solve the Poisson equation in the media. This obstacle can be bypassed using MD methods that treat the local polarization charge density as a dynamic variable, but such approaches require access to a true free energy functional, one that evaluates to the equilibrium electrostatic energy at its minimum. In this Letter, we derive the needed functional. As an application, we develop a Car-Parrinello MD method for the simulation of free charges present near a spherical emulsion droplet separating two immiscible liquids with different dielectric constants. Our results show the presence of nonmonotonic ionic profiles in the dielectric with a lower dielectric constant.