Comparing isotropic solids and fluids at either imposed volume or pressure, we investigate various correlations of the instantaneous pressure and its ideal and excess contributions. Focusing on the compression modulus K, it is emphasized that the stress fluctuation representation of the elastic moduli may be obtained directly (without a microscopic displacement field) by comparing the stress fluctuations in conjugated ensembles. This is made manifest by computing the Rowlinson stress fluctuation expression Krow of the compression modulus for NPT-ensembles. It is shown theoretically and numerically that Krow|P = Pid(2 − Pid/K) with Pid being the ideal pressure contribution.

Among the fundamental properties of any equilibrium system are its elastic moduli characterizing the fluctuations of its extensive and/or conjugated intensive variables.1–7 The isothermal compression modulus K of an isotropic solid or fluid may thus be obtained in the NPT-ensemble at imposed particle number N, pressure P, and temperature T from the fluctuations

$\delta \hat{V}= \hat{V}- V$
$δV̂=V̂−V$ of the instantaneous volume
$\hat{V}$
$V̂$
around its mean value
$V = \langle \hat{V}\rangle$
$V=⟨V̂⟩$
according to the strain fluctuation relation2

$$K = K_\mathrm{vol}|_P\equiv \mbox{k_{\rm B}T}V / \langle \delta \hat{V}^2 \rangle |_P$$
$K=K vol |P≡kBTV/⟨δV̂2⟩|P$
(1)

with kB being Boltzmann's constant. Equivalently, K may be obtained in a canonical NVT-ensemble using Rowlinson's stress fluctuation relation4,8,9

$$K = K_\mathrm{row}|_V\equiv P + \eta _\mathrm{B}- \beta V \big\langle \delta \hat{P}_\mathrm{ex}^2 \big\rangle \big|_V$$
$K=K row |V≡P+ηB−βVδP̂ ex 2|V$
(2)

with β = 1/kBT being the inverse temperature,

$\hat{P}_\mathrm{ex}$
$P̂ ex$ the instantaneous excess pressure contribution, and ηB a Born-Lamé coefficient9 which for pairwise additive potentials becomes a simple sum of moments of derivatives of the potential with respect to the particle distance.8 In this Communication, we emphasize that the stress fluctuation representation of the elastic moduli4–6 may be obtained directly from the well-known transformation rules between conjugated ensembles.10 Focusing on the compression modulus K this is made manifest by computing Rowlinson's expression Krow deliberately for NPT-ensembles where the volume is allowed to freely fluctuate. We show that

$$K_\mathrm{row}|_P= P_\mathrm{id}\ (2 - P_\mathrm{id}/K)$$
$K row |P=P id (2−P id /K)$
(3)

with Pid being the ideal pressure contribution. We demonstrate first Eq. (3) by considering theoretically the fluctuations of the instantaneous normal pressure

$\hat{P}$
$P̂$ and its ideal and excess contributions
$\hat{P}_\mathrm{id}$
$P̂ id$
and
$\hat{P}_\mathrm{ex}$
$P̂ ex$
in both conjugated ensembles. These correlations are then checked numerically by means of Monte Carlo (MC) simulation of simple coarse-grained model systems.

As discussed in the literature,2,8,10 a simple average

$A = \langle \hat{A} \rangle$
$A=⟨Â⟩$ of an observable
${\cal A}$
$A$
does not depend on the chosen ensemble, at least not if the system is large enough (V → ∞). A correlation function
$\langle \delta \hat{A}\delta \hat{B}\rangle$
$⟨δÂδB̂⟩$
of two observables
${\cal A}$
$A$
and
${\cal B}$
$B$
may differ, however, depending on whether V or P are imposed. As shown by Lebowitz, Percus, and Verlet,10 one verifies that

$$\langle \delta \hat{A}\delta \hat{B}\rangle |_{V} = \langle \delta \hat{A}\delta \hat{B}\rangle|_{P} - \frac{K}{\beta V} \ \frac{\partial A}{\partial P} \frac{\partial B}{\partial P},$$
$⟨δÂδB̂⟩|V=⟨δÂδB̂⟩|P−KβV∂A∂P∂B∂P,$
(4)

where K = −VP/∂V2 has been used. For

$\hat{A}=\hat{B}=\hat{P}$
$Â=B̂=P̂$⁠, this implies the transformation

$$\beta V \langle \delta \hat{P}^2 \rangle |_V = \beta V \langle \delta \hat{P}^2 \rangle |_P - K,$$
$βV⟨δP̂2⟩|V=βV⟨δP̂2⟩|P−K,$
(5)

i.e., the compression modulus K may be obtained from the difference of the pressure fluctuations in both ensembles. Interestingly, the numerically more convenient Rowlinson expression Krow for NVT-ensembles can be derived from Eq. (5)9 without using a microscopic displacement field (only possible for solids)5 and avoiding the volume rescaling trick used originally for liquids.4

There is a considerable freedom for defining the instantaneous pressure

$\hat{P}= \hat{P}_\mathrm{id}+ \hat{P}_\mathrm{ex}$
$P̂=P̂ id +P̂ ex$ as long as its average P = Pid + Pex does not change.8 It is convenient for the subsequent derivations and the presented MC simulations to define the instantaneous ideal pressure
$\hat{P}_\mathrm{id}$
$P̂ id$
by

$$\hat{P}_\mathrm{id}= \mbox{k_{\rm B}T}N /\hat{V}\ \ \mbox{(MC-gauge)}$$
$P̂ id =kBTN/V̂(MC-gauge)$
(6)

and the instantaneous excess pressure

$\hat{P}_\mathrm{ex}$
$P̂ ex$ by the Kirkwood expression.8,9 Within this “MC-gauge” the thermal momentum fluctuations are assumed to be integrated out and the (effective) Hamiltonian
${\cal H}_s$
$Hs$
of a state s of the system may be written

$${\cal H}_s(\hat{V}) = - \mbox{k_{\rm B}T}N \log (\hat{V}) + {\cal U}_s(\hat{V}) + \text{const.}$$
$Hs(V̂)=−kBTNlog(V̂)+Us(V̂)+const.$
(7)

with

${\cal U}_s$
$Us$ being the total excess potential energy.

In the following, the concise notation

$\eta \equiv \beta V \langle \delta \hat{P}^2 \rangle$
$η≡βV⟨δP̂2⟩$ is used. An immediate consequence of the MC-gauge is, of course, that the fluctuations of
$\hat{P}_\mathrm{id}$
$P̂ id$
vanish for the NVT-ensemble and that, hence,

$$\eta |_V= \beta V \big\langle \delta \hat{P}_\mathrm{ex}^2 \big\rangle \big|_V.$$
$η|V=βVδP̂ ex 2|V.$
(8)

Since K > 0 for a stable system, Eq. (5) implies η|P > η|V. Depending on the disorder, η|V is, however, not a negligible contribution as assumed (implicitly) by Born.1 For solids it measures the effect of non-affinedisplacements under an imposed macroscopic strain.6,7,9

The second moment of any intensive variable computed in an ensemble, where its mean value is imposed, is obtained readily by integration by parts. Using Eq. (7), this shows that

$$\eta |_P= V \langle {\cal H}_s^{\prime \prime }(\hat{V})\rangle |_P = P_\mathrm{id}+ V \langle {\cal U}_s^{\prime \prime }(\hat{V})\rangle |_P,$$
$η|P=V⟨Hs″(V̂)⟩|P=P id +V⟨Us″(V̂)⟩|P,$
(9)

where a prime denotes a derivative with respect to the indicated variable. Albeit the indicated averages are taken over all states s and all volumes

$\hat{V}$
$V̂$ at imposed P, being simple averagesthey can also be evaluated for sufficiently large systems in the NVT-ensemble yielding identical results. Denoting the last term in Eq. (9) by ηA, ex, one can show that it is equivalent for pair interaction potentials to the already mentioned Born-Lamé coefficient:9 ηA, ex = ηB + Pex. Substituting Eqs. (8) and (9) into the Legendre transform Eq. (5), this confirms Eq. (2).

