Search

Supercapacitors are promising electrochemical energy storage devices due to their prominent performance in rapid charging/discharging rates, long cycle life, stability, etc. Experimental measurement and theoretical prediction on charging timescale for supercapacitors often have large difference. This work develops a matched asymptotic expansion method to derive the charging dynamics of supercapacitors with porous electrodes, in which the supercapacitors are described by the stack-electrode model. Coupling leading-order solutions between every two stacks by continuity of ionic concentration and fluxes leads to an ODE system, which is a generalized equivalent circuit model for zeta potentials, with the potential-dependent nonlinear capacitance and resistance determined by physical parameters of electrolytes, e.g., specific counterion valences for asymmetric electrolytes. Linearized stability analysis on the ODE system after projection is developed to theoretically characterize the charging timescale. The derived asymptotic solutions are numerically verified. Further numerical investigations on the biexponential charging timescales demonstrate that the proposed generalized equivalent circuit model, as well as companion linearized stability analysis, can faithfully capture the charging dynamics of symmetric/asymmetric electrolytes in supercapacitors with porous electrodes., Comment: 26 pages, 6 figures

In this paper we propose and analyze a finite difference numerical scheme for the Poisson-Nernst-Planck equation (PNP) system. To understand the energy structure of the PNP model, we make use of the Energetic Variational Approach (EnVarA), so that the PNP system could be reformulated as a non-constant mobility H − 1 H^{-1} gradient flow, with singular logarithmic energy potentials involved. To ensure the unique solvability and energy stability, the mobility function is explicitly treated, while both the logarithmic and the electric potential diffusion terms are treated implicitly, due to the convex nature of these two energy functional parts. The positivity-preserving property for both concentrations, n n and p p , is established at a theoretical level. This is based on the subtle fact that the singular nature of the logarithmic term around the value of 0 0 prevents the numerical solution reaching the singular value, so that the numerical scheme is always well-defined. In addition, an optimal rate convergence analysis is provided in this work, in which many highly non-standard estimates have to be involved, due to the nonlinear parabolic coefficients. The higher order asymptotic expansion (up to third order temporal accuracy and fourth order spatial accuracy), the rough error estimate (to establish the ℓ ∞ \ell ^\infty bound for n n and p p ), and the refined error estimate have to be carried out to accomplish such a convergence result. In our knowledge, this work will be the first to combine the following three theoretical properties for a numerical scheme for the PNP system: (i) unique solvability and positivity, (ii) energy stability, and (iii) optimal rate convergence. A few numerical results are also presented in this article, which demonstrates the robustness of the proposed numerical scheme.

This work proposes a new variational, thermodynamically consistent model to predict thermal electrokinetics in electric double layer capacitors (EDLCs) by using an energetic variational approach. The least action principle and maximum dissipation principle from the non-equilibrium thermodynamics are employed to develop modified Nernst-Planck equations for non-isothermal ion transport with temperature inhomogeneity. Laws of thermodynamics are employed to derive a temperature evolution equation with heat sources due to thermal pressure and electrostatic interactions. Numerical simulations successfully predict temperature oscillation in the charging-discharging processes of EDLCs, indicating that the developed model is able to capture reversible and irreversible heat generations. The impact of ionic sizes and scan rate of surface potential on ion transport, heat generation, and charge current is systematically assessed in cyclic voltammetry simulations. It is found that the thermal electrokinetics in EDLCs cannot follow the surface potential with fast scan rates, showing delayed dynamics with hysteresis diagrams. Our work thus provides a useful tool for physics-based prediction of thermal electrokinetics in EDLCs.

Charge dynamics play essential role in many practical applications such as semiconductors, electrochemical devices and transmembrane ion channels. A Maxwell-Ampère Nernst-Planck (MANP) model that describes charge dynamics via concentrations and the electric displacement is able to take effects beyond mean-field approximations into account. To obtain physically faithful numerical solutions, we develop a structure-preserving numerical method for the MANP model whose solution has several physical properties of importance. By the Slotboom transform with entropic-mean approximations, a positivity preserving scheme with Scharfetter-Gummel fluxes is derived for the generalized Nernst-Planck equations. To deal with the curl-free constraint, the dielectric displacement from the Maxwell-Ampère equation is further updated with a local relaxation algorithm of linear computational complexity. We prove that the proposed numerical method unconditionally preserves the mass conservation and the solution positivity at the discrete level, and satisfies the discrete energy dissipation law with a time-step restriction. Numerical experiments verify that our numerical method has expected accuracy and structure-preserving properties. Applications to ion transport with large convection, arising from boundary-layer electric field and Born solvation interactions, further demonstrate that the MANP formulation with the proposed numerical scheme has attractive performance and can effectively describe charge dynamics with large convection of high numerical cell Péclet numbers.

