Pressure-energy correlations in liquids. I. Results from computer simulations

Physicists are familiar with the idea of thermal fluctuations in equilibrium. They also know how to extract useful information from them, using linear response theory.^1–4 These methods started with Einstein’s observation that the specific heat in the canonical ensemble is determined by the magnitude of energy fluctuations. In any thermodynamic system some variables are fixed and some fluctuate. The magnitude of the variances of the latter, as well as their mutual covariances, determines the thermodynamic “response” parameters.¹ For example, in the canonical $(NVT)$ ensemble, pressure $p$ and energy $E$ fluctuate; the magnitude of pressure fluctuations is related to the isothermal bulk modulus $KT≡−V(∂p∕∂V)T$⁠, that of the energy fluctuations to the specific heat at constant volume $cV≡T(∂S∕∂T)V$⁠, while the covariance $⟨ΔpΔE⟩$ is related⁴ to the thermal pressure coefficient $βV≡(∂p∕∂T)V$⁠. If the latter is nonzero, it implies a degree of correlation between pressure and energy fluctuations. There is no obvious reason to suspect any particularly strong correlation, and to the best of our knowledge none has ever been reported. However, in the course of investigating the physics of highly viscous liquids by computer simulation, we noted strong correlations between pressure and energy equilibrium fluctuations in several model liquids, also in the high temperature, low-viscosity state. These included the most studied of all computer liquids, the Lennard-Jones system. Surprisingly, these strong correlations survive crystallization, and they are also present in the glass phase. “Strong” here and henceforth means a correlation coefficient of order 0.9 or larger. In this paper we examine several model liquids and detail which systems exhibit strong correlations and which do not. In the companion paper⁵ (referred to as Paper II) we present a detailed analysis of the correlations for the single-component Lennard-Jones (SCLJ) system, and discuss some consequences.

Specifically, the fluctuations that are in many cases strongly correlated are those of the configurational parts of pressure and energy. The (instantaneous) pressure $p$ and energy $E$ (the parts in addition to the ideal gas terms) have contributions both from particle momenta and positions as follows:

where $ri$ is the position of the $ith$ particle. Note that $W$ has dimension energy. For a pair interaction we have

where $rij$ is the distance between particles $i$ and $j$ and $v(r)$ is the pair potential. The expression for the virial [Eq. (2)] becomes⁴ 

where for convenience we define

Figure 1(a) shows normalized instantaneous values of $p$ and $E$⁠, shifted and scaled to have zero mean and unit variance, as a function of time for the standard SCLJ liquid, while Fig. 1(b) shows the corresponding fluctuations of $W$ and $U$⁠. We quantify the degree of correlation by the standard correlation coefficient $R$⁠, defined by

Here the angle brackets ⟨ ⟩ denote thermal averages while $Δ$ denotes deviation from the average value of the given quantity. The correlation coefficient is ensemble dependent, but our main focus—the $R→1$ limit—is not. Most of the simulations reported below were carried out in the NVT ensemble. Another important characteristic quantity is the “slope” $γ$⁠, which we define as the ratio of standard deviations as follows:

Considering the “total” quantities, $p$ and $E$ [Fig. 1(a)], there is some correlation; the correlation coefficient is $0.70$⁠. For the configurational parts, $W$ and $U$⁠, on the other hand [Fig. 1(b)], the degree of correlation is much higher, $R=0.94$ in this case. Another way to exhibit the correlation is a scatter plot of $W$ against $U$⁠, as shown in Fig. 2(a).

Is this correlation surprising? Actually, there are some interatomic potentials for which there is a 100% correlation between virial and potential energy. If we have a pair potential of the form $v(r)∝r−n$⁠, an inverse power law, then $w(r)=−nv(r)$ and $Wpair=(n∕3)Upair$⁠, holds exactly. In this case the correlation is 100% and $γ=n∕3$⁠.

Conversely, suppose a system is known to be governed by a pair potential and that there is 100% correlation between $W$ and $U$⁠. We can write both $U$ and $W$ at any given time $t$ as integrals over the instantaneous radial distribution function defined⁴ as

from which

and

Here the factor of $12$ is to avoid double counting, and $ρ=N∕V$ is the number density. 100% correlation means that $W(t)=γU(t)$ holds for arbitrary $g(r,t)$ (a possible additive constant could be absorbed into the definition of $U$⁠). In particular, we could consider $g(r,t)=δ(r−r0)$⁠.⁶ Substituting this into the above expressions, the integrals go away and we find $w(r0)=−3γv(r0)$⁠. Since $r0$ was arbitrary, $v′(r)=−3γv(r)∕r$⁠, which has the solution $v(r)∝r−3γ$⁠. This connection between an inverse power-law potential and perfect correlations suggests that strong correlations can be attributed to an effective inverse power-law potential, with exponent given by three times the observed value of $γ$⁠. This will be detailed in Paper II, which shows that while this explanation is basically correct, matters are somewhat more complicated than this. For instance, the fixed volume condition, under which the strong correlations are observed, imposes certain constraints on $g(r,t)$⁠.

The celebrated Lennard-Jones potential is given⁷ by

One might think that in the case of the Lennard-Jones potential the fluctuations are dominated by the repulsive $r−12$ term, but this naive guess leads to a slope of 4, rather than the 6.3 seen in Fig. 2(a). Nevertheless the observed correlation and the above mentioned association with inverse power-law potentials suggest that an effective inverse power-law description (involving short distances), with a more careful identification of the exponent, may apply. In fact, the presence of the second, attractive, term increases the steepness of the repulsive part, thus increasing the slope of the correlation, or equivalently the effective inverse power-law exponent (Fig. 3). Note the distinction between repulsive term and repulsive part of the potential: The latter is the region where $v(r)$ has a negative slope; thus, the region $r<rm$ (⁠$rm$ being the distance where the pair potential has its minimum, $21∕6σ$ for $vLJ$⁠). This region involves both the repulsive and attractive terms (see Fig. 3, which also illustrates the approximation of the repulsive part by a power law with exponent 18). The same division was made by Weeks, Chandler, and Andersen in their noted paper of 1971,⁸ in which they showed that the thermodynamic and structural properties of the Lennard-Jones fluid were dominated by the repulsive part at high temperatures for all densities, and also at low temperatures for high densities. Ben-Amotz and Stell⁹ noted that the repulsive core of the Lennard-Jones potential may be approximated by an inverse power law with $n∼18–20$⁠. The approximation by an inverse power law may be directly checked by computing the potential and virial with an inverse power-law potential for configurations drawn from actual simulations using the Lennard-Jones potential. The agreement (apart from additive constants) is good; see Paper II.

Consider now a system with different types of pair interactions, for example, a binary Lennard-Jones system with $AA$⁠, $BB$⁠, and $AB$ interactions, or a hydrogen-bonding system modeled via both Lennard-Jones and Coulomb interactions. We can write arbitrary deviations of $U$ and $W$ from their mean values, denoted $ΔU$ and $ΔW$⁠, as a sum over types (indexed by $t$⁠; sums over pairs of a given type are implicitly understood) as follows:

Now, supposing there is near-perfect correlation for the individual terms with corresponding slopes $γt$⁠, we can rewrite $ΔW$ as

If the $γt$ are all more or less equal to a single value $γ$⁠, then this can be factored out and we get $ΔW≃γΔU$⁠. Thus the existence of different Lennard-Jones interactions in the same system does not destroy the correlation, since they have $γt∼6$⁠. On the other hand the slope for Coulomb interaction, which as an inverse power law has perfect $W$⁠, $U$-correlations, is $1∕3$⁠, so we cannot expect overall strong correlation in this case [Fig. 2(b)]. Indeed such reasoning also accounts for the reduction of correlation when the total pressure and energy are considered: $ΔE=ΔU+ΔK$⁠, while (for a large atomic system) $VΔp=γΔU+(2∕3)ΔK$⁠. The fact that $γ$ is (for the Lennard-Jones potential) quite different from $2∕3$ implies that the $p$⁠, $E$-correlation is significantly weaker than that of $W$⁠, $U$ (Fig. 1). Even in cases of unequal slopes, however, there can be circumstances under which one kind of term, and therefore one slope, dominates the fluctuations. In this case strong correlations will be observed. Examples include the high-temperature limits of hydrogen-bonded liquids (Sec. III D) and the time-averaged (total) energy and pressure in viscous liquids (Paper II).

Some of the results detailed below were published previously in letter form;¹⁰ the aim of the present contribution is to make a comprehensive report covering more systems, while Paper II contains a detailed analysis and discusses applications. In the following section, we describe the systems simulated. In Sec. III we present the results for all the systems investigated, in particular, the degree of correlation (correlation coefficient $R$⁠) and the slope. Section IV gives a summary.

A range of simulation methods, thermodynamic ensembles, and computational codes were used. One reason for this was to eliminate the possibility that strong correlations are an artifact of using a particular ensemble or code. In addition, no code can simulate the full range of systems presented. Most of the data we present are from molecular dynamics (MD) simulations, although some are from Monte Carlo¹¹ (MC) and event-driven¹² (ED) simulations. Most of the MD simulations (and of course all MC simulations), had fixed temperature (⁠$NVT$⁠), while some had fixed total energy $(NVE)$⁠. Three MD codes were used: GROMACS (GRO),^13,14 ASAP (ASAP),¹⁵ and DIGITALMATERIAL (DM).¹⁶ Homemade (HM) codes were used for the MC and ED simulations.

We now list the 13 systems studied, giving each a code name for future reference. The systems include monatomic systems interacting with pair potentials, binary atomic systems interacting with pair potentials, molecular systems consisting of Lennard-Jones particles joined rigidly together in a fixed configuration (here the Lennard-Jones interaction models the van der Waals forces), molecular systems that have Coulomb as well as Lennard-Jones interactions, metallic systems with a many-body potential, and a binary system interacting with a discontinuous “square-well” potential. Included with each system is a list specifying which simulation method(s), which ensemble(s), and which code(s) were used [semicolons separate the method(s) from the ensemble(s) and the ensemble(s) from the code(s)]. Details of the potentials are given in Appendix  A.

