import numpy as np
from matplotlib import pyplot as plt
plt.rcParams['figure.dpi'] = 200
import pandas as pd
from ase.io import read, write
from ase.geometry.analysis import get_rdf
import sportran as st
plt.style.use('AIP.mplstyle')
plt.rcParams['text.latex.preamble'] = r'''\usepackage{amsfonts}
\usepackage{bm}
\usepackage{mathptm}'''
import pyanalisi as pa
from glob import glob
from scipy.optimize import curve_fit
from pathlib import Path
Path('images').mkdir(exist_ok = True)
gofr = {}
for spec in ['Na', 'Cl', 'NaCl']:
gofr[spec] = {}
for func in ['PBE', 'PBE0']:
gofr[spec][func] = {}
for i, f in enumerate(glob(f'./RPDF/{func}/gofr_{spec}_*')):
gofr[spec][func][i] = np.loadtxt(f)
fs = plt.rcParams['figure.figsize']
fig, axes = plt.subplots(ncols = 3, figsize = (2*fs[0], 2*fs[0]/3/1.618), gridspec_kw = {'wspace':0.2})
w = 1
for i, (ax, spec) in enumerate(zip(axes, gofr)):
handles = []
labels = []
for func, c in zip(gofr[spec], ['tab:blue', 'tab:red']):
x = np.mean([gofr[spec][func][i][:, 0] for i in gofr[spec][func]], axis = 0)
y = np.mean([gofr[spec][func][i][:, 1] for i in gofr[spec][func]], axis = 0)
sh = len(list(gofr[spec][func].keys()))
y_std = np.std([gofr[spec][func][i][:, 1] for i in gofr[spec][func]], axis = 0)/np.sqrt(sh)
pl, = ax.plot(x, pd.Series(y).rolling(w).mean(),
color = c,
label = func)
fi = ax.fill_between(x, pd.Series(y-y_std).rolling(w).mean(), pd.Series(y+y_std).rolling(w).mean(),
color = c,
alpha = 0.5)
handles.append((pl, fi))
labels.append(func)
ax.set_xlim(0, 6)
ax.set_ylim(0,)
ax.set_yticks([0, 1, 2])
ax.set_xlabel('Distance $r$ ($\mathrm{\AA}$)')
ax.set_ylabel('$g_{{\mathrm{{{}}}\\textnormal{{-}}b}}(r)$'.format(spec))
if i == 1:
ax.legend(handles, labels)
ax.set_ylim(0, 2.5)
fig.savefig('images/gofr.pdf', dpi = 1200, bbox_inches = 'tight')
def do_GK(log_arr, header, header_list, n_corr_GK, n_blocks_GK, skip = 1):
log_pa = pa.ReadLog(log_arr, header)
block_size = (log_pa.getNtimesteps() - n_corr_GK) // n_blocks_GK
res = []
print('Computing the GK integral for a trajectory with {:d} steps'.format(log_pa.getNtimesteps()))
print('Trajectory divided into {:d} blocks of {:d} steps'.format(n_blocks_GK, block_size))
print('Correlation length = {} steps'.format(n_corr_GK))
gk = pa.GreenKubo(log_pa, # Log object
'', # Name of dump file
skip, # Skip in doing the trajectory average
header_list, # header list
False, # if True, dump a file
n_corr_GK, # Max length of correlation function
4, # Number of threads
False, # if true subtract the average from the correlation function before integrating
0, # calculate the average starting from this time
1, # number of segments (it does impact only cache performance)
False, # do a benchmark
1, # benchmark parameter start
100 # benchmark parameter stop
)
gk.reset(block_size)
for ib in range(n_blocks_GK):
print(ib)
gk.calculate(ib*block_size)
res.append(np.array(gk, copy=True))
return np.array(res)
atoms = {}
atoms['PBE'] = read('trajectories/NaCl_bip.PBE.extxyz', index = '-28000:')
atoms['PBE0'] = read('trajectories/NaCl_bip.PBE0.extxyz', index = '-14000:')
rdf_peaks = {}
spread = {}
for func in atoms:
spread[func] = np.array([a.info['spread'] for a in atoms[func]])
rdf_peaks[func] = []
for a in atoms[func]:
g, r = get_rdf(a, rmax = 6, elements=[0, 11], nbins = 30)
peak = np.argmax(g)
rdf_peaks[func].append(r[peak])
nb = {'PBE0': 13, 'PBE': 25}
ncorr = {'PBE0': 1000, 'PBE': 1000}
GK = {}
for func in atoms:
log = np.array([spread[func], spread[func], spread[func], rdf_peaks[func], rdf_peaks[func], rdf_peaks[func]]).T
log_arr = np.copy((log-log.mean(axis=0))/log.std(axis=0), order = 'C')
header = ['spread[1]', 'spread[2]', 'spread[3]', 'dist[1]', 'dist[2]', 'dist[3]']
header_list = ['spread[1]', 'dist[1]']
GK[func] = do_GK(log_arr,
header,
header_list,
ncorr[func],
nb[func],
skip = 1)
colors = {'PBE': 'tab:blue', 'PBE0': 'tab:red'}
fs = plt.rcParams['figure.figsize']
fig, axes = plt.subplots(nrows = 2, sharex = True, gridspec_kw = {'hspace': 0}, figsize = (fs[0], 1.*fs[1]))
normalize = False
for func, ax in zip(GK, axes):
labels = []
handles = []
time = np.arange(GK[func].shape[1])/1000
y_ = GK[func][:,:,0]
if normalize:
y_ = np.array([y__/y__[0] for y__ in y_])
y = y_.mean(axis=0)
y_std = y_.std(axis=0)/np.sqrt(GK[func].shape[0]*2)
pl, = ax.plot(time, y, color = colors[func], zorder = 1)
fi = ax.fill_between(time, y-y_std, y+y_std, color = colors[func], alpha = 0.3, zorder = 1)
labels.append(r'$\varsigma$-$\varsigma$')
handles.append((pl, fi))
y_ = GK[func][:,:,3]
if normalize:
y_ = np.array([y__/y__[0] for y__ in y_])
y = y_.mean(axis=0)
y_std = y_.std(axis=0)/np.sqrt(GK[func].shape[0]*2)
pl, = ax.plot(time, y, label = 'min dist', ls = '--', color = colors[func], zorder = 2)
fi = ax.fill_between(time, y-y_std, y+y_std, color = colors[func], alpha = 0.3, zorder = 2)
labels.append(r'$r_{\mathrm{peak}}$-$r_{\mathrm{peak}}$')
# labels.append(r'$d_{\mathrm{min}}$-$d_{\mathrm{min}}$')
handles.append((pl, fi))
y_ = np.vstack([GK[func][:,:,1], GK[func][:,:,2]])
if normalize:
y_ = np.array([y__/y__[0] for y__ in y_])
y = y_.mean(axis=0)
y_std = y_.std(axis=0)/np.sqrt(GK[func].shape[0]*2)
pl, = ax.plot(time, y, label = 'mixed', ls = 'dotted', color = colors[func], zorder = 3)
fi = ax.fill_between(time, y-y_std, y+y_std, color = colors[func], alpha = 0.3, zorder = 3)
labels.append(r'$\varsigma$-$r_{\mathrm{peak}}$')
handles.append((pl, fi))
ax.set_xlim(0, 1)
ax.axhline(0, color = 'black', ls = 'solid', lw = 0.5, zorder = 0)
ax.set_ylim(-0.3, 1.3)
ax.set_xticks([0, 0.5, 1])
ax.set_yticks([0, 1])
ax.text(0.15, 0.8, func, transform = ax.transAxes)
ax.legend(handles, labels, ncol=2)
ax.set_xlabel('Time $t$ (ps)')
fig.text(0.065, 0.5, r'Correlation function', rotation = 90, ha = 'center', va = 'center')
fig.savefig('images/correlation_function_spread-peak.pdf', dpi = 1200, bbox_inches = 'tight')
nb = {'PBE0': 6, 'PBE': 10}
ncorr = {'PBE0': 2000, 'PBE':2000}
GK_output = {}
for func in ['PBE', 'PBE0']:
GK_output[func] = {}
for key in ['polarization', 'OS', 'bipolaron', 'topological']:
print(func, key)
log_arr = np.loadtxt(f'fluxes/charge_flux.{func}.{key}.dat')
header = ['c_flux[1]', 'c_flux[2]', 'c_flux[3]']
header_list = ['c_flux[1]']
GK_output[func][key] = do_GK(log_arr,
header,
header_list,
ncorr[func],
nb[func])
## 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
ALAT = 13.0e-10 # meters
TEMP = 1250 # K
ang2m = 1e-10
DT = 1 # fs
conv = (ELCH*ang2m)**2 * DT*1e9 * (1/ALAT**3/TEMP/k_BOLTZ) / 100 * 1e6
fs = plt.rcParams['figure.figsize']
fig, axes = plt.subplots(ncols = 2, figsize = (fs[0], .9*fs[1]), gridspec_kw = {'wspace': 0})
# 'polarization', 'OS', 'bipolaron', 'topological'
lbl = {'bipolaron': r'Bipolaron',
'topological': r'Total',
'OS': r'Ions'}
colors = {'bipolaron': r'tab:green',
'topological': r'black',
'OS': r'tab:orange'}
colors = {
'PBE': {'bipolaron': r'tab:green',
'topological': r'tab:gray',
'OS': r'tab:orange'},
'PBE0': {'bipolaron': r'tab:green',
'topological': r'tab:gray',
'OS': r'tab:orange'}
}
low_cut = 0.8
for ax, func in zip(axes, ['PBE', 'PBE0']):
handles = []
labels = []
for key in ['topological',
'bipolaron',
'OS']:
# Contribution of the different fluxes
ytot = GK_output[func][key][:, :, -1]
y = ytot.mean(axis=0)
ye = ytot.std(axis=0)/np.sqrt(ytot.shape[0])
conv = 6.509739091807338e-08
x = np.arange(y.size)*DT
pl, = ax.plot(x/1000, x*y*conv, zorder = 2, color = colors[func][key])
pl2 = ax.fill_between(x/1000, x*(y-ye)*conv, x*(y+ye)*conv,
color = pl.get_color(),
alpha = 0.5, zorder = 0)
x_fit = x[x/1000>low_cut]/1000
y_fit = x[x/1000>low_cut]*y[x/1000>low_cut]*conv
popt, pcov = curve_fit(lambda t, m, q: m*t+q,
x_fit, y_fit)
m, q = popt
pl3, = ax.plot(x/1000, m*x/1000+q, ls = '--', color = pl.get_color())
handles.append((pl, pl2, pl3))
labels.append(lbl[key])
if func == 'PBE':
low_cut = 1.5
else:
low_cut = 1.5
ax.set_xlim(0,2)
ax.set_ylim(0,30)
ax.set_yticks([0, 10, 20, 30])
ax.set_xlabel('Time $t$ (ps)')
ax.legend(handles, labels,
title = func,
fontsize = 6,
title_fontsize = 6,
ncol = 1, loc = 'upper left')
ax = axes[0]
ax.set_xticks([0, 1])
ax.set_ylabel('Mean square displaced dipole\n' +\
r'$\langle|\bm{\Delta\mu}|^2\rangle/(6k_BTL^3)$ ($\mathrm{S\,cm^{-1}ps}$)')
ax = axes[1]
ax.set_ylabel('')
ax.set_xticks([0, 1, 2])
# ax.set_yticks([])
ax.set_yticklabels([])
fig.savefig('images/MSDD.pdf', dpi = 1200, bbox_inches = 'tight')
local_J = np.load('fluxes/local_flux.npy', allow_pickle = True).item()
bipolaron_J = {'PBE': np.loadtxt('fluxes/charge_flux_bipolaron.PBE.dat'),
'PBE0': np.loadtxt('fluxes/charge_flux_bipolaron.PBE0.dat')}
alat = 13.0 # Angstrom
temperature = 1300 # Kelvin
pstar = None
pstar = 2
fstar = 10
sigma = {}
for func in local_J:
sigma[func] = {}
for cut in local_J[func]:
sigma[func][cut] = {}
J_ = st.ElectricCurrent(local_J[func][cut],
DT_FS = 1,
UNITS = 'metal',
TEMPERATURE = temperature,
VOLUME = alat**3)
J_ = J_.resample(fstar_THz=fstar)
J_.cepstral_analysis(manual_cutoffK=pstar)
# print(func, cut, 'loc', J_.cepf.cutoffK)
sigma[func][cut]['local'] = np.array([J_.kappa, J_.kappa_std])
J_ = st.ElectricCurrent(local_J[func][cut] + bipolaron_J[func],
DT_FS = 1,
UNITS = 'metal',
TEMPERATURE = temperature,
VOLUME = alat**3)
J_ = J_.resample(fstar_THz=fstar)
# pstar = None
J_.cepstral_analysis(manual_cutoffK=pstar)
print(func, cut, 'tot', J_.cepf.cutoffK)
sigma[func][cut]['total'] = np.array([J_.kappa, J_.kappa_std])
J_ = st.ElectricCurrent(bipolaron_J[func],
DT_FS = 1,
UNITS = 'metal',
TEMPERATURE = temperature,
VOLUME = alat**3)
J_ = J_.resample(fstar_THz=fstar)
J_.cepstral_analysis(manual_cutoffK=pstar)
print(func, cut, 'bip', J_.cepf.cutoffK)
sigma[func][cut]['lone pair'] = np.array([J_.kappa, J_.kappa_std])
fs = plt.rcParams['figure.figsize']
fig, axes = plt.subplots(nrows = 2,
figsize = (fs[0], fs[1]*1.1),
gridspec_kw = {'hspace': 0}, sharex = True)
ms = 4.5
for func, c, ax in zip(sigma, ['tab:blue', 'tab:red'], axes):
s0 = np.array([sigma[func][cut]['total'][0] for cut in sigma[func]])/100
s1 = np.array([sigma[func][cut]['lone pair'][0] for cut in sigma[func]])/100
s2 = np.array([sigma[func][cut]['local'][0] for cut in sigma[func]])/100
ds0 = np.array([sigma[func][cut]['total'][1] for cut in sigma[func]])/100
ds1 = np.array([sigma[func][cut]['lone pair'][1] for cut in sigma[func]])/100
ds2 = np.array([sigma[func][cut]['local'][1] for cut in sigma[func]])/100
x = np.array([cut for cut in sigma[func]])
ev = 1
b = 0
ax.errorbar(x[b:-1:ev], s0[b:-1:ev], ds0[b:-1:ev],
lw = 0.3,
elinewidth = 1,
ls = '-',
marker = 'o',
markersize = ms,
color = c,
label = '$\overline{\sigma}$',
capsize = 0,
capthick = 1)
eb = ax.errorbar(x[b:-1:ev], (s1+s2)[b:-1:ev], (ds1+ds2)[b:-1:ev],
lw = 0.3,
elinewidth = 1,
markersize = ms,
markeredgewidth=1,
marker = 'v',
markerfacecolor='white',
color = c,
label = '$\sigma_{b} + \sigma_{\mathrm{loc}}$',
capsize = 0,
capthick = 1)
ax.legend(title = func, ncol = 2, handletextpad = 0.4, loc = 'upper left')
ax.axvline(1, color = 'black', ls = 'dotted', lw = 0.5)
ax.set_xlabel(r'Local flux cutoff $\lambda/\varsigma_{\mathrm{avg}}$')
ax.set_xlim(0,2)
print(func, s0[-1], ds0[-1])
axes[0].set_ylim(9,20)
axes[1].set_ylim(0,2.8)
fig.text(0.045, 0.5, r'Local electrical conductivity (S/cm)', rotation = 90, ha = 'center', va = 'center')
fig.savefig('images/conductivity.pdf', dpi = 1200, bbox_inches = 'tight')