Example 1: Cepstral Analysis of solid amorphous Silica¶
This example shows the basic usage of sportran to compute the thermal conductivity of a classical MD simulation of a-SiO\(_2\).
[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. with LAMMPS format).
For other input formats see Example #3.
[3]:
jfile = st.i_o.TableFile('./data/Silica/Silica.dat', group_vectors=True)
# Solid Silica - BKS potential, melted and quenched
# 216 atoms, T~1000K, dens~2.295g/cm^3
# NVE, dt = 1.0 fs, 100 ps, print_step = 1.0 fs
# Temperature = 983.172635 K, Volume = 3130.431110818 A^3
# LAMMPS metal units
Temp c_flux1[1] c_flux1[2] c_flux1[3]
#####################################
all_ckeys = [('Temp', [0]), ('flux1', array([1, 2, 3]))]
#####################################
Data length = 100001
[4]:
jfile.read_datalines(start_step=0, NSTEPS=0, select_ckeys=['flux1', 'Temp'], max_vector_dim=3)
ckey = [('Temp', [0]), ('flux1', array([1, 2, 3]))]
step = 100000 - 100.00% completed
( 100000 ) steps read.
DONE. Elapsed time: 0.9626550674438477 seconds
[4]:
{'flux1': array([[ -265.30586 , 1520.6107 , 67.461829],
[ -168.68352 , 1377.4459 , 101.82146 ],
[ -93.688306, 1180.375 , 117.20939 ],
...,
[ 1226.9778 , 212.0939 , -1126.4643 ],
[ 1223.8753 , 186.93836 , -881.39541 ],
[ 1232.7723 , 141.30647 , -620.41895 ]]),
'Temp': array([[ 998.48171],
[1003.699 ],
[1003.8906 ],
...,
[ 967.21723],
[ 978.47566],
[ 985.41455]])}
2. Heat Current¶
Define a HeatCurrent from the trajectory, with the correct parameters.
[5]:
DT_FS = 1.0 # time step [fs]
TEMPERATURE = np.mean(jfile.data['Temp']) # temperature [K] (983.173 K)
VOLUME = 3130.431110818 # volume [A^3]
print(DT_FS, TEMPERATURE, VOLUME)
j = st.HeatCurrent(jfile.data['flux1'], UNITS='metal', DT_FS=DT_FS, TEMPERATURE=TEMPERATURE, VOLUME=VOLUME)
1.0 983.1726353043 3130.431110818
Using single component code.
[6]:
# trajectory
f = plt.figure()
ax = plt.plot(j.timeseries()/1000., j.traj);
plt.xlim([0, 1.0])
plt.xlabel(r'$t$ [ps]')
plt.ylabel(r'$J$ [eV A/ps]');
Compute the Power Spectral Density (\(\rm PSD\)) and filter it for visualization.
[7]:
# Periodogram with given filtering window width
ax = j.plot_periodogram(PSD_FILTER_W=0.5, kappa_units=True)
print(j.Nyquist_f_THz)
plt.xlim([0, 50])
ax[0].set_ylim([0, 150]);
ax[1].set_ylim([12, 18]);
500.0
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 = 28.0
jf, ax = j.resample(fstar_THz=FSTAR_THZ, plot=True, freq_units='thz')
plt.xlim([0, 80])
ax[1].set_ylim([12,18]);
Using single component code.
-----------------------------------------------------
RESAMPLE TIME SERIES
-----------------------------------------------------
Original Nyquist freq f_Ny = 500.00000 THz
Resampling freq f* = 27.77778 THz
Sampling time TSKIP = 18 steps
= 18.000 fs
Original n. of frequencies = 50001
Resampled n. of frequencies = 2778
PSD @cutoff (pre-filter&sample) ~ 2666288.02378
(post-filter&sample) ~ 2206043.36451
log(PSD) @cutoff (pre-filter&sample) ~ 14.61451
(post-filter&sample) ~ 14.40954
min(PSD) (pre-filter&sample) = 4.03008
min(PSD) (post-filter&sample) = 60168.84968
% of original PSD Power f<f* (pre-filter&sample) = 77.164 %
-----------------------------------------------------
[9]:
ax = jf.plot_periodogram(PSD_FILTER_W=0.1)
ax[1].set_ylim([12, 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) = 33 (auto, AIC_Kmin = 33, corr_factor = 1.0)
L_0* = 13.114928 +/- 0.097614
S_0* = 717769.108506 +/- 70064.257396
-----------------------------------------------------
kappa* = 2.205099 +/- 0.215248 W/m/K
-----------------------------------------------------
[11]:
# Cepstral Coefficients
print('c_k = ', jf.cepf.logpsdK)
ax = jf.plot_ck()
ax.set_xlim([0, 50])
ax.set_ylim([-0.5, 0.5])
ax.grid();
c_k = [ 1.51720808e+01 -7.09450459e-01 -1.75761630e-01 ... 6.61581529e-04
-1.41143322e-02 5.27000116e-03]
[12]:
# AIC function
f = plt.figure()
plt.plot(jf.cepf.aic, '.-', c=c[0])
plt.xlim([0, 150])
plt.ylim([min(jf.cepf.aic[:150]) - 10, max(jf.cepf.aic[50:150]) + 10]);
print('K of AIC_min = {:d}'.format(jf.cepf.aic_Kmin))
print('AIC_min = {:f}'.format(jf.cepf.aic_min))
K of AIC_min = 33
AIC_min = 2948.593166
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, 150])
ax.set_ylim([12.5, 14.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 = 33
AIC_min = 2948.593166
[14]:
# kappa as a function of cutoff K
ax = jf.plot_kappa_Pstar()
ax.set_xlim([0, 150])
ax.set_ylim([0, 5.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 = 33
AIC_min = 2948.593166
Print the results :)
[15]:
results = jf.cepstral_log
print(results)
-----------------------------------------------------
CEPSTRAL ANALYSIS
-----------------------------------------------------
cutoffK = (P*-1) = 33 (auto, AIC_Kmin = 33, corr_factor = 1.0)
L_0* = 13.114928 +/- 0.097614
S_0* = 717769.108506 +/- 70064.257396
-----------------------------------------------------
kappa* = 2.205099 +/- 0.215248 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., 50.])
ax[1].set_ylim([12,18])
ax[0].legend(['original', 'resampled', 'cepstrum-filtered'])
ax[1].legend(['original', 'resampled', 'cepstrum-filtered']);
[ ]: