We extend the correlation functions obtained by molecular dynamics (MD) simulation for a molten salt modeled as a superposition of the Lennard-Jones (LJ) and Coulomb potentials using the hybrid closure method, which employs the Ornstein–Zernike (OZ) theory coupled with a closure relation. An appropriate distance for switching the short-range MD part and the long-range OZ part is determined by monitoring the isothermal compressibility, excess internal energy, and pressure. The Kobryn–Gusarov–Kovalenko (KGK) closure relation is mainly employed for the hybrid closure method (MD–KGK hybrid closure). The hybrid closure with either the hypernetted chain (HNC) or Kovalenko–Hirata (KH) closure was also tested to confirm that the performance was almost equivalent to one another among the MD–HNC, MD–KH, and MD–KGK methods. The bridge function for the model molten salt is extracted using the MD–KGK hybrid closure method. At a high-density state, the bridge function shows a steep increase in the repulsive core region, as is often observed for simple fluids, whereas when the density is relatively low, the bridge function for the cation–anion pair shows a downward-sloping behavior. Furthermore, the accuracies of excess internal energy, pressure, and isothermal compressibility were also examined for the HNC, KH, and KGK approximations. For molten salt systems, these approximations exhibited a similar behavior to those for monatomic LJ fluids, especially in the high-density state. The analysis of the integrand for excess internal energy and pressure is also discussed.

The integral equation theory (IET) is one of the most powerful methods to describe the structure and thermodynamics of solutions.1 The application of the IETs to solution systems has been extended even to the hydration of proteins and molecular recognition.2–7 In many cases, these applications can be looked at in terms of the solvation of the solute. In the discussion, the distribution function of solvents around the solute, solvation free energy, partial molar volume, and so on play a central role.1 When the IETs are employed for these calculations, an infinite dilution limit is usually assumed. In this case, first, the IET for evaluating the correlation function between the solvents (or the solvent susceptibility response function) is solved, followed by solving IET between the solute and the solvent. The solvent susceptibility response function can be regarded as input data used when solving the IET between the solute and the solvent.1 

The accuracy of the calculation of solvation by the IETs depends on the approximation employed.8–14 The more accurate the approximation employed is, the more accurate the obtained result is. The calculation of the solvation is quantitatively affected not only by the accuracy of the approximation employed for the IET calculation between the solute and solvent but also by the quality of the solvent susceptibility response function.15,16 As a matter of course, using an accurate solvent susceptibility response function leads to an accurate calculation of the solvation. This study focuses on a method to accurately obtain a solvent susceptibility response function.

The improvement of the accuracy of the IETs is usually performed by bridge correction.17 This is based on the idea that by using an accurate bridge function, numerical results obtained from IETs become accurate. From this viewpoint, various approximated bridge functions have been proposed to date (for example, Refs. 10, and 1825). On the other hand, in the calculation of the solvation, a complicated shape is quite typically seen on the solute molecules rather than the solvent ones. In fact, the solvent molecule is not so complicated in shape in many cases. It would be possible to include even water in a typical example of such a solvent molecule. For a solvent molecule whose shape is relatively simple, molecular dynamics (MD) simulation can be used to evaluate the solvent susceptibility response function. For systems interacting through a spherically symmetric potential such as a monatomic Lennard-Jones (LJ) fluid, for example, a method has been proposed that extends the radial distribution function (RDF) obtained from an MD simulation to a longer-region with the help of the IETs.26–30 Hereafter, the method for this extension is referred to as the hybrid closure method. Our group has also studied the hybrid closure method for LJ fluids, taking the application of the method to the solvation in mind.15,16

When the IETs are applied to systems consisting of complicated polyatomic molecules such as proteins, a versatile potential model is often used [for example, Assisted Model Building with Energy Refinement (AMBER),31 Chemistry at Harvard Molecular Mechanics (CHARMM),32 Optimized Potentials for Liquid Simulations (OPLS),33 and Groningen Molecular Simulation (GROMOS)34]. These versatile potential models are expressed by a superposition of the LJ and Coulomb potentials. Hereafter, such a model is referred to as a “Coulomb system.” In order to apply the IETs to a realistic model, it is necessary to evaluate an accurate solvent susceptibility response function for the Coulomb system.

The objective of this study is to accurately evaluate the solvent susceptibility response function for the Coulomb system by using the above-mentioned hybrid closure method. In this study, we employ a modeled molten salt as the simplest Coulomb system. This article discusses the hybrid closure method, which is a method to extend the RDF of modeled molten salt obtained by an MD simulation to a long-ranged space by incorporating it into the closure relation of the IET. While the hybrid closure method has been studied for single-component systems,15,16,26–30 this study considers a two-component system. Furthermore, we also discuss the accuracy of the hypernetted chain (HNC), Kovalenko–Hirata (KH), and Kobryn–Gusarov–Kovalenko (KGK) closures, which are existing closure approximations, when they are applied to the modeled molten salt.

In this study, we employ a simple model of molten salts. The interaction potential we assume is the superposition of the LJ and Coulomb potentials. It is written as
uαγr=4εαγσαγr12σαγr6+qαqγ4πε0r,
(1)
where σ, ɛ, and q are the size parameter, the energy parameter, and the charge, respectively. The symbol ɛ0 denotes the permittivity of a vacuum. The subscripts α and γ are the indices to identify the species, i.e., either the cation or the anion. Hereafter, the cation and anion are denoted by “C” and “A,” respectively. Although the parameter σαγ is the mixed LJ parameter, it is practically identical both to σC and σA because we assume σC = σA in this study. Here, σC and σA indicate the size of the cation and the anion, respectively. Similarly, the mixed LJ parameter, ɛαγ, is practically identical both to ɛC and ɛA because ɛC = ɛA is assumed (see Sec. II D for details).
This study focuses on the extension of the RDF obtained from MD simulations with the aid of Ornstein–Zernike theory, which is one of the IETs. Here, we describe the OZ and closure equations for systems consisting of cations and anions. For such two-component systems, the OZ equation reads
ĥαγk=ĉαγk+μ=C,Aĉαμkρμĥμγ(k),
(2)
where h, c, and ρ represent the total correlation function, direct correlation function, and density, respectively. The caret indicates the Fourier transform. The wave number is represented by k.
The closure equations to extend the RDFs of MD simulations are as follows. Mainly, this study employs the hybrid closure between the MD simulation and the KGK closure equation.35 It is written as
1+hαγr=gαγMD(r)forrr00forr>r0anddαγr11+dαγrforr>r0anddαγr>1,
(3)
where
dαγr=βuαγr+hαγrcαγr.
(4)
In Eq. (4), β=kBT1, where kB and T are the Boltzmann constant and temperature, respectively. Hereafter, we refer to this hybrid closure as the MD–KGK hybrid closure. In Eq. (3), the term gαγMD(r) represents the RDF obtained from the MD simulations, which is given by
gαγMDr=VΔNγ4πr2ΔrNγ.
(5)
Equation (5) assumes that the molecule α is located at the origin, and r means the separation from the origin (i.e., molecule α in this case). The symbols V, Nγ, and ΔNγ represent the system volume, the number of molecules γ in the system, and the number of molecules γ contained in the spherical shell of the radial width of Δr. The angle bracket in Eq. (5) indicates the ensemble average. In the MD–KGK hybrid closure method, Eqs. (2) and (3) are solved simultaneously. In this study, we basically use the MD–KGK hybrid closure unless otherwise stated.
Other options are the hybrid between the MD simulation and the HNC closure equation1,17 and that between the MD and the KH closure equation.1,36 The former is referred to as the MD–HNC hybrid closure, and the latter as the MD–KH hybrid closure. The MD–HNC closure is written as
1+hαγr=gαγMD(r)forrr0expdαγrforr>r0,
(6)
where as the MD–KH closure is given by
1+hαγr=gαγMD(r)forrr0expdαγrforr>r0anddαγr01+dαγrforr>r0anddαγr>0.
(7)
The RDF gαγr is related to hαγr as
gαγr=1+hαγr.
(8)
The excess internal energy per molecule uex was evaluated from
uex=2πρα=C,Aγ=C,Axαxγ0drr2uαγ(r)gαγr,
(9)
where ρ = ρC + ρA, xC = ρC/ρ, and xA = ρA/ρ. The pressure p was obtained via
p=ρkBT23πα=C,Aγ=C,Aραργ0drr3duαγrdrgαγr.
(10)
To estimate the isothermal compressibility, the Kirkwood–Buff theory37 was employed. It is written as
κT=ζkBTη,
(11)
where η and ζ for the current system are expressed by
η=ρC+ρA+ρCρA(GCC+GAA2GCA),
(12)
and
ζ=1+ρCGCC+ρAGAA+ρCρA(GCCGAAGCA2),
(13)
respectively. The term Gαγ is defined by
Gαγ=4π0drr2gαγr1.
(14)
In terms of thermodynamic accuracy, the HNC, KH, and KGK closures are also examined for comparison in this study. Here, these closure equations are also described for completeness. The HNC closure reads1,17
1+hαγr=expdαγr.
(15)
The KH1,36 and KGK35 closures are
1+hαγr=expdαγrdαγr0,1+dαγrdαγr>0,
(16)
and
1+hαγr=0,dαγr1,1+dαγr,dαγr>1,
(17)
respectively. Solving these closure equations together with Eq. (2), the numerical solutions under the bare HNC, bare KH, and bare KGK approximations are obtained.

In this study, σC = σA = 3.405 Å and ɛC = ɛA = 0.2380 kcal/mol were assumed both for the monatomic cation “C” and anion “A” throughout this study. Although the Lorentz–Berthelot mixing rule was applied to evaluate the mixed LJ parameters σαγ and ɛαγ, they are practically identical to σC = σA and ɛC = ɛA, as mentioned in Sec. II A. Four different charges were employed for qC = −qA, which were 0.3e, 0.6e, 0.8e, and 1.0e. The densities assumed in this study are ρCσC3=ρAσA3=0.05, 0.15, 0.25, 0.35, and 0.45 (this means ρC = ρA since σC = σA is assumed). The temperature considered herein is as follows: for the model of qC = −qA = 0.3e, 0.6e, 0.8e, and 1.0e, kBT/ɛC = kBT/ɛA = 10.0, 15.0, 30.0, and 45.0, respectively. For later convenience, we introduce the notation ρC*, ρA*, and T* as ρC*=ρCσC3, ρA*=ρAσA3, and T* = kBT/ɛC = kBT/ɛA, respectively. The details of the model parameters and thermodynamic conditions are summarized in Table I. Hereafter, the models of qC = −qA = 0.3e, 0.6e, 0.8e, and 1.0e are referred to as the model types 1, 2, 3, and 4, respectively.

TABLE I.

The model parameters and thermodynamic conditions employed in this study. The symbols ρC*, ρA*, and T* are defined as ρC*=ρCσC3, ρA*=ρAσA3, and T* = kBT/ɛC = kBT/ɛA, respectively.