We focus now on stress fluctuations in the NPT-ensemble. By comparing with Eq. (5), one sees that if the Rowlinson formula Krow is applied at imposed P, this must yield

$$K_\mathrm{row}|_P= \beta V\big\langle \delta \hat{P}_\mathrm{id}^2 \big\rangle \big|_P + 2 \beta V\langle \delta \hat{P}_\mathrm{id}\delta \hat{P}_\mathrm{ex}\rangle |_P.$$
$K row |P=βVδP̂ id 2|P+2βV⟨δP̂ id δP̂ ex ⟩|P.$
(10)

Interestingly, Eq. (10) does not completely vanish for finite T as does the corresponding stress fluctuation expression for the shear modulus G at imposed shear stress τ.9 As a next step, we demonstrate the relations

\begin{eqnarray}\eta _\mathrm{id}|_P\equiv \beta V \big\langle \delta \hat{P}_\mathrm{id}^2 \big\rangle \big|_P & = & P_\mathrm{id}^2/K ,\end{eqnarray}
$η id |P≡βVδP̂ id 2|P=P id 2/K,$
(11)
\begin{eqnarray}\eta _\mathrm{mix}|_P\equiv \beta V \langle \delta \hat{P}_\mathrm{id}\hat{P}_\mathrm{ex}\rangle |_P & = & P_\mathrm{id}\ (1 - P_\mathrm{id}/K)\end{eqnarray}
$η mix |P≡βV⟨δP̂ id P̂ ex ⟩|P=P id (1−P id /K)$
(12)

from which Eq. (3) is then directly obtained by substitution into Eq. (10). Returning to the general transformation relation Eq. (4), we first note that the l.h.s. must vanish if at least one of the observables is a function of

$\hat{V}$
$V̂$⁠. With
$\hat{A}= \hat{B}= 1/\hat{V}$
$Â=B̂=1/V̂$
, it follows that

$$\langle \delta (1/\hat{V})^2 \rangle |_P = \frac{K}{\beta V} \ \left(\frac{\partial \langle 1/\hat{V}\rangle }{\partial P}\right)^2 \approx \frac{1}{\beta K V^3}$$
$⟨δ(1/V̂)2⟩|P=KβV∂⟨1/V̂⟩∂P2≈1βKV3$
(13)

making the steepest-descent approximation

$\langle 1/\hat{V}\rangle \approx 1/V$
$⟨1/V̂⟩≈1/V$ for simple averages and using finally V/K = −∂V/∂P.11 Remembering Eq. (6), this implies Eq. (11). With
$\hat{A}= \hat{P}_\mathrm{id}\break= k_{\rm B}TN/\hat{V}$
$Â=P̂ id =kBTN/V̂$
and
$\hat{B}= \hat{P}$
$B̂=P̂$
, one obtains similarly

$$\beta V \langle \delta \hat{P}_\mathrm{id}\delta \hat{P}\rangle |_P = \mbox{k_{\rm B}T}N K \ \frac{\partial \langle 1/\hat{V}\rangle }{\partial P} \approx P_\mathrm{id}$$
$βV⟨δP̂ id δP̂⟩|P=kBTNK∂⟨1/V̂⟩∂P≈P id$
(14)

to leading order for V → ∞. This relation implies finally the claimed correlation between ideal and excess pressure fluctuations, Eq. (12), using

$\hat{P}= \hat{P}_\mathrm{id}+ \hat{P}_\mathrm{ex}$
$P̂=P̂ id +P̂ ex$ and the already demonstrated Eq. (11).12 Please note that in Ref. 9, ideal and excess pressure fluctuations have incorrectly been assumed to be uncorrelated.

The numerical results reported here to check our predictions have been obtained by MC simulation of (i) one-dimensional (1D) nets with permanent cross-links and (ii) two-dimensional (2D) glass-forming liquids. Periodic boundary conditions are used and the pressure P is first imposed using a standard MC barostat.8,9 After equilibrating and sampling in the NPT-ensemble, the volume is fixed,

