Ionic Electrical Conductivity of Molten Sodium Chloride at 1100 K (OpenMM)¶
Warning
This example notebook is work in progress. There are still some issues with the MD results obtained with OpenMM, which are discussed in the notebook.
This notebook shows how to post-process trajectories from OpenMM simulations
to calculate the ionic electrical conductivity.
The OpenMM trajectories are converted to NPZ files within the Jupyter Notebooks of the simulation,
making the approach here easily adaptable to other codes or physical systems.
All OpenMM simulation notebooks can be found in the directory docs/data/openmm_salt
in STACIE’s source repository.
The required theoretical background is explained the
Ionic Electrical Conductivity section.
The MD simulations are performed using the Born-Huggins-Mayer-Tosi-Fumi potential,
which is a popular choice for molten salts. [TF64]
This potential does not use mixing rules and it is not natively implemented in OpenMM,
but it can be incorporated using the CustomNonbondedForce and some creativity,
see docs/data/openmm_salt/bhmft.py in the Git repository.
The molten salt was simulated with a 3D periodic box of 1728 ions (864 Na\(^+\) and 864 Cl\(^-\)).
The time step in all simulations was 5 fs.
Following the Recommendations for MD Simulations, an initial block size of 10 steps (50 fs) was used. The initial production runs consist of 800 blocks, i.e. 8000 steps (40 ps) of simulation time.
The spectrum of the charge current has a minimum at zero frequency, where its value is the ionic conductivity. At higher frequencies, the spectrum rises, because the ions are temporarily confined by their neighbors and oscillate in place, while on longer timescales they escape and diffuse. The Debye model describes the onset of this rise and is used throughout this notebook. It has \(P=3\) parameters. Only the amount of data differs between the analyses below: the same 100 trajectories are analyzed after one, two and three simulation stages.
In total 100 NVE production runs were performed. For each run, the system was first equilibrated in the NVT and later NPT ensemble. According to the section How to Prepare Sufficient Inputs for STACIE?, 100 runs should be sufficient to obtain a relative error on the ionic conductivity of about 1%:
where \(P\) is the number of parameters in the model and \(M = 100 \times 3\) is the number of independent input sequences. (100 trajectories, 3 Cartesian components of the charge current per trajectory) The guideline of \(400 P\) blocks would call for 1200 blocks with this model, somewhat more than the 800 blocks of the initial production runs.
Note
The results in this example were obtained using OpenMM 8.2.0. Minor differences may arise when using a different version of OpenMM, or even the same version compiled with a different compiler.
Library Imports and Configuration¶
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from path import Path
import scipy.constants as sc
from scipy.stats import chi2
from stacie import (
DebyeModel,
UnitConfig,
compute_spectrum,
estimate_acint,
plot_fitted_spectrum,
plot_extras,
)
from utils import plot_instantaneous_percentiles, plot_cumulative_temperature_histogram
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", "./")) / "openmm_salt/output/"
Analysis of the NpT Equilibration Runs¶
To validate that the equilibration runs have reached to proper temperature distribution, the following cell implements a plot of the percentiles (over the 100 trajectories) of a thermodynamic quantity (temperature or volume).
def plot_openmm_percentiles(
ensemble: str,
field: str,
unitstr: str,
unit: float = 1,
npart: int = 1,
ntraj: int = 100,
expected: None = None,
ymin: float | None = None,
ymax: float | None = None,
):
"""Plot the temperature of the NpT equilibration runs."""
time = None
natom = None
sequences = []
time = None
for itraj in range(ntraj):
row = []
if itraj == 0:
time = []
for ipart in range(npart):
path_npz = DATA_ROOT / f"sim{itraj:04d}_part{ipart:02d}_{ensemble}_traj.npz"
if not path_npz.exists():
print(f"File {path_npz} not found, skipping.")
row = None
break
data = np.load(path_npz)
natom = len(data["atnums"])
if itraj == 0:
time.append(data["time"])
row.append(data[field])
if row is None:
continue
if itraj == 0:
time = np.concatenate(time)
row = np.concatenate(row)
sequences.append(row)
sequences = np.array(sequences)
percents = np.array([95, 80, 50, 20, 5])
if field == "temperature":
temp_d = 1100
ndof = 3 * natom - 3
expected = chi2.ppf(percents / 100, ndof) * temp_d / ndof
ymin = chi2.ppf(0.01, ndof) * temp_d / ndof
ymax = chi2.ppf(0.99, ndof) * temp_d / ndof
else:
expected = None
time_unit = 1e-12
num = f"percentiles_{field}_{ensemble}"
plt.close(num)
_, ax = plt.subplots(num=num)
plot_instantaneous_percentiles(
ax,
time / time_unit,
sequences / unit,
percents,
None if expected is None else expected / unit,
ymin,
ymax,
)
ax.set_title(f"{field.title()} percentiles during the {ensemble.upper()} run")
ax.set_xlabel("Time [ps]")
ax.set_ylabel(f"{field.title()} [{unitstr}]")
The following cell plots the temperature percentiles for the NVT and NPT equilibration runs.
plot_openmm_percentiles("nvt", "temperature", "K")
plot_openmm_percentiles("npt", "temperature", "K")
The percentiles look good for the equilibration runs: they quickly reach their theoretical values (black dotted lines) and then fluctuate around them. The following cell plots the temperature percentiles for the initial NVE production runs.
plot_openmm_percentiles("nve", "temperature", "K")
This is clearly not the correct temperature distribution! There is a known problem with restart files in the NVE ensemble in OpenMM. Due to a bug, it tends to lower the temperature of the system. More details on this issue can be found here: https://github.com/openmm/openmm/issues/4948
The following cell plots the percentiles of the volume (over the 100 trajectories), which we can only use to validate that the volume distribution converges. However, we cannot trivially compare these percentiles to an expected distribution.
plot_openmm_percentiles("npt", "volume", "nm$^3$", unit=1e-27, expected="last")
This plot is not completely satisfactory either, as it suggests that the volume fluctuations of the 100 runs exhibit synchronized fluctuations, while they should be independent. (They use different random seeds for their MC Barostat.)
Reusable Code for the Analysis of the Production Runs¶
The analyze function applies the same analysis with STACIE
to different amounts of data (initial and extended production runs).
After the analysis, it generates screen output and figures
as discussed in the minimal example.
BOLTZMANN_CONSTANT = sc.value("Boltzmann constant") # J/K
def analyze(npart: int = 1, ntraj: int = 100, neff_max: int | None = 1000) -> float:
"""Analyze MD trajectories to compute the ionic conductivity.
Parameters
----------
npart
The number of parts in the simulation to load.
The default value of 1 corresponds to only loading the initial production runs.
ntraj
The number of trajectories to load.
neff_max
The maximum number of effective points in the fit, see :func:`stacie.estimate_acint`.
Returns
-------
acint
The estimated ionic conductivity, mainly used for regression testing.
"""
# Get the time step from the first NPZ file.
time = np.load(DATA_ROOT / "sim0000_part00_nve_traj.npz")["time"]
timestep = time[1] - time[0]
def iter_sequences():
"""A generator that only loads one MD trajectory at a time in memory."""
for itraj in range(ntraj):
paths_npz = [
DATA_ROOT / f"sim{itraj:04d}_part{ipart:02d}_nve_traj.npz"
for ipart in range(npart)
]
if not all(path_npz.exists() for path_npz in paths_npz):
print(f"Some of {paths_npz} not found, skipping.")
continue
dipole = []
for path_npz in paths_npz:
data = np.load(path_npz)
dipole.append(data["dipole"])
dipole = np.concatenate(dipole, axis=1)
data = np.load(paths_npz[0])
prefactor = 1.0 / (
data["volume"][0] * data["temperature"].mean() * BOLTZMANN_CONSTANT
)
# The finite difference is equivalent to a block-averaged charge current.
current = np.diff(dipole, axis=1) / timestep
yield prefactor, current
# Configure units for output
uc = UnitConfig(
acint_symbol=r"\sigma",
acint_unit_str=r"S/m",
acint_fmt=".1f",
time_unit=1e-15,
time_unit_str="fs",
time_fmt=".3f",
freq_unit=1e12,
freq_unit_str="THz",
)
# Perform the analysis with STACIE
spectrum = compute_spectrum(
iter_sequences(),
timestep=timestep,
prefactors=None,
include_zero_freq=False,
)
result = estimate_acint(
spectrum, DebyeModel(), neff_max=neff_max, verbose=True, uc=uc
)
# Plot some basic analysis figures.
prefix = "conductivity"
plt.close(f"{prefix}_fitted")
_, ax = plt.subplots(num=f"{prefix}_fitted")
plot_fitted_spectrum(ax, uc, result)
plt.close(f"{prefix}_extras")
_, axs = plt.subplots(2, 2, num=f"{prefix}_extras")
plot_extras(axs, uc, result)
# Return the ionic conductivity.
return result.acint
Analysis of the Initial Production Simulation¶
The following cell computes the ionic conductivity of the molten salt at 1100 K, from the initial production runs (8000 steps each).
conductivity_1 = analyze()
CUTOFF FREQUENCY SCAN cv2l(125%)
neff criterion fcut [THz]
--------- ---------- ----------
15.0 inf 3.78e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-1.41687554e-04, 4.32857722e-01, 2.56728397e+00]))
16.0 inf 4.02e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-1.81902689e-04, 4.35714404e-01, 2.56446750e+00]))
17.1 inf 4.28e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-2.50166855e-04, 4.37129989e-01, 2.56312018e+00]))
18.2 inf 4.56e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-2.12201358e-04, 4.39264241e-01, 2.56094796e+00]))
19.4 inf 4.85e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-2.11285519e-04, 4.40372341e-01, 2.55983894e+00]))
20.7 inf 5.17e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-1.94887911e-04, 4.42707229e-01, 2.55748766e+00]))
22.0 inf 5.50e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-1.57378947e-04, 4.46166592e-01, 2.55399079e+00]))
23.5 inf 5.86e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
25.0 inf 6.23e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
26.7 inf 6.64e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
28.4 inf 7.06e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
30.3 inf 7.52e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
32.3 inf 8.00e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
34.4 inf 8.52e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
36.7 -38.3 9.07e-01
39.1 -60.3 9.65e-01
41.6 -71.2 1.03e+00
44.3 -77.8 1.09e+00
47.2 -82.5 1.16e+00
50.3 -85.9 1.24e+00
53.6 -88.4 1.32e+00
57.1 -90.4 1.40e+00
60.8 -92.0 1.50e+00
64.7 -93.3 1.59e+00
68.9 -94.3 1.69e+00
73.4 -95.2 1.80e+00
78.2 -96.0 1.92e+00
83.3 -96.6 2.04e+00
88.7 -96.5 2.18e+00
94.4 -94.6 2.32e+00
100.5 -89.0 2.47e+00
107.1 -76.4 2.62e+00
114.0 -52.2 2.79e+00
121.4 -9.3 2.97e+00
129.2 61.9 3.17e+00
Cutoff criterion exceeds incumbent + margin: -96.6 + 100.0.
INPUT TIME SERIES
Time step: 50.000 fs
Simulation time: 39950.000 fs
Maximum degrees of freedom: 600.0
MAIN RESULTS
Autocorrelation integral: 342.4 ± 4.9 S/m
Integrated correlation time: 24.786 ± 0.356 fs
SANITY CHECKS (weighted averages over cutoff grid)
Effective number of points: 82.8 (ideally > 60)
Regression cost Z-score: -1.5 (ideally < 2)
Cutoff criterion Z-score: 1.0 (ideally < 2)
MODEL debye() | CUTOFF CRITERION cv2l(125%)
Number of parameters: 3
Average cutoff frequency: 2.03e+00 THz
Exponential correlation time: 94.115 ± 10.971 fs
RECOMMENDED SIMULATION SETTINGS (EXPONENTIAL CORR. TIME)
Block time: < 29.567 ± 3.447 fs
Simulation time: > 5913.404 ± 689.357 fs
The Debye model describes the spectrum well up to a cutoff frequency of about 2 THz. The fit uses 83 effective points, which is sufficient for a three-parameter model. The relative error on the conductivity is 1.4%, which is of the same order of magnitude as the coarse estimate of 0.75%.
The Debye model also provides a relaxation time of about 94 fs. Its relative error (12%) is much larger than that of the conductivity (1.4%), so the curvature of the dip is only moderately constrained by the initial runs. The cutoff weight in the extra plots is concentrated around 2 THz and drops to zero well before the lowest and highest cutoff frequencies are reached, so the cutoff scan is not limited by the grid. A finer frequency grid is still useful to constrain the model more tightly. Hence, we extended the production runs by 8000 additional steps, of which the results are discussed in the following subsection.
Analysis of the Extended Production Simulation (8000 + 8000 steps)¶
We simply call the same analyze() function, but now with npart=2,
which loads the initial production runs and their first extension.
conductivity_2 = analyze(npart=2)
CUTOFF FREQUENCY SCAN cv2l(125%)
neff criterion fcut [THz]
--------- ---------- ----------
15.0 inf 1.89e-01 (No correlation time estimate available.)
16.0 inf 2.01e-01 (No correlation time estimate available.)
17.1 inf 2.14e-01 (No correlation time estimate available.)
18.2 inf 2.28e-01 (No correlation time estimate available.)
19.4 inf 2.43e-01 (No correlation time estimate available.)
20.7 inf 2.58e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
22.0 inf 2.75e-01 (No correlation time estimate available.)
23.5 inf 2.93e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-3.73792913e-06, 4.33897008e-01, 2.56610673e+00]))
25.0 inf 3.11e-01 (No correlation time estimate available.)
26.7 inf 3.32e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-3.28037401e-05, 4.35641126e-01, 2.56439168e+00]))
28.4 inf 3.53e-01 (No correlation time estimate available.)
30.3 inf 3.76e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-1.08351006e-04, 4.39075697e-01, 2.56103265e+00]))
32.3 inf 4.00e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-2.07045291e-04, 4.40510955e-01, 2.55969609e+00]))
34.4 inf 4.26e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-3.55238609e-04, 4.42706286e-01, 2.55764895e+00]))
36.7 inf 4.53e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-4.63517892e-04, 4.44506824e-01, 2.55595669e+00]))
39.1 inf 4.82e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-5.28159961e-04, 4.46790058e-01, 2.55373810e+00]))
41.6 inf 5.14e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-5.21024915e-04, 4.48943271e-01, 2.55157775e+00]))
44.3 inf 5.47e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-3.58399964e-04, 4.50864573e-01, 2.54949383e+00]))
47.2 inf 5.82e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-1.59800057e-06, 4.53043327e-01, 2.54695827e+00]))
50.3 inf 6.19e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
53.6 inf 6.59e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
57.1 inf 7.02e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
60.8 inf 7.47e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
64.7 inf 7.95e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
68.9 inf 8.47e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
73.4 -27.3 9.01e-01
78.2 -51.9 9.59e-01
83.3 -65.0 1.02e+00
88.7 -73.4 1.09e+00
94.4 -79.2 1.16e+00
100.5 -83.5 1.23e+00
107.1 -86.7 1.31e+00
114.0 -89.0 1.40e+00
121.4 -90.5 1.49e+00
129.2 -91.6 1.58e+00
137.6 -92.5 1.68e+00
146.5 -93.4 1.79e+00
156.0 -94.2 1.91e+00
166.1 -94.2 2.03e+00
176.8 -92.0 2.16e+00
188.3 -85.5 2.30e+00
200.4 -71.5 2.45e+00
213.4 -44.4 2.61e+00
227.2 3.8 2.78e+00
241.9 85.7 2.96e+00
Cutoff criterion exceeds incumbent + margin: -94.2 + 100.0.
INPUT TIME SERIES
Time step: 50.000 fs
Simulation time: 79950.000 fs
Maximum degrees of freedom: 600.0
MAIN RESULTS
Autocorrelation integral: 345.8 ± 3.5 S/m
Integrated correlation time: 25.056 ± 0.250 fs
SANITY CHECKS (weighted averages over cutoff grid)
Effective number of points: 156.5 (ideally > 60)
Regression cost Z-score: 0.5 (ideally < 2)
Cutoff criterion Z-score: 2.7 (ideally < 2)
MODEL debye() | CUTOFF CRITERION cv2l(125%)
Number of parameters: 3
Average cutoff frequency: 1.91e+00 THz
Exponential correlation time: 88.607 ± 8.504 fs
RECOMMENDED SIMULATION SETTINGS (EXPONENTIAL CORR. TIME)
Block time: < 27.837 ± 2.671 fs
Simulation time: > 5567.323 ± 534.299 fs
With twice as much data, the number of effective points in the fit also roughly doubles (157), and the uncertainty on the conductivity decreases from 4.9 to 3.5 S/m. The result is consistent with the one from the initial runs. The cutoff criterion Z-score is slightly above the ideal limit of 2, but this does not affect the conclusions, because the conductivity hardly changes over the cutoff scan (extra plot at the top right). To obtain more robust results, we extended the simulations once more. We added 184000 more steps, resulting in a simulation time of 1 ns for each of the 100 trajectories.
Analysis of the Extended Production Simulation (8000 + 8000 + 184000 steps)¶
We simply call the same analyze() function, but now with npart=3,
which loads the initial production runs and their first and second extensions.
conductivity_3 = analyze(npart=3, neff_max=4000)
CUTOFF FREQUENCY SCAN cv2l(125%)
neff criterion fcut [THz]
--------- ---------- ----------
15.0 inf 1.51e-02 (No correlation time estimate available.)
16.0 inf 1.61e-02 (No correlation time estimate available.)
17.1 inf 1.71e-02 (No correlation time estimate available.)
18.2 inf 1.82e-02 (No correlation time estimate available.)
19.4 inf 1.94e-02 (No correlation time estimate available.)
20.7 inf 2.06e-02 (No correlation time estimate available.)
22.0 inf 2.20e-02 (No correlation time estimate available.)
23.5 inf 2.34e-02 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-2.74279924e-05, 4.30590654e-01, 2.56943677e+00]))
25.0 inf 2.49e-02 (No correlation time estimate available.)
26.7 inf 2.65e-02 (No correlation time estimate available.)
28.4 inf 2.82e-02 (No correlation time estimate available.)
30.3 inf 3.00e-02 (No correlation time estimate available.)
32.3 inf 3.20e-02 (No correlation time estimate available.)
34.4 inf 3.40e-02 (No correlation time estimate available.)
36.7 inf 3.62e-02 (No correlation time estimate available.)
39.1 inf 3.86e-02 (No correlation time estimate available.)
41.6 inf 4.11e-02 (No correlation time estimate available.)
44.3 inf 4.37e-02 (No correlation time estimate available.)
47.2 inf 4.65e-02 (No correlation time estimate available.)
50.3 inf 4.95e-02 (No correlation time estimate available.)
53.6 inf 5.27e-02 (No correlation time estimate available.)
57.1 inf 5.61e-02 (No correlation time estimate available.)
60.8 inf 5.97e-02 (No correlation time estimate available.)
64.7 inf 6.36e-02 (No correlation time estimate available.)
68.9 inf 6.77e-02 (No correlation time estimate available.)
73.4 inf 7.21e-02 (No correlation time estimate available.)
78.2 inf 7.67e-02 (No correlation time estimate available.)
83.3 inf 8.17e-02 (No correlation time estimate available.)
88.7 inf 8.69e-02 (No correlation time estimate available.)
94.4 inf 9.25e-02 (No correlation time estimate available.)
100.5 inf 9.85e-02 (No correlation time estimate available.)
107.1 inf 1.05e-01 (No correlation time estimate available.)
114.0 inf 1.12e-01 (No correlation time estimate available.)
121.4 inf 1.19e-01 (No correlation time estimate available.)
129.2 inf 1.26e-01 (No correlation time estimate available.)
137.6 inf 1.35e-01 (No correlation time estimate available.)
146.5 inf 1.43e-01 (No correlation time estimate available.)
156.0 inf 1.53e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-1.54359370e-06, 4.41366618e-01, 2.55863493e+00]))
166.1 inf 1.62e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-1.49142797e-06, 4.41582231e-01, 2.55841926e+00]))
176.8 inf 1.73e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-2.91330612e-06, 4.41639799e-01, 2.55836311e+00]))
188.3 inf 1.84e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-8.02534064e-06, 4.41998992e-01, 2.55800903e+00]))
200.4 inf 1.96e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-1.49369521e-05, 4.42334093e-01, 2.55768084e+00]))
213.4 inf 2.09e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-2.27159654e-05, 4.42747757e-01, 2.55727496e+00]))
227.2 inf 2.22e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-2.87488075e-05, 4.43190437e-01, 2.55683831e+00]))
241.9 inf 2.36e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-4.04863273e-05, 4.43778776e-01, 2.55626171e+00]))
257.5 inf 2.52e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-4.33777791e-05, 4.44331986e-01, 2.55571139e+00]))
274.2 inf 2.68e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-4.06989167e-05, 4.44918800e-01, 2.55512190e+00]))
291.9 inf 2.85e-01 (opt: Hessian matrix has non-positive eigenvalues: evals=array([-2.40531661e-05, 4.45570259e-01, 2.55445379e+00]))
310.7 inf 3.03e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
330.8 inf 3.23e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
352.2 inf 3.44e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
374.9 inf 3.66e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
399.1 inf 3.90e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
424.9 inf 4.15e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
452.3 inf 4.41e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
481.5 inf 4.70e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
512.6 inf 5.00e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
545.7 inf 5.32e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
580.9 inf 5.67e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
618.4 inf 6.03e-01 (rel.err. tau_exp > 1.0e+02 x rel.err. ac integral)
658.4 -6.9 6.42e-01
700.9 -29.1 6.84e-01
746.1 -45.6 7.28e-01
794.2 -57.4 7.75e-01
845.5 -66.2 8.25e-01
900.1 -72.5 8.78e-01
958.1 -77.1 9.35e-01
1020.0 -80.3 9.95e-01
1085.8 -82.8 1.06e+00
1155.8 -84.8 1.13e+00
1230.4 -86.7 1.20e+00
1309.8 -88.4 1.28e+00
1394.3 -90.1 1.36e+00
1484.3 -91.7 1.45e+00
1580.0 -93.3 1.54e+00
1682.0 -94.9 1.64e+00
1790.5 -96.1 1.75e+00
1906.0 -96.0 1.86e+00
2028.9 -92.1 1.98e+00
2159.8 -78.5 2.11e+00
2299.1 -42.5 2.24e+00
2447.5 41.0 2.39e+00
Cutoff criterion exceeds incumbent + margin: -96.1 + 100.0.
INPUT TIME SERIES
Time step: 50.000 fs
Simulation time: 999950.000 fs
Maximum degrees of freedom: 600.0
MAIN RESULTS
Autocorrelation integral: 350.6 ± 1.0 S/m
Integrated correlation time: 25.286 ± 0.074 fs
SANITY CHECKS (weighted averages over cutoff grid)
Effective number of points: 1815.4 (ideally > 60)
Regression cost Z-score: -0.3 (ideally < 2)
Cutoff criterion Z-score: 3.0 (ideally < 2)
MODEL debye() | CUTOFF CRITERION cv2l(125%)
Number of parameters: 3
Average cutoff frequency: 1.77e+00 THz
Exponential correlation time: 92.272 ± 2.979 fs
RECOMMENDED SIMULATION SETTINGS (EXPONENTIAL CORR. TIME)
Block time: < 28.988 ± 0.936 fs
Simulation time: > 5797.593 ± 187.160 fs
The full production runs give a very detailed spectrum.
With the default upper limit of 1000 effective points (neff_max),
the cutoff scan is stopped by this limit rather than by the data.
Hence, the call above raises neff_max to 4000,
which allows the model to be fitted to higher frequencies.
The scan now ends by itself, at about 2400 effective points,
because the cutoff criterion rises well beyond its minimum.
The average cutoff frequency is about 1.8 THz.
The relative error estimate on the conductivity is 0.29%
(with the default limit, it would be 0.33%).
The relaxation time of about 92 fs agrees with the values from the shorter runs,
and its relative error improves from 9% (default limit) to 3%.
This is the main benefit of the higher limit:
the conductivity is well determined because it is the value of the spectrum at zero frequency,
while the relaxation time depends on the curvature of the dip,
which is better constrained by the additional frequencies.
The price is that the cutoff criterion Z-score rises to about 3, slightly above the ideal limit of 2, as for the extended runs above. This shows that the Debye model starts to deviate from the spectrum at the higher frequencies, which is expected because it only describes the onset of the resonance peak. The conductivity changes by less than one standard error when the limit is raised, so this deviation does not affect the conclusions.
Density¶
To enable a proper comparison with the experimental and other simulation results, we also need to estimate the density of the system. This is done by averaging the density over the NpT trajectories from the production runs.
def estimate_density(ntraj: int = 100):
densities = []
molar_vols = []
masses = {11: 22.990, 17: 35.45} # g/mol
avogadro = 6.02214076e23 # 1/mol
for itraj in range(ntraj):
path_npz = DATA_ROOT / f"sim{itraj:04d}_part00_npt_traj.npz"
if not path_npz.exists():
print(f"File {path_npz} not found, skipping.")
continue
data = np.load(path_npz)
mass = sum(masses[atnum] for atnum in data["atnums"]) / avogadro
volume = data["volume"] * 10**6 # from m³ to cm³
densities.append(mass / volume)
molar_vols.append(2 * avogadro * volume / len(data["atnums"]) / 2)
density = np.mean(densities)
print(f"Mass density: {density:.3f} ± {np.std(densities):.3f} g/cm³")
print(f"Molar volume: {np.mean(molar_vols):.4f} ± {np.std(molar_vols):.4f} cm³/mol")
return density
density = estimate_density()
Mass density: 1.454 ± 0.014 g/cm³
Molar volume: 20.1042 ± 0.1977 cm³/mol
Comparison to Literature Results¶
Transport properties for this system are challenging to compute accurately. Consequently, simulation results from the literature may exhibit some variation. While the results should be broadly comparable to some extent, deviations may arise due to the differences in post-processing techniques, and the absence of reported error bars in some studies. Furthermore, in [WSLY14] smaller simulation cells were used (512 ions instead of 1728), which may also contribute to discrepancies.
The table below lists the STACIE results discussed above and some results from the literature.
Ensemble |
Simulated time [ns] |
Density [g/cm3] |
Conductivity [S/m] |
Reference |
|---|---|---|---|---|
NpT+NVE |
4 |
1.454 ± 0.014 |
342.4 ± 4.9 |
init debye |
NpT+NVE |
8 |
1.454 ± 0.014 |
345.8 ± 3.5 |
ext1 debye |
NpT+NVE |
100 |
1.454 ± 0.014 |
350.7 ± 1.0 |
ext2 debye |
NpT+NVT |
6 |
1.456 |
348 ± 7 |
[WDZ+20] |
NpT+NVT |
> 5 |
1.444 |
≈ 310 |
[WSLY14] |
Experiment |
N.A. |
1.542 ± 0.006 |
366 ± 3 |
The comparison shows that the results obtained with STACIE align reasonably well with the literature. In terms of statistical efficiency, STACIE achieves comparable or smaller error bars for about the same simulation time. The deviation from experiment is attributed to the approximations in the NaCl potential. [WDZ+20]
Finally, this example shows how the quality of the result improves with the amount of data. With a model that suits the spectrum, the effective number of points in the fit grows with the simulation time. As a result, the uncertainty on the conductivity decreases steadily: from 4.9 S/m (4 ns) to 3.5 S/m (8 ns) and 1.0 S/m (100 ns). This follows the inverse square root scaling with simulation time, which is expected when the number of effective points grows in proportion to the simulation time.
Technical Details of the Analysis of the Literature Data¶
References for the experimental data:
The following cell converts a molar ionic conductivity from the literature back to a conductivity.
def convert_molar_conductivity():
"""Convert a specific conductance to a conductivity."""
# Parameters taken from Wang 2020 (https://doi.org/10.1063/5.0023225)
# and immediately converted to SI units
molar_conductivity = 140 * 1e-4 # S m²/mol
molar_conductivity_std = 3 * 1e-4 # S m²/mol
density = 1.456 * 1e3 # kg/m³
molar_mass = (22.990 + 35.45) * 1e-3 # kg/mol
molar_volume = molar_mass / density # m³/mol
conductivity = molar_conductivity / molar_volume
conductivity_std = molar_conductivity_std / molar_volume
print("Conductivity [S/m]", conductivity)
print("Conductivity std [S/m]", conductivity_std)
convert_molar_conductivity()
Conductivity [S/m] 348.8021902806297
Conductivity std [S/m] 7.474332648870638
Validation of the Production Runs¶
To further establish that our NVE runs together represent the NpT ensemble, the following two cells perform additional validation checks.
A plot of the conserved quantity of the separate NVE runs, to detect any drift.
The distribution of the instantaneous temperature, which should match the desired NpT distribution. For each individual NVE run and for the combined NVE runs, cumulative distributions are plotted. The function also plots the expected cumulative distribution of the NpT ensemble.
def plot_total_energy(npart: int = 3, ntraj: int = 100):
time = None
energies = []
for itraj in range(ntraj):
if itraj == 0:
time = []
energies_traj = []
for ipart in range(npart):
path_npz = DATA_ROOT / f"sim{itraj:04d}_part{ipart:02d}_nve_traj.npz"
if not path_npz.exists():
print(f"File {path_npz} not found, skipping.")
continue
data = np.load(path_npz)
if itraj == 0:
time.append(data["time"])
energies_traj.append(data["total_energy"])
if itraj == 0:
time = np.concatenate(time)
energies.append(np.concatenate(energies_traj))
num = "total_energy"
plt.close(num)
_, ax = plt.subplots(num=num)
for energies_traj in energies:
plt.plot(time, energies_traj)
plt.title("Total energy of the NVE production runs")
plt.xlabel("Time [ps]")
plt.ylabel("Total energy [kJ/mol]")
plot_total_energy()
There is no noticeable drift in the total energy of the NVE runs. Apart from the usual (and acceptable) numerical noise, the total energy is conserved perfectly.
def plot_temperature_production(npart: int = 3, ntraj: int = 100):
"""Plot cumulative distributions of the instantaneous temperature."""
# Load the temperature data from the NVE production runs.
natom = None
temps = []
for itraj in range(ntraj):
temps.append([])
for ipart in range(npart):
path_npz = DATA_ROOT / f"sim{itraj:04d}_part{ipart:02d}_nve_traj.npz"
if not path_npz.exists():
print(f"File {path_npz} not found, skipping.")
continue
data = np.load(path_npz)
natom = len(data["atnums"])
temps[-1].append(data["temperature"])
temps = np.array([np.concatenate(t) for t in temps])
# Plot the instantaneous and desired temperature distribution.
plt.close("tempprod")
_, ax = plt.subplots(num="tempprod")
ndof = 3 * natom - 3
temp_d = 1100
plot_cumulative_temperature_histogram(ax, temps, temp_d, ndof, "K")
plot_temperature_production()
Alas, as mentioned above, there is still a small mismatch between the obtained and expected NVT temperature distributions. This notebook will be updated after OpenMM issue #4948 has been resolved.
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(conductivity_1 - 342) > 5:
raise ValueError(f"wrong conductivity (production): {conductivity_1:.0f}")
if abs(conductivity_2 - 346) > 4:
raise ValueError(f"wrong conductivity (production): {conductivity_2:.0f}")
if abs(conductivity_3 - 351) > 2:
raise ValueError(f"wrong conductivity (production): {conductivity_3:.0f}")
if abs(density - 1.449) > 0.02:
raise ValueError(f"wrong density (production): {density:.3f}")