(191ad) Faster Computation of Virial Coefficients Via Mayer-Sampling Monte Carlo, Integral-Equation Theories, and Graphics Processing Units

The ability of a molecular model to reproduce the pressure isotherms and critical point of the fluid it describes is an important measure of its efficacy, and analysis of the model's performance can help improve its parameterization. Estimation of these quantities from the model's truncated virial equation of state (VEOS) is trivial, but computation of the high-order virial coefficients B_n required to extend the range of conditions over which the VEOS can be applied is not. Mayer-sampling Monte Carlo (MSMC) has enabled calculation of higher-order B_n than previously possible for a variety of potential models, including up to B₈ for the Lennard-Jones potential.¹ However, even higher-order B_n are required to apply the VEOS at densities approaching the critical density, and calculation of these quantities remains infeasible: the higher the order of B_n, the more operations there are per MSMC step, and the more steps must be taken to achieve a desired precision. Eventually, at an order depending upon the model's complexity, MSMC calculations become prohibitively long.

For spherically symmetric potentials like the Lennard-Jones potential, the virial coefficients B_n can be decomposed into two parts²: one which can be formulated in terms of convolution integrals and thus be readily computed by fast Fourier transforms, and another which cannot and thus must be computed by a method such as quadrature or MSMC. Among the most promising and interesting decompositions are those resulting from approximations defined by the Percus-Yevick and hypernetted-chain integral-equation theories. Because these approximations are recursive in nature, and so become successively worse as the order of B_n increases, we have developed an approach in which we parse the density-expansion coefficients of the direct correlation function c_n into approximation and correction, rather than parsing B_n. The functions c_n lie one layer of integration beneath the virial coefficients and are utilized within the Ornstein-Zernike equation and the integral-equation theory to compute c_n+1. A fully accurate c_n, rather than an approximate one, yields a more accurate approximate c_n+1 and thus reduces the work that must be done by MSMC to compute the correction.

In addition to optimizing the quantity that we compute by MSMC, we have found that we can dramatically improve the calculation speed by utilizing graphics processing units (GPUs) rather than central processing units (CPUs). GPUs have found use in computational science where problems can be broken up into very many small tasks, each needing a small amount of memory. Because MSMC only considers the interactions of a few molecules, the simulation itself can serve as the small task. Thousands of such simulations can be run simultaneously on a GPU and their results averaged together. Using GPUs, we have computed up to B₇ for the Lennard-Jones potential as well as TraPPE ethane and propane. These simulations can run hundreds of times faster than equivalent calculations on a CPU.

1. A. J. Schultz and D. A. Kofke, "Sixth, seventh and eighth virial coefficients of the Lennard-Jones model," Mol. Phys. 107 (21), 2309-2318 (2009).

2. Kippi M. Dyer, John S. Perkyns, and B. Montgomery Pettitt, "A reexamination of virial coefficients of the Lennard-Jones fluid," Theor. Chem. Acc. 105 (3), 244-251 (2001).