A non-Bornian analysis of the Gibbs energy of hydration for organic ions

Abstract

Recently, a non-Bornian model was successfully applied to evaluate the Gibbs energy of hydration (ΔG^°_hyd) for spherical ions (mainly inorganic ions). In this model, the long-range, Born-type electrostatic ion–solvent interaction is not explicitly included in the calculation of ΔG^°_hyd, since its contribution is small, whereas the short-range interaction, including Coulomb, polarization, and charge-transfer interactions, is considered as the dominant factor that determines the ΔG^°_hyd of ions. The ΔG^°_hyd scaled by the surface area of an ion can be given by a quadratic function of the surface field strength (E) of the ion. In this study, the non-Bornian model was further applied to organic ions with charged groups. Using the Gaussian 09 program package, the geometries of ions in vacuum were optimized at the B3LYP/6-311++G(2d,p) level, and the partial atomic charges were computed in the Mulliken, Merz-Kollman (MK), natural population analysis (NPA), Hirshfeld, and ChelpG methods. Introducing a new subprogram, we could estimate local electric fields on the ion surface (van der Waals surface or solvent-accessible surface (SAS)). This enabled us to perform regression analyses based on the non-Bornian model, by using the experimental values of ΔG^°_hyd for 109 ions. When the NPA-SAS combination was chosen, the best regression result was obtained, giving the mean absolute error of 4.3 kcal mol⁻¹. The non-Bornian model would provide a simple and relatively accurate way of determining ΔG^°_hyd of ions.

---

1. Introduction

The Gibbs energy of ion hydration (ΔG^°_hyd), which is defined as the energy required for the transfer of an ion from vacuum to water, has received broad attention due to its close relationship to various chemical processes of ions in aqueous solution (e.g., acid–base reactions, redox reactions, ion transport, enzyme reactions, protein folding, etc.). As early as 1920, Born¹ proposed a continuum electrostatic model for the evaluation of ΔG^°_hyd‡² for spherical ions. Since the Born model overestimated the values of ΔG^°_hyd systematically, Abraham and Liszi^3–6 and Kornyshev and Volkov⁷ proposed some modified models, taking into account the dielectric saturation, i.e. a lowering of the permittivity of solvent adjacent to an ion due to its high electric field. In these modified Born models,^3–7 the relative permittivity of solvent (ε) is usually assumed to be very low (∼2) for the primary solvent layer (cf. ε = ∼80 for bulk water). This suggests that the contribution from the primary solvent layer should be more significant than that from the outer bulk solvent. The continuum electrostatic model has been developed into the generalized Born model,^8,9 which treats a solute molecule as a discrete set of overlapping charged spheres imbedded in a polarizable dielectric continuum.

In computational methods such as molecular dynamics (MD) simulation^10,11 and Monte Carlo (MC) simulation,^10,12 the solvent is treated explicitly, i.e. taking into account the detailed structure of solvent molecules. Using these simulations combined with a thermodynamic integration technique^13,14 or a free energy perturbation technique,¹⁵ ΔG^°_hyd can be calculated efficiently. On the other hand, a statistical mechanical integral equation method called the reference interaction site model (RISM)^16,17 has also been utilized for calculating ΔG^°_hyd of ions.^18–20

Previously, our group^21,22 proposed a non-Bornian model for evaluation of the ion-transfer or resolvation energy of ions between two immiscible liquids (i.e., water and nitrobenzene). In the model, the short-range interactions between an ion and its adjacent solvent molecules were taken into account, while the long-range, Born-type electrostatic interaction for the outer bulk solvents was not explicitly included in the calculation. In the recent paper,²³ we have successfully employed the non-Bornian model to evaluate ΔG^°_hyd for spherical ions, including 85 cations (with the charge number, z = 1 to 4) and 53 anions (with z = −1 and −2). In the proposed model,²³ the ΔG^°_hyd scaled by the solvent-accessible surface^24,25 (SAS) area of an ion is simply given by a quadratic function of the electric field strength (E) at the SAS. Three coefficients of the quadratic function of E should be associated with the short-range ion–solvent interaction, which includes Coulomb, polarization, and charge-transfer interactions. Although the coefficients could not accurately be determined in a theoretical manner, they were estimated by multivariate regression analyses for the above spherical cations or anions. The obtained regression equations were found to be useful for accurate prediction of ΔG^°_hyd. The performance of the non-Bornian model was compared with those of other existing models including Born,¹ non-local electrostatic,⁷ and Poisson–Langevin²⁶ models, showing that the non-Bornian model is one of the highest performance models for the prediction of ΔG^°_hyd for spherical ions.

So far, theoretical studies of ΔG^°_hyd have been mainly confined to those for the ions that are spherical or may be assumed as such (i.e., mainly for inorganic ions). Recently, however, experimental values of ΔG^°_hyd for organic ions have been evaluated by using thermochemical cycles that involve pK_a values, gas-phase acidities, and ΔG^°_hyd values of the neutral species and hydrogen ion.^27–29 Then, theoretical models for prediction of ΔG^°_hyd for organic ions have been proposed by previous authors.^28–32

In this study, we have extended the non-Bornian model to organic ions with charged groups. In the extended model, a quantum mechanical calculation by Gaussian 09 (ref. 33 and 34) is used to evaluate partial atomic charges located at nuclear positions in an organic ion. Using a newly developed subprogram, local E values have been obtained at arbitrary points on the molecular surface (or SAS) of the ion. Since the E value is regarded as constant at each local point, the non-Bornian approach can be applied to the respective minute surfaces of the organic ion. The developed novel and simple algorithm has been found to be promising for relatively accurate evaluation of ΔG^°_hyd for organic ions.

2. Theoretical model

The non-Bornian model for the prediction of ΔG^°_hyd for spherical ions has been described in the previous paper.²³ In the model as well as the Born model, it is postulated that a spherical ion has a uniform electric field strength (E) on the ion surface. In this study, we extend the model to organic ions with charged groups, which have non-uniform distribution of E on their surfaces. The proposed model is shown in Fig. 1. In the following, we will present a general formulation for ΔG^°_hyd, which can be applied to either inorganic (spherical) or non-spherical organic ions.

Fig. 1 The proposed hydration model for an organic ion. The dashed line shows the SAS. For details, see the text.

In the non-Bornian approach, the ΔG^°_hyd of an ion is regarded as the energy required for the formation of a nano-sized ion–solvent interface around the ion. The contribution from the long-range electrostatic interaction of an ion with water molecules in the outer solvation shell is not explicitly included in the theoretical framework. By generalizing the Uhlig formula,³⁵ the ΔG^°_hyd of an ion is here expressed as

where S₀ is a closed surface enclosing the ion (the positioning of the surface will be discussed below), dS is an infinitesimal area on S₀, and γ₁₂ is the surface tension of the solvent (2) on the charged ion (1), which is dependent on the local position of the infinitesimal surface. Then, γ₁₂ is related to the (hypothetical) surface tension of the ion (γ₁) in vacuum and that of the solvent (γ₂) with its own saturated vapor, by the Dupré equation:³⁶

γ₁₂ = γ₁ + γ₂ − W₁₂ (2)

where W₁₂ is the work of adhesion, which has the opposite sign of the energy of short-range ion–solvent interaction. For simplicity, it is here assumed that the short-range energy is given by multiplying the energy of the ion–solvent 1 : 1 interaction (U_SR) by the number (N) of primary solvent molecules in the unit surface area of an ion, though this may not be fully valid for small ions. Then, eqn (2) is given by

γ₁₂ = γ₁ + γ₂ + NU_SR (3)

Substituting this equation into eqn (1) yields

Referring to the energy decomposition analysis in ab initio self-consistent field molecular orbital studies,^37–39 U_SR can be partitioned into Coulomb (COU), polarization (POL), charge-transfer (CT), and exchange (EX) terms:

U_SR = U_COU + U_POL + U_CT + U_EX (5)

As described in our previous papers,^21–23 COU, POL, and CT terms can be given by the function of E, i.e., U_COU = −μ < cos θ > E, U_POL = −(1/2)αE², and U_CT = −ζ₀ − ζ₁E − ζ₂E², respectively. Substituting these functions into eqn (5), we obtain

U_SR = −ζ₀ + U_EX − μ < cos θ > E − ζ₁E − (1/2)αE² − ζ₂E² (6)

where μ and α are the dipole moment and electronic polarizability of the solvent molecule, respectively, θ is the angle between the dipole axis and the line connecting the point dipole and point charge, < > indicates the ensemble average, and the CT coefficients of ζ₀, ζ₁, and ζ₂ are influenced by various molecular properties of the solvent and ion, including their electron-donating or -accepting abilities.

Further, substituting eqn (6) into eqn (4), we obtain

with

A = Nζ₀ − NU_EX − γ₁ − γ₂ (8)

B = Nμ < cos θ > + Nζ₁ (9)

C = (1/2)Nα + Nζ₂ (10)

As seen in eqn (9) and (10), coefficient B is due to the COU and CT interactions, while coefficient C is due to the POL and CT interactions.

In order to calculate ΔG^°_hyd based on eqn (7), the closed surface S₀ has been divided into a finite number of minute surface elements (whose area is shown by S_i with i being a finite positive integer). Then it has been assumed that the electric field strength at an arbitrary minute surface has a specific value of E_i (defined as the vertical component). If we also assume that coefficients A, B, and C are independent of either S_i or E_i, the value of ΔG^°_hyd in eqn (7) can be obtained as

ΔG^°_hyd = ∑S_i(−A − BE_i − CE_i²) = −A∑S_i − B∑S_iE_i − C∑S_iE_i² (11)

As shown in eqn (8)–(10), coefficients A, B, and C are defined in terms of several parameters for short-range ion–solvent interactions. However, it is rather difficult to estimate their theoretical values in an accurate manner. This is because the present non-Bornian model is based on some daring assumptions, e.g., neglect of the interaction between the first and outer solvation layers and of the lateral dipole–dipole interactions among the primary solvent molecules. It is also not easy to estimate the CT-interaction parameters, ζ₀, ζ₁, and ζ₂. Consequently, coefficients A, B, and C have been determined by a regression analysis with the literature data of ΔG^°_hyd for organic ions.²⁸

3. Calculation methods

3.1. Minute surface area and electric field strength

In this study, we carried out a Gaussian quantum chemical calculation^33,34 (for details, vide infra) and then obtained atomic coordinate data and partial atomic charges for the optimized structures of organic ions in vacuum. A Microsoft Visual Basic 6.0 program was written for calculating S_i and E_i of each minute surface of an ion (the source code is provided in the ESI†). This program can load a Gaussian logfile that includes atomic coordinate data and partial atomic charges of the ion fully optimized (when optimization is insufficient, further Gaussian calculation is requested). In a polar coordinate system, the program creates the van der Waals (vdW) surface of an ion molecule by using vdW radii (a_vdW) of atoms. The authorized values of a_vdW have been adopted from the literatures,^40,41 being shown in Table 1 as well as in the source code of the program. Furthermore, if necessary, the program provides the corresponding SAS being defined as the locus of the center of the probe sphere (water molecule) as it rolls along the vdW surface; accordingly, the SAS is apart from the vdW surface by Δa (= 0.14 nm; the radius of water molecule; see Fig. 2b). In the first step of the program, a spherical surface is created around each atomic nucleus, for the formation of a molecular surface (i.e., the vdW surface or SAS). The spherical surface is divided into 10386 of minute surfaces of an equal area (S_i). The spherical surface formed around an atom intersects those around adjacent atom(s). Such portions of the spherical surface as included in the region of adjacent atom(s) are omitted, and the remaining part of the spherical surface is used to create the molecular surface.

Table 1 The van der Waals radii (a_vdW) of atoms

Atom	avdW (Å)	Ref.	Atom	avdW (Å)	Ref.
H	1.20	40	Rb	3.04	41
He	1.40	40	Sr	2.63	41
Li	1.82	40	Y	2.20	41
Be	1.38	41	Zr	1.96	41
B	1.20	41	Nb	1.86	41
C	1.70	40	Mo	1.76	41
N	1.55	40	Tc	1.73	41
O	1.52	40	Ru	1.81	41
F	1.47	40	Rh	1.75	41
Ne	1.54	40	Pd	1.63	40
Na	2.27	40	Ag	1.77	40
Mg	1.73	40	Cd	1.58	40
Al	1.75	41	In	1.93	40
Si	2.10	40	Sn	2.17	40
P	1.80	40	Sb	2.03	41
S	1.80	40	Te	2.06	40
Cl	1.75	40	I	1.98	40
Ar	1.88	40	Xe	2.16	40
K	2.75	40	Cs	3.27	41
Ca	2.41	41	Ba	2.71	41
Sc	1.98	41	La	2.29	41
Ti	1.80	41	Hf	1.94	41
V	1.72	41	Ta	1.87	41
Cr	1.67	41	W	1.77	41
Mn	1.66	41	Re	1.75	41
Fe	1.65	41	Os	1.83	41
Co	1.64	41	Ir	1.77	41
Ni	1.63	40	Pt	1.72	40
Cu	1.40	40	Au	1.86	41
Zn	1.39	40	Hg	1.55	40
Ga	1.87	40	Tl	1.96	40
Ge	1.77	41	Pb	2.02	40
As	1.85	40	Bi	2.17	41
Se	1.90	40	Th	2.20	41
Br	1.85	40	U	1.86	41
Kr	2.02	40

---

Fig. 2 Calculation of the electric field (E_i) at a minute surface of (a) the van der Waals surface and (b) the SAS for a two-atom ion. Outlined white arrows show the electric fields created from the respective point charges (q₁ and q₂), which give the electric field (black arrow) at the minute surface as the vector sum. For further details, see the text.

Because electric fields satisfy the superposition principle, the electric field strength at any point on a molecular surface can be obtained from the vector sum of the electric fields created from the point charges of all constituent atoms. This is exemplified in the scheme of Fig. 2, in which the ion is simplified as a two-atom ion. A contribution to the electric field at a minute surface from each point charge (q_j, where j indicates a specific atom) can be given by Gauss' law: E_ij = q_j/(4πε₀R_ij²) where ε₀ and R_ij are, respectively, the permittivity of vacuum and the distance from atom j to minute surface i (i.e., R_ij = a_vdW or a_vdW + Δa). In the present program, the electric field (E_i) on a minute surface is provided as the component vertical to the molecular surface (see again Fig. 2). An execution of the program provides two files in txt and RES formats. The former file describes Cartesian coordinates in three dimensions for all points indicating the minute surfaces and the corresponding S_i and E_i values, while the latter one describes the calculation of numerical sums (i.e., ∑S_i, ∑S_iE_i, and ∑S_iE_i²) appearing in eqn (11). Using these calculation values and the literature values of ΔG^°_hyd, we have performed regression analysis to obtain coefficients A, B, and C in eqn (11).

3.2. Quantum chemical calculation

