Application of Screening Functions as Cutoff-Based Alternatives to

Subscriber access provided by READING UNIV

Article

Application of Screening Functions as Cut-off Based Alternative to Ewald Summation in Molecular Dynamics Simulations Using Polarizable Force Fields Jenel Vatamanu, Oleg Borodin, and Dmitry Bedrov J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b01043 • Publication Date (Web): 02 Jan 2018 Downloaded from http://pubs.acs.org on January 3, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Application of Screening Functions as Cut-off Based Alternative to Ewald Summation in Molecular Dynamics Simulations Using Polarizable Force Fields Jenel Vatamanu1,2,*, Oleg Borodin2, Dmitry Bedrov1 *

[email protected]

1

Department of Materials Science and Engineering, University of Utah, 122 South Central Campus Dr., Salt Lake

City, Utah, 84112, USA 2

Electrochemistry Branch, Sensors and Electron Devices Directorate, Army Research Laboratory, 2800 Powder

Mill Rd., Adelphi, MD, 20703, USA

Abstract The range-dependent screening of the charge – charge, charge – induced dipole and induced dipole – induced dipole interactions was examined for a variety of liquids modeled using polarizable force fields. A cut-off based method for calculation of the electrostatic interactions in molecular dynamics (MD) is presented as an alternative to Ewald-type summation for simulations of the disordered materials modeled using many-body polarizable force fields with permanent charges and induced point dipoles. The proposed approach was tested on bulk water, room temperature ionic liquids, and solutions of ions in polar solvents. The smooth, short range and atom type independent screening functions for interactions between the charges and induced dipoles were obtained using the force matching approach. An excellent agreement for both the magnitude and directionality of forces, structural and dynamic properties was found in MD simulations utilizing the developed screening functions compared to those with Ewald summation. While similar in shape and range the charge-charge screening functions were somewhat dependent on the material chemistry. The charge - induced dipole and induced dipoleinduced dipole screening functions, on the other hand, were found to be nearly universal for the tested materials.

1 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table of Content Figure (TOC):


Page 2 of 34

Page 3 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction Molecular dynamics (MD) simulations are commonly used for prediction of the material structural and dynamic properties as well as obtaining molecular scale insight into mechanisms of ion and solvent transport. In classical MD simulations, the interatomic forces can be categorized as short-range van der Waals interactions (typically diminishing within 10-12 Å) and the long-range electrostatic interactions originated from interactions between partial atomic charges and (in polarizable force fields) induced dipoles. The computationally most demanding part of the MD simulations is calculation of the long-range forces. To emulate the bulk-like macroscopic sample using a smaller simulation cell accessible to classical MD simulations, the original simulation cell (with its geometry defined by = , , ,

= a , a , a is replicated in each direction, generating the so called ‘supercell’. Taking the periodic images of the primary cell into account, the electrostatic interaction energy V of a system with N-point charges qi is given by

= ∑k* ∑,

(1)

,

where index k labels the summation over the periodic images, and k=nL with n being a vector of integers. The star symbol in the summation over k indicates that the self-interaction is excluded in the primary cell, i.e., i≠j when k=0. The rij,k represents the distance between the charge i in the primary cell and the charge j in the kth periodic image of the cell. The factor ½ in Eq. 1 removes the double counting of each pair of charges. The computation of the long-range interactions involving the summation over the infinite number of periodic cells (Eq. 1) is typically done using the Ewald summation technique.1 The Ewald summation approach can be conceptually understood as follows: at the position of each point charge, a Gaussian charge equal and opposite in sign to that of the point-charge, is added. The interactions between the added Gaussian-distributed screening functions is evaluated in the Fourier (or reciprocal) space while interactions of the point charges screened with the Gaussian functions become short-ranged and, hence, are evaluated in real space. In the Ewald summation, the total potential energy of the system consisting of a collection of point charges qi that satisfies the electroneutrality condition, ∑ ! = 0, becomes1, 2: = ∑)*)+*),-, ,; 0

∑+*),-, ,; 0

#$%&'$ +

1 C 34 56 7 /9:7 ∑DE ;∑ 2 67

#$%'$ + GH, I#JK#L$M


! exp −@AB ; −

:

√1

(2)

∑ ! −


According to Eq. 2, the slowly convergent sum Σk*Σijqiqj/rij,k is divided into two components: i) a real space component space component 34 56 7 /9:7 , 67

#$%&'$ computable within a short-range cut-off radius and ii) a reciprocal

34 56 7 /9: 7 ;∑ 67

! exp −@AB ;

which, due to the Gaussian envelope

converges in the Fourier space. Note, that the second term in Eq. 2 contains the self- and

intra-molecular interactions that are subtracted by the third and fourth terms of Eq. 2. In this equation, α is the so-called Ewald splitting parameter and it is related to the width of the assumed Gaussian charges. 3, 4

The last term GH, I#JK#L$M depends on the geometry (or periodicity) of the summation and on the

dipole moment of the simulation cell.2, 3, 5, 6 It can be shown that for 3D systems and spherical symmetry of the summation the last term J becomes zero if the metallic (or tinfoil) boundaries are assumed. For details on the GH, I#JK#L$M term see previous works.7-11 As commonly utilized for bulk 3D-periodic systems, metallic (infinite dielectric constant of surrounding medium) boundary conditions were assumed in this work for the periodic directions. The main computational burden of the Ewald method is the reciprocal term. Optimized implementations of the Ewald scheme can reduce the overall computational cost to O(N3/2). 12 However, this scaling becomes expensive for large systems (i.e., systems with more than 20,000 atoms). The approximation of the reciprocal part with the smooth particle mesh Ewald (SPME)4 method, on the other hand, can significantly decrease the computational time to O(NlogN). In SPME techniques, an evenly spaced grid of charges is fit to the actual charge distributions in the system and the reciprocal part of the energy is evaluated based on the generated charge grid in the reciprocal space. The electrostatic field at atoms’ locations is back interpolated from the regular grid transformed in real space by an inverse fast Fourier transform. The reciprocal part contribution to atomic forces is often computed from the derivatives of the spline functions utilized to back interpolate the fields. The speed-up of the SPME arises from the utilization of the fast NlogN Fourier algorithm, possible on evenly spaced grids. However, the SPME can still become computationally prohibitive for the larger scale simulations mainly because it is difficult to efficiently parallelize the 3D fast Fourier transform. For example, there are phenomena, such as grain formation, grain domain evolution, crack formation and propagation, polymorphism during crystallization, dislocation dynamics, etc. for which simulations with millions or even billions might be necessary.13,

14

Also, the sampling of rare events such as homogeneous nucleation of crystals,

deposition/electrodeposition, and dynamics in glassy systems often requires very long trajectories (multiple microseconds). For such large-scale simulations further improvements even for highly optimized SPME approach are desired to reduce computational cost.


Page 4 of 34

Page 5 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The utilization of cut-off based (COB) alternatives to Ewald-type summations, which approximate the energies/forces only within a spherical cut-off and eliminate the reciprocal space calculations, is attractive due to lower computational cost and straightforward parallelization. A number of COB schemes has been previously proposed: i) truncation/shift of the energies/forces based on either atom or neutral group cut-off15-19 (extended to dipoles20), ii) reaction-field methods21-23, iii) the Wolf summation and its variants24-26, iv) isotropic periodic sum27 for 3D, 2D and 1D periodicity, v) iterative Boltzmann inversion28-30, vi) local molecular31 field theory or symmetry-preserving32 mean-field theories,33-35 and vii) force-matching short-range pairwise forces36, 37 or pair-type independent screening functions based on Izvekov and Voth36 scheme to coarse-grain in “interaction space”. exception of the mean field approaches designed for inhomogeneous systems (e.g. ref.

32

36, 38

With the

and ref 35) most

COB approaches are primarily applicable for bulk systems where structural heterogeneities are on the scale of cut-off radius used in the simulation. For the condensed phase ionic systems, where ions are present at rather high concentrations, there are physical grounds for considering the shorter-range effective interactions. Because of the electrostatic screening by the local environment (other neighboring ions, surrounding polar solvents, etc.), the total/effective forces/energies in the condensed phase decay faster vs. distance than the corresponding electrostatic interactions in vacuum. Therefore, choosing the effective short-range forces/energies that resemble the actual physical screening in the system is conceptually possible. While there are instances of the electrostatic interactions screened out to truly short-ranged terms,39 in general, additional techniques are required to impose the decay of electrostatic forces to zero at the arbitrary predefined desired cut-off distance. The force matching approach is well-suited for finding such short range effective interactions. This method requires a set of atomic forces be available from the high-fidelity simulations (i.e., with full Ewald summation or ab-initio MD). Short-ranged forces are then fit to the atomic forces calculated using a high-fidelity model. The screened short-range forces are typically represented numerically and depend on the chemical environment of the interacting species. Two methods are commonly used for the force matching. The first approach is to fit a type :N

dependent numerical function described by coefficients %6

for all non-bonded αβ atom pair types :N

defined by the site/atom of types α and β. The grid spacing $6 over which the coefficients %6 ($6 ) are calculated encompasses interatomic separations ranging from close contacts to the predefined cut-off distance. The fit minimizes the following objective function: Q,R,S,TUVW,

O = ∑Q,R,S ∑ P

XYZ

− ∑

:N

I%6 , $

,


(3)


Page 6 of 34

:N

where the g(f,r) is an interpolation function (spline, polynomial or linear) of forces from %6

at distances

r between the discrete grid points. In Eq. 3, the superscript x,y,z labels the Cartesian component of the XYZ

force vector Fi. The quantity ∑

:N