ModelσC(=σA) (Å)ɛC(=ɛA) (kcal/mol)qC (e)qA (e)ρC*(=ρA*)T*
Type 1 3.405 0.2380 0.3 −0.3 0.05–0.45 10.0 
Type 2 3.405 0.2380 0.6 −0.6 0.05–0.45 15.0 
Type 3 3.405 0.2380 0.8 −0.8 0.05–0.45 30.0 
Type 4 3.405 0.2380 1.0 −1.0 0.05–0.45 45.0 
ModelσC(=σA) (Å)ɛC(=ɛA) (kcal/mol)qC (e)qA (e)ρC*(=ρA*)T*
Type 1 3.405 0.2380 0.3 −0.3 0.05–0.45 10.0 
Type 2 3.405 0.2380 0.6 −0.6 0.05–0.45 15.0 
Type 3 3.405 0.2380 0.8 −0.8 0.05–0.45 30.0 
Type 4 3.405 0.2380 1.0 −1.0 0.05–0.45 45.0 

The grid spacing and the number of grids for the OZ calculation were 0.05 Å and 16384, respectively. The modified method of direct inversion in iterative subspace (MDIIS)38 was used for numerical convergence of the theories, where we set the numerical tolerance at 10−9. The offset of gαγr less than 2 × 10−6 was cutoff after a numerical convergence was achieved, which was typically observed inside the repulsive core.

Concerning the MD simulations, canonical (NVT) ensemble simulations were performed to obtain the RDF, excess internal energy uex, and pressure p. In a cubic cell, 2662 cations and the same number of anions were placed, assuming a periodic boundary condition. A timestep of 1 fs was used. After the system was equilibrated for 10 ps, a production run of 4 ns was performed. In the NVT ensemble MD simulations, the temperature was controlled by a Nose–Hoover chain thermostat.39,40 Isothermal–isobaric (NPT) ensemble MD simulations were also performed to evaluate the isothermal compressibility κT, where the length of equilibration was 10 ps and that of the production run was 3 ns. The Nose–Hoover chain thermostat and Andersen barostat were applied to the NPT ensemble MD simulation,39,40 where the pressure p obtained from the NVT ensemble MD simulation was assumed. In these MD simulations, the LJ potential was truncated at 7σC, while the Coulomb potential was handled by the Ewald method.

The MD–KGK hybrid closure involves the parameter r0 as in Eq. (3). Other hybrid closures [Eqs. (6) and (7)] also involve it. On one hand, the too small value of r0 in Eq. (3) gives rise to just a KGK closure almost everywhere. On the other hand, computation using a too large value of r0 would be significantly contaminated by a statistical error in the RDF obtained from the MD simulation. Therefore, it is necessary to optimize the r0 value. In the previous studies on monatomic LJ fluids,15,16 this optimization was mainly judged in terms of the r0-dependence of the isothermal compressibility κT. In this study, we add two more criteria: the r0-dependence of the excess internal energy uex and that of the pressure p.

Examples of the dependence of κT on r0 are shown in Fig. 1. Figure 1(a) corresponds to the results obtained for model type 1 at ρC*=ρA*=0.45 and T* = 10.0, whereas (b) for model type 3 at ρC*=ρA*=0.45 and T* = 30.0. The data on the isothermal compressibility obtained from the NPT–MD simulation are also shown in Fig. 1. Similar data to Fig. 1 obtained at other conditions can be found in Fig. S1 in the supplementary material. In Figs. 2 and 3, the r0-dependence of uex and that of p, respectively, are shown. The panels (a) in Figs. 2 and 3 show the results for model type 1 at ρC*=ρA*=0.45 and T* = 10.0, while the panels (b) correspond to those for model type 3 at ρC*=ρA*=0.45 and T* = 30.0. The results of uex and p obtained from the NVT–MD simulation are also shown in Figs. 2 and 3, respectively. Similar data to Figs. 2 and 3 obtained at other conditions can be found in Figs. S2 and S3 in the supplementary material.

FIG. 1.

Dependence of κT on r0: (a) for model type 1 at ρC*=ρA*=0.45 and T* = 10.0; (b) for model type 3 at ρC*=ρA*=0.45 and T* = 30.0. The data on the isothermal compressibility obtained from the NPT–MD simulation are also shown.

FIG. 1.

Dependence of κT on r0: (a) for model type 1 at ρC*=ρA*=0.45 and T* = 10.0; (b) for model type 3 at ρC*=ρA*=0.45 and T* = 30.0. The data on the isothermal compressibility obtained from the NPT–MD simulation are also shown.

Close modal
FIG. 2.

Dependence of uex on r0: (a) for model type 1 at ρC*=ρA*=0.45 and T* = 10.0; (b) for model type 3 at ρC*=ρA*=0.45 and T* = 30.0. The data on the excess internal energy obtained from the NVT–MD simulation are also shown.

FIG. 2.

Dependence of uex on r0: (a) for model type 1 at ρC*=ρA*=0.45 and T* = 10.0; (b) for model type 3 at ρC*=ρA*=0.45 and T* = 30.0. The data on the excess internal energy obtained from the NVT–MD simulation are also shown.

Close modal
FIG. 3.

Dependence of p on r0: (a) for model type 1 at ρC*=ρA*=0.45 and T* = 10.0; (b) for model type 3 at ρC*=ρA*=0.45 and T* = 30.0. The data on the pressure obtained from the NVT–MD simulation are also shown.

FIG. 3.

Dependence of p on r0: (a) for model type 1 at ρC*=ρA*=0.45 and T* = 10.0; (b) for model type 3 at ρC*=ρA*=0.45 and T* = 30.0. The data on the pressure obtained from the NVT–MD simulation are also shown.

Close modal

