Entropic Origin of Ionic Interactions in Polar Solvents

Entropic

Origin

of Ionic

Interactions

in Polar

Solvents

Published

as part of The Journal

of Physical

Chemistry

virtual

special

issue “Pablo

G. Debenedetti

Festschrift”.

Samuel

Varner,

‡

Christopher

Balzer,

‡

and Zhen-Gang

Wang

Cite This:

J. Phys.

Chem.

2023,

127, 4328−4337

Read

Online

ACCESS

Metrics

& More

Article

Recommendations

sı

Supporting

Information

ABSTRACT:

Implicit

solvent models that reduce solvent degrees of

freedom

into effective interaction

potentials

are widely used in the study

of soft materials

and biophysical

systems.

For electrolyte

and

polyelectrolyte

solutions,

coarse-graining

the solvent degrees of freedom

into an effective

dielectric

constant

embeds entropic

contributions

into

the temperature

dependence

of the dielectric

constant.

Properly

accounting

for this electrostatic

entropy is essential

to discern whether

a free energy change is enthalpically

or entropically

driven. We address

the entropic

origin of electrostatic

interactions

in a dipolar solvent and

provide a clarified physical picture of the solvent dielectric

response.

calculate

the potential

of mean force (PMF) between

two oppositely

chargedionsinadipolarsolventusingmolecular

dynamics

anddipolarself-consistent

fieldtheory.Wefindwithbothtechniques

that

thePMFisdominated

bytheentropygainfromthedipolerelease,owingtothediminished

orientational

polarization

ofthesolvent.

We also find that the relative contribution

of the entropy to the free energy change is nonmonotonic

with temperature.

We expect

that our conclusions

are applicable

to a broad range of problems

involving

ionic interactions

in polar solvents.

■

INTRODUCTION

A wide variety of simulation

and theoretical

approaches

utilize

implicit solvent models, where the solvent degrees of freedom

arelumpedintoeffectiveinteractions.

1,2

Treating

thesolventas

a background

medium

can significantly

reduce the computa-

tional cost.

−

In doing so, however,

the solvent degrees of

freedom

become

hidden in effective

interaction

potentials,

which are often specified

in an approximate

manner. Common

methods

inbiological

simulations

includetheaccessible

surface

area (ASA) method and continuum

electrostatic

methods

such

as the generalized

Born model.

−

Dielectric

materials

with polar molecules

respond

to an

electric field through reorientation

of the dipoles.

Generally,

the presence

of an electric field will cause the dipoles to align

and increase the order in the system. The free energy change

during dipole reorganization

is thus composed

of both

energetic

and entropic

contributions.

The energy comes

from the electrostatic

interactions

of the dipoles with the

electricfieldandwitheachother,while theentropyarisesfrom

the changes in the orientation

of the dipoles. For a uniform

dielectric

material,

the electrostatic

entropy

contribution

encapsulated

in the temperature

dependence

of the dielectric

constant.

{

(1)

where

is the temperature,

is the electrostatic

entropy

change due to the application

of an electric field,

is the

Helmholtz

free energy change,

is the dielectric

constant,

the system volume, and

is the electric field. Capturing

the

entropic

contribution

relies on knowing

the temperature

dependence

of the dielectric

constant.

In solutions

containing

ions and charged macromolecules,

the presence

of charged species generates

the electric field that

polarizes

the solvent.

Recently,

this phenomenon

has been

used to explain the apparent

discrepancy

between

experi-

ments

−

and coarse-grained

molecular

dynamics

simula-

tions

−

in describing

the driving force for polyelectrolyte

complex

coacervation.

Chen and Wang showed

that the

Coulomb

potential

used in coarse-grained

implicit

solvent

models inherently

includes

an entropic

contribution,

which

they term the electrostatic

entropy.

By correctly

accounting

for this electrostatic

entropy

contribution

through

the

temperature

dependence

of the dielectric

constant,

they were

able to predict entropy

driven coacervation

from implicit

solvent molecular

dynamics,

in agreement

with experimental

observations.

Further,

they rationalized

this entropic

driving

force as arising from the solvent reorganization

using the

example

of two oppositely

charged ions forming an ion pair.

Received:

January 27, 2023

Revised:

March 30, 2023

Published:

May 9,

2023

Article

pubs.acs.org/JPCB

The Authors.

Published

American

Chemical

Society

4328

https://doi.org/10.1021/acs.jpcb.3c00588

J. Phys.

Chem.

2023,

127, 4328

−

4337

The entropic

contributions

to the potential

of mean force have

also been addressed

in other studies which used molecular

dynamics

and the extended

reference

interaction

site model

(RISM).

However,

these studies were restricted

to water,

which has specific properties

and interactions

with different

types of ions.

The effective

interaction

potential,

or potential

of mean

force (PMF),

between

ions in solution

is a result of the

combined

effects of direct ion

−

ion

interactions

and

interactions

of ions with the solvent. At the most basic level,

one can assume the solvent has a uniform

temperature-

dependent

dielectric

constant.

In the Debye approximation,

for the process of bringing

two ions from infinity to a distance

, we have

q q

)