$V=\hat{V}$
$V=V̂$⁠, and various simple means and fluctuations are obtained in the NVT-ensemble at the same state point.13 For the 1D nets, we assume ideal harmonic springs, U = ∑lkl(xlRl)2/2, with xl being the distance between the connected particles, the reference length Rl of the springs being set to unity and the spring constants kl being taken randomly from a uniform distribution of half-width δk centered around a mean value also set to unity. Only simple networks are presented here where two particles i − 1 and i along the chain are connected by one spring l = i, i.e., at zero temperatures all forces fl along the chain become identical. This implies K ∼ 1/⟨1/kl⟩. The compression modulus decreases thus strongly with δk as indicated by the bold line in the inset of Fig. 1. Our 2D systems are polydisperse Lennard-Jones (pLJ) beads7 kept at a constant pressure P = 2 as described in Ref. 9.

FIG. 1.

Compression modulus K computed using the rescaled volume fluctuations Kvol|P (filled spheres), the Rowlinson stress fluctuation formula Krow|V (crosses), the difference between the total pressure fluctuations in both ensembles (squares) and the fluctuations of the inverse volume

$P_\mathrm{id}^2/\eta _\mathrm{id}|_P$
$P id 2/η id |P$ (large spheres). (Main panel) The upper data refer to systems of glass-forming 2D pLJ beads at P = 2,9 the lower data to simple 1D nets of harmonic springs at P = 0. (Inset) Compression modulus Kvs. polydispersity δk of the spring constants for 1D nets with T = 0.01 and P = 0.

FIG. 1.

Compression modulus K computed using the rescaled volume fluctuations Kvol|P (filled spheres), the Rowlinson stress fluctuation formula Krow|V (crosses), the difference between the total pressure fluctuations in both ensembles (squares) and the fluctuations of the inverse volume

$P_\mathrm{id}^2/\eta _\mathrm{id}|_P$
$P id 2/η id |P$ (large spheres). (Main panel) The upper data refer to systems of glass-forming 2D pLJ beads at P = 2,9 the lower data to simple 1D nets of harmonic springs at P = 0. (Inset) Compression modulus Kvs. polydispersity δk of the spring constants for 1D nets with T = 0.01 and P = 0.

Close modal

As shown in Fig. 1, the compression modulus K may be determined using the volume fluctuations in the NPT-ensemble, Eq. (1), or using Rowlinson's stress fluctuation formula, Eq. (2), for the NVT-ensemble. The same values of K are obtained from the Legendre transform for the pressure fluctuations, Eq. (5), and from the ideal pressure fluctuations ηid|P, Eq. (11), which thus confirms both relations. As seen in the inset of Fig. 1, the compression modulus of the 1D nets decreases with δk. Also indicated is the “affine” contribution η|P to K, measuring the mean spring constant ⟨kl⟩ = 1, and the “non-affine” contribution η|V which is seen to increase with δk. The decrease of K is thus due to the increase of the non-affine contribution.

As shown in the inset of Fig. 2, we have also checked the correlations between the ideal and the excess pressure fluctuations ηmix|P. To make both models comparable, the reduced correlation function y = ηmix|P/Pid is traced as a function of the reduced ideal pressure x = Pid/K with K as determined independently above. A perfect data collapse on the prediction y = 1 − x (bold line) is observed for all systems. The main panel of Fig. 2 shows finally the scaling of the Rowlinson formula computed in the NPT-ensemble. As before a scaling collapse of the data is achieved by plotting y = Krow|P/Pidvs.x. The bold line indicates our key prediction, Eq. (3). Interestingly, the latter result does not depend on the MC-gauge which has been used above to simplify the derivation of Eq. (3). Please note that it is not possible to increase x beyond unity for our liquid systems (KPid) and the deviations from the low-temperature plateau y = 2 are thus necessarily small. The additional Pid/K correction has thus been overlooked in our previous publication.9

