(472c) An Efficient Implementation of the Test Particle Insertion Method for High-Throughput Screening of Gas Absorption in Ionic Liquids

Ionic Liquids (ILs) are commonly regarded as a promising class of materials for use as selective gas entrainers due to their customizable structure and properties and their extremely low vapor pressure. Despite years of experimental data collection, only a relatively small number of systems have been experimentally investigated compared to the vast number of potential ILs. Molecular simulations can be used to guide experimental investigation toward attractive solvent candidates, but molecular simulations of ILs, especially free energy calculations required to compute gas isotherms, are computationally expensive. A simulation method that can quickly and efficiently screen Henry’s law constants of many gases in ILs could be used to direct more in-depth computational investigation of the solubility of gases in ILs beyond Henry’s constants.

The test particle insertion method, also known as the Widom insertion method, is well suited for high-throughput screening simulations of solvent-solute systems. This is because the configurations of the solvent are not perturbed by the insertion of the solute, thereby making the insertions highly parallelizable and enabling independent generation of libraries of different solvent trajectories. These trajectories can be reused to compute solubilities of many different solutes in parallel. Since Henry’s constant is obtained from the limit of infinite dilution, the Widom insertion simulations used to compute it are performed with pure solvent. This means that for a given solvent and temperature, a single pure solvent trajectory (collection of solvent configurations) can be used as the collection of solvent configurations to compute the Henry’s constants of any number of solutes. The simulation used to sample the solvent configurations can also be used to compute other, simpler thermophysical properties of the pure solvent.

Unfortunately, the larger a test particle molecule is, the harder it is to find statistically important insertion positions due to the high probability of high-energy overlap with solvent molecules. We have developed and implemented a generalized, embarrassingly parallel version of the Widom insertion method within the open source Monte Carlo simulation engine Cassandra that mitigates this challenge with configurational biasing of the test particle insertions (CBMC) as well as rapid detection and termination of CBMC trials with high-energy overlap. While this is not the first implementation of configurationally biased Widom insertions, it has the advantage of compatibility with practically any molecule type. The rapid overlap detection relies on a cell list with tiny cells that is used exclusively to check for intermolecular overlap prior to computing the energy of the trial. Intermolecular overlap is flagged when two atoms belonging to different molecules (one of which is the test particle) are closer than the overlap radius, which is either uniform for the whole simulation, specific to the atom type pair and based on a target maximum atom pair energy, or specific to the pair of in-molecule atom IDs and estimated by a prior, shorter simulation. This escapes wastefully computing energies for many positions for which the energy is unimportant.

Since molecular dynamics (MD) simulations are often faster than Monte Carlo simulations for generating decorrelated solvent configurations, especially for these bulky ILs, we also implemented a trajectory reader that allows Widom insertions to be performed on trajectory snapshots generated by other simulation engines.

We used Cassandra to perform Widom insertions on MD trajectories to study solubility of hydrofluorocarbons R-32 (difluoromethane) and R-125 (pentafluoroethane) in ionic liquids, primarily focusing on the Henry’s constants but also computing an absorption isotherm. These results have good agreement with those computed with a free energy-based MD-only method using the same force fields.

We identified more than twenty ILs for which we have good force fields and for which there is experimental Henry’s constant data for gaseous solutes of interest in the ILThermo database. Solutes of interest for this study are gases with fewer than three carbon atoms, since the test particle insertions would probably be too inefficient for larger gases. We are using these methods to conduct initial high-throughput computation of Henry’s constants of more than fifteen solutes of interest for each of these ILs. The results will be compared with experimental data where experimental data is available.