The Gaussian 09 program package^33,34 with the B3LYP hybrid density functional theory^42–45 was used to optimize the ground-state geometry of organic ions in vacuum. Recently, a variety of quantum chemical studies have been performed on the solvent effect on ion solvation (e.g., see ref. 28 and 30–32), however we carried out quantum chemical calculations for ions isolated in vacuum, in view of figuring out ΔG^°_hyd from the intrinsic properties of ions (i.e., their size, shape, charge distribution, etc.). At the B3LYP/6-311++G(2d,p) level, five different quantum chemical charges were computed, including Mulliken,⁴⁶ Merz–Kollman (MK),^47,48 natural population analysis (NPA),⁴⁹ Hirshfeld,⁵⁰ and ChelpG⁵¹ charges.

4. Results and discussion

4.1. Regression analysis for 109 ions

In this study, a non-Bornian analysis was carried out for 109 ions in total, which include 52 cations and 57 anions. They are listed in Table 2, being mainly organic ions. The ΔG^°_hyd values (in kcal mol⁻¹) for all the ions were reported by Kelly et al.,²⁸ which were determined by using ΔG^°_hyd(H⁺) = −265.9 kcal mol⁻¹ (= −1113 kJ mol⁻¹).⁵²

Table 2 Estimation of the Gibbs hydration energies at 298 K for 52 cations and 57 anions using the non-Bornian model

Cation	ΔG^°hyd (kcal mol^−1)			Anion	ΔG^°hyd (kcal mol^−1)
	Exp.^a	Non-Bornian^b	Diff.^c		Exp.^a	Non-Bornian^b	Diff.^c
