Example 2: Cepstral Analysis of liquid NaCl¶
This example shows the basic usage of sportran to compute the thermal conductivity of a classical MD simulation of molten NaCl, which requires the multi-component theory.
For reference: Bertossa et al., PRL 122, 255901 (2019)
[1]:
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
try:
import sportran as st
except ImportError:
from sys import path
path.append('..')
import sportran as st
c = plt.rcParams['axes.prop_cycle'].by_key()['color']
[2]:
%matplotlib notebook
1. Load trajectory¶
Read the heat current from a simple column-formatted file. The desired columns are selected based on their header (e.g. in LAMMPS format).
For other input formats see the corresponding example.
[3]:
jfile = st.i_o.TableFile('./data/NaCl/NaCl.dat', group_vectors=True)
# Molten NaCl - Fumi-Tosi potential
# https://journals.aps.org/prl/supplemental/10.1103/PhysRevLett.122.255901/Bertossa_et_al_Supplemental_Material.pdf
# 1728 atoms, T~1400K
# NVE, dt = 1.0 fs, 100 ps, print_step = 5.0 fs
# Volume = 65013.301261 A^3
# LAMMPS metal units
Temp c_flux[1] c_flux[2] c_flux[3] c_vcm[1][1] c_vcm[1][2] c_vcm[1][3]
#####################################
all_ckeys = [('Temp', [0]), ('flux', array([1, 2, 3])), ('vcm[1]', array([4, 5, 6]))]
#####################################
Data length = 20000
We need to load the energy flux (flux) and one mass flux (vcm[1], that is the velocity of the center of mass of one species).
[4]:
jfile.read_datalines(start_step=0, NSTEPS=0, select_ckeys=['Temp', 'flux', 'vcm[1]'])
ckey = [('Temp', [0]), ('flux', array([1, 2, 3])), ('vcm[1]', array([4, 5, 6]))]
( 20000 ) steps read.
DONE. Elapsed time: 0.2213594913482666 seconds
[4]:
{'Temp': array([[1442.7319],
[1440.806 ],
[1440.6008],
...,
[1402.7391],
[1392.6822],
[1385.4031]]),
'flux': array([[ 250.86549 , 20.619423 , 200.115 ],
[ 196.22265 , 82.667342 , 284.3325 ],
[ 126.39441 , 160.75472 , 340.36769 ],
...,
[ 179.91856 , 18.612706 , -132.65623 ],
[ 204.71193 , -0.46643315, -204.0165 ],
[ 241.23318 , -18.295461 , -252.46475 ]]),
'vcm[1]': array([[-0.15991832, -0.07137043, 0.02068792],
[-0.13755206, -0.07100293, -0.01127988],
[-0.10615044, -0.06238124, -0.04156812],
...,
[-0.0919399 , -0.08477829, 0.06001139],
[-0.13384949, -0.1147453 , 0.08932355],
[-0.18385053, -0.1369343 , 0.1143406 ]])}
2. Heat Current¶
Define a HeatCurrent from the trajectory, with the correct parameters. The difference with respect to the single-component case is only here: we need to load a list of currents instead of a single one.
[5]:
DT_FS = 5.0 # time step [fs]
TEMPERATURE = np.mean(jfile.data['Temp']) # temperature [K]
VOLUME = 40.21**3 # volume [A^3]
print('T = {:f} K'.format(TEMPERATURE))
print('V = {:f} A^3'.format(VOLUME))
j = st.HeatCurrent([jfile.data['flux'], jfile.data['vcm[1]']], UNITS='metal', DT_FS=DT_FS,
TEMPERATURE=TEMPERATURE, VOLUME=VOLUME)
T = 1399.347781 K
V = 65013.301261 A^3
Using multicomponent code.
[6]:
# trajectory
f, ax = plt.subplots(2, sharex=True)
ax[0].plot(j.timeseries()/1000., j.traj);
ax[1].plot(j.timeseries()/1000., j.otherMD[0].traj);
plt.xlim([0, 1.0])
plt.xlabel(r'$t$ [ps]')
ax[0].set_ylabel(r'$J^0$ [eV A/ps]');
ax[1].set_ylabel(r'$J^1$ [A/ps]');
Compute the Reduced periodogram \(\bar{\mathcal{S}}^0_k\) and filter it for visualization. You can notice the difference with respect to the energy-flux periodogram \(\mathcal{S}^0_k\).
[7]:
# Periodogram with given filtering window width
ax = j.plot_periodogram(PSD_FILTER_W=0.4, kappa_units=True, label=r'$\bar{\mathcal{S}}^0_k$')
print(j.Nyquist_f_THz)
# compare with the spectrum of the energy flux
jen = st.HeatCurrent(jfile.data['flux'], UNITS='metal', DT_FS=DT_FS,
TEMPERATURE=TEMPERATURE, VOLUME=VOLUME)
ax = jen.plot_periodogram(axes=ax, PSD_FILTER_W=0.4, kappa_units=True, label=r'$\mathcal{S}^0_k$')
plt.xlim([0, 20])
ax[0].set_ylim([0, 0.8]);
ax[1].set_ylim([7, 18]);
ax[0].legend(); ax[1].legend();
100.0
Using single component code.
3. Resampling¶
If the Nyquist frequency is very high (i.e. the sampling time is very small), such that the log-spectrum goes to low values, you may want resample your time series, in order to reduce the maximum frequency to \(f^*\). Before performing that operation, the time series is automatically filtered to reduce the amount of aliasing introduced. Ideally you do not want to go too low in \(f^*\), because the statistical error may increase. In an intermediate region the results should not change.
To perform resampling you can choose the resampling frequency \(f^*\) (fstar_THz) or the resampling step (TSKIP). If you choose \(f^*\), the code will try to find the closest value allowed. The resulting \(\rm PSD\) can be visualized to ensure that the low-frequency region was not affected by this operations
[8]:
FSTAR_THZ = 14.0
jf, ax = j.resample(fstar_THz=FSTAR_THZ, plot=True, freq_units='thz')
plt.xlim([0, 20])
ax[1].set_ylim([7, 18]);
Using multicomponent code.
-----------------------------------------------------
RESAMPLE TIME SERIES
-----------------------------------------------------
Original Nyquist freq f_Ny = 100.00000 THz
Resampling freq f* = 14.28571 THz
Sampling time TSKIP = 7 steps
= 35.000 fs
Original n. of frequencies = 10001
Resampled n. of frequencies = 1429
PSD @cutoff (pre-filter&sample) ~ 392810.55342
(post-filter&sample) ~ 367183.46153
log(PSD) @cutoff (pre-filter&sample) ~ 12.70308
(post-filter&sample) ~ 12.42357
min(PSD) (pre-filter&sample) = 0.31536
min(PSD) (post-filter&sample) = 22166.11934
% of original PSD Power f<f* (pre-filter&sample) = 96.679 %
-----------------------------------------------------
[9]:
ax = jf.plot_periodogram(PSD_FILTER_W=0.1)
ax[1].set_ylim([7, 18]);
4. Cepstral Analysis¶
Perform Cepstral Analysis. The code will perform the following operations: 1. compute the parameters describing the theoretical distribution of the \(\rm PSD\) 2. compute the Cepstral coefficients by Fourier transforming the \(\log(\rm{PSD})\) 3. apply the Akaike Information Criterion, to choose the number of cepstral coefficients \(P^*\) 4. return the resulting transport coefficient (thermal conductivity) \(\kappa\)
[10]:
jf.cepstral_analysis()
-----------------------------------------------------
CEPSTRAL ANALYSIS
-----------------------------------------------------
cutoffK = (P*-1) = 3 (auto, AIC_Kmin = 3, corr_factor = 1.0)
L_0* = 15.158757 +/- 0.056227
S_0* = 6824108.702608 +/- 383697.095268
-----------------------------------------------------
kappa* = 0.498310 +/- 0.028018 W/m/K
-----------------------------------------------------
[11]:
# Cepstral Coefficients
print('c_k = ', jf.cepf.logpsdK)
ax = jf.plot_ck()
ax.set_xlim([0, 25])
ax.set_ylim([-0.2, 1.0])
ax.grid();
c_k = [ 1.43138519e+01 6.58383813e-01 -6.48049965e-03 ... 4.76081379e-03
2.28241054e-02 3.95167451e-02]
[12]:
# AIC function
f = plt.figure()
plt.plot(jf.cepf.aic, '.-', c=c[0])
plt.xlim([0, 50])
plt.ylim([1400, 1500]);
print('K of AIC_min = {:d}'.format(jf.cepf.aic_Kmin))
print('AIC_min = {:f}'.format(jf.cepf.aic_min))
K of AIC_min = 3
AIC_min = 1410.555757
Plot the thermal conductivity \(\kappa\) as a function of the cutoff \(P^*\)
[13]:
# L_0 as a function of cutoff K
ax = jf.plot_L0_Pstar()
ax.set_xlim([0, 50])
ax.set_ylim([14.75, 16.5]);
print('K of AIC_min = {:d}'.format(jf.cepf.aic_Kmin))
print('AIC_min = {:f}'.format(jf.cepf.aic_min))
K of AIC_min = 3
AIC_min = 1410.555757
[14]:
# kappa as a function of cutoff K
ax = jf.plot_kappa_Pstar()
ax.set_xlim([0, 50])
ax.set_ylim([0, 1.0]);
print('K of AIC_min = {:d}'.format(jf.cepf.aic_Kmin))
print('AIC_min = {:f}'.format(jf.cepf.aic_min))
K of AIC_min = 3
AIC_min = 1410.555757
Print the results :)
[15]:
results = jf.cepstral_log
print(results)
-----------------------------------------------------
CEPSTRAL ANALYSIS
-----------------------------------------------------
cutoffK = (P*-1) = 3 (auto, AIC_Kmin = 3, corr_factor = 1.0)
L_0* = 15.158757 +/- 0.056227
S_0* = 6824108.702608 +/- 383697.095268
-----------------------------------------------------
kappa* = 0.498310 +/- 0.028018 W/m/K
-----------------------------------------------------
You can now visualize the filtered PSD…
[16]:
# filtered log-PSD
ax = j.plot_periodogram(0.5, kappa_units=True)
ax = jf.plot_periodogram(0.5, axes=ax, kappa_units=True)
ax = jf.plot_cepstral_spectrum(axes=ax, kappa_units=True)
ax[0].axvline(x = jf.Nyquist_f_THz, ls='--', c='r')
ax[1].axvline(x = jf.Nyquist_f_THz, ls='--', c='r')
plt.xlim([0, 20])
ax[0].set_ylim([0, 0.8]);
ax[1].set_ylim([7, 18]);
ax[0].legend(['original', 'resampled', 'cepstrum-filtered'])
ax[1].legend(['original', 'resampled', 'cepstrum-filtered']);
[ ]: