An angle dependent site-renormalized theory for the conformations of n-butane in a simple fluid

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.


Introduction
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 sitesite model fluids is the Reference Interaction Site Model (RISM) theory [1].Another case is the diagrammatically proper interaction site model (PISM) formalism of Chandler, Silbey, and Ladanyi [2] 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 [3], which are a straightforward conceptual extension [4] of the RISM equations to the 3 spatial dimensions of the orientationdependent 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 [5], molecular recognition [6], ionic solvation and transport [7].
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 effects 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 [8], Karplus [9] and Rossky [10], 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 [11] 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 [12].An especially encouraging recent development [13] 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 [14] for single-center molecular models.The method retains the general topological form of the PISM formalism of Chandler, Silbey, and Ladanyi [2].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.

A general Ornstein-Zernike method and self-consistent S bonds
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 [17].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, h vv (r) = g vv (r) − 1 is the total correlation function, g vv (r) is the radial distribution function, and c vv (r), the total correlation function, is defined by equation (1).The solvent closure equation used here is where u vv (r) is the pair potential between solvent atoms, vv (r) is the 4-point, h-bonded bridge diagram, and t vv (r) = h vv (r) − c vv (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 [18].
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 solutesolvent 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 [19], 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 c iv (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 ρ , 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 U iv (r, φ b ) 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 Ω = π 0 2π 0 2π 0 sin θdθdφdφ b .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 [20] between the Stell [21] [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 [24], or, more directly, by the equivalent method of Rushbrooke and Scoins [25] used in demonstrating the term by term equivalence of the pressure integral of g(r) to the virial series.Consider the definition [17], where A is an arbitrary constant, and c i (φ b ) 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 c i (φ b ) 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 )dφ b ≡ 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 u iv (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 selfconsistent generation of the g iv (r, φ b ) functions, we will investigate the numerical consequences of equation (11) in combination with equations (3) and (4).

A site-renormalized method for n-butane in LJ fluids
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 ≡ r 12 , 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 l 3 and l 4 are fixed, while l 5 is fixed along θ 5 but rotates through the dihedral angle, φ 5 , and, where necessary, U 1v = U 1v − u(φ 5 ).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 U iv 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 t l iv (r, φ b ) in turn has a unique topological representation, such that these terms may be separated sitewise, That is, the subset of terms in j =i t o jv (|r − l j |) are those diagrams in t l iv which are simple coordinate transforms of t o jv , 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 τ l iv (r, φ b ) 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 t o , cl;000 iv (r) = exp (−βu o iv (r) + t o iv (r)) exp(−β ūl iv (r, φ b ) + τ l;000 (r)) − 1 where the prime on the sum in the definition of hl;000 (k) means that the s 45 term is not included, c l;000 iv (r) = cl;000 (r) − Φ l;000 (r), τ l,000 iv (r) = h l,000 (r) − cl,000 (r), the renormalization functions are defined by and Again, the ij = 45 term is not included in the primed sum.The c(k), h(k) are the usual zeroth order Hankel transforms of their respective r-space functions.The s ij (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 s 45 (k) pair, the terminus site-site intramolecular distributions.We exclude the s 45 (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 [8], the coordinate transformation s(φ b ) → s 45 (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, Φ l iv , 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 s 45 (r) is uniquely defined, can be included in the OZ equations, and Φ l 4v and Φ l 5v have the same form as for the interior Φ l iv terms.

Results
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 [17], the Lennard-Jones parameters used here for CCl 4 are CCl4 /k b = 373.2K, where K is units Kelvin, and σ CCl4 = 5.27 Å, while the unique sites (C i the interior, or methylene groups, and C e the exterior or methyl groups) for the n-butane model are given by i /k b = 57.51K, σ i = 3.983 Å, e /k b = 91.19K, and σ e = 3.861 Å.For the dihedral potential of the n-butane model, we use the Scheraga parameters [16] 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 * = k b T / = 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 /k b = 91.19K 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.6K.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 CCl 4 solvent with the solute fixed in the cis and trans conformations, using standard simulation methods [26] 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 solutesolvent and solute thermodynamic results.The results using the HNCH2 approximation for the solvent are discussed in previous work [18].The excess solute-solvent internal energy, βU ex /N , is calculated from the usual molecular expression, though we take advantage of the fact that the site-site definitions (24) 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, x g /x t , 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 x g /x t indicates that the structural results for this function are in strong quantitative agreement with the simulation results of Rebertus,et. al. [16] 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 [8], 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 [8] 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, g i (r, φ b ) and g e (r, φ b ), respectively, are given in figures 3 and 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 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 10 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.

Figure 1 .
Figure 1.S(φ b), the intramolecular distribution function for the dihedral angle of nbutane in CCl4.The solid line is the simulation result of Rebertus, et.al.[16], the dashed line is the result for this work, the dotted line is the ideal gas result, and the dot-dashed line is the XRISM result.

Figure 2 .
Figure 2. S(φ b ), the intramolecular distribution function for the dihedral angle of nbutane in methane.The solid line is the ideal gas distribution, the dashed line is the integral equation result for this work.

Figure 3 .
Figure 3. gi(r, φ b ), the solute-solvent averaged radial distribution function as a function of the n-butane dihedral angle, φ b , between the interior methylene groups and the CCl4 solvent groups at ρ = 0.0063 molecules/ Å3 and T = 298 K.The distances on the x, y axes are in Å.

Figure 4 .
Figure 4. ge(r, φ b ), the solute-solvent averaged radial distribution function as a function of the n-butane dihedral angle, φ b , between the exterior methyl groups and the CCl4 solvent groups at ρ = 0.0063 molecules/ Å3 and T = 298 K.The axes are as in figure 3.

Figure 5 .
Figure5.gi(r, φ b = 0), the solute-solvent averaged radial distribution function with the solute angle in the cis conformation, between the interior methylene groups and the CCl4 solvent group.The solid line is our simulation with the solute fixed, and the dashed line is the result of the integral equation theory.

Figure 6 .
Figure 6.gi(r, φ b = π), the solute-solvent averaged radial distribution function for the trans conformer, between the methylene groups and the CCl4 solvent groups.The line types are as in figure 5.

Figure 7 .
Figure 7. ge(r, φ b = 0), the solute-solvent averaged radial distribution function for the cis conformer, between the exterior methyl groups and the CCl4 solvent groups.The solid and dashed lines are the simulation and the integral equation theory, as in figure 5.The dot-dashed line is the result from the integral equation theory holding the solute fixed in the cis conformation.

Figure 8 .
Figure 8. ge(r, φ b ), the solute-solvent averaged radial distribution function for the trans conformer, between the exterior methyl groups and the CCl4 solvent groups.The line types are as in figure 6.

Figure 9 .
Figure 9. g(r, φ b ), the solute-solvent averaged radial distribution function for the cis conformer, between the solute groups and the methane solvent groups.The solid line here is for the exterior methyl groups, the dashed line is for the interior methylene groups.

Figure 10 .
Figure 10.g(r, φ b ), the solute-solvent averaged radial distribution function for the trans conformer, between the solute groups and the methane solvent groups.The line types are as in figure 9.

Table 1 .
Various thermodynamic quantities calculated in this work.