Effect of truncating electrostatic interactions on predicting thermodynamic properties of water–methanol systems

ABSTRACT

The combination of the TraPPE and OPLS/2016 force fields with five water models, TIP3P, SPC/E, OPC, TIP4P/2005 and TIP4P/EW was used to compute mixing enthalpies, excess chemical potentials and activity coefficients of water and methanol. Excess chemical potentials and activity coefficients were computed in an expanded version of the NPT ensemble. We found the best agreement between experimental data for all the computed properties of water–methanol mixtures for the TIP4P/2005-TraPPE force fields. The performance of the spherical cutoff methods in MC and MD simulations was compared to the Ewald summation. The radial distribution functions obtained from the Ewald summation and the Damped-Shifted Force (DSF) method were in excellent agreement. Numerical artifacts appeared at the cutoff radius when the original Wolf method was used to calculate the electrostatic interactions. The calculated excess mixing enthalpies, excess chemical potentials, and activity coefficients of water and methanol obtained from the Wolf method were in good agreement with the DSF method. Our simulation results show that the numerical artifacts of the original Wolf method have little effect for energy calculations in aqueous methanol mixtures.

KEYWORDS:

• Water
• water–methanol mixtures
• wolf method
• electrostatics
• chemical potential
• activity coefficients
• excess enthalpy of mixing

1. Introduction

Water is one of the most important molecules in science [1,2], central to life, and probably the most studied molecule in the field of molecular simulation [3–7]. From a chemistry perspective, it is a simple molecule formed by two hydrogen atoms and one oxygen atom. The liquid phase has unique properties and complex behaviour [5,8–10]. To date, our understanding of intermolecular and intramolecular interactions of water is far from complete [4,5]. To reproduce different thermophysical properties of water, numerous theoretical and computational models for water were developed and published [3,6,11–13]. In the last decades, due to exponential increase in computational power, molecular simulations are used to compute many properties of water [4–6] in biomolecular [14,15], chemical and industrial applications [12,13,16]. The performance of different water models used to describe the intermolecular interactions are reported extensively in literature [4,17], and is central for predicting properties of water and reproducing experimental data. Since water models are usually fitted to a limited set of experimental data, no water model can simultaneously reproduce all thermophysical properties in good agreement with experiments [4–6]. Therefore, reproducing thermophysical properties of water depends strongly on the choice of experimental data used to fit the intermolecular interaction parameters [17]. Deviations from experimental data may also arise because of inherent limitations in the molecule models of water [3,5,6,18].

Water is a flexible and polarisable molecule [11,19,20]. This means that the electronic structure of the water molecule undergoes deformation due to the electric field induced by the surrounding water or other polar molecules [21]. To account for polarisation effects in the recent years, polarisable force fields for water have been developed [11,20–28]. Some properties of water, such as vapour pressure, critical properties, dielectric constant, and virial coefficient are most accurately predicated by considering polarisation effects [20]. However, the performance of different polarisable force fields to compute chemical potentials and activity coefficients of water is not fully investigated in literature.

Although the water molecule is flexible and polarisable [11], including polarisation effects in MC simulations significantly increases the computational costs [29–31]. In many studies, water is considered as an explicit solvent/medium and is not the main focus. Therefore, the majority of water models in literature are rigid with point charges, and polarisation is ignored [4–6,11]. These simplified models for water are computationally advantageous and the reproduced bulk properties of water are usually in good agreement with experiments at ambient conditions [4,11]. The most popular classical force fields of water include three-point site interaction models [32–34], four-point interaction site models [5], five-point interaction site models [35,36], and six-point interaction site models [8,37]. Many studies compared the performance of different interaction potentials for predicting thermophysical properties of pure water [5,6].

Aqueous mixtures of methanol are investigated frequently in academia [38–57] and are of practical importance in industrial applications [58–66]. Depending on the application, different mixture properties are required for process design. To compute the mixture properties at different conditions, different force field combinations of water and methanol are considered in few molecular simulation studies [7,38,39,41–45,53,57]. To the best of our knowledge, chemical potentials and activity coefficients of water and methanol are not reported in molecular simulation studies for different/recent force field combinations of water and methanol. Because of the low vapour pressures of water and methanol at ambient conditions, it is also experimentally challenging to determine the activity coefficients.

The most commonly used method to compute the chemical potential in molecular simulation is the well-known Widom's Test Particle Insertion (WTPI) method [67–69]. In this method, the chemical potential of a component is determined by sampling the interaction energy of a test molecule, inserted at a randomly selected position. Successful sampling of the chemical potential depends strongly on frequent occurrence of spontaneous cavities in the system large enough to accommodate the test molecules [70–72]. Due to the high density of aqueous methanol mixtures at ambient conditions, these spontaneous cavities occur rarely, rendering the WTPI method practically useless [70]. Therefore, advanced simulation techniques are required to compute chemical potentials and activity coefficients in aqueous methanol mixtures [57,67,68,70,71,73–79]. Here, we investigate the performance of different combinations of water and methanol force fields for the computation of chemical potentials and activity coefficients, using state of the art molecular simulation techniques, based on the idea of expanded ensembles [76–78,80]