CU: Pure liquid Cu simulated using the many-body potential derived from effective medium theory (EMT);^17,18 (MD; $NVE$⁠; ASAP).

DB: Asymmetric “dumbbell” molecules¹⁹ consisting of two unlike Lennard-Jones spheres connected by a rigid bond; (MD; $NVT$⁠; GRO).

DZ: The potential introduced by Dzugutov²⁰ as a candidate for a monatomic glass-forming system. Its distinguishing feature is a peak in $v(r)$ around $1.5σ$⁠, after which it decays exponentially to zero at a finite value of $r$⁠; (MD; $NVT$⁠, $NVE$⁠; DM).

EXP: A system interacting with a pair potential with exponential repulsion and a van der Waals attraction; (MC; $NVT$⁠; HM).

KABLJ: The Kob–Andersen binary Lennard-Jones liquid;²¹ (MD; $NVT$⁠, $NVE$⁠; GRO, DM).

METH: The GROMOS (Ref. 22) three-site model for methanol; (MD; $NVT$⁠; GRO).

MGCU: A model of the metallic alloy $Mg85Cu15$ using an EMT-based potential;²³ (MD; $NVE$⁠; ASAP).

OTP: A three-site model of the fragile glass-former orthoterphenyl (OTP);²⁴ (MD; $NVT$⁠; GRO).

SCLJ: The standard single-component Lennard-Jonessystem with the interaction given in Eq. (11); (MD, MC; $NVT$⁠, $NVE$⁠; GRO, DM).

SPC/E: The SPC/E model of water;²⁵ (MD; $NVT$⁠; GRO).

SQW: A binary model with a pair interaction consisting of an infinitely hard core and an attractive square well;^12,26 (ED; $NVE$⁠; HM).

TIP5P: A five-site model for liquid water, which reproduces the density anomaly;²⁷ (MD; $NVT$⁠; GRO).

TOL: A seven-site united-atom model of toluene; (MD; $NVT$⁠; GRO).

The number of particles (atoms or molecules) was in the range 500–2000. Particular simulation parameters (⁠$N,ρ,T$⁠, duration of simulation) are given when appropriate in Sec. III.

SCLJ is the system we have most completely investigated. $W,U$-plots are shown for a range of thermodynamic state points in Fig. 4. Here the ensemble was $NVT$ with $N=864$⁠, and each simulation consisted of a $10ns$ run taken after $10ns$ of equilibration; for all SCLJ results so-called “argon” units are used (⁠$σ=0.34nm$⁠, $ϵ=0.997kJ∕mol$⁠). Each elongated oval in Fig. 4 is a collection of $W,U$ pairs for a given state point. Varying temperature at fixed density moves the oval parallel to itself, following an almost straight line as indicated by the dashed lines. Different densities correspond to different lines, with almost the same slope. In a system with a pure inverse power-law interaction, the correlation would be exact, and moreover the data for all densities would fall on the same straight line [see the discussion immediately after Eq. (5)]. Our data, on the other hand, show a distinct dependence on volume, but for a given volume, because of the strong correlation, the variation in $W$ is almost completely determined by that of $U$⁠.

Values of correlation coefficient $R$ for the state points of Fig. 4 are listed in Table I, along with the slope $γ$⁠. In Fig. 5 we show the temperature dependence of both $R$ and $γ$ for different densities. Lines have been drawn to indicate isochores and one isobar $(p=0)$⁠. Note that when we talk of an isobar here, we mean a set of $NVT$ ensembles with $V,T$ chosen so that the thermal average of $p$ takes on a given value, rather than fixed-pressure ensembles. This figure makes it clear that for fixed density, $R$ increases as $T$ increases, while it also increases with density for fixed temperature; the slope slowly decreases in these circumstances. In fact, it eventually reaches 4, the value expected for a pure $r−12$ interaction (e.g., at $ρ=34.6mol∕l$⁠, $T=1000K$⁠, $γ=4.61$⁠, see Ref. 10). This is consistent with the idea that the repulsive part, characterized by an effective inverse power law, dominates the fluctuations: Increasing either temperature or density increases the frequency of short-distance encounters while reducing the typical distances of such encounters. On the other hand, along an isobar, these two effects work against each other, since as $T$ increases, the density decreases. The density effect “wins” in this case, which is equivalent to a statement about the temperature and volume derivative of $R$⁠: Our simulations imply that

which is equivalent to

where $αp≡(∂V∕∂T)p∕V$ is the thermal expansivity at constant pressure and $ρ$ is the particle density. This can be recast in terms of logarithmic derivatives [valid whenever $(∂R∕∂ρ)T>0$] as follows:

Thus what we observe in the simulations, namely, that the correlation becomes stronger as temperature is reduced at fixed pressure, is a priori more to be expected when the thermal expansivity is large [since then the right hand side of Eq. (16) is large]. This has particular relevance in the context of supercooled liquids, which we discuss in Paper II, because these are usually studied by lowering temperature at fixed pressure. On the other hand if the expansivity becomes small, as for example, when a liquid passes through the glass transition, inequality (16) is a priori less likely to be satisfied. We have, in fact, observed this in a simulation of OTP: Upon cooling through the (computer) glass transition, the correlation became weaker with further lowering of temperature at constant pressure.

Remarkably, the correlation persists when the system has crystallized, as seen in the data for the highest density—the occurrence of the first-order phase transition can be inferred from the gap between the data for 90 and $110K$⁠, but the data fall on the same line above and below the transition. One would not expect the dynamical fluctuations of a crystal, which are usually assumed to be well described by a harmonic approximation, to resemble those of the high-temperature liquid. In fact, for a one-dimensional crystal of particles interacting with a harmonic potential $v(r)=12k(r−rm)2$⁠, it is easy to show (Paper II) that there is a negative correlation with slope equal to $−2∕3$⁠. To investigate whether the harmonic approximation ever becomes relevant for the correlations, we prepared a perfect fcc crystal of SCLJ particles at zero temperature and simulated it at increasing temperatures, from $0.02to90K$ in argon units, along a constant density path. The results are shown in Fig. 6. Clearly the correlation is maintained right down to zero temperature. The harmonic approximation is therefore useless for dealing with the pressure fluctuations even as $T→0$⁠, because the slope is far from $−2∕3$⁠. The reason for this is that the dominant contribution to the virial fluctuations comes from the third-order term, as shown in Paper II.

Before presenting data for all the systems studied, it is useful to see what it means for the correlation not to hold. In this subsection we consider the Dzugutov system,²⁰ whose potential contains a peak at the second-neighbor distance (Fig. 7, see Appendix  A for details) whose presence might be expected to interfere with the effectiveness of an inverse power-law description. In the next subsection we show how in a specific model of water the lack of correlation can be explicitly seen to be the result of competing interactions. Figure 8 shows $W,U$-plots for the Dzugutov system for two nearby temperatures at the same density. The ovals are much less elongated than was the case for SCLJ, indicating a significantly weaker correlation—the correlation coefficients here are 0.585 and 0.604, respectively. In Paper II it is shown explicitly that the weak correlation is due to contributions arising from the second peak. Note that the major axes of the ovals are not aligned with the line joining the state points, given by the mean values of $W$ and $U$⁠, here identifiable as the intersection of the dashed and straight lines. On the other hand, the lines of best fit from linear regression, indicated by the dashed lines in each case, do coincide with the line connecting state points. This holds generally, a fact which follows from statistical mechanics (Appendix  B). The interesting thing is rather that the major axes point in different directions, whereas in the SCLJ case they are also aligned with the state-point line. The linear-regression slope, being equal to $⟨ΔUΔW⟩∕⟨(ΔU)2⟩$⁠, treats $W$ and $U$ in an asymmetric manner by involving $⟨(ΔU)2⟩$⁠, but not $⟨(ΔW)2⟩$⁠. This is because a particular choice of independent and dependent variables is made. If instead we plotted $U$ against $W$⁠, we would expect the slope to be simply the inverse of the slope in the $W,U$-plot, but, in fact, the new slope is $⟨ΔUΔW⟩∕⟨(ΔW)2⟩$⁠. This equals the inverse of the original slope only in the case of perfect correlation, where $⟨ΔUΔW⟩2=⟨(ΔW)2⟩⟨(ΔU)2⟩$⁠. For our purposes a more symmetric estimate of the slope is desired, one which agrees with the linear-regression slope in the limit of perfect correlation. We use simply the ratio of standard deviations $⟨(ΔW)2⟩∕⟨(ΔU)2⟩$ [Eq. (7)]. This slope was used to plot the dashed line in Fig. 2(a) and the full lines in Fig. 8, where it clearly represents the orientation of the data better.²⁸ 

As we shall see in the next section, the systems that show little correlation include several, which involve both van der Waals and hydrogen bonding, modeled by Lennard-Joness and Coulomb interactions, respectively. As noted already, the latter, being a pure inverse power law $(n=1)$⁠, by itself exhibits perfect correlation with slope $γ=1∕3$⁠, while the Lennard-Jones part has near perfect correlation. However, the significant difference in slopes means that no strong correlation is seen for the full interaction. To check explicitly that this is the reason the correlation is destroyed we have calculated the correlation coefficients for the Lennard-Jones and Coulomb parts separately in a model of water. Water is chosen because the density of hydrogen bonds is quite high. Simulations were done with the TIP5P model of water,²⁷ which has the feature that the density maximum is reasonably well reproduced. This existence of the density maximum is, in fact, related to pressure and energy becoming uncorrelated, as we shall see.

