In [1]:
#### Uncomment the line below to install thermocepstrum
#! pip install thermocepstrum
In [2]:
import numpy as np
import thermocepstrum as tc

from matplotlib import pyplot as plt
plt.style.use('good_looking.mplstyle')

from scipy.optimize import curve_fit
In [3]:
! mkdir -p Images
In [4]:
## Some constants

# Hartree atomic time unit in seconds
TAUPS = 2.4189e-5
# 1 bohr in Angstrom
BOHR = 0.529177 
# elementary charge in Coulomb
ELCH = 1.60217662e-19
# Boltzmann constant in J/K
k_BOLTZ = 1.38064852e-23 
In [5]:
# define a linear function f(x)=m*x+q
linear = lambda x, m, q: m*x + q

$\mathrm{H_3^+}$¶

Displaced dipoles¶

In [6]:
# Load data
mu_x, mu_y, mu_z = np.loadtxt('Data/displaced_dipole.H3+.dat').T

# Plot
from matplotlib.ticker import MaxNLocator

fig, ax = plt.subplots()

x = np.arange(len(mu_x))

ax.yaxis.set_major_locator(MaxNLocator(integer=True))
ax.set_xticks([x[0], x[-1]/3, 2*x[-1]/3, x[-1]])
ax.set_xticklabels(['I', 'II', 'III', 'I'])

lw = 2.5
ax.plot(x, mu_x, label = r'$x$', lw = lw)
ax.plot(x, mu_y, label = r'$y$', lw = lw)
ax.plot(x, mu_z, label = r'$z$', lw = lw)

delta = x[-1]-x[0]
ax.set_xlim(x[0] - 0.05*delta, x[-1] + 0.05*delta)

M = np.max([mu_x, mu_y, mu_z])
m = np.min([mu_x, mu_y, mu_z])
delta = M - m
ax.set_ylim(m - 0.05*delta, M + 0.05*delta)

ax.set_xlabel('Path')
ax.set_ylabel(r'$\Delta \mu /(Le)$')

ax.legend();

fig.savefig('Images/figure2c.png')
No description has been provided for this image

$\mathrm{K_3Cl}$¶

Displaced dipoles¶

In [7]:
# Load data
x, y, z, mu_x, mu_y, mu_z = np.loadtxt('Data/displaced_dipole.K3Cl.dat').T

# Plot
from matplotlib.ticker import MaxNLocator

fig, ax = plt.subplots()

ax.yaxis.set_major_locator(MaxNLocator(integer=True))
ax.set_xticks([x[0], x[-1]])
ax.set_xticklabels(['I', 'I'])


# x = np.arange(len(mu_x))
lw = 2.5
ax.plot(x, mu_x, label = r'$x$', lw = lw)
ax.plot(y, mu_y, label = r'$y$', lw = lw)
ax.plot(z, mu_z, label = r'$z$', lw = lw)

delta = x[-1]-x[0]
ax.set_xlim(x[0] - 0.05*delta, x[-1] + 0.05*delta)

M = np.max([mu_x, mu_y, mu_z])
m = np.min([mu_x, mu_y, mu_z])
delta = M - m
ax.set_ylim(m - 0.05*delta, M + 0.05*delta)

ax.set_xlabel('Path')
ax.set_ylabel(r'$\Delta \mu /(Le)$')

ax.legend();

fig.savefig('Images/figure3c.png')
No description has been provided for this image

$\mathrm{K_{33}Cl_{31}}$¶

Electrical conductivity computed with the Einstein-Helfand method¶

Load core + Wannier centers data¶

In [8]:
# Load file with conductivity
sigma_file = np.loadtxt('Data/sigma.K-KCl.core_and_wcs.dat').T

# Define timestep
dt = TAUPS*15*20*1000
# Define simulation cell linear size
ALAT = 14.070816430e-10 # meters
# Compute average temperature
TEMP_h = np.mean(np.loadtxt('Data/current.K-KCl.core_and_wcs.dat', skiprows=1).T[1])
# Conversion factor from angstroms to meters 
ang2m = 1e-10
# Conversion factor to obtain the conductivity in S/cm
conv = (ELCH*ang2m)**2 * dt*1e9 * (1/ALAT**3/TEMP_h/k_BOLTZ) / 100


sigma_core_wcs = np.zeros((2, len(sigma_file[0])-2))
sigma_core_wcs[0] = sigma_file[8, 2:]*conv
sigma_core_wcs[1] = np.sqrt(sigma_file[9, 2:])*conv

