key: cord-0827006-5t7t8oba
authors: Muralidharan, Ajay; Yethiraj, Arun
title: Fast estimation of ion-pairing for screening electrolytes: A cluster can approximate a bulk liquid
date: 2021-10-26
journal: J Chem Phys
DOI: 10.1063/5.0077013
sha: b092e92af4ce8265af2e9db2a004ed78b6ee4ed6
doc_id: 827006
cord_uid: 5t7t8oba

The propensity for ion-pairing can often dictate the thermodynamic and kinetic properties of electrolyte solutions. Fast and accurate estimates of ion-pairing can thus be extremely valuable for supplementing design and screening efforts for novel electrolytes. Here, we introduce an efficient cluster model to estimate the local ion-pair potential-of-mean-force (PMF) between ionic solutes in electrolytes. The model incorporates an ion-pair and a few layers of explicit solvent in a gas-phase cluster and leverages an enhanced sampling approach to achieve high efficiency and accuracy. We employ harmonic restraints to prevent solvent escape from the cluster and restrict sampling of large inter ion distances. We develop a Cluster Ion-Pair Sampling (CLIPS) tool that implements our cluster model and demonstrate its potential utility for screening simple and poly-electrolyte systems.

Ion-ion correlations and pair potentials-of-mean-force (PMF) are fundamental elements that enter the theory for ion-pairing in electrolyte solutions. [1] [2] [3] The propensity for ion-pairing can often dictate the conductivity and viscosity of electrolyte solutions. 4 In battery electrolytes, for example, weak correlations between lithium ion (Li + ) and the anion is often a prerequisite for achieving both a high Li + conductivity and a high transference number. 5 This is because long-lived charge-neutral aggregates do not contribute to useful transport in batteries. In addition, high energetic barriers for dissociation of ions can also lead to lowered salt solubility leading to poor performance. Hence, quick and accurate estimates of ionpairing and dissociation behavior between ionic species can be an immensely useful criteria for assessing electrolytes.

Extensive theoretical and computational efforts have been directed towards estimating ion-ion PMF's in a variety of contexts. [6] [7] [8] [9] [10] [11] [12] [13] [14] However, obtaining accurate and converged PMF's with computational models that include bulk explicit solvent can be computationally intensive. Because of the difficulty of these molecular calculations, a cluster model for estimating ion pair PMF's that reliably captures the important local features would be extremely valuable. Here, we utilize a simplified cluster model and perform free energy calculations along a chosen inter-molecular distance coordinate (r) to estimate dissociation barriers and short range PMF behavior between ion-pairs. (Figure 1 ) Li + with arbitrary anion and solvent (green triangles) is pictured above. We perform enhanced sampling along the distance coordinate, r, between the ions. A one sided harmonic upper wall restricts sampling of r above 0.65 nm and a solvent barrier prevents solvent diffusion away from the ions.

for the best compromise between efficiency and accuracy. The first choice is the level of chemical detail to be included for reasonable estimates of the PMF. For the case of ionic solutes, including explicit solvent has been shown to be essential for capturing specific ion effects [15] [16] [17] and other local effects. 6 This is because the physical process of dissociation of a contact ion pair (CP) to form a solvent shared/separated ion pair (SSIP/SIP) also involves the migration of solvent from outer shells into the inner solvation shell of the ions. Here, we incorporate a few layers of solvent (depending on solvent identity and density) to solvate the ions and fill the space between the ions up to a predetermined distance set to 0.65 nm. Beyond 0.65 nm the ion-pair PMF behavior is relatively flat for simple univalent ions (vide infra). Additionally, we apply a solvent barrier to prevent solvent diffusion away from the ions. We achieve this using a harmonic restraint on arXiv:2110.14036v1 [cond-mat.soft] 26 Oct 2021 the number of solvents within a spherical radius, R sol , relative to the center of distance between the ions (Figure 1 ). When R sol is empirically set to 2 nm, we show that it sufficiently captures the local PMF behavior for small univalent ions relevant to battery electrolytes. The second choice is the level of theory required to capture the physics of ion dissociation. Ab initio molecular dynamics based on modern wave function or density functional theory approaches are accurate but computationally intensive. Classical non polarizable force-fields based on partial charges derived from quantum mechanics are more efficient and adequate for screening purposes. For instance, previous estimates of the PMF for NaCl dissociation show that force-field models capture the physics of binding accurately compared to AIMD calculations. 18 The method developed here is applicable for any choice of (classical or ab intio) interaction model and is expected to be computationally feasible in all cases.

The final choice is the method used for molecular dynamics sampling. Simulations based on plain vanilla molecular dynamics are useful when the barriers associated with association are small. However, enhanced sampling methods are required for evaluating PMF's in a general situation with larger barriers. Here, we choose a state-of-the-art, robust and flexible approach that incorporates an on-the-fly probability enhanced sampling (OPES) scheme. 19 The OPES sampling scheme is able to sample all configurations necessary for characterizing the binding behavior, starting from a contact pair to a solvent separated pair in a back and forth manner until convergence of free energy. To further accelerate convergence, we restrict sampling of non-local distances (≥0.65nm) by utilizing a one sided harmonic upper wall. These practical choices lead to an efficient computational model for estimating local ion-pair PMF's. In the last few decades, modern computational approaches have played a significant role in engineering materials with desirable properties. Recent efforts such as the Materials Initiative 20,21 seek to establish an infrastructure to accelerate advanced materials discovery by leveraging the use of computational and experimental capabilities in an integrated approach to materials engineering. In the sphere of electrolyte design, the Electrolyte Genome project 22 enabled an automated high throughput computation of properties such as oxidationreduction potential, salt solubility, and electrochemical solution stability for assessing and screening novel battery electrolytes. Further, screening electrolytes for lithium ion (Li + ) battery applications requires the simultaneous optimization of a variety of design parameters. Thus, a need of the hour is the development of robust, efficient models and tools to compute design parameters for supplementing such design and screening efforts.

