In the late 1970s, Hubbard and Onsager predicted that adding salt to a polar solution would result in a reduced dielectric permittivity that arises from the unexpected tendency of solvent dipoles to align opposite to the applied field. Here we develop a novel non-equilibrium molecular dynamics simulation approach to determine this decrement accurately. Using a thermodynamic consistent all-atom force field we show that for an aqueous solution containing sodium chloride around 4.8 mol/l, this effect accounts for 12% of the total dielectric permittivity. The dielectric decrement can be strikingly different if a less accurate force field for the ions is used. Using the widespread GROMOS parameters, we observe in fact an increment of the dielectric permittivity rather than a decrement, caused by ion pairing and introduced by a too low dispersion force.

It is well known that addition of salt to a polar liquid changes its dielectric properties, in particular by lowering its static dielectric permittivity.1 Much less is known about the actual physical mechanisms which contribute to the decrement. Hubbard and Onsager2 first proposed that the motion of ions, dragged by an external electric field, should induce a depolarization of the solvent molecules, by orienting their dipole moments opposite to the electric field. This effect would be genuinely kinetic in nature, not having any static equivalent.3–5 However, the existence of this effect has never been proven experimentally, because it is masked by the decrement caused by dielectric saturation. In other words, solvent molecules in the vicinity of an ion are strongly polarized by its local field,6–8 and do not take part in the usual dielectric relaxation process, therefore contributing to the reduction of the dielectric permittivity. Molecular dynamics simulations represent an ideal tool to investigate this effect, which was first addressed in simulations by Chandra,9 but so far a coherent framework for the calculation of the Hubbard-Onsager decrement has not been proposed suited for equilibrium and out-of-equilibrium simulations. We use the out-of-equilibrium approach to estimate with high accuracy the static limit of the kinetic decrement, showing that there is indeed a sizable kinetic contribution, even though the Hubbard-Onsager continuum theory overestimates it. In addition, we scrutinize directly the impact of long-living ion pairs on the dielectric spectrum of the solution. Finally, we show that the ionic current correlation results in a dielectric increment. This kinetic contribution is often disregarded, but can counterbalance the Hubbard-Onsager kinetic decrement, even in the absence of ion pairing.

Several approaches have been proposed to date to compute the dielectric spectrum ε(ω) in presence of free ions (e.g., Refs. 10,11). The common denominator of their theoretical background is the separate coupling of the dipole vector of neutral moieties and of the current of charged ones to the electrostatic scalar and vector potentials, respectively. As only the curl-free part of the electric field survives in classical simulations, all these approaches are equivalent12 to calculating the dielectric spectrum from that of the conductivity σ(ω) as ε(ω) = 1 + 4πiσ(ω)/ω, where the quantity 4πiσ(ω)/ω is sometimes called the generalized dielectric function and denoted by Σ(ω). The conductivity spectrum, in turn, can be computed in the linear response regime from the time correlation current13 σ(ω) = β⟨jT(t)jT(0)⟩ω/3V of the total electric current jT, generated by both partial charges on solvent molecules and by ionic charges. The symbols ⟨…⟩ and ⟨…⟩ω denote the canonical average and its Fourier-Laplace transform,

$\langle f \rangle _\omega = \int _0^\infty \langle f \rangle \exp (i \omega t) \mathrm{d}t$
fω=0fexp(iωt)dt⁠, respectively.

The total current can be split arbitrarily into different contributions. For aqueous solutions of atomic ions, however, it is natural to decompose it into solvent contributions, jw, and ionic contributions, ji. The two currents stem from motion of the partial charges of water molecules and from the ionic charges, respectively. This splitting has a particularly important meaning in the context of dynamical contributions. In the linear response relation Δja(ω)∝⟨ja(t)jb(0)⟩ω = Cab(ω) (where subscript a and b are completely general), the current jb is the derivative of the dipole moment that couples to the perturbing field in the Hamiltonian. Hence, for example, in an aqueous solution of ions, the correlation Cww represents the response of water in a fictitious case where only water molecules – and not the ions – are coupled to the electric field. Similar considerations apply to the other possible permutations Cwi, Ciw, and Cii. Evidently, the zero-frequency limit of Δεiw(ω) = 4πβiCwi(ω)/(3ωV) is precisely the first dynamic mechanism identified by Hubbard and Onsager, namely, the depolarization of water molecules due to the motion of ions. The symmetry between Cwi and Ciw is nothing else but an instance of Onsager reciprocal relations, and shows directly that the second dynamical mechanism identified by Hubbard and Onsager (namely, a reduction in ion conductivity due to the rotation of water molecules) leads to the same contribution to the static permittivity. The total kinetic decrement is therefore ΔεHO = 2Δεwi(0). The static dielectric permittivity can be written as ε = 1 + ΔεwT + ΔεiT = 1 + Δεww + Δεii + 2Δεwi, where here and in the following the absence of a frequency dependency implies that the limit to zero frequency has to be taken.