The isothermal compressibility κT of the hybrid MD–KGK closure is almost constant when r0/σC < 0.7 in Fig. 1, since in this region the approximation is nearly equivalent to the bare KGK closure. In the region of 0.8 < r0/σC < 2.0, the curves of the hybrid MD–KGK closure in Fig. 1 exhibit some oscillations. When the r0/σC value becomes larger than 2.0, the curves in Fig. 1 show a plateau again. In the region of r0/σC > 3.0, the curve deviates from this plateau because of the statistical error involved in the data of MD simulations. In the plateau of 2.0 < r0/σC < 3.0, the κT values of MD–KGK closure basically become close to those of the NPT-MD simulations, though there is such a tendency that the κT value of MD–KGK closure is slightly smaller than that of the NPT-MD simulation (see also Fig. S1).

Concerning the excess internal energy uex (Fig. 2) and pressure p (Fig. 3), the results of hybrid MD–KGK closure generally deviate from those of the NVT-MD simulations in the region of r0/σC < 0.7. This indicates that the values of uex and p obtained from the bare KGK closure deviate from the MD simulation values (the accuracy of thermodynamics based on bare HNC, bare KH, and bare KGK closures are discussed in Sec. III C). When the values of r0/σC are greater than 1.5, both uex and p values become almost constant irrespective of r0/σC. The uex and p values in the region of r0/σC > 1.5 are almost coincident with those of the NVT-MD simulations, though the former is slightly smaller than the latter. This tendency is also observed at other conditions (Figs. S2 and S3). The underestimation of thermodynamic quantities in the hybrid closure calculation may partly be attributed to the relatively small system size in the MD simulation (i.e., the thermodynamics obtained from NPT-MD and/or NVT-MD simulations may have been somewhat contaminated by a system size effect).

From the results shown in Figs. 13, it is possible to conclude that an appropriate value of r0 should be determined in the region of 2.0 < r0/σC < 3.0. Moreover, we examine the difference among the MD–HNC, MD–KH, and MD–KGK hybrid closures. A comparison among these approximations was made in Fig. 4, where Figs. 4(a)4(c) show the r0-dependence of κT, that of uex, and that of p, respectively, obtained for model type 1 at ρC*=ρA*=0.45 and T* = 10.0. Similar data to Fig. 4 obtained at other conditions can be found in Figs. S4–S6 in the supplementary material. By definition, the results in the region of r0/σC < 0.7 are nearly equivalent to those obtained by one of the bare HNC, KH, and KGK closures. Reflecting this fact, the curves in Fig. 4 in this region deviate from one another. On the other hand, the curves in the region of r0/σC > 2.0 coincide with one another. This indicates that, as far as the numerical value of r0 is given appropriately, the MD–HNC, MD–KH, and MD–KGK hybrid closures are almost equivalent to one another.

FIG. 4.

Comparison of r0-dependences of thermodynamic quantities among MD–HNC, MD–KH, and MD–KGK hybrid closures for model type 1 at ρC*=ρA*=0.45 and T* = 10.0: (a) r0-dependences of κT, (b) r0-dependences of uex, and (c) r0-dependences of p.

FIG. 4.

Comparison of r0-dependences of thermodynamic quantities among MD–HNC, MD–KH, and MD–KGK hybrid closures for model type 1 at ρC*=ρA*=0.45 and T* = 10.0: (a) r0-dependences of κT, (b) r0-dependences of uex, and (c) r0-dependences of p.

Close modal

The differences among the three hybrid closures (MD–HNC, MD–KH, and MD–KGK) were also examined in terms of the direct correlation functions. In evaluating the direct correlation function, we chose the value of r0/σC = 2.5. Examples of cCCr and cCAr are shown in Fig. 5, which was obtained for model type 1 at ρC*=ρA*=0.45 and T* = 10.0. Similar data to Fig. 5 obtained at other conditions can be found in Fig. S7 in the supplementary material. In both curves of cCCr and cCAr, practically no differences are observed among the three hybrid closures. Similar coincidences among these three hybrid closures were also obtained for monatomic LJ fluids.16 Therefore, only the MD–KGK closure is employed for the hybrid closure method hereafter in this study.

FIG. 5.

Comparison of direct correlation functions cCCr and cCAr among MD–HNC, MD–KH, and MD–KGK hybrid closures for model type 1 at ρC*=ρA*=0.45 and T* = 10.0.

FIG. 5.

Comparison of direct correlation functions cCCr and cCAr among MD–HNC, MD–KH, and MD–KGK hybrid closures for model type 1 at ρC*=ρA*=0.45 and T* = 10.0.

Close modal
In the remaining part of this article, we employ r0/σC = 2.5 for the hybrid closure method. Using the MD–KGK hybrid closure, we extracted the bridge function for the current model of Coulomb systems. The bridge function is estimated from
bαγr=lngαγr+βuαγrhαγr+cαγr.
(18)
Examples of the obtained bridge functions are shown in Fig. 6. Figure 6(a) shows both bCCr and bCAr obtained for model type 1 at ρC*=ρA*=0.05 and T* = 10.0, whereas (b) those for model type 1 at ρC*=ρA*=0.45 and T* = 10.0. Figures 6(c) and 6(d) correspond to the results for model type 3 at ρC*=ρA*=0.05 and T* = 30.0 and to those at ρC*=ρA*=0.45 and T* = 30.0, respectively.
FIG. 6.

Bridge functions bCCr and bCAr extracted by the MD-KGK hybrid closure: (a) for model type 1 at ρC*=ρA*=0.05 and T* = 10.0; (b) for model type 1 at ρC*=ρA*=0.45 and T* = 10.0; (c) for model type 3 at ρC*=ρA*=0.05 and T* = 30.0; (d) for model type 3 at ρC*=ρA*=0.45 and T* = 30.0. The bridge function of a single-component LJ fluid is also shown, where bLJr is the one extracted by the MD–KGK hybrid closure method and bVMr is the Verlet modified (VM) bridge function.20 The bridge function of an LJ fluid was evaluated by assuming the LJ fluid density as ρ = ρC + ρA (see the text for details).