Water fluctuates in a hydrophobic confinement, forming multiple dry and wet hydration states through evaporation and condensation. Transitions between such states are critical to both thermodynamics and kinetics of solute molecular processes, such as protein folding and protein–ligand binding and unbinding. To efficiently predict such dry–wet transition paths, we develop a hybrid approach that combines a variational implicit solvation model, a generalized string method for minimum free-energy paths, and the level-set numerical implementation. This approach is applied to three molecular systems: two hydrophobic plates, a carbon nanotube, and a synthetic host molecule Cucurbit[7]uril. Without an explicit description of individual water molecules, our mesoscale approach effectively captures multiple dry and wet hydration states, multiple dry–wet transition paths, such as those geometrically symmetric and asymmetric paths, and transition states, providing activation energy barriers between different states. Further analysis shows that energy barriers depend on mesoscopic lengths, such as the separation distance between the two plates and the cross section diameter of the nanotube, and that the electrostatic interactions strongly influence the dry–wet transitions. With the inclusion of solute atomic motion, general collective variables as reaction coordinates, and the finite-temperature string method, together with an improved treatment of continuum electrostatics, our approach can be further developed to sample an ensemble of transition paths, providing more accurate predictions of the transition kinetics.

In a continuum model of the solvation of charged molecules in an aqueous solvent, the classical Poisson-Boltzmann (PB) theory is generalized to include the solute point charges and the dielectric boundary that separates the high-dielectric solvent from the low-dielectric solutes. With such a setting, we construct an effective electrostatic free-energy functional of ionic concentrations, where the solute point charges are regularized by a reaction field. We prove that such a functional admits a unique minimizer in a class of admissible ionic concentrations and that the corresponding electrostatic potential is the unique solution to the boundary-value problem of the dielectric-boundary PB equation. The negative first variation of this minimum free energy with respect to variations of the dielectric boundary defines the normal component of the dielectric boundary force. Together with the solute-solvent interfacial tension and van der Waals interaction forces, such boundary force drives an underlying charged molecular system to a stable equilibrium, as described by a variational implicit-solvent model. We develop an $L^2$-theory for the continuity and differentiability of solutions to elliptic interface problems with respect to boundary variations, and derive an explicit formula of the dielectric boundary force. With a continuum description, our result of the dielectric boundary force confirms a molecular-level prediction that the electrostatic force points from the high-dielectric and polarizable aqueous solvent to the charged molecules. Our method of analysis is general as it does not rely on any variational principles., 52 pages, 2 figures

In this work, we design and analyze a conservative, positivity preserving, and free energy dissipating finite difference method for the multi-dimensional nonlocal Fokker–Planck (FP) equation. Based on a non-logarithmic Landau transformation, a central-differencing spatial discretization using harmonic-mean approximations is developed. Both forward and backward Euler discretizations in time are employed to derive an explicit scheme and a linearized semi-implicit scheme, respectively. Three desired properties that are possessed by analytical solutions: i) mass conservation, ii) free-energy dissipation, and iii) positivity, are proved to be maintained at discrete level. Remarkably, numerical analysis demonstrates that the semi-implicit time discretization ensures the property of positivity preserving unconditionally. Due to the advantages brought by the harmonic-mean approximations, an estimate on the upper bound of the condition number of the resulting coefficient matrix is further established for the semi-implicit scheme. Extensive numerical tests are performed to validate aforementioned properties numerically.

We develop a computational model for the solvation of charged molecules in an aqueous solvent. Such solvent is modeled as an incompressible fluid, and its dynamics is described by the Stokes equations. All the surface tension force, van der Waals dispersive force, electrostatic force, hydrostatic pressure, and viscous force are balanced on the solute-solvent interface, giving rise to the traction boundary conditions on such an interface. We use the level-set method for the solute-solvent interface motion. To describe the hydrophobic interaction that is characterized by the volume change of the solute region, we design special numerical boundary conditions for the Stokes equations. We then reformulate and discretize the Stokes equations with a second-order finite difference scheme using a MAC grid. At virtual nodes near the solute-solvent interface, we use interpolation to discretize the boundary conditions and formulate the difference equations at such nodes. The resulting linear system is reduced by the Schur complement and then solved by the least-squares method with various techniques of acceleration. Numerical tests show that our method has a second-order convergence in the maximum norm for both the velocity and pressure up to the boundary. We apply our methods to some model molecular systems. With the comparison of our approach with a well established, static solvation theory and method, our dynamic model and numerical techniques accurately reproduce the solvation free energies, and capture the dry and wet molecular states. Our work provides a possibility of describing solvent fluctuations through the solvent fluid flow.

Metal–Intermetallic Laminate (MIL) composites are laminate structures that are fabricated by optimizing unique benefits of constituent components to have attractive physical and mechanical features, such as high strength, stiffness, and toughness. The synthesis of MIL composites involves annealing reaction of multiple metallic elements at high temperatures. In this work, we propose, analyze, and implement a second-order semi-implicit scheme for the annealing of a ternary metallic reaction system. The robustness of such a numerical scheme is demonstrated by various numerical experiments, as well as the expected accuracy tests. At the theoretical level, we provide a detailed convergence analysis, which confirms the second-order accuracy in both the temporal and spatial discretization. Moreover, this numerical scheme is used to study the annealing process of the Al–Fe–Ni ternary system. The computation reproduces the formation of the Al-rich and Ni-rich layers in the annealing process. We present the computational results as diffusion paths in a ternary phase diagram, with a nice agreement with that of experimental data. We also study the growth kinetics of the two layers by calculating the rate constant and kinetic exponent of the reaction. The morphology of interfaces between different layers is thoroughly investigated as well. The computational results indicate that the numerical scheme is an effective, useful tool for predicting the microstructure evolution in the annealing process of a ternary reaction system.