We have investigated the dielectric spectrum of aqueous solutions of sodium chloride, performing extensive MD simulations. The systems were composed of 1621 water molecules (SPC/E potential14) and 160 ion pairs, corresponding to a salt concentration of about 4.7–4.9 mol/l, depending on the average volume yielded by the use of different potentials for the ions. Distinct simulations were performed using the GROMOS15 and the Kirkwood–Buff force field (KBFF)16 force fields (see Table I for the parameters). The ensemble of choice is the isothermal/isobaric one at 298 K and 1 atm, realized using the Nosè–Hoover17,18 and Parrinello–Rahman coupling schemes.19 The electrostatic interaction has been computed using the Particle Mesh Ewald method in its smooth variant20 with tin-foil boundary conditions, 4th order spline interpolation on a grid with 0.12 nm spacing, and a relative strength of the interaction at 0.9 nm of 10−5. Both the Lennard-Jones and short-range part of the electrostatic potentials were switched smoothly to zero using a fourth-degree polynomial in the range from 0.9 to 1.2 nm. Several 8 ns long runs with different initial conditions were performed, up to a cumulative time of about 4 μs for each system. Water and ionic currents were saved to disk every time step (1 fs) for off-line analysis. Current correlations were calculated on the complete 8 ns datasets using a fast Fourier transform based approach.21–23 The correlations computed on each dataset were eventually averaged, and spectra were calculated using delays of up to 2 ns of the averaged correlations, as at longer delay times the noise in the correlation degraded the quality of the spectrum noticeably. No tapering or padding procedure has been used in order to avoid any bias in the spectral features, especially in those at low frequencies.

Table I.

Lennard Jones interaction parameters.a

ModelAtomε (kJ/mol)σLJ (nm)Charge (e)
KBFF Na 0.320a 0.245 
  Cl 0.470 0.440 −1 
GROMOS Na 0.062 0.258 
  Cl 0.446 0.445 −1 
SPC/E 0.650 0.316 −0.8476 
  0.0 0.0 0.4238 
ModelAtomε (kJ/mol)σLJ (nm)Charge (e)
KBFF Na 0.320a 0.245 
  Cl 0.470 0.440 −1 
GROMOS Na 0.062 0.258 
  Cl 0.446 0.445 −1 
SPC/E 0.650 0.316 −0.8476 
  0.0 0.0 0.4238 
a

Geometric mixing rules for both σLJ and ε are employed.

b

The interaction energy between oxygen and sodium atoms has to be scaled by a factor 0.75 in KBFF.

At first, we investigated the solution using the GROMOS force field for Na and Cl ions,15 obtaining the spectra reported in Fig. 1. The spectrum has been split into the contributions coming from water CwT and from the ions CiT, where the index T again denotes the total current. In addition, the positive part of the cross correlations Cwi has also been reported. A clear, main relaxation peak of water (left panel of Fig. 1) appears at ω ≃ 78 GHz (τ = 14.6 ps), but the ionic contribution shows the appearance of an important, unexpected relaxation process at ω ≃ 14 GHz (τ = 70 ps, d.c. contribution subtracted). A simultaneous fit of the real and imaginary part of the spectra to the phenomenological Cole-Cole functional form augmented by a conductivity term,

\begin{equation}\epsilon (\omega )=\epsilon _\infty + \frac{\Delta \epsilon }{1+(i\omega \tau )^{\alpha }}+\frac{4\pi i}{\omega }\sigma ,\end{equation}
ε(ω)=ε+Δε1+(iωτ)α+4πiωσ,
(1)

led to the best fit parameters (amplitude Δε, relaxation time τ exponent α and conductivity σ) reported in Table II (cf. the experimental results in Ref. 24, supplementary material sample 12a for concentration c = 5.0105 mol/l and T = 298.15 K: Δε = 43.48, τ = 6.29, α = 0.801, ε = 5.65).

FIG. 1.

Contributions to the dielectric spectrum in the GROMOS potential case. The water (left panel) an ions (right panel) contributions (squares), and Cole-Cole fit to their real (dashed line) and imaginary (dotted-dashed) parts. In addition, the real part of the ion-ion self-correlation is shown in the left panel (circles) and the water-ion cross correlation is shown in the right panel (triangles).

FIG. 1.

Contributions to the dielectric spectrum in the GROMOS potential case. The water (left panel) an ions (right panel) contributions (squares), and Cole-Cole fit to their real (dashed line) and imaginary (dotted-dashed) parts. In addition, the real part of the ion-ion self-correlation is shown in the left panel (circles) and the water-ion cross correlation is shown in the right panel (triangles).

Close modal
Table II.

Several spectral properties.

 GROMOSKBFFKBFF+pairs
ΔεwT 42.1 ± 0.3 27.4 ± 0.1 39.9 ± 0.3 
ΔεiT 18.0 ± 2.0 … 36.0 ± 4.0 
Δεwi 6.0a −1.7 ± 0.1a 5.9b 
ΔεHO, cont −6.2 −8.0 −8.6 
τwT/ps 14.6 ± 0.2 10.6 ± 0.1 17.1 ± 0.3 
τiT/ps 70.0 ± 1 … 110 ± 10 
αwT 0.820 ± 0.005 0.830 ± 0.003 0.800 ± 0.005 
αiT 0.91 ± 0.01 … 0.94 ± 0.01 
σ/(S/m) 5.1 ± 0.1 10.2 ± 0.1 5.4 ± 0.1 
 GROMOSKBFFKBFF+pairs
ΔεwT 42.1 ± 0.3 27.4 ± 0.1 39.9 ± 0.3 
ΔεiT 18.0 ± 2.0 … 36.0 ± 4.0 
Δεwi 6.0a −1.7 ± 0.1a 5.9b 
ΔεHO, cont −6.2 −8.0 −8.6 
τwT/ps 14.6 ± 0.2 10.6 ± 0.1 17.1 ± 0.3 
τiT/ps 70.0 ± 1 … 110 ± 10 
αwT 0.820 ± 0.005 0.830 ± 0.003 0.800 ± 0.005 
αiT 0.91 ± 0.01 … 0.94 ± 0.01 
σ/(S/m) 5.1 ± 0.1 10.2 ± 0.1 5.4 ± 0.1 
a

From NEMD simulation.

b

Value of Σ(ω) at the smallest frequency.

The contribution coming from the cross-correlation Cwi is large and positive, in open contrast to the prediction of a kinetic decrement by the Hubbard-Onsager mechanism. In the continuum approximation, i.e., in the limit of infinite dilution and large ionic radius, this decrement can be estimated25 as ΔεHO, cont = −8πτwTσ(ε0 − ε)/3ε0, where σ, ε0, and ε are the static conductivity, static permittivity, and infinite frequency dielectric permittivity of the solution, and τwT is the (Debye-like) relaxation time of the solvent. Note that the single Debye relaxation in the Hubbard-Onsager model is an approximation, as the solvent relaxation in this case would be better described by either a superposition of Debye processes or by a Cole-Cole function.24,26

In the ionic contribution to the spectrum (right panel of Fig. 1), however, a feature which is compatible with a Cole-Cole process with a characteristic time of about 70 ps is present. This suggests that ion pairing is taking place on this relatively large timescale which is consistent with earlier reports that the effective attraction between solvated Na and Cl ions in the GROMOS force field is too strong.27 Sodium chloride, instead, is a strong electrolyte which is generally believed not to lead to ion pairing on such a timescale.28 

To test our hypothesis, we simulated the NaCl solution using the KBFF,16 which is more accurate than the GROMOS one in reproducing structural and energetic properties of solvated sodium chloride. The water contribution is qualitatively similar to the GROMOS case, but the ion contribution differs dramatically. The prominent ion relaxation process at about 14 GHz (70 ps), which characterized the GROMOS case, disappears completely, and the zero-frequency contribution becomes so small that it can hardly be estimated (see Fig. 2).

FIG. 2.