FIG. 6.

Bridge functions bCCr and bCAr extracted by the MD-KGK hybrid closure: (a) for model type 1 at ρC*=ρA*=0.05 and T* = 10.0; (b) for model type 1 at ρC*=ρA*=0.45 and T* = 10.0; (c) for model type 3 at ρC*=ρA*=0.05 and T* = 30.0; (d) for model type 3 at ρC*=ρA*=0.45 and T* = 30.0. The bridge function of a single-component LJ fluid is also shown, where bLJr is the one extracted by the MD–KGK hybrid closure method and bVMr is the Verlet modified (VM) bridge function.20 The bridge function of an LJ fluid was evaluated by assuming the LJ fluid density as ρ = ρC + ρA (see the text for details).

Close modal

For comparison, the bridge function of a single-component LJ fluid is also shown in Fig. 6, which was obtained as follows. An MD simulation was executed to obtain an RDF of a single-component LJ fluid, where the density of the LJ fluid was assumed to be ρ = ρC + ρA. For example, in the cases of Figs. 6(a) and 6(b), the LJ fluid densities were assumed to be ρσ3 = 0.10 and ρσ3 = 0.90, respectively. Then, using the obtained RDF of the MD simulation, the MD–KGK hybrid method was applied to obtain the total and direct correlation functions. Finally, using a similar equation to Eq. (18), the bridge function of the LJ fluid was estimated. Therefore, the obtained bridge function is denoted by bLJr. For further comparison, the bridge function bVMr of the Verlet modified (VM) closure approximation for LJ fluids20 was also plotted in Fig. 6, where the LJ fluid density was assumed to be ρ = ρC + ρA, as before. The formulas for the VM closure can also be found in Supplementary explanation 1 in the supplementary material. Similar data to Fig. 6 obtained at other conditions can be found in Fig. S8 in the supplementary material.

The curve of the bridge function for the cation–anion (CA) pair, bCAr, lies over that for the cation–cation pair, bCCr, in Fig. 6, irrespective of the model types and densities. The VM bridge function bVMr is almost coincident with bLJr. The curve of bLJr is located between those of bCCr and bCAr. These results indicate that the bridge function for the Coulomb systems is significantly different from that for the LJ system. When the density is relatively high, as in Figs. 6(b) and 6(d), both bCCr and bCAr show a steep increase with r in the region of r/σC < 1.0. Such a steep increase in the bridge function is known to be observed, for example, for LJ fluids.18,20 On the other hand, at a relatively low-density state, only bCCr shows an increasing behavior with r in the region of r/σC < 1.0. The value of bCAr in Fig. 6(c) decreases with increasing r when density is relatively low. Such a downward-sloping behavior as observed in the curve of bCAr had not been reported for LJ fluids.18,20 Even if the density of the LJ fluid made coincident with the total density of the Coulomb system, the bridge functions of the LJ fluids cannot fully describe those of the Coulomb systems. This tendency is prominent, especially for low-density states, i.e., the bridge function for the cation–anion pair is qualitatively different from that for the LJ fluids. When the density becomes relatively high, the bridge function of the Coulomb systems may be qualitatively described by that of the LJ systems in terms of a steep increase in the small r-region. Nevertheless, the bridge functions of the Coulomb systems and those of the LJ systems are quantitatively different from each other, even at high-density regions.

In this section, existing closures are examined in terms of the accuracy of thermodynamics. The excess internal energy uex, pressure p, and isothermal compressibility κT were evaluated using the bare HNC, KH, and KGK closure relations. Figures 7(a)7(c) show the dependences of uex, p, and κT, respectively, on the cation density obtained for model type 3 at T* = 30.0. For comparison, the results of the MD simulation and the MD–KGK hybrid closure are also shown in Fig. 7. Similar data to Fig. 7 obtained at other conditions can be found in Figs. S9–S11 in the supplementary material.

FIG. 7.

Dependences of (a) uex, (b) p, and (c) κT on the cation density for model type 3 at T* = 30.0. The data of the MD simulation and that of the MD–KGK hybrid closure are also shown.

FIG. 7.

Dependences of (a) uex, (b) p, and (c) κT on the cation density for model type 3 at T* = 30.0. The data of the MD simulation and that of the MD–KGK hybrid closure are also shown.

Close modal

The HNC closure provided relatively accurate uex values, though at high density regions, it tended to overestimate uex values [Fig. 7(a)] compared to those by the MD and MD–KGK methods (the MD results and the MD–KGK ones were coincident with each other). The KH closure overestimated uex at the density region we examined. The uex values obtained by the KGK closure in relatively low-density states were close to those obtained by the KH closure, whereas they approached those obtained by the HNC closure in relatively high-density regions. The fact that the values of uex in the high density region are overestimated by the HNC, KH, and KGK closures is qualitatively similar to the results for the LJ fluids.17–19,41

The pressure p obtained by the HNC closure was significantly overestimated in the high density region in comparison with the MD and/or MD–KGK results [Fig. 7(b)]. The overestimation was more prominent when the KH closure was used. The p values obtained by the KGK closure were underestimated in the intermediate-density region, whereas in the high-density region, they were somewhat overestimated. These results qualitatively agree with the behavior observed for LJ fluids.17–19,41

Concerning the isothermal compressibility κT, the result of the MD–KGK hybrid closure seemed to be slightly smaller than that of the MD simulation in some cases, but generally their agreement was excellent [Fig. 7(c)]. When the density was low, the HNC closure seemed to be relatively accurate for κT, whereas in the high-density region, κT values were overestimated by the HNC. This tendency is similar to that observed for LJ fluids.17–19,41 Both the KH and KGK closures underestimated the κT values when the density is relatively low. The isothermal compressibility obtained by these two closures became close to that obtained by the MD simulation when the density increased. Namely, the KH and KGK closures were relatively accurate in terms of the isothermal compressibility in the high-density region.