time = np.arange(1, len(sigma_core_wcs[0])+1)*dt
fitto = int(4.5/6*len(time))

sigma_core_wcs[0] *= time
sigma_core_wcs[1] *= time

(m, q), _ = curve_fit(linear, time[fitto:], sigma_core_wcs[0][fitto:],
                      sigma = (sigma_core_wcs[1][fitto:]))
sig_core_wcs = m
q_core_wcs = q

# Estimate error on sigma as the semi-dispersion of the maximum and minimum fit
(m, q), _ = curve_fit(linear, time[fitto:], sigma_core_wcs[0][fitto:] + sigma_core_wcs[1][fitto:])
dsig1 = m
(m, q), _ = curve_fit(linear, time[fitto:], sigma_core_wcs[0][fitto:] - sigma_core_wcs[1][fitto:])
dsig2 = m
dsig_core_wcs = (dsig1 - dsig2)/2

print('sigma[cores+WCs] = {:4.2f} +\- {:4.2f} S/cm'.format(sig_core_wcs, dsig_core_wcs))
sigma[cores+WCs] = 15.90 +\- 3.16 S/cm

Load ionic charges + lone pair data¶

In [13]:
# Load file with conductivity
sigma_file = np.loadtxt('Data/sigma.K-KCl.nominal_charges_and_lone_pair.dat').T

# Define timestep
dt = TAUPS*15*20*1000
# Define simulation cell linear size
ALAT = 14.070816430e-10 # meters
# Compute average temperature
TEMP_h = np.mean(np.loadtxt('Data/current.K-KCl.nominal_charges_and_lone_pair.dat', skiprows=1).T[1])
# Conversion factor from angstroms to meters 
ang2m = 1e-10
# Conversion factor to obtain the conductivity in S/cm
conv = (ELCH*ang2m)**2 * dt * 1e9 * (1/ALAT**3/TEMP_h/k_BOLTZ) / 100

# Load data
sigma_ions_lp = np.zeros((2, len(sigma_file[0])-2))
sigma_ions_lp[0] = sigma_file[8, 2:]*conv
sigma_ions_lp[1] = np.sqrt(sigma_file[9, 2:])*conv

# 
time = np.arange(1, len(sigma_ions_lp[0])+1)*dt
fitto = int(4.5/6*len(time))

sigma_ions_lp[0] *= time
sigma_ions_lp[1] *= time

(m, q), _ = curve_fit(linear, time[fitto:], sigma_ions_lp[0][fitto:],
                      sigma = (sigma_ions_lp[1][fitto:]))
sig_ions_lp = m
q_ions_lp = q

# Estimate error on sigma as the semi-dispersion of the maximum and minimum fit
(m, q), _ = curve_fit(linear, time[fitto:], sigma_ions_lp[0][fitto:] + sigma_ions_lp[1][fitto:])
dsig1 = m
(m, q), _ = curve_fit(linear, time[fitto:], sigma_ions_lp[0][fitto:] - sigma_ions_lp[1][fitto:])
dsig2 = m
dsig_ions_lp = (dsig1 - dsig2)/2

print('sigma[ionic charges + lone pair] = {:4.2f} +\- {:4.2f} S/cm'.format(sig_ions_lp, dsig_ions_lp))
sigma[ionic charges + lone pair] = 14.74 +\- 2.55 S/cm

Load cross-term data¶

In [14]:
# Load file with conductivity
sigma_file = np.loadtxt('Data/sigma.K-KCl.cross.dat').T

# Define timestep
dt = TAUPS*15*20*1000
# Define simulation cell linear size
ALAT = 14.070816430e-10 # meters
# Compute average temperature
TEMP_h = np.mean(np.loadtxt('Data/current.K-KCl.nominal_charges_and_lone_pair.dat', skiprows=1).T[1])
# Conversion factor from angstroms to meters 
ang2m = 1e-10
# Conversion factor to obtain the conductivity in S/cm
conv = (ELCH*ang2m)**2 * dt*1e9 * (1/ALAT**3/TEMP_h/k_BOLTZ) / 100

sigma_cross = np.zeros((2, len(sigma_file[0])-2))
sigma_cross[0] = sigma_file[12, 2:]*conv
sigma_cross[1] = np.sqrt(sigma_file[13, 2:])*conv

time = np.arange(1, len(sigma_cross[0])+1)*dt
fitto = int(4.5/6*len(time))

sigma_cross[0] *= time
sigma_cross[1] *= time