To this end, we create a tool titled "CLIPS" (short for CLuster Ion-Pair Sampling) that estimates metrics related to the binding of ionic species using our simplified cluster model ( Figure 1 ). This tool is made available for free on Github (https://github.com/ajaymur91/ CLIPS.git). We achieve a remarkable efficiency (≈ 5 minutes per PMF calculation with modest resources: 50% of cpu resources on a single AMD Ryzen 5 3600) by combining a reduced gas-phase cluster model (containing an ion pair and a few layers of solvent) with an enhanced sampling scheme that allows efficient crossing of free energy barriers. We demonstrate how this tool can generate useful metrics related to ion correlations for screening candidates for battery electrolytes as well as poly-electrolyte solutions.

On-the-fly Probability Enhanced Sampling Method (OPES) 19 is a state-of-the-art method to sample systems that present with large energy barriers along slow reaction coordinates or collective variables (CV's). We use the PLUMED 23-25 plugin with Gromacs 2019.6 for OPES simulations. We perform 5 ns of OPES at 313K in an equilibrated box of N s solvent molecules with 1 ion pair. We apply a history dependent bias, V n (r), at each step to the distance (r) between the cation (e.g. Li + ) and a distinguished atom of the anion (e.g. S of TFSI).

where P n (r) is an estimate (at step n) of the unbiased probability distribution of r, and and Z n are the regularization and normalization terms over the explored CV-space. 19 We set the maximum barrier height (∆E) or the free energy barrier to overcome to ∆E = 100 kJ/mol where = e −β∆E/(1−1/γ) and γ = β∆E. This choice of ∆E is based on sampling observed in trial OPES simulations. We also limit the exploration of the cation-anion distance (r) by introducing a large upper wall bias according to:

where κ w = 200000 kJ/nm and r 0 = 0.65 nm. We apply a solvent barrier using harmonic restraints relative to the center of distance (c) between the ions to prevent solvent diffusion away from the ions (Figure 1) .

where κ s = 1000 kJ/nm. where N is the number of solvent molecules within a spherical radius R sol . We choose a radius, R sol = 2 nm, for all cases in our analysis. and N s = 30 or 60 for the case of ethylene carbonate and water respectively. In general, a solvent molecule is considered to lie inside a spherical volume if a distinguished atom (i) of the solvent lies within the sphere. However, in practice we use a smooth switching function (s ic ) that goes from 0 to 1 based on the distance of the distinguished solvent atom i from c.

We make standard choices of n = 6, m = 2n for rational switching functions. R sol = 2 nm. Finally, we obtain an estimate for the unbiased distribution for the cation-anion distance, r, by using a reweighting scheme 26 :

Here, 0 represents an ensemble average over an unbiased trajectory, while V T represents a biased trajectory. V T is the sum total of the deposited bias after time t, including bias due to OPES, upper wall and solvent barrier (eq. 1,2,3). From the unbiased distribution, r 0 , we obtain the unbiased free energy as − log P (r). Further, we add a Jacobian correction 27,28 of 2kT ln r to the distance dependent free energies estimated from the unbiased probability distribution (eq. 6).

We use Gromacs 2019.6 [29] [30] [31] [32] [33] [34] [35] with PLUMED 23-25 plugin to perform OPES simulations of NaCl (in H 2 O) and some common Li + battery electrolytes. Lithium bis(trifluoromethane) sulfonimide (LiTFSI) and Lithium trifluoromethanesulfonate (LiOTf) in ethylene carbonate (EC) are ideal test cases because they display distinct ion pairing behavior despite the anions being structurally similar. 6 Force field parameters for Li + , OTf − , TFSI − and ethylene carbonate are taken from our previous work. 6 OPLSAA parameters 36 are used for Na + and Cl − with SPC water model. 37 All bonds are constrained using the LINCS algorithm. Electrostatic interactions are calculated using the particle mesh ewald technique with a 1.2 nm cutoff. A cut-off Lennard-Jones interaction is implemented for the non-electrostatic interactions with a 1.2 nm cutoff. Simulation systems undergo energy minimization using the steepest descent minimization algorithm. After energy minimization, 5 ns of OPES is performed with the stochastic dynamics 38 integrator at 313K in a simulation box with a few layers of solvent and 1 ion pair. Our implementation is made available on Zenodo 39 and Github.

We perform a 1µs NVT simulation (Figure 2 After energy minimization using a steepest descent algorithm, we carry out an initial 1µs of NVT equilibration with the Nose-Hoover thermostat 42,43 at a density of 1320 kg/m 3 followed by a 1µs molecular dynamics trajectory that is used for further analysis. We perform cluster simulations to characterize the PMF between Li + and bis(nonenylmalonato)-borate (BMB − ) anion in EC solvent. We label the partial charge distribution from the model of Wang et al. 40 as Q1 (inset, Figure 6 ). We also create an alternate hypothetical delocalized charge model for BMB − anion labeled as Q2. Note that we obtain Q2 model (from Q1) by manual redistribution of atomic partial charges (inset Figure 6 ).

Our cluster model with enhanced sampling (Figure 1 ) approximates bulk solvent ion-pairing potential of mean force (PMF). The cluster simulations begin with the ions in the contact pair (CP) state. Dissociation barriers are crossed as OPES biases (eq. 1) are added, and the solvent shared ion pair (SSIP) and solvent separated ion pair (SIP) states are explored revealing the local free energy landscape. By fast sampling back and forth, we obtain very good estimates of the dissociation and association barriers and free energy differences between the CP and SIP states. An animation of our model in action is provided on figshare (https://doi.org/10.6084/m9. figshare.16755277.v3). The upper wall (blue vertical lines in the figures) harmonic bias prevents the exploration of distances beyond 0.65 nm and is sufficient to capture the interesting details of the local free energy landscape for small ions.

When compared to bulk solvent simulations, the cluster model shows good agreement for the PMF for NaCl salt in water ( Figure 3 ) and LiTFSI and Li-OTf in EC (Figure 4) . In all cases, the agreement is excellent for distances less than the upper wall at 0.65nm. Beyond that distance the bulk solvent free energy is relatively flat compared to smaller distances and hence not targeted by our cluster model. For NaCl, the barrier between the CIP and SSIP is ∼ 4.55kT which is accurately reproduced by the cluster model. For Li + -OTf − that barrier is much higher ∼ 9.21kT , however, reproduced faithfully by the cluster model as well. The model is therefore able to capture the essential physics of ion-pairing for distinct cases of free energy barriers at a very modest computational expense.

In light of the successful treatment of PMF's for simple salts with our cluster model, we shift our focus to the case of poly-electrolytes. Recent experiments 44 have sought to develop a polymer solution electrolyte based on poly(bis(nonenylmalonato)-borate) [Li(pBMB)] that is electrochemically stable over a wide potential window in organic carbonate solvents with high Li + transference for batteries. However, in practice the poly-electrolyte [Li(pBMB)] appears to make a gel in EC solvent the mechanism for which is not clear. Figure 5 depicts the PMF between Li + and BMB − , and shows that there is a deep attractive well which can result in strong ion pairing. This calculation was based on the force field obtained from Wang et al. 40 which we label as Q1. From molecular dynamics simulations of the solution, we find that in this case there is an aggregation of ions that persists over the duration of the simulation. Figure 6 (left panel) shows the final configuration from molecular dynamics simulations at experimentally relevant 44 concentration of 0.3M Li(pBMB) in EC solvent. The degree of polymerization of the poly-anion is 5. The simulations reveal that the system based on the Q1 (strong interaction) model for BMB − formed non-homogeneous micro porous structures with self diffusivities of 0.24 x 10 −6 cm 2 /s for Li + and 0.05 x 10 −6 cm 2 /s for monomer beads of BMB − .

The undesirable gelation behavior of the polyelectrolyte can be reduced dramatically with subtle changes in the charge distribution of the anion. To test the effect of charge delocalization we adjust the partial charges in the Q1 force-field by hand. In the resulting Q2 model, the partial charges on O b and O c oxygen atoms are reduced by 8% and 15% (inset) respectively compared to the Q1 model. The charges on boron and carbonyl carbon are adjusted to keep the net charge of Q2 model to be -1. With these slight adjustments the PMF well depth decreases by more than a factor of two (see Figure 5 ) in cluster simulations. This difference has a dramatic influence on the micro-structure, and the weak association (for Q2 model) between Li + and BMB − in the solution phase leading to homogeneous micro-structure with better Li + transport (see Figure 6 ). The Q2 based polyelectrolyte model produces improved self diffusivities of 5.22 x 10 −6 cm 2 /s for Li + and 0.57 x 10 −6 cm 2 /s for monomer beads of BMB − .

Ion-pairing can determine the viscosity and ionic transport coefficients in electrolyte solutions. An intrinsic propensity for ion pairing is the PMF between ions in dilute solution. We introduce a robust and efficient cluster model to estimate this PMF. To achieve high efficiency, the model combines enhanced sampling OPES with a cluster of explicit solvent molecules. The predictions of the model are in excellent agreement with bulk solvent simulations.

We use the model to predict the PMF between Li + and the complex anion BMB − in EC. The calculations show a deep attractive well that should lead to ion-pairing in dilute solution. From molecular dynamics simulations of the bulk system we find that there is an aggregation of ion pairs and a gel-like micro-structure. We find that subtle changes in the charge distribution of the anion (without changing the overall charge of course) results in a large difference in the PMF, and the micro-structure of the bulk liquid, which becomes homogeneous. The ion self-diffusion coefficients increase by an order of magni-tude. The calculations therefore provide a design principle for anion selection. Judicious modification of the anion chemistry to make the charge more diffuse might result in better structural and transport properties for the poly-electrolyte.

In concentrated solutions, especially with divalent ions, aggregation of ions might play an important role. In this case, in addition to the pair PMF the free energy of multiple ions could be of importance. The cluster model can be extended to such systems with larger clusters of solvent molecules, and an appropriate choice of collective variables, i.e., order parameters. Investigating many body effects is a possible future direction of this work.

An interesting aspect to consider is the the number of solvent molecules required to reproduce the PMF obtained from bulk solvent simulations 6 . In all our calculations, we empirically set the number of solvent molecules to completely cover Li + and the space between the ions at all ionic separations up to 0.65 nm. It might he possible to optimize this because in practice we are not interested in the long-range behavior of the PMF, but rather the short-ranged behavior where ion-pairing occurs. Another possibility is to use very few solvent molecules, and a reaction field.

The method is easily extensible to more realistic interactions between the ions and solvent molecules. Our CLIPS tool performs calculations at the classical force-field level because highly automated implementations of accurate force-fields are available (OPLSAA Ligpargen, 36,45 CHARMM 46 ) with parameters derived from quantum mechanics. However, our model is general and independent of the level of theory used for performing the simulations. For instance, if force-fields are 

Ion pair potentialsof-mean-force in water

Generalized transition state theory in terms of the potential of mean force

Marcus theory of ion-pairing

Conductivity and viscosity of solutions of licf3 so 3, li (cf 3 so 2) 2 n, and their mixtures

Ion correlations and their impact on transport in polymer-based electrolytes

Why lithium ions stick to some anions and not others

Water-mediated ion pairing: Occurrence and relevance

Potential of mean force between oppositely charged nanoparticles: A comprehensive comparison between poisson-boltzmann theory and monte carlo simulations

Driving force for the complexation of charged polypeptides

Ion pairing in molecular simulations of aqueous alkali halide solutions

Simulation study of ion pairing in concentrated aqueous salt solutions with a polarizable force field

Calculation of effective interaction potentials from radial distribution functions: A reverse monte carlo approach

Molecular dynamics-potential of mean force calculations as a tool for understanding ion permeation and selectivity in narrow channels

Potential of mean force between two molecular ions in a polar molecular solvent: A study by the three-dimensional reference interaction site model

Quasi-chemical theory for anion hydration and specific ion effects: Cl − (aq) vs. f − (aq)

Quasi-chemical theory with cluster sampling from ab initio molecular dynamics: Fluoride (f-) anion hydration

Hydration mimicry by membrane ion channels

Dissociation of nacl in water from ab initio molecular dynamics simulations

Rethinking metadynamics: from bias potentials to probability distributions

US), Materials genome initiative for global competitiveness (Executive Office of the President

Commentary: The materials project: A materials genome approach to accelerating materials innovation

The electrolyte genome project: A big data approach in battery materials discovery

Plumed: A portable plugin for free-energy calculations with molecular dynamics

Plumed 2: New feathers for an old bird

Promoting transparency and reproducibility in enhanced molecular simulations

A time-independent free energy estimator for metadynamics

Free energy calculations

A comparison of methods to compute the potential of mean force

GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers

GRO-MACS: A message-passing parallel molecular dynamics implementation

GRO-MACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation

GROMACS 2019.4 Source code

Tackling Exascale Software Challenges in Molecular Dynamics Simulations with GROMACS

GROMACS 4.5: A high-throughput and highly parallel open source molecular simulation toolkit

GROMACS: Fast, flexible, and free

Molecular Force Field for Ionic Liquids Composed of Triflate or Bistriflylimide Anions

Interaction models for water in relation to protein hydration

Efficient algorithms for langevin and dpd dynamics

Cluster ion-pair sampling (clips) tool for fast estimation of ion-pairing

Atomistic insight into orthoborate-based ionic liquids: force field development and evaluation

Optimization of the opls-aa force field for long hydrocarbons

A molecular dynamics method for simulations in the canonical ensemble

Canonical dynamics: Equilibrium phase-space distributions

Electrochemically Stable, High Transference Number Lithium Bis(malonato)borate Polymer Solution Electrolytes

Ligpargen web server: an automatic opls-aa parameter generator for organic ligands

Charmm general force field: A force field for drug-like molecules compatible with the charmm all-atom additive biological force fields

Appendix A: Table of contents graphic

This research was supported in part by UW-Madison Department of Chemistry PHOENIX & STARLING research cluster through National Science Foundation Grant CHE-0840494. This study was also partially supported by US Department of Energy, Basic Energy Sciences Contract DE-SC0017877. A.M. is a Hirschfelder fellow.