Here, we focus on the excess internal energy uex and pressure p to determine the sources of their errors. For this purpose, we examine the corresponding integrands. Figures 8 and 9 show examples of the differences in the integrands for uex and p, respectively, obtained for model type 3 at ρC*=ρA*=0.45 and T* = 30.0. The terms “diff. in integrand” at the right axes of ordinate in Figs. 8 and 9, which mean the differences in the integrands, are defined as 1/(σC2εC)r2uαγ(r)gαγOZrgαγMDr and 1/(σC2εC)r3duαγr/drgαγOZrgαγMDr, respectively, where gαγOZr and gαγMDr represent the RDFs obtained by the OZ and MD methods, respectively. By definition, if the “diff. in integrand” is positive, the difference contributes to the overestimation of uex or p, and vice versa. The panels (a) and (b) in Figs. 8 and 9 show the results under the HNC closure. The panels (c) and (d) correspond to the results of the KH closure, whereas (e) and (f) correspond to those of the KGK closure. Furthermore, panels (a), (c), and (e) were obtained for the cation–cation (CC) pair, while (b), (d), and (f) were obtained for the cation–anion (CA) pair. Figures 8 and 9 also show the RDFs, which correspond to the left axes of the ordinate. Similar data to Figs. 8 and 9 obtained at other conditions can be found in Figs. S12–S15 and S16–S19, respectively, in the supplementary material.

FIG. 8.

Differences of the integrands for uex, obtained for model type 3 at ρC*=ρA*=0.45 and T* = 30.0: (a) for cation–cation pair under the HNC approximation; (b) for cation–anion pair under the HNC approximation; (c) for cation–cation pair under the KH approximation; (d) for cation–anion pair under the KH approximation; (e) for cation–cation pair under the KGK approximation; (f) for cation–anion pair under the KGK approximation. The term “diff. in integrand” at the right axis of the ordinate is defined as 1/(σC2εC)r2uαγ(r)gαγOZrgαγMDr. The left axis of the ordinate corresponds to the RDFs.

FIG. 8.

Differences of the integrands for uex, obtained for model type 3 at ρC*=ρA*=0.45 and T* = 30.0: (a) for cation–cation pair under the HNC approximation; (b) for cation–anion pair under the HNC approximation; (c) for cation–cation pair under the KH approximation; (d) for cation–anion pair under the KH approximation; (e) for cation–cation pair under the KGK approximation; (f) for cation–anion pair under the KGK approximation. The term “diff. in integrand” at the right axis of the ordinate is defined as 1/(σC2εC)r2uαγ(r)gαγOZrgαγMDr. The left axis of the ordinate corresponds to the RDFs.

Close modal
FIG. 9.

Differences of the integrands for p, obtained for model type 3 at ρC*=ρA*=0.45 and T* = 30.0: (a) for cation–cation pair under the HNC approximation; (b) for cation–anion pair under the HNC approximation; (c) for cation–cation pair under the KH approximation; (d) for cation–anion pair under the KH approximation; (e) for cation–cation pair under the KGK approximation; (f) for cation–anion pair under the KGK approximation. The term “diff. in integrand” at the right axis of the ordinate is defined as 1/(σC2εC)r3duαγr/drgαγOZrgαγMDr. The left axis of the ordinate corresponds to the RDFs.

FIG. 9.

Differences of the integrands for p, obtained for model type 3 at ρC*=ρA*=0.45 and T* = 30.0: (a) for cation–cation pair under the HNC approximation; (b) for cation–anion pair under the HNC approximation; (c) for cation–cation pair under the KH approximation; (d) for cation–anion pair under the KH approximation; (e) for cation–cation pair under the KGK approximation; (f) for cation–anion pair under the KGK approximation. The term “diff. in integrand” at the right axis of the ordinate is defined as 1/(σC2εC)r3duαγr/drgαγOZrgαγMDr. The left axis of the ordinate corresponds to the RDFs.

Close modal

Under the HNC closure, the magnitude of the difference in the integrand for uex became largest at r/σC ≈ 0.8–0.9, which corresponded to the first rising (FR) region of the RDF [Figs. 8(a) and 8(b)]. For the cation–cation pair under either the KH or KGK closure, the largest magnitude of the difference in the integrand for uex was observed at the FR region of the RDF [Figs. 8(c) and 8(e)]. These facts are similar to the results that were reported for LJ fluids.18,19,41 In the case of the cation–anion pair under the KH or KGK closure [Figs. 8(d) and 8(f)], the largest difference in the integrand for uex was observed over the first peak of the RDF, which was wider than the case of the cation–cation pair. It had been reported that for LJ fluids, those differences are significant only in the FR region of the RDF.18,19,41 In contradistinction to the case of LJ fluids, it is striking that the differences in the integrand for uex were significant even at r/σC ≈ 4.0–5.0 in Fig. 8. This is because the interaction potential includes a long-ranged Coulomb term. The difference of the integrand for uex seems to converge beyond r/σC ≈ 7.0. Furthermore, it is to be mentioned that in the region of r/σC > 2.0, the overestimation and the underestimation of the integrand for uex considerably cancel with each other, which is expected from a damped oscillation seen in the curves of the “diff. in integrand.”

The differences in the integrand for p exhibited similar behaviors to those for uex (Fig. 9). The differences in the integrand for p were largest in the FR region and/or over the first peak of the RDF. Because the Coulomb potential is a long-ranged interaction, those differences in the integrand remained significant at r/σC ≈ 4.0–5.0, which almost converged at r/σC ≈ 7.0. In the region of r/σC > 2.0, large canceling between the overestimation and underestimation would take place in the integrand for p because these overestimations and underestimations appear repeatedly. As is similar to the case of uex, it was reported previously18,19,41 that the differences in the integrand for p in LJ fluids are significant only at the FR region of the RDF, which reflects the fact that the LJ potential is relatively short-ranged.

