Effect of Hydrogen Bonds on the Dielectric Properties of Interfacial
Water
Sleeba Varghese,
†
Sridhar Kumar Kannam,
‡,∥
Jesper Schmidt Hansen,
§
and Sarith P. Sathian*
,†
†
Department of Applied Mechanics, Indian Institute of Technology Madras, Chennai 600036, India
‡
Faculty of Science, Engineering and Technology, Swinburne University of Technology, Melbourne, Victoria 3122, Australia
§
Department of Science and Environment, Roskilde University, Roskilde 4000, Denmark
∥
School of Applied Sciences, RMIT University, Melbourne, Victoria 3001, Australia
ABSTRACT: The dielectric constant for water is reduced
under confinement. Although this phenomenon is well known,
the underlying physical mechanism for the reduction is still in
debate. In this w ork, we investigate t he effect of the
orientation of hydrogen bonds on the dielectric properties
of confined water using molecular dynamics simulations. We
find a reduced rotational diffusion coefficient for water
molecules close to the solid surface. The reduced rotational
diffusion arises due to the hindered rotation away from the
plane parallel to the channel walls. The suppressed rotation in
turn affects the orientational polarization of water, leading to a
low value for the dielectric constant at the interface. We
attribute the constrained out-of-plane rotation to originate
from a higher density of planar hydrogen bonds formed by the interfacial water molecules.
■
INTRODUCTION
Confinement induces significant changes in the structural and
dynamical properties of water with most perturbations being
extremely local and observable only within a few molecular
layers close to the surface.
1−3
However, recent studies indicate
that the effect of confinement on the dielectric properties of
water extends well beyond the interfacial regions.
4,5
The
extremely low dielectric constant of nanoconfined water also
suggests a far more significant influence by the interfaces on
the dielectric properties, when the confinement width reduces
to nanometric scales.
5−10
Despite these data available on the
dielectric permittivity of water, knowledge on the physical
origin of the anomalous dielectric behavior of confined water is
still not fully complete. Understanding the underlying physical
mechanism of the dielectric response is essential in the
solvation dynamics of aqueous solutions
11−13
that are primarily
governed by the bulk dielectric properties of the solvent as well
as in the development of nanoscale electromechanical
devices.
14,15
Experimental investigations of the problem are very
complicated. This has led to the development of theoretical
and computer simulation techniques to calculate the molecular
dielectric response in highly confined geometries. An
appropriate theoretical model should at least be able to
reproduce the experimental dielectric constant of homoge-
neous water. The Kirkwood−Fro
hlich fluctuation formula
16,17
has shown satisfactory results on the static dielectric constant
of homogeneous water for different temperature ranges.
18,19
However, this approach has significant limitations when
applied to confined systems. For instance, while calculating
permittivity of an interfacial layer, the Kirkwood−Fro
hlich
fluctuation formula neglects the dipolar contribution of water
molecules away from the layer but is significant due to the
substantial ordering of water molecules near a solid sur-
face.
20,21
Ballenegger and Hansen
22
addressed these issues and
derived a correct fluctuation formula to calculate the dielectric
permittivity of confined polar fluids. They found that under
confinement, the dielectric constant is a space-dependent
dielectric response function with components parallel and
perpendicular to the surface. Later, Bonthuis et al.
23
modified
the fluctuation formula to calculate the dielectric properties of
molecules with higher multipole moments (e.g., water). Both
studies observed the emergence of singularities in the
perpendicular dielectric response function at the molecular
layers close to the solid surface. Further investigations into the
interfacial dielectric properties revealed the existence of an
abnormally high polarization field close to the surface, which
overscreens any external electric field.
22,23
The consequence of
this overscreening phenomenon is the negative values observed
in the dielectric permittivity profile of water close to the solid
surface, which ultimately leads to a low value of the effective
dielectric constant in the interface regions.
24
Experimental and simulation studies show a substantial
decrement in the perpendicular component of dielectric
Received: February 24, 2019
Revised: May 22, 2019
Published: May 23, 2019
Article
pubs.acs.org/Langmuir
Cite This: Langmuir 2019, 35, 8159−8166
© 2019 American Chemical Society 8159 DOI: 10.1021/acs.langmuir.9b00543
Langmuir 2019, 35, 8159−8166
Downloaded by ROYAL LBRY ROSKILDE UNIV LBRY at 05:03:28:501 on June 19, 2019
from https://pubs.acs.org/doi/10.1021/acs.langmuir.9b00543.
permittivity for a few layers of water molecules near the solid
surface.
23−27
Even though this interfacial region with the
reduced dielectric response extends only up to a few
nanometers, under confinement, its influence on the overall
dielectric permittivity becomes quite significant. Molecular
dynamics simulation studies reported an extremely low
effective dielectric constant for water under nanoconfined
geometries.
6,8
A recent experiment on the water con fined
between hexagonal boron nitride layers substantiated this low
dielectric response under nanoconfinement and it was
observed that water regains its bulk permittivity only when
channel widths are higher than 100 nm.
9
This reduced dielectric response under confinement is
modeled on the basis of the capacitor model approach, which
considers the overall permittivity as a series of low-permittivity
interfacial regions and a nominal permittivity bulk region.
9,28
Although the capacitor model explains the overall permittivity
of a confined system, the physical origin of the reduced
dielectric response is still to be pondered. Understanding the
physical mechanism of this reduced dielectric response at the
interfaces requires investigation into the distinctive structural
features of water close to solid surfaces, which can be
rationalized in terms of intermolecular hydrogen bonding of
water. Hence, like most anomalous properties, any deviation in
the dielectric behavior of water, we here conjecture, should
also arise from the restruct uring of its hydrogen-bo nd
network.
29−31
Despitethelargeamountofknowledge
gathered, the influence of hydrogen bonds on the dielectric
permittivity of water is still unclear.
26,27,31,32
In this work, we investigate the orientation of hydrogen
bonds and its correlation with the dielectric permittivity of
confined water. We find a higher density of hydrogen bonds
paral lel to the walls at the solid−water interface when
compared with the bulk liquid. This enhanced planar hydrogen
bond network constrains the rotation of water molecules to the
planes parallel to the channel walls. This extensive network of
planar hydrogen bonds indicates that the reorientation of water
molecules in the direction perpendicular to the channel walls is
energetically unfavorable. The consequence will be low
polarizability of the interfacial water molecules under an
external electric field. The next section provides a detailed
description of the molecular dynamics simulation methods
employed, which is followed by a section on the Results and
Discussion.
■
METHODS
We perform molecular dynamics simulations of a graphene nano-
channel enclosing SPC/E
33
water using the large atomic/molecular
massively parallel simulator
34
package. The bond and angular
vibrations of each water molec ule are constrained using the
SHAKE
35
algorithm. The SPC/E water model has a dipole moment
of 2.35 D, in agreement with the experimental value of 2.9 ± 0.6 D.
36
Graphene sheets are placed in the x−y plane with periodicity in both
directions, and z is the perpendicular or confined direction, as shown
in Figure 1. Liquid-surface chemical reactions and atomic and
electronic polarizability effects are neglected. The graphene walls are
kept as rigid, uncharged, and electrically transparent throughout the
simulations. Each hydrogen atom carries a partial charge of 0.4238e
and oxygen atoms carry a partial charge of −0.8476e, where e is the
elementary charge. The repulsive and attractive parts of van der Waals
interaction are modeled using the pairwise additive Lennard-Jones
(LJ) potential with a cutoff length of 10 Å. The van der Waals
parameters for water−carbon (WC) interaction (ϵ
WC
= 0.09369 kcal/
mol, σ
WC
= 3.19 Å) are taken from Werder et al.
37
The Lennard-Jones
parameters for oxygen−oxygen (OO) interaction are ϵ
OO
= 0.15535
kcal/mol and σ
OO
= 3.166 Å, and Lennard-Jones interaction
coefficients are zero for hydrogen atoms. The short-range Coulombic
interactions with the charged sites are also modeled with the same
truncation cutoff of 10 Å. Long-range electrostatic interactions are
calculated using the Ewald algorithm with a particle−particle-
particle−mesh solver
38
of LAMMPS
34,39
with a relative root mean
square error in the per-atom force calculations below 1 × 10
−4
.
Nonperiodicity in the z-direction is accommodated in the Ewald
algorithm by inserting empty volumes above the channel walls such
that the extended confined dimension is 3 times the actual channel
size and a correction term in the calculation of long-range Coulombic
forces.
34,39
The corrected Ewald algorithm is known as EW3DC.
39
The simulation box dimensions for the graphene channel in the x−y
plane are L
x
= 34 Å and L
y
= 32 Å. Channel width is kept at 40 Å.
Water up to 5 Å from the graphene surface (0−5 and 35−40 Å, which
includes the density peaks near the wall) is considered as interfacial
water.
3
The equilibrium density for water inside a nanochannel is
obtained by immersing our graphene sheets inside a 50 × 50 × 60 Å
3
water reservoir at the ambient conditions of 1.01 bar pressure and
temperature of 300 K. We attained a density of 1.02 g/cm
3
at the
center of the channel, which agrees well with the value reported for
channel widths of similar dimensions.
3
After achieving the density, the
system is equilibrated in the NVT ensemble at 300 K with a Nose
−
Hoover
40,41
thermostat. To calculate the perpendicular component of
the dielectric constant, an electric field of strength E
ext
= 0.1 V/Å is
applied in the direction perpendicular to the walls. The total
simulation time is 12 ns. The first 6 ns is performed without an
electric field and the rest with the field. A time step of 1 fs is used
throughout the simulation. Five independent simulations with
different initial configurations are performed to calculate the statistical
errors.
■
RESULTS AND DISCUSSION
Dielectric Response under Confinement. The static
dielectric constant for a rigid nonpolarizable water model using
Figure 1. (a) Schematic representation of the confined water system (b) Top view of the graphene channel.
Langmuir Article
DOI: 10.1021/acs.langmuir.9b00543
Langmuir 2019, 35, 8159−8166
8160
Ewald’s summation under conducting boundary conditions is
given by
18
M
kTV
1
4
3
2
B
π
ϵ= +
⟨⟩
(1)
where M(t)=∑
i=1
N
μ
i
(t) is the total dipole moment of the
system and μ
i
is the dipole moment of the ith molecule, and k
B
,
T, and V represent Boltzmann’s constant, temperature, and
volume of the system, respectively. We obtain the static
dielectric constant of isotropic, homogeneous water as 72 ± 3,
which is in good agreement with the available data in the
literature.
19
Close to the confining walls, nonzero polarization exists even
at a vanishing electric field. Hence, the net polarization density
upon the application of an external field E
ext
is given by
22
PP P
E 0
ext
Δ
=⟨ ⟩ −⟨ ⟩
(2)
where ⟨...⟩
E
ext
and ⟨...⟩
0
denote the ensemble averages with and
without the applied electric field, respectively. Due to the
planar symmetry of the confined system in the (x, y) plane, the
average polarization has no components parallel to the surface
(i.e., ⟨P⟩ = (0, 0, P(z))). Under constant displacement
boundary conditions the net polarization density from eq 2,
relates to the inverse perpendicular dielectric response function
(ϵ
⊥
−1
(z)) as
22,23,42
z
Pz
E
() 1
()
1
0ext
ϵ=−
Δ
ϵ
⊥
−
(3)
where E
ext
is the external electric field applied perpendicular to
the walls and ϵ
0
is the vacuum dielectric permittivity.
The induced electric field due to the polarization of water in
response to the external electric field E
ext
applied across z is
given by
43
Ez
zz
()
()d
z
p
0
∫
ρ
=
′′
ϵ
−∞
(4)
where ρ(z′) is the charge density distribution perpendicular to
the surface. According to the multipole expansion approach the
polarization field is related to the polarization density as
23
P
zEz() ()
0p
=−ϵ
(5)
Therefore, from eqs 4 and 5 we get
P
zzz() ( )d
z
∫
ρ=−
′
′
−∞
(6)
The local charge density ρ(z)isdefined as
z
LL
qz z()
1
()
xy
i
N
i
i
∑
ρ
δ=⟨ −⟩
(7)
where q
i
is the charge of ith atom and the ensemble average is
replaced by a sample average in a bin, with z as the mid-point
of the slab/bin.
Figure 2a plots the net polarization density calculated from
eqs 2 and 6. We observe the net polarization density profile
(ΔP) to be symmetric even though the polarization density
profiles with (P
E
ext
) and without (P
0
) electric field are anti-
symmetric around the channel center. This behavior arises due
to the alignment of water dipoles in the direction of the applied
electric field (hydrogen atoms point toward the top wall and
oxygen atoms point toward the bottom wall unde r the
influence of an external field). Close to the walls, the
perpendicular dielectric permittivity, given in Figure 2b,
shows large amplitude variations similar to the net polarization
density profile. However, from the inset plot of Figure 2b, we
can observe that the dielectric response reaches the bulk value
ϵ
⊥
−1
≈ 1/71 = 0.014 at the center of the channel matching the
previous observations in the case of confined water.
22,44
Similar to eq 2, the net electric field induced due to the
polarization response under the external electric field E
ext
is
given by
Ez E z E z() () ()
Epp, p,0
ext
Δ
=−
(8)
From Figure 3a, we observe that the magnitude of polarization
field response ΔE
p
at the interfaces is much higher than the
magnitude of the applied field. This results in the over-
screening of the external field by the polarization field. The
electric field screened by the polarization field can be measured
using a dimensionless quantity called the screening factor
defined as
24
Ez E
S
F()/
pext
≡−Δ
(9)
Using eqs 5 and 9, we can modify eq 3 for the perpendicular
dielectric response ϵ
⊥
−1
and write it in terms of the screening
factor as
z
Ez
E
z() 1
()
1SF()
1
p
ext
ϵ=+
Δ
=−
⊥
−
(10)
From Figure 3b, we observe that SF > 1 close to the walls,
affirming the overscreening phenomenon, which results in
negative values of permittivity, as defined by eq 10, for water
near the graphene walls, and concurs with the previous
observations on confined fluids.
22−24,42
Close to the surface,
the permittivity 1/ϵ
⊥
−1
(z) diverges whenever the polarization
field completely screens the external field, i.e., SF = 1,
Figure 2. (a) Polarization density profile with (red) and without
(green) an external field for water confined in a 40 Å channel. The
dashed line represents the variation of the net polarization density
(due to E
ext
= 0.1 V/Å) along the channel width. (b) Inverse dielectric
response of water confined in the 40 Å channel under an external
perpendicular electric field of E
ext
= 0.1 V/Å, calculated from eq 3.
Figure 3. (a) Change in the polarization field and (b) variation of the
screening factor, along the channel width under E
ext
= 0.1 V/Å.
Langmuir Article
DOI: 10.1021/acs.langmuir.9b00543
Langmuir 2019, 35, 8159−8166
8161
indicating the need of a nonlocal description of the dielectric
response especially near the solid−water interfaces, as reported
by Schaaf et al.
45
The effective dielectric constant (ϵ
⊥
eff
)atdifferent regions in
the channel is calculated from the formula given by Schlaich et
al.
44
zz L
L
()d
1
1
L
L
1eff
eff
1
2
Ä
Ç
Å
Å
Å
Å
Å
Å
Å
Å
Å
Å
É
Ö
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
Ñ
∫
ϵ=
ϵ
−+
⊥
−
⊥
⊥
(11)
where L
⊥
eff
is the effective width of the region over which the
effective dielectric constant ϵ
⊥
eff
acts and L = L
2
− L
1
is the
distance enclosed by the region. The difference between the
effective width and actual width of the region is
LL
eff
δ
=−
⊥
(12)
As mentioned earlier (Figure 2b; inset plot), the center of the
channel regains the bulk permittivity of water; hence, with ϵ
⊥
eff
= ϵ
bulk
for L
1
=15ÅtoL
2
= 25 Å and using eqs 11 and 12,we
can estimate the value of δ. Now using this value of δ and eq
11, we calculate the effective dielectric constant at different
regions of water, which is given in Table 1. We observe a
reduced effective dielectric permittivity for confined water, as
reported in previous studies,
9,28
with the least value of
permittivity shown by the interface regions. The low dielectric
constant at the interface indicates that the interfacial water
molecules restrain themselves from aligning along the direction
of the external field. We investigate the reason for this
constrained alignment in the following sections.
Orientation of Water Molecules. Since a detailed
analysis on the orientation of water molecules near planar
surfaces has already been performed,
29,46,47
we herein consider
only the key orientational features that are relevant in the study
of the anisotropic dielectric behavior of water under confine-
ment. To re-examine the orientational preference of water
molecules under confinement, we chose two di fferent angles
formed by a water molecule with the positive direction of the
z-axis, viz., the angle ϕ
μ
, formed by the molecular dipole
moment of water μ
i
and the z-axis, and the angle ϕ
OH
, formed
by the OH bonds of water and the z-axis.
Figure 4 shows the probability density distribution of the
two angles defined above for the interfacial and bulk regions.
Water dipoles and OH bonds in the bulk region have random
orientation, whereas in the interfacial regions, they show
preferential orientation. The peak at cos ϕ
μ
= 0 indicates that
molecular dipoles align parallel to the confining surfaces. The
OH bond distribution shows the presence of two additional
peaks at cos ϕ
OH
= 1 and −1 besides the one at cos ϕ
OH
=0.
This indicates that OH bonds have two preferred orientations:
the first one, with OH bonds nearly parallel to the surface and
the other with OH bonds parallel or antiparallel to the
direction of z-axis.
To substantiate our observations, we also present the
probability distribution of angles along the channel in Figure 5.
As mentioned above, the molecular dipoles show preferential
orientation only close to the surface (Figure 5a). The higher
intensity of cos ϕ
OH
= −1 near the bottom graphene layer and
cos ϕ
OH
= 1 near the top graphene layer indicates presence of
hydrogen atoms directly pointing to the confining walls. The
normally oriented hydrogen atoms, as observed in our study,
are a peculiar aspect of water−hydrophobic interactions.
29
From the orientational analysis, we observe the interfacial
water molecules to most likely orient along the planes parallel
to the walls of the confinement, which results in a substantial
restructuring of the hydrogen bond network at the interfaces.
This preferential orientation of the water molecules is also the
reason for the nonzero polarization observed (Figure 2a) close
to the channel walls even at no external field, suggesting a
spontaneous ferroelectric behavior for the interfacial water
molecules.
48
Reorientation Times of Water Molecules under
Confinement. We analyze the rotational and reorientational
dyna mics of water from the time correlation function C
l
defined as
49
Ct P tuu() ( () (0))
ll
≡⟨ · ⟩
μμ
(13)
where P
l
is the lth Legendre polynomial and u
μ
is a unit vector
parallel to the dipole moment of the molecule. Figure 6a,b
represents the results of ln(C
l
(t)) for l = 1 and 2, at different
regions of the confined water system.
The reorientational time τ
l
is obtained from a linear fitof
ln(C
l
(t)) to
y
xb
1
l
τ
=− +
(14)
where b is some constant. For isotropic diffusive motion, the
rotational diffusion coefficient relates to the reorientation time
as
49,50
Table 1. Effective Out-of-Plane Dielectric Constant at
Different Regions for a 40 Å Channel
a
region (Å) ϵ
⊥
eff
interface [0−5 and 35−40 Å] 5.099 ± 0.443
bulk [15−25 Å] 71
total [0−40 Å] 15.44 ± 2.18
a
Values in square brackets indicate the distance enclosed by the
respective region.
Figure 4. Probability distribution of angular orientation of (a) dipole
moments and (b) OH bonds for bulk and interfacial regions.
Figure 5. Probability distribution of angular orientation of (a) dipole
moments and (b) OH bonds with respect to the distance from
graphene walls.
Langmuir Article
DOI: 10.1021/acs.langmuir.9b00543
Langmuir 2019, 35, 8159−8166
8162
D
ll
1
(1)
l
r
τ
=
+
(15)
The average values obtained for τ
1
, τ
2
, and D
r
are provided in
Table 2. For water in unconfined condition and the bulk region
of the confined system, the values of reorientation times
obtained are in close agreement with the reported values.
51
This contrasts with the values close to the walls, where
molecules reorient themselves on a longer timescale, consistent
with the observations found for water near hydrophobic
substrates.
7,52−54
We conjecture this to arise from the
anisotropic rotation of water molecules around the system
axes near the surface. Since the rotational diffusion coefficient
is inversely proportional to the reorientation time, D
r
decreases
as we move closer to the walls.
Now, we focus on the analysis of reorientation/rotation of
water molecules along different axes. The correlation function
C
l
′ along x, y, and z axes is given by
Ct P tuu() ( () (0))
ll,,
′
≡⟨ ·
⟩
μα μα
(16)
where u
μ,α
is a unit vector parallel to the α-component of the
dipole moment of a molecule and α = x, y,orz. Figure 7a,b
shows ln(C
1
′(t)) along x, y, and z axes for the bulk and
interface, which is fitted to eq 14 to obtain the reorientation
time.
Figure 7c shows reorientation times along different axes for
bulk and interfacial water molecules. In the bulk region, we
observe that τ
xx
≈ τ
yy
≈ τ
zz
whereas in the interfacial region, τ
xx
≈ τ
yy
≫ τ
zz
. This means that in the interface water takes a
much longer time to rotate out of the x−y plane than doing a
rotation in the x−y plane. We also observe that τ
xx,interface
,
τ
yy,interface
> τ
xx,bulk
, τ
yy,bulk
, and τ
zz,interface
≪ τ
zz,bulk
. This suggests
that a significant suppression in the dipole rotation is exhibited
in the direction perpendicular to the confining channel.
Orientation of Hydrogen Bonds under Confinement.
Molecules are considered to be hydrogen-bonded if r
OO′
< 3.5
Å and θ
HOO′
<30°, where r
OO′
is the distance between the
donor and acceptor oxygen atoms, O and O′, and θ
HOO′
is the
angle between the OH bond and the OO′ vector.
51,55,56
The
orientation of hydrogen bonds is analyzed by measuring the
angle made by the hydrogen bonds with the confinement axis.
We here conjecture that the estimation of the planar hydrogen
bonds formed by a water molecule can provide information on
the restructuring of the hydrogen bond network under
confinement. Figure 8 illustrates the criteria to consider a
hydrogen bond as planar.
Figure 9a demonstrates the distribution of hydrogen bonds
per molecule lying in planes parallel to the confining walls (x−
y plane). Planar hydrogen bonds restrict the rotation of water
molecules away from the planes parallel to channel walls. The
interface shows a higher density of planar hydrogen bonds
compared with the bulk. Hence, the reorientation of water
molecules becomes suppressed in the direction perpendicular
to the planes parallel to the confining walls. This explains the
Figure 6. Logarithm of the time correlation functions when (a) l =1
and (b) l = 2 for confined water. The dashed lines show the fitted
curves to the logarithm of time correlation functions.
Table 2. Correlation Times and Rotational Diffusion
Coefficient for Confined and Unconfined Water
region (Å) τ
1
(ps) τ
2
(ps) D
r
(ps
−1
)
interface 5.2 ± 0.7 2.7 ± 0.5 0.06−0.09
bulk 4.5 ± 0.1 1.8 ± 0.1 0.09−0.11
total 4.4 ± 0.1 1.91 ± 0.04 0.09−0.11
liq. wat 4.23 ± 0.09 1.78 ± 0.03 0.09−0.12
liq. wat (Expt) 2−7.5 1.7−2.6
Figure 7. Logarithm of the time correlation functions along x, y, and z
axes of the system for (a) bulk and (b) interface regions. The dashed
lines show the linear fit. (c) Values of reorientation time along
different axes for water under confinement.
Figure 8. Schematic representation of the criteria to consider a planar
hydrogen bond. A maximum of 10° shift from the x−y plane for the
hydrogen bond is tolerated while considering the planarity of H
bonds.
Figure 9. (a) Distribution of planar hydrogen bonds per molecule for
water and (b) probability density distribution of angular orientation of
hydrogen bonds confined in a 40 Å channel. ϕ represents the angle
between the hydrogen bond and z-axis.
Langmuir Article
DOI: 10.1021/acs.langmuir.9b00543
Langmuir 2019, 35, 8159−8166
8163