(m, q), _ = curve_fit(linear, time[fitto:], sigma_cross[0][fitto:],
                      sigma = (sigma_cross[1][fitto:]))
sig_cross = m
q_cross = q

# Estimate error on sigma as the semi-dispersion of the maximum and minimum fit
(m, q), _ = curve_fit(linear, time[fitto:], sigma_cross[0][fitto:] + sigma_cross[1][fitto:])
dsig1 = m
(m, q), _ = curve_fit(linear, time[fitto:], sigma_cross[0][fitto:] - sigma_cross[1][fitto:])
dsig2 = m
dsig_cross = (dsig1 - dsig2)/2

print('sigma[cross-term] = {:4.2f} +\- {:4.2f} S/cm'.format(sig_cross, dsig_cross))
sigma[cross-term] = 0.28 +\- 2.19 S/cm

Load ionic contribution data¶

In [15]:
# Load file with conductivity
sigma_file = np.loadtxt('Data/sigma.K-KCl.ionic.dat').T

# Define timestep
dt = TAUPS*15*20*1000
# Define simulation cell linear size
ALAT = 14.070816430e-10 # meters
# Compute average temperature
TEMP_h = np.mean(np.loadtxt('Data/current.K-KCl.nominal_charges_and_lone_pair.dat', skiprows=1).T[1])
# Conversion factor from angstroms to meters 
ang2m = 1e-10
# Conversion factor to obtain the conductivity in S/cm
conv = (ELCH*ang2m)**2 * dt*1e9 * (1/ALAT**3/TEMP_h/k_BOLTZ) / 100


sigma_ions = np.zeros((2, len(sigma_file[0])-2))
sigma_ions[0] = sigma_file[8, 2:]*conv
sigma_ions[1] = np.sqrt(sigma_file[9, 2:])*conv

time = np.arange(1, len(sigma_ions[0])+1)*dt
fitto = int(4.5/6*len(time))

sigma_ions[0] *= time
sigma_ions[1] *= time

(m, q), _ = curve_fit(linear, time[fitto:], sigma_ions[0][fitto:],
                      sigma = (sigma_ions[1][fitto:]))
sig_ions = m
q_ions = q

# Estimate error on sigma as the semi-dispersion of the maximum and minimum fit
(m, q), _ = curve_fit(linear, time[fitto:], sigma_ions[0][fitto:] + sigma_ions[1][fitto:])
dsig1 = m
(m, q), _ = curve_fit(linear, time[fitto:], sigma_ions[0][fitto:] - sigma_ions[1][fitto:])
dsig2 = m
dsig_ions = (dsig1 - dsig2)/2

print('sigma[ionic contribution] = {:4.2f} +\- {:4.2f} S/cm'.format(sig_ions, dsig_ions))
sigma[ionic contribution] = 3.59 +\- 0.57 S/cm

Plot Figure 8¶

In [16]:
fig, ax = plt.subplots()

#######################################################
s = sigma_core_wcs
m = sig_core_wcs
q = q_core_wcs

time = np.arange(len(s[0]))*dt
pf = ax.plot(time/1000, s[0]/1000, lw=0)
pffill = ax.fill_between(time/1000, (s[0] - s[1])/1000,
                (s[0] + s[1])/1000, color = pf[0].get_color(), alpha = 0.5)
pffit = ax.plot(time/1000, linear(time, m, q)/1000,
        ls = '--',
        lw = 3,
        color = pf[0].get_color())
#######################################################
s = sigma_ions_lp
m = sig_ions_lp
q = q_ions_lp

time = np.arange(len(s[0]))*dt
p = ax.plot(time/1000, s[0]/1000, lw=0)
pfill = ax.fill_between(time/1000, (s[0] - s[1])/1000,
                (s[0] + s[1])/1000, color = p[0].get_color(), alpha = 0.5)
pfit = ax.plot(time/1000, linear(time, m, q)/1000,
        ls = '--',
        lw = 3,
        color = p[0].get_color())
#######################################################
s = sigma_ions
m = sig_ions
q = q_ions

time = np.arange(len(s[0]))*dt
pi = ax.plot(time/1000, s[0]/1000, lw=0)
pifill = ax.fill_between(time/1000, (s[0] - s[1])/1000,
                (s[0] + s[1])/1000, color = pi[0].get_color(), alpha = 0.5)
pifit = ax.plot(time/1000, linear(time, m, q)/1000,
        ls = '--',
        lw = 3,
        color = pi[0].get_color())