Contributions to the dielectric spectrum in the KBFF potential case. The water (left panel) an ions (right panel) contributions (squares), and Cole-Cole fit to their real (dashed line) and imaginary (dotted-dashed) parts. I addition, the real part of the ion-ion self-correlation is shown in the left panel (circles) and the negative part of the water-ion cross correlation is shown in the right panel (triangles).

FIG. 2.

Contributions to the dielectric spectrum in the KBFF potential case. The water (left panel) an ions (right panel) contributions (squares), and Cole-Cole fit to their real (dashed line) and imaginary (dotted-dashed) parts. I addition, the real part of the ion-ion self-correlation is shown in the left panel (circles) and the negative part of the water-ion cross correlation is shown in the right panel (triangles).

Close modal

To check explicitly that the phenomenon underlying the relaxation observed in the GROMOS case is ion pairing, we performed another series of simulations using the KBFF parameters. This time, however, we created artificial ion pairs out of 50% of the ions in solution by introducing a stiff bond between the Na and Cl ions with a bond length corresponding to the first peak in the radial distribution function of the non-bonded case. The effect of persistent pairing on the dielectric relaxation of ions modeled using KBFF can be appreciated in Fig. 3. A relaxation process appears again at low frequencies (ω ≃ 9 GHz, τ = 110 ps), qualitatively very similar to the one observed using GROMOS force field, a clear indication that the process has to be attributed to ion pairing. This also shows, in particular, that it is the contact ion pair mechanism which is involved in the present case, as the first peak of the radial distribution function is located slightly closer than 0.3 nm, and therefore a water molecule cannot be accommodated in between the two ions. In this respect, KBFF proves to be superior to GROMOS not only in describing the structural properties, but also the dynamical ones, as it correctly reproduces the absence of ion pairs that are stable on time scales longer than 10–100 ps.

FIG. 3.

Contributions to the dielectric spectrum in the KBFF potential case, when ion pairs are enforced. The water (left panel) an ions (right panel) contributions (squares), and Cole-Cole fit to their real (dashed line) and imaginary (dotted-dashed) parts. I addition, the real part of the ion-ion self-correlation is shown in the left panel (circles) and the water-ion cross correlation is shown in the right panel (triangles).

FIG. 3.

Contributions to the dielectric spectrum in the KBFF potential case, when ion pairs are enforced. The water (left panel) an ions (right panel) contributions (squares), and Cole-Cole fit to their real (dashed line) and imaginary (dotted-dashed) parts. I addition, the real part of the ion-ion self-correlation is shown in the left panel (circles) and the water-ion cross correlation is shown in the right panel (triangles).

Close modal

Nevertheless, the failure of the GROMOS force field in reproducing the spectrum of aqueous NaCl solutions has provided new important knowledge. On the one hand, this is the first direct observation of the relaxation process associated with ion pairing, with a characteristic time of about 100 ps, well separated from the water relaxation time (approximately 10 ps). On the other hand, these results show that kinetic contributions can lead not just to the Hubbard-Onsager decrement, but also to kinetic increments, which arise from the completely different mechanism of pairing. It is interesting to investigate the origin of the excessive pairing in the GROMOS force field. The small Lennard-Jones interaction energy of sodium (0.062 kJ/mol) suggests a weak attraction between ions. Paradoxically, this results in excessive pairing, as also the hard-core repulsion of the ion is small, allowing two oppositely charged particles to come so close that strong electrostatic attraction sets in, favoring the pairing.

The reliability of the KBFF results offers the opportunity to test quantitatively the Hubbard-Onsager continuum theory2,29 and the Hubbard-Colonomos-Wolynes molecular theory25 for the kinetic decrement. Obtaining a precise estimate of the kinetic decrement from the analysis of the spectrum Cwi(ω), however, can be computationally quite demanding, especially if the static contribution is small. In the present case, in fact, the fluctuations of the spectrum of the cross-correlation function at low frequencies did not allow us to perform a reliable extrapolation to ω = 0. The interpretation in terms of linear response of the correlation Cwi suggests, however, that it is possible to compute the decrement using a non-equilibrium MD (NEMD) approach. This consists in applying a constant (fictitious) electric field acting only on ions, and in measuring the induced polarization in water by the ion flow. This method realizes directly the thought experiment proposed by Hubbard and co-workers25 to explain the first dynamic mechanism for the decrement. This NEMD approach has obvious advantages over the linear-response approach: relatively high fields can be used to increase the signal-to-noise ratio, and no extrapolation procedure is needed. Indeed, a 500 ps long run with an electric field of about 0.25 V/nm is enough to determine the dynamic contribution Δεwi = −1.7 with an accuracy of 5%. This term is only half of the contribution to the total dielectric decrement, which therefore amounts to −3.4. Two separate simulations with a different applied voltage were used to control the validity of linear response assumption.