In this work, we numerically study linear stability of multiple steady-state solutions to a type of steric Poisson--Nernst--Planck (PNP) equations with Dirichlet boundary conditions, which are applicable to ion channels. With numerically found multiple steady-state solutions, we obtain $S$-shaped current-voltage and current-concentration curves, showing hysteretic response of ion conductance to voltages and boundary concentrations with memory effects. Boundary value problems are proposed to locate bifurcation points and predict the local bifurcation diagram near bifurcation points on the $S$-shaped curves. Numerical approaches for linear stability analysis are developed to understand the stability of the steady-state solutions that are only numerically available. Finite difference schemes are proposed to solve a derived eigenvalue problem involving differential operators. The linear stability analysis reveals that the $S$-shaped curves have two linearly stable branches of different conductance levels and one linearly unstable intermediate branch, exhibiting classical bistable hysteresis. As predicted in the linear stability analysis, transition dynamics, from a steady-state solution on the unstable branch to a one on the stable branches, are led by perturbations associated to the mode of the dominant eigenvalue. Further numerical tests demonstrate that the finite difference schemes proposed in the linear stability analysis are second-order accurate. Numerical approaches developed in this work can be applied to study linear stability of a class of time-dependent problems around their steady-state solutions that are computed numerically.

In this work, we develop novel structure-preserving numerical schemes for a class of nonlinear Fokker--Planck equations with nonlocal interactions. Such equations can cover many cases of importance, such as porous medium equations with external potentials, optimal transport problems, and aggregation-diffusion models. Based on the Energetic Variational Approach, a trajectory equation is first derived by using the balance between the maximal dissipation principle and least action principle. By a convex-splitting technique, we propose energy dissipating numerical schemes for the trajectory equation. Rigorous numerical analysis reveals that the nonlinear numerical schemes are uniquely solvable, naturally respect mass conservation and positivity at fully discrete level, and preserve steady states. Under certain smoothness assumptions, the numerical schemes are shown to be second order accurate in space and first order accurate in time. Extensive numerical simulations are performed to demonstrate several valuable features of the proposed schemes. In addition to the preservation of physical structures, such as positivity, mass conservation, discrete energy dissipation, blue and steady states, numerical simulations further reveal that our numerical schemes are capable of solving \emph{degenerate} cases of the Fokker--Planck equations effectively and robustly. It is shown that the developed numerical schemes have convergence order even in degenerate cases with the presence of solutions having compact support, can accurately and robustly compute the waiting time of free boundaries without any oscillation, and can approximate blow-up singularity up to machine precision.

In this work, we consider numerical methods for the Poisson–Nernst–Planck–Cahn–Hilliard (PNPCH) equations with steric interactions, which correspond to a non-constant mobility H − 1 gradient flow of a free-energy functional that consists of electrostatic free energies, steric interaction energies of short range, entropic contribution of ions, and concentration gradient energies. We propose a novel energy stable numerical scheme that respects mass conservation and positivity at the discrete level. Existence and uniqueness of the solution to the proposed nonlinear scheme are established by showing that the solution is a unique minimizer of a convex functional over a closed, convex domain. The positivity of numerical solutions is further theoretically justified by the singularity of the entropy terms, which prevents the minimizer from approaching zero concentrations. A further numerical analysis proves discrete free-energy dissipation. Extensive numerical tests are performed to validate that the numerical scheme is first-order accurate in time and second-order accurate in space, and is capable of preserving the desired properties, such as mass conservation, positivity, and free energy dissipation, at the discrete level. Moreover, the PNPCH equations and the proposed scheme are applied to study charge dynamics and self-assembled nanopatterns in highly concentrated electrolytes that are widely used in electrochemical energy devices. Numerical results demonstrate that the PNPCH equations and our numerical scheme are able to capture nanostructures, such as lamellar patterns and labyrinthine patterns in electric double layers and the bulk, and multiple time relaxation with multiple time scales. The multiple time relaxation dynamics with metastability take long time to reach an equilibrium, highlighting the need for robust, energy stable numerical schemes that allow large time stepping. In addition, we numerically characterize the interplay between cross steric interactions of short range and the concentration gradient regularization, and their impact on the development of nanostructures in the equilibrium state.

Ion transport, often described by the Poisson--Nernst--Planck (PNP) equations, is ubiquitous in electrochemical devices and many biological processes of significance. In this work, we develop conservative, positivity-preserving, energy dissipating, and implicit finite difference schemes for solving the multi-dimensional PNP equations with multiple ionic species. A central-differencing discretization based on harmonic-mean approximations is employed for the Nernst--Planck (NP) equations. The backward Euler discretization in time is employed to derive a fully implicit nonlinear system, which is efficiently solved by a newly proposed Newton's method. The improved computational efficiency of the Newton's method originates from the usage of the electrostatic potential as the iteration variable, rather than the unknowns of the nonlinear system that involves both the potential and concentration of multiple ionic species. Numerical analysis proves that the numerical schemes respect three desired analytical properties (conservation, positivity preserving, and energy dissipation) fully discretely. Based on advantages brought by the harmonic-mean approximations, we are able to establish estimate on the upper bound of condition numbers of coefficient matrices in linear systems that are solved iteratively. The solvability and stability of the linearized problem in the Newton's method are rigorously established as well. Numerical tests are performed to confirm the anticipated numerical accuracy, computational efficiency, and structure-preserving properties of the developed schemes. Adaptive time stepping is implemented for further efficiency improvement. Finally, the proposed numerical approaches are applied to characterize ion transport subject to a sinusoidal applied potential.

The charge transport in an environment with inhomogeneous dielectric permittivity is ubiquitous in many areas such as electrochemical energy devices and biophysical systems. We theoretically study ...

The Poisson--Nernst--Planck (PNP) equations have been widely applied to describe ionic transport in ion channels, nanofluidic devices, and many electrochemical systems. Despite their wide applications, the PNP equations fail in predicting dynamics and equilibrium states of ionic concentrations in confined environments, due to the ignorance of the excluded volume effect. In this work, a simple but effective modified PNP (MPNP) model with the excluded volume effect is derived, based on a modification of diffusion coefficients of ions. At the steady state, a modified Poisson--Boltzmann (MPB) equation is obtained with the help of the Lambert-W special function. The existence and uniqueness of a weak solution to the MPB equation are established. Further analysis on the limit of weak and strong electrostatic potential leads to two modified Debye screening lengths, respectively. A numerical scheme that conserves total ionic concentration and satisfies energy dissipation is developed for the MPNP model. Numerical analysis is performed to prove that our scheme respects ionic mass conservation and satisfies a corresponding discrete free energy dissipation law. Positivity of numerical solutions is also discussed and numerically investigated. Numerical tests are conducted to demonstrate that the scheme is of second-order accurate in spatial discretization and has expected properties. Extensive numerical simulations reveal that the excluded volume effect has pronounced impacts on the dynamics of ionic concentration and flux. In addition, the effect of volume exclusion on the timescales of charge diffusion is systematically investigated by studying the evolution of free energies and diffuse charges.

Ligand-receptor binding and unbinding are fundamental biomolecular processes and particularly essential to drug efficacy. Environmental water fluctuations, however, impact the corresponding thermodynamics and kinetics and thereby challenge theoretical descriptions. Here, we devise a holistic, implicit-solvent, multi-method approach to predict the (un)binding kinetics for a generic ligand-pocket model. We use the variational implicit-solvent model (VISM) to calculate the solute-solvent interfacial structures and the corresponding free energies, and combine the VISM with the string method to obtain the minimum energy paths and transition states between the various metastable ('dry' and 'wet') hydration states. The resulting dry-wet transition rates are then used in a spatially-dependent multi-state continuous-time Markov chain Brownian dynamics simulations, and the related Fokker-Planck equation calculations, of the ligand stochastic motion, providing the mean first-passage times for binding and unbinding. We find the hydration transitions to significantly slow down the binding process, in semi-quantitative agreement with existing explicit-water simulations, but significantly accelerate the unbinding process. Moreover, our methods allow the characterization of non-equilibrium hydration states of pocket and ligand during the ligand movement, for which we find substantial memory and hysteresis effects for binding versus unbinding. Our study thus provides a significant step forward towards efficient, physics-based interpretation and predictions of the complex kinetics in realistic ligand-receptor systems., 6 pages, 5 figures, accepted for publication in Proc. Natl. Acad. Sci. (PNAS)

Combining classical density functional theory with modified interfacial statistical associating fluid theory, we investigate the surface electrostatic force (SEF) between two negatively and similarly charged planar surfaces immersed in electrolyte solution consisting of monomer anion and dimer cation. The dimer cation is comprised of two hard spheres (HSs) tangentially bonded together and freely rotating on the top of each other; one of the two HS sites is positively charged and the other is neutral. New findings about novel properties of the SEF caused by the dimer counter-ion are summarized below. (i) In the presence of only univalent counter-ion, the SEF always becomes more and more repulsive with increasing of the dimer counter-ion neutral site size, whether the bulk electrolyte solution is dense or not, and the surface charge strength is small, moderate, or large; repulsion strength of the SEF always reduces with the bulk concentration, and at equal surface distance rapidly increases with the surface charge strength. (ii) In the presence of bivalent counter-ion of monomer or dimer, rather complex SEF behavior emerges. Specifically, (a) at low surface charge strength, the dimer counter-ion with small neutral site size induces one obvious like-charge attraction (LCA) at lower bulk concentrations; whereas at higher bulk concentrations, a repulsive SEF emerges instead, whose strength monotonously decreases with the neutral site size. (b) For high surface charge strengths, the SEF curves exhibit oscillatory structure with the LCA feature; moreover, the dimer counter-ion always induces larger LCA position and strength than the monomer counter-ion does; the LCA strength decreases with the counter-ion neutral site size, and the LCA position always reduces with increase of the surface charge strength. (c) Dependence of the LCA strength on the bulk concentration is conditional on the surface charge strengths considered; to speak specifically, for very high surface charge strength, the LCA strength reduces with the bulk concentration, whereas for not so high value, the opposite is true. The above observations can be analyzed from fluctuation and balance of the three factors: osmotic pressure difference between inside and outside the plates; hard sphere repulsion in congested space, and “hydrogen bond” formation.