(2)

q q

)

(3)

q q

)

(4)

where

and

are the charges on the ions,

is the internal

energy,

is a dimensionless

measure

of the dipole

strength,

is the dipole moment,

= 1/

, and

is the

molecular

volume

of the solvent.

Thus, the electrostatic

entropy will dominate

the PMF when

> 1. A similar result

can be obtained

using a more complete

dielectric

theory, such

as Onsager’s

theory.

29,30

While this is a satisfying

result, it is

phenomenological

and does not provide

a clear molecular

picture, since the solvent degrees of freedom

are not explicitly

included.

Inthisstudy,weaddresstheconceptofelectrostatic

entropy

and demonstrate

its generality

by studying

two ions in a

dipolar

fluid. Using dipolar

self-consistent

field theory

(DSCFT)

and molecular

dynamics

simulations,

we analyze

the PMF between

two oppositely

charged,

monovalent

ions

immersed

in a dipolar solvent, explicitly

accounting

for the

solvent degrees of freedom.

We separate

the PMF into its

entropic

and energetic

contributions

to determine

the

conditions

where the electrostatic

entropy dominates.

Finally,

we connect the entropic

driving force to the release of dipoles

as the ions approach

one another,

quantified

through

the

decrease

in the solvent polarization.

While solvent reorganiza-

tion has long been recognized

as a source of entropy

physical processes,

especially

in the context of binding of small

molecules

to proteins

and multivalent

ions to charged

macromolecules

in water,

−

the phenomenon

is often

presented

as solvent-

or system-specific,

where other effects

suchashydrogen

bondingmaydominate.

Ourresultshighlight

the generality

of the electrostatic

entropy and the importance

of properly

accounting

for the solvent degrees of freedom

studying

charge-containing

systems.

■

METHODS

Simulation

Model.

To model two ions in a dipolar fluid,

we use a Stockmayer

fluid model based on work by Shock et

al.

Thesolventparticlespossesspermanent

pointdipoles

their center of mass, while the ions are described

as point

charges

with no dipole (Figure

S1). The nonelectrostatic

nonbonded

potential

energyforallbeadtypesisdescribed

bya

truncated

and shifted Lennard-Jones

(LJ) potential,

