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.

## I. INTRODUCTION

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 18–25). 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.

## II. THEORY AND COMPUTATION

### A. Interaction potential

*σ*,

*ɛ*, 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).

### B. Hybrid closure method

*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*.

^{35}It is written as

*k*

_{B}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\alpha \gamma MD(r)$ represents the RDF obtained from the MD simulations, which is given by

*α*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.

^{1,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

### C. Thermodynamic quantities

*u*

^{ex}was evaluated from

*ρ*=

*ρ*

_{C}+

*ρ*

_{A},

*x*

_{C}=

*ρ*

_{C}/

*ρ*, and

*x*

_{A}=

*ρ*

_{A}/

*ρ*. The pressure

*p*was obtained via

^{37}was employed. It is written as

*η*and

*ζ*for the current system are expressed by

*G*

_{αγ}is defined by

^{1,17}

^{1,36}and KGK

^{35}closures are

*bare*HNC,

*bare*KH, and

*bare*KGK approximations are obtained.

### D. Molecular model and computational details

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 *q*_{C} = −*q*_{A}, which were 0.3*e*, 0.6*e*, 0.8*e*, and 1.0*e*. The densities assumed in this study are $\rho C\sigma C3=\rho A\sigma 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 *q*_{C} = −*q*_{A} = 0.3*e*, 0.6*e*, 0.8*e*, and 1.0*e*, *k*_{B}*T*/*ɛ*_{C} = *k*_{B}*T*/*ɛ*_{A} = 10.0, 15.0, 30.0, and 45.0, respectively. For later convenience, we introduce the notation $\rho C*$, $\rho A*$, and *T** as $\rho C*=\rho C\sigma C3$, $\rho A*=\rho A\sigma A3$, and *T** = *k*_{B}*T*/*ɛ*_{C} = *k*_{B}*T*/*ɛ*_{A}, respectively. The details of the model parameters and thermodynamic conditions are summarized in Table I. Hereafter, the models of *q*_{C} = −*q*_{A} = 0.3*e*, 0.6*e*, 0.8*e*, and 1.0*e* are referred to as the model types 1, 2, 3, and 4, respectively.

Model . | σ_{C}(=σ_{A}) (Å)
. | ɛ_{C}(=ɛ_{A}) (kcal/mol)
. | q_{C} (e)
. | q_{A} (e)
. | $\rho C*(=\rho 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)
. | q_{C} (e)
. | q_{A} (e)
. | $\rho C*(=\rho 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\alpha \gamma 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 *u*^{ex}, 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.

## III. RESULTS AND DISCUSSION

### A. Dependence of thermodynamic quantities on *r*_{0}

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

Examples of the dependence of *κ*_{T} on *r*_{0} are shown in Fig. 1. Figure 1(a) corresponds to the results obtained for model type 1 at $\rho C*=\rho A*=0.45$ and *T** = 10.0, whereas (b) for model type 3 at $\rho C*=\rho 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 *r*_{0}-dependence of *u*^{ex} and that of *p*, respectively, are shown. The panels (a) in Figs. 2 and 3 show the results for model type 1 at $\rho C*=\rho A*=0.45$ and *T** = 10.0, while the panels (b) correspond to those for model type 3 at $\rho C*=\rho A*=0.45$ and *T** = 30.0. The results of *u*^{ex} 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.

The isothermal compressibility *κ*_{T} of the hybrid MD–KGK closure is almost constant when *r*_{0}/*σ*_{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 < *r*_{0}/*σ*_{C} < 2.0, the curves of the hybrid MD–KGK closure in Fig. 1 exhibit some oscillations. When the *r*_{0}/*σ*_{C} value becomes larger than 2.0, the curves in Fig. 1 show a plateau again. In the region of *r*_{0}/*σ*_{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 < *r*_{0}/*σ*_{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 *u*^{ex} (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 *r*_{0}/*σ*_{C} < 0.7. This indicates that the values of *u*^{ex} 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 *r*_{0}/*σ*_{C} are greater than 1.5, both *u*^{ex} and *p* values become almost constant irrespective of *r*_{0}/*σ*_{C}. The *u*^{ex} and *p* values in the region of *r*_{0}/*σ*_{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. 1–3, it is possible to conclude that an appropriate value of *r*_{0} should be determined in the region of 2.0 < *r*_{0}/*σ*_{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 *r*_{0}-dependence of *κ*_{T}, that of *u*^{ex}, and that of *p*, respectively, obtained for model type 1 at $\rho C*=\rho 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 *r*_{0}/*σ*_{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 *r*_{0}/*σ*_{C} > 2.0 coincide with one another. This indicates that, as far as the numerical value of *r*_{0} is given appropriately, the MD–HNC, MD–KH, and MD–KGK hybrid closures are almost equivalent to one another.

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 *r*_{0}/*σ*_{C} = 2.5. Examples of $cCCr$ and $cCAr$ are shown in Fig. 5, which was obtained for model type 1 at $\rho C*=\rho 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.

### B. Extracted bridge function

*r*

_{0}/

*σ*

_{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

*T** = 10.0, whereas (b) those for model type 1 at $\rho C*=\rho A*=0.45$ and

*T** = 10.0. Figures 6(c) and 6(d) correspond to the results for model type 3 at $\rho C*=\rho A*=0.05$ and

*T** = 30.0 and to those at $\rho C*=\rho A*=0.45$ and

*T** = 30.0, respectively.

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 fluids^{20} 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.

### C. Thermodynamics of bare HNC, KH, and KGK closures

In this section, existing closures are examined in terms of the accuracy of thermodynamics. The excess internal energy *u*^{ex}, 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 *u*^{ex}, *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.

The HNC closure provided relatively accurate *u*^{ex} values, though at high density regions, it tended to overestimate *u*^{ex} 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 *u*^{ex} at the density region we examined. The *u*^{ex} 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 *u*^{ex} 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 *u*^{ex} 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 *u*^{ex} and *p*, respectively, obtained for model type 3 at $\rho C*=\rho 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/(\sigma C2\epsilon C)r2u\alpha \gamma (r)g\alpha \gamma OZr\u2212g\alpha \gamma MDr$ and $\u22121/(\sigma C2\epsilon C)r3du\alpha \gamma r/drg\alpha \gamma OZr\u2212g\alpha \gamma MDr$, respectively, where $g\alpha \gamma OZr$ and $g\alpha \gamma 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 *u*^{ex} 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.

Under the HNC closure, the magnitude of the difference in the integrand for *u*^{ex} 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 *u*^{ex} 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 *u*^{ex} 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 *u*^{ex} 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 *u*^{ex} 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 *u*^{ex} 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 *u*^{ex} (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 *u*^{ex}, it was reported previously^{18,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.

### D. First peak of the radial distribution function

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 *u*^{ex} 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 $\rho C*=\rho 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.

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 *u*^{ex} and *p*, as seen in Sec. III C.

## IV. CONCLUSION

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 *r*_{0} at which the short-range MD part and the long-range OZ part were switched was 2.0 < *r*_{0}/*σ*_{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 *u*^{ex}, pressure *p*, and isothermal compressibility *κ*_{T}. In high-density regions, these closures tended to overestimate *u*^{ex} 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 *u*^{ex} 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 *u*^{ex} 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 *u*^{ex} 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.

## SUPPLEMENTARY MATERIAL

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

## ACKNOWLEDGMENTS

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).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

## REFERENCES

*Theory of Simple Liquids*

*Molecular Theory of Solutions*

*Computer Simulation of Liquids*

*Understanding Molecular Simulation: From Algorithms to Applications*