430SS MIL composites were initially fabricated under semi-solid reaction at 655 °C and then annealed at high temperature (1000 °C) to obtain the ductile Fe-rich intermetallic phases. The microstructure evolution and micro-hardness distribution across the intermetallic/metal interfaces were investigated, as well as the experimentally and computationally evaluating the growth kinetics of the Fe-rich intermetallic layers. A multilayer structure (namely Al-rich layer, Fe-rich layer and B2 phase layer) with a complex mix of intermetallic phases (FeAl 2 , Cr 5 Al 8 and FeAl) form in the high temperature annealed 430SS MIL composites. The FeAl phase with a B2 structure formed in the intermetallic layer provides a more gradual transition of mechanical properties across the intermetallic/metal layer. The growth mechanism of the Fe-rich layer is interface controlled, while the growth of the B2 phase layer is diffusion controlled based on both the experimental and numerical analysis. A relationship between the processing-microstructure-properties of the MIL composites has been established and applied to predict and tailor the microhardness distribution across the multi-layered structure of the 430SS MIL composites.

We construct a mean-field variational model to study how the dependence of dielectric coefficient (i.e., relative permittivity) on local ionic concentrations affects the electrostatic interaction in an ionic solution near a charged surface. The electrostatic free-energy functional of ionic concentrations, which is the key object in our model, consists mainly of the electrostatic potential energy and the ionic ideal-gas entropy. The electrostatic potential is determined by Poisson’s equation in which the dielectric coefficient depends on the sum of concentrations of individual ionic species. This dependence is assumed to be qualitatively the same as that on the salt concentration for which experimental data are available and analytical forms can be obtained by the data fitting. We derive the first and second variations of the free-energy functional, obtain the generalized Boltzmann distributions, and show that the free-energy functional is in general nonconvex. To validate our mathematical analysis, we numerically minimize our electrostatic free-energy functional for a radially symmetric charged system. Our extensive computations reveal several features that are significantly different from a system modeled with a dielectric coefficient independent of ionic concentration. These include the non-monotonicity of ionic concentrations, the ionic depletion near a charged surface that has been previously predicted by a one-dimensional model, and the enhancement of such depletion due to the increase of surface charges or bulk ionic concentrations.

By using classical density functional theory we investigate influences of atomic scale strip shape surface charge separation on the effective electrostatic force (EEF) between two overall neutral surfaces immersed in 1:1 electrolyte solution. Main novel findings are confirmed, and summarized as follows: (i) At lower bulk concentration, the EEF is always attractive for asymmetrical surface charge pattern and repulsive for symmetrical surface charge pattern. Strength of the EEF is positively correlated with the strip shape domain width and surface charge density intensity, and generally smaller than that of the EEF between two equally and oppositely or similarly charged planar surfaces with same surface charge intensity. (ii) With increasing of the bulk concentration, the asymmetrical strip shape surface charge pattern induces two repulsion peaks, strengths of the repulsion peaks around one times and double the ion diameter surface separation are negatively and positively correlated with the strip shape domain width, respectively; the second repulsion peak position reduces with the domain charge strength. (iii) At higher bulk concentration, the symmetrical strip shape charge pattern induces two attraction wells of the EEF curve; one positioned at wall separation larger than one times and smaller than double the ion size, the other approximately 2.75 times the ion size. (iv) Considering that systems favor low-energy configurations, the present calculations provide clear evidence that effective repulsion between colloids with non-uniform surface charge distribution or charge fluctuation is less than that with uniform surface charge distribution or with fixed charges, respectively. All of the observations can be explained self-consistently by using concepts such as ion saturation adsorption capacity, arrangement compactness of the adsorption layer, formation of quasi-hydrogen-bonding, and charge reversal.

Classical density functional theory is combined with the extended primitive model and solvent primitive model to investigate how differential capacitance Cd of electrical double-layer formed inside cylindrical pore is influenced by solvent granularity, bulk concentration, counter-ion diameter, and neutral and non-hard sphere (HS) potential utail_αβ(r) between ion species, between solvent and counter-ion (and co-ion). Several main conclusions are drawn. (i) The Cd−surface charge strength (|σ|) curve generally rises as a result of consideration of the solvent granularity. (ii) An interionic attractive utail_αβ(r) helps in raising the Cd curve whereas a repulsive utail_αβ(r) tends to lower the Cd curve. Moreover, a secondary Cd peak appears when |σ| increases to an appropriately large value, and becomes increasingly evident with the strength of the attractive utail_αβ(r) and eventually disappears as the attractive utail_αβ(r)reduces in attraction or changes into repulsion at all. (iii) The Cd − |σ| curve is raised by increasing the repulsion between the solvent and counter-ion; however, both counter-ion diameter and bulk concentration cause obvious changes of the Cd − |σ∗| curve morphology. Concretely speaking, the Cd peak height and the peak position rise with the counter-ion size decreasing. Moreover, changing the utail_αβ(r) between solvent and counter-ion from attraction to repulsion facilitates a transition from bell-shaped curve to camel-shaped curve. (iv) Effects of all of the factors considered in influencing the Cd curve diminish increasingly with |σ|. All of the above observations can be explained by considering the interionic depletion potential induced by the solvent and its changes with the system parameters.