The other kinetic contribution to the static permittivity is due to the ion-ion current correlation, but unlike the water-ion cross term this cannot be computed using a NEMD approach, since the ionic cloud is infinitely polarizable. As it is seen from the spectrum of CiT in Fig. 2, the scatter of the data at low frequencies is rather large. This makes a precise extrapolation to zero frequency very difficult, but a rough estimate suggests Δεii ≃ 2−5. Such a contribution would in large part counterbalance the total kinetic decrement 2Δεwi.

Coming back to the ion-water contribution, we now compare the simulation results with the available theoretical models. The continuum theory2,29 requires as input parameters the main dielectric relaxation time and the conductivity of the solution. Using the values computed by the MD simulation, the continuum model predicts for the KBFF case a considerably larger decrement ΔεHO, mol ≃ −8 (to be compared with 2Δwi = −3.4). The situation does not improve if the molecular theory developed by Hubbard, Colonomos, and Wolynes25 is used. For a system with many ions in solution the molecular theory estimate of the dielectric decrement is

\begin{equation}\Delta \epsilon _{\text{HO},mol}= -8\pi \frac{\tau _{wT}}{3 e} \left\langle \sum _{ij} R_i \sigma _i \bm{\mu }_j \cdot \mathbf {r}_{ij} / r_{ij}^3 \right \rangle ,\end{equation}
ΔεHO,mol=8πτwT3eijRiσiμj·rij/rij3,
(2)

where μj is the dipole moment of the jth water molecule, located at the relative position rij from the ith ion of radius Ri and specific conductance σi (i.e., the contribution to the total conductance stemming from the single ion current). Applying Eq. (2) to the KBFF case leads to ΔεHO, cont ≃ −19, which is in absolute value even larger than ΔεHO, cont. This is because ΔεHO, mol is very sensitive to the local structure of water around each ion, due to the 1/rij term in the statistical average.

In conclusion, we made a coherent derivation of the Green–Kubo formulas for the two mechanisms identified by Hubbard and Onsager as responsible for the kinetic dielectric decrement, and computed the different (ionic, dipolar, cross) contributions to the dielectric spectrum of sodium chloride aqueous solutions. Using a highly accurate non-equilibrium approach, we showed that the phenomenon of kinetic decrement can account for a considerable fraction of the static dielectric permittivity in salt solutions. However, both the continuum and molecular theories for the dielectric decrement fail describing the decrement at the present concentrations. The reason for this failure has to be sought in the inadequacy of the infinite dilution approximation. The high salt density screens the electrostatic interaction, therefore reducing the efficiency of the Hubbard-Onsager mechanism. This has been confirmed by a supplementary simulation at lower concentration (60 ion pairs and same number of water molecules), for which the kinetic contribution, Δεwi = −2.3, is higher (in absolute value), despite the lower conductivity of the solution. In the limit of a continuum charged background, in fact, no Hubbard-Onsager mechanism can take place.

M.S. acknowledges support from the European Community's Seventh Framework Programme (FP7-PEOPLE-2012-IEF) funded under Grant No. 331932 SIDIS. S.S.K. acknowledges support from RFBR Grant Nos. mol-a 1202-31-374 and mol-a-ved 12-02-33106, from the Ministry of Science and Education of RF 2.609.2011 and, from Austrian Science Fund (FWF): START-Projekt Y 627-N27. A.A. and C.H. acknowledge support by the cluster of excellence SimTech of the University of Stuttgart. The authors thank Friedrich Kremer, Christian Schröder, and Othmar Steinhauser for useful discussions.

1.
Broadband Dielectric Spectroscopy
, edited by
F.
Kremer
and
A.
Schönhals
(
Springer
,
Berlin
,
2003
).
2.
J. B.
Hubbard
and
L. J.
Onsager
,
J. Chem. Phys.
67
,
4850
(
1977
).
3.
P. G.
Wolynes
,
Ann. Rev. Phys. Chem.
31
,
345
(
1980
).
4.
U.
Kaatze
,
J. Phys. Chem.
91
,
3111
(
1987
).
5.
U.
Kaatze
,
J. Solution Chem.
26
,
1049
(
1997
).
6.
J.
Aragones
,
L.
MacDowell
,
J.
Siepmann
, and
C.
Vega
,
Phys. Rev. Lett.
107
,
155702
(
2011
).
7.
A.
Levy
,
D.
Andelman
, and
H.
Orland
,
Phys. Rev. Lett.
108
,
227801
(
2012
).
8.
A.
Levy
,
D.
Andelman
, and
H.
Orland
,
J. Chem. Phys.
139
,
164909
(
2013
).
9.
A.
Chandra
,
J. Chem. Phys.
113
,
903
(
2000
).
10.
J. M.
Caillol
,
D.
Levesque
, and
J. J.
Weis
,
J. Chem. Phys.
85
,
6645
(
1986
).
11.
C.
Schröder
,
T.
Rudas
, and
O.
Steinhauser
,
J. Chem. Phys.
125
,
244506
(
2006
).
12.
M.
Sega
,
S. S.
Kantorovich
,
A.
Arnold
, and
C.
Holm
, “
On the calculation of the dielectric properties of liquid ionic systems
,” in
Recent Advances in Broadband Dielectric Spectroscopy
,
NATO Science Peace Series B
, edited by
Y. P.
Kalmykov
(
Springer
,
Dordrecht
,
2013
), Chap. 8, pp.
103
122
.
13.
R.
Kubo
,
J. Phys. Soc. Jpn.
12
,
570
(
1957
).
14.
H. J. C.
Berendsen
,
J. R.
Grigera
, and
T. P.
Straatsma
,
J. Phys. Chem.
91
,
6269
(
1987
).
15.
W. R. P.
Scott
,
P. H.
Hünenberger
,
I. G.
Tironi
,
A. E.
Mark
,
S. R.
Billeter
,
J.
Fennen
,
A. E.
Torda
,
T.
Huber
,
P.
Krüger
, and
W. F.
van Gunsteren
,
J. Phys. Chem. A
103
,
3596
(
1999
).
16.
S.
Weerasinghe
and
P. E.
Smith
,
J. Chem. Phys.
119
,
11342
(
2003
).
17.
S.
Nosé
,
Mol. Phys.
52
,
255
(
1984
).
18.
W. G.
Hoover
,
Phys. Rev. A
31
,
1695
(
1985
).
19.
M.
Parrinello
and
A.
Rahman
,
J. Appl. Phys.
52
,
7182
(
1981
).
20.
U.
Essmann
,
L.
Perera
,
M. L.
Berkowitz
,
T.
Darden
,
H.
Lee
, and
L.
Pedersen
,
J. Chem. Phys.
103
,
8577
(
1995
).
21.
R. P.
Futrelle
,
Chem. Phys. Lett.
12
,
285
(
1971
).
22.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
, 1st ed. (
Clarendon Press
,
Oxford
,
1987
).
23.
W. H.
Press
,
S. A.
Teukolsky
,
W. T.
Vetterling
, and
B. P.
Flannery
,
Numerical Recipes: The Art of Scientific Computing
, 2nd ed. (
Cambridge University Press
,
Cambridge
,
1992
).
24.
R.
Buchner
,
G. T.
Hefter
, and
P. M.
May
,
J. Phys. Chem. A
103
,
1
(
1999
).
25.
J. B.
Hubbard
,
P.
Colonomos
, and
P. G.
Wolynes
,
J. Chem. Phys.
71
,
2652
(
1979
).
26.
E.
Levy
,
A.
Puzenko
,
U.
Kaatze
,
P. B.
Ishai
, and
Y.
Feldman
,
J. Chem. Phys.
136
,
114503
(
2012
).
27.
B.
Hess
,
C.
Holm
, and
N.
van der Vegt
,
J. Chem. Phys.
124
,
164509
(
2006
).
28.
R.
Buchner
,
Pure Appl. Chem.
80
,
1239
(
2008
).
29.
J. B.
Hubbard
,
J. Chem. Phys.
68
,
1649
(
1978
).