5614
Ind. Eng. Chem. Res. 2006, 45, 5614-5618
Do Inverse Monte Carlo Algorithms Yield Thermodynamically Consistent Interaction Potentials? Sandeep Jain, Shekhar Garde, and Sanat K. Kumar* The Howard P. Isermann Department of Chemical and Biological Engineering, Rensselaer Polytechnic Institute, Troy, New York 12180
We numerically verify a statistical mechanics theorem which shows that there is a one-to-one equivalence between the structure of a liquid (i.e., the pair correlation function) and its pairwise additive intermolecular potential. Specifically, we show for three systems interacting with simple spherically symmetric pairwise additive potentials that inverse Monte Carlo (IMC) simulations can obtain the underlying potentials by only using the target pair correlation functions. The convergence of potentials obtained by the standard IMC procedure is, however, extremely slow. Interestingly, we find that the repulsive part of the potential converges rapidly, consistent with the well-accepted notion that it essentially determines the structure of condensed liquids. We show that additional information about the system, such as thermodynamic properties (e.g., average energy and or pressure) can be included in a modified IMC procedure. Because internal energy and pressure are primarily sensitive to the attractive part of the potential, the convergence to the true potential is improved by an order of magnitude. Although the improved convergence is a technical advance, no new information is obtained on the final converged potential by this approach, as expected by the Henderson theorem. 1. Introduction The development of potentials (or force fields) for the description of molecular systems is crucial if simulations need to be employed as predictive tools. For computationally efficient studies of molecular systems, the form of the Hamiltonian is empirically chosen. Frequently, it is the simplest possible form that can reproduce as many structural, dynamic, and thermodynamic properties of the system. The parameters of the force field are, thus, chosen to “best fit” available experimental data. In this sense, force fields are “effective” potentials that attempt to represent many body systems governed by nonpairwise additive interactions by simpler pairwise additive (frequently spherically symmetric) potentials which are useful over a broad range of phase space. To achieve further increases in computational efficiency, there has been a recent drive toward developing coarse-grained descriptions, where “nonessential” parts of the moieties of interest are not modeled. The interaction potentials between the resulting coarse-grained beads are then chosen to accurately represent the system structure when it is examined at the same level of coarse graining.1-12 Since only the system structure is reproduced by the coarse-grained model, it then follows that, in general, the thermodynamics (e.g., pressure-volumetemperature (PVT) equation of state) of coarse-grained models do not track that of the underlying system.13 Corrections can be made locally to the coarse-grained potentials to fix, for example, pressure or other thermodynamic properties, but this then represents a second level of model fitting that is statepoint dependent. The ideal goal in this context would be to develop a self-consistent procedure to generate a transferable potential that reproduces the structure and the thermodynamics of a given system, at any arbitrary level of coarse graining. To date, the development of coarse-grained potentials is generally achieved through the use of reverse or inverse Boltzmann methods.7,8,14 In a typical algorithm, a trial potential * Corresponding author. E-mail:
[email protected].
is used to generate the structure of a coarse-grained system. Comparison of this structure to the actual structure of the fluid of interest (in its coarse-grained form) then allows one to obtain better estimates of the potential. This procedure is then repeated till the simulated structure converges to the structure of the parent fluid in its coarse-grained form. Similar methods are also applied when structural information is available from experiments (e.g., neutron or other scattering methods).7 Several important fundamental questions arise during such a development, especially regarding the uniqueness and convergence of the obtained potentials. In this context, it is worthwhile to go back to the simplest systems and study these issues and the factors that control them, especially the convergence to a target potential. By now, several workers have shown that, at a given state point, a pairwise additive potential is uniquely linked to its pair correlation function.9,14 In other words, knowledge of the radial distribution function uniquely defines the underlying pair potential. Other works, notably by Gray and Gubbins15 and Zwicker and Lovett,9 assert that the specification of a two-body density (along with proper closures for integral equations) completely specify the thermodynamic state of the system. Thus, in principle, the multibody Hamiltonian can be completely specified by knowledge of the one particle and the pair density functions. While this statement of uniqueness is, thus, reassuring, nevertheless, implementation of these ideas require closure relationships, which themselves depend on the form of the higher-order multibody interactions. The major point being made here is that uniqueness relationships between liquid structure (and thermodynamics) and the system Hamiltonian have been proven. However, the relevance of these formal derivations to the problem of uniquely determining the effective interactions between coarse-grained potentials in a simulation is an issue of considerable interest to us and to the larger community. Precisely to this point, the recent calculations by Rutledge for a simple pairwise additive Lennard-Jones (LJ) system showed that, while the radial distribution function converges rather quickly in the inverse monte carlo (IMC) method, the
10.1021/ie060042h CCC: $33.50 © 2006 American Chemical Society Published on Web 07/12/2006
Ind. Eng. Chem. Res., Vol. 45, No. 16, 2006 5615
exact form of the potential was not reproduced over the simulation times employed.13 This work highlights the difficulty in computationally realizing Henderson’s theorem.14 In light of this discussion, and seeking to critically evaluate the numerical realization of Henderson’s theorem,14 here we examine the problem of reproducing three different spherically symmetric pairwise additive potentials. For these systems, we obtain target radial distribution functions from independent, long canonical ensemble Metropolis Monte Carlo simulations. We then perform IMC simulations to verify Henderson’s theorem as it applies to numerical simulations with noise. Specifically, we focus on the convergence of the potential obtained from the IMC procedure to the expected underlying potential. We find that the right form of the potential is reproduced in each case, even though the convergence is unexpectedly sluggish in cases where only structural information, i.e., g(r), is included. We then present a physically motivated procedure to speed up the convergence by systematically adding thermodynamic information to the standard structure-based IMC procedure. Specifically, knowledge of the average internal energy or the system pressure can be included in the algorithm, and the inclusion of this additional information speeds up the convergence of the obtained potential. It is important to reiterate that we recover the input potentials in all cases, even though the rate of convergence is strongly dependent on the amount of information that is used in the IMC procedure. We discuss the physical origins of this speedup and provide recommendations toward further improvement of our methods to the development of “effective” coarse-grained potentials. 2. Methods 2.A. Obtaining Target Pair Distribution Functions and Thermodynamic Properties. Canonical Monte Carlo (MC) simulations16 were carried out for three different systems of particles interacting with spherically symmetric, pairwise additive potentials. The potentials used in three different cases are as follows: (a) a shifted-force Lennard-Jones (LJ)16 potential with a cutoff radius of Rcut ) 2.5σ, (b) a Weeks-ChandlerAnderson potential corresponding to that LJ potential,17 and (c) a square-well potential with a well width of 0.25σ and a well depth of . All systems were simulated at a reduced density of Fσ3 ) 0.8 and a reduced temperature of kT/ ) 1.2. These simulations were used to extract the pair distribution functions of interest in each case. In addition, the internal energy and pressure were also tabulated. We used 100 particles and conducted 1 000 000 cycles to obtain structural and thermodynamic information from these simulations. (Each cycle corresponds to an attempted move of each particle in the system.) 2.B. Inverse Monte Carlo Simulations. Given target distribution functions, we need to choose the size of the system as well as the cutoff radius for the potential to be obtained. Whereas a detailed study of the dependence of the final results on both of these factors can be made, for simplicity, here we use the same box size and the potential cutoff as those used in the target systems. (See more on this topic below.) Further, the temperature and density in the IMC simulations were identical to those of the target systemssthis is a requirement, since the theorem of Henderson14 only holds when the original system and the model are at the same thermodynamic state. In our first implementation, which is the same as the standard IMC procedure, we only used the target radial distribution function, gtarget(r), as input at a given thermodynamic state point. The potential of mean force [ui)0(r) ) -kT ln g(r)] was taken
as the starting guess potential for each system. The trial potential was updated every 10 000 MC cycles using the prescription
(
ui+1(r) ) ui(r) + fkBT ln
gi(r)
gtarget(r)
)
(1)
gi(r) is the radial distribution function obtained during the ith step of the update, gtarget(r) is the target radial distribution function, kBT is the product of Boltzmann’s constant and the absolute temperature, and f is a numerical factor (