The steady-state Poisson-Nernst-Planck (ssPNP) equations are an effective model for the description of ionic transport in ion channels. It is observed that an ion channel exhibits voltage-dependent switching between open and closed states. Different conductance states of a channel imply that the ssPNP equations probably have multiple solutions with different level of currents. We propose numerical approaches to study multiple solutions to the ssPNP equations with multiple ionic species. To find complete current-voltage (I-V ) and current-concentration (I-C) curves, we reformulate the ssPNP equations into four different boundary value problems (BVPs). Numerical continuation approaches are developed to provide good initial guesses for iteratively solving algebraic equations resulting from discretization. Numerical continuations on V , I, and boundary concentrations result in S-shaped and double S-shaped (I-V and I-C) curves for the ssPNP equations with multiple species of ions. There are five solutions to the ssPNP equations with five ionic species, when an applied voltage is given in certain intervals. Remarkably, the current through ion channels responds hysteretically to varying applied voltages and boundary concentrations, showing a memory effect. In addition, we propose a useful computational approach to locate turning points of an I-V curve. With obtained locations, we are able to determine critical threshold values for hysteresis to occur and the interval for V in which the ssPNP equations have multiple solutions. Our numerical results indicate that the developed numerical approaches have a promising potential in studying hysteretic conductance states of ion channels.

© 2015 American Chemical Society. Protein-ligand binding is a key biological process at the molecular level. The identification and characterization of small-molecule binding sites on therapeutically relevant proteins have tremendous implications for target evaluation and rational drug design. In this work, we used the recently developed level-set variational implicit-solvent model (VISM) with the Coulomb field approximation (CFA) to locate and characterize potential protein-small-molecule binding sites. We applied our method to a data set of 515 protein-ligand complexes and found that 96.9% of the cocrystallized ligands bind to the VISM-CFA-identified pockets and that 71.8% of the identified pockets are occupied by cocrystallized ligands. For 228 tight-binding protein-ligand complexes (i.e, complexes with experimental pKdvalues larger than 6), 99.1% of the cocrystallized ligands are in the VISM-CFA-identified pockets. In addition, it was found that the ligand binding orientations are consistent with the hydrophilic and hydrophobic descriptions provided by VISM. Quantitative characterization of binding pockets with topological and physicochemical parameters was used to assess the "ligandability" of the pockets. The results illustrate the key interactions between ligands and receptors and can be very informative for rational drug design.

Mean-field methods, such as the Poisson-Boltzmann equation (PBE), are often used to calculate the electrostatic properties of molecular systems. In the past two decades, an enhancement of the PBE, the size-modified Poisson-Boltzmann equation (SMPBE), has been reported. Here, the PBE and the SMPBE are reevaluated for realistic molecular systems, namely, lipid bilayers, under eight different sets of input parameters. The SMPBE appears to reproduce the molecular dynamics simulation results better than the PBE only under specific parameter sets, but in general, it performs no better than the Stern layer correction of the PBE. These results emphasize the need for careful discussions of the accuracy of mean-field calculations on realistic systems with respect to the choice of parameters and call for reconsideration of the cost-efficiency and the significance of the current SMPBE formulation.

To study ion transport in electrolyte solutions, we propose numerical methods for a modified Poisson–Nernst–Planck model with ionic steric effects (SPNP). Positivity preserving schemes based on harmonic-mean approximations are proposed on a nonuniform mesh for the spatial discretization of the SPNP equations. Both explicit and semi-implicit discretization are considered in time. Numerical analysis shows that explicit forward Euler and semi-implicit trapezoidal discretization lead to schemes that maintain fully discrete solution positivity under a constraint on a mesh ratio, while the semi-implicit backward Euler discretization maintain fully discrete solution positivity unconditionally. We further study the condition numbers of the matrix associated with the semi-implicit backward Euler discretization, and establish an upper bound on condition numbers, indicating that the developed discretization based on harmonic-mean approximations effectively solves the issue of the presence of large condition numbers when using the Slotboom-type variables. Further numerical simulations confirm the analysis results on accuracy, positivity, and bounded condition numbers. The proposed schemes are also applied to study practical applications, such as the impact of self and cross steric interactions on ion distribution and rectifying behavior in a slit-shaped nanopore with surface charges. Possible extensions of the numerical methods to other modified PNP models with ion correlations and steric effects are also discussed.