{

(

)

1 /6

(5)

where

for all pairs, and

is the distance

between

beads

and

. With a cutoff of 2

1/6

, the LJ potential

is purely repulsive

and

is relatively

unimportant.

For all of

our systems,

we set

= 1.

The electrostatic

interactions

are composed

of charge

−

charge, charge

−

dipole,

and dipole

−

dipole

interactions.

The

standard

Coulomb

potential

describes

the charge

−

charge

interactions

between

the two ions,

q q

(

)

(6)

where

is vacuum permittivity

and

is the charge on each

ion. The charge

−

dipole

and dipole

−

dipole

interactions

are,

respectfully,

given as

(

)

(

)

(7)

where

is the center of mass distance

between

the charge

and the solvent dipole, and

(

)

(

)(

)

(8)

where

μμ

isthecenterofmassdistancebetween

pointdipoles

and

.We reiterate that onlythe ionscarry charge. Alllong-

range electrostatic

interactions

are computed

using a standard

Ewald summation.

The solvent dipoles reorient due to the torque arising from

charge

−

dipole

or dipole

−

dipole

interactions.

We use a

Langevin

thermostat

that takes into account

the angular

degrees of freedom

of the solvent. The LJ size parameter

used to describe

the spherical

diameter

required

to update the

angular velocity via the solvent bead’s moment

of inertia.

Throughout

the work, we consider

a coarse-grained

model

solvent based on water.

Namely,

the number density is

0.03344/Å

= 300 K, with diameter

= 3.0 Å and mass

18.015

mol

. Water has a gas-phase

dipole moment

of 1.85 D;

however,

we vary the magnitude

of the dipole moment

highlight

the role of the dipoles, ranging from

= 0

−

2 D. For

simplicity,

we treat the anion and cation as monovalent

ions

with the same size parameter

as the solvent,

−

For all simulations,

we use reduced units with the length scale

= 3 Å, energy scale 1

= 2.49 kJ/mol, and mass scale 1

18.015 g/mol that give a corresponding

time scale

The PMF between

two ions is calculated

using the adaptive

biasing force (ABF) method.

36,37

The ion separation

distance

isdividedinto8windows

intherangeof0.75

to8.5

. Ineach

window,

the system was equilibrated

for 5

timesteps

(

= 0.005

) and production

of 10

timesteps.

Each simulation

consists

of 5000 solvent

particles,

corresponding

to a

simulation

box length of 17.69

. We use the GPU

38,39

and

Colvars

packages

LAMMPS

for simulations

and

OVITO

for

The Journal

of Physical

Chemistry

pubs.acs.org/JPCB

Article

https://doi.org/10.1021/acs.jpcb.3c00588

J. Phys.

Chem.

2023,

127, 4328

−

4337

4329

visualizations.

Example

LAMMPS

and Colvars

scripts for

calculating

the PMF are available

at https://github.com/

chrisbalzer/Stockmayer-Two-Ions.

Field-Theoretic

Model.

In recent years, several groups

have developed

statistical

field theories that explicitly

account

for solvent polarization

and the dielectric

response.

−

Here,

we adapt a dipolar self-consistent

field theory that was

previously

developed

and used to study ion solvation

energy

and electron transfer reorganization

energy.

−

Alternatively,

we could have used the Ornstein

−

Zernike

integral equation

theory with the hypernetted-chain

(HNC) approximation,

the extended

RISM equation

which have both been used to

study the potential

of mean force between

ions at infinite

dilution.

27,51

−

The purpose of the theoretical

model is not to

compare

or validate the simulation,

but rather to provide an

alternative

approach

for studying

the ionic interaction

dipolar fluids to emphasize

the generality

of the presented

behavior.

We believe any theory which can capture the solvent

orientational

polarization

inthevicinityofions cancapturethe

presented

results, which includes allof the mentioned

theories.

Weconsider

twoionsatfixedseparation

immersed

inadipolar

solvent. The solvent dipole is composed

of a permanent

dipole

andaninduceddipole

. Theinduceddipoleisrelatedtothe

electronic

degrees

of freedom.

This contribution

is not

essential

for this work, but we include it for completeness

it does not add much complexity

to the theory or calculations.

The ions are modeled

as Gaussian

smeared

charges inside a

spherical

solute cavity, where the solvent is excluded.

The use

of smeared

charges is for convenience,

to avoid the diverging

self-energy.

The cavities

are modeled

using a spherically

symmetric

cavity function.

−

The ions are treated as fixed

external charges to the solvent in this field theory. See Figure

S2foranexampleofthesystemsetup.Themicroscopic

charge

density in the system is given by the contributions

from the

two ions and the solvent dipoles.

(

)

(

)

(

)

(

)

(9)

Here,

(

)

(

)

(

)

is the total charge density of the

ions, modeled

as Gaussian

smeared

charges,

{

z e

(

)

exp

3/ 2

(10)

where

are the ion valencies,

are the radii of the charge

spread,

are the ion positions,

and

is the elementary

charge.

(

)

and

(

)

are the charge densities

due to the permanent

and induced

dipoles of the solvent, respectively.

The charge

density due to theorientational

and electronic

contributions

the solvent can be expressed

in terms of their dipole moments

(

)

(

)

(11)

(

)

(

)

(12)

where

is the number

of solvent molecules,

is the

permanent

dipole moment

onmolecule

, and

is theinduced

dipole moment

on molecule

. In the energy, we only consider

the Coulomb

interactions

between

all charges, and a harmonic

penalty associated

with the induced dipoles

(

)

(

)

(

)

| |

(13)

where

is the polarizability.

We work in the grand canonical

ensemble

with volume

, temperature

, and solvent chemical

potential

. We assume the system is incompressible.

The

grand canonical

partition

function

is given by

(

)

(

)

[

]

(14)

where

is the solvent density operator,

is the volume

fraction occupied

by the ions, and

is a factor similar to the

cube of the thermal wavelength,

which has no thermodynamic

consequence.

The solvent

density

operator

is given by

(

)

(

)

. The local volume fraction due to the

two ions, is modeled

by the superposition

of two cavity

functions,

(

)

tanh

(15)

where

is the diameter

of the two ions,

is a positive

parameter

for shifting the boundary

of the cavity, and

is a

positive

parameter

for tuning the width of the interface

between

the ion and the solvent. In practice,

any reasonable

choice of

and

will yield the same qualitative

results. We

choose

= 0.01 for a rapid but continuous

decay near the ion

boundary,

and we choose

= 0.95 such that the ion volume

fraction practically

decays to zero by

We take advantage

of the Fourier representation

of a delta

functional

to introduce

the incompressibility

field

(

)

(

)

exp

(

)

(

)

{

}

[

] =

[

]

(16)

To turn the partition

function

into a field theory, we introduce

coarse-grained

density fields

and

through the identities

(

)

(

)

exp

(

)

(

)

{

}

[

]

[

]

(17)

(

)

(

)

exp

(

)

(

)

{

}

[

]

[

]

(18)

where

and

are auxiliary

fields introduced

through the

Fourier representation

of the delta functional.

In eqs 16

−

18,

The Journal

of Physical

Chemistry

pubs.acs.org/JPCB

Article

https://doi.org/10.1021/acs.jpcb.3c00588

J. Phys.

Chem.

2023,

127, 4328

−

4337

4330

the notation

[···]

denotes functional

integration.

Applying

these identities

to the partition

function

results in

}

exp

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

| |

[

] +

[

]

[

]

(19)

where

(

)

(

)

(

)

(

)

. Through

these identity

transformations,

the particle

−

particle

interactions

are turned

into particles

interacting

with fluctuating

fields. The particle

degrees of freedom

can be easily integrated

out to yield the

single-particle

partition

function,

[

], given by

{

}

ivw

sin(

)

exp

(

)

(

)

3 /2

(20)

The transformed

grand canonical

partition

function

is now

written as

(21)

where the Hamiltonian,

[

], is a

functional

of the field variables

only and is given by

(

)

(

)

(

)

(

)

(

)

(

)

(

) 1

(

)

[

]

[

]

(22)

The integrals

in the partition

function

cannot be evaluated

closed form, so we use the saddle-point

configuration

as an

approximation

to the full partition

function.

The saddle point

is found by extremizing

the Hamiltonian

with respect to all the

fieldvariables.

Thegrandfreeenergyisgivenby

−

≈

−

Ξ*

, where

Ξ*

is the saddle-point

contribution

to the

partition

function.

Extremizing

the free energy results in a set

of self-consistent

equations

that can be solved to find the

equilibrium

field configurations,

and therefore

observables

such as the free energy and the electric field. Since the saddle-

pointsforthe

fieldslieontheimaginary

axis,toavoidtheuse

of imaginary

numbers,

we define real potentials

−

and

βφ

(

at equilibrium).

The set of self-

consistent

equations

obtained

fromextremizing

withrespect

(and

, and

respectively,

is given by

(

)

sinh(

(

) )

(

)

(

)

(

/ 2)

(

)

(23)

(

)

(

))

(

)

(

)

·[

] =

(24)

(

)

(

)

(

)

(

/ 2)

(

)

(25)

(

)

(

)

[

]

(26)

where

(

) = [1/tanh(

)

−

] sinh(

. Additionally,

have defined the orientational

electric susceptibility,

, and

the electronic

electric susceptibility,

. For convenience,

The

Poisson equation,

eq 24, is written in a way that separates

the

ion charge density from the bound solvent charge density

(reflected

and

). Equation

24 can be rewritten

(

)

(27)

where

−

[1 +

(

) +

(

)]

∇

(

)

is the electric

displacement.

Conveniently,

the electric displacement

can be

determined

solely from the location

of the free charges.

simplify

the solution

to eq 27 by assuming

that the total

electric displacement

is the superposition

of the electric

displacement

due to each of the ions. (The superposition

principle

is used as an approximation

to avoid solving the full

Poisson

equation.

In general it does not hold exactly for

nonlinear

dielectrics.

Theapproximation

becomes

lessaccurate

for strong interactions,

and in particular

when the ions are at

close approach.

) The displacement

at a distance

from an

isolated Gaussian

smeared

charge

{

(

)

erf

exp

(28)

where

−

, and

is the unit vector in the direction

The equilibrium

electric field is found by iterating

(

) =

−

(

)

∇

(

)

usingeqs23, 25, and26tocalculate

thedielectric

function

(

) =

[1 +

(

) +

(

)]. For all self-consistent

calculations,

we iterate until the maximum

difference

in the

electric field between

two consecutive

steps is less than 10

−

The equilibrium

free energy can be simplified

by plugging

the saddle-point

equations

into the Hamiltonian,

}

(

)

(

))

(

)

(

)

sinh(

(

) )

(

)

(

)

ln(

(

))

{

[

]

(29)

where

(

) = 1

−

(

) is the coarse-grained

solvent volume

fraction.

The system is rotationally

symmetric

about the axis

connecting

the ions, allowing

the integral in eq 29 to be

calculated

in cylindrical

coordinates

with variation

restricted

the

-plane. We use a large enough domain to avoid cutting

off the long-range

tail of the electric field from the ions.

When calculating

the PMF, the quantity

of interest is the

difference

between

the free energy at separation

and that at

infinite separation.

We call this difference

ΔΩ

(

)

−

(

∞

). The entropic

contribution

ΔΩ

can be calculated

with the derivative

{

(30)

Note that the change in grand free energy is equal to the

change in Helmholtz

free energy for this process

as the

chemical

potential

and particle number do not change when

bringing

the ions together

at a fixed

and

due to the

The Journal

of Physical

Chemistry

pubs.acs.org/JPCB

Article

https://doi.org/10.1021/acs.jpcb.3c00588

J. Phys.

Chem.

2023,

127, 4328

−

4337

4331

incompressibility

constraint.

Mathematically,

we can write the

following

(

)

(

)

(

)

[

]

(31)

where the subscript

∞

denotes

a quantity

at infinite ion

separation.

We then decompose

ΔΩ

into entropic

and

energetic

contributions

by calculating

with eq 30 and

internal energy with

ΔΩ

■

RESULTS

AND

DISCUSSION

We start with a discussion

of the behavior

of the PMF for

bringing

two oppositely

charged monovalent

ions together

in a

dipolar solvent; this is shown in Figure 1. Both panels show

that increasing

the dipole moment

of the solvent decreases

the

attraction

between

the two ions. This behavior

can be easily

understood

as due to the increase

in the effective

dielectric

constant

with the dipole moment.

The diverging

behavior

below

= 1

in Figure 1b is due to the diverging

LJ potential

from simulation.

Above

= 1

, we see close agreement

between

theDSCFTandsimulation

resultsatlowto moderate

dipole moment;

a direct comparison

is given in Figure S3 in

the Supporting

Information.

At higher dipole moments

(insets

of Figure 1), oscillations

appear in the PMFs for simulations,

Figure

PMFsforvariousdipolemoments

calculated

from(a)DSCFTand(b)molecular

dynamics

simulation

with

=3Å,

=300K,

−

= 30 Å

. The insets of both panels show PMFs for larger dipole moments.

For reference,

= 1.85D corresponds

to the gas-phase

dipole

moment

of water.

Figure

PMFs decomposed

into their energetic

and entropic

contributions

with solvent dipoles of (a, b)

= 0 D and (c, d)

= 1 D. The PMFs

are calculated

from DSCFT in (a) and (c) and MD simulation

in (b) and (d). Other parameters

are the same as in Figure 1.

The Journal

of Physical

Chemistry

pubs.acs.org/JPCB

Article

https://doi.org/10.1021/acs.jpcb.3c00588

J. Phys.

Chem.

2023,

127, 4328

−

4337

4332