In Sec. III C, it has been revealed that the first peak of the RDF, especially its FR region, is important in terms of the accuracy of uex and p for the current Coulomb systems. Therefore, in this section, the first peak of the RDF is examined. Figure 10 shows the first peaks of the RDFs obtained for model type 3 at ρC*=ρA*=0.45 and T* = 30.0, as examples. The figure shows the RDFs obtained by the HNC, KH, and KGK closures, together with those of the MD simulation. Figures 10(a) and 10(b) show gCCr and gCAr, respectively. Similar data to Fig. 10 obtained at other conditions can be found in Fig. S20 in the supplementary material.

FIG. 10.

The first peaks of the RDFs obtained for model type 3 at ρC*=ρA*=0.45 and T* = 30.0; the panels (a) and (b) show gCCr and gCAr, respectively.

FIG. 10.

The first peaks of the RDFs obtained for model type 3 at ρC*=ρA*=0.45 and T* = 30.0; the panels (a) and (b) show gCCr and gCAr, respectively.

Close modal

Concerning the FR region of the RDF, the HNC closure resulted in a leftward shift in comparison with the MD simulation for both gCCr and gCAr. The leftward shift of the FR region became more prominent when the KH closure was used. In the curve of the RDF under the KGK closure, a sudden rise was observed at a point located at the right of the rising point of the MD curve. The slope of the FR region under the KGK closure was steeper than that of the MD simulation. The curve of gCAr obtained by the KGK closure approached that of the KH closure in the region further than the rising point. These results qualitatively agree with those obtained for LJ fluids.41 For gCCr, the result of the KGK closure was somewhat deviated from that of the KH closure over the first peak, though the former curve roughly followed the latter one. These deviations in the RDF curves obtained by the HNC, KH, or KGK closure from those of the MD simulations cause errors in the thermodynamics such as uex and p, as seen in Sec. III C.

We have attempted to extend the correlation functions obtained by the MD simulation using the OZ theory for a simple molten salt that was modeled as a superposition of LJ and Coulomb potentials. The extension method employed the hybrid closure, where the short-range part of the RDF was given by the data of an MD simulation and its long-range part was described by an existing closure relation coupled with the OZ equation. This study mainly used the MD–KGK hybrid closure, while the MD–HNC and MD–KH hybrid closures were also tested to confirm that these three hybrid closures gave almost equivalent results to one another. Judging from the behavior of the isothermal compressibility, internal energy, and pressure, we have found that the appropriate distance r0 at which the short-range MD part and the long-range OZ part were switched was 2.0 < r0/σC < 3.0. This conclusion is similar to that obtained for LJ fluids in the previous studies.15,16 The bridge function was extracted using the MD–KGK hybrid closure method. The bridge function for a cation–cation pair bCCr and that for a cation–anion pair bCAr were significantly different from each other. When the density was relatively high, both bridge functions showed steep increases with r in the region of r/σC < 1.0, where the value of bCAr was greater than that of bCCr. At a relatively low-density state, bCCr showed an increase with r in the repulsive region (r/σC < 1.0), whereas there was a case where bCAr decreased with r in the region.

Furthermore, we have examined the accuracy of the HNC, KH, and KGK closures in terms of the excess internal energy uex, pressure p, and isothermal compressibility κT. In high-density regions, these closures tended to overestimate uex and p. The HNC closure overestimated κT values when the density became high. The value of p obtained by KGK closure was underestimated in the intermediate-density region. These results are qualitatively similar to those observed for LJ fluids.17–19,41 However, the errors in the integrands for uex and p were found to be considerably different from those of LJ fluids. In LJ fluids, those errors are limited almost only to the FR region (say, r/σC ≈ 0.8–0.9) of the RDF,18,19,41 whereas those of the Coulomb systems were significant even at r/σC ≈ 4.0–5.0. Nevertheless, in terms of the values of uex and p, the overestimation and the underestimation of the integrands would considerably cancel out each other in the region of r/σC > 2.0. The most important region for accurately estimating thermodynamic quantities such as uex and p is the first peak of the RDF, including the FR region. The errors in the FR region of the RDFs obtained by the HNC, KH, and KGK closures were qualitatively similar to the behaviors for LJ fluids. Namely, the FR region tended to shift toward the left when we used these closure relations in comparison with the MD simulation.

See the supplementary material for additional explanations and figures displaying the data other than those shown in the main text.

This work was supported by the JSPS KAKENHI (Grant Nos. 16K05659 and 23K04666). The computation was performed using the Research Center for Computational Science, Okazaki, Japan (Project Nos. 20-IMS-C031, 21-IMS-C027, 22-IMS-C207, and 23-IMS-C148).

The authors have no conflicts to disclose.

