Bulk Viscosity of a Lennard-Jones Liquid Near the Triple Point (LAMMPS)¶
This example demonstrates how to compute the bulk viscosity of a Lennard-Jones liquid near its triple point using LAMMPS. It uses the same production runs and conventions as in the Shear viscosity example. The required theoretical background is explained the section Bulk Viscosity. In essence, it is computed in the same way as the shear viscosity, except that the isotropic pressure fluctuations are used as input.
Note
The results in this example were obtained using LAMMPS 29 Aug 2024 Update 3. Minor differences may arise when using a different version of LAMMPS, or even the same version compiled with a different compiler.
Library Imports and Matplotlib Configuration¶
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from path import Path
from yaml import safe_load
from stacie import (
UnitConfig,
compute_spectrum,
estimate_acint,
LorentzModel,
plot_fitted_spectrum,
plot_extras,
)
mpl.rc_file("matplotlibrc")
%config InlineBackend.figure_formats = ["svg"]
# You normally do not need to change this path.
# It only needs to be overridden when building the documentation.
DATA_ROOT = Path(os.getenv("DATA_ROOT", "./")) / "lammps_lj3d/sims/"
Analysis of the Production Simulations¶
The following code cells define analysis functions used below.
get_piso: Computes the isotropic pressure from the diagonal components of the time-dependent pressure tensor (\(P_{xx}\), \(P_{yy}\), and \(P_{zz}\)), as explained in the bulk viscosity theory section.estimate_bulk_viscosity: Computes the bulk viscosity, visualizes the results, and provides recommendations for data reduction (block averaging) and simulation time, as explained in the following two sections of the documentation:
def estimate_bulk_viscosity(name, p_iso, av_temperature, volume, timestep):
# Compute spectrum of the isotropic pressure fluctuations.
# Note that the Boltzmann constant is 1 in reduced LJ units.
uc = UnitConfig(
acint_fmt=".3f",
acint_symbol="η_b",
acint_unit_str="η*",
freq_unit_str="1/τ*",
time_fmt=".3f",
time_unit_str="τ*",
)
spectrum = compute_spectrum(
p_iso,
prefactors=volume / av_temperature,
timestep=timestep,
include_zero_freq=False,
)
# Estimate the bulk viscosity from the spectrum.
result = estimate_acint(spectrum, LorentzModel(), verbose=True, uc=uc)
# Plot some basic analysis figures.
plt.close(f"{name}_spectrum")
_, ax = plt.subplots(num=f"{name}_spectrum")
plot_fitted_spectrum(ax, uc, result)
plt.close(f"{name}_extras")
_, axs = plt.subplots(2, 2, num=f"{name}_extras")
plot_extras(axs, uc, result)
# Return the bulk viscosity
return result.acint
Note
When computing bulk viscosity, the include_zero_freq argument in
the compute_spectrum function must be set to False,
as the average pressure is nonzero.
This ensures the DC component is excluded from the spectrum.
See the bulk viscosity theory section for more details.
def demo_production(npart: int = 3, ntraj: int = 100):
"""
Perform the analysis of the production runs.
Parameters
----------
npart
Number of parts in the production runs.
For the initial production runs, this is 1.
ntraj
Number of trajectories in the production runs.
Returns
-------
eta_bulk
The estimated bulk viscosity.
"""
# Load the configuration from the YAML file.
with open(DATA_ROOT / "replica_0000_part_00/info.yaml") as fh:
info = safe_load(fh)
# Load trajectory data.
thermos = []
p_isos = []
for itraj in range(ntraj):
thermos.append([])
p_isos.append([])
for ipart in range(npart):
prod_dir = DATA_ROOT / f"replica_{itraj:04d}_part_{ipart:02d}/"
thermos[-1].append(np.loadtxt(prod_dir / "nve_thermo.txt"))
# The average over columns 2, 3 and 4 of parts corresponds
# to the time-dependent isotropic pressure.
p_comps = np.loadtxt(prod_dir / "nve_pressure_blav.txt")
p_isos[-1].append(p_comps[:, 1:4].mean(axis=1))
thermos = [np.concatenate(parts).T for parts in thermos]
p_iso = [np.concatenate(parts) for parts in p_isos]
# Compute the average temperature
av_temperature = np.mean([thermo[1] for thermo in thermos])
# Compute the bulk viscosity
return estimate_bulk_viscosity(
f"part{npart}",
p_iso,
av_temperature,
info["volume"],
info["timestep"] * info["block_size"],
)
eta_bulk_production = demo_production(3)
CUTOFF FREQUENCY SCAN cv2l(125%)
neff criterion fcut [1/τ*]
--------- ---------- ----------
15.0 inf 5.03e-02 (No correlation time estimate available.)
16.0 inf 5.36e-02 (No correlation time estimate available.)
17.1 inf 5.71e-02 (No correlation time estimate available.)
18.2 inf 6.07e-02 (No correlation time estimate available.)
19.4 inf 6.46e-02 (No correlation time estimate available.)
20.7 inf 6.88e-02 (No correlation time estimate available.)
22.0 inf 7.33e-02 (No correlation time estimate available.)
23.5 inf 7.80e-02 (No correlation time estimate available.)
25.0 inf 8.30e-02 (No correlation time estimate available.)
26.7 inf 8.84e-02 (No correlation time estimate available.)
28.4 inf 9.41e-02 (No correlation time estimate available.)
30.3 inf 1.00e-01 (No correlation time estimate available.)
32.3 inf 1.07e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-1.73755917e-03, 4.17952194e-01, 2.58378537e+00]))
34.4 inf 1.13e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-1.33790325e-03, 4.17271199e-01, 2.58406670e+00]))
36.7 inf 1.21e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-2.33922935e-04, 4.16524612e-01, 2.58370931e+00]))
39.1 inf 1.29e-01 (Variance of the correlation time estimate is too large.)
41.6 inf 1.37e-01 (Variance of the correlation time estimate is too large.)
44.3 inf 1.46e-01 (Variance of the correlation time estimate is too large.)
47.2 inf 1.55e-01 (Variance of the correlation time estimate is too large.)
50.3 inf 1.65e-01 (Variance of the correlation time estimate is too large.)
53.6 inf 1.76e-01 (Variance of the correlation time estimate is too large.)
57.1 inf 1.87e-01 (Variance of the correlation time estimate is too large.)
60.8 inf 1.99e-01 (Variance of the correlation time estimate is too large.)
64.7 inf 2.12e-01 (Variance of the correlation time estimate is too large.)
68.9 inf 2.26e-01 (Variance of the correlation time estimate is too large.)
73.4 inf 2.40e-01 (Variance of the correlation time estimate is too large.)
78.2 inf 2.56e-01 (Variance of the correlation time estimate is too large.)
83.3 inf 2.72e-01 (Variance of the correlation time estimate is too large.)
88.7 inf 2.90e-01 (Variance of the correlation time estimate is too large.)
94.4 inf 3.08e-01 (Variance of the correlation time estimate is too large.)
100.5 inf 3.28e-01 (Variance of the correlation time estimate is too large.)
107.1 inf 3.49e-01 (Variance of the correlation time estimate is too large.)
114.0 inf 3.72e-01 (Variance of the correlation time estimate is too large.)
121.4 inf 3.96e-01 (Variance of the correlation time estimate is too large.)
129.2 inf 4.22e-01 (Variance of the correlation time estimate is too large.)
137.6 inf 4.49e-01 (Variance of the correlation time estimate is too large.)
146.5 inf 4.78e-01 (Variance of the correlation time estimate is too large.)
156.0 3.6 5.09e-01
166.1 3.3 5.41e-01
176.8 2.9 5.76e-01
188.3 2.5 6.13e-01
200.4 2.0 6.53e-01
213.4 1.7 6.95e-01
227.2 1.6 7.40e-01
241.9 2.0 7.88e-01
257.5 2.8 8.38e-01
274.2 4.3 8.92e-01
291.9 6.6 9.50e-01
310.7 9.9 1.01e+00
330.8 14.3 1.08e+00
352.2 20.1 1.15e+00
374.9 27.6 1.22e+00
399.1 37.2 1.30e+00
424.9 49.2 1.38e+00
452.3 64.0 1.47e+00
481.5 82.1 1.57e+00
512.6 104.5 1.67e+00
Cutoff criterion exceeds incumbent + margin: 1.6 + 100.0.
INPUT TIME SERIES
Time step: 0.030 τ*
Simulation time: 300.000 τ*
Maximum degrees of freedom: 200.0
MAIN RESULTS
Autocorrelation integral: 1.191 ± 0.021 η*
Integrated correlation time: 0.050 ± 0.001 τ*
SANITY CHECKS (weighted averages over cutoff grid)
Effective number of points: 214.6 (ideally > 60)
Regression cost Z-score: 0.1 (ideally < 2)
Cutoff criterion Z-score: 2.5 (ideally < 2)
MODEL lorentz(0.1) | CUTOFF CRITERION cv2l(125%)
Number of parameters: 3
Average cutoff frequency: 6.99e-01 1/τ*
Exponential correlation time: 0.504 ± 0.040 τ*
RECOMMENDED SIMULATION SETTINGS (EXPONENTIAL CORR. TIME)
Block time: 0.158 ± 0.022 τ*
Simulation time: 31.694 ± 0.316 τ*
The cutoff criterion Z-score is relatively high, around 2. This suggests that the fits on the two halves deviate more from each other than what would be expected from the uncertainty of the spectrum. There are multiple potential explanations for this observation:
One potential explanation is that the isotropic pressure fluctuations are not perfectly Gaussian. This is expected for a Lennard-Jones fluid, as expansion of the system will result in slightly lower restoring forces than compression. Such a slightly non-Gaussian distribution of the pressure fluctuations can result in a distribution of spectral data that deviates from the Gamma distribution employed by STACIE.
Another potential cause is that there is not yet sufficient data to fix the cutoff frequency. This can be addressed by generating more trajectory data, which will make it easier to determine the suitable range of cutoff frequencies. We have not further expanded the production runs in this example, to keep the computational cost low. Furthermore, as shown in the comparison below, we already obtained a good agreement with the literature results and relatively small uncertainties.
Comparison to Literature Results¶
Computational estimates of the bulk viscosity of a Lennard-Jones fluid can be found in [MLK04b]. Since the simulation settings (\(r_\text{cut}^{*}=2.5\), \(N=1372\), \(T^*=0.722\) and \(\rho^{*}=0.8442\)) are identical to those used in this notebook, the reported values should be directly comparable.
Method |
Simulation time [τ*] |
Bulk viscosity [η\(_b\)*] |
Reference |
|---|---|---|---|
EMD NVE (STACIE) |
10800 |
1.158 ± 0.030 |
(here) extension 1 |
EMD NVE (STACIE) |
30000 |
1.191 ± 0.021 |
(here) extension 2 |
EMD NVE (Helfand-Einstein) |
300000 |
1.186 ± 0.084 |
[MLK04b] |
This comparison demonstrates that STACIE accurately reproduces bulk viscosity results while achieving lower statistical uncertainty with significantly less data than existing methods.
Note that the results for only the initial NVE production run are not included because the sanity checks indicated that the data was not sufficient.
Regression Tests¶
If you are experimenting with this notebook, you can ignore any exceptions below. The tests are only meant to pass for the notebook in its original form.
if abs(eta_bulk_production - 1.195) > 0.1:
raise ValueError(f"wrong viscosity (production): {eta_bulk_production:.3e}")