a Experimental data taken from ref. 28.b Theoretical values obtained from eqn (13) with the NPA charges and the E values at the SAS.c Differences between the experimental and theoretical values.
H3O^+	−110.3	−105.8	−4.5	OH^−	−104.7	−107.1	2.4
CH3OH2^+	−93.0	−83.1	−9.9	HOO^−	−97.3	−91.7	−5.6
CH3CH2OH2^+	−88.4	−73.5	−14.9	O2^−	−83.3	−89.1	5.8
(CH3)2OH^+	−79.7	−71.1	−8.6	HS^−	−72.1	−88.2	16.1
(C2H5)2OH^+	−71.5	−61.8	−9.7	HC2^−	−76.5	−81.6	5.1
(CH3)2C=OH^+	−77.1	−67.5	−9.6	CN^−	−70.2	−87.7	17.5
CH3C(OH)C6H5^+	−64.5	−59.1	−5.4	CH3O^−	−95.0	−88.4	−6.6
NH4^+	−85.2	−97.3	12.1	CH3CH2O^−	−90.7	−80.2	−10.5
CH3NH3^+	−76.4	−80.4	4.0	CH3CH2CH2O^−	−88.3	−80.4	−7.9
CH3(CH2)2NH3^+	−71.5	−72.7	1.2	(CH3)2CHO^−	−86.3	−76.1	−10.2
(CH3)2CHNH3^+	−69.6	−69.1	−0.5	CH3CH2CHOCH3^−	−84.2	−75.7	−8.5
(CH3)3CNH3^+	−67.3	−65.6	−1.7	(CH3)3CO^−	−82.3	−74.1	−8.2
c-C6H11NH3^+	−68.7	−65.7	−3.0	H2C=CHCH2O^−	−86.6	−78.9	−7.7
H2C=CHCH2NH3^+	−72.0	−71.9	−0.1	CH3OCH2CH2O^−	−89.4	−81.7	−7.7
(CH3)2NH2^+	−68.6	−69.9	1.3	HOCH2CH2O^−	−85.3	−79.7	−5.6
(C2H5)2NH2^+	−63.4	−60.4	−3.0	C6H5CH2O^−	−85.1	−77.1	−8.0
(n-C3H7)2NH2^+	−60.5	−60.1	−0.4	CF3CH2O^−	−77.5	−71.4	−6.1
(H2C=CHCH2)2NH2^+	−61.6	−59.6	−2.0	CH(CF3)2O^−	−65.5	−63.0	−2.5
(CH3)3NH^+	−61.1	−64.0	2.9	CH3OO^−	−93.2	−86.1	−7.1
(C2H5)3NH^+	−54.6	−56.0	1.4	CH3CH2OO^−	−89.2	−85.2	−4.0
(n-C3H7)3NH^+	−50.9	−55.9	5.0	HCO2^−	−76.2	−79.1	2.9
C6H5NH3^+	−72.4	−70.7	−1.7	CH3CO2^−	−77.6	−77.6	0.0
o-CH3C6H4NH3^+	−70.3	−68.0	−2.3	CH3CH2CO2^−	−76.2	−76.3	0.1
m-CH3C6H4NH3^+	−69.6	−70.3	0.7	CH3(CH2)4CO2^−	−74.6	−81.2	6.6
p-CH3C6H4NH3^+	−69.8	−70.6	0.8	H2C=CHCO2^−	−74.0	−74.9	0.9
m-NH2C6H4NH3^+	−65.8	−69.8	4.0	CH3COCO2^−	−68.5	−72.3	3.8
C6H5NH2CH3^+	−62.6	−63.6	1.0	CH2ClCO2^−	−69.7	−72.4	2.7
C6H5NH2CH2CH3^+	−62.2	−59.7	−2.5	CHCl2CO2^−	−62.3	−68.5	6.2
C6H5NH(CH3)2^+	−57.2	−60.0	2.8	CF3CO2^−	−59.3	−64.4	5.1
p-CH3C6H4NH(CH3)2^+	−55.9	−60.5	4.6	C6H5CO2^−	−71.2	−73.1	1.9
C6H5NH(CH2CH3)2^+	−54.0	−56.8	2.8	C6H5O^−	−71.9	−69.4	−2.5
C10H7NH3^+	−67.4	−70.3	2.9	o-CH3C6H4O^−	−70.2	−67.4	−2.8
C2H4NH2^+	−70.9	−73.6	2.7	m-CH3C6H4O^−	−71.1	−70.2	−0.9
C3H6NH2^+	−67.7	−68.1	0.4	p-CH3C6H4O^−	−72.0	−69.9	−2.1
C4H8NH2^+	−66.0	−65.3	−0.7	m-HOC6H4O^−	−73.8	−68.3	−5.5
C5H10NH2^+	−64.2	−63.1	−1.1	p-HOC6H4O^−	−77.6	−69.1	−8.5
C6H12NH2^+	−63.3	−61.6	−1.7	o-NO2C6H4O^−	−60.1	−66.3	6.2
C4H4NH2^+	−61.4	−68.1	6.7	m-NO2C6H4O^−	−61.9	−63.8	1.9
PyridineH^+	−61.1	−64.8	3.7	p-NO2C6H4O^−	−57.8	−63.1	5.3
C9H7NH^+	−56.0	−60.0	4.0	o-ClC6H4O^−	−66.1	−66.5	0.4
C4H8NHNH2^+	−66.0	−63.4	−2.6	p-ClC6H4O^−	−66.0	−66.7	0.7
CH3CNH^+	−75.3	−75.1	−0.2	CH2CHO^−	−76.5	−76.4	−0.1
H2NNH3^+	−84.6	−85.9	1.3	CH3C(O)CH2^−	−76.2	−73.6	−2.6
p-CH3OC6H4NH3^+	−71.2	−74.4	3.2	CH3CH2C(O)CHCH3^−	−73.7	−67.9	−5.8
p-NO2C6H4NH3^+	−75.9	−83.9	8.0	NCNH^−	−72.2	−78.5	6.3
C4H8ONH2^+	−69.6	−68.6	−1.0	CH2CN^−	−66.6	−74.2	7.6
CH3CONH3^+	−73.9	−77.6	3.7	C6H5NH^−	−62.9	−66.1	3.2
C6H5CONH3^+	−67.2	−70.4	3.2	p-NO2C6H5NH^−	−57.4	−62.3	4.9
(CH3)2SH^+	−64.5	−68.3	3.8	(C6H5)2N^−	−54.6	−62.5	7.9
(CH3)2SOH^+	−67.7	−70.7	3.0	CH3CONH^−	−80.2	−78.5	−1.7
m-ClC6H4NH3^+	−74.7	−72.1	−2.6	CH2NO2^−	−76.5	−75.9	−0.6
p-ClC6H4NH3^+	−74.1	−72.3	−1.8	CH3S^−	−73.8	−76.0	2.2
				CH3CH2S^−	−71.8	−70.9	−0.9
				CH3CH2CH2S^−	−70.5	−71.1	0.6
				C6H5S^−	−63.4	−65.3	1.9
				CH3S(O)CH2^−	−67.7	−74.5	6.8
				CCl3^−	−54.1	−61.6	7.5

---

In the above, we describe that coefficients A, B, and C in eqn (8)–(10) are assumed to be independent of E_i, i.e., the magnitude of the electric field strength at minute ion surfaces. However, these coefficients would be different between cations and anions, because the orientation of primary water molecules on the respective ion surfaces may be dependent largely on the sign of E_i (i.e., positive/negative). Fig. 3 exemplifies the distributions of E_i for (a) 4-nitrophenoxide anion and (b) 4-methoxyanilinium cation, which have been calculated with the NPA–SAS combination (vide infra). As seen in (a), the anion shows negative E_i values at all surface positions. On the other hand, as seen in (b), the cation shows positive E_i values at most surface positions, but shows slightly negative values (as indicated by pale blue) at around the oxygen atom of the methoxy group. In general, a cation can partially have anionic surface, and vice versa for an anion. Anyway, since the orientation of water molecules on positive and negative surfaces may be different, we used two different values for coefficients B and C, depending on the sign of E_i, and then rewrote eqn (11) as

where B₊ and C₊ are the coefficients, B and C, for positive E_i values, and where B₋ and C₋ are those for negative E_i values. This equation enables us to analyze both cations and anions simultaneously.

Fig. 3 Distributions of E_i (in V nm⁻¹) at the SAS for (a) 4-nitrophenoxide anion and (b) 4-methoxyanilinium cation. In this calculation, the partial atomic charges have been obtained by the NPA method.

Multivariate regression analyses were performed on the experimental values²⁸ of ΔG^°_hyd to determine the coefficients in eqn (12) (i.e., A, B₊, C₊, B₋, and C₋). The values of ∑S_i, , , , and required for the regression analyses were obtained from the quantum chemical calculation described above. In the present study, we investigated the performance of regression analysis by using all combinations of five different partial atomic charges (Mulliken, MK, NPA, Hirshfeld, and ChelpG) and two different molecular surfaces (vdW surface and SAS). As shown in Fig. 4, we obtained the highest performance when using a combination of NPA charge and SAS. The regression equation obtained is

where ΔG^°_hyd is in kcal mol⁻¹, S_i in nm², and E_i in V nm⁻¹. In this analysis with the NPA–SAS combination, the square of the correlation coefficient (R²) is 0.9943; the mean absolute error (, where n is the number of data; x_j is the difference between the experimental and theoretical values) and the root mean square error () for ΔG^°_hyd are 4.3 kcal mol⁻¹ and 5.5 kcal mol⁻¹, respectively. Fig. 5 shows the comparison of the experimental and theoretical values of ΔG^°_hyd. The theoretical as well as experimental values are shown in Table 2. Thus, eqn (13) has been found to be promising for relatively accurate evaluation of ΔG^°_hyd for organic ions other than those analyzed here.

Fig. 4 Comparison of the performances of non-Bornian analyses made by using all combinations of five different partial atomic charges (Mulliken, MK, NPA, Hirshfeld, and ChelpG) and two different molecular surfaces (vdW surface and SAS). For each combination, the values of MAE and RMSE (left y-axis, bars) and R² (right y-axis, circles) are shown.

Fig. 5 Scatter plots of the theoretical values of ΔG^°_hyd obtained using the non-Bornian model (i.e., eqn (13)) against the experimental values of ΔG^°_hyd for 109 ions (shown in Table 2). The combination of NPA–SAS has been employed. R² = 0.9943.

As shown in Fig. 4, generally higher regression performance was achieved in the use of SAS, rather than in the use of vdW surface, particularly when the Mulliken and NPA methods were used to evaluate partial atomic charges. This is in harmony with the previous regression results obtained for spherical ions.²³ It is generally considered that the electrostatic interaction, including COU and POL interactions, makes a dominant contribution to ΔG^°_hyd; the SAS seems to be more suitable for evaluating such electrostatic interactions.

There is no “correct” method for assigning partial atomic charges, because they are not quantum mechanical observable.⁵³ However, the NPA method (with SAS) has achieved the highest performance among the methods tested. This method generally overestimates ionic character, giving larger magnitude of atomic charges than the Mulliken method.⁵⁴ This may give the NPA method an advantage for evaluating the electrostatic ion–solvent interaction as the main contribution to ΔG^°_hyd. As also seen in Fig. 4, the Mulliken method gave no good results for the present regression analysis. The Mulliken method, which is the default method in Gaussian 09,^33,34 has been widely used for determining partial atomic charges, however its weakness is well known. The method tends to yield unnatural values at large basis sets,⁵⁵ such as B3LYP/6-311++G(2d,p) employed in this study.

4.2. Comparison with other previous models

In the present non-Bornian analysis, the lowest value of MAE (alias mean unsigned error; MUE) is 4.3 kcal mol⁻¹. The data set used for the regression analysis is similar to that used by Cramer and Truhlar's (C–T) group,²⁸ albeit our data set does not include water-clustered ions. The C–T group used their SM6 continuum solvation model (employing a generalized Born model) to analyze the data set including 112 unclustered ions, plus 31 clustered ions. The MAE value obtained when using the most suitable basis set(s) is 4.5 kcal mol⁻¹, which is equivalent to that obtained in this study. The C–T group, however, refined their model on several occasions to present the newest model (called SM12) giving a lower MAE value of 2.9 kcal mol⁻¹ for 112 unclustered and clustered ions.³² Besides, Curutchet et al.³⁰ employed the optimized Miertus–Scrocco–Tomasi (MST) continuum model to obtain the MAE value of 3.6 kcal mol⁻¹ (calculated from Table 2 in ref. 30) for a data set of 47 ions. Lee et al.²⁹ used the modified solvation free energy density (SFED) model to obtain a rather low MAE value of 1.7 kcal mol⁻¹ for a data set of 90 ions. Da Silva et al.³¹ proposed a method explicitly representing the solvation shell in continuum solvent calculations, and then obtained the MAE value of 2.1 and 2.8 kcal mol⁻¹ for 30 cations and 30 anions, respectively.