The synthetic host cucurbit[7]uril (CB[7]) binds aromatic guests or metal complexes with ultrahigh affinity compared with that typically displayed in protein-ligand binding. Due to its small size, CB[7] serves as an ideal receptor-ligand system for developing computational methods for molecular recognition. Here, we apply the recently developed variational implicit-solvent model (VISM), numerically evaluated by the level-set method, to study hydration effects in the high-affinity binding of the B2 bicyclo[2.2.2]octane derivative to CB[7]. For the unbound host, we find that the host cavity favors the hydrated state over the dry state due to electrostatic effects. For the guest binding, we find reasonable agreement to experimental binding affinities. Dissection of the individual VISM free-energy contributions shows that the major driving forces are water-mediated hydrophobic interactions and the intrinsic (vacuum) host-guest van der Waals interactions. These findings are in line with recent experiments and molecular dynamics simulations with explicit solvent. It is expected that the level-set VISM, with further refinement on the electrostatic descriptions, can efficiently predict molecular binding and recognition in a wide range of future applications.

Diffusion is often the rate-determining step in many biological processes. Currently, the two main computational methods for studying diffusion are stochastic methods, such as Brownian dynamics, and continuum methods, such as the finite element method. A previous study introduced a new hybrid diffusion method that couples the strengths of each of these two methods, but was limited by the lack of interactions among the particles; the force on each particle had to be from an external field. This study further develops the method to allow charged particles. The method is derived for a general multidimensional system and is presented using a basic test case for a one-dimensional linear system with one charged species and a radially symmetric system with three charged species.

We design and implement numerical methods for the incompressible Stokes solvent flow and solute---solvent interface motion for nonpolar molecules in aqueous solvent. The balance of viscous force, surface tension, and van der Waals type dispersive force leads to a traction boundary condition on the solute---solvent interface. To allow the change of solute volume, we design special numerical boundary conditions on the boundary of a computational domain through a consistency condition. We use a finite difference ghost fluid scheme to discretize the Stokes equation with such boundary conditions. The method is tested to have a second-order accuracy. We combine this ghost fluid method with the level-set method to simulate the motion of the solute---solvent interface that is governed by the solvent fluid velocity. Numerical examples show that our method can predict accurately the blow up time for a test example of curvature flow and reproduce the polymodal (e.g., dry and wet) states of hydration of some simple model molecular systems.

The solute-solvent interface that separates biological molecules from their surrounding aqueous solvent characterizes the conformation and dynamics of such molecules. In this work, we construct a solvent fluid dielectric boundary model for the solvation of charged molecules and apply it to study the stability of a model cylindrical solute-solvent interface. The motion of the solute-solvent interface is defined to be the same as that of solvent fluid at the interface. The solvent fluid is assumed to be incompressible and is described by the Stokes equation. The solute is modeled simply by the ideal-gas law. All the viscous force, hydrostatic pressure, solute-solvent van der Waals interaction, surface tension, and electrostatic force are balanced at the solute-solvent interface. We model the electrostatics by Poisson's equation in which the solute-solvent interface is treated as a dielectric boundary that separates the low-dielectric solute from the high-dielectric solvent. For a cylindrical geometry, we find multiple cylindrically shaped equilibrium interfaces that describe polymodal (e.g., dry and wet) states of hydration of an underlying molecular system. These steady-state solutions exhibit bifurcation behavior with respect to the charge density. For their linearized systems, we use the projection method to solve the fluid equation and find the dispersion relation. Our asymptotic analysis shows that, for large wavenumbers, the decay rate is proportional to wavenumber with the proportionality half of the ratio of surface tension to solvent viscosity, indicating that the solvent viscosity does affect the stability of a solute-solvent interface. Consequences of our analysis in the context of biomolecular interactions are discussed.

In this paper, a linearly semi-implicit compact scheme is developed for the Burgers-Huxley equation. The equation is decomposed into two subproblems, i.e. a Burgers equation and a nonlinear ODE, by the operator splitting technique. The Burgers equation is solved by a linearly self-starting compact scheme which is fourth-order accurate in space and second-order accurate in time. The nonlinear ODE is discretized by a third-order semi-implicit Runge-Kutta method, which possesses good numerical stability with low computational cost. The numerical experiments show that the scheme provides the expected convergence order. Finally, several experiments are conducted to simulate the solutions of the Burgers-Huxley equation to validate our numerical method.

The numerical simulation of coupled nonlinear Schrodinger equations on unbounded domains is considered in this paper. By using the operator splitting technique, the original problem is decomposed into linear and nonlinear subproblems in a small time step. The linear subproblem turns out to be two decoupled linear Schrodinger equations on unbounded domains, where artificial boundaries are introduced to truncate the unbounded physical domains into finite ones. Local absorbing boundary conditions are imposed on the artificial boundaries. On the other hand, the coupled nonlinear subproblem is an ODE system, which can be solved exactly. To demonstrate the effectiveness of our method, some comparisons in terms of accuracy and computational cost are made between the PML approach and our method in numerical examples.

