|Home | About | Journals | Submit | Contact Us | Français|
The angular dependent site-renormalized integral equation theory is developed to compute the dihedral conformation distribution and intermolecular pair distributions of n-butane at infinite dilution in a Lennard-Jones solvent. The equations take advantage of the topological diagrammatic expansion of the full angular dependent molecular system by resumming the series in conjunction with the intramolecular degree of freedom. To first order in an angular basis set, the numerical results of these site-renormalized equations are a systematic quantitative improvement over previous methods. In particular, the thermodynamics and conformational distribution of the solute are essentially indistinguishable from simulation.
In contrast to the study of simple liquids, a significant challenge in the study of the statistical mechanics of molecular liquids is the intricate multidimensional coupling of the pair distribution functions of molecules. For single-center molecular models based upon familiar multipole methods, analytic methods for accounting for the translation plus angular dimensions of the pair distribution functions of rigid molecules have been formally realized for quite some time. Utilitarian concerns have limited our understanding of the equivalent functions for multi-center site-site molecular model fluids. In most cases we have only the study of the radially symmetric pair distribution functions. The most well known example of this class of theories of the structure and thermodynamics of site-site model fluids is the Reference Interaction Site Model (RISM) theory . Another case is the diagrammatically proper interaction site model (PISM) formalism of Chandler, Silbey, and Ladanyi 
However, with ever-increasing computational resources available, as well as continued interest for understanding the detailed solvation and conformations of large, biologically active molecules, a welcome and increasingly useful trend in recent studies has developed, in which many methods are being developed for studying higher dimension projections, or averages, of the pair functions of site-site model fluids. In particular, as demonstrated especially by a significant body of work by Professor Hirata and his co-workers, the so-called 3D-RISM theories , which are a straightforward conceptual extension  of the RISM equations to the 3 spatial dimensions of the orientation-dependent pair functions for a single, arbitrarily shaped solute, have been used in many diverse and challenging application areas, including self-consistent solvation models for electronic structure calculations , molecular recognition , ionic solvation and transport .
A related question is that of the intramolecular conformations and the structure of molecules in solution. Again, the internal dimensionality of most molecules, especially polymers and biomolecules, delineates the range of the methods developed and available to study the effectss of solvation on systems of interest. Here, as well, the RISM class of theories, and extensions to it, have been the modern tool of choice for investigating site-site model systems yielding the popular mean field electrostatics such as Poisson-Boltzmann as the zeroth order term in the expansion. While there are conceptually many such extensions, examples include the methods of Pratt , Karplus  and Rossky , which have been found to be especially useful as a tool for use in compliment to simulation studies, as well as the methods developed by Eu  for the broad study of the conformations of polymers.
The increasing use of 3D-RISM methods, as well as our own continuing interests in improving the quantitative predictions and understanding of the theory of liquids, have recently led us to reexamine the integral equation methods for site-site model systems in some detail . An especially encouraging recent development  was, through means of topological resummation or renormalization for each site in the diagrammatic expansion of the pair functions, the derivation of an exact set of integral equations for site-site models which is formally equivalent to the rotational invariant methods of Blum and Torruella  for single-center molecular models. The method retains the general topological form of the PISM formalism of Chandler, Silbey, and Ladanyi . Thus, we have an exact generalization of the diagrammatically proper site-site interaction site model theory which includes the full angular dimensional details of arbitrary site-site rigid molecules. In addition to the formal development, for diatomic models, the quantitative predictions of the theory using a minimal, isotropic basis, were a significant improvement over existing methods, and, since full molecular dimensionality is intrinsic to the theory, the 3D structural projections (as well as higher order projections) are a natural feature of the theory.
This confluence of interests, together with this special tribute to Fumio Hirata, suggests an extension to intramolecular degrees of freedom. We derive and apply the method to a simple test case with ample theoretic background and simulation results, and for which the geometry of the system is simple yet illustrative. In particular, we extend the site-renormalized methods to the study of n-butane at infinite dilution in the Lennard-Jones fluid. Below, we detail the extension of the angle dependent site-renormalized theory to n-butane in the Lennard-Jones fluid, a formally exact method for self-consistently calculating the conformation distribution of the solute, numerical results, and conclusions.
Here and throughout, we study a model system of n-butane at infinite dilution in a Lennard-Jones solvent [8,10,15,16]. Our discussion assumes throughout that n-butane is represented by four unique Lennard-Jones sites, and all intramolecular degrees of freedom are assumed rigid except for the dihedral angle. Thus, the total solute-solvent pair potential has 4 intermolecular terms, plus the dihedral angle-dependent interaction implicit in the geometry. At infinite dilution, the solvent distributions are decoupled and uniquely determined by the standard integral equation methods for a single component fluid . The solvent pair distribution functions are determined by the familiar Ornstein-Zernike(OZ) equation,
where the vv subscript denotes solvent-solvent tags, ρv is the number density of the solvent, hvv(r) = gvv(r) − 1 is the total correlation function, gvv(r) is the radial distribution function, and cvv(r), the total correlation function, is defined by equation (1). The solvent closure equation used here is
where uvv(r) is the pair potential between solvent atoms, β−1 = kbT, kb is the Boltzmann constant, T is the absolute temperature, is the 4-point, h-bonded bridge diagram, and tvv(r) = hvv(r) − cvv(r) is the indirect correlation function. This approximation scheme, labeled previously as the HNCH2 approximation, is formally exact to order ρ2, and has been found to be a useful compromise between formal exactness and numerical complexity for the Lennard-Jones system .
In order to proceed, we must first detail the angular conventions used throughout. Consider the interaction between the 4-site n-butane solute and the spherical Lennard-Jones solvent. If we choose a spherical polar coordinate system relative to the solute and fixed between an arbitrary site within the solute and the atomic solvent, there are 3 unique angles in the system, the polar and azimuthal angles of the solvent atom, and the dihedral angle of the solute. Thus, the solute-solvent functions are, e.g., functions of the set (r, b) = (r, θ, , b), where r is the scalar distance between the origin site in the solute molecule and the solvent atom, and we label the dihedral angle b and reserve as the standard azimuthal angle for the solvent. Since there are 4 sites in the solute molecule, there are 4 simple choices of origin, and thus 4 unique site-site solute-solvent distributions. These site-wise solute-solvent pair functions are determined by a simple extension of the standard molecular conventions , with the OZ equations for each site origin defined as
where the iv subscripts label a given site i on the solute and v for the solvent. Because we are at infinite dilution, and since the solvent atom is spherically symmetric, the typical orientation average associated with the 3rd molecule in the OZ equation does not appear here, though the orientation dependence is still coupled through the vector dependence of r in civ(r, b).
In order to define the closure for this system, it is first necessary to choose a normalization convention for the intramolecular potential contribution to the system. We must first define how we will include the intramolecular conformational distribution, s(b). There are at least two distinct ways of conceptualizing the problem. First, we may choose to assign the intramolecular potential exclusively to single-particle properties, akin to the separation typically invoked in inhomogeneous fluids. That is, we could separate the problem such that the 2-particle density function is defined by , with the single particle densities in turn defined self-consistently as , where ρi is the solute density. In this case, the closure does not include the dihedral potential directly, and we would have, for the hypernetted-chain (HNC) approximation,
where is the total intermolecular potential between all 4 sites and the solvent atom in the iv coordinate system, excluding the dihedral contribution. In this convention, the Ω−1 normalization in the usual site-site radial pair functions,
is an extension of the standard case for rigid molecules, with . However, we may also include the intramolecular potential directly into the closure, so that
In this case, the normalization convention is now
where u(b) is the dihedral potential, and δw(b) is the excess potential of mean force on the dihedral. This construction for Ω would be necessary to avoid over counting the intramolecular contribution to the pair functions in equation (6). It would be equivalent to choosing the 2-body density to be . These distinctions between normalization are a useful, alternate description of the correspondence  between the Stell  [equations (4) and (5)] and Wertheim [22,23] [equations (6) and (7)) conventions in reactive fluids. At infinite dilution, since the OZ equation is independent of the particular choice, we use the first convention of equations (4) and (5).
To proceed, we require a self-consistent definition for s(b). As with the related problems of inhomogeneous and reactive fluids, there are several equivalent means to determine s(b). Here, we use a topological identity which may be derived by either the functional integration methods used by Morita and Hiroike , or, more directly, by the equivalent method of Rushbrooke and Scoins  used in demonstrating the term by term equivalence of the pressure integral of g(r) to the virial series. Consider the definition ,
where A is an arbitrary constant, and is the solute single particle equivalent of the total direct correlation function defined by the functional derivative with respect to the density of the solvent,
Following Rushbrooke and Scoins, we recognize that the diagrammatic expansion of is simply related to the pair distribution function in the solute-solvent system
where r is the gradient operator. The arbitrary constant A is determined by the requirement that ∫ s(b)db 1, and we have that
The normalization constant, ω, is defined
Note that, as with other sum rules, there is a required agreement between all iv coordinate systems, i.e. s(b) must be the same regardless of the site-site coordinate system chosen to perform the integrals. This is easily satisfied formally by using the usual site-site representation
where uiv(r) is the potential between site i and the solvent atom when site i is at the origin. The site-site representation for this integral is exact in this case, due to the fact that the integrals must be independent of the choice of origin. For now, we will leave to later work a more general derivation and investigation of these equations, which should be generally applicable wherever the inhomogeneity, or, equivalently, the intramolecular potential, is a functional form sufficiently localized that the origin is freely chosen. Below, in combination with a site-renormalized approach to the self-consistent generation of the giv(r, b) functions, we will investigate the numerical consequences of equation (11) in combination with equations (3) and (4).
In our previous work, we demonstrated for molecular fluids that the exact molecular OZ and closure relationships may be easily related to the PISM formalism, which itself provides a subset of correct terms. For the purpose at hand, the chief advantage of so doing is that, if we expand the molecular site-site pair functions in a spherical basis set, then the first, isotropic basis approximation in the renormalized equations is a significant quantitative and qualitative improvement over the RISM and PISM methods. Since we consider here a molecular solute at infinite dilution in an atomic solvent, the diagrammatic complexity is sufficiently reduced so that a brief, illustrative derivation is easy to present. As mentioned above, we consider n-butane as the united atom limit, with 4 Lennard-Jones sites. In this system, we label the solvent atom as particle 2, the two interior methylene sites of the n-butane molecule as 1 and 3, and the exterior methyl sites as 4 and 5. This system has two unique site-site solute-solvent coordinate systems. In the first system we place the origin coincident with site 1, with r r12, fix the z-axis along the 1–3 bond, and fix the 1–4 bond in space at 14 = 0. In this case, the dihedral reduces to b = 5, and the total potential is defined
where l3 and l4 are fixed, while l5 is fixed along θ5 but rotates through the dihedral angle, 5, and, where necessary, . Thus, the 3 unique angles in this system are (θ, , b) = (θ2, 2, 5). The second unique site-site coordinate system is coincident with one of the exterior methyl group sites of the molecule. Given the coordinate system above for the interior sites, then, keeping the methyl site labeled 4 fixed in orientation, then for any angle 5 the exterior site-site coordinate system is given by simply translating site 4 to the origin, and similarly translating sites 1,3, and 5.
Using the interior coordinate system as an example, we recognize that the potential separates naturally into 2 contributions, such that
and we again use the labeling of Chandler, Silbey, and Ladanyi in recognition of the fact that the o, or none, functions depend upon only r, while the l, or left, functions depend upon both r and b.
Given this separation of for any of the site-site coordinate systems, the closure equations also separate naturally into o, l pairs,
and a similar separation holds for the molecular OZ, equation (3). A site-renormalized form of these equations is determined by recognizing that in turn has a unique topological representation, such that these terms may be separated sitewise,
That is, the subset of terms in are those diagrams in which are simple coordinate transforms of , and thus depend only upon the position of the j ≠ i site. Further, these diagrams are simply those which, in the isotropic basis, are generated by the S-bond convolutions in the PISM equations, by definition. The diagrams are then those which depend simultaneously upon the position of more than one of the j ≠ i sites in the solute. Recognizing this diagrammatic distinction, we re-sum the left closure and OZ equations (in the isotropic basis) with respect to to,
where the prime on the sum in the definition of l;000(k) means that the s45 term is not included,
the renormalization functions are defined by
Again, the ij = 45 term is not included in the primed sum. The (k), (k) are the usual zeroth order Hankel transforms of their respective r-space functions. The sij(k) are the intramolecular site-site distribution functions, with the usual rigid constraint definition
for all i, j pairs in the solute except for the s45(k) pair, the terminus site-site intramolecular distributions. We exclude the s45(r) functions from the renormalized OZ equation, and thus use a different definition of the Φ functions in reference to the different generating topologies. This is for the simple fact that, as discussed by Pratt and Chandler , the coordinate transformation s(b) → s45(r) is dependent upon integrating over b numerically in any case. Thus, unlike the rigid bonds, there is no numerical advantage to directly pulling these terms through the molecular equations. As such, we simply integrate b in r-space and keep the renormalization terms, , in the form appropriate. This has the additional property that all b terms are chirality-preserving. We further note that in the case that the b angle is kept fixed, and all intramolecular degrees of freedom held rigid, then s45(r) is uniquely defined, can be included in the OZ equations, and and have the same form as for the interior terms.
The potential model used here is the standard 4-site Lennard-Jones model system used by several other groups [10,16]. In combination with the standard Lorentz-Berthelot mixing rules , the Lennard-Jones parameters used here for CCl4 are εCCl4/kb = 373.2 K, where K is units Kelvin, and σCCl4 = 5.27 Å, while the unique sites (Ci the interior, or methylene groups, and Ce the exterior or methyl groups) for the n-butane model are given by εi/kb = 57.51 K, σi = 3.983 Å, εe/kb = 91.19 K, and σe = 3.861 Å. For the dihedral potential of the n-butane model, we use the Scheraga parameters  for the potential form
with torsional parameters
The solvent phase point typically investigated for this system is at T = 298 K and ρCCl4 = 0.0063 Molecules/Å3. In reduced units, this corresponds to a Lennard-Jones liquid at T* = kbT/ε = 0.796 and ρ* = ρσ3 = 0.922. For a useful contrast, both numerically and physically, we also calculate the results of using solvent parameters corresponding to methane, where we take the simple approximation that εCH4/kb = 91.19 K and σCH4 = 3.861, i.e. that to a first approximation solvent methane is equivalent to the methyl groups of the solute. For this choice of parameters, the same Lennard-Jones liquid in reduced units is equivalent to ρCH4 = 0.016 Molecules/Å3 and T = 72.6 K. We note that quantum effects for this phase point, especially with respect to the dihedral rotation, may be discernible, but, with that caveat, the choice of methane should still be reasonably illustrative for the model. Finally, in order to compare the results for the intermolecular pair functions, we ran 2 molecular dynamics simulations for the CCl4 solvent with the solute fixed in the cis and trans conformations, using standard simulation methods  and 1536 solvent atoms. Numerical results for the integral equations were calculated using standard methods on a radial grid of 2048 points and spacing dr = 0.0514 Å. The angular integrations were all calculated using a regular trapezoid rule method, and satisfactory numerical convergence for all angular integrals investigated was found using (θ, , b) = (16, 32, 256) points.
Our numerical results are summarized in table 1 and below. We report here only the solute-solvent and solute thermodynamic results. The results using the HNCH2 approximation for the solvent are discussed in previous work . The excess solute-solvent internal energy, βUex/N, is calculated from the usual molecular expression, though we take advantage of the fact that the site-site definitions
are exact here by construction. Note, the factor of 2 here for the solute-solvent energy, since there are 2 equivalent solute-solvent terms in the total excess energy for the mixture. In table 1, the cis and trans labels refer to the results from our own simulations with the solute fixed in the cis and trans conformations, respectively. The excess rotational energy of the dihedral angle is, similarly,
Finally, the equilibrium constant for the dihedral rotation, or, equivalently, the ratio of gauche to trans conformers for the system, xg/xt, is defined
We note that, in all cases, the thermodynamic results for the present theory are numerically within the error of the simulation results. The intermolecular distribution functions given below show that this result is predominately due to the very close agreement of the predicted intermolecular pair functions to the simulation results outside the first solvation shell.
The most intriguing result of this work is the intramolecular distribution function, s(b). The ratio xg/xt indicates that the structural results for this function are in strong quantitative agreement with the simulation results of Rebertus, et. al.  This is confirmed in figure 1. The distribution function is essentially indistinguishable from the simulation results, and is a significant improvement over the results from the XRISM theory. We note that the construction of equation (11) is, in an important sense, an exact definition of the basic results of the theory of Pratt and Chandler , in that the packing forces of all solvent atoms in the system on the solute molecule are aggregated in their effect upon the conformational distribution. The close agreement of the theory presented here, given the form of the integrand in equation (11), indicates that the intermediate and long-range structure of the fluid, in aggregate, is an essential part of the force influencing the conformational distribution of the molecule. We further illuminate the effects of the solvent through changing the solvent to methane, the results for which are plotted in figure 2. The effect of the smaller solvent molecule is essentially negligible, and is dominated mostly by the temperature change, though there is still a qualitative solvent effect.
To investigate the intermolecular pair functions, we follow Pratt and Chandler  and calculate the site-site solute-solvent radial projections as a function of the dihedral angle, g(r, b), defined as
The complete functions for the interior methylene and exterior methyl groups, gi(r, b) and ge(r, b), respectively, are given in figures 3 and and4.4. We compare these functions for g(r, b = 0, π) against our simulation results for the fixed solutes in figures 5–8. There are two primary results that we take from these plots. First, in the first solvation shell, the integral equation theory appears to overemphasize the solvent affect on the interior methylene groups of the butane, at the expense of the solvation of the exterior methyl groups. Second, the secondary solvation structure is extremely well-represented by the integral equation theory as constructed. In order to further illustrate the situation, in figures 7 and and8,8, we show the results of holding the solute fixed in the cis and trans conformations, respectively, in the integral equation theory. From this, we conclude that the construction of the integral equations, in particular the choice of using only the linear terms in the construction of the renormalization functions, overemphasizes the contribution of the interior methylene groups with respect to the solute-solvent interaction. Further, the sharp secondary shoulder in the exterior distribution function is partially an artifact of holding the solute fixed, since this feature is qualitatively reproduced in the theory by holding the solute fixed. The overall effect of the theory as constructed appears to be to push the first solvation shell out some distance from the exact result for the model.
It is useful to investigate whether these effects are a result of the theory, or whether the combination of model and theory is to blame. We do this by changing the solvent relative to the solute, making the solvent size and interaction equivalent to the exterior methyl groups, which is basically equivalent to making the solvent equivalent to liquid methane. The results for these calculations are shown in figures 9 and and1010 for the cis and trans conformers, respectively. The key feature here is that, even though the solute is not held fixed in these calculations, the relative size change makes the solvation structure significantly more detailed, and the theory produces these results for this system without modification. This indicates that the relative contribution of the renormalization functions in the theory is sensitive to the details of the system under investigation. Given the structure of the theory, and the results here, it would seem that this particular choice of renormalization function should be most useful for the study of solutes that are large compared to the solvent constituents.
In this work, we presented an angle dependent site-renormalized integral equation theory for calculating the intramolecular and intermolecular pair distributions of n-butane at infinite dilution in a Lennard-Jones solvent. The equations were derived from the diagrammatic expansion of the full angular dependent molecular system of equations by resumming the series self-consistently. To first order, the numerical results of these site-renormalized equations were shown to be a quantitative method for predicting the intramolecular conformational distributions and solute thermodynamics. The intermolecular pair distributions between the solute and the solvent molecules were qualitatively reasonable.
We found that the thermodynamics and conformational distributions of the solute were essentially indistinguishable from the simulation results. The intermolecular pair distributions were also shown to be in qualitative agreement with simulation, but especially better with regard to the intermediate and long-range solvation structure of the liquid. These results were for a simple test case. Further work in this area will focus on developing the analysis of intramolecular distributions for larger, more complex molecules, especially polar and charged systems, as well as the contributions from higher order terms in the angular basis set expansions.
We gratefully acknowledge the support of the several agencies responsible for funding this work. K.M.D, J.S.P., and B.M.P acknowledge the support of the Robert A. Welch Foundation and N.I.H. The work of G.S. was supported by the Division of Chemical Sciences, Office of Basic Energy Sciences, Office of Energy Research, U.S. Department of Energy. The work of K.M.D. done at Stony Brook was supported by the National Science Foundation NIRT Award 29963.