#######################################################
s = sigma_cross
m = sig_cross
q = q_cross

time = np.arange(len(s[0]))*dt
pc = ax.plot(time/1000, s[0]/1000, lw=0)
pcfill = ax.fill_between(time/1000, (s[0] - s[1])/1000,
                (s[0] + s[1])/1000, color = pc[0].get_color(), alpha = 0.5)
pcfit = ax.plot(time/1000, linear(time, m, q)/1000,
        ls = '--',
        lw = 3,
        color = pc[0].get_color())
#######################################################
ax.legend([(pffill, pffit[0]), (pfill, pfit[0]), (pifill, pifit[0]), (pcfill, pcfit[0])],
              [r'From Eq. (2)', r'From Eq. (4)', r'Ionic contribution', r'Cross-term'],
              loc = 'upper left')

ax.set_xlabel(r'Time (ps)')
ax.set_ylabel(r'$\langle | \Delta \mu | ^2 \rangle / (6 k_B T L^3)$ (S cm$^{-1}$ ps)')

# ax.tick_params(axis='both', which='major', labelsize = fs)
ax.set_xticks(range(4))
ax.set_yticks(range(0, 80, 20))
ax.set_xlim(0,3.4)
ax.set_ylim(-5,75)

fig.savefig('Images/figure8.png')
No description has been provided for this image

Compute the Electrical conductivity via the Cepstral Method¶

For details on the method please see:

R. Bertossa, F. Grasselli, L. Ercole, and S. Baroni, Phys. Rev. Lett. 122, 255901.

For details on the thermocepstrum code, please see

L. Ercole, R. Bertossa, S. Bisacchi, and S. Baroni, ThermoCepstrum: a code to estimate transport coefficients from the cepstral analysis of a multi-variate current stationary time series, https://github.com/lorisercole/thermocepstrum (2017–2020)

Load cores + Wannier centers data¶

In [10]:
wind = 0.1

dt = TAUPS*15*20*1000
ALAT = 14.070816430e-10 # metres
VOL = 14.070816430**3
fstar = 35

current_file = np.loadtxt('Data/current.K-KCl.core_and_wcs.dat', skiprows=1).T[1:]
TEMP = np.mean(current_file[0])

conv = np.sqrt(TEMP)/10

j_core_wcs = current_file[1:].T * conv

hc_core_wcs = tc.HeatCurrent(j_core_wcs, 'metal', dt, TEMP, VOL, wind)
hcr_core_wcs = hc_core_wcs.resample_current(fstar_THz = fstar, plot = True, PSD_FILTER_W = wind)[0]
hcr_core_wcs.cepstral_analysis(Kmin_corrfactor = 0.8)
Using single component code.
Using single component code.
-----------------------------------------------------
  RESAMPLE TIME SERIES
-----------------------------------------------------
 Original Nyquist freq  f_Ny =      68.90184 THz
 Resampling freq          f* =      34.45092 THz
 Sampling time         TSKIP =             2 steps
                             =        14.513 fs
 Original  n. of frequencies =          4751
 Resampled n. of frequencies =          2376
 PSD      @cutoff  (pre-filter) = 280147.34383
                  (post-filter) = 345764.48312
 log(PSD) @cutoff  (pre-filter) =     12.38206
                  (post-filter) =     12.59417
 min(PSD)          (pre-filter) =     27.23741
 min(PSD)         (post-filter) =  13450.20312
 % of original PSD Power f<f* (pre-filter)  = 96.757350
-----------------------------------------------------

-----------------------------------------------------
  CEPSTRAL ANALYSIS
-----------------------------------------------------
  AIC_Kmin  = 7  (P* = 8, corr_factor = 0.800000)
  L_0*   =          15.611947 +/-   0.049943
  S_0*   =     8718189.758674 +/- 435414.050016
-----------------------------------------------------
  kappa* =          16.167553 +/-   0.807459  W/mK
-----------------------------------------------------

No description has been provided for this image

Load ionic charges + lone pair data¶

In [12]:
wind=0.1

dt = TAUPS*15*20*1000
ALAT = 14.070816430e-10 # metres
VOL = 14.070816430**3
fstar = 35

current_file = np.loadtxt('Data/current.K-KCl.nominal_charges_and_lone_pair.dat', skiprows=1).T[1:]
TEMP = np.mean(current_file[0])

conv = np.sqrt(TEMP)/10

j_ions_lp = current_file[1:].T * conv

