Chemical potential calculations of charged species using expanded ensemble simulations: A study to assess the approximations of recent developments in the SAFT-γ equation of state framework

Molecular simulation techniques form a powerful set of tools that can be used to assess the accuracy of fluid theories. Thermodynamic modelling of electrolyte solutions through statistical mechanics theories... [ view full abstract ]

Molecular simulation techniques form a powerful set of tools that can be used to assess the accuracy of fluid theories. Thermodynamic modelling of electrolyte solutions through statistical mechanics theories proves particularly challenging when considering organic electrolyte molecules that are characterized by complex non-spherical geometries.The SAFT-γ Mie equation of state¹ follows a group-contribution approach that allows the calculation of thermodynamic properties of heteronuclear models of molecules, composed of transferable functional groups.

Within the SAFT-γ framework the Mean Spherical Approximation (MSA) approach² is employed to account for the interactions between ions in the presence of the solvent dielectric medium, while the polar-polar and ion-polar interactions are implicitly treated with a separate Born solvation energy term³. However, the existing MSA theory is only applicable to spherical ions. Recently, Lazarou et al.⁴ have proposed two alternative approaches for calculating the Coulombic contribution to the free energy, but the accuracy of the theory needs to be assessed against molecular simulation data.

In this work we focused on calculating the excess chemical potential of model charged molecules. The simplest method to achieve this is by performing Widom insertions⁵.In this case, the molecules – or particles – of interest are inserted randomly in an equilibrated trajectory of N particles in order to sample the energy difference they cause. Nevertheless, charged species interact through long-range Coulombic interactions, and the atoms tend to form clusters in order to minimize the energy of the fluid. Thus, the N ensemble is rarely a subset the N+1 system and this leads to random insertions being inefficient or even wrong⁶. To overcome this problem it is common to introduce subdivisions ofthe region between the two states.

In our simulations the Van der Waals interactions are treated following the simple scheme of the Widom insertions, while the Coulombic contribution to the chemical potential is calculated using the Expanded Ensemble Transition Matrix (EETM) method^7,8. This method is based on an expanded ensemble simulation where the Coulombic potential of the test particles is switched on gradually with weighted Monte Carlo moves. The weights introduce a bias in order to ensure that all intermediate sub-ensembles are sampled appropriately. During the process the unweighted probabilities are averaged and saved in a “transition” matrix.

Finally,our simulation results are compared against the SAFT-γ Mie theoretical calculations in order to gain a better understanding of the effect of the approximations involved.

¹V. Papaioannou, T. Lafitte, C. Avendaño, C.S. Adjiman, G. Jackson, E.A. Müller, and A. Galindo, J. Chem.Phys. 140, 054107 (2014).

²L. Blum, Mol. Phys. 30, 1529 (1975).

³J. M. A. Schreckenberg, S. Dufal, J. Haslam, C.S. Adjiman, and G. Jackson, Mol. Phys. 112, 2339 (2014).

⁴G. Lazarou, A. Galindo, C.S. Adjiman, and G. Jackson (Unpublished).

⁵J. Vrabec, A. Lotfi, and J. Fischer, Mol. Phys. 76, 1319 (1992).

⁶B.D.A. Kofke and P.T. Cummings, Mol. Phys. 92, 973 (1997).

⁷A.S. Paluch, D.L. Mobley, and E.J. Maginn, J. Chem. Theory Comput., 7, 2910 (2011).

⁸F. Drunsel and J. Gross, Fluid Phase Equilib. 429, 205 (2016).