We focus on popular rigid non-polarisable water models of water. Popular force fields used for water include: TIP3P [33], SPC/E [34], OPC [81], TIP4P/2005 [17], and TIP4P/EW [82]. These models are frequently used in benchmark studies [3–6]. In this study, we have chosen for the TraPPE [83] and OPLS/2016 [84] force fields as the TraPPE force field is widely used, and the OPLS/2016 force field is recently published. The OPLS/2016 force field is claimed to be a more accurate force field for methanol [84]. The force fields for methanol are rigid and non-polarisable. These models are different in bond geometries and/or point charge distributions, and predicting different properties of the mixture depends on different combinations of these force fields [39,40]. The potential energy of the system of rigid molecules with point charges is the sum of pairwise interaction potentials consisting of Lennard–Jones (LJ) and Coulombic interactions:(1) $E=\frac{1}{2}\sum _i\sum _{j\ne i}\left[4{\epsilon }_{ij}\left[{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^{12}-{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^6\right]+\frac{1}{4\pi {\epsilon }_0}\frac{q_i{q}_j}{r_{ij}}\right],$(1) where $r_{ij}$ where the double summation is over all the interaction sites. $r_{ij}$ is the distance between atoms i and j, ${\epsilon }_0$ is the dielectric constant, $q_i$ is the partial charge of atom i, ${\sigma }_{ij}$ and ${\epsilon }_{ij}$ are the LJ parameters between atoms i and j, obtained from ${\sigma }_i$, ${\sigma }_j$, ${\epsilon }_i$ and ${\epsilon }_j$. Proper treatment of LJ and Coulombic interactions are essential for an accurate description of molecular structure, and thermodynamic properties of water and methanol [69,85–93].

In molecular simulation, periodic boundary conditions are applied to simulate bulk phases and compute structural and thermodynamic properties [69]. In systems with periodic boundary conditions, the long-range electrostatic interactions decay slowly with $r^{-1}$, and the treatment of the electrostatic interactions becomes computationally demanding. Direct summation of Coulombic interactions is conditionally convergent which means that the results depend on the order of summation [86,94]. Various molecular simulation studies investigated the accuracy and scalability of different methods for treating long-range Coulombic interactions [87,95,96]. These methods can be divided into Ewald-based and cutoff-based methods. In this paper, we investigate to what extent different electrostatic methods influence thermophysical properties of water–methanol mixtures, and whether numerical artifacts are observed. A brief description on electrostatic methods is provided below.

The Ewald summation [69,86,87,92,97] is the most widely used and accepted method to compute electrostatic interactions in molecular simulation and scales as $\mathcal{O}\left(N^2\right)$. In the Ewald summation, the electrostatic interactions are split into effective short-range interactions, evaluated in real space, and a Fourier series to account for the long-range contribution of the electrostatic interactions, evaluated in reciprocal space [69,86,88,92]. The high computational costs of the Ewald summation have lead to efforts to develop faster alternatives [88,90,93,94,96,98]. In the past decades, more Ewald-based algorithms have been developed to reduce the time complexity to $\mathcal{O}(N\log N)$ by optimising the reciprocal space summation [99,100]. Particle mesh algorithms such as Staggered Mesh Ewald method (StEM) [99], Particle-Particle Particle-Mesh (PPPM) [101], and Particle-Mesh Ewald (PME) [100] are based on the Ewald method and scale as $\mathcal{O}(N\log N)$ [87,102]. These algorithms allow efficient parallel implementations, and are especially advantageous to use in MD simulations [87]. When the Ewald summation is used in Monte Carlo (MC) simulations, one only need consider the charged atoms that are changed in a trial move. This can be achieved by storing the Fourier part in memory in an efficient way [94]. For a faster computation of the electrostatic interactions, truncation or spherical cutoff-based methods are proposed as an alternative to the Ewald-type methods [89,90,96,103–105].

Cutoff-based methods are based on the idea that the effective electrostatic potential of condensed phases has rather short-ranged behaviour [89,90,93,94,96,98,106,107], and the effect of the long range interactions beyond the cutoff radius becomes negligible due to the screening of charges [87,89,98,104,108]. In the gas phase, the screening of charges is weak compared to the liquid phase [89,90,93,98]. Compared to the Ewald-type methods, the cutoff-based methods are much simpler to implement and scale as $\left(\mathcal{O}\right(N\left)\right)$. The performance of the cutoff-based methods are promising especially for bulk homogeneous systems, while in systems with an interface, severe numerical problems may be expected [88]. A common criticism about cutoff-based methods is that numerical artifacts may arise due to the truncation of long-range interactions [87,88,95,109–111]. These numerical artifacts can affect the structure of the liquid. For example, two atoms of opposite charges (positive-negative) prefer to be located within the cutoff radius, while two atoms with the same charge (positive-positive or negative-negative) are preferentially located outside the cutoff radius. As the concept of using a cutoff radius for electrostatic interactions is due to simulation efficiency reasons rather than a physical effect, changes in the liquid structure due to an imposed cutoff radius are non-physical and should be avoided. To observe whether anomalies occur in the liquid structure, the Radial Distribution Function (RDF) of the liquid can be studied. Depending on the extent of the numerical artifacts, energy calculations are affected. Mark et al. reported artificial structuring of water when using atom-based cutoff methods developed by Steinbach et al. [104] to compute the electrostatic interactions [108]. The slow convergence issue appeared especially in truncation and shifting schemes in which the truncation sphere inside the cutoff was not electroneutral [87,90,91]. Wolf et al. realised that the error in computing the electrostatic interactions was related to the net charge inside the cutoff sphere and enforced charge neutrality within the cutoff radius [88,90]. By using a charge-neutralised damped pair potential, fast convergence is achieved and effects of possible artifacts due to the truncation of electrostatic interactions are minimised [88–91,98].

Using the so-called Wolf method, the damped shifted pairwise electrostatic potential for a system of N charges is obtained from [90,93,94,98](2) $\begin{array}{cc}{E}^{\mathrm{Wolf}}& =\frac{1}{2}\underset{r_{ij}<R_{\mathrm{c}}}{\sum _{i=1}^N\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^N}{q}_i{q}_j\left[\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^N{q}_i^2,\end{array}$(2) where $r_{ij}$ is the distance between the partial charges i and j, $\mathrm{erfc}\left(x\right)=1-\mathrm{erf}\left(x\right)$ is the complementary error function. In principle, other damping functions may be used [85,90]. $R_{\mathrm{c}}$ is the cutoff radius, α is the Wolf damping parameter, and its value determines how fast the complementary error function approaches zero with increasing $r_{ij}$. The second term on the right-hand side of Equation (2(2) $\begin{array}{cc}{E}^{\mathrm{Wolf}}& =\frac{1}{2}\underset{r_{ij}<R_{\mathrm{c}}}{\sum _{i=1}^N\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^N}{q}_i{q}_j\left[\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^N{q}_i^2,\end{array}$(2) ) is the energy associated with the so-called self term [90,98]. The Wolf method has been used in different studies including water and ionic systems [96,112–114]. For a nice and readable derivation of the Wolf method for molecular systems, the reader is referred to the paper by Waibel and Gross [98]. Applying the Wolf method to a molecular system, the charge-charge interactions within each molecule should be excluded from the total sum of the electrostatic energy [98] because these are not part of intermolecular interactions. Initially, we consider a system without intramolecular Coulombic interactions, e.g. for rigid molecules. The total electrostatic energy equals [93,94,115]:(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) In Equation (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) ), $N_{\mathrm{m}}$ is the number of molecules, $N_{\mathrm{a}}^i$ is the number of atoms in molecule i, indices i and j are used to count the number of molecules, indices a and b are used to count atoms within molecules. $q_{ia}$ is the partial charge of atom a in molecule i, $r_{iajb}$ is the distance between interaction sites a and b. The first and second terms on the right-hand side of Equation (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) ) are used to compute the pairwise electrostatic interactions between all partial charges in the system, which are screened electrostatic interactions between different molecules and electrostatic interactions inside molecules, respectively. For rigid molecules, we are only interested in intermolecular interactions, and therefore we have to exclude electrostatic interactions within molecules, i.e. charge-charge interactions between the atoms within the molecules. This is achieved by the third term on the right-hand side of Equation (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) ), the so-called exclusion term, similar to the exclusion term in the Ewald summation [116]. The fourth term on the right-hand side of Equation (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) ) is the self-interaction term which is similar to that of an atomic system equation (2(2) $\begin{array}{cc}{E}^{\mathrm{Wolf}}& =\frac{1}{2}\underset{r_{ij}<R_{\mathrm{c}}}{\sum _{i=1}^N\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^N}{q}_i{q}_j\left[\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^N{q}_i^2,\end{array}$(2) ). For non-rigid molecules, the intramolecular Coulombic interactions should be added to Equation (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) ) by:(4) $E_{\mathrm{Intra}}=\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}{c}_{ab}\frac{q_{ia}{q}_{ib}}{r_{iaib}}$(4) in which $c_{ab}$ is a scaling parameter for electrostatic interactions between atom a and atom b. $c_{ab}$ equals zero if atoms a and b do not have any intramolecular electrostatic interactions. Several force fields have intramolecular interactions that are scaled [117–121], so the value of $c_{ab}$ can be non-integer. Note that a cutoff radius is not applied for the intramolecular interactions of Equation (4(4) $E_{\mathrm{Intra}}=\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}{c}_{ab}\frac{q_{ia}{q}_{ib}}{r_{iaib}}$(4) ), and that these interactions are not screened with a $\mathrm{erfc}\left(\alpha r\right)/r$ term. The reason is that for a system consisting of an isolated molecule with intramolecular interactions but no intermolecular interactions, the correct result is obtained, i.e. a direct pairwise $1/r$ summation over the intramolecular Coulombic interactions. It is important to note that in sharp contrast to the terms for intramolecular Coulombic interactions and intramolecular Coulombic exclusions proposed by Gross et al. [98], in this work no cutoff is imposed for these interactions. It is expected that this avoids potential artifacts for systems including molecules larger than the cutoff radius.

Fennell and Gezelter [89] found that the damped shifted potential proposed by Wolf et al. (Equation (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) )), results in force discontinuity at the cutoff radius in MD simulations. This is undesirable as it may lead to energy drifts [69]. To remove this discontinuity, the electrostatic potential proposed by Fennell and Gezelter (the so-called Damped-Shifted Force (DSF) method) is calculated using(5) $\begin{array}{cc}{E}^{\mathrm{DSF}}& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{R}_{\mathrm{c}}^2\right)}{R_{\mathrm{c}}}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(5) Similar to Equation (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) ), the intramolecular Coulombic interactions are not included here, and these interactions need to be taken into account according to Equation (4(4) $E_{\mathrm{Intra}}=\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}{c}_{ab}\frac{q_{ia}{q}_{ib}}{r_{iaib}}$(4) ). For simplicity, Equations (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) )–(5(5) $\begin{array}{cc}{E}^{\mathrm{DSF}}& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{R}_{\mathrm{c}}^2\right)}{R_{\mathrm{c}}}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(5) ) are formulated for single-component systems. Extending these equations to multicomponent systems is trivial: (1) The summation in the first terms of Equations (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) ) and (5(5) $\begin{array}{cc}{E}^{\mathrm{DSF}}& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{R}_{\mathrm{c}}^2\right)}{R_{\mathrm{c}}}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(5) ) should be over all intermolecular interactions between all the atoms of all molecules; (2) the summations in other terms of Equations (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) ) and (5(5) $\begin{array}{cc}{E}^{\mathrm{DSF}}& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{R}_{\mathrm{c}}^2\right)}{R_{\mathrm{c}}}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(5) ) (and the summation in Equation (4(4) $E_{\mathrm{Intra}}=\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}{c}_{ab}\frac{q_{ia}{q}_{ib}}{r_{iaib}}$(4) )) should be over all molecules. The extra term on the right-hand side of Equation (5(5) $\begin{array}{cc}{E}^{\mathrm{DSF}}& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{R}_{\mathrm{c}}^2\right)}{R_{\mathrm{c}}}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(5) ) compared to Equation (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) ) makes both the potential and its first derivative (the force) continuous at the cutoff radius [89]. Although the continuity of the force at the cutoff radius is of primary interest in MD simulations, the electrostatic potential derived in Equation (5(5) $\begin{array}{cc}{E}^{\mathrm{DSF}}& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{R}_{\mathrm{c}}^2\right)}{R_{\mathrm{c}}}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(5) ) is slightly different compared to the original Wolf method. This small difference in the electrostatic interaction potential can potentially change the computed structure and/or other properties of the system. In different studies of ionic liquids [89,100,122–124], it was found that the electrostatic energies and forces obtained using the DSF method are in excellent agreement with the conventional Ewald/smooth PME method.

Instead of using an atom-based spherical cutoff method, one could also apply a group-based or charge-group cutoff method [87,89,104]. In this method, atoms are assigned to charge-neutral groups. If the distance between the geometric centres of the charge-neutral groups is smaller than a certain cutoff, all atoms belonging to the charge-neutral groups interact, otherwise no atomic interaction between the two groups is taken into account. Therefore the computational time is reduced by ignoring the pairwise electrostatic interactions between charge-neutral groups which are further away from the cutoff radius. Another advantage is that the leading $1/r$ term for Coulombic interactions reduces to a higher order term that decays faster [104]. Gross et al. [98] found no differences between a group-based cutoff radius and an atom-based cutoff radius for the Wolf method, and therefore we have not considered this further in this paper.

This paper is organised as follows. In Section 2, expressions to compute the following thermodynamic properties: excess mixing enthalpy, chemical potentials and activity coefficients are described. Simulation details, force field parameters, and scaling of the intramolecular interactions are described in Section 3. Our simulation results are presented in Section 4. It is shown that the liquid structure obtained using the Wolf method, Equation (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) ), and the DSF method, Equation (5(5) $\begin{array}{cc}{E}^{\mathrm{DSF}}& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{R}_{\mathrm{c}}^2\right)}{R_{\mathrm{c}}}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(5) ), is different especially near the cutoff radius. However, computing the excess chemical potentials of water and methanol obtained using the Wolf method and the DSF method yield very similar results. It is also shown that the activity coefficients of water and methanol in water–methanol mixtures computed based on TIP4P/2005 and TraPPE force fields show the best agreement with experimental data. Our conclusions are summarised in Section 5.

2. Theory

The enthalpy of a system (pure component or a mixture) can be directly computed from simulations in the NPT ensemble [69,86](6) $H={⟨U⟩}_{NPT}+P{⟨V⟩}_{NPT},$(6) ${⟨U⟩}_{NPT}$ is the ensemble average internal energy of the system, and ${⟨V⟩}_{NPT}$ is the ensemble average volume in the NPT ensemble. For the water–methanol mixture, the excess enthalpy of the mixing is obtained from(7) $h_{\mathrm{mix}}^{\mathrm{ex}}=h_{\mathrm{mix}}-(1-x_{\mathrm{MeOH}})h_{{\mathrm{H}}_2\mathrm{O}}-x_{\mathrm{MeOH}}{h}_{\mathrm{MeOH}},$(7) $h^{\mathrm{ex}}$ is the excess enthalpy of mixing with respect to pure liquid, $h_{\mathrm{mix}}$ is the enthalpy of the water–methanol mixture, $h_{\mathrm{MeOH}}$ is the enthalpy of pure methanol, $h_{{\mathrm{H}}_2\mathrm{O}}$ is the enthalpy of pure water, and $x_{\mathrm{MeOH}}$ is the mole fraction of methanol. The most common method to compute the chemical potential of a component is the WTPI method [67,68]. Because of inefficient sampling of the chemical potentials of water and methanol using the WTPI method at low temperatures (298 K) [70,72], the Continuous Fractional Component NPT (CFCNPT) ensemble is used [73] to compute the chemical potentials. For the general case of a multicomponent mixture in this ensemble, the NPT ensemble is expanded with a fractional molecule of component A, as opposed to all other molecules which are referred to as ‘whole molecules’. For the complete expression of the partition function of this expanded NPT ensemble, the reader is referred to Ref. [73]. The interaction of the fractional molecule is scaled with a coupling parameter ${\lambda }_{\mathrm{A}}\in [0,1]$. For ${\lambda }_{\mathrm{A}}=0$, the fractional molecule does not interact with other molecules in the systems, and for ${\lambda }_{\mathrm{A}}=1$ the fractional molecule is fully interacting with the rest of the molecules in the system, and thus acts as a ‘whole’ molecule. In our previous work, we have shown that the chemical potential of component A in the CFCNPT ensemble is obtained according to [72,73](8) ${\mu }_{\mathrm{A}}=-\frac{1}{\beta }\ln \left(\frac{{⟨V⟩}_{\mathrm{CFC}\mathit{N}\mathit{P}\mathit{T}}}{{\mathrm{\Lambda }}_{\mathrm{A}}^3\left(N_{\mathrm{A}}+1\right)}\right)-\frac{1}{\beta }\ln \left(\frac{p({\lambda }_{\mathrm{A}}\uparrow 1)}{p({\lambda }_{\mathrm{A}}\downarrow 0)}\right),$(8) ${\mu }_{\mathrm{A}}$ is the chemical potential of component A, $\beta =1/\left(k_{\mathrm{B}}T\right)$, and $k_{\mathrm{B}}$ is the Boltzmann constant. ${⟨V⟩}_{\mathrm{CFC}\mathit{N}\mathit{P}\mathit{T}}$ is the ensemble average volume in the CFCNPT ensemble, ${\mathrm{\Lambda }}_{\mathrm{A}}$ is the thermal wavelength, $N_{\mathrm{A}}$ is the number of whole molecules of type A, $p({\lambda }_{\mathrm{A}}\downarrow 0)$ is the probability of ${\lambda }_{\mathrm{A}}$ approaching zero, and $p({\lambda }_{\mathrm{A}}\uparrow 1)$ is the probability of ${\lambda }_{\mathrm{A}}$ approaching one. The second term in Equation (8(8) ${\mu }_{\mathrm{A}}=-\frac{1}{\beta }\ln \left(\frac{{⟨V⟩}_{\mathrm{CFC}\mathit{N}\mathit{P}\mathit{T}}}{{\mathrm{\Lambda }}_{\mathrm{A}}^3\left(N_{\mathrm{A}}+1\right)}\right)-\frac{1}{\beta }\ln \left(\frac{p({\lambda }_{\mathrm{A}}\uparrow 1)}{p({\lambda }_{\mathrm{A}}\downarrow 0)}\right),$(8) ) is the excess chemical potential of component A (with respect to the ideal gas phase). In [73], we show that the WTPI method fails to compute the excess chemical potential of pure water and methanol at ambient conditions. The activity coefficient of component A, ${\gamma }_{\mathrm{A}}$, in a mixture is obtained according to [125,126](9) ${\gamma }_{\mathrm{A}}=\frac{{\rho }_{\mathrm{A}}}{x_{\mathrm{A}}{\rho }_{0\mathrm{A}}}\cdot \exp \left[\beta \left({\mu }_{\mathrm{A}}^{\mathrm{ex}}-{\mu }_{0\mathrm{A}}^{\mathrm{ex}}\right)\right],$(9) ${\rho }_{\mathrm{A}}$ and ${\rho }_{0\mathrm{A}}$ are the number densities of component A in the mixture and the reference number density of the pure solvent, respectively. ${\mu }_{\mathrm{A}}^{\mathrm{ex}}$ is the excess chemical potential of component A in a mixture with mole fraction of $x_{\mathrm{A}}$, and ${\mu }_{0\mathrm{A}}^{\mathrm{ex}}$ is the excess chemical potential of A in the pure fluid A.

3. Simulation details

3.1. Monte Carlo simulations