hc_ions_lp = tc.HeatCurrent(j_ions_lp, 'metal', dt, TEMP, VOL, wind)
hcr_ions_lp = hc_ions_lp.resample_current(fstar_THz = fstar, plot = True, PSD_FILTER_W = wind)[0]
hcr_ions_lp.cepstral_analysis(Kmin_corrfactor = 0.8)
Using single component code.
Using single component code.
-----------------------------------------------------
  RESAMPLE TIME SERIES
-----------------------------------------------------
 Original Nyquist freq  f_Ny =      68.90184 THz
 Resampling freq          f* =      34.45092 THz
 Sampling time         TSKIP =             2 steps
                             =        14.513 fs
 Original  n. of frequencies =          4751
 Resampled n. of frequencies =          2376
 PSD      @cutoff  (pre-filter) = 293071.92439
                  (post-filter) = 375387.08603
 log(PSD) @cutoff  (pre-filter) =     12.36241
                  (post-filter) =     12.67744
 min(PSD)          (pre-filter) =     26.69397
 min(PSD)         (post-filter) =  24379.45734
 % of original PSD Power f<f* (pre-filter)  = 96.697306
-----------------------------------------------------

-----------------------------------------------------
  CEPSTRAL ANALYSIS
-----------------------------------------------------
  AIC_Kmin  = 7  (P* = 8, corr_factor = 0.800000)
  L_0*   =          15.598159 +/-   0.049943
  S_0*   =     8598806.522775 +/- 429451.672540
-----------------------------------------------------
  kappa* =          15.946161 +/-   0.796402  W/mK
-----------------------------------------------------

No description has been provided for this image

Plot Figure 9¶

In [17]:
fig, ax = plt.subplots()

alpha = 0.3
fmax = fstar
ymax=25

f = hcr_core_wcs
conv = f.kappa_scale/2
psd = ax.plot(f.freqs_THz, conv * f.psd,
              alpha = alpha, zorder = 0, lw=0)
fpsd = ax.plot(f.freqs_THz, conv * f.fpsd,
               color = psd[0].get_color(), zorder = 1)
ceps = ax.plot(f.freqs_THz, conv * f.dct.psd, lw = 2,
               label="From Eq. (2)",
               color = psd[0].get_color(), zorder = 2)
###########################################################################################
f = hcr_ions_lp
conv = f.kappa_scale/2
psd = ax.plot(f.freqs_THz, conv * f.psd,
#               label = "Bare PSD",
              alpha = alpha, zorder = 0, lw=0)
fpsd = ax.plot(f.freqs_THz, conv * f.fpsd,
#                label = "Window filter",
               color = psd[0].get_color(), zorder = 1)
ceps = ax.plot(f.freqs_THz, conv * f.dct.psd, lw = 2,
               label= "From Eq. (4)", #"Cepstral estimate",
               color = psd[0].get_color(), zorder = 2)


ax.set_xlim(0, 8)
ax.set_ylim(0, ymax)

ax.set_xlabel(r'Frequency (THz)')
ax.set_ylabel(r'Conductivity  (S/cm)')

ax.legend()

fig.savefig('Images/figure9.png')
No description has been provided for this image

Energy band-gap¶

In [57]:
fig, ax = plt.subplots()

TEMP = np.loadtxt('Data/current.K-KCl.core_and_wcs.dat', skiprows = 1).T[1].mean()

gap1 = np.loadtxt('Data/gap_HOMOLUMO.K33Cl31.dat')
gap2 = np.loadtxt('Data/gap_HOMO-1LUMO.K33Cl31.dat')

dt = TAUPS*15*20 
time = np.arange(len(gap1))*dt

lw = 0.5

gap1p = ax.plot(time, gap1, label = 'HOMO/LUMO',
                zorder = 0,
                lw = lw);
gap2p = ax.plot(time, gap1+gap2, label = 'HOMO-1/LUMO',
                zorder = 1,
                color='orange',
                lw = lw);

ax.set_xlabel('Time (ps)');
ax.set_ylabel(r'Energy (eV)');

kT = 8.617333262145e-5 * TEMP
ax.axhline(kT,
           lw = 3, color = 'red')

stoic_gap = 2.59 # eV
q1 = 2.01 #2.0619
q2 = 3.00 #2.9374

gapSp = ax.axhline(stoic_gap,
                   lw = 3, color = 'darkgreen', zorder=2)


ax.set_ylim(0,4)
ax.tick_params(axis='both', which='major')

fig.savefig('Images/figure4.png')
No description has been provided for this image