- Email: [email protected]
Correction terms for the solvation free energy functional of three-dimensional reference interaction site model based on the reference-modified density functional theory
Journal of Molecular Liquids 291 (2019) 111160
Contents lists available at ScienceDirect
Journal of Molecular Liquids journal homepage: www.elsevier.com/locate/molliq
Correction terms for the solvation free energy functional of three-dimensional reference interaction site model based on the reference-modified density functional theory夽 Yutaka Maruyama Architecture Development Team, FLAGSHIP 2020 Project, RIKEN Center for Computational Science, Kobe 650-0047, Japan
A R T I C L E
I N F O
Article history: Received 28 February 2019 Received in revised form 13 May 2019 Accepted 8 June 2019 Available online 26 June 2019 Keywords: Solvation free energy 3D-RISM theory Classical DFT Water solvent Organic molecules
A B S T R A C T Recently, we derived the reference-modified density functional theory (RMDFT) to calculate the solvation free energy (SFE) with high accuracy, introducing a reference system to the DFT for polyatomic molecular liquids. The RMDFT functional consists of the DFT using second-order approximation (DFT(D)) functional and the correction terms. The DFT(D) functional are equivalent to the Singer-Chandler formula for the hypernetted chain closure equation. The contribution of the correction terms has a strong correlation with the partial molar volume (PMV) of solute. Thus, the RMDFT functional has almost the same accuracy as the PMV correction. Furthermore, we combined the Singer-Chandler formula for the Kovalenko-Hirata closure and the correction terms to avoid inconsistency between the closure equation and the SFE functional. It become clear that the inconsistency is unimportant, whereas the combined functional works well for the amino acid side chain analogues molecules. © 2019 Elsevier B.V. All rights reserved.
1. Introduction The solvation free energy (SFE) is an important property for investigation of the thermodynamic stabilities of systems, in addition to chemical reaction, protein folding, and protein-ligand binding. One way to calculate the SFE is to employ the three-dimensional reference interaction site model (3D-RISM) theory, which is a statistical mechanics theory for molecular liquids [1-3]. To increase the speed of 3D-RISM calculations, a highly efficient multigrid algorithm was developed [4, 5]. In addition, the 3D-RISM program is implemented using graphics processing unit (GPU) [6] and massively parallel supercomputer [7]. It is also implemented in commercial programs such as the Amsterdam Density Functional (ADF) modeling suite [8] and the Amber package [9]. Therefore, the 3D-RISM theory has been applied to various biological systems, such as, for example, ion binding protein [10, 11], staphylococcal nuclease [12], aquaporin [13-15], rosette nanotubes [16-20], amyloidb (Ab42) [21-23], photoactive yellow protein [24], DNA [25, 26], and telomere [27, 28].
夽 The article was accepted for publication under the “VSI: EMLG/JMLG 2018 in Nagoya”. E-mail address: [email protected] (Y. Maruyama).
https://doi.org/10.1016/j.molliq.2019.111160 0167-7322/© 2019 Elsevier B.V. All rights reserved.
The Singer-Chandler formula [29, 30] is commonly used when calculating the SFE using the RISM/3D-RISM theory. However, Ten-no pointed out that the absolute values of the SFE artificially increase with the number of sites of solute molecules, a theoretical drawback of the RISM theory [31, 32]. Therefore, some researchers proposed introducing a phenomenological correction in SFE calculation methods in which the partial molar volume of the solute molecule cancels the artificial solute-size dependence of the SFE [33-35]. Similarly, a cavity correction [36] and a pressure correction are also proposed for the SFE [37-39]. However, which of these correction method is the best remains controversial. Furthermore, the differences between isochoric and isobaric conditions and ensembleindependence of the SFE are discussed in detail [40-42]. Recently, Miyata and Thapa pointed out that one of the problems on the RISM theory should be attributed to the theoretical drawback of the hypernetted-chain (HNC) type approximation itself rather than to the simple extension of the HNC-type closure. They showed that the SFE of a Lennard-Jones (LJ) solute in an LJ solvent depends on solute size by using the Ornstein-Zernike (OZ) integral equation with the HNC and Kovalenko-Hirata (KH) closures. SFE values obtained from the HNC approximation and the KH approximation increase monotonically as the solute size increases, whereas the SFE determined using molecular dynamics (MD) simulations decrease [43]. They actually applied an approximate bridge function adjusting the s parameter in LJ potential of
2
Y. Maruyama / Journal of Molecular Liquids 291 (2019) 111160
the HNC-type closure equation to obtain a more accurate thermodynamic quantity [44]. Thus, not only the RISM theory for molecular liquids but also the OZ theory for simple liquids results in solute-size dependence and significant overestimation of the SFE using the HNC or KH closure. Inspired by their study, we derived a new SFE functional using another approach. We proposed a reference-modified density functional theory (RMDFT) to systematically develop a free-energy density functional model [45-47]. We employ a hardsphere (HS) fluid as the reference system instead of an ideal polyatomic molecular gas to derive the SFE functional of a solute molecule in water. The SFE functional derived from the RMDFT approach improved the overestimation of the SFE originating from the HNC-type approximation. The RMDFT functional improved the absolute values of the SFE for neutral amino acid side chain analog molecules and for 504 small organic molecules [45, 46]. Furthermore, we applied the RMDFT functional to investigate the Lennard-Jones solution [47], the structural stability of protein during folding process [48], and the high predictive ability for temperature and pressure dependences of the SFE [49]. The RMDFT functional is composed of two parts: the secondorder approximation DFT (DFT(D)) functional and the correction terms (CT). The DFT(D) functional is equivalent to the SingerChandler formula for HNC closure (HNC functional) [50, 51]. We actually compared the DFT(D) functional and the HNC functional in the SFE calculation for neutral amino acid side chain analogues molecules (Table 1). Despite the deviation from the experimental values, all values of the DFT functional and the HNC functional agree with an accuracy of 0.1 kcal/mol or less. If the DFT functional contains a higher order approximation, deviation of the DFT functional from the HNC functional will occur. However, the higher order terms do not contribute significantly to improve the calculation value [47]. For closing the 3D-RISM equation, we chose the KH closure equation [3, 52] because of its advantage of rapid convergence, whereas the HNC closure shows an extremely poor convergence. The DFT(D) functional is equivalent to the HNC functional, but not SingerChandler formula for the KH closure (KH functional). In other words, when using RMDFT functional with KH closure equation, there is an inconsistency between the closure equation and the SFE functional because the Singer-Chandler formula has a different form for each closure equation. In this study, we combine the KH functional and the CT based RMDFT for consistency. We calculate and verify the SFE values for
Table 1 Solvation free energy values of amino acid side chain analogues in water. Amino acid
Analog solute
Experimental
DFT(D)
HNC
Leu Ile Val Ala Phe Cys Met Thr Ser Trp Tyr Gln Asn Hid Hie
Isobutane n-Butane Propane Methane Toluene Methanethiol Methyl ethyl sulfide Ethanol Methanol 3-Methylindole p-Cresol Propionamide Acetamide 4-Methylimidazole 4-Methylimidazole
2.3 2.2 2.0 1.9 −0.8 −1.2 −1.5 −4.9 −5.1 −5.9 −6.1 −9.4 −9.7 −10.3 −10.3
17.2 17.2 14.3 8.1 16.5 6.1 13.2 4.6 0.7 17.1 11.9 4.2 0.4 4.0 3.5
17.2 17.2 14.3 8.1 16.5 6.1 13.2 4.6 0.7 17.1 11.9 4.2 0.4 4.0 3.5
Energy unit is kcal/mol. Experimental data were taken from the literature [78] . We used hypernetted-chain closure in the calculations.
amino acid side chain analogues and 504 organic molecules using the KH functional with the CT. We also investigate how the CT of RMDFT functional works.
2. Methods 2.1. RMDFT functional We derived RMDFT by introducing a reference system to the DFT for polyatomic molecular liquids. Here, we consider the SFE functional for water solvent. Then, we incorporate the reference HS only in water oxygen because water hydrogen is buried in the oxygen atom. The RMDFT functional of the SFE for water solvent is expressed as [45, 46] 1 qc drhc (r) b c q ex qc dr dr Ccc + (|r − r |)hc (r) b c c 1 ex qc qc dr dr Ccc + (|r − r |)hc (r)hc (r ) 2b c c dF HS [qO ] (hO (r) + 1) − lOHS , + DF HS [qO ] − qO dr d(qO hO (r))
DlRMDFT = −
(1)
where ex Cab (|r
− r |) =
HS C¯ OO (|r − r |) − COO (|r − r |) (a = b = O) C¯ ab (|r − r |) (otherwise)
IM C¯ ab (|r − r |) = Cab (|r − r |) + Cab (|r − r |) −1 dab sin(kLab ) C˜ ab (k) = − dab q + (1 − dab )q . q kLab
(2)
Here, b = 1/(kB T) is the inverse of the Boltzmann constant times temperature, qc denotes the average number density of the solvent site, c, q is the average number density of the water solvent, hc (r) is the 3D total correlation function of the solvent site c around solute, Cab (r) is the site-site direct correlation function of pure solvent, IM Cab (r) is the intramolecular direct correlation function [53], Lab is the bond length between a and b sites, O denotes oxygen site of water solvent, CHS OO (r) correspond to the site-site direct correlation function of the reference HS fluid, FHS [qO ] is the excess intrinsic free energy functional for HS fluid, and l HS O is the excess chemical potential of the reference HS fluid. We obtain hc (r) by 3D-RISM calculations, Cab (r) by solvent-solvent RISM calculations, and CHS OO (r) by hardsphere fluids theory. The RMDFT functional is decomposed into the DFT(D) functional, 1 qc drhc (r) b c q qc dr dr C¯ cc (|r − r |)hc (r) + b c c 1 qc qc dr dr C¯ cc (|r − r |)hc (r)hc (r ) + 2b c
DlDFT = −
c
(3)
Y. Maruyama / Journal of Molecular Liquids 291 (2019) 111160
and the CT, qO HS dr dr COO (|r − r |)hO (r) b q2 HS − O dr dr COO (|r − r |)hO ( r)hO (r ) 2b dF HS [qO ] (hO (r) + 1) − lOHS . + DF HS [qO ] − qO dr d(qO hO (r))
DlCT = −
(4)
As mentioned above, the DFT(D) functional, which is second-order approximation, is equivalent to the HNC functional (Eq. (15)). We require a concrete expression for the excess intrinsic free energy functional for the reference HS system, FHS [qO ], to calculate the RMDFT functional or the CT. Therefore, we require an HS fluids theory, for example, the fundamental measure theory [54, 55] or its modified versions [56]. Here, we employ the effective density approximation (EDA) theory [57], and the excess intrinsic free energy functional for the reference HS system is F HS [qO ] = qO
dr(hO (r) + 1)fHS qeff O (r|qO )
(5)
1 g(4 − 3g) b (1 − g)2 (6)
where dHS is the diameter of the reference HS and qeff O (r|qO ) is the effective density, dr WOO (|r − r |)hO (r ). qeff O (r|qO ) =qO + qO ˜ OO (k) = −2bf (qO ) + (2bf (qO ))2 − 4qO bf (qO )C˜ HS (k) / W HS HS HS OO 2qO bfHS (qO )
(7)
Similarly, the excess chemical potential is expressed as
dr WOO (|r − r |).
(8)
Then, the CT is expressed as
(9)
2.2. 3D-RISM theory
cc (r) ∗ wvv c c (r) + qc Hc c (r) .
exp(wc ) − 1 (wc < 0) . wc (wc ≥ 0)
wc = −buc (r) + hc (r) − cc (r)
(12)
The 3D distribution function gc (r) is defined from hc (r) by, gc (r) = hc (r) + 1.
(13)
1 ∂ uc (r; k) 1 qc dk dr gc (r; k). b c ∂k 0
(14)
The coupling parameter, k, gradually switches on the interaction potential from k = 0 (no interaction) to k = 1 (full interaction). uc (r; k) varies according to k, and gc (r; k) denotes the distribution function under uc (r; k). To evaluate this formula with high accuracy, we need to perform 3D-RISM calculations at least 40 times. Singer and Chandler derived the closed form to avoid the necessity of numerical coupling parameter integrations. They analytically integrated over k using RISM and HNC closure equations and showed one dimensional Singer-Chandler formula [29]. The Singer-Chandler formula can be simply extended to three dimensions. Then, the HNC functional is expressed as DlHNC =
1 qc b c
dr
1 2 1 hc (r) − cc (r) − cc (r)hc (r) . 2 2
DlKH =
(15)
1 2 1 1 hc (r)H(−hc (r)) − cc (r) − cc (r)hc (r) , qc dr b c 2 2 (16)
where H denotes Heaviside step function. The Singer-Chandler formula is equivalent to Kirkwood charging formula in each closure equation. Here, we define the sum of the KH functional and the correction terms as KH + CT functional, DlKH + CT = DlKH + DlCT .
To calculate Dl using RMDFT functional, we require hc (r). We employ the 3D-RISM theory to obtain hc (r) under the interaction potential acting on the solvent site, c, of position r, uc (r). For a solutesolvent system at infinite dilution, the 3D-RISM equation is written as
c
hc (r) =
The KH functional is similarly derived and expressed as [30]
1 HS (|r − r |)hO (r) qO dr dr COO b 1 2 HS − dr dr COO (|r − r |)hO (r)hO (r ) q 2b O (qO ) dr dr WOO (|r − r |) + q2O fHS
− q2O dr dr fHS qeff O (r |qO ) WOO (|r − r |)(hO (r) + 1)(hO (r ) + 1).
(11)
or KH closure equation,
DlCT = −
hc (r) =
hc (r) = exp(−buc (r) + hc (r) − cc (r)) − 1,
Dl =
pq(dHS )3 , g= 6
ex = fHS (qO ) + qO fHS (qO ) lHS
where cc (r) is the 3D direct correlation function of the solvent site c around the solute, the asterisk denotes a convolution integral in the real space, wvv c c (r) is the site-site intramolecular correlation function of the solvent, and Hc c (r) is the site-site total correlation functions of pure solvent. The site-site correlation functions of the solvent are obtained in advance from the 1D-RISM theory for pure solvent. The 3D-RISM equation is complemented by HNC closure equation,
One of the methods for obtaining the SFE is to calculate the following Kirkwood charging formula [59],
where fHS (n) is Carnahan-Starling equation [58], fHS (q) =
3
(10)
(17)
Similarly, the sum of the HNC functional and the CT can be defined, but this is numerically equivalent to the RMDFT functional. A modification to Singer-Chandler formula can be obtained in the so-called Gaussian fluctuations (GF) approximation [60], DlGF =
1 1 qc dr −cc (r) − cc (r)hc (r) . b c 2
(18)
4
Y. Maruyama / Journal of Molecular Liquids 291 (2019) 111160
Despite dropping the first term of the Singer-Chandler formula, the GF approximation showed improvement for the solvation thermodynamics in aqueous solution [61, 62]. The partial molar volume [63, 64] is expressed as ⎛ ⎞ jT ⎝ V= 1− qc drcc (r)⎠ , b c
(19)
where j T is the pure solvent isothermal compressibility. V or quantity related to its nonpolar component is used for the partial molar volume corrections. 2.3. Computational details We used two data sets to verify the KH + CT functional for the SFE. The first set consisted of amino acid side chain analogues molecules. We used the OPLS-AA parameter to compare the previously reported computational results [65, 66]. We obtained the input data (structures and partial charges) from the Supporting Information section of the paper by Shirts et al. [65]. The second dataset comprised 504 organic molecules that have been used in several studies on integral equations and DFT [34, 36, 37, 67]. The atomic partial charges and geometries were obtained from the Supporting Information section of the paper by Mobley et al. [68]. In the RISM calculations of Cab (r) and Hab (r), we used the TIP3P model with an additional parameter (s = 0.4 Å and 4 = 0.046 kcal/mol) for water molecule [69, 70]. The number density of water molecules and the temperature were 0.03334 molecule/Å3 and 298.15 K, respectively. We used KH closure in the RISM/3D-RISM calculations except for Table 1 where we used HNC closure. We employed the EDA theory for calculation of CHS OO (r). The diameters of the HS were 2.88 Å for RMDFT functional and 2.85 Å for KH + CT functional. We determined these diameters to optimize the SFE of small hydrophobic solutes, methane, propane, and isobutane. The grid spacing was set to 0.01 Å and the number of grids was 16,384 for the EDA and the 1D-RISM calculations. The 3D-RISM integral equations were solved for a grid of 2563 points in a cubic cell with a size of 64 Å3 on GPU [6]. A grid space of 0.25 Å was sufficient to calculate the SFE without significant numerical errors. 3. Results and discussion In Fig. 1 and Table 2, the SFE values of the neutral amino acid side chain analogues using KH + CT functional are compared with those determined by experiment and using RMDFT functional. The mean absolute deviation (MAD) of the KH + CT results from the experimental values is 0.54 kcal/mol, whereas the MAD of the RMDFT values is 0.48 kcal/mol. Similarly, the root mean square deviation (RMSD) between the KH + CT results and the experimental data is 0.47 kcal/mol, whereas the RMSD of the RMDFT values is 0.41 kcal/mol. Compared with RMDFT, the KH + CT is less accurate but acceptable for practical use. As can be seen in the figure, the features of the amino acid side chain analogues molecules are reproduced. These molecules are ordered on the basis of the experimental values of the SFE. Discrepancies in their ordering between the SFE values and the experimental data were found between Cys and Met and between Ser and Trp. However, similar discrepancies with the experimental data were found in the energy representation method [66]. Fig. 2 shows the correlation between the experimental data and the uncorrected calculated SFEs. SFEs of both the DFT(D) functional and the KH functional have only a weak correlation with the experimental data. We obtained R = 0.408 and s = 7.68 kcal/mol for DFT(D) functional, where R is the correlation coefficient and s is the standard deviation of the residual. For the KH functional, we obtained
Fig. 1. Solvation free energy (SFE) values of neutral amino acid side chain analogues. SFE values calculated using the KH + CT functional are compared with the experimental values [78] and the computational values determined using the RMDFT functional.
R = 0.479 and s = 6.76 kcal/mol. Note that the KH functional is not better, even though its correlation coefficient was slightly higher. Because the same hc (r) and cc (r) were used in these calculations, the difference was caused by the first term of the KH functional. As can be seen in the figure, SFEs are too positive compared to the experimental data, which indicates that both functionals predict these molecules to be too hydrophobic. Fig. 3 shows the correlation between the SFE values of the DFT(D) functional and the KH functional for 504 small organic molecules. The dotted line indicates a linear regression, and the equation is y = 0.909x − 1.78. We obtained R = 0.996 and s = 0.35 kcal/mol, showing a strong correlation. In addition, it can be seen that the KH functional estimates a lower SFE overall. Even if there is a deviation in absolute value in the KH functional and the DFT(D) functional, features of individual molecules are reproduced. Next, we show the correlations between the dimensionless PMV and the differences between the DFT(D) functional and KH functional and the experimental data in Fig. 4. The correlation coefficient between qV and its nonpolar component for 504 organic molecules is 0.9998, which is almost the same. As is well known, there are strong correlations between the PMV (or its nonpolar component) and the energy differences. In fact, it was observed that R = −0.981 and s = 9.18 kcal/mol for the DFT(D) functional, and R = −0.971 Table 2 Solvation free energy values with correction terms of amino acid side chain analogues in water. Amino acid
Analog solute
Experimental
RMDFT
KH + CT
Leu Ile Val Ala Phe Cys Met Thr Ser Trp Tyr Gln Asn Hid Hie
Isobutane n-Butane Propane Methane Toluene Methanethiol Methyl ethyl sulfide Ethanol Methanol 3-Methylindole p-Cresol Propionamide Acetamide 4-Methylimidazole 4-Methylimidazole
2.3 2.2 2.0 1.9 −0.8 −1.2 −1.5 −4.9 −5.1 −5.9 −6.1 −9.4 −9.7 −10.3 −10.3
2.5 2.4 2.2 1.7 −1.5 −2.6 −2.0 −4.7 −5.7 −4.6 −6.0 −8.4 −9.5 −9.8 −10.3
2.6 2.5 2.2 1.7 −1.4 −2.7 −2.0 −5.3 −6.3 −4.7 −6.3 −8.9 −10.2 −10.3 −10.8
Energy unit is kcal/mol. Experimental data were taken from the literature [78] . We used Kovalenko-Hirata closure in the calculations.
Y. Maruyama / Journal of Molecular Liquids 291 (2019) 111160
a
5
b
Fig. 2. Comparison of solvation free energy values of 504 small organic molecules obtained using the DFT(D) functional (left) and the KH functional (right) and the experimentally determined values.
and s = 8.24 kcal/mol for the KH functional. (To be precise, those are negative correlations, but the negative sign is only a matter of how differences are taken.) The linear regressions for differences of the DFT(D) functional and the KH functional are y = −5.04x and y = −4.30x, respectively. The inclination coefficients corresponded to the parameter a of the PMV correction. Note that the difference between the experimental SFE and the SFE calculated using the GF formula is also strongly correlated with the dimensionless partial molar volume [33]. Our main results are shown in Fig. 5, which contains the correlations between qV and the CT values of the RMDFT functional (CT(DFT(D))) and the KH + CT functional (CT(KH)). The linear
Fig. 3. Correlation between the solvation free energy values of the DFT(D) functional and the KH functional for 504 small organic molecules.
regressions for the RMDFT functional and the KH + CT functional were y = −5.06x and y = −4.34x, respectively. Both had almost the same inclination values as in Fig. 4. This indicates that the CT based on the RMDFT is suitable for the correction of the SFE. It was also observed that R = −0.998 and s = 9.13 kcal/mol for the RMDFT functional, and R = −0.998 and s = 8.05 kcal/mol for the KH + CT functional. Both CTs have strong correlation with PMV. (Those are negative correlations, but negative sign is only a matter of how a coefficient is taken.) This shows that both CTs have the same effect as the PMV correction. The inclination coefficients also corresponded to the parameter a of the PMV correction. In the RMDFT functional, the diameter of the reference HS plays the role of this parameter. In short, the CT of the RMDFT functional provides the contribution of cavity formation by the solute depending on the size of the reference HS. Fig. 6 shows comparisons of the calculation results and the experimental data for the SFE values of 504 small organic molecules in water [68, 71]. Both SFE results agree well with the experimental data, especially for hydrophobic solutes that have positive or small SFE values. However, there are differences in the degree of variation between the RMDFT and KH + CT functionals. The MADs of the RMDFT and KH + CT functionals from experimental data were 1.12 and 1.22 kcal/mol, respectively. Similarly, the RMSDs for RMDFT and KH + CT functionals were 1.47 and 2.85 kcal/mol, respectively. Although the effects of both CTs are similar, the result of the RMDFT functional was better than that of the KH + CT functional. In contrast, the correlation coefficient and the standard deviation of RMDFT results from experimental data were R = 0.907 and s = 1.48 kcal/mol, and R = 0.912 and s = 1.61 kcal/mol for KH + CT functional. The RMDFT is somewhat worse from the correlation coefficient viewpoint. Fig. 7 shows the correlation between the SFE values of the RMDFT and the KH + CT functionals for the 504 small organic molecules. The linear regression was y = 1.08x − 0.191 and the inclination exceeded 1, in contrast to the results shown in Fig. 3. However, the reciprocal of the slope of the correlation between the DFT(D) functional and the KH functional was 1.10(= 1/0.909), which was the same degree of deviation. This is essentially the difference between the HNC functional (which is equivalent to DFT(D)) and KH functional. The difference between the HNC functional and
6
Y. Maruyama / Journal of Molecular Liquids 291 (2019) 111160
a
b
Fig. 4. Correlation between the dimensionless partial molar volume and differences between the experimental data and the calculation values using the DFT(D) functional (left) and the KH functional (right).
the KH functional is only the first term of the integral (see Eqs. (15) and (16)). In other words, the deviation of the two SingerChandler formulas is caused by this term. In contrast to the inclination, we observed R = 0.998, and s = 0.351 kcal/mol, indicating a strong correlation.
4. Conclusion We divided the RMDFT functional into the DFT(D) and CT parts and investigated their contributions. We also investigated its usefulness by combining the KH functional and the CT to confirm consistency with the closure equation. The CT estimates the correction of the cavity formation of the solute by the reference HS. The contribution of the CT increases with the diameter of the HS of the reference system, and it strongly correlated with
a
the PMV. As the result, the CT based on the RMDFT showed the same effect as PMV correction, although the RMDFT functional does not contain the parameter b of the PMV correction. On the other hand, the KH + CT functional showed an accuracy comparable to that of the RMDFT functional. The result shows that consistency between the SFE functional and the closure equation is unnecessary. SFE calculations and application by the 3D-RISM theory is still under development. Recently, Kobryn et al. derived a new KGK closure for the 3D-RISM theory that is particularly suitable for polar and charged molecules in electrolyte solutions. The KGK closure levels out the density distribution function to zero inside the repulsive core to avoid the problem of the HNC-type approximation [72]. Sosnin et al. demonstrated that the bioconcentration factor can be effectively predicted by the combination of the 3D-RISM theory and machine learning [73].
b
Fig. 5. Correlation between the dimensionless partial molar volume and the values of correction terms for the DFT(D) functional (left) and for the KH functional (right).
Y. Maruyama / Journal of Molecular Liquids 291 (2019) 111160
a
7
b
Fig. 6. Comparison of the solvation free energy values of 504 small organic molecules obtained using the RMDFT (left) and the KH + CT (right) functionals and the experimentally determined values.
Improvement of 3D distribution function is essential to improve the accuracy of SFE of hydrophilic molecules because the contribution of the Coulomb interaction between solute and solvent is large. One way to do this is to introduce a bridge correction term in the HNC-type closure equation to overcome its softness. Such a term is obtained as a functional derivative of the SFE correction term [74, 75]. That is, the bridge correction term for the closure equation is in correspondence with the correction term of the SFE functional [76, 77]. In fact, we derived the HS bridge correction term based on RMDFT and applied it to Lennard-Jones fluids [47]. In addition, we can systematically improve the SFE functional and the closure equation by using RMDFT. However, it is currently difficult to converge 3D-RISM calculation with the RMDFT closure equation including the HS bridge correction term.
We may need to develop a new convergence scheme or convergence acceleration method for the 3D-RISM calculation with the RMDFT closure equation.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27]
Fig. 7. Correlation between the solvation free energy values of the RMDFT and KH + CT functionals for 504 small organic molecules.
[28] [29]
D. Beglov, B. Roux, J. Phys. Chem. B 101 (1997) 7821–7826. A. Kovalenko, F. Hirata, Chem. Phys. Lett. 290 (1998) 237–244. A. Kovalenko, F. Hirata, J. Phys. Chem. B 103 (1999) 7942–7957. S. Gusarov, B.S.B. Pujari, A. Kovalenko, J. Comput. Chem. 33 (2012) 1478–1494. V.P. Sergiievskyi, M.V. Fedorov, J. Chem. Theory Comput. 8 (2012) 2062–2070. Y. Maruyama, F. Hirata, J. Chem. Theory Comput. 8 (2012) 3015–3021. Y. Maruyama, N. Yoshida, H. Tadano, D. Takahashi, M. Sato, F. Hirata, J. Comput. Chem. 35 (2014) 1347–1355. D. Casanova, S. Gusarov, A. Kovalenko, T. Ziegler, J. Chem. Theory Comput. 3 (2007) 458–476. T. Luchko, S. Gusarov, D.R. Roe, C. Simmerling, D.A. Case, J. Tuszynski, A. Kovalenko, J. Chem. Theory Comput. 6 (2010) 607–624. N. Yoshida, S. Phongphanphanee, Y. Maruyama, T. Imai, F. Hirata, J. Am. Chem. Soc. 128 (2006) 12042–12043. N. Yoshida, S. Phongphanphanee, F. Hirata, J. Phys. Chem. B 111 (2007) 4588–4595. T. Yamazaki, T. Imai, F. Hirata, A. Kovalenko, J. Phys. Chem. B 111 (2007) 1206–1212. S. Phongphanphanee, N. Yoshida, F. Hirata, J. Am. Chem. Soc. 130 (2008) 1540–1541. S. Phongphanphanee, N. Yoshida, F. Hirata, J. Mol. Liq. 147 (2009) 107–111. S. Phongphanphanee, N. Yoshida, F. Hirata, J. Phys. Chem. B 114 (2010) 7967–7973. R.S. Johnson, T. Yamazaki, A. Kovalenko, H. Fenniri, J. Am. Chem. Soc. 129 (2007) 5735–5743. G. Tikhomirov, T. Yamazaki, A. Kovalenko, H. Fenniri, Langmuir 24 (2008) 4447–4450. T. Yamazaki, H. Fenniri, A. Kovalenko, ChemPhysChem 11 (2010) 361–367. G. Borzsonyi, R.S. Johnson, A.J. Myles, J.-Y. Cho, T. Yamazaki, R.L. Beingessner, A. Kovalenko, H. Fenniri, Chem. Commun. 46 (2010) 6527–6529. T. Yamazaki, H. Fenniri, J. Mol. Liq. 217 (2016) 70–74. T. Yamazaki, N. Blinov, D. Wishart, A. Kovalenko, Biophys. J. 95 (2008) 4540–4548. S.-H. Chong, S. Ham, Proc. Natl. Acad. Sci. 109 (2012) 7636–7641. S.-H. Chong, S. Ham, Angew. Chem. 126 (2014) 1–5. D.J. Sindhikara, N. Yoshida, M. Kataoka, F. Hirata, J. Mol. Liq. 164 (2011) 120–122. Y. Maruyama, N. Yoshida, F. Hirata, J. Phys. Chem. B 114 (2010) 6464–6471. G.M. Giambasu, ¸ M.K. Gebala, M.T. Panteva, T. Luchko, D.A. Case, D.M. York, Nucleic Acids Res. 43 (2015) 8405–8415. Y. Maruyama, T. Matsushita, R. Ueoka, F. Hirata, J. Phys. Chem. B 115 (2011) 2408–2416. J. Lee, H. Im, S.-H. Chong, S. Ham, Chem. Phys. Lett. 693 (2018) 216–221. S. Singer, D. Chandler, Mol. Phys. 55 (1985) 621–625.
8
Y. Maruyama / Journal of Molecular Liquids 291 (2019) 111160 [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53]
A. Kovalenko, F. Hirata, J. Chem. Phys. 110 (1999) 10095–10112. S. Ten-No, J. Chem. Phys. 115 (2001) 3724. K. Sato, H. Chuman, S. Ten-No, J. Phys. Chem. B 109 (2005) 17290–17295. D.S. Palmer, A.I. Frolov, E.L. Ratkova, M.V. Fedorov, J. Phys. Condens. Matter 22 (2010) 492101. D.S. Palmer, A.I. Frolov, E.L. Ratkova, M.V. Fedorov, Mol. Pharm. 8 (2011) 1423–1429. A.I. Frolov, E.L. Ratkova, D.S. Palmer, M.V. Fedorov, J. Phys. Chem. B 115 (2011) 6011–6022. J.-F. Truchon, B.M. Pettitt, P. Labute, J. Chem. Theory Comput. 10 (2014) 934–941. V.P. Sergiievskyi, G. Jeanmairet, M. Levesque, D. Borgis, J. Phys. Chem. Lett. 5 (2014) 1935–1942. V.P. Sergiievskyi, G. Jeanmairet, M. Levesque, D. Borgis, J. Chem. Phys. 143 (2015) 184116–184117. M. Misin, D.S. Palmer, M.V. Fedorov, J. Phys. Chem. B 120 (2016) 5724–5731. M. Kinoshita, Y. Harano, R. Akiyama, J. Chem. Phys. 125 (2006) 244504–244507. S.-H. Chong, S. Ham, Chem. Phys. Lett. 535 (2012) 152–156. S.-H. Chong, S. Ham, J. Chem. Theory Comput. 11 (2015) 378–380. T. Miyata, J. Thapa, Chem. Phys. Lett. 604 (2014) 122–126. T. Miyata, Y. Ebato, J. Mol. Liq. 217 (2016) 75–82. T. Sumi, A. Mitsutake, Y. Maruyama, J. Comput. Chem. 36 (2015) 1359–1369. T. Sumi, A. Mitsutake, Y. Maruyama, J. Comput. Chem. (2015) 2009–2011. T. Sumi, Y. Maruyama, A. Mitsutake, K. Koga, J. Chem. Phys. 144 (2016) 224104–224116. Y. Maruyama, A. Mitsutake, J. Phys. Chem. B 121 (2017) 9881–9885. T. Sumi, Y. Maruyama, A. Mitsutake, K. Mochizuki, K. Koga, J. Comput. Chem. 39 (2018) 202–217. J.P. Hansen, I.R. McDonald, Theory of Simple Liquids, Third ed. ed., Elsevier/Academic Press, London, 2006. D. Henderson (Ed.), Fundamentals of Inhomogeneous Fluids, Marcel Dekker, New York., 1992, A. Kovalenko, F. Hirata, J. Chem. Phys. 112 (2000) 10391–10402. D. Chandler, J. Mccoy, S. Singer, J. Chem. Phys. 85 (1986) 5971–5976.
[54] Y. Rosenfeld, Phys. Rev. Lett. 63 (1989) 980–983. [55] R. Roth, R. Evans, A. Lang, G. Kahl, J. Phys. Condens. Matter 14 (2002) 12063–12078. [56] Y.-X. Yu, J. Wu, J. Chem. Phys. 117 (2002) 10156–10164. [57] T. Sumi, H. Sekino, J. Phys. Soc. Jpn. 77 (2008) 034605. [58] N.F. Carnahan, K.E. Starling, J. Chem. Phys. 51 (1969) 635–636. [59] A. Ben-Naim, Molecular Theory of Solutions, Oxford University Press, New York, 2006. [60] D. Chandler, Y. Singh, D.M. Richardson, J. Chem. Phys. 81 (1984) 1975. [61] D. Chandler, T. Ichiye, J. Phys. Chem. 92 (1988) 5257. [62] P.H. Lee, G.M. Maggiora, J. Phys. Chem. 97 (1993) 10175. [63] T. Imai, M. Kinoshita, F. Hirata, J. Chem. Phys. 112 (2000) 9469–9478. [64] T. Imai, Y. Harano, A. Kovalenko, F. Hirata, Biopolymers 59 (2001) 512–519. [65] M.R. Shirts, J.W. Pitera, W.C. Swope, V.S. Pande, J. Chem. Phys. 119 (2003) 5740–5761. [66] Y. Karino, M.V. Fedorov, N. Matubayasi, Chem. Phys. Lett. 496 (2010) 351–355. [67] Y. Liu, J. Fu, J. Wu, J. Phys. Chem. Lett. 4 (2013) 3687–3691. [68] D.L. Mobley, C.I. Bayly, M.D. Cooper, M.R. Shirts, K.A. Dill, J. Chem. Theory Comput. 5 (2009) 350–358. [69] W. Jorgensen, J. Chandrasekhar, J. Madura, R. Impey, M. Klein, J. Chem. Phys. 79 (1983) 926–935. [70] B. Pettitt, P. Rossky, J. Chem. Phys. 77 (1982) 1451–1457. [71] R.C. Rizzo, T. Aynechi, D.A. Case, I.D. Kuntz, J. Chem. Theory Comput. 2 (2006) 128–139. [72] A.E. Kobryn, S. Gusarov, A. Kovalenko, J. Phys. Condens. Matter 28 (2016) 404003–404024. [73] S. Sosnin, M. Misin, D.S. Palmer, M.V. Fedorov, J. Phys. Condens. Matter 30 (2018) 32LT03. [74] Y. Liu, S. Zhao, J. Wu, J. Chem. Theory Comput. 9 (2013) 1896–1908. [75] Y. Liu, J. Wu, J. Chem. Phys. 140 (2014) 084103–084108. [76] A. Kovalenko, F. Hirata, J. Chem. Phys. 113 (2000) 2793–2805. [77] D. Yokogawa, J. Chem. Theory Comput. 14 (2018) 3272–3278. [78] R. Wolfenden, L. Andersson, P.M. Cullis, C.C. Southgate, Biochemistry 20 (1981) 849–855.
Contents lists available at ScienceDirect
Journal of Molecular Liquids journal homepage: www.elsevier.com/locate/molliq
Correction terms for the solvation free energy functional of three-dimensional reference interaction site model based on the reference-modified density functional theory夽 Yutaka Maruyama Architecture Development Team, FLAGSHIP 2020 Project, RIKEN Center for Computational Science, Kobe 650-0047, Japan
A R T I C L E
I N F O
Article history: Received 28 February 2019 Received in revised form 13 May 2019 Accepted 8 June 2019 Available online 26 June 2019 Keywords: Solvation free energy 3D-RISM theory Classical DFT Water solvent Organic molecules
A B S T R A C T Recently, we derived the reference-modified density functional theory (RMDFT) to calculate the solvation free energy (SFE) with high accuracy, introducing a reference system to the DFT for polyatomic molecular liquids. The RMDFT functional consists of the DFT using second-order approximation (DFT(D)) functional and the correction terms. The DFT(D) functional are equivalent to the Singer-Chandler formula for the hypernetted chain closure equation. The contribution of the correction terms has a strong correlation with the partial molar volume (PMV) of solute. Thus, the RMDFT functional has almost the same accuracy as the PMV correction. Furthermore, we combined the Singer-Chandler formula for the Kovalenko-Hirata closure and the correction terms to avoid inconsistency between the closure equation and the SFE functional. It become clear that the inconsistency is unimportant, whereas the combined functional works well for the amino acid side chain analogues molecules. © 2019 Elsevier B.V. All rights reserved.
1. Introduction The solvation free energy (SFE) is an important property for investigation of the thermodynamic stabilities of systems, in addition to chemical reaction, protein folding, and protein-ligand binding. One way to calculate the SFE is to employ the three-dimensional reference interaction site model (3D-RISM) theory, which is a statistical mechanics theory for molecular liquids [1-3]. To increase the speed of 3D-RISM calculations, a highly efficient multigrid algorithm was developed [4, 5]. In addition, the 3D-RISM program is implemented using graphics processing unit (GPU) [6] and massively parallel supercomputer [7]. It is also implemented in commercial programs such as the Amsterdam Density Functional (ADF) modeling suite [8] and the Amber package [9]. Therefore, the 3D-RISM theory has been applied to various biological systems, such as, for example, ion binding protein [10, 11], staphylococcal nuclease [12], aquaporin [13-15], rosette nanotubes [16-20], amyloidb (Ab42) [21-23], photoactive yellow protein [24], DNA [25, 26], and telomere [27, 28].
夽 The article was accepted for publication under the “VSI: EMLG/JMLG 2018 in Nagoya”. E-mail address: [email protected] (Y. Maruyama).
https://doi.org/10.1016/j.molliq.2019.111160 0167-7322/© 2019 Elsevier B.V. All rights reserved.
The Singer-Chandler formula [29, 30] is commonly used when calculating the SFE using the RISM/3D-RISM theory. However, Ten-no pointed out that the absolute values of the SFE artificially increase with the number of sites of solute molecules, a theoretical drawback of the RISM theory [31, 32]. Therefore, some researchers proposed introducing a phenomenological correction in SFE calculation methods in which the partial molar volume of the solute molecule cancels the artificial solute-size dependence of the SFE [33-35]. Similarly, a cavity correction [36] and a pressure correction are also proposed for the SFE [37-39]. However, which of these correction method is the best remains controversial. Furthermore, the differences between isochoric and isobaric conditions and ensembleindependence of the SFE are discussed in detail [40-42]. Recently, Miyata and Thapa pointed out that one of the problems on the RISM theory should be attributed to the theoretical drawback of the hypernetted-chain (HNC) type approximation itself rather than to the simple extension of the HNC-type closure. They showed that the SFE of a Lennard-Jones (LJ) solute in an LJ solvent depends on solute size by using the Ornstein-Zernike (OZ) integral equation with the HNC and Kovalenko-Hirata (KH) closures. SFE values obtained from the HNC approximation and the KH approximation increase monotonically as the solute size increases, whereas the SFE determined using molecular dynamics (MD) simulations decrease [43]. They actually applied an approximate bridge function adjusting the s parameter in LJ potential of
2
Y. Maruyama / Journal of Molecular Liquids 291 (2019) 111160
the HNC-type closure equation to obtain a more accurate thermodynamic quantity [44]. Thus, not only the RISM theory for molecular liquids but also the OZ theory for simple liquids results in solute-size dependence and significant overestimation of the SFE using the HNC or KH closure. Inspired by their study, we derived a new SFE functional using another approach. We proposed a reference-modified density functional theory (RMDFT) to systematically develop a free-energy density functional model [45-47]. We employ a hardsphere (HS) fluid as the reference system instead of an ideal polyatomic molecular gas to derive the SFE functional of a solute molecule in water. The SFE functional derived from the RMDFT approach improved the overestimation of the SFE originating from the HNC-type approximation. The RMDFT functional improved the absolute values of the SFE for neutral amino acid side chain analog molecules and for 504 small organic molecules [45, 46]. Furthermore, we applied the RMDFT functional to investigate the Lennard-Jones solution [47], the structural stability of protein during folding process [48], and the high predictive ability for temperature and pressure dependences of the SFE [49]. The RMDFT functional is composed of two parts: the secondorder approximation DFT (DFT(D)) functional and the correction terms (CT). The DFT(D) functional is equivalent to the SingerChandler formula for HNC closure (HNC functional) [50, 51]. We actually compared the DFT(D) functional and the HNC functional in the SFE calculation for neutral amino acid side chain analogues molecules (Table 1). Despite the deviation from the experimental values, all values of the DFT functional and the HNC functional agree with an accuracy of 0.1 kcal/mol or less. If the DFT functional contains a higher order approximation, deviation of the DFT functional from the HNC functional will occur. However, the higher order terms do not contribute significantly to improve the calculation value [47]. For closing the 3D-RISM equation, we chose the KH closure equation [3, 52] because of its advantage of rapid convergence, whereas the HNC closure shows an extremely poor convergence. The DFT(D) functional is equivalent to the HNC functional, but not SingerChandler formula for the KH closure (KH functional). In other words, when using RMDFT functional with KH closure equation, there is an inconsistency between the closure equation and the SFE functional because the Singer-Chandler formula has a different form for each closure equation. In this study, we combine the KH functional and the CT based RMDFT for consistency. We calculate and verify the SFE values for
Table 1 Solvation free energy values of amino acid side chain analogues in water. Amino acid
Analog solute
Experimental
DFT(D)
HNC
Leu Ile Val Ala Phe Cys Met Thr Ser Trp Tyr Gln Asn Hid Hie
Isobutane n-Butane Propane Methane Toluene Methanethiol Methyl ethyl sulfide Ethanol Methanol 3-Methylindole p-Cresol Propionamide Acetamide 4-Methylimidazole 4-Methylimidazole
2.3 2.2 2.0 1.9 −0.8 −1.2 −1.5 −4.9 −5.1 −5.9 −6.1 −9.4 −9.7 −10.3 −10.3
17.2 17.2 14.3 8.1 16.5 6.1 13.2 4.6 0.7 17.1 11.9 4.2 0.4 4.0 3.5
17.2 17.2 14.3 8.1 16.5 6.1 13.2 4.6 0.7 17.1 11.9 4.2 0.4 4.0 3.5
Energy unit is kcal/mol. Experimental data were taken from the literature [78] . We used hypernetted-chain closure in the calculations.
amino acid side chain analogues and 504 organic molecules using the KH functional with the CT. We also investigate how the CT of RMDFT functional works.
2. Methods 2.1. RMDFT functional We derived RMDFT by introducing a reference system to the DFT for polyatomic molecular liquids. Here, we consider the SFE functional for water solvent. Then, we incorporate the reference HS only in water oxygen because water hydrogen is buried in the oxygen atom. The RMDFT functional of the SFE for water solvent is expressed as [45, 46] 1 qc drhc (r) b c q ex qc dr dr Ccc + (|r − r |)hc (r) b c c 1 ex qc qc dr dr Ccc + (|r − r |)hc (r)hc (r ) 2b c c dF HS [qO ] (hO (r) + 1) − lOHS , + DF HS [qO ] − qO dr d(qO hO (r))
DlRMDFT = −
(1)
where ex Cab (|r
− r |) =
HS C¯ OO (|r − r |) − COO (|r − r |) (a = b = O) C¯ ab (|r − r |) (otherwise)
IM C¯ ab (|r − r |) = Cab (|r − r |) + Cab (|r − r |) −1 dab sin(kLab ) C˜ ab (k) = − dab q + (1 − dab )q . q kLab
(2)
Here, b = 1/(kB T) is the inverse of the Boltzmann constant times temperature, qc denotes the average number density of the solvent site, c, q is the average number density of the water solvent, hc (r) is the 3D total correlation function of the solvent site c around solute, Cab (r) is the site-site direct correlation function of pure solvent, IM Cab (r) is the intramolecular direct correlation function [53], Lab is the bond length between a and b sites, O denotes oxygen site of water solvent, CHS OO (r) correspond to the site-site direct correlation function of the reference HS fluid, FHS [qO ] is the excess intrinsic free energy functional for HS fluid, and l HS O is the excess chemical potential of the reference HS fluid. We obtain hc (r) by 3D-RISM calculations, Cab (r) by solvent-solvent RISM calculations, and CHS OO (r) by hardsphere fluids theory. The RMDFT functional is decomposed into the DFT(D) functional, 1 qc drhc (r) b c q qc dr dr C¯ cc (|r − r |)hc (r) + b c c 1 qc qc dr dr C¯ cc (|r − r |)hc (r)hc (r ) + 2b c
DlDFT = −
c
(3)
Y. Maruyama / Journal of Molecular Liquids 291 (2019) 111160
and the CT, qO HS dr dr COO (|r − r |)hO (r) b q2 HS − O dr dr COO (|r − r |)hO ( r)hO (r ) 2b dF HS [qO ] (hO (r) + 1) − lOHS . + DF HS [qO ] − qO dr d(qO hO (r))
DlCT = −
(4)
As mentioned above, the DFT(D) functional, which is second-order approximation, is equivalent to the HNC functional (Eq. (15)). We require a concrete expression for the excess intrinsic free energy functional for the reference HS system, FHS [qO ], to calculate the RMDFT functional or the CT. Therefore, we require an HS fluids theory, for example, the fundamental measure theory [54, 55] or its modified versions [56]. Here, we employ the effective density approximation (EDA) theory [57], and the excess intrinsic free energy functional for the reference HS system is F HS [qO ] = qO
dr(hO (r) + 1)fHS qeff O (r|qO )
(5)
1 g(4 − 3g) b (1 − g)2 (6)
where dHS is the diameter of the reference HS and qeff O (r|qO ) is the effective density, dr WOO (|r − r |)hO (r ). qeff O (r|qO ) =qO + qO ˜ OO (k) = −2bf (qO ) + (2bf (qO ))2 − 4qO bf (qO )C˜ HS (k) / W HS HS HS OO 2qO bfHS (qO )
(7)
Similarly, the excess chemical potential is expressed as
dr WOO (|r − r |).
(8)
Then, the CT is expressed as
(9)
2.2. 3D-RISM theory
cc (r) ∗ wvv c c (r) + qc Hc c (r) .
exp(wc ) − 1 (wc < 0) . wc (wc ≥ 0)
wc = −buc (r) + hc (r) − cc (r)
(12)
The 3D distribution function gc (r) is defined from hc (r) by, gc (r) = hc (r) + 1.
(13)
1 ∂ uc (r; k) 1 qc dk dr gc (r; k). b c ∂k 0
(14)
The coupling parameter, k, gradually switches on the interaction potential from k = 0 (no interaction) to k = 1 (full interaction). uc (r; k) varies according to k, and gc (r; k) denotes the distribution function under uc (r; k). To evaluate this formula with high accuracy, we need to perform 3D-RISM calculations at least 40 times. Singer and Chandler derived the closed form to avoid the necessity of numerical coupling parameter integrations. They analytically integrated over k using RISM and HNC closure equations and showed one dimensional Singer-Chandler formula [29]. The Singer-Chandler formula can be simply extended to three dimensions. Then, the HNC functional is expressed as DlHNC =
1 qc b c
dr
1 2 1 hc (r) − cc (r) − cc (r)hc (r) . 2 2
DlKH =
(15)
1 2 1 1 hc (r)H(−hc (r)) − cc (r) − cc (r)hc (r) , qc dr b c 2 2 (16)
where H denotes Heaviside step function. The Singer-Chandler formula is equivalent to Kirkwood charging formula in each closure equation. Here, we define the sum of the KH functional and the correction terms as KH + CT functional, DlKH + CT = DlKH + DlCT .
To calculate Dl using RMDFT functional, we require hc (r). We employ the 3D-RISM theory to obtain hc (r) under the interaction potential acting on the solvent site, c, of position r, uc (r). For a solutesolvent system at infinite dilution, the 3D-RISM equation is written as
c
hc (r) =
The KH functional is similarly derived and expressed as [30]
1 HS (|r − r |)hO (r) qO dr dr COO b 1 2 HS − dr dr COO (|r − r |)hO (r)hO (r ) q 2b O (qO ) dr dr WOO (|r − r |) + q2O fHS
− q2O dr dr fHS qeff O (r |qO ) WOO (|r − r |)(hO (r) + 1)(hO (r ) + 1).
(11)
or KH closure equation,
DlCT = −
hc (r) =
hc (r) = exp(−buc (r) + hc (r) − cc (r)) − 1,
Dl =
pq(dHS )3 , g= 6
ex = fHS (qO ) + qO fHS (qO ) lHS
where cc (r) is the 3D direct correlation function of the solvent site c around the solute, the asterisk denotes a convolution integral in the real space, wvv c c (r) is the site-site intramolecular correlation function of the solvent, and Hc c (r) is the site-site total correlation functions of pure solvent. The site-site correlation functions of the solvent are obtained in advance from the 1D-RISM theory for pure solvent. The 3D-RISM equation is complemented by HNC closure equation,
One of the methods for obtaining the SFE is to calculate the following Kirkwood charging formula [59],
where fHS (n) is Carnahan-Starling equation [58], fHS (q) =
3
(10)
(17)
Similarly, the sum of the HNC functional and the CT can be defined, but this is numerically equivalent to the RMDFT functional. A modification to Singer-Chandler formula can be obtained in the so-called Gaussian fluctuations (GF) approximation [60], DlGF =
1 1 qc dr −cc (r) − cc (r)hc (r) . b c 2
(18)
4
Y. Maruyama / Journal of Molecular Liquids 291 (2019) 111160
Despite dropping the first term of the Singer-Chandler formula, the GF approximation showed improvement for the solvation thermodynamics in aqueous solution [61, 62]. The partial molar volume [63, 64] is expressed as ⎛ ⎞ jT ⎝ V= 1− qc drcc (r)⎠ , b c
(19)
where j T is the pure solvent isothermal compressibility. V or quantity related to its nonpolar component is used for the partial molar volume corrections. 2.3. Computational details We used two data sets to verify the KH + CT functional for the SFE. The first set consisted of amino acid side chain analogues molecules. We used the OPLS-AA parameter to compare the previously reported computational results [65, 66]. We obtained the input data (structures and partial charges) from the Supporting Information section of the paper by Shirts et al. [65]. The second dataset comprised 504 organic molecules that have been used in several studies on integral equations and DFT [34, 36, 37, 67]. The atomic partial charges and geometries were obtained from the Supporting Information section of the paper by Mobley et al. [68]. In the RISM calculations of Cab (r) and Hab (r), we used the TIP3P model with an additional parameter (s = 0.4 Å and 4 = 0.046 kcal/mol) for water molecule [69, 70]. The number density of water molecules and the temperature were 0.03334 molecule/Å3 and 298.15 K, respectively. We used KH closure in the RISM/3D-RISM calculations except for Table 1 where we used HNC closure. We employed the EDA theory for calculation of CHS OO (r). The diameters of the HS were 2.88 Å for RMDFT functional and 2.85 Å for KH + CT functional. We determined these diameters to optimize the SFE of small hydrophobic solutes, methane, propane, and isobutane. The grid spacing was set to 0.01 Å and the number of grids was 16,384 for the EDA and the 1D-RISM calculations. The 3D-RISM integral equations were solved for a grid of 2563 points in a cubic cell with a size of 64 Å3 on GPU [6]. A grid space of 0.25 Å was sufficient to calculate the SFE without significant numerical errors. 3. Results and discussion In Fig. 1 and Table 2, the SFE values of the neutral amino acid side chain analogues using KH + CT functional are compared with those determined by experiment and using RMDFT functional. The mean absolute deviation (MAD) of the KH + CT results from the experimental values is 0.54 kcal/mol, whereas the MAD of the RMDFT values is 0.48 kcal/mol. Similarly, the root mean square deviation (RMSD) between the KH + CT results and the experimental data is 0.47 kcal/mol, whereas the RMSD of the RMDFT values is 0.41 kcal/mol. Compared with RMDFT, the KH + CT is less accurate but acceptable for practical use. As can be seen in the figure, the features of the amino acid side chain analogues molecules are reproduced. These molecules are ordered on the basis of the experimental values of the SFE. Discrepancies in their ordering between the SFE values and the experimental data were found between Cys and Met and between Ser and Trp. However, similar discrepancies with the experimental data were found in the energy representation method [66]. Fig. 2 shows the correlation between the experimental data and the uncorrected calculated SFEs. SFEs of both the DFT(D) functional and the KH functional have only a weak correlation with the experimental data. We obtained R = 0.408 and s = 7.68 kcal/mol for DFT(D) functional, where R is the correlation coefficient and s is the standard deviation of the residual. For the KH functional, we obtained
Fig. 1. Solvation free energy (SFE) values of neutral amino acid side chain analogues. SFE values calculated using the KH + CT functional are compared with the experimental values [78] and the computational values determined using the RMDFT functional.
R = 0.479 and s = 6.76 kcal/mol. Note that the KH functional is not better, even though its correlation coefficient was slightly higher. Because the same hc (r) and cc (r) were used in these calculations, the difference was caused by the first term of the KH functional. As can be seen in the figure, SFEs are too positive compared to the experimental data, which indicates that both functionals predict these molecules to be too hydrophobic. Fig. 3 shows the correlation between the SFE values of the DFT(D) functional and the KH functional for 504 small organic molecules. The dotted line indicates a linear regression, and the equation is y = 0.909x − 1.78. We obtained R = 0.996 and s = 0.35 kcal/mol, showing a strong correlation. In addition, it can be seen that the KH functional estimates a lower SFE overall. Even if there is a deviation in absolute value in the KH functional and the DFT(D) functional, features of individual molecules are reproduced. Next, we show the correlations between the dimensionless PMV and the differences between the DFT(D) functional and KH functional and the experimental data in Fig. 4. The correlation coefficient between qV and its nonpolar component for 504 organic molecules is 0.9998, which is almost the same. As is well known, there are strong correlations between the PMV (or its nonpolar component) and the energy differences. In fact, it was observed that R = −0.981 and s = 9.18 kcal/mol for the DFT(D) functional, and R = −0.971 Table 2 Solvation free energy values with correction terms of amino acid side chain analogues in water. Amino acid
Analog solute
Experimental
RMDFT
KH + CT
Leu Ile Val Ala Phe Cys Met Thr Ser Trp Tyr Gln Asn Hid Hie
Isobutane n-Butane Propane Methane Toluene Methanethiol Methyl ethyl sulfide Ethanol Methanol 3-Methylindole p-Cresol Propionamide Acetamide 4-Methylimidazole 4-Methylimidazole
2.3 2.2 2.0 1.9 −0.8 −1.2 −1.5 −4.9 −5.1 −5.9 −6.1 −9.4 −9.7 −10.3 −10.3
2.5 2.4 2.2 1.7 −1.5 −2.6 −2.0 −4.7 −5.7 −4.6 −6.0 −8.4 −9.5 −9.8 −10.3
2.6 2.5 2.2 1.7 −1.4 −2.7 −2.0 −5.3 −6.3 −4.7 −6.3 −8.9 −10.2 −10.3 −10.8
Energy unit is kcal/mol. Experimental data were taken from the literature [78] . We used Kovalenko-Hirata closure in the calculations.
Y. Maruyama / Journal of Molecular Liquids 291 (2019) 111160
a
5
b
Fig. 2. Comparison of solvation free energy values of 504 small organic molecules obtained using the DFT(D) functional (left) and the KH functional (right) and the experimentally determined values.
and s = 8.24 kcal/mol for the KH functional. (To be precise, those are negative correlations, but the negative sign is only a matter of how differences are taken.) The linear regressions for differences of the DFT(D) functional and the KH functional are y = −5.04x and y = −4.30x, respectively. The inclination coefficients corresponded to the parameter a of the PMV correction. Note that the difference between the experimental SFE and the SFE calculated using the GF formula is also strongly correlated with the dimensionless partial molar volume [33]. Our main results are shown in Fig. 5, which contains the correlations between qV and the CT values of the RMDFT functional (CT(DFT(D))) and the KH + CT functional (CT(KH)). The linear
Fig. 3. Correlation between the solvation free energy values of the DFT(D) functional and the KH functional for 504 small organic molecules.
regressions for the RMDFT functional and the KH + CT functional were y = −5.06x and y = −4.34x, respectively. Both had almost the same inclination values as in Fig. 4. This indicates that the CT based on the RMDFT is suitable for the correction of the SFE. It was also observed that R = −0.998 and s = 9.13 kcal/mol for the RMDFT functional, and R = −0.998 and s = 8.05 kcal/mol for the KH + CT functional. Both CTs have strong correlation with PMV. (Those are negative correlations, but negative sign is only a matter of how a coefficient is taken.) This shows that both CTs have the same effect as the PMV correction. The inclination coefficients also corresponded to the parameter a of the PMV correction. In the RMDFT functional, the diameter of the reference HS plays the role of this parameter. In short, the CT of the RMDFT functional provides the contribution of cavity formation by the solute depending on the size of the reference HS. Fig. 6 shows comparisons of the calculation results and the experimental data for the SFE values of 504 small organic molecules in water [68, 71]. Both SFE results agree well with the experimental data, especially for hydrophobic solutes that have positive or small SFE values. However, there are differences in the degree of variation between the RMDFT and KH + CT functionals. The MADs of the RMDFT and KH + CT functionals from experimental data were 1.12 and 1.22 kcal/mol, respectively. Similarly, the RMSDs for RMDFT and KH + CT functionals were 1.47 and 2.85 kcal/mol, respectively. Although the effects of both CTs are similar, the result of the RMDFT functional was better than that of the KH + CT functional. In contrast, the correlation coefficient and the standard deviation of RMDFT results from experimental data were R = 0.907 and s = 1.48 kcal/mol, and R = 0.912 and s = 1.61 kcal/mol for KH + CT functional. The RMDFT is somewhat worse from the correlation coefficient viewpoint. Fig. 7 shows the correlation between the SFE values of the RMDFT and the KH + CT functionals for the 504 small organic molecules. The linear regression was y = 1.08x − 0.191 and the inclination exceeded 1, in contrast to the results shown in Fig. 3. However, the reciprocal of the slope of the correlation between the DFT(D) functional and the KH functional was 1.10(= 1/0.909), which was the same degree of deviation. This is essentially the difference between the HNC functional (which is equivalent to DFT(D)) and KH functional. The difference between the HNC functional and
6
Y. Maruyama / Journal of Molecular Liquids 291 (2019) 111160
a
b
Fig. 4. Correlation between the dimensionless partial molar volume and differences between the experimental data and the calculation values using the DFT(D) functional (left) and the KH functional (right).
the KH functional is only the first term of the integral (see Eqs. (15) and (16)). In other words, the deviation of the two SingerChandler formulas is caused by this term. In contrast to the inclination, we observed R = 0.998, and s = 0.351 kcal/mol, indicating a strong correlation.
4. Conclusion We divided the RMDFT functional into the DFT(D) and CT parts and investigated their contributions. We also investigated its usefulness by combining the KH functional and the CT to confirm consistency with the closure equation. The CT estimates the correction of the cavity formation of the solute by the reference HS. The contribution of the CT increases with the diameter of the HS of the reference system, and it strongly correlated with
a
the PMV. As the result, the CT based on the RMDFT showed the same effect as PMV correction, although the RMDFT functional does not contain the parameter b of the PMV correction. On the other hand, the KH + CT functional showed an accuracy comparable to that of the RMDFT functional. The result shows that consistency between the SFE functional and the closure equation is unnecessary. SFE calculations and application by the 3D-RISM theory is still under development. Recently, Kobryn et al. derived a new KGK closure for the 3D-RISM theory that is particularly suitable for polar and charged molecules in electrolyte solutions. The KGK closure levels out the density distribution function to zero inside the repulsive core to avoid the problem of the HNC-type approximation [72]. Sosnin et al. demonstrated that the bioconcentration factor can be effectively predicted by the combination of the 3D-RISM theory and machine learning [73].
b
Fig. 5. Correlation between the dimensionless partial molar volume and the values of correction terms for the DFT(D) functional (left) and for the KH functional (right).
Y. Maruyama / Journal of Molecular Liquids 291 (2019) 111160
a
7
b
Fig. 6. Comparison of the solvation free energy values of 504 small organic molecules obtained using the RMDFT (left) and the KH + CT (right) functionals and the experimentally determined values.
Improvement of 3D distribution function is essential to improve the accuracy of SFE of hydrophilic molecules because the contribution of the Coulomb interaction between solute and solvent is large. One way to do this is to introduce a bridge correction term in the HNC-type closure equation to overcome its softness. Such a term is obtained as a functional derivative of the SFE correction term [74, 75]. That is, the bridge correction term for the closure equation is in correspondence with the correction term of the SFE functional [76, 77]. In fact, we derived the HS bridge correction term based on RMDFT and applied it to Lennard-Jones fluids [47]. In addition, we can systematically improve the SFE functional and the closure equation by using RMDFT. However, it is currently difficult to converge 3D-RISM calculation with the RMDFT closure equation including the HS bridge correction term.
We may need to develop a new convergence scheme or convergence acceleration method for the 3D-RISM calculation with the RMDFT closure equation.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27]
Fig. 7. Correlation between the solvation free energy values of the RMDFT and KH + CT functionals for 504 small organic molecules.
[28] [29]
D. Beglov, B. Roux, J. Phys. Chem. B 101 (1997) 7821–7826. A. Kovalenko, F. Hirata, Chem. Phys. Lett. 290 (1998) 237–244. A. Kovalenko, F. Hirata, J. Phys. Chem. B 103 (1999) 7942–7957. S. Gusarov, B.S.B. Pujari, A. Kovalenko, J. Comput. Chem. 33 (2012) 1478–1494. V.P. Sergiievskyi, M.V. Fedorov, J. Chem. Theory Comput. 8 (2012) 2062–2070. Y. Maruyama, F. Hirata, J. Chem. Theory Comput. 8 (2012) 3015–3021. Y. Maruyama, N. Yoshida, H. Tadano, D. Takahashi, M. Sato, F. Hirata, J. Comput. Chem. 35 (2014) 1347–1355. D. Casanova, S. Gusarov, A. Kovalenko, T. Ziegler, J. Chem. Theory Comput. 3 (2007) 458–476. T. Luchko, S. Gusarov, D.R. Roe, C. Simmerling, D.A. Case, J. Tuszynski, A. Kovalenko, J. Chem. Theory Comput. 6 (2010) 607–624. N. Yoshida, S. Phongphanphanee, Y. Maruyama, T. Imai, F. Hirata, J. Am. Chem. Soc. 128 (2006) 12042–12043. N. Yoshida, S. Phongphanphanee, F. Hirata, J. Phys. Chem. B 111 (2007) 4588–4595. T. Yamazaki, T. Imai, F. Hirata, A. Kovalenko, J. Phys. Chem. B 111 (2007) 1206–1212. S. Phongphanphanee, N. Yoshida, F. Hirata, J. Am. Chem. Soc. 130 (2008) 1540–1541. S. Phongphanphanee, N. Yoshida, F. Hirata, J. Mol. Liq. 147 (2009) 107–111. S. Phongphanphanee, N. Yoshida, F. Hirata, J. Phys. Chem. B 114 (2010) 7967–7973. R.S. Johnson, T. Yamazaki, A. Kovalenko, H. Fenniri, J. Am. Chem. Soc. 129 (2007) 5735–5743. G. Tikhomirov, T. Yamazaki, A. Kovalenko, H. Fenniri, Langmuir 24 (2008) 4447–4450. T. Yamazaki, H. Fenniri, A. Kovalenko, ChemPhysChem 11 (2010) 361–367. G. Borzsonyi, R.S. Johnson, A.J. Myles, J.-Y. Cho, T. Yamazaki, R.L. Beingessner, A. Kovalenko, H. Fenniri, Chem. Commun. 46 (2010) 6527–6529. T. Yamazaki, H. Fenniri, J. Mol. Liq. 217 (2016) 70–74. T. Yamazaki, N. Blinov, D. Wishart, A. Kovalenko, Biophys. J. 95 (2008) 4540–4548. S.-H. Chong, S. Ham, Proc. Natl. Acad. Sci. 109 (2012) 7636–7641. S.-H. Chong, S. Ham, Angew. Chem. 126 (2014) 1–5. D.J. Sindhikara, N. Yoshida, M. Kataoka, F. Hirata, J. Mol. Liq. 164 (2011) 120–122. Y. Maruyama, N. Yoshida, F. Hirata, J. Phys. Chem. B 114 (2010) 6464–6471. G.M. Giambasu, ¸ M.K. Gebala, M.T. Panteva, T. Luchko, D.A. Case, D.M. York, Nucleic Acids Res. 43 (2015) 8405–8415. Y. Maruyama, T. Matsushita, R. Ueoka, F. Hirata, J. Phys. Chem. B 115 (2011) 2408–2416. J. Lee, H. Im, S.-H. Chong, S. Ham, Chem. Phys. Lett. 693 (2018) 216–221. S. Singer, D. Chandler, Mol. Phys. 55 (1985) 621–625.
8
Y. Maruyama / Journal of Molecular Liquids 291 (2019) 111160 [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53]
A. Kovalenko, F. Hirata, J. Chem. Phys. 110 (1999) 10095–10112. S. Ten-No, J. Chem. Phys. 115 (2001) 3724. K. Sato, H. Chuman, S. Ten-No, J. Phys. Chem. B 109 (2005) 17290–17295. D.S. Palmer, A.I. Frolov, E.L. Ratkova, M.V. Fedorov, J. Phys. Condens. Matter 22 (2010) 492101. D.S. Palmer, A.I. Frolov, E.L. Ratkova, M.V. Fedorov, Mol. Pharm. 8 (2011) 1423–1429. A.I. Frolov, E.L. Ratkova, D.S. Palmer, M.V. Fedorov, J. Phys. Chem. B 115 (2011) 6011–6022. J.-F. Truchon, B.M. Pettitt, P. Labute, J. Chem. Theory Comput. 10 (2014) 934–941. V.P. Sergiievskyi, G. Jeanmairet, M. Levesque, D. Borgis, J. Phys. Chem. Lett. 5 (2014) 1935–1942. V.P. Sergiievskyi, G. Jeanmairet, M. Levesque, D. Borgis, J. Chem. Phys. 143 (2015) 184116–184117. M. Misin, D.S. Palmer, M.V. Fedorov, J. Phys. Chem. B 120 (2016) 5724–5731. M. Kinoshita, Y. Harano, R. Akiyama, J. Chem. Phys. 125 (2006) 244504–244507. S.-H. Chong, S. Ham, Chem. Phys. Lett. 535 (2012) 152–156. S.-H. Chong, S. Ham, J. Chem. Theory Comput. 11 (2015) 378–380. T. Miyata, J. Thapa, Chem. Phys. Lett. 604 (2014) 122–126. T. Miyata, Y. Ebato, J. Mol. Liq. 217 (2016) 75–82. T. Sumi, A. Mitsutake, Y. Maruyama, J. Comput. Chem. 36 (2015) 1359–1369. T. Sumi, A. Mitsutake, Y. Maruyama, J. Comput. Chem. (2015) 2009–2011. T. Sumi, Y. Maruyama, A. Mitsutake, K. Koga, J. Chem. Phys. 144 (2016) 224104–224116. Y. Maruyama, A. Mitsutake, J. Phys. Chem. B 121 (2017) 9881–9885. T. Sumi, Y. Maruyama, A. Mitsutake, K. Mochizuki, K. Koga, J. Comput. Chem. 39 (2018) 202–217. J.P. Hansen, I.R. McDonald, Theory of Simple Liquids, Third ed. ed., Elsevier/Academic Press, London, 2006. D. Henderson (Ed.), Fundamentals of Inhomogeneous Fluids, Marcel Dekker, New York., 1992, A. Kovalenko, F. Hirata, J. Chem. Phys. 112 (2000) 10391–10402. D. Chandler, J. Mccoy, S. Singer, J. Chem. Phys. 85 (1986) 5971–5976.
[54] Y. Rosenfeld, Phys. Rev. Lett. 63 (1989) 980–983. [55] R. Roth, R. Evans, A. Lang, G. Kahl, J. Phys. Condens. Matter 14 (2002) 12063–12078. [56] Y.-X. Yu, J. Wu, J. Chem. Phys. 117 (2002) 10156–10164. [57] T. Sumi, H. Sekino, J. Phys. Soc. Jpn. 77 (2008) 034605. [58] N.F. Carnahan, K.E. Starling, J. Chem. Phys. 51 (1969) 635–636. [59] A. Ben-Naim, Molecular Theory of Solutions, Oxford University Press, New York, 2006. [60] D. Chandler, Y. Singh, D.M. Richardson, J. Chem. Phys. 81 (1984) 1975. [61] D. Chandler, T. Ichiye, J. Phys. Chem. 92 (1988) 5257. [62] P.H. Lee, G.M. Maggiora, J. Phys. Chem. 97 (1993) 10175. [63] T. Imai, M. Kinoshita, F. Hirata, J. Chem. Phys. 112 (2000) 9469–9478. [64] T. Imai, Y. Harano, A. Kovalenko, F. Hirata, Biopolymers 59 (2001) 512–519. [65] M.R. Shirts, J.W. Pitera, W.C. Swope, V.S. Pande, J. Chem. Phys. 119 (2003) 5740–5761. [66] Y. Karino, M.V. Fedorov, N. Matubayasi, Chem. Phys. Lett. 496 (2010) 351–355. [67] Y. Liu, J. Fu, J. Wu, J. Phys. Chem. Lett. 4 (2013) 3687–3691. [68] D.L. Mobley, C.I. Bayly, M.D. Cooper, M.R. Shirts, K.A. Dill, J. Chem. Theory Comput. 5 (2009) 350–358. [69] W. Jorgensen, J. Chandrasekhar, J. Madura, R. Impey, M. Klein, J. Chem. Phys. 79 (1983) 926–935. [70] B. Pettitt, P. Rossky, J. Chem. Phys. 77 (1982) 1451–1457. [71] R.C. Rizzo, T. Aynechi, D.A. Case, I.D. Kuntz, J. Chem. Theory Comput. 2 (2006) 128–139. [72] A.E. Kobryn, S. Gusarov, A. Kovalenko, J. Phys. Condens. Matter 28 (2016) 404003–404024. [73] S. Sosnin, M. Misin, D.S. Palmer, M.V. Fedorov, J. Phys. Condens. Matter 30 (2018) 32LT03. [74] Y. Liu, S. Zhao, J. Wu, J. Chem. Theory Comput. 9 (2013) 1896–1908. [75] Y. Liu, J. Wu, J. Chem. Phys. 140 (2014) 084103–084108. [76] A. Kovalenko, F. Hirata, J. Chem. Phys. 113 (2000) 2793–2805. [77] D. Yokogawa, J. Chem. Theory Comput. 14 (2018) 3272–3278. [78] R. Wolfenden, L. Andersson, P.M. Cullis, C.C. Southgate, Biochemistry 20 (1981) 849–855.