© 2015 Wiley Periodicals, Inc. We introduce a software package for the analysis of biomolecular solvation. The package collects computer codes that implement numerical methods for a variational implicit-solvent model (VISM). The input of the package includes the atomic data of biomolecules under consideration and the macroscopic parameters such as solute-solvent surface tension, bulk solvent density and ionic concentrations, and the dielectric coefficients. The output includes estimated solvation free energies and optimal macroscopic solute-solvent interfaces that are obtained by minimizing the VISM solvation free-energy functional among all possible solute-solvent interfaces enclosing the solute atoms. We review the VISM with various descriptions of electrostatics. We also review our numerical methods that consist mainly of the level-set method for relaxing the VISM free-energy functional and a compact coupling interface method for the dielectric Poisson-Boltzmann equation. Such numerical methods and algorithms constitute the central modules of the software package. We detail the structure of the package, format of input and output files, workflow of the codes, and the postprocessing of output data. Our demo application to a host-guest system illustrates how to use the package to perform solvation analysis for biomolecules, including ligand-receptor binding systems. The package is simple and flexible with respect to minimum adjustable parameters and a wide range of applications. Future extensions of the package use can include the efficient identification of ligand binding pockets on protein surfaces.

Near a charged surface, counterions of different valences and sizes cluster; and their concentration profiles stratify. At a distance from such a surface larger than the Debye length, the electric field is screened by counterions. Recent studies by a variational mean-field approach that includes ionic size effects and by Monte Carlo simulations both suggest that the counterion stratification is determined by the ionic valence-to-volume ratios. Central in the mean-field approach is a free-energy functional of ionic concentrations in which the ionic size effects are included through the entropic effect of solvent molecules. The corresponding equilibrium conditions define the generalized Boltzmann distributions relating the ionic concentrations to the electrostatic potential. This paper presents a detailed analysis and numerical calculations of such a free-energy functional to understand the dependence of the ionic charge density on the electrostatic potential through the generalized Boltzmann distributions, the role of ionic valence-to-volume ratios in the counterion stratification, and the modification of Debye length due to the effect of ionic sizes.

We incorporate the Poisson–Boltzmann (PB) theory of electrostatics into our variational implicit-solvent model (VISM) for the solvation of charged molecules in an aqueous solvent. In order to numerically relax the VISM free-energy functional by our level-set method, we develop highly accurate methods for solving the dielectric PB equation and for computing the dielectric boundary force. We also apply our VISM-PB theory to analyze the solvent potentials of mean force and the effect of charges on the hydrophobic hydration for some selected molecular systems. These include some single ions, two charged particles, two charged plates, and the host–guest system Cucurbit[7]uril and Bicyclo[2.2.2]octane. Our computational results show that VISM with PB theory can capture well the sensitive response of capillary evaporation to the charge in hydrophobic confinement and the polymodal hydration behavior and can provide accurate estimates of binding affinity of the host–guest system. We finally discuss several issues for further improvement of VISM.

The interplay between geometry and electrostatics contributes significantly to hydrophobic interactions of biomolecules in an aqueous solution. With an implicit solvent, such a system can be described macroscopically by the dielectric boundary that separates the high-dielectric solvent from low-dielectric solutes. This work concerns the motion of a model cylindrical dielectric boundary as the steepest descent of a free-energy functional that consists of both the surface and electrostatic energies. The effective dielectric boundary force is defined and an explicit formula of the force is obtained. It is found that such a force always points from the solvent region to solute region. In the case that the interior of a cylinder is of a lower dielectric, the motion of the dielectric boundary is initially driven dominantly by the surface force but is then driven inward quickly to the cylindrical axis by both the surface and electrostatic forces. In the case that the interior of a cylinder is of a higher dielectric, the competition between the geometrical and electrostatic contributions leads to the existence of equilibrium boundaries that are circular cylinders. Linear stability analysis is presented to show that such an equilibrium is only stable for a perturbation with a wavenumber larger than a critical value. Numerical simulations are reported for both of the cases, confirming the analysis on the role of each component of the driving force. Implications of the mathematical findings to the understanding of charged molecular systems are discussed.

Competitive adsorption of counterions of multiple species to charged surfaces is studied by a size-effect-included mean-field theory and Monte Carlo (MC) simulations. The mean-field electrostatic free-energy functional of ionic concentrations, constrained by Poisson's equation, is numerically minimized by an augmented Lagrangian multiplier method. Unrestricted primitive models and canonical ensemble MC simulations with the Metropolis criterion are used to predict the ionic distributions around a charged surface. It is found that, for a low surface charge density, the adsorption of ions with a higher valence is preferable, agreeing with existing studies. For a highly charged surface, both the mean-field theory and the MC simulations demonstrate that the counterions bind tightly around the charged surface, resulting in a stratification of counterions of different species. The competition between mixed entropy and electrostatic energetics leads to a compromise that the ionic species with a higher valence-to-volume ratio has a larger probability to form the first layer of stratification. In particular, the MC simulations confirm the crucial role of ionic valence-to-volume ratios in the competitive adsorption to charged surfaces that had been previously predicted by the mean-field theory. The charge inversion for ionic systems with salt is predicted by the MC simulations but not by the mean-field theory. This work provides a better understanding of competitive adsorption of counterions to charged surfaces and calls for further studies on the ionic size effect with application to large-scale biomolecular modeling.