Figure 9 shows the correlation coefficients and slopes for a range of temperatures; the correlation is almost nonexistent, passing through zero around where the density attains its maximum value. We have separately determined the correlation coefficient of the Lennard-Jones part of the interaction; it ranges from 0.9992 at $−25°C$ to 0.9977 at $75°C$⁠, even larger than we have seen in the SCLJ system. The reason for this is that the (attractive) Coulomb interaction forces the centers of the Lennard-Jones interaction closer together than they would be otherwise; thus, the relevant fluctuations are occurring higher up the repulsive part of the Lennard-Jones pair potential. Correspondingly the slope from this interaction ranges between 4.45 and 4.54, closer to the high-$T$⁠, high density limit of 4 than was the case for the SCLJ system. This is confirmed by inspection of the oxygen-oxygen radial distribution function in Ref. 27 where it can be seen that the first peak lies entirely to the left of the $vLJ=0$ distance $σ=0.312nm$⁠. Finally note that the near coincidence between the vanishing of the correlation coefficient and the density maximum, which is close to the experimental value of $4°C$⁠, is not accidental: The correlation coefficient is proportional to the configurational part of the thermal pressure coefficient $βV$ (Paper II). The full $βV$ vanishes exactly at the density maximum $(4°C)$⁠, but the absence of the kinetic term means that the correlation coefficient vanishes at a slightly higher temperature $(∼12°C)$⁠.

In Fig. 10 we summarize the results for the various systems. Here we plot the numerator of Eq. (6) against the denominator, including factors of $1∕(kBTV)$ in both cases to make an intensive quantity with units of pressure. Since $R$ cannot be greater than unity, no points can appear above the diagonal. Being exactly on the diagonal indicates perfect correlation $(R=1)$⁠, while being significantly below indicates poor correlation. Different types of symbols indicate different systems, as well as different densities for the same system, while symbols of the same type correspond to different temperatures.

All of the simple liquids, SCLJ, KABLJ, EXP, DB, and TOL, show strong correlations, while METH, SPC/E, and TIP5P show little correlation. Values of $R$ and $γ$ at selected state points for all systems are listed in Table II. What determines the degree of correlation? The Dzugutov and TIP5P cases have already been discussed. The poor correlation for METH and SPC/E is (presumably) because these models, like TIP5P, involve both Lennard-Jones and Coulomb interactions. In systems with multiple Lennard-Jones species, but without any Coulomb interaction, there is overall a strong correlation because the slope is almost independent of the parameters for a given kind of pair.

As the temperature is increased, the data for the most poorly correlated systems, which are all hydrogen-bonding organic molecules, slowly approach the perfect-correlation line. This is consistent with the idea that this correlation is observed when fluctuations of both $W$ and $U$ are dominated by close encounters of pairs of neighboring atoms; at higher temperature there are increasingly many such encounters, which therefore come to increasingly dominate the fluctuations. Also because the Lennard-Jones potential rises much more steeply than the Coulomb potential, the latter becomes less important as short distances dominate more. Although not obvious in the plot, we find that increasing the density at fixed temperature generally increases the degree of correlation, as found in the SCLJ case; this is also consistent with the increasing relevance of close encounters or collisions.

A system quite different from the others presented so far is the binary square-well system, SQW, with a discontinuous potential combining hard-core repulsion and a narrow attractive well (Fig. 11; see Appendix  A for details). Such a potential models attractive colloidal systems,¹² one of whose interesting features, predicted from simulations and theory, is the existence of two different glass phases, termed the “repulsive” and “attractive” glasses.²⁹ The discontinuous potential not only makes the simulations substantially different from a technical point of view, but there are also conceptual differences—in particular, the instantaneous virial is a sum of delta functions in time. The dynamical algorithm employed in the simulations is ED, where events involve a change in the relative velocity of a pair of particles due to hitting the hard-core inner wall of the potential or crossing the potential step. The algorithm must detect the next event, advance time by the appropriate amount, and adjust the velocities of all particles appropriately. There are four kinds of events (illustrated in Fig. 11): (1) “collisions,” involving the inner repulsive core; (2) “bounces,” involving bouncing off the outer (attractive) wall of the potential; (3) “trapping,” involving the separation going below the range of the outer wall; and (4) “escapes,” involving the separation increasing beyond the outer wall. To obtain meaningful values of the virial a certain amount of time averaging must be done—we can no longer consider truly instantaneous quantities. As shown in Appendix  C the time-averaged virial may be written as the following sum over all events that take place in the averaging interval $tav$⁠:

Here $re$ and $ve$ refer to the relative position and velocity for the pair of particles participating in event $e$⁠, while $Δ$ indicates the change taking place in that event. Positive contributions to $W¯$ come from collisions; the three other event types involve the outer wall, which, as is easy to see, always gives a negative contribution. The default $tav$ in the simulation was 0.05. Strong correlations emerge only at longer averaging times, however. An appropriate $tav$ may be chosen by considering the correlation functions (auto- and cross-) for virial and potential energy, plotted in Fig. 12, where the “instantaneous” values $E(t)$ and $W(t)$ correspond to averaging over 0.05 time units. We choose $tav≃τα∕10$⁠, where $τα$ is the relaxation time determined from the cross-correlation function $⟨ΔU(0)ΔW(t)⟩$⁠. Data for a few state points are shown in Table II. Remarkably, this system, so different from the continuous potential systems, also exhibits strong $W,U$-correlations, with $R=0.94$ in the case $T=1.0$⁠, $ϕ=0.595$ (something already hinted at in Fig. 12 in the fact that the curves coincide). There is a notable difference, however, compared to continuous systems: The correlation is negative.

The reason for the negative correlation is that at high density, most of the contributions to the virial are from collisions: A particle will collide with neighbor 1, recoil, and then collide with neighbor 2 before there is a chance to make a bounce event involving neighbor 1. The number of collisions that occur in a given time interval is proportional to the number of bound pairs that is exactly anticorrelated with the energy. The effective slope $γ$ has a large (negative) value of order $−50$⁠, which does not seem to depend strongly on temperature. This example is interesting because it shows that strong pressure-energy correlations can appear in a wider range of systems that might first have been guessed. Note, however, that the ordinary hard-sphere system cannot display such correlations, since potential energy does not exist as a dynamical variable in this system, i.e., it is identically zero. The idea of correlations emerging when quantities are averaged over a suitable time interval is one we shall meet again in Paper II in the context of viscous liquids.

We have demonstrated a class of model liquids whose equilibrium thermal fluctuations of virial and potential energy are strongly correlated. We have presented detailed investigations of the presence or absence of such correlations in various liquids, with extra detail presented for the standard SCLJ case. One notable aspect is how widespread these correlations are, appearing not just in Lennard-Jones potentials or in potentials that closely resemble the Lennard-Jones one, but also in systems involving many-body potentials (CU and MGCU) and discontinuous potentials (SQW). We have seen how the presence of different types of terms in the potential, such as Lennard-Jones and Coulomb interactions, spoil the correlations, even though each by itself would give rise to a strongly correlating system. Most simulations were carried out in the NVT ensemble; R is not ensemble independent, but the $R→1$ limit is.

Several of the hydrocarbon liquids studied here were simulated using simplified “united-atom” models where groups such as methyl groups or even benzene rings were represented by Lennard-Jones spheres. These could also be studied using more realistic “all-atom” models, where every atom (including hydrogen atoms) is included. It would be worth investigating whether the strength of the correlations is reduced by the associated Coulomb terms in such cases.

In Paper II we provide a detailed analysis for the SCLJ case, including consideration of contributions beyond the effective inverse power-law approximation and the $T→0$ limit of the crystal. There we also discuss some consequences, including a direct experimental verification of the correlations for supercritical argon and consequences of strong pressure-energy correlations in highly viscous liquids and biomembranes.

Useful discussions with Søren Toxværd are gratefully acknowledged. Center for viscous liquid dynamics “Glass and Time” is sponsored by The Danish National Research Foundation.

Here we give more detailed information about the interatomic potentials used. These details have been published elsewhere as indicated, except for the case of EXP and TOL.

CU: Pure liquid Cu simulated using the many-body potential derived from EMT.^17,18 This is similar to the embedded atom method of Daw and Baskes,³⁰ where the energy of a given atom $i$⁠, $Ei$⁠, is some nonlinear function (the “embedding function”) of the electron density due to the neighboring atoms. In the EMT, it is given as the energy of an atom in an equivalent reference system, the “effective medium,” plus a correction term, $Ei=EC,i(ni)+12[∑j≠ivij(rij)−∑j≠irefvij(rij)]$⁠. Specifically, the reference system is chosen as a fcc crystal of the given kind of atom, and “equivalent” means that the electron density is used to choose the lattice constant of the crystal. The correction term is an ordinary pair potential involving a simple exponential, but notice that the corresponding sum in the reference system is subtracted (guaranteeing that the correct reference energy is given when the configuration is in fact, the reference configuration). The parameters were $E0=−3.510eV$⁠, $s0=1.413Å$⁠, $V0=2.476eV$⁠, $η2=3.122Å−1$⁠, $κ=5.178$⁠, $λ=3.602$⁠, and $n0=0.0614Å−3$⁠.

DB: Asymmetric dumbbell molecules,¹⁹ consisting of two unlike Lennard-Jones spheres, labelled P and M, connected by a rigid bond. The parameters were $ϵp=5.726kJ∕mol$⁠, $σp=0.4963nm$⁠, $mp=77.106u$⁠, $ϵm=0.66944kJ∕mol$⁠, $σm=0.3910nm$⁠, and $mm=15.035u$⁠; the bond length was $d=0.29nm$⁠. Cross interactions, $ϵpm$ and $σpm$⁠, were set equal to the geometric and arithmetic means of the $p$ and $m$ parameters, respectively (Lorentz–Berthelot mixing rule³¹).