Different water–methanol mixtures with compositions ranging from $x_{\mathrm{MeOH}}=0$ to $x_{\mathrm{MeOH}}=1$ are simulated at P=1 bar and T=298 K, both in the conventional NPT ensemble [69,86] and the CFCNPT ensemble [73]. All MC simulations were performed using our in-house code which is verified to produce the same results as the RASPA software package [127,128] in various works [93,115]. For water, the TIP3P [32], SPC/E [34], OPC [81], TIP4P/2005 [17], and TIP4P/EW [82] force fields, and for methanol the TraPPE [83] and OPLS/2016 [84] force fields are used. All force field parameters are provided in Tables 1 and 2. All molecules are rigid and the interactions between the molecules only consist of LJ and Coulombic interactions. Periodic boundary conditions were used. LJ potentials were truncated but not shifted, and analytic tail corrections and the Lorentz–Berthelot mixing rules were applied [69,86]. For every mixture composition and every water–methanol force field combination, the Wolf and the DSF methods, Equations (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) ) and (5(5) $\begin{array}{cc}{E}^{\mathrm{DSF}}& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{R}_{\mathrm{c}}^2\right)}{R_{\mathrm{c}}}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(5) ), were both used to compute the electrostatic interactions.

Table 1. Force field parameters for water used in this study.

Force field	TIP3P [32]	SPCE/E [34]	OPC [81]	TIP4P/2005 [17]	TIP4P/EW [82]
${\epsilon }_{\mathrm{O}\mathrm{O}}$	76.500	78.175	107.086	93.196	81.899
${\sigma }_{\mathrm{O}\mathrm{O}}$	3.1506	3.1660	3.1666	3.1589	3.1644
$q_{\mathrm{O}}$	−0.8340	−0.8476	−1.3582	–	–
$q_{\mathrm{H}}$	0.41700	0.42380	0.67910	0.55640	0.52422
$q_{\mathrm{L}}$	–	–	–	−1.11280	−1.04844
$r_{\mathrm{OH}}$	0.9572	1.0000	0.8724	0.9572	0.9572
$r_{\mathrm{OL}}$	–	–	0.1594	0.1546	0.1250

Notes: All molecules are considered rigid. ε is reported in units of K, σ and r are reported in units of Å.

Table 2. Force field parameters for methanol used in this study.

Force field	OPLS/2016 [84]	TraPPE [83]
${\epsilon }_{\mathrm{O}\mathrm{O}}$	97.775	93.000
${\sigma }_{\mathrm{O}\mathrm{O}}$	3.1659	3.0200
${\epsilon }_{{\mathrm{CH}}_3{\mathrm{CH}}_3}$	110.450	98.000
${\sigma }_{{\mathrm{CH}}_3{\mathrm{CH}}_3}$	3.6449	3.75
$q_{\mathrm{O}}$	−0.6544	−0.70000
$q_{\mathrm{H}}$	0.49980	0.43500
$q_{{\mathrm{CH}}_3}$	0.1546	0.2650
$r_{\mathrm{OH}}$	0.9450	0.9450
$r_{{\mathrm{CH}}_3\mathrm{O}}$	1.43	1.43

Notes: All molecules are considered rigid. ε is reported in units of K, σ and r are reported in units of Å.

To obtain the parameters for the Wolf method for pure water and pure methanol, independent simulations were performed in the NVT ensemble, using SPC/E [34] and OPLS/2016 [84] force fields, close to the experimental densities at T=298 K [129,130]. For dense liquids such as water and methanol at ambient conditions, it is sufficient to plot Figure 1 for a single configuration [93]. For the single equilibrated configuration, the electrostatic energies were calculated for different values of cutoff radii ranging from $R_{\mathrm{c}}=10$ Å to $R_{\mathrm{c}}=15$ Å, as a function of α. The relative difference in electrostatic energies, for water and methanol, were compared to the results obtained from the Ewald summation, and the results are shown in Figure 1. It is shown in Figure 1 that the relative difference between the electrostatic energies between these methods is within 0.5% for the cutoff radii ranging from 10 Å to 15 Å, and α ranging from 0.1 ${\text{Å}}^{-1}$ to 0.15 ${\text{Å}}^{-1}$. This means that the results obtained from the Wolf method in this (α,$R_{\mathrm{c}}$) range, are consistent with the energetics from the Ewald summation. For all simulations, Wolf parameters, $R_{\mathrm{c}}$ and α were set to 14 Å and 0.12 ${\text{Å}}^{-1}$, respectively, and a cutoff radius of 14 Å was used for LJ interactions. The same parameters were used for the DSF method (Equation (5(5) $\begin{array}{cc}{E}^{\mathrm{DSF}}& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{R}_{\mathrm{c}}^2\right)}{R_{\mathrm{c}}}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(5) )).

Figure 1. (Colour online) Relative differences in computed electrostatic energies between the Wolf method, Equation (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) ), and the Ewald summation for (a) water and (b) methanol. The parameters for the Ewald summation are calculated based on relative precision of ${10}^{-6}$ [93]. The SPC/E [34] and OPLS/2016 [84] force fields were used to obtain the densities of water and methanol at T=298 K and P=1 bar. Individual configurations were obtained at constant densities of 1000 $\mathrm{kg}{\mathrm{m}}^{-3}$ and 748 $\mathrm{kg}{\mathrm{m}}^{-3}$ for water and methanol, respectively.

Simulations in the NPT ensemble were performed to compute the excess mixing enthalpies of water–methanol mixtures (with respect to pure liquid) based on Equation (7(7) $h_{\mathrm{mix}}^{\mathrm{ex}}=h_{\mathrm{mix}}-(1-x_{\mathrm{MeOH}})h_{{\mathrm{H}}_2\mathrm{O}}-x_{\mathrm{MeOH}}{h}_{\mathrm{MeOH}},$(7) ). Simulations in the CFCNPT ensemble were performed to compute the excess chemical potentials of water and methanol (with respect to ideal gas phase) using Equation (8(8) ${\mu }_{\mathrm{A}}=-\frac{1}{\beta }\ln \left(\frac{{⟨V⟩}_{\mathrm{CFC}\mathit{N}\mathit{P}\mathit{T}}}{{\mathrm{\Lambda }}_{\mathrm{A}}^3\left(N_{\mathrm{A}}+1\right)}\right)-\frac{1}{\beta }\ln \left(\frac{p({\lambda }_{\mathrm{A}}\uparrow 1)}{p({\lambda }_{\mathrm{A}}\downarrow 0)}\right),$(8) ). The activity coefficients of water and methanol were obtained using Equation (9(9) ${\gamma }_{\mathrm{A}}=\frac{{\rho }_{\mathrm{A}}}{x_{\mathrm{A}}{\rho }_{0\mathrm{A}}}\cdot \exp \left[\beta \left({\mu }_{\mathrm{A}}^{\mathrm{ex}}-{\mu }_{0\mathrm{A}}^{\mathrm{ex}}\right)\right],$(9) ). Each simulation in the NPT ensemble was carried out with ${10}^5$ equilibration cycles and $4\times {10}^6$ production cycles. In each cycle, the number of MC steps equals the total number of molecules, and trial moves were selected with the following probabilities: 1% volume changes, 49.5% translations, and 49.5% rotations. Simulations in the CFCNPT ensemble were performed to compute the chemical potentials and activity coefficients of water and methanol at different mixture compositions. To scale the interactions of the fractional molecule in the CFCNPT ensemble, the LJ and Coulombic interactions are decoupled. It is specifically important to switch on the repulsive LJ interactions before the Coulombic interactions, to protect the charges from overlapping [131–136]. Charge-overlaps can potentially lead to huge electrostatic potentials, inaccuracies and numerical instabilities [131–135]. First, the LJ interactions of the fractional molecule are switched on at $\lambda =0$ and are fully interacting with the surrounding molecules when λ reaches a certain predefined threshold value of ${\lambda }^{\ast }$. Second, the Coulombic interactions of the fractional molecule are switched on at $\lambda ={\lambda }^{\ast }$ and are fully interacting with the surrounding molecules when $\lambda =1$. In Figure 2, it is shown how the LJ and Coulombic interactions are decoupled and scaled with separate coupling parameters ${\lambda }_{\mathrm{LJ}}\in [0,1]$ and ${\lambda }_{\mathrm{Coul}}\in [0,1]$ respectively. The precise functional forms for the scaling of the LJ and electrostatic interactions are discussed in Section 3.1.1. Different choices are possible for ${\lambda }^{\ast }$ depending on the system. Here, we selected ${\lambda }^{\ast }=0.8$ and did not attempt to choose the value of ${\lambda }^{\ast }$ which leads to the most efficient simulation. Note that only the efficiency of the simulation depends on ${\lambda }^{\ast }$ and not the result.

Figure 2. (Colour online) Coupling parameter λ to scale the interactions of fractional molecules. ${\lambda }_{\mathrm{LJ}}\in [0,1]$ is the coupling parameter used to scale the LJ interactions of the fractional molecule (Equation (10(10) $\begin{array}{rl}{u}_{\mathrm{LJ}}& \left(r,{\lambda }_{\mathrm{LJ}}\right)={\lambda }_{\mathrm{LJ}}4\epsilon \\ & \left(\frac{1}{{\left[\frac{1}{2}{\left(1-{\lambda }_{\mathrm{LJ}}\right)}^2+{\left(\frac{r}{\sigma }\right)}^6\right]}^2}-\frac{1}{\left[\frac{1}{2}{\left(1-{\lambda }_{\mathrm{LJ}}\right)}^2+{\left(\frac{r}{\sigma }\right)}^6\right]}\right),\end{array}$(10) )). At $\lambda ={\lambda }^{\ast }$, the Coulombic interactions are switched on. ${\lambda }_{\mathrm{Coul}}\in [0,1]$ is the coupling parameter used to scale the Coulombic interactions of the fractional molecule (Equation (11(11) $\begin{array}{cc}{E}_{\mathrm{Coul}}^{\mathrm{DSF}}\left(r,{\lambda }_{\mathrm{Coul}}\right)& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{\lambda }_{\mathrm{Coul}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha \left[r_{iajb}+r_{\ast }\right]\right)}{r_{iajb}+r_{\ast }}-\frac{\mathrm{erfc}\left(\alpha \left[R_{\mathrm{c}}+r_{\ast }\right]\right)}{R_{\mathrm{c}}+r_{\ast }}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha \left[R_{\mathrm{c}}+r_{\ast }\right]\right)}{{\left(R_{\mathrm{c}}+r_{\ast }\right)}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{\left[R_{\mathrm{c}}+r_{\ast }\right]}^2\right)}{R_{\mathrm{c}}+r_{\ast }}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{\lambda }_{\mathrm{Coul}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha \left[r_{iaib}+r_{\ast }\right]\right)}{r_{iaib}+r_{\ast }}-\frac{\mathrm{erfc}\left(\alpha \left[R_{\mathrm{c}}+r_{\ast }\right]\right)}{R_{\mathrm{c}}+r_{\ast }}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{{\lambda }_{\mathrm{Coul}}{q}_{ia}{q}_{ib}}{r_{iaib}+r_{\ast }}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{\lambda }_{\mathrm{Coul}}{q}_{ia}^2,\end{array}$(11) )).

Beside thermalisation trial moves (translations, rotations, volume changes etc.) [69], three types of trial moves involving the fractional molecule are used to facilitate gradual particle insertions/removals [73]: (1) Changes in λ: changing the coupling parameter (λ) of the fractional molecule while keeping the positions of all molecules including the fractional molecule fixed. The trial move is automatically rejected if $\lambda <0$ or $\lambda >1$ [72,73,115]; (2) Reinsertions: the fractional molecule is moved to a randomly selected position while keeping the value of λ and the positions of all the whole molecules fixed; (3) Identity changes: the fractional molecule is changed into a whole molecule of the same molecule type, and a randomly selected whole molecule of the same type is changed into a fractional molecule. The positions of all molecules and the value of λ are kept fixed in this trial move. These three trial moves are accepted or rejected based on Metropolis acceptance rules [69,115]. Trial moves of type (2) have a high acceptance probability at low values of λ, and trial moves of type (3) has a high acceptance probability only at high values of λ [73,115]. It is possible to define a hybrid trial move in which trial moves of type (2) are only selected at low values of λ, and trial moves of type (3) are only selected at high values of λ. Therefore, trial moves (2) and (3) are only selected when the acceptance probabilities are high. In Refs [73,115], it is shown that such a hybrid trial move obeys detailed balance.

Each simulation in the CFCNPT ensemble was carried out with ${10}^5$ equilibration cycles and $4\times {10}^6$ production cycles. To facilitate the sampling of λ, a weight function ($W\left(\lambda \right)$) was used to make the sampled probability of λ flat [77,78]. During equilibration, the Wang-Landau algorithm was used to construct the weight function [137,138]. In each MC step, trial moves were selected with the following probabilities: $1\text{\%}$ volume changes, $35\text{\%}$ translations, $30\text{\%}$ rotations, $17\text{\%}$ λ changes, $8.5\text{\%}$ reinsertions, and $8.5\text{\%}$ identity changes.

3.1.1. Scaling of the Lennard–Jones and electrostatic interactions

In CFCNPT simulations, the LJ interactions of the fractional molecule were scaled as follows [77]:(10) $\begin{array}{rl}{u}_{\mathrm{LJ}}& \left(r,{\lambda }_{\mathrm{LJ}}\right)={\lambda }_{\mathrm{LJ}}4\epsilon \\ & \left(\frac{1}{{\left[\frac{1}{2}{\left(1-{\lambda }_{\mathrm{LJ}}\right)}^2+{\left(\frac{r}{\sigma }\right)}^6\right]}^2}-\frac{1}{\left[\frac{1}{2}{\left(1-{\lambda }_{\mathrm{LJ}}\right)}^2+{\left(\frac{r}{\sigma }\right)}^6\right]}\right),\end{array}$(10) ${\lambda }_{\mathrm{LJ}}$ is the coupling parameter for LJ interactions as shown in Figure 2. The Coulombic interactions for the DSF potential are scaled as(11) $\begin{array}{cc}{E}_{\mathrm{Coul}}^{\mathrm{DSF}}\left(r,{\lambda }_{\mathrm{Coul}}\right)& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{\lambda }_{\mathrm{Coul}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha \left[r_{iajb}+r_{\ast }\right]\right)}{r_{iajb}+r_{\ast }}-\frac{\mathrm{erfc}\left(\alpha \left[R_{\mathrm{c}}+r_{\ast }\right]\right)}{R_{\mathrm{c}}+r_{\ast }}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha \left[R_{\mathrm{c}}+r_{\ast }\right]\right)}{{\left(R_{\mathrm{c}}+r_{\ast }\right)}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{\left[R_{\mathrm{c}}+r_{\ast }\right]}^2\right)}{R_{\mathrm{c}}+r_{\ast }}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{\lambda }_{\mathrm{Coul}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha \left[r_{iaib}+r_{\ast }\right]\right)}{r_{iaib}+r_{\ast }}-\frac{\mathrm{erfc}\left(\alpha \left[R_{\mathrm{c}}+r_{\ast }\right]\right)}{R_{\mathrm{c}}+r_{\ast }}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{{\lambda }_{\mathrm{Coul}}{q}_{ia}{q}_{ib}}{r_{iaib}+r_{\ast }}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{\lambda }_{\mathrm{Coul}}{q}_{ia}^2,\end{array}$(11) where $r_{\ast }=A(1-{\lambda }_{\mathrm{Coul}}{)}^2$, and $A=\frac{1}{2}\text{Å}$. ${\lambda }_{\mathrm{Coul}}$ is the couplingparameter for Coulombic interactions as shown in Figure 2. By ignoring the second line in Equation (11(11) $\begin{array}{cc}{E}_{\mathrm{Coul}}^{\mathrm{DSF}}\left(r,{\lambda }_{\mathrm{Coul}}\right)& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{\lambda }_{\mathrm{Coul}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha \left[r_{iajb}+r_{\ast }\right]\right)}{r_{iajb}+r_{\ast }}-\frac{\mathrm{erfc}\left(\alpha \left[R_{\mathrm{c}}+r_{\ast }\right]\right)}{R_{\mathrm{c}}+r_{\ast }}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha \left[R_{\mathrm{c}}+r_{\ast }\right]\right)}{{\left(R_{\mathrm{c}}+r_{\ast }\right)}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{\left[R_{\mathrm{c}}+r_{\ast }\right]}^2\right)}{R_{\mathrm{c}}+r_{\ast }}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{\lambda }_{\mathrm{Coul}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha \left[r_{iaib}+r_{\ast }\right]\right)}{r_{iaib}+r_{\ast }}-\frac{\mathrm{erfc}\left(\alpha \left[R_{\mathrm{c}}+r_{\ast }\right]\right)}{R_{\mathrm{c}}+r_{\ast }}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{{\lambda }_{\mathrm{Coul}}{q}_{ia}{q}_{ib}}{r_{iaib}+r_{\ast }}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{\lambda }_{\mathrm{Coul}}{q}_{ia}^2,\end{array}$(11) ), a similar expression is obtained for the scaled Coulombic interaction of the fractional molecule using the Wolf method (Equation (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) )).