Thus, the MAE values reported by the previous authors are comparable with or somewhat lower than that obtained in this study. However, no rigorous comparison can be made, because the data sets employed in the previous studies are different with each other, though there are many ions in common. Nevertheless, we believe the performance of our non-Bornian model is not so bad, probably comparable with other previous models. In most of the previous models, a number of parameters have been adjusted so as to optimize the agreement with experimental values; some sort of parameters, such as “atomic surface tension”, have been optimized for each atom or atomic group. In our non-Bornian model, however, there are only five adjusting parameters (i.e., A, B₊, C₊, B₋, and C₋ in eqn (12)), which are common to all atoms or atomic groups. It should be stressed that this very simple model can predict ΔG^°_hyd for various organic ions with fair accuracy. Using a higher-level quantum chemical calculation, the performance of the non-Bornian model would be further improved. In practice, a higher-level calculation at the M06-2X/6-311++G(2d,p) level slightly improved the result shown in eqn (13) (with MAE = 4.1 kcal mol⁻¹; RMSE = 5.3 kcal mol⁻¹; detailed data not shown), though the calculation time was 1.5 to 3 times longer than the above-shown calculation at the B3LYP/6-311++G(2d,p) level. The use of MP2 and CCSD, which were reported as more reliable than B3LYP,⁵⁶ might provide a better or worse result than the present analysis.

We would like to add that in recent years, some other models with few adjusting parameters have been proposed for the evaluation of hydration energies of monovalent spherical ions and noble gas atoms⁵⁷ and neutral and ionic solutes.⁵⁸ Though no rigorous comparison can be made because of the difference in the solutes to be tested, a comparatively low value of MUE (= 2.4 kcal mol⁻¹) has been reported for ionic solutes.⁵⁸

5. Conclusions

It has been found that the non-Bornian model is also available for the evaluation of ΔG^°_hyd for organic ions having charged groups. Using the Gaussian 09 program package with a newly developed subprogram, the local electric field strength (E_i) on the surface of an organic ion can be estimated from partial atomic charges of the ion. Using the simple, semi-empirical equation (eqn (12)) with the estimated E_i values allows us to perform a satisfactory regression analysis for the experimental data of ΔG^°_hyd for organic ions. When the NPA method is used for the evaluation of partial atomic charges and the E_i values at the SAS of an ion are employed, the most excellent match between the experimental and theoretical values is achieved over 109 ions. The MAE value is 4.3 kcal mol⁻¹, being comparable with or somewhat higher than those reported by the previous authors.^28–32 However, considering the smaller number (i.e., five) of adjusting parameters and their independence from atoms or groups, the non-Bornian model would be useful as a simple tool to estimate ΔG^°_hyd of ions in a sufficiently accurate manner.

References