DZ: A monatomic liquid introduced by Dzugutov as a candidate for a monatomic glass-forming system.²⁰ The potential is a sum of two parts, $v(r)=v1(r)+v2(r)$⁠, with $v1(r)=A(r−m−B)exp(c∕(r−a))$ for $r<a$ and zero otherwise, and $v2(r)=Bexp(d∕(r−b))$ for $r<b$⁠, zero otherwise. The parameters are chosen to match the location and curvature of the Lennard-Jones potential: $m=16$⁠, $A=5.82$⁠, $c=1.1$⁠, $a=1.87$⁠, $B=1.28$⁠, $d=0.27$⁠, and $b=1.94$⁠.

EXP: A system interacting with a pair potential with exponential repulsion $U(r)=(ϵ∕8)[6exp(−14(r∕σ−1))−14(σ∕r)6]$⁠. Note that the attractive term has the same form as the Lennard-Jones potential.

KABLJ: The Kob–Andersen binary Lennard-Jones liquid,²¹ a mixture of two kinds of particles $A$ and $B$⁠, with $A$ making 80% of the total number. The energy and length parameters are $ϵAA=1.0$⁠, $ϵBB=0.5$⁠, $ϵAB=1.5$⁠, $σAA=1.0$⁠, $σBB=0.88$⁠, and $σAB=0.8$⁠. The masses are both equal to unity. We use the standard density $ρ=1.2σAA−3$⁠.

METH: The GROMOS three-site model for methanol.³² The sites represent the methyl $(M≡CH3)$ group $(m=15.035u)$⁠, the O atom $(m=15.999u)$⁠, and the O-bonded H atom $(m=1.008u)$⁠. The charges for Coulomb interactions are, respectively, $0.176e$⁠, $−0.574e$⁠, and $0.398e$⁠. The $M$ and O groups additionally interact via Lennard-Jones forces, with parameters $ϵMM=0.9444kJ∕mol$⁠, $ϵOO=0.8496kJ∕mol$⁠, $ϵMO=0.9770kJ∕mol$⁠, $σMM=0.3646nm$⁠, $σOO=0.2955nm$⁠, and $σMO=0.3235nm$⁠. Lennard-Jones interactions are smoothly cutoff between 0.9 and $1.1nm$⁠. The $M–O$ and O–H distances are fixed at 0.136 and $0.1nm$⁠, respectively, while the $M–O–H$ bond angle is fixed at 108.53°.

MGCU: A model of the metallic alloy $Mg85Cu15$⁠, simulated by EMT with parameters as in Ref. 23. In this potential there are seven parameters for each element. However, some of the Cu parameters were allowed to vary from their original values in the process of optimizing the potential for the Mg–Cu system. The parameters for Cu were $E0=−3.510eV$⁠, $s0=1.413Å$⁠, $V0=1.994eV$⁠, $η2=3.040Å−1$⁠, $κ=4.944$⁠, $λ=3.694$⁠, and $n0=0.0637Å−3$⁠, while those for Mg were $E0=−1.487eV$⁠, $s0=1.766Å$⁠, $V0=2.230eV$⁠, $η2=2.541Å−1$⁠, $κ=4.435$⁠, $λ=3.293$⁠, and $n0=0.0355Å−3$⁠. It should be noted that there is an error in Ref. 23: The parameter $s0$ for Cu is given in units of bohr instead of Å.

OTP: The Lewis–Wahnström three-site model of OTP (Ref. 24) consisting of three identical Lennard-Jones spheres located at the apices $A$⁠, $B$⁠, and $C$ of an isosceles triangle. Sides $AB$ and $BC$ are $0.4830nm$ long, while the $ABC$ angle is 75°. The Lennard-Jones interaction parameters are $ϵ=4.989kJ∕mol$ and $σ=0.483nm$⁠, while the mass of each sphere, not specified in Ref. 24, was taken as one-third of the mass of an OTP molecule, $m=76.768u$⁠.

SCLJ: The standard single-component Lennard-Jones system with potential given by Eq. (11).

SPC/E: The SPC/E model of water,²⁵ in which each molecule consists of three rigidly bonded point masses, with an OH distance of $0.1nm$ and the HOH angle equal to the tetrahedral angle. Charges on O and each H are equal to $−0.8476e$ and $+0.4238e$⁠, respectively. O atoms interact with each other via a Lennard-Jones potential with $ϵ=2.601kJ∕mol$ and $σ=0.3166nm$⁠.

SQW: A binary model with a pair interaction consisting of an infinitely hard core and an attractive square well:^12,26 $vij(r)=∞$⁠, $r<σij$⁠, $vij(r)=−u0$⁠, $σij<r<σij(1+ϵ)$⁠, $vij(r)=0$⁠, and $r>σij(1+ϵ)$⁠. The radius parameters are $σAA=1.2$⁠, $σBB=1$⁠, and $σAB=1.1$⁠, while $ϵ=0.03$ and $u0=1$⁠. The composition was equimolar, and the masses of both particles were equal to unity.