3.2. Molecular dynamics simulations

All MD simulations were performed with the LAMMPS software package [139]. Periodic boundary conditions were used. All molecules were kept rigid using the SHAKE algorithm [86]. LJ potentials were truncated and analytic tail corrections were applied to compute the energy and pressure of the system [86]. The Lorentz–Berthelot mixing rules were used for non-bonded LJ interactions [69,86]. To compute densities and enthalpies of water–methanol mixtures, MD simulations are performed in the NPT ensemble. The Nosé–Hoover thermostat and barostat are used in all MD simulations performed in this work [86]. To calculate the RDFs, the ensemble average densities obtained from NPT simulations were used to fix the densities in the NVT ensemble simulations. The temperature of the system is regulated by using the Nosé–Hoover thermostat [86]. The length of each simulation for computing thermodynamic properties, i.e. densities and enthalpies, and RDFs in the NPT and NVT ensembles are 5 ns and 10 ns, respectively. A time step of 1 fs is used to integrate the equations of motion. The specifications of the force fields used in MD simulations are the same as the specifications used in MC simulations. For the Ewald summation method, long-range electrostatic interactions are computed with a relative precision of ${10}^{-6}$ [86].

4. Results

4.1. Electrostatics

The simplest way to characterise the structure of the liquid phase is by calculating its RDF [69,86,140]. To investigate how different methods for handling the electrostatics may influence the structure of the liquid, a binary mixture of water–methanol (50–50%) was considered as a representative case. The RDF obtained using the Ewald summation, from MD simulations, was considered as a reference case. For this comparison, the TIP4P/2005 [17] and TraPPE [83] force fields were considered. The densities of the binary water–methanol mixtures were obtained from independent MC and MD simulations, based on different treatments of the long-range electrostatic interactions. The relative difference in densities obtained from MD and MC simulations was below 0.2%, see Tables S1, S5 and S15 of the Supporting Information. In Figure 3, three different RDFs are shown for water–water, water–methanol, and methanol–methanol, in a 50–50% water–methanol mixture. The oxygen atoms in water and methanol represent the molecules. Excellent agreement between the RDFs in Figure 3 shows that all electrostatic methods, both in MD and MC simulations can capture the same probability distributions for the first and second coordination shells. The small difference observed at the first coordination shells in Figure 3(a–c) is practically negligible. This can be due to the difference in the computed densities obtained from independent simulations, and/or statistical noise in the simulations [86]. As shown in Figure 3, all the RDFs have converged to unity after distance of 10 Å. For systematic discussion of the RDFs obtained from experiments and molecular simulations, the reader is referred to Refs. [39,41].

Figure 3. (Colour online) Radial distribution functions of water–methanol mixtures (50–50%), at T=298 K and P=1 bar, for: (a) water–water (b) water–methanol (c) methanol–methanol. The TIP4P/2005 [17] and TraPPE [83] force fields were used to compute the density of water–methanol mixtures in MD and MC simulations. The relative difference in densities obtained from MD and MC simulations was 0.2%. To compute the long-range electrostatic interactions, the Ewald and DSF methods were used in MD simulations. In the MC simulations, the Wolf and DSF methods (Equations (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) ) and (5(5) $\begin{array}{cc}{E}^{\mathrm{DSF}}& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{R}_{\mathrm{c}}^2\right)}{R_{\mathrm{c}}}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(5) ))) were used.

In all the RDFs obtained using the Wolf method, a numerical artifact is observed at $R_{\mathrm{c}}=14$ Å. This Gaussian-shaped artifact is most noticeable for $g\left({\mathrm{O}}_{{\mathrm{H}}_{2}O}-{\mathrm{O}}_{{\mathrm{H}}_{2}O}\right)$ in Figure 3(a), with less than 5% deviation from 1. The artificial structuring at $R_{\mathrm{c}}=14$ Å indicates a non-physical behaviour at the cutoff [104,108]. This is due to the discontinuity in intermolecular electrostatic interactions between the molecules inside and outside the cutoff sphere.

To investigate to what extent the numerical artifact of the Wolf method affects the computed thermodynamic properties, excess mixing enthalpies of water–methanol mixtures obtained from MD and MC simulations were compared based on different treatments of electrostatics: the Ewald, DSF and the Wolf methods. The water–methanol mixtures were defined based on the TIP4P/2005 [17] and TraPPE [83] force fields. The results are compared in Figure 4. From the MD simulation results, it is clear that the Ewald and DSF methods yield identical excess mixing enthalpies at different mole fractions of methanol, $x_{\mathrm{MeOH}}$. Therefore, the excess mixing enthalpies based on the DSF method, from MD simulations, were considered as reference. At $x_{\mathrm{MeOH}}=0.9$, excess mixing enthalpies obtained from MD and MC simulations are equal within statistical uncertainty. At mole fractions $x_{\mathrm{MeOH}}=0.5$ and $x_{\mathrm{MeOH}}=0.7$, the computed excess mixing enthalpies using the Wolf method, from MC simulations, are equal, within statistical uncertainty, to those obtained using the DSF method, from MD simulations. At mole fraction $x_{\mathrm{MeOH}}=0.3$, the results using the DSF method, from MC simulations are in excellent agreement with the MD simulation results. Furthermore, the MD and MC simulations yield marginally different results at $x_{\mathrm{MeOH}}=0.1$ and identical results at $x_{\mathrm{MeOH}}=0.9$. This means that there is no clear distinction between the Wolf and DSF methods in calculating excess mixing enthalpies of water–methanol mixtures.

Figure 4. (Colour online) Excess enthalpies of mixing for water–methanol mixtures based on the TIP4P/2005 [17] and TraPPE [83] force fields at T=298 K and P=1 bar. To compute the electrostatic energies, the DSF and Ewald [92] methods were used in MD simulations. In MC simulations, the Wolf and DSF methods (Equations (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) ) and (5(5) $\begin{array}{cc}{E}^{\mathrm{DSF}}& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{R}_{\mathrm{c}}^2\right)}{R_{\mathrm{c}}}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(5) )) were used to treat the electrostatic interactions. The solid line indicates experimental values for the excess mixing enthalpy [49]. Dotted lines are a guide to the eye. Error bars are smaller than symbol sizes. Raw data are listed in Tables S1, S5 and S15 of the Supporting Information.

To further investigate the numerical artifact of the Wolf method, excess chemical potentials and activity coefficients of water and methanol are computed using the Wolf and DSF methods, in MC simulations. The experimental values for activity coefficients of water and methanol are taken from Ref. [141]. Experimental excess chemical potentials of water and methanol, ${\mu }^{\mathrm{ex}}$, at different mole fractions were calculated using Equation (8(8) ${\mu }_{\mathrm{A}}=-\frac{1}{\beta }\ln \left(\frac{{⟨V⟩}_{\mathrm{CFC}\mathit{N}\mathit{P}\mathit{T}}}{{\mathrm{\Lambda }}_{\mathrm{A}}^3\left(N_{\mathrm{A}}+1\right)}\right)-\frac{1}{\beta }\ln \left(\frac{p({\lambda }_{\mathrm{A}}\uparrow 1)}{p({\lambda }_{\mathrm{A}}\downarrow 0)}\right),$(8) ). The experimental density data are provided in Ref. [56]. The excess chemical potentials of pure water and pure methanol, ${\mu }_0^{\mathrm{ex}}$, were computed from empirical equations of state [142–144], using the REFPROP software [145]. For this comparison, the TIP4P/2005 [17] and TraPPE [83] force fields were used for simulations in the CFCNPT ensemble. The computed chemical potentials and activity coefficients of water and methanol are shown in Figures 5 and 6, respectively. It is important to note that at low concentrations of water and methanol, i.e. $x_{\mathrm{MeOH}}=0.1$ and 0.9, larger error bars are observed for the excess chemical potentials. This is due to the smaller number of molecules of one of the species (either water or methanol), limiting the number of identity changes of the fractional molecule. Therefore, the sampling of the excess chemical potential becomes more difficult. Since activity coefficients are computed directly from the excess chemical potentials (Equation (9(9) ${\gamma }_{\mathrm{A}}=\frac{{\rho }_{\mathrm{A}}}{x_{\mathrm{A}}{\rho }_{0\mathrm{A}}}\cdot \exp \left[\beta \left({\mu }_{\mathrm{A}}^{\mathrm{ex}}-{\mu }_{0\mathrm{A}}^{\mathrm{ex}}\right)\right],$(9) )), only at low concentrations of water and methanol, the values of the activity coefficients display scatter. Considering larger error bars at low mole fractions, it can be seen that the results from the Wolf and DSF methods are in excellent agreement. This suggests that the artifact of the Wolf method observed in the RDFs has a minor effect on thermodynamic properties of water–methanol mixtures. For the rest of the paper, the results based on the DSF method in MC simulations are presented, and the raw data corresponding the DSF method are presented in Tables S2 to S11 of The Supporting Information. The results obtained based on the Wolf method, Equation (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) ) are not considered further in Section 4 due to this artifact. The corresponding properties computed using the Wolf method are listed in Tables S12 to S21 of the Supporting Information.