FIG. 2.

Characterization of stress fluctuations in the NPT-ensemble. Large spheres refer to 2D pLJ beads for P = 2, all other symbols to 1D nets for different δk and P as indicated. (Main panel) Rescaled Rowlinson formula Krow|P/Pid as a function of the reduced ideal pressure x = Pid/K. The bold line represents our key prediction, Eq. (3), on which all data points collapse. (Inset) Similar scaling for the reduced correlation function ηmix|P/Pid confirming Eq. (12).

FIG. 2.

Characterization of stress fluctuations in the NPT-ensemble. Large spheres refer to 2D pLJ beads for P = 2, all other symbols to 1D nets for different δk and P as indicated. (Main panel) Rescaled Rowlinson formula Krow|P/Pid as a function of the reduced ideal pressure x = Pid/K. The bold line represents our key prediction, Eq. (3), on which all data points collapse. (Inset) Similar scaling for the reduced correlation function ηmix|P/Pid confirming Eq. (12).

Close modal

Emphasizing the underlying Legendre transform, Eq. (5), of the stress fluctuation formalism, we have investigated here the well-known Rowlinson stress fluctuation expression Krow for the compression modulus, Eq. (2), using deliberately the NPT-ensemble. Correcting several statements made in Ref. 9, it has been demonstrated theoretically and numerically that Eq. (3) holds. The latter result, as the other correlation relations indicated in the paper, may allow to readily calibrate (correctness, convergence, and precision) various barostats commonly used.8

P.P. thanks the Région Alsace and the IRTG Soft Matter and F.W. the DAAD for funding. We are indebted to A. Blumen (Freiburg) and A. Johner (Strasbourg) for helpful discussions.

1.
M.
Born
and
K.
Huang
,
Dynamical Theory of Crystal Lattices
(
Clarendon Press
,
Oxford
,
1954
).
2.
H. B.
Callen
,
Thermodynamics and an Introduction to Thermostatistics
(
Wiley
,
New York
,
1985
).
3.
P. M.
Chaikin
and
T. C.
Lubensky
,
Principles of Condensed Matter Physics
(
Cambridge University Press
,
1995
).
4.
J. S.
Rowlinson
,
Liquids and Liquid Mixtures
(
Butterworths Scientific Publications
,
London
,
1959
).
5.
D. R.
Squire
,
A. C.
Holt
, and
W. G.
Hoover
,
Physica
42
,
388
(
1969
).
6.
J. F.
Lutsko
,
J. Appl. Phys.
65
,
2991
(
1989
).
7.
J. P.
Wittmer
,
A.
Tanguy
,
J.-L.
Barrat
, and
L.
Lewis
,
Europhys. Lett.
57
,
423
(
2002
).
8.
M.
Allen
and
D.
Tildesley
,
Computer Simulation of Liquids
(
Oxford University Press
,
Oxford
,
1994
).
9.
J. P.
Wittmer
,
H.
Xu
,
P.
Polińska
,
F.
Weysser
, and
J.
Baschnagel
,
J. Chem. Phys.
138
,
12A533
(
2013
).
10.
J. L.
Lebowitz
,
J. K.
Percus
, and
L.
Verlet
,
Phys. Rev.
153
,
250
(
1967
).
11.
That Eq. (13) becomes exact for V → ∞ can be also seen by using that the distribution of
$\hat{V}$
$V̂$
is Gaussian.
12.
For the excess pressure fluctuations a similar relation
$$\beta V \langle \delta \hat{P}_\mathrm{ex}^2 \rangle |_P = \eta _\mathrm{A,ex}- P_\mathrm{id}(1 - P_\mathrm{id}/K)$$
$βV⟨δP̂ ex 2⟩|P=ηA, ex −P id (1−P id /K)$
is obtained using Eq. (9) together with Eqs. (11) and (12). This relation has been also confirmed numerically.
13.
In d = 1 “volume” corresponds to the linear length of the system and in d = 2 to its surface. Pressure and elastic moduli take units of energy per d-dimensional volume.