TIP5P: In this water model²⁷ there are five sites associated with a single water molecule. One for the O atom, one for each H, and two to locate the centers of negative charge corresponding to the electron lone pairs on the O. The OH bond length and the HOH bond angle are fixed at the gas-phase experimental values $rOH=0.09572nm$ and $θHOH=104.52°$⁠. The negative charge sites are located symmetrically along the lone-pair directions at distance $rOL=0.07nm$ and with an intervening angle $θLOL=109.47°$⁠. A charge of $+0.241e$ is located on each hydrogen site, while charges of equal magnitude and opposite sign are placed on the lone-pair sites. O atoms on different molecules interact via the Lennard-Jones potential with $σO=0.312nm$ and $ϵO=0.669kJ∕mol$⁠.

TOL: A seven site united-atom model of toluene, consisting of six “ring” C atoms and a methyl group (H atoms are not explicitly represented). In order to handle the constraints more easily, only three mass points were used; one at the ring C attached to the methyl group $(m=40.065u)$⁠, and one at each of the two “meta” C atoms $(m=26.038)$ (note that with this mass distribution, the moment of inertia is not reproduced correctly). Parameters were derived from the information in Ref. 33: $ϵring=0.4602kJ∕mol$⁠, $ϵmethyl=0.6694kJ∕mol$⁠, $σring=0.375nm$⁠, and $σmethyl=0.391nm$⁠. The Lorentz–Berthelot rule was used for cross interactions.³³ 

If $A$ is a dynamical quantity that depends only on the configurational degrees of freedom, then its average in the canonical ensemble $(NVT)$ is given by (where, for convenience, we use a discrete-state notation, with $Ai$ referring to the value of $A$ in the $ith$ microstate, etc.)

where $β=1∕kBT$ and $Q$ is the configurational partition function. Then the inverse temperature derivative of $⟨A⟩$ can be written in terms of equilibrium fluctuations as follows:

Now taking $A=W$ and $A=U$ successively we find that

This last expression is precisely the formula for the slope obtained by linear regression when plotting $W$ against $U$⁠.

Consider now volume derivatives. Because volume dependence comes in through the microstate values, $Ai$ and $Ui$⁠, and the volume derivatives of these are not necessarily related in a simple way, the corresponding expression is not as simple: The derivative of $⟨W⟩$ with respect to volume at fixed temperature is given by

where $KT$ is the isothermal bulk modulus. The derivative of $U$ can be obtained by writing pressure as the derivative of Helmholtz free energy $F$ (⁠$K$ is kinetic energy) as follows:

Then using the thermodynamic identity $(∂S∕∂V)T=(∂⟨p⟩∕∂T)V≡βV$⁠, we obtain the ratio of volume derivatives of $⟨W⟩$ and $⟨U⟩$ as follows:

which becomes $−KT∕(TβV)$ when the pressure is small compared to the bulk modulus. As discussed in Paper II, $βV$ can be expressed in terms of $⟨ΔUΔW⟩$ again, but the fluctuation expression for $KT$ is more complicated. Thus we cannot simply identify the lines of constant $T$⁠, varying $V$⁠, on a $⟨W⟩$⁠, $⟨U⟩$-plot, as we could with lines of fixed $V$⁠, varying $T$⁠, by examining the fluctuations at fixed $V,T$⁠.

Here we derive the expression for the time-averaged virial, Eq. (17), as a convenience for the reader. The idea is to replace the step $u0$ in the potential with a finite slope $u0∕δ$ over a range $δ$⁠, and take the limit $δ→0$⁠. We start by replacing a two-body interaction in three dimensions with the equivalent one-dimensional, one-body problem using the radial separation $r$ and the reduced mass $mr$⁠. Let the potential step be at $r=rs$ and define $x=r−rs$ (see Fig. 13). We consider an “escape event” over a positive step, so that an initial (relative) velocity $v0$ becomes a final velocity $v1$ and $r$ goes from a value less than $r0$ to a value greater than $r0+δ$⁠. The resulting formula also applies for the other kinds of events.

The contribution to the time integral of the virial from this event is given by

where $F$ is the (constant) force in the region $0<x<δ$ and $tδ$ is the time taken for the “particle” (the radial separation) to cross this region. The trajectory $x(t)$ is given by the standard formula for uniform acceleration

which by setting $x(tδ)=δ$ gives the following expression for $tδ$⁠:

here we have taken the negative root, appropriate for a positive $u0$ (we want the smallest positive $tδ$⁠). Returning to $Δ$⁠, it can be split into two parts as follows:

Consider the second term: Using the expression for $x(t)$ from Eq. (C2), we see that the result of the integral will involve a term proportional to $tδ2$ and one proportional to $tδ3$⁠. Using Eq. (C3) to replace $tδ∝δ$⁠, and noting the $δ$ in the denominator, the terms will have linear and quadratic dependences on $δ$⁠, respectively. Thus they will vanish in the limit $δ→0$⁠. On the other hand, the first term gives