Figure 5. (Colour online) Excess chemical potentials of: (a) water, (b) methanol, with respect to the ideal gas phase, in water–methanol mixtures obtained from MC simulations in the CFCNPT ensemble [73], at T=298 K and P=1 bar. The Wolf and the DSF methods (Equations (3) and (5)) were used to calculate the electrostatic interactions. The TIP4P/2005 [17] and TraPPE [83] force fields were used. Error bars are smaller than symbol sizes. Raw data are listed in Tables S5 and S15 of the Supporting Information.

Figure 6. (Colour online) Activity coefficients of: (a) water, (b) methanol in water–methanol mixtures obtained from MC simulations in the CFCNPT ensemble, at T=298 K and P=1 bar. The Wolf [90] and the DSF [89] methods were used to calculate the electrostatic interactions. The TIP4P/2005 [17] and TraPPE [83] force fields were used. The line indicates experimental values for the activity coefficients [55]. Raw data are listed in Tables S5 and S15 of the Supporting Information.

4.2. Thermodynamic properties of water–methanol mixtures

Few molecular simulation studies investigated the properties of water–methanol mixtures using the OPLS/2016 force field [84] and popular rigid, non-polarisable force fields for water [38–40]. In some cases, non-Lorentz–Berthelot mixing rules [40] were applied to improve the predicted thermodynamic properties of water–methanol compared to experiments. To the best of our knowledge, a comparative analysis of these force field combinations to compute the activity coefficients and excess chemical potentials of water and methanol is missing. To compute the chemical potentials and activity coefficients of water and methanol, TIP3P [32], SPC/E [34], OPC [81], TIP4P/2005 [17], TIP4PEW [82] force fields for water, and the OPLS/2016 [84] and TraPPE [83] force fields for methanol are considered.

The excess mixing enthalpies of water–methanol mixtures were computed for all water–methanol force field combinations, using the DSF method (Equation (5(5) $\begin{array}{cc}{E}^{\mathrm{DSF}}& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{R}_{\mathrm{c}}^2\right)}{R_{\mathrm{c}}}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(5) )) to calculate the electrostatic interactions. The results are shown in Figure 7 as a function of $x_{\mathrm{MeOH}}$, and the raw data are listed in Tables S2 to S11 of the Supporting Information. From Figure 7, it is clear that computing the excess mixing enthalpies for water–methanol mixtures using the TraPPE force field for methanol provides considerably better results compared to the OPLS/2016 force field. The sign of the excess mixing enthalpy is predicted correctly and its parabolic shape is reproduced with partial success for all four-site water force fields [39], for $x_{\mathrm{MeOH}}>0.5$. None of the water–methanol force field combinations can precisely reproduce experimental excess enthalpies and the location of the minimum for $x_{\mathrm{MeOH}}<0.5$. Different experimental studies suggest that the unique thermodynamic behaviour of water–methanol mixtures arises from incomplete mixing of the species at molecular level [41,146,147]. Segregation of water and methanol, and formation of clusters, is reported in neutron diffraction experiments [41,146,147] and molecular simulation studies [39] for the whole concentration range. In aqueous methanol solutions, it is observed that hydrophobic methyl groups tend to cluster together, while the hydrophilic hydroxyl groups are pushed further apart and oriented more towards water-rich regions [41]. This leads to a reduction in the extent of the methanol–methanol hydrogen bonding network compared to pure methanol and an addition in the extent of the water–methanol hydrogen bonding network [41,146]. In contrast to methanol, no significant change in the local structure of water is observed in neutron diffraction studies [146]. These observations suggest that water–methanol hydrogen bonding network has a strong influence on the behaviour of the excess properties of water–methanol mixtures. Based on the significant deviation between simulation results and experiments, for $x_{\mathrm{MeOH}}<0.5$, it can be concluded that the selected force field combinations, cannot reproduce the actual clustering/orientation of methanol molecules in aqueous mixtures.

Figure 7. (Colour online) Excess mixing enthalpies for water–methanol mixtures defined by: (a) TraPPE [83] and (b) OPLS/2016 [84] force fields at T=298 K and P=1 bar. The TIP3P [33], SPC/E [34], OPC [81], TIP4P/2005 [17], TIP4P/EW [82] force fields were considered for water. The DSF method (Equation (5(5) $\begin{array}{cc}{E}^{\mathrm{DSF}}& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{R}_{\mathrm{c}}^2\right)}{R_{\mathrm{c}}}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(5) )) was used to treat the electrostatic interactions. The solid line indicates experimental values for the excess mixing enthalpy [49]. Dotted lines are a guide to the eye. Raw data are listed in Tables S2-S11 of the Supporting Information.

In this work, the TIP4P/2005-TraPPE potential outperforms the other force field combinations in predicting the excess mixing enthalpy and its shape. For water–methanol mixtures defined by the OPLS/2016 force field, the sign of the excess mixing enthalpy is not reproduced, except partially for the OPC-OPLS/2016 potential. A comparative analysis of the TraPPE and OPLS/2016 force field parameters, in Table 2, shows that partial charges on the oxygen and methyl sites are larger for the TraPPE force field. Similarly, in Table 1, it is shown that the four-site water models have larger partial charges on the oxygen or dummy site. Clearly, increasing partial charges plays an important role in producing the excess molar enthalpies of water–methanol mixtures closer to the experimental values. This is in agreement with the work of Dopazo-Paz et al. [39].

Excess chemical potentials, with respect to the ideal gas phase, and activity coefficients of water and methanol, were calculated using Equations (8(8) ${\mu }_{\mathrm{A}}=-\frac{1}{\beta }\ln \left(\frac{{⟨V⟩}_{\mathrm{CFC}\mathit{N}\mathit{P}\mathit{T}}}{{\mathrm{\Lambda }}_{\mathrm{A}}^3\left(N_{\mathrm{A}}+1\right)}\right)-\frac{1}{\beta }\ln \left(\frac{p({\lambda }_{\mathrm{A}}\uparrow 1)}{p({\lambda }_{\mathrm{A}}\downarrow 0)}\right),$(8) ) and (9(9) ${\gamma }_{\mathrm{A}}=\frac{{\rho }_{\mathrm{A}}}{x_{\mathrm{A}}{\rho }_{0\mathrm{A}}}\cdot \exp \left[\beta \left({\mu }_{\mathrm{A}}^{\mathrm{ex}}-{\mu }_{0\mathrm{A}}^{\mathrm{ex}}\right)\right],$(9) ). The results are shown in Figures 8 and 9, and the raw data are listed in Tables S2 to S11 of the Supporting Information. In Figure 8, the computed activity coefficients of water and methanol for all force field combinations were plotted as a function of $x_{\mathrm{MeOH}}$. Overall, activity coefficients of methanol obtained based on the OPLS/2016 and TraPPE force fields are in good agreement with the experiments for $x_{\mathrm{MeOH}}>0.5$. For the OPLS/2016 force field, significant deviation from experiments is observed for $x_{\mathrm{MeOH}}<0.5$. In contrast to the OPLS/2016 force field, the computed activity coefficients of methanol for the TraPPE force field are considerably closer to experimental results. For water, the predicted activity coefficients are in good agreement for $x_{\mathrm{MeOH}}<0.3$. The predicted activity coefficients of water in mixtures defined by the OPLS/2016 force field deviate significantly from experimental data for $x_{\mathrm{MeOH}}>0.3$, except for the TIP3P force field. It is clear from Figure 8 that the activity coefficients of different water models, obtained in combination with the TraPPE force field, are in better agreement with the experiments. In Figure 8, it can be seen that the TIP4P2005-TraPPE potential outperforms other force field combinations to predict the activity coefficients closest to the experimental values. Among water–methanol mixtures defined by the OPLS/2016 force field, the activity coefficients obtained from TIP3P-OPLS/2016 force fields deviate less from the experiments.