1. M. Born, Z. Phys., 1920, 1, 45–48.
2. N. Bazhin, ISRN Thermodyn., 2012, 204104.
3. M. H. Abraham and J. Liszi, J. Chem. Soc., Faraday Trans. 1, 1978, 74, 1604–1614.
4. M. H. Abraham and J. Liszi, J. Chem. Soc., Faraday Trans. 1, 1978, 74, 2858–2867.
5. M. H. Abraham, J. Liszi and L. Mészáros, J. Chem. Phys., 1979, 70, 2491–2496.
6. M. H. Abraham and J. Liszi, J. Inorg. Nucl. Chem., 1981, 43, 143–151.
7. A. A. Kornyshev and A. G. Volkov, J. Electroanal. Chem., 1984, 180, 363–381.
8. R. Constanciel and R. Contreras, Theor. Chim. Acta, 1984, 65, 1–11.
9. W. Still, A. Tempczyk, R. Hawley and T. Hendrickson, J. Am. Chem. Soc., 1990, 112, 6127–6129.
10. T. P. Lybrand, I. Ghosh and J. A. McCammon, J. Am. Chem. Soc., 1985, 107, 7793–7794.
11. M. Soniat and S. W. Rick, J. Chem. Phys., 2012, 137, 044511.
12. M. Migliore, G. Corongiu, E. Clementi and G. C. Lie, J. Chem. Phys., 1988, 88, 7766–7771.
13. T. P. Straatsma and H. J. C. Berendsen, J. Chem. Phys., 1988, 89, 5876–5886.
14. K. Leung, S. B. Rempe and O. A. von. Lilienfeld, J. Chem. Phys., 2009, 130, 204507.
15. J. Carisson and J. Åqvist, J. Phys. Chem. B, 2009, 113, 10255–10260.
16. D. Chandler and H. C. Andersen, J. Chem. Phys., 1972, 57, 1930–1937.
17. F. Hirata and P. Rossky, Chem. Phys. Lett., 1981, 83, 329–334.
18. B. Roux, H. A. Yu and M. Karplus, J. Phys. Chem., 1990, 94, 4683–4688.
19. S.-H. Chong and F. Hirata, J. Phys. Chem. B, 1997, 101, 3209–3220.
20. M. V. Fedotova and S. E. Kruchinin, Russ. Chem. Bull., 2011, 60, 223–228.
21. T. Osakai and K. Ebina, J. Phys. Chem. B, 1998, 102, 5691–5698 , erratum sheets are available from http://www2.kobe-u.ac.jp/%7Eosakai/OldPapers.html.
22. W. Murakami, K. Eda, M. Yamamoto and T. Osakai, J. Electroanal. Chem., 2013, 704, 38–43.
23. W. Murakami, K. Eda, M. Yamamoto and T. Osakai, Bull. Chem. Soc. Jpn., 2014, 87, 403–411.
24. B. Lee and F. M. Richards, J. Mol. Biol., 1971, 55, 379–400.
25. A. Shrake and J. A. Rupley, J. Mol. Biol., 1973, 79, 351–371.
26. P. Koehl, H. Orland and M. Delarue, J. Phys. Chem. B, 2009, 113, 5694–5697.
27. J. R. Pliego, Jr and J. M. Riveros, Phys. Chem. Chem. Phys., 2002, 4, 1622–1627.
28. C. P. Kelly, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2006, 110, 16066–16081.
29. S. Lee, K.-H. Cho, C. J. Lee, G. E. Kim, C. H. Na, Y. In and K. T. No, J. Chem. Inf. Model., 2011, 51, 105–114.
30. C. Curutchet, A. Bidon-Chanal, I. Soteras, M. Orozco and F. J. Luque, J. Phys. Chem. B, 2005, 109, 3565–3574.
31. E. F. da Silva, H. F. Svendsen and K. M. Merz, J. Phys. Chem. A, 2009, 113, 6404–6409.
32. A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Chem. Theory Comput., 2013, 9, 609–620.
33. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, N. J. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, Ö. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, and D. J. Fox, Gaussian 09, Revision C.01, Gaussian, Inc., Wallingford, CT, 2009.
34. R. Dennington, T. Keith and J. Millam, GaussView, Version 5.0.9, Semichem Inc., Shawnee Mission, KS, 2009.
35. H. H. Uhlig, J. Phys. Chem., 1937, 41, 1215–1225.
36. J. N. Israelachvili, Intermolecular and Surface Forces, Academic Press, London, 3rd edn, 2011, ch. 17.
37. K. Kitaura and K. Morokuma, Int. J. Quantum Chem., 1976, 10, 325–340.
38. H. Umeyama and K. Morokuma, J. Am. Chem. Soc., 1977, 99, 1316–1332.
39. A. Karpfen and P. Schuster, The Chemical Physics of Solvation, Part A, ed. R. R. Dogonadze, E. Kálmán, A. A. Kornyshev and J. Ulstrup, Elsevier, Amsterdam, 1985, ch. 7.
40. A. Bondi, J. Phys. Chem., 1964, 68, 441–451.
41. S. S. Batsanov, Inorg. Mater., 2001, 37, 871–885.
42. A. D. Becke, Phys. Rev. A, 1988, 38, 3098–3100.
43. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
44. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
45. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627.
46. R. S. Mulliken, J. Chem. Phys., 1955, 23, 1833–1840.
47. U. C. Singh and P. A. Kollman, J. Comput. Chem., 1984, 5, 129–145.
48. B. H. Besler, K. M. Merz and P. A. Kollman, J. Comput. Chem., 1990, 11, 431–439.
49. A. E. Reed, L. A. Curtiss and F. Weinhold, Chem. Rev., 1988, 88, 899–926.
50. F. L. Hershel, Theor. Chim. Acta, 1977, 44, 129–138.
51. C. M. Breneman and K. B. Wiberg, J. Comput. Chem., 1990, 11, 361–373.
52. M. D. Tissandier, K. A. Cowen, W. Y. Feng, E. Gundlach, M. H. Cohen, A. D. Earhart, J. V. Coe and T. R. Tuttle, J. Phys. Chem. A, 1998, 102, 7787–7794.
53. J. W. Storer, D. J. Giesen, C. J. Cramer and D. G. Truhlar, J. Comput.-Aided Mol. Des., 1995, 9, 87–110.
54. T. Kar, J. G. Angyan and A. B. Sannigrahi, J. Phys. Chem. A, 2000, 104, 9953–9963.
55. D. C. Young, Computational Chemistry: A Practical Guide for Applying Techniques to Real-World Problems, Wiley, London, 2001, ch. 12.
56. M. T. Nguyen, S. Creve and L. G. Vanquickenborne, J. Phys. Chem., 1996, 100, 18422–18425.
57. T. T. Duignan, D. F. Parsons and B. W. Ninham, J. Phys. Chem. B, 2013, 117, 9421–9429.
58. A. Pomogaeva and D. M. Chipman, J. Chem. Theory Comput., 2014, 10, 211–219.

---

Footnotes

† Electronic supplementary information (ESI) available: The source code (c4ra02422b1.txt) and Readme (c4ra02422b5.txt) of the program for the calculation of E_i and S_i and three Visual Basic files (ProgramNonBorn.vbp; ProgramNonBorn.frm; ProgramNonBorn.vbw; being uploaded as c4ra02422b2.txt, c4ra02422b3.txt, and c4ra02422b4.txt, respectively) for the execution of the program. See DOI: 10.1039/c4ra02422b

‡ It is conventionally accepted that the Born formula describes the Gibbs energy, however it has recently been claimed that the Born formula describes the enthalpy. See ref. 2.

---