I%6 , $ is the total force on atom i evaluated using the fitted

:N

short-ranged force coefficients %6 . In the case of linear interpolation the function g is: :N

:N

5_

:N

I [%6 , $ \ = %6 ]6 + %6^ ]6^ where ]6^ =

_`a 5_

:N

In this approach the number of coefficients %6

, ]6 = 1 − ]6^

(4)

increases as a square of the number of atom types present

in the system. For example, a system containing atoms of 10 different types will require 55 distinct 'c

types of pairwise combinations. Using 100 grid points to represent each % :N

coefficients %6

:N

will result in 55,000

necessary to describe all interaction in the system. For that many variables, a fitting

procedure that involves direct matrix inversion algorithm or conjugate gradient methods quickly becomes computationally inefficient, not to mention the complexity of the logistics for documenting and disseminating of a “force field” with such a large number of coefficients that might work only for one specific system. The second approach to force matching electrostatic interactions is to fit the atom pair type independent screening functions sk as was originally proposed by Izvekov et al.

36, 38

for molecular

simulations employing the non-polarizable force fields containing only permanent atomic charges. In this case, the objective function for minimization is: O

Q,R,S,TUVW,

= ∑Q,R,S ∑ P

XYZ

− ! ∑

! Id e6 , $

(5)

The coefficients sk are type independent because the pre-factors qiqj are not included in the fit (compare with Eqs. 3 and 5). As the screening function Id is the same for all pair types the number of coefficients to be fitted is equal to the number of points representing the grid rk, i.e. about a hundred. In this work, we extend the approach by Izvekov and Voth to MD simulations using polarizable force fields. Polarizable force fields gain popularity due to increases in computational power, improved transferability, ease of parameterization using the gas-phase cluster energies obtained from quantum chemistry, and often lead to more accurate description of ionic structure and transport for electrolytes containing small or divalent ions. For example, the structure and dynamics of electrolytes doped with Li+ salts was accurately predicted in simulations using polarizable force fields.40-42 Three types of polarizable models are commonly using in classical MD simulations: i) charge equilibration,43, 44 ii) particle on the 6 ACS Paragon Plus Environment

Page 7 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


spring or Drude oscillators,45 and iii) induced atomic dipoles.40-42,

46-48

Representation of atomic

polarizabilities with induced dipoles has proven to be rather robust for very diverse set of systems. For a more detailed introduction and discussion on this topic see e.g., ref. [49] and references therein. Here we propose and benchmark the technique to extend the screening function COB method to the polarizable force fields with induced dipoles. The remainder of the paper is organized as follows. The mathematical formalism for the computation of the screening functions for charge-induced dipole and induced dipoleinduced dipole interactions is presented in the next section. Next, the simulation and fitting protocol will be described. The transferability of the short-range charge-dipole and dipole-dipole screening function coefficients is discussed for bulk water, room temperature ionic liquids, and solutions of ions in polar solvents. Finally, a benchmark of the performance of MD simulations employing screening functions is given.

2. Methodology In MD simulations with polarizable force fields, the total electrostatic energy U of a molecular system containing point charges and point induced dipoles can be written as,

U = U C − U thole +

where

αi

1 µiµi ∑ − E14 2 i αi

(6)

is the polarizability assigned for atom i and

µi = (µix , µiy , µiz ) the corresponding induced

dipoles, the UC is the electrostatic energy due to the interactions between charges and dipoles in the systems, the

Uthole

is the Thole50,

51

dumping over short distances of the dipole-dipole interactions to

avoid the phenomenon of “polarization catastrophe” than 4α α /g to each other. The

52

possible if the point induced dipoles are closer

1 µ iµi ∑ is the polarization energy and the E14 is a damping 2 i αi

energy term of interactions between atoms within the same molecule separated by three consecutive chemical bonds. In this work we employed the APPLE&P force field42, 53 in which only the charge-dipole 1-4 interactions were damped by 20% for RTILs and Li-salts in carbonate solvents and by 50% for solutions of dimethoxymethane (DME) in water. The energy Uc is given by:



UC =

where

Page 8 of 34

1 1 Qi Q j ∑ 2 i, j ,k* R ij + kL

(7)

Qi is a generalized charge formally defined as, Qi = qi + µi o ∇i

In Eq. 7 the qi is the electrostatic scalar charge and µ i is the induced dipole on atom i. The term Uc can be expediently computed with the Ewald summation modified for mixture of point charges and point dipoles as described e.g., in ref. [54]. The induced dipoles are calculated iteratively by minimizing the total energy given by Eq. 6, with respect to the Cartesian components of the dipoles,

hi hjk,l,m

= 0.

Next we show how to rearrange the Uc in a form that makes apparent the definition of type pair independent short-range screening functions. As the charge-charge, charge-dipole and dipole-dipole energy/forces decay differently with distance, the Uc is decomposed into three components as follows,

U c = U q −q + U q − µ + U µ − µ where

(8)

Uq−q is the charge-charge energy given by U q −q =

qi q j 1 , ∑ 2 i , j ,k * r ijk

(9)

U q −µ is the charge-dipole energy defined as U q−µ =

and

(

)

1 1 qi (µ jrijk ) − q j (µ i rijk ) 3 , ∑ 2 i , j ,k * rijk

(10)

U µ −µ is the dipole-dipole energy,

U µ −µ =

 µ iµ j 1 d  + [(µ i rijk )(µ jrijk )] ∑ 3 2 i , j ,k *  rijk drijk

 1   r3  ijk

 1 µ µ    = ∑  i j − 3 [(µ i rijk )(µ jrijk )]    2 i , j ,k *  r 3  rijk5   ijk 

(11)

Taking into account that the forces on atom i are defined from the electrostatic energy U as Fi=-∇ ∇iU, it is straightforward to show that Eqs. 9-11 lead to the following expressions for the charge-charge, chargedipole, and dipole-dipole components of forces acting on atom i:


Page 9 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fq −q (ri ) = − ∑

qi q j rijk

(12)

rijk2 rijk

j ,k *

(

Fq−µ (ri ) = −∑ qi (µ jrijk ) − q j (µ i rijk ) j ,k*

Fµ − µ (ri ) = − ∑ (µ i µ j ) j ,k *

d  1 drijk  rijk3

 1  3  ijk  rijk

 rijk  1  − ∑ (qi µ j − q j µ i ) 3 r r  ijk j ,k*  ijk

) drd

 rijk d  + [(µ i rijk )(µ jrijk )] r drijk  ijk

 3   r5  ijk

   

(13)

 rijk  3  − [µ i (µ jrijk ) + µ j (µ i rijk )] 5 r r  ijk  ijk

   (14) 

Based on the structure of Eqs. 12-14, it is possible to define type independent screening functions s(r), as follows:

Fq −q (ri ) = −∑ qi q j j

rij rij

s qq ( rij )

(15)

(

)(

Fq − µ (ri ) = −∑ qi (µ jrij ) − q j (µ i rij ) s qµ (rij ) j

) r − ∑ (q µ r '

ij

i

ij

Fµ − µ (ri ) = −∑ (µ i µ j )s µµ (rij ) j

rij rij

− q j µ i )s qµ (rij )

(16)

j

− [(µ i rij )(µ jrij )] s µµ (rij )

qµ

qq

(

j

)r '

ij

rij

− [µ i (µ jrij ) + µ j (µ i rij )]s µµ (rij ) (17)

µµ

In Eqs. 15-17, s (rij ) , s (rij ) and s (rij ) are the so-called charge-charge, charge-dipole and dipoledipole screening functions respectively. The corresponding expressions for the energies are obtained from forces Eq. 15-17 as: c

U q −q

1 = ∑ qi q j ∫ s qq (rij ) dr 2 i, j r

Uq−µ =

U µ −µ

(

(18)

)

1 qi (µ jrij ) − q j (µ irij ) s qµ (rij ) ∑ 2 i, j

(19)

c  1  = ∑  µ i µ j ∫ s µµ ( rij )dr + [(µ i rij )(µ jrij )]s µµ ( rij )  (20) 2 i, j  r 



Page 10 of 34

To implement this method in a computer code one has to supply (either numerically on a grid or as a c

qq

function) the following functions: i) s (r ) and

∫s

qq

qµ ( r )dr for the charge-charge interactions, ii) s (r )

r

and its derivative

ds qµ (r )

dr

µµ for the charge-dipole interactions, and iii) the s (r ) , together with its integral

c

∫s

µµ

( r )dr and its derivative

r

s µµ and

ds µµ

dr

ds µµ ( r )

dr

qq

for the dipole-dipole interactions. Note that the s , s c

qq

µµ

,

ds qµ

dr

,

c

∫

∫s

r

r

are fitted from force-matching (Eqs. 3-5) while the functions s qq ( r ) dr and

can be straightforwardly calculated by integrating the fitted s and s

qµ

µµ

( r )dr

.

It is important to note that the summation over periodic cells k in Eqs. 9-14 is no longer needed because the screening functions s are chosen to be short-range and decay to zero at the cut-off distance. Utilizing the short-range screening functions s provides (in Eqs 15-20) a smooth truncation of forces to zero and a shift of energies to zero at the cut-off, in a way that atomic forces calculated with this COB scheme match the forces of the high fidelity model such as Ewald summation on average. Based on the above definitions for the screening functions s, it can be shown that in the limit of small interatomic separation where the electrostatic interactions can be assumed as not screened, the sqq, sqµ and sµµ become qµ

s (r → 0) = rs (rij → 0) = − qq

r 2 s µµ (rij → 0)

3

=

1 r2

(21)

It is straightforward to verify that substituting Eq. 21 in Eqs. 15-20 generates the unscreened electrostatics in vacuum for large r (Eqs. 9-14) and therefore the screening functions s reproduce the correct asymptotic behavior in the low charge density limit. In order to fit the sqq from atomic forces, we used the method given in ref. [38] for linear αβ interpolation functions I [%6 , $ \ as in Eq. 4. The charge-dipole and dipole-dipole terms s

qµ

and s

µµ

were fitted utilizing two functions, s and its derivative s’=ds/dr under the constraints that both s and s’ became zero at cut-off. From the fit, a set of coefficients sk on the grid with distance spacing rk, is obtained. Note that from the obtained sets of points (sk, rk) polynomial dependences s=s(r), also called screening polynomials, can be further fitted. 10 ACS Paragon Plus Environment

Page 11 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


As an example, shown below are expressions of the charge-charge, charge-dipole, and dipoledipole screening polynomials s(r) fitted for the RPOL polarizable water model and a 10Å cut-off, are given in supporting information. Note that the interactions at distances within close contacts were not sampled well (if at all) during simulations due to their rare occurrence, and therefore these polynomials should be modified at short distances r within close-contacts, for example to the corresponding unscreened form Eq. 21. 36, 38 In simulations utilizing such screening function approach the computational speed-up arises from i) elimination of the reciprocal part of the Ewald summation, ii) force evaluation within a short-range cutoff, iii) easier and efficient parallelization in the codes that employ spatial decomposition and hence limiting computations only to the short-range forces, and iv) faster convergence of the self-consistent iterations that computes the induced dipoles.

Figure 1. (a) The screening functions sqq, sqµ, sµµ computed from the RPOL water model. The dashed line shows the corresponding unscreened terms (1/r2) from bare electrostatics. (b) The dependence of sqq on distance for various systems. The units of the y–axis are 1/Å2.

3. Simulation and Screening Function Fitting Protocols The simulations were performed utilizing the in-house developed MD code that incorporates electronic polarization via point induced dipoles. The simulations were conducted in cubic cells with periodic boundary condition in all three directions. The bulk densities of the studied systems were 11 ACS Paragon Plus Environment


determined from preliminary NPT trajectories performed at atmospheric pressure and corresponding predefined temperatures. The production runs were done in the NVT ensemble, at the density predicted by NPT simulations. The equations of motion were integrated with the reversible multiple time-scale integrator scheme.

56

55

and outer time step of 3fs. The temperature was controlled utilizing the Nosé-Hoover

The details of simulated systems are summarized in Table 1 and S1 from the supporting

information. A three-step procedure was utilized as follows: i) First, we generated several nanoseconds (6-20 depending on the system and temperature) of simulation trajectories utilizing full 3D Ewald method. At this stage, the reciprocal part of Ewald was handled with the more expedient SPME technique.57 The configurations were saved every 300fs. From these configurations we recomputed the forces, energies and the induced dipoles with higher accuracy, i.e, using classical Ewald summation with 15-17Å (depending on cell dimensions) real space cut-off and 12x12x12 k-space grid for the reciprocal space part. The non-SPM version of Ewald used for the higher precision force/energy/dipoles updates ensured that the atomic forces utilized in the force matching fitting strictly obeyed the Newton third law. ii) The screening coefficients sqq(rk), sqµ(rk), sµµ(rk) were fit over a distance grid rk with k=1..P, using the least square approach and an objective to match the atomic forces from previous step, as described in the previous sections and in refs. [36, 38, 58] From the fitted discrete numerical dependence s=s(rk), k=1..P , polynomial functions s=s(r) were further fitted. Note, that when implementing Eqs. S1(SI)-S3(SI) (given in supporting information) in computer code all the shown decimals must be explicitly written, i.e., they are double precision numbers. iii) Finally, MD simulations were re-run utilizing the COB scheme with the obtained short-range screening functions s(r) instead of the Ewald summation. The accuracy of simulations employing COB approach was assessed by comparing the results obtained from simulations utilizing the Ewald method for electrostatic calculations.


Page 12 of 34

Page 13 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


System

Force field

Temperature (K) 298 298 1250 298

Thole parameter (Å) 0.5 0.5 0.5 0.5

298

0.4

H2O/DME

Ref. 59 Ref.60 Ref. 61 Ref. 61 for NaCl Ref. 60 for water Ref. 62

C2mim-TFSI Pyr13-TFSI 1M LiPF6/EC:DMC(3:7)

Ref. 53 Ref.53 Ref.53, 63

393 393 363

0.2 0.2 0.2

1M LiPF6/TMS:DMC(1:2)

Ref.53, 64

363

0.2

LiTFSI/DMMSA

Ref.53

363

0.2

H2O DG H2O RPOL NaCl melt H2O/NaCl(10/1)

Number of molecules 1434 1434 1000 pairs 1290 water 72 NaCl 600 water 24 DME 70 pairs 70 pairs 114 EC 256 DMC 31 LiPF6 120 TMS 240 DMC 36 LiPF6 238 DMSA 59 LiTFSI

Table 1. A summary of the investigated systems.

To test the proposed COB approach we considered the following systems: bulk water, solution of NaCl and dimethoxymethane (DME) in water, room temperature ionic liquids of 1-ethyl-3methylimidazolium bis(trifluoromethylsulfonyl)imide bis(trifluoromethylsulfonyl)imide

(c2mimTFSI), N-propyl-N-methylpyrrolidinium

(pyr13TFSI), and 1-n-butyl-3-methylimidazolium tetracyanoborate

(c4mimB(CN)4), melt of NaCl, and electrolytes containing lithium hexafluorophosphate (LiPF6) or lithium bis(trifluoromethylsulfonyl)imide (LiTFSI) salts in polar solvents such as tetramethylene sulfone (TMS),

ethylene

carbonate

(EC)/dimethyl

carbonate

(DMC)

mixture,

and

N,N-

dimethylmethanesulfonamide (DMMSA) solvents (see Table S1). Ideally, one would like to have the screening functions to be system independent. However, as we will show in the next section, there is slight system dependence. As the charge-charge interactions have the longer-range scaling with the distance and therefore will have the most impact, the sqq were refitted for each system. However, because the charge-dipole and dipole-dipole interactions are shorter-range than the charge-charge interactions, the sqµ and sµµ were computed only for liquid water and their transferability was tested for other systems. Such an approach to calculate the sqq for all systems while probing the transferability sqµ and sµµ is a reasonable trade-off between the accuracy and the difficulty to extract these coefficients.



4. Quality of the Screening Function Fitting In this section, we discuss several measures that estimate the quality of the fits of the proposed screening functions based on the force matching approach as well as the dependence of these screening functions on such parameters as the imposed cut-off radius, system composition, temperature, etc. In the next section, we present results of MD simulations employing these screening functions and compare structural and dynamic property predictions with those obtained from high fidelity simulations (i.e., the full Ewald summation).

Figure 2. The histograms of (a) Uqq, (b) Uqq,COB - Uqq ,Ewald, (c) Uqµ, and (d) Uqµ,COB - Uqµ,Ewald obtained using screening functions with 10, 12, 14, and 16Å cut-off distance for molten NaCl at 1250K. The histograms were computed from the energies assigned to individual atoms. In these histograms the unit of energy is 1 kJ/mol for individual atoms. The histograms are normalized to 1. The pairwise energies were assigned equally on both individual atoms i, j constituting the pair ij. The sum of all individual atom energies gives twice the total energy, which is in agreement with the ½ prefactor in the right-hand-side of Eq. 1.

Figure 1a shows the fitted sqq(rk), rksqµ(rk), and -⅓rk2sµµ(rk) functions obtained for the case of bulk water using the three-site RPOL polarizable water model.60 Supporting information shows additional results for the 4-site polarizable (DC) water model given in ref. 59. For comparison, the unscreened scaling of 1/r2 is also shown. As expected, at short distances, the screening functions are similar to the unscreened terms and follow a linear dependence with the slope of -2 on the log-log scale. On the other hand, with 14 ACS Paragon Plus Environment

Page 14 of 34

Page 15 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


approaching the cut-off separations, the screening functions decay much faster than the unscreened 1/r2 factor (Fig. 1a) such that at cut-off (in this case Rc=10 Å) they become zero. Note, that if the Rc is changed to a different value the sharp decay of the screening functions will shift to the cut-off radius, essentially smoothly truncating all interactions outside of the Rc to zero. Figure 1b shows the influence of the chemical details of each system on the sqq obtained using the force matching approach. differences in s

qq

The

become noticeable at distances larger than 6Å. In other words, the specificity of the

system has a more pronounced influence on sqq at distances closer to Rc and it determines the shape of the sqq decay to zero. One the other hand, the charge-dipole sqµ and dipole-dipole sµµ screening coefficients for different systems were found to be very similar therefore suggesting their transferability from one system to another, which we have utilized in this work. Next we assess the accuracy of the fitting protocol and the ability of the obtained screened interactions to represent characteristics of the high fidelity electrostatic interactions by comparing several properties calculated with the screening functions COB and with the Ewald method for the configurations (atomic positions) used in the force matching fitting. The induced dipoles were recomputed by iterations for each utilized method. The compared properties (labeled with “X”) are the energy (U), the force (F), and the induced dipoles (µ µ), for the charge-charge, charge-dipole, and dipole-dipole components. For the vector quantities F and µ, the statistics were analyzed on the norm ||X||=(Xx2+Xy2+Xz2)1/2 and on the appended set of all three x,y,z components, i.e., the union set Fx∪Fy∪Fz. We verified that each Cartesian component Fx, Fy, and Fz generates the same statistics as the Fx∪Fy∪Fz, which is expected for isotropic systems. Three types of the quantitative measures of the COB accuracy were studied for the mentioned above properties. The first measure is based on the slope, intercept and regression coefficient of the linear fit of XCOB =f(XEwald) dependence. Note that a slope equal to one and an intercept equal to zero correspond to ideal fit. The parameters of the linear fit XCOB =f(XEwald) as well as the sum of deviations of absolute values of atomic forces, energies, and induced dipoles are shown in the Tables S2-S4 shown in Supporting Information. Table S2 compares these properties for different systems, Table S3 shows the role of the cut off distance used for the screening coefficients, and Table S4 compares different methods as well as the transferability of the screening function from one water model to another. The second measure of the screening function performance is based on comparing the histograms (distributions) computed with COB and with Ewald approaches for X=[U,F,µ µ,] as well as the differences XEwald – XCOB. These histograms are shown in Figures 2 and 3 for NaCl melt (force-field token from ref. 61) and bulk water, respectively, and in the SI for additional systems. Finally, the third measure is based on the preservation of the direction of the forces and induced dipoles. For this purpose histograms of the 15 ACS Paragon Plus Environment


µEwald,µ µCOB) are shown in Fig. 4. Note that the statistics for properties “X” and cos(FEwald,FCOB) and cos(µ their distributions shown in Tables S2-S4 (see SI) and in Figures 2-4 and Figures S2-S3 (shown in SI) were collected from thousands of configurations each containing several thousand atoms.

Figure 3. The histograms of (a) Uqq, (b) Uqµ, (c) Uµµ, (d) Fqq, (e) Fqµ, (f) Fµµ, (g) |Fqq|, (h) |Fqµ|, (i) |Fµµ|, (j) µ, (k) |µ µ| and (l) induced energies, calculated with screening functions over 10 Å cut-off (red lines) and Ewald (black lines), for the RPOL water system. The blue lines represent the histograms of the differences between the quantities calculated with the screening function and with Ewald. The histograms for the energies and forces were computed from the corresponding values assigned on individual atoms. The numbers in y-axis are probabilities and the number in x-axes represent energies (panels a-c and l), forces (panels d-i) and induced dipoles (panels j-k). Note that in panels (d) and (g) for Fqq and |F|qq the red and black lines overlap in the shown scale and only black line is visible. The unit of energy is 1kJ/mol per individual atom. The unit of force is 1 kJ/mol/Å per individual atom. The unit of dipole is 1Debye. These histograms are normalized to integrate to 1.

Comparison of COB with Ewald predictions. We begin with analyzing the accuracy of the description of forces by fitted screening functions. Shown in Tables S2 are the slopes, intercepts and regression coefficient obtained from the linear fit of FCOB =f(FEwald) as well as the deviation of Σ|F| from Σ|FEwald| for several representative systems. For water the regression slope for the charge-charge forces Fqq is 0.9996, while for the charge-dipole and dipole-dipole forces the slope drops to 0.986 and 0.971,


Page 16 of 34

Page 17 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


respectively. In other words, the faster the decay vs. distance of the fitted interactions, the larger the error of the forces recomputed with the short-ranged screening functions. Note that the slopes for dipole-dipole forces are similar to the slopes fitted from dipole-dipole energies (Tables S2). For the screening functions computed with Rc=10 Å, the fitted slope of Fqq varies as a function of system from 0.990 to 0.999 as follows: 0.9900 for NaCl melt, 0.9944 for pyr13TFSI, 0.9980 for 1M NaCl, 0.9988 for DME in water (force field token from ref. 62), 0.9970 for Li-salts in carbonate solvents, and 0.9990 for bulk water using the DG59 model. (Table S2). Table S3 shows the regression characteristics obtained for bulk water represented by the RPOL water and various Rc distances utilized for the screening functions. Increasing the Rc used for sqq increased the Fqq slopes from 0.9966 (Rc=5Å) to 0.9996 (Rc=12Å) and to 0.9993 (for Rc=16Å). Note that increasing the cut off above 12 Å does not improve the Fqq slope, however, it slightly improves the regression coefficients. The regression coefficient for Fqq increased systematically from 0.9965 (Rc=5Å) to 0.9997 (Rc=10Å) and to 0.999945 (Rc=16Å) (see Table S4). The procedure to fit the sqq coefficients converged with less than 100 (independent) configurations for all studied systems. The faster decaying (vs. distance) sqµ(r) and sµµ(r) can be more difficult to converge with an unbiased fit. A biased fit, giving more weight to coefficients near Rc, generated short-ranged sqµ(r) and sµµ(r) for all studied systems, however the predicted slopes, intercepts and regression coefficients for FCOB=f(FEwald) did not improve compared to the results obtained with transferred sqµ(r) and sµµ(r) from those obtained for water models. Based on that, the following approach was pursued in this work: the sqq screening functions for longer range and higher magnitude forces were refitted for each system, while the shorter-range sqµ and sµµ were obtained with the unbiased quickly converging fit for water and then transferred to other systems. As we will show below, the screening functions sqµ and sµµ from water fits are largely transferable to other systems. For example, for the RPOL water model with its own sqµ the slope of Fqµ, COB=f(Fqµ,Ewald) is 0.986. Some systems with sqµ transferred from RPOL water generated comparable or even better slopes than the fits for the original RPOL water with its own fitted sqµ. For example, the slopes of the Fqµ, COB=f(Fqµ,Ewald) linear fits using water sqµ are: 1.006 for NaCl melt, 0.990 for NaCl in water , 1.063 for DME in water, 0.970 for pyr13TFSI, and 0.988 for LiTFSI/DMMSA64 system. Similarly, the even faster decaying dipole-dipole terms were also obtained from fits for water system and then transferred to other systems. The transferability of the sµµ screening function from water model to other systems can be clearly seen from Table S2. A lower accuracy description for energies (as compared to forces) is expected from the screening functions approach because the fit is performed on atomic forces and not on energies. Indeed, as shown in



Page 18 of 34

Tables S2 and S3, the intercepts from the linear fits of Uqq,COB =f(Uqq,Ewald) are deviating from zero and the combined energy terms Σ|Uqq,COB| show larger systematic error than the force terms Σ|Fqq,COB|. The differences between the histograms of Uqq,COB and Uqq,Ewald are illustrated in Figure 2 for NaCl melt at 1250K. Specifically, the Uqq, COB histograms computed with COB approach are systematically shifted by as much as 50% from the Uqq,Ewald histograms computed with Ewald. Increasing the Rc utilized to fit the screening functions from 10 to 16 Å moves the Uqq, COB histograms closer to the Uqq,Ewald distributions. However, even for the largest cut-off utilized (Rc=16 Å) a large gap between the two histograms still remains (Figure 2a). Other ionic liquid systems presented in SI show similar trends for Uqq,

COB

histograms. In contrast to ionic systems discussed above, the liquids comprised of neutral polar molecules yield the significantly more accurate Uqq. This is demonstrated by the values of slopes/intercepts/Σ|Uqq| in Table S2 for water. For a 10Å cut-off used to fit the screening functions, the slope of Uqq, COB =f(Uqq,Ewald) is 0.957 compared to 0.55-0.7 for NaCl melt or RTILs. For polar molecules, a relatively small cut-off for screening functions fit is sufficient to generate averaged Uqq,COB within 1% from the Ewald values. As shown in Table S3 for the RPOL water, the increase of the Rc from 10 to 16 Å improved the Σ|Uqq| from being 5% off from the Ewald prediction to only 0.6% deviation. The Uqq,COB histograms for the studied polar systems are in a significantly better agreement/overlap with the Uqq,Ewald histograms than it is the case for ionic liquids (compare Figures 3a and 2a). Importantly, the histograms of Uqq computed using COB approach for polar systems are not shifted from the Ewald Uqq to the large extent as observed for ionic liquids (Figure 3a). While the histograms of the energy differences Uqq,COB – Uqq-Ewald for NaCl is broad and shifted from zero even at the largest cut-off (Fig. 2d), the same histograms for water are centered very close to zero (compare Figs. 2d and 3a). In summary, the above analysis indicates that for the charge-charge interactions, the screening functions will work reasonably well for systems comprised of polar molecules, but will result in systematic shift (truncation) of predicted energies for ionic liquids. Such behavior for Uqq underlines the limitation of this approach for ionic systems, i.e., it gives poor prediction of the energy and thermodynamic properties related to the energy, e.g. heat of vaporization. Nevertheless, the forces acting on atoms are described very well and, therefore, one can expect the ionic liquid structure and dynamics to be captured well using the COB approach as discussed in the next section. Next, we examine the charge-dipole and dipole-dipole energy components which are important because the induced dipoles are the derivative of the energy with respect to components µx, µy, µz. Note that, for any given atomic configuration, large errors for Uqq shown above for ionic systems will have no impact on the accuracy of the induced dipoles since Uqq does not depend on µ. As shown in Figs. 2 and 3, 18 ACS Paragon Plus Environment

Page 19 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the histograms of Uqµ and Uµµ computed from COB approach (Eqs. 18-20) are in a significantly better agreement with those computed using Ewald. Also, the histograms of differences Uqµ,COB - Uqµ,Ewald and Uµµ,COB - Eµµ,Ewald are centered very close to zero both for ionic liquids (Figs. 2and S2) and for neutral polar molecules (Figs 3b-c), indicating that the average values of Uqµ and Uµµ computed using COB approach agree well with Ewald values. Tables S2-S4 show that Σ| Uqµ | and Σ| Uµµ | are represented more accurately using COB approach

than Σ|Uqq| contributions. The histograms of the induced energies n),op-, = ∑

qr qr :

computed with

COB and Ewald also are in good agreement, i.e., they almost overlap (Fig. 3l). As expected from trends observed in Tables S2-S4 and Figs. 3-4 for Uqµ and Uµµ, the histograms of the induced dipoles calculated with screening function (µCOB) also agree remarkably well with those computed with Ewald (µEwald), see Figs. 3j-k. Furthermore, the histograms of the differences µCOB-µ µEwald, shown as dashed blue lines, are almost as sharp as the histograms for forces (Figures 3d-i) for which the actual fit (force matching) was performed. The slope of the fitted linear dependence of µCOB=f(µ µEwald) varies from 0.988 to 1.005 indicating high accuracy of representing induced dipoles within COB approach because they are determined by electric field at each atom that is closely related to the forces used to fit the screening functions. Comparison with other COB methods. In order to assess the accuracy of the screening functions obtained using the force matching approach against other COB approximations, we calculated in Table S4 the fitted slopes, the intercepts, and the regression coefficients of XCOB =f(XEwald) for the RPOL water using the following COB approximations: i) its own RPOL screening functions, ii) screening functions transferred from DG water model, iii) a simple truncation of energies and forces given by Eqs. 9-14, and iv) a shift of energies and a smooth truncation of energies and forces to zero at Rc as following:

n

= ∑,

n

j

= ∑, s

n

jj

= ∑, s

[1 − Yp \

qt Art u

qr qt u

+

−

(22)

qr Art u

qr Art qt Art z

y

v s1 − 3 [Yp \ + 2 [Yp \ v

y

9

v s1 − 4 [ \ + 3 [ \ v Yp

Yp

(23)

(24)

Note that the tapering polynomials applied to energies in Eqs 22-24 will also smoothly tapper the resulted forces F=-∇ ∇U to zero at cut-off. Two cut-off distances were considered: Rc=10 and Rc=16Å. As expected, the precision increases with the increase of the cut-off for all methods. In the case of a simple truncation, 19 ACS Paragon Plus Environment


Page 20 of 34

Fqq,COB= f(Fqq,Ewald) are spread out with somewhat small regression coefficients of 0.9674 for 10 Å and 0.987 for 16 Å cut-off. In contrast, the COB schemes studied here generated regression coefficients for forces Fqq larger than 0.9995. For energies Uqq,COB= f(Uqq,Ewald), simple truncation generated regression coefficients as small as 0.06, as expected. The screening functions or the smooth truncation of forces by Eq. 22-24 dramatically improved the Σ|Uqq| within a few percentages from the Ewald values (see Table S4); also, the regression coefficients improved from 0.06 (for simple cut off) to above 0.97 (for screening functions). It is remarkable that the COB scheme (Eq. 22-24) provides such excellent description of the slopes, intercepts, regression coefficients and Σ|Uqq| compared to Ewald values, indicating that carefully tuned COB approach can be a viable alternative to the Ewald summation for simulation of bulk systems. 15

Table S4 shows that the screening functions presented in this work clearly outperform the COB scheme of Eqs. 22-24 for all three charge-charge, charge-dipole and dipole-dipoles components. The predicted slopes of XCOB= f(XEwald) drop as follows: for Uqq from 0.959 (with screening) to 0.918 (with Eqs. 22), for Uqµ from 0.945 (with screening) to 0.784 (with Eq. 23), and for Uµµ from 0.965 (with screening) to 0.763 (with Eqs. 24). Similar improvements are observed for forces, i.e. the slopes decrease as follows, for Fqq from 0.9992 (with screening) to 0.9779 (the force from Eq. 22), for Fqµ from 0.986 (with screening) to 0.87 (the force from Eq. 23), and for Fµµ from 0.989 (with screening) to 0.886 (the force from Eq. 24). The regression coefficients and the Σ|Xqq,

qµ, µµ

| in Table S4 also show better results

for the screening functions then for the COB scheme by Eqs. 22-24. Directionality of the forces and induced dipoles. It is important to accurately capture induced dipole directionality in addition to its magnitude.65 Deviations of the induced dipole and total force vector directions obtained from the COB vs. Ewald are examined by considering cos(∢(FEwald,FCOB)) and cos(∢(µ µEwald,µ µCOB)) distributions. The resulting histograms were then integrated from -1 (i.e., antiparallel orientations of FEwald and F COB) to a given value x (of the cos(∢(FEwald,F COB)) ), Q

|} = ~5[Histogram of &Je [∢ , \] [&Je∢ , \]

(25)

The cumulative integral of histograms I(x) was normalized to 1. The values of I(x) at a given value of x ( i.e., cos(∢(FEwald,FCOB))) represent the percentage of force vectors computed with the screening functions that are ‘misaligned’ from those computed using Ewald by more than the particular value of x=cos(∢(FEwald,FCOB)). For the random distribution of angles between FEwald and FCOB , the I(x) increases linearly from 0 to 100% as x increases from -1 to 1. For identical vectors, I(x) is zero in the interval x∈[1,1) and 100% at x=1. Figure 4 shows that the I(x) for the RPOL water has very small values for most of 20 ACS Paragon Plus Environment

Page 21 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the x interval and then sharply increases at x>0.99, essentially indicating a very large extent of co-aligned orientations of Ewald and COB computed forces and induced dipoles. For any given angle, the force directionality is best preserved for induced dipoles µ and Fqq and slightly less accurately captured for Fqµ and Fµµ forces. For example, let us consider the RPOL water and select a value of x=0.99838 which corresponds to an angle of 3.2 degrees between FEwald and FCOB vectors. The percentage of vectors of FCOB that deviates from collinearity from FEwald with more than 3.2 degrees (quantified by I(0.99838) in Figure 5) is less than 0.15 % for µ, 1% for Fqq, and 4-5% for Fqµ and Fµµ . Note that for a random distribution of vectors (angles), the percentage of angles with larger deviation than this x would be 99.838%. Therefore the screening functions method is quite successful in preserving the force and induced dipole directionality for the RPOL water. The collinearity of FEwald and FCOB increases systematically with the increase of cut-off radius used in COB fitting (see Figures 4 and 5). An important question is: to what extent other COB methods preserve the directionality of the forces, i.e., how much improvement the proposed screening functions really bring compared to other truncation approaches? To address this question in Figure 5a for the RPOL water we show the deviation from collinearity (vs. Ewald) for the forces and induced dipoles calculated with i) simple truncation, ii) tapered energies Eqs. 22-24, iii) screening functions transferred from DG-water model, and iv) RPOL model fitted screening functions. As expected, the simple truncation approach generates very large I(0.99838), i.e., almost no preservation of force directionality. The smooth truncation approach of Eqs. 22-24 decreases the deviation (for x=0.99838) to 2% which is a dramatic improvement, while the force matching fitted screening functions decrease the deviation to as low as 1%. However, the shifted force scheme by Eqs. 22-24 did poorly in preserving the orientation of the dipoles, i.e., around 15% of dipoles were deviating from those predicted by Ewald by more than 3.2 degrees. In contrast, the screening functions resulted in only 0.25% induced dipole vectors that deviated more than 3.2 degrees from the Ewald predictions. Therefore, as compared to other simpler COB schemes, the screening functions achieve a significant improvement to the directionality of the induced diploes (see Figure 5 and S3). Furthermore, the preservation of directionality for Fqµ and Fµµ is increased substantially if the COB based on Eqs. 22-24 is replaced with the screening functions. Therefore the screening functions capture the directionality of all key vector quantities in polarizable systems better than the alternative approaches.



Figure 4. The cumulative distributions I(x) (in percent, %) of the deviations of vector directions using COS approach from Ewald for the force components (a-c) and induced dipoles (d), where x=cos(∢(FEwald,FCOB)) for the charge-charge (a), charge dipole (b) and dipole-dipole (c) force components, and x=cos(∢(µ µEwald,µ µCOB)) (d) for the induced dipoles. Functions correspond to different cut-offs for the RPOL water model.

Figure 5. The dependence of I(0.99838) as a function of (a) cut-off radius and (b) applied COB method. See Figure S3 from SI for more details regarding the panel (b).


Page 22 of 34

Page 23 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5. Molecular Dynamics Simulations Using the Screening Function COB Approach Accurate prediction of the liquid structural and dynamic properties is central to many MD simulations. Ability of the COB approach to accurately approximate the structure and dynamics was examined in MD simulations performed in NVT ensemble using densities obtained from Ewald simulations. The velocity autocorrelation functions were calculated in the NVE ensemble. Structural properties. The structure of bulk water is quantified by the nearest neighbors distribution y

and defined by the so called tetrahedral order parameter qt : ! = 1 − ∑y ∑96^&Je6 + 1/3 , where the summation over i, k indices uniquely cycles through all triangles formed by the four nearest neighbors water molecules i, k around a centrally placed water molecule j. The parameter qt was calculated as follows: for each oxygen atom of water the closest four other oxygen were identified and for such configuration of five oxygen atoms the qt was computed. The qt is normalized such that for a perfectly tetrahedral network with cosθijk=-1/3 the qt has a value of one, while for a random distribution (which cancels ) the ensemble averaged is zero.

Figure 6. A comparison of structural properties predicted from simulations with Ewald summation and COB approach with 10Å cut-off screening functions for the RPOL water model: (a) tetrahedral order parameter, (b) the distributions of induced permanent induced permanent cos(∢(µ µ ,µ µ )) as a function of cos(∢(µ ,µ )), and (c) rdfs of atomic pairs. The histograms in panels (a) and µ µ (b) are normalized to integrate to 1. For a random distribution the histogram in (b) will be a constant line.

Shown in Figure 6a are the histograms of collected for individual atoms from simulations of bulk water using the RPOL water model and computed with COB Rc=10 Å and Ewald approaches. Additional histograms of are presented in SI for the DG water model and aqueous solution of NaCl. Note that differences in qt distributions are expected to have a dramatic impact on liquid properties because the closest neighbors have the strongest interaction with the central molecule. As shown in Figure 23 ACS Paragon Plus Environment


6a the distributions of calculated with COB Rc=10 Å are in excellent agreement with the Ewald predictions. Also, we verified that decreasing the cut-off of the sqµ and sµµ to as low as 6Å, while maintaining the cut-off for sqq at 10Å, results in a similar accuracy as in Figure 6a indicating that the charge-dipole and dipole-dipole components can be treated as even shorter-range interactions than the typical cut-off radius for the van der Waals interactions. Similar conclusion is observed if we use other descriptors to characterize the reproducibility of water structure. Figure 6b and 6c compare the distribution of the angle between vectors defining the direction of induced and permanent dipoles on each water molecule and the regular radial distribution function, respectively. For both properties, the predictions by COB and Ewald are almost identical.

Figure 7. The 3D-space distribution function of water for (a) H and (b) O sites around a centrally considered water molecule (molecular frame of reference defined in the main text). The white wireframes and solid (a) red and (b) blue colored lobes represent areas in space of high probability for (a) H and (b) O relative to the water molecule. The overlap between the white wireframe and the colored red/blue lobes shows that COB methods predict the same space distribution of atoms as Ewald. Sdf for more complex systems (RTILs, Li-salts in solvent) are shown in Figures S5 and S6 from Supporting Information.

The arrangement of neighboring atoms can also be characterize with the 3D spatial distribution functions (sdf). In order to calculate sdf, a molecular frame was defined based on three vectors. For water, two of the vectors defined the HOH plane (one vector connected O with the middle of the HH vector, and the second vector connecting H atoms) and the third one was perpendicular to the HOH plane. Then the distribution of other atoms was binned depending on the atom’s location around centrally placed water molecule in the defined reference frame. As shown in Figure 7, the regions of higher density for the neighboring H and O atoms (from other molecules) lie preferentially along the O-H bond of the reference (central) molecule. A second lobe of H/O density is located on a plane perpendicular to the HOH of the central molecule on the backside of the reference molecule resulting in a tetrahedral arrangement of nearest neighboring O atoms. Figure 7 also shows a lobe of high density of the second coordination shell 24 ACS Paragon Plus Environment

Page 24 of 34

Page 25 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in the areas extended above and beneath the HOH central plane. As it can be seen from Figure 7, the 3D sdf calculated with COB Rc=10 Å agrees very well with the Ewald summation predictions. Supplemental information S5 and S6-B shows similar excellent agreement for sdf for more complex RTILs and Libased electrolytes. The radial distribution functions for molten NaCl obtained with screening functions utilizing cut off radii ranging between 10 Å and 16 Å are shown in Figure S4-B (SI). This case is interesting because the Uqq calculated with screening functions (see Figures 2a and 2b) were off from the Ewald values by as much as 50%, and the increase of cut off from 10 to 16Å did not improve the Uqq-COB (see the intercept/slopes in the Table S2). However, despite such large systematic errors in Uqq-COB , the ensemble averaged rdfs shown in Figure S4-B (SI) indicate no noticeable difference with the cut-off increase from 10 to 16 Å and they are almost identical to those from MD simulation using Ewald. In the case of NaCl melt (and other RTILs), an accurate fit of the longer range (and most dominant) term Fqq (sqq) was sufficient to generate virtually identical liquid structure as in simulations with Ewald. The shorter range screening functions sqµ and sµµ were sufficiently accurate to generate a similar extent of agreement in rdfs as for RTILs, NaCl and water/DME. In other words, sqµ and sµµ are transferable from one system to another without significant loss in accuracy.

Next, we examine a more complex system of the Li-salt/organic solvent electrolytes, specifically, LiPF6 dissolved in EC:DMC(3:7) mixture (force-field token from ref.

63

) and the LiPF6 in

TMS:DMC(1:2) mixture. In order to probe the screening function method for electrolytes, we want to understand how well simulations with the screening functions reproduce the solvation structure of Li+ ions. As Li+ is the key element in Li-ion-batteries, its solvation and coordination by the solvent are particularly important for the ion transport mechanisms and the formation of solid electrolyte interphases at electrode surfaces. Note that the structure of these electrolytes in bulk phase has been extensively studied using APPLE&P force field in previous works of Borodin et al.40-42, 53 In this paper we compare the radial distribution functions and the coordination numbers for the following pairs: Li-O(carbonyl), LiO(ether) and Li-F. Figure 8 confirms that the first solvation shell of Li+ is dominated by Li-F and LiO(carbonyl) large peaks in rdf at around 1.85 and 2Å, respectively, and very little by Li-O(ether) pairs, as expected. As shown in Figure 8, the coordination numbers (CN) of Li are 2 with O(DMC) at Li-O(DMC) distances between 2 to 4Å, 0.8-0.9 with O(EC) at distances Li-O(EC) between 2.3 to 4 Å and 1.1-1.2 with P(of PF6) at Li-P distances between 3.5 to 4 Å, or a total coordination number of Li of roughly 4. Figure 8 shows that rdf and CN of COB-10 (COB simulations with screening functions fitted to Rc=10 Å)



practically overlap with the Ewald values. Similar agreement between rdfs predicted from simulations using the COB and Ewald methods were found for most other systems as shown in Figure 9 and S4-S7 for other systems and atomic pairs.

Figure 8.The radial distribution functions and coordination numbers for LiPF6 salt in aprotic organic solvent utilized as electrolytes in Li-batteries (a,c) LiPF6 in EC:DMC and (b,d) LiPF6 in TMS:DMC, calculated with screening functions (red lines) and Ewald (black lines) methods.

The worst predictions of rdf were found for aqueous solutions of NaCl, specifically, for ion-ion correlations. As shown in Figure 9g-i, the location of the maxima/minima in rdf and the overall shape of rdf are qualitatively preserved for all pairs. However, spurious oscillations near Rc distance are observed particularly when an Rc =10 Å was used for screening functions. Increasing the screening to the cut-off above 14Å noticeably improves the rdf and removes the longer range spurious oscillations. Note that, depending on the utilized cut-off, poor predictions of rdf with COB methods are not surprising. For dilute ionic solutions the cut-off should be large enough to encompass sufficient ionic pairs to screen out the electrostatic field. In other words, the choice of cut-off radius should be based on the ionic screening length. For RTILs the screening length is small and the 10Å cut-off radius is sufficient. For solutions of ions, on the other hand, the screening length is large and it increases with the decrease of ion 26 ACS Paragon Plus Environment

Page 26 of 34

Page 27 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


concentration, therefore larger cut-offs are required. However, for systems consisting of Li-salts in aprotic solvents the screening functions with 10 Å cut-off radius generated a much better agreement with Ewald than the 1M solution of NaCl did (see Figures 9g-i).

Figure 9. The upper panels are the radial distribution functions in RTILs a) centers of mass in c2mimTFSI, (b) N-N groups in c2mimTFSI, (c) C-alkyl with N in pyr13TFSI. The middle row of panels show the radial distribution functions for ionic pairs of Li-salts in organic solvents relevant for Li-batteries: (d) LiPF6 in EC:DMC, (e) LiPF6 in TMS:DMC and (f) LiTFSI in DMMSA. The bottom panels show the radial distribution functions for ionic pairs of 1M aqueous NaCl as follows: (g) Cl-Cl, (h) Na-Cl and (i) Na-Na pairs. The panels (g)-(h) also illustrate the improvement of screening function predictions for solutions of ions if the radius utilized for screening functions is increased up to 16Å. Additional rdf for solvent-solvent and solvent-solute pairs shown in Figures S4-S7 Supporting Information show significantly better agreement COB vs Ewald rdf for pairs involving species in higher concentration.

We also verified the histograms of the angle between the permanent and induced dipoles cos(∢(µ µpermanent,µ µinduced)) and compared these histograms for trajectories calculated with Ewald and with screening functions. The µpermanent was calculated with respect to the center of mass of ions/molecules. The histograms of cos(∢(µ µpermanent,µ µinduced)) are shown in Figures S4-S6 in Supporting Information. There is a 27 ACS Paragon Plus Environment


Page 28 of 34

good preservation of the orientation of the induced dipole relative to the permanent dipole, i.e. Ewald vs COB-10 histograms almost overlap (see Figures S4-S6). Dynamical properties. Dynamics at both short and long timescale was probed by the velocity autocorrelation functions, mean square displacements (MSD) and diffusion coefficients (D). The diffusion coefficients were calculated based on the Einstein equation, i.e., from the slope of the linear fit of MSD in the limit of long times. Figure 10a shows the velocity autocorrelation functions for the subpicosecond short timescale regime. The COB approximation is quite accurate as the observed correlations are virtually identical to those predicted from simulations using Ewald. Figures S8 and S9 from SI show for multiple systems that the MSD and diffusion coefficients computed from simulations with screening functions are very similar to those obtained from the Ewald method. For example, the COB-10 simulations predict for water at 298K a diffusion coefficient of 2.56⋅10-9 m2/s, which is very close to 2.53⋅10-9 m2/s given by simulations with Ewald. For both Li+ and PF6- ions in EC:DMC solvent, the Ewald and COB-10 predict very similar diffusion coefficients of around 0.9⋅10-10 m2/s at 363K.

42

Somewhat larger discrepancies in the diffusion coefficients were observed for water/DME system at 298K, i.e., for the water solvent the COB-10 over-predicted the D to be 1.52⋅10-10 m2/s (the Ewald value is 1.45⋅10-10 m2/s). For c2mimTFSI RTIL at 363K the COB-10 generated dynamics as follows: the diffusion coefficient of c2mim ion was 1.98⋅10-10 m2/s (1.91⋅10-10 m2/s with Ewald) and for TFSI anion the diffusion coefficient was 1.09⋅10-10 m2/s (1.05⋅10-10 m2/s with Ewald), i.e. an excellent agreement. The influence of the cut-off distance on the diffusion coefficient was investigated for RPOL water using the following mixed cut-off range screening functions: the longer range charge-charge screening functions were calculated with a typical 10 Å cut-off radius while the shorter range sqµ and sµµ had cutoffs of 9, 8, 7, 6 and 5 Å. As shown in Figure S9-a, the MSDs were not affected by lowering the cut-off of sqµ and sµµ and the fitted diffusion coefficients only varied within 4%. As the screening functions for systems comprised of polar molecules (as well as RTILs) yielded good predictions (vs. Ewald) of the force/energy histograms and of rdfs, it is then not surprising that the dynamics of these compounds obtained from simulations using the COB approach is also accurate. It may be, however, somewhat surprising that the MSDs and D of aqueous NaCl mixture computed with COB-10 shown in Figure S9 (SI) agree rather well with the Ewald predictions despite that Na and Cl had larger disagreement and even spurious effects in rdfs for longer-range correlations shown in Figures 9g-i. For aqueous NaCl, increasing the screening function cut-off from 10 to 16 Å was required to eliminate the spurious effects in rdf; however this increase in the cut-off had a rather small impact on diffusion coefficients, i.e. over this range the COB diffusion coefficients change from 1.6⋅10-10 to 1.7⋅10-10m2/s for Na (D=1.6⋅10-10m2/s with


Page 29 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Ewald) and from 2.2⋅10-10 to 2.4⋅10-10m2/s for Cl (D=2.1⋅10-10 m2/s with Ewald). Finally, the coefficients of all studied systems were represented as DEwald=f(DCOB) are shown in Figure 10b. The good alignment of the points (DEwal,DCOB) on a log-log scale with a line of slope 1 shows great agreement between simulations using Ewald and COB methods.

Figure 10. (a) The velocity autocorrelation functions for the polarizable DG model of water as obtained from NVE simulations using Ewald (solid lines) and COB approach with screening functions (dashed lines). (b) The diffusion coefficients DCOB=f(DEwald) for multiple systems studied in this paper. The shown line is an ideal correlation with a slope of 1. Note that from the fit DCOB=f(DEwald) with the data represented in log-log scale, a slope of 1.0025 is indicating excellent agreement with only very slightly increasing dynamics for COB methods versus Ewald.

6. Conclusions Electrostatic interactions are screened in liquids via reorganization of induced dipoles, permanent dipoles and ion positions. The first term represents electronic contribution to dielectric constant measured at high frequencies, while the second and third terms yield the lower frequency or relaxed dielectric constant. They are critical for fundamental understanding of ion and molecular solvation as screening is largely responsible for salt dissociation in liquids and reduced interaction between ions in ionic liquids. In this contribution, the screening functions were investigated separately for the permanent charge – charge, charge – induced dipole and induced dipole – induced dipole contributions obtained from MD simulations with polarizable force fields. While the shape of the screening functions were qualitatively similar for all systems reflecting the increased screening as the charge separation increases to 10 – 15 Å, the screening functions exhibited some differences. Ability to accurately incorporate the long – range interactions into the short-range screening allowed us to suggest a cut-off based alternative to Ewald method. It was investigated for a variety of materials modeled using the many-body polarizable force fields with atombased induced dipoles.



The smooth truncation of forces to zero at the cut-off distance is achieved with short-range atom type-independent screening functions fitted to reproduce the atomic forces in Ewald-based simulation. The key point of the presented scheme is to split the interactions in the charge-charge, charge-dipoles and dipoles-dipoles components and to fit pair type- independent screening functions for each component. Only relatively short MD simulation trajectories were required in order to obtain the screening functions for the longer range charge-charge interactions. The shorter range charge-dipole and dipole-dipole screening functions are nearly universal and represented by polynomials, suggesting that they could be applied to MD simulations prediction of the structural and dynamics properties of other bulk materials with high fidelity. This is clearly an improvement compared to previous66 efforts to utilize distance dependent dielectric constant to effectively model charge – induced dipole interactions and screen them after the first peak in the radial distribution function.

Supplementary Information Shown in Supporting information are additional Tables and Figures showing Ewald versus COB histogram distributions of energy and forces, histograms of cos(∢(µ µpermanent,µ µinduced)), radial distribution functions, 3D spatial distribution functions.

Acknowledgment Authors gratefully acknowledge the support from project sponsored by the Army Research Laboratory under Cooperative Agreement Number W911NF-12-2-0023. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of ARL or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein. JV acknowledges funding support through ORAU (contract # W911NF-16-2-0107) to U.S. Army Research Laboratory.


Page 30 of 34

Page 31 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


References (1) Ewald, P. P. The Computation of Optical and Electrostatic Lattice Potentials. Ann. Physik 1921, 64, 253 (2) Simulation of Electrostatic Systems in Periodic Boundary Conditions. I. Lattice Sums and Dielectric Constants. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 1980, 373, 27-56. (3) Yeh, I.-C.; Berkowitz, M. L. Ewald Summation for Systems with Slab Geometry. J. Chem. Phys. 1999, 111, 3155-3162. (4) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577-8593. (5) Pan, C.; Hu, Z. Rigorous Error Bounds for Ewald Summation of Electrostatics at Planar Interfaces. J. Chem. Theory Comput. 2014, 10, 534-542. (6) Herce, H. D.; Garcia, A. E.; Darden, T. The Electrostatic Surface Term: (I) Periodic Systems. J. Chem. Phys. 2007, 126, 124106. (7) Pan, C.; Yi, S.; Hu, Z. The Effect of Electrostatic Boundaries in Molecular Simulations: Symmetry Matters. Phys. Chem. Chem. Phys. 2017, 19, 4861-4876. (8) Hu, Z. Infinite Boundary Terms of Ewald Sums and Pairwise Interactions for Electrostatics in Bulk and at Interfaces. J. Chem. Theory Comput. 2014, 10, 5254-5264. (9) Smith, E. R. Electrostatic Potentials in Systems Periodic in One, Two, and Three Dimensions. J. Chem. Phys. 2008, 128, 174104. (10) Ballenegger, V. Communication: On the Origin of the Surface Term in the Ewald Formula. J. Chem. Phys. 2014, 140, 161102. (11) Simulation of Electrostatic Systems in Periodic Boundary Conditions. II. Equivalence of Boundary Conditions. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 1980, 373, 57-66. (12) Perram, J. W.; Petersen, H. G.; De Leeuw, S. W. An Algorithm for the Simulation of Condensed Matter Which Grows as the 3/2 Power of the Number of Particles. Mol. Phys. 1988, 65, 875-893. (13) Kadau, K.; Germann, T. C.; Lomdahl, P. S. Molecular Dynamics Comes of Age: 320 Billion Atom Simulation on Bluegene/L. International Journal of Modern Physics C 2006, 17, 1755-1761. (14) Abraham, F. F.; Walkup, R.; Gao, H.; Duchaineau, M.; Diaz De La Rubia, T.; Seager, M. Simulating Materials Failure by Using up to One Billion Atoms and the World's Fastest Computer: Work-Hardening. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 5783-5787. (15) Fennell, C. J.; Gezelter, J. D. Is the Ewald Summation Still Necessary? Pairwise Alternatives to the Accepted Standard for Long-Range Electrostatics. J. Chem. Phys. 2006, 124, 234104. (16) Steinbach, P. J.; Brooks, B. R. New Spherical-Cutoff Methods for Long-Range Forces in Macromolecular Simulation. J. Comput. Chem. 1994, 15, 667-683. (17) Nguyen, T. D.; Carrillo, J.-M. Y.; Dobrynin, A. V.; Brown, W. M. A Case Study of Truncated Electrostatics for Simulation of Polyelectrolyte Brushes on Gpu Accelerators. J. Chem. Theory Comput. 2013, 9, 73-83. (18) Orsi, M. Comparative Assessment of the Elba Coarse-Grained Model for Water. Mol. Phys. 2013, 112, 1566-1576. (19) Ma, K.; Forsman, J.; Woodward, C. E. Influence of Ion Pairing in Ionic Liquids on Electrical Double Layer Structures and Surface Force Using Classical Density Functional Approach. J. Chem. Phys. 2015, 142, 174704. (20) Toukmaji, A.; Sagui, C.; Board, J.; Darden, T. Efficient Particle-Mesh Ewald Based Approach to Fixed and Induced Dipolar Interactions. J. Chem. Phys. 2000, 113, 10913-10927.



(21) Barker, J. A.; Watts, R. O. Monte Carlo Studies of the Dielectric Properties of Water-Like Models. Mol. Phys. 1973, 26, 789-792. (22) Onsager, L. Electric Moments of Molecules in Liquids. J. Am. Chem. Soc. 1936, 58, 1486-1493. (23) van Gunsteren, W. F.; Berendsen, H. J. C.; Rullmann, J. A. C. Inclusion of Reaction Fields in Molecular Dynamics. Application to Liquid Water. Faraday Discussions of the Chemical Society 1978, 66, 58-70. (24) Wolf, D.; Keblinski, P.; Phillpot, S. R.; Eggebrecht, J. Exact Method for the Simulation of Coulombic Systems by Spherically Truncated, Pairwise R−1 SummaIon. J. Chem. Phys. 1999, 110, 82548282. (25) Zahn, D.; Schilling, B.; Kast, S. M. Enhancement of the Wolf Damped Coulomb Potential: Static, Dynamic, and Dielectric Properties of Liquid Water from Molecular Simulation. J. Phys. Chem. B 2002, 106, 10725-10732. (26) Fanourgakis, G. S. An Extension of Wolf’s Method for the Treatment of Electrostatic Interactions: Application to Liquid Water and Aqueous Solutions. J. Phys. Chem. B 2015, 119, 1974-1985. (27) Wu, X.; Brooks, B. R. Isotropic Periodic Sum: A Method for the Calculation of Long-Range Interactions. J. Chem. Phys. 2005, 122, 044107. (28) Hadley, K. R.; McCabe, C. On the Investigation of Coarse-Grained Models for Water: Balancing Computational Efficiency and the Retention of Structural Properties. J. Phys. Chem. B 2010, 114, 45904599. (29) Reith, D.; Pütz, M.; Müller-Plathe, F. Deriving Effective Mesoscale Potentials from Atomistic Simulations. J. Comput. Chem. 2003, 24, 1624-1636. (30) Moore, T. C.; Iacovella, C. R.; McCabe, C. Derivation of Coarse-Grained Potentials via Multistate Iterative Boltzmann Inversion. J. Chem. Phys. 2014, 140, 224104. (31) Jocelyn, M. R.; John, D. W. Local Molecular Field Theory for the Treatment of Electrostatics. J. Phys.: Condens. Matter 2008, 20, 494206. (32) Hu, Z. Symmetry-Preserving Mean Field Theory for Electrostatics at Interfaces. Chem. Commun. 2014, 50, 14397-14400. (33) Yi, S.-S.; Pan, C.; Hu, Z.-H. Accurate Treatments of Electrostatics for Computer Simulations of Biological Systems: A Brief Survey of Developments and Existing Problems. Chinese Physics B 2015, 24, 120201. (34) Yi, S.; Pan, C.; Hu, L.; Hu, Z. On the Connections and Differences among Three Mean-Field Approximations: A Stringent Test. Phys. Chem. Chem. Phys. 2017, 19, 18514-18518. (35) Remsing, R. C.; Liu, S.; Weeks, J. D. Long-Ranged Contributions to Solvation Free Energies from Theory and Short-Ranged Models. Proceedings of the National Academy of Sciences of the United States of America 2016, 113, 2819-2826. (36) Izvekov, S.; Swanson, J. M. J.; Voth, G. A. Coarse-Graining in Interaction Space: A Systematic Approach for Replacing Long-Range Electrostatics with Short-Range Potentials. J. Phys. Chem. B 2008, 112, 4711-4724. (37) Lafond, P. G.; Izvekov, S. Multiscale Coarse-Graining of Polarizable Models through ForceMatched Dipole Fluctuations. J. Chem. Theory Comput. 2016, 12, 5737-5750. (38) Izvekov, S.; Parrinello, M.; Burnham, C. J.; Voth, G. A. Effective Force Fields for Condensed Phase Systems from Ab Initio Molecular Dynamics Simulation: A New Method for Force-Matching. J. Chem. Phys. 2004, 120, 10896-10913. (39) Kondrat, S.; Georgi, N.; Fedorov, M. V.; Kornyshev, A. A. A Superionic State in Nano-Porous Double-Layer Capacitors: Insights from Monte Carlo Simulations. Phys. Chem. Chem. Phys. 2011, 13, 11359-11366. (40) Borodin, O.; Smith, G. D. Development of Many−Body Polarizable Force Fields for Li-Battery Components: 1. Ether, Alkane, and Carbonate-Based Solvents. J. Phys. Chem. B 2006, 110, 6279-6292. 32 ACS Paragon Plus Environment

Page 32 of 34

Page 33 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(41) Borodin, O.; Smith, G. D. Development of Many−Body Polarizable Force Fields for Li-Battery Applications: 2. Litfsi-Doped Oligoether, Polyether, and Carbonate-Based Electrolytes. J. Phys. Chem. B 2006, 110, 6293-6299. (42) Borodin, O.; Smith, G. D.; Henderson, W. Li+ Cation Environment, Transport, and Mechanical Properties of the Litfsi Doped N-Methyl-N-Alkylpyrrolidinium+TFSI- Ionic Liquids. J. Phys. Chem. B 2006, 110, 16879-16886. (43) Rick, S. W.; Stuart, S. J.; Berne, B. J. Dynamical Fluctuating Charge Force Fields: Application to Liquid Water. J. Chem. Phys. 1994, 101, 6141-6156. (44) Ribeiro, M. C. C. Ionic Dynamics in the Glass-Forming Liquid Ca0.4k0.6no31.4 a Molecular Dynamics Study with a Polarizable Model. Phys. Rev. B 2001, 63, 094205. (45) Schröder, C.; Steinhauser, O. Simulating Polarizable Molecular Ionic Liquids with Drude Oscillators. J. Chem. Phys. 2010, 133, 154511. (46) Yan, T.; Burnham, C. J.; Del Pópolo, M. G.; Voth, G. A. Molecular Dynamics Simulation of Ionic Liquids: The Effect of Electronic Polarizability. J. Phys. Chem. B 2004, 108, 11877-11881. (47) Chang, T.-M.; Dang, L. X. Computational Studies of Structures and Dynamics of 1,3Dimethylimidazolim Salt Liquids and Their Interfaces Using Polarizable Potential Models. The Journal of Physical Chemistry A 2009, 113, 2127-2135. (48) Gu, Y.; Yan, T. Thole Model for Ionic Liquid Polarizability. The Journal of Physical Chemistry A 2013, 117, 219-227. (49) Yu, H.; van Gunsteren, W. F. Accounting for Polarization in Molecular Simulation. Comput. Phys. Commun. 2005, 172, 69-85. (50) Thole, B. T. Molecular Polarizabilities Calculated with a Modified Dipole Interaction. Chem. Phys. 1981, 59, 341-350. (51) Taylor, T.; Schmollngruber, M.; Schröder, C.; Steinhauser, O. The Effect of Thole Functions on the Simulation of Ionic Liquids with Point Induced Dipoles at Various Densities. J. Chem. Phys. 2013, 138, 204119. (52) Elking, D.; Darden, T.; Woods, R. J. Gaussian Induced Dipole Polarization Model. J. Comput. Chem. 2007, 28, 1261-1274. (53) Borodin, O. Polarizable Force Field Development and Molecular Dynamics Simulations of Ionic Liquids. J. Phys. Chem. B 2009, 113, 11463-11478. (54) Aguado, A.; Madden, P. A. Ewald Summation of Electrostatic Multipole Interactions up to the Quadrupolar Level. J. Chem. Phys. 2003, 119, 7471-7483. (55) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Explicit Reversible Integrators for Extended Systems Dynamics. Mol. Phys. 1996, 87, 1117-1157. (56) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695-1697. (57) Wang, W.; Skeel, R. D. Fast Evaluation of Polarizable Forces. J. Chem. Phys. 2005, 123, 164107. (58) Iuchi, S.; Izvekov, S.; Voth, G. A. Are Many-Body Electronic Polarization Effects Important in Liquid Water? J. Chem. Phys. 2007, 126, 124505. (59) Dang, L. X.; Chang, T.-M. Molecular Dynamics Study of Water Clusters, Liquid, and Liquid–Vapor Interface of Water with Many-Body Potentials. J. Chem. Phys. 1997, 106, 8149-8159. (60) Dang, L. X. The Nonadditive Intermolecular Potential for Water Revised. J. Chem. Phys. 1992, 97, 2659-2660. (61) Smith, D. E.; Dang, L. X. Computer Simulations of NaCl Association in Polarizable Water. J. Chem. Phys. 1994, 100, 3757-3766. (62) Starovoytov, O. N.; Borodin, O.; Bedrov, D.; Smith, G. D. Development of a Polarizable Force Field for Molecular Dynamics Simulations of Poly (Ethylene Oxide) in Aqueous Solution. J. Chem. Theory Comput. 2011, 7, 1902-1915. 33 ACS Paragon Plus Environment


(63) Vatamanu, J.; Borodin, O.; Smith, G. D. Molecular Dynamics Simulation Studies of the Structure of a Mixed Carbonate/LiPF6 Electrolyte near Graphite Surface as a Function of Electrode Potential. J. Phys. Chem. C 2011, 116, 1114-1121. (64) Xing, L.; Vatamanu, J.; Borodin, O.; Smith, G. D.; Bedrov, D. Electrode/Electrolyte Interface in Sulfolane-Based Electrolytes for Li Ion Batteries: A Molecular Dynamics Simulation Study. J. Phys. Chem. C 2012, 116, 23871-23881. (65) Bedrov, D.; Borodin, O.; Li, Z.; Smith, G. D. Influence of Polarization on Structural, Thermodynamic, and Dynamic Properties of Ionic Liquids Obtained from Molecular Dynamics Simulations. J. Phys. Chem. B 2010, 114, 4984-4997. (66) Borodin, O.; Smith, G. D.; Douglas, R. Force Field Development and Md Simulations of Poly(Ethylene Oxide)/Libf4 Polymer Electrolytes. J. Phys. Chem. B 2003, 107, 6824-6837.


Page 34 of 34

Application of Screening Functions as Cutoff-Based Alternatives to

Recommend Documents