Tatsuhiko Miyata: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (equal); Methodology (lead); Resources (lead); Software (equal); Supervision (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). Yu Funahara: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal). Seiya Omori: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal). Taro Shinjo: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
Molecular Theory of Solvation
, edited by
F.
Hirata
(
Kluwer
,
Dordrecht
,
2003
).
2.
T.
Imai
,
A.
Kovalenko
, and
F.
Hirata
,
Chem. Phys. Lett.
395
,
1
(
2004
).
3.
T.
Imai
,
R.
Hiraoka
,
A.
Kovalenko
, and
F.
Hirata
,
J. Am. Chem. Soc.
127
,
15334
(
2005
).
4.
N.
Yoshida
,
S.
Phongphanphanee
,
Y.
Maruyama
,
T.
Imai
, and
F.
Hirata
,
J. Am. Chem. Soc.
128
,
12042
(
2006
).
5.
S.
Phongphanphanee
,
N.
Yoshida
, and
F.
Hirata
,
J. Am. Chem. Soc.
130
,
1540
(
2008
).
6.
Y.
Kiyota
,
R.
Hiraoka
,
N.
Yoshida
,
Y.
Maruyama
,
T.
Imai
, and
F.
Hirata
,
J. Am. Chem. Soc.
131
,
3852
(
2009
).
7.
T.
Imai
,
K.
Oda
,
A.
Kovalenko
,
F.
Hirata
, and
A.
Kidera
,
J. Am. Chem. Soc.
131
,
12430
(
2009
).
8.
T.
Miyata
and
J.
Thapa
,
Chem. Phys. Lett.
604
,
122
(
2014
).
9.
T.
Miyata
and
Y.
Ebato
,
J. Mol. Liq.
245
,
2
(
2017
).
10.
T.
Miyata
,
Bull. Chem. Soc. Jpn.
90
,
1095
(
2017
).
11.
T.
Miyata
and
N.
Yabuki
,
AIP Adv.
9
,
025310
(
2019
).
12.
T.
Miyata
,
Chem. Phys. Lett.
755
,
137777
(
2020
).
13.
T.
Miyata
,
N.
Yabuki
, and
J.
Leung
,
Chem. Lett.
49
,
1372
(
2020
).
14.
T.
Miyata
and
Y.
Hikasa
,
AIP Adv.
12
,
085206
(
2022
).
15.
T.
Miyata
,
Y.
Ogasawara
,
T.
Fujii
,
D.
Yano
, and
Y.
Ebato
,
J. Mol. Liq.
290
,
111167
(
2019
).
16.
T.
Miyata
,
S.
Nishida
, and
Y.
Ogasawara
,
AIP Adv.
11
,
025026
(
2021
).
17.
J.-P.
Hansen
and
I. R.
McDonald
,
Theory of Simple Liquids
, 3rd ed. (
Elsevier; Academic Press
,
London
,
2006
).
18.
Y.
Rosenfeld
and
N. W.
Ashcroft
,
Phys. Rev. A
20
,
1208
(
1979
).
19.
T.
Miyata
and
Y.
Ebato
,
J. Mol. Liq.
217
,
75
(
2016
).
20.
N.
Choudhury
and
S. K.
Ghosh
,
J. Chem. Phys.
116
,
8517
(
2002
).
21.
F. J.
Rogers
and
D. A.
Young
,
Phys. Rev. A
30
,
999
(
1984
).
22.
J. S.
Hoye
and
G.
Stell
,
Mol. Phys.
52
,
1071
(
1984
).
23.
G.
Zerah
and
J.-P.
Hansen
,
J. Chem. Phys.
84
,
2336
(
1986
).
24.
C.
Caccamo
,
G.
Pellicane
, and
E.
Enciso
,
Phys. Rev. E
56
,
6954
(
1997
).
25.
J. A.
Anta
,
E.
Lomba
,
M.
Lombardero
, and
C.
Martin
,
J. Chem. Phys.
105
,
4265
(
1996
).
26.
L.
Verlet
,
Phys. Rev.
165
,
201
(
1968
).
27.
S. M.
Foils
,
N. W.
Ashcroft
, and
L.
Reatto
,
J. Chem. Phys.
81
,
6140
(
1984
).
28.
E.
Lomba
,
M.
Alvarez
,
G.
Stell
, and
J. A.
Anta
,
J. Chem. Phys.
97
,
4349
(
1992
).
29.
S.
Kambayashi
and
Y.
Hiwatari
,
J. Non-Cryst. Solids
156–158
,
80
(
1993
).
30.
S.
Kambayashi
and
J.
Chihara
,
Phys. Rev. E
50
,
1317
(
1994
).
31.
S. J.
Weiner
,
P. A.
Kollman
,
D. A.
Case
,
U. C.
Singh
,
C.
Ghio
,
G.
Alagona
,
S.
Profeta
, Jr.
, and
P.
Weiner
,
J. Am. Chem. Soc.
106
,
765
(
1984
).
32.
B. R.
Brooks
,
R. E.
Bruccoleri
,
B. D.
Olafson
,
D. J.
States
,
S.
Swaminathan
, and
M.
Karplus
,
J. Comput. Chem.
4
,
187
(
1983
).
33.
W. L.
Jorgensen
,
D. S.
Maxwell
, and
J.
Tirado-Rives
,
J. Am. Chem. Soc.
118
,
11225
(
1996
).
34.
X.
Daura
,
A. E.
Mark
, and
W. F.
van Gunsteren
,
J. Comput. Chem.
19
,
535
(
1998
).
35.
A. E.
Kobryn
,
S.
Gusarov
, and
A.
Kovalenko
,
J. Phys.: Condens. Matter
28
,
404003
(
2016
).
36.
A.
Kovalenko
and
F.
Hirata
,
J. Chem. Phys.
110
,
10095
(
1999
).
37.
A.
Ben-Naim
,
Molecular Theory of Solutions
(
Oxford University Press
,
Oxford
,
2006
).
38.
A.
Kovalenko
,
S.
Ten-no
, and
F.
Hirata
,
J. Comput. Chem.
20
,
928
(
1999
).
39.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Clarendon
,
Oxford
,
1987
).
40.
D.
Frenkel
and
B.
Smit
,
Understanding Molecular Simulation: From Algorithms to Applications
(
Academic
,
San Diego
,
2002
).
41.
T.
Miyata
and
K.
Tange
,
Chem. Phys. Lett.
700
,
88
(
2018
).

Supplementary Material