Figure 8. (Colour online) Activity coefficients of water and methanol in water–methanol mixtures for different combinations of water–methanol force fields, at T=298 K and P=1 bar. In subfigures (a) and (b); the TraPPE force field was used for methanol and in subfigures (c) and (d); the OPLS/2016 force field was used for methanol. The TIP3P [33], SPC/E [34], OPC [81], TIP4P/2005 [17], TIP4P/EW [82] force fields were considered for water. The DSF method (Equation (5(5) $\begin{array}{cc}{E}^{\mathrm{DSF}}& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{R}_{\mathrm{c}}^2\right)}{R_{\mathrm{c}}}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(5) )) was used to treat the electrostatic interactions. The solid lines indicate experimental values for the activity coefficients [55]. Dotted lines are a guide to the eye. Raw data are listed in Tables S2–S11 of the Supporting Information.

Figure 9. (Colour online) Excess chemical potentials of water and methanol, with respect to the ideal gas phase, for different combinations of water–methanol force fields, at T=298 K and P=1 bar. In subfigures (a) and (b); the TraPPE force field was used for methanol and in subfigures (c) and (d); the OPLS/2016 force field was used for methanol. The TIP3P [33], SPC/E [34], OPC [81], TIP4P/2005 [17], TIP4P/EW [82] force fields were considered for water. Error bars are smaller than symbol sizes. The solid lines indicate experimental values for the chemical potentials [142–145]. Dotted lines are a guide to the eye. Raw data are listed in Tables S2–S11 of the Supporting Information.

As the excess chemical potential and activity coefficient are related (Equation (9(9) ${\gamma }_{\mathrm{A}}=\frac{{\rho }_{\mathrm{A}}}{x_{\mathrm{A}}{\rho }_{0\mathrm{A}}}\cdot \exp \left[\beta \left({\mu }_{\mathrm{A}}^{\mathrm{ex}}-{\mu }_{0\mathrm{A}}^{\mathrm{ex}}\right)\right],$(9) )), it is important that both properties agree well with experiments. It is expected that the performance of a force field combination should be the same in predicting both the excess chemical potentials and activity coefficients of water and methanol. Since water and methanol force fields are not fitted to experimental chemical potentials, some deviation is expected depending on the force field [70]. In Figure 9, it is shown that all force fields predict the excess chemical potentials with some deviation/shift with respect to the experimental values. Since computing the activity coefficients depends only on the difference between the excess chemical potentials, see Equation (9(9) ${\gamma }_{\mathrm{A}}=\frac{{\rho }_{\mathrm{A}}}{x_{\mathrm{A}}{\rho }_{0\mathrm{A}}}\cdot \exp \left[\beta \left({\mu }_{\mathrm{A}}^{\mathrm{ex}}-{\mu }_{0\mathrm{A}}^{\mathrm{ex}}\right)\right],$(9) ), a constant shift between the predicted excess chemical potentials and the experimental data does not introduce an error in computing the activity coefficients. This is because term $({\mu }_{\mathrm{A}}^{\mathrm{ex}}-{\mu }_{0\mathrm{A}}^{\mathrm{ex}})$ in Equation (9(9) ${\gamma }_{\mathrm{A}}=\frac{{\rho }_{\mathrm{A}}}{x_{\mathrm{A}}{\rho }_{0\mathrm{A}}}\cdot \exp \left[\beta \left({\mu }_{\mathrm{A}}^{\mathrm{ex}}-{\mu }_{0\mathrm{A}}^{\mathrm{ex}}\right)\right],$(9) ) remains constant. This is especially the case for methanol when $x_{\mathrm{MeOH}}>0.5$. For $x_{\mathrm{MeOH}}<0.5$, the calculated excess chemical potentials deviate significantly from the experiments, for methanol OPLS/2016. Clearly, this leads to a considerable error in computing the activity coefficients. It can be seen in Figure 9 that the TIP4P/2005 and TIP3P show the best performance combined with the TraPPE and OPLS/2016 force fields, respectively. For pure components, the excess chemical potential of pure water predicted by the TIP3P force field has the best agreement with the empirical equation of state [144]. The excess chemical potential of water predicated by the OPC force field deviates the most from experimental data. For pure methanol, the computed excess chemical potential using the TraPPE force field agrees best with the experimental equation of state [142,143].

5. Conclusions

To predict and reproduce the thermodynamic properties of water, a large number of force fields have been published in literature. The most popular force fields for water are rigid non-polarisable potentials. The performance of these force fields depend on the experimental data used for fitting the force field parameters. These force fields can be combined with other force fields to calculate thermodynamic properties of aqueous solutions. In this work, different force field combinations for water–methanol mixtures were considered to compute excess mixing enthalpies, excess chemical potentials and activity coefficients of water and methanol. In MC simulations, spherical cutoff-based methods are computationally more efficient compared to the Ewald-type methods since computation of electrostatic interactions is reduced to the molecules inside the cutoff sphere. To investigate the accuracy of two spherical cutoff-based methods, i.e. the Wolf and DSF methods (Equations (3(3) $\begin{array}{cc}{E}^{\mathrm{Wolf}}=& \frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & +\frac{1}{2}\underset{r_{iaib}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(3) ) and (5(5) $\begin{array}{cc}{E}^{\mathrm{DSF}}& =\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{j=1}{j\ne i}}^{N_{\mathrm{m}}}\sum _{b=1}^{N_{\mathrm{a}}^j}}{q}_{ia}{q}_{jb}\left[\frac{\mathrm{erfc}\left(\alpha r_{iajb}\right)}{r_{iajb}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right.\\ & +\left.\left(\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}^2}+\frac{2\alpha }{\sqrt{\pi }}\frac{\exp \left(-{\alpha }^2{R}_{\mathrm{c}}^2\right)}{R_{\mathrm{c}}}\right)\left(r_{iajb}-R_{\mathrm{c}}\right)\right]\\ & +\frac{1}{2}\underset{r_{iajb}<R_{\mathrm{c}}}{\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}}{q}_{ia}{q}_{ib}\left[\frac{\mathrm{erfc}\left(\alpha r_{iaib}\right)}{r_{iaib}}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\right]\\ & -\frac{1}{2}\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}\sum _{\genfrac{}{}{0ex}{}{b=1}{b\ne a}}^{N_{\mathrm{a}}^i}\frac{q_{ia}{q}_{ib}}{r_{iaib}}\\ & -\left[\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{2R_{\mathrm{c}}}+\frac{\alpha }{\sqrt{\pi }}\right]\sum _{i=1}^{N_{\mathrm{m}}}\sum _{a=1}^{N_{\mathrm{a}}^i}{q}_{ia}^2.\end{array}$(5) )), RDFs and excess mixing enthalpies of aqueous solutions of methanol were computed and compared to the results obtained from the Ewald summation method. The RDFs and excess mixing enthalpies obtained from the Ewald summation and DSF methods were in excellent agreement. We observed numerical artifacts at the cutoff radius in RDFs in simulations using the Wolf method. Based on the RDFs, it can be concluded that some orientational correlation exists between between the molecules inside and outside the cutoff sphere. This may imply that the dielectric properties of water are not treated correctly [111]. However, the good agreement between the excess mixing enthalpies, activity coefficients, and chemical potentials computed based on the DSF and Wolf methods suggests that these numerical artifacts have a small effect on thermodynamic properties. By using the DSF method, we investigated the performance of the TraPPE and OPLS/2016 force fields for methanol combined with five water models: TIP3P, SPC/E, OPC, TIP4P/2005, and TIP4P/EW. For these force field combinations, we computed excess mixing enthalpies, chemical potentials and activity coefficients of water and methanol. All these properties are reproduced better by the TraPPE force field compared to the OPLS/2016 force field. The predicted properties of water–methanol mixture defined by TIP4P/2005-TraPPE force fields show the best agreement with experimental data compared to other force fields combinations.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

A. Rahbari http://orcid.org/0000-0002-6474-3028

R. Hens http://orcid.org/0000-0002-6147-0749

S. H. Jamali http://orcid.org/0000-0002-4198-0901

M. Ramdin http://orcid.org/0000-0002-8476-7035

D. Dubbeldam http://orcid.org/0000-0002-4382-1509

T. J. H. Vlugt http://orcid.org/0000-0003-3059-8712

Additional information

Funding

This work was sponsored by NWO Exacte Wetenschappen (Physical Sciences) for the use of supercomputer facilities, with financial support from the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (Netherlands Organization for Scientific Research, NWO). TJHV acknowledges NWO-CW for a VICI grant.