In [5]:

plt.plot(time/1000, elc[:,6]*DT*10)
plt.fill_between(time/1000, (elc[:,6] + np.sqrt(elc[:,7]))*DT*10, (elc[:,6] - np.sqrt(elc[:,7]))*DT*10, alpha=0.4, label='Green-Kubo')
plt.plot(time/1000, elc[:,8]*DT*10)
plt.fill_between(time/1000, (elc[:,8] + np.sqrt(elc[:,9]))*DT*10, (elc[:,8] - np.sqrt(elc[:,9]))*DT*10, alpha=0.4, label='Helfand-Einstein')
plt.xlim(0,2.)
plt.ylim(0,1)
plt.xlabel('$\mathcal{T}$ [ps]')
plt.ylabel('time integral [$\\times 10^{2} \;\AA^2$/ps]')
plt.tick_params('both', bottom=True, top=False, left=True, right=False, direction='in')
plt.legend(frameon=False, loc='lower right')
plt.yticks([0,0.5,1.0])
plt.xticks([0, 1, 2])
plt.savefig('GK-HE_comparison.pdf', bbox_inches = 'tight', pad_inches = 0)

Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode

No description has been provided for this image

In [6]:

fig,ax = plt.subplots()
from matplotlib.ticker import FuncFormatter, MultipleLocator

x=np.arange(100)*0.2
ax.plot(x, 2*np.sinc(x/np.pi), label='$\\tilde{\Theta}^{\mathcal{T}}_{GK}(\omega\mathcal{T})$')
ax.plot(x, np.sinc(x/2/np.pi)**2, label='$\\tilde{\Theta}^{\mathcal{T}}_{HE}(\omega\mathcal{T})$')
ax.plot(x, 0*x, ':k')
ax.set_xticks([0,2*np.pi,4*np.pi])
ax.tick_params('both', bottom=True, top=False, left=True, right=False, direction='in')
ax.xaxis.set_major_formatter(FuncFormatter(
   lambda val,pos: '{:.0g}$\pi$'.format(val/np.pi) if val !=0 else '0'
))
ax.xaxis.set_major_locator(MultipleLocator(base=2*np.pi))
ax.set_xlabel("$\omega\mathcal{T}$")
ax.set_yticks([0,1,2])
lista = ["$0$", "$\mathcal{T}$", "$\mathcal{2T}$"]
ax.set_yticklabels(lista)
ax.set_xlim(0,19.8)
ax.set_ylabel("$\\tilde{\Theta}_X$")
ax.legend(frameon=False)
fig.savefig('Theta_GK-HE.pdf', bbox_inches = 'tight', pad_inches = 0)

Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol 2 from STIXGeneral
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol 2 from STIXGeneral
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol 2 from STIXGeneral
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol 2 from STIXGeneral
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode

In [7]:

msd_elc = elc[:,8]*6.0*time*DT
msd_elq = elq[:,8]*6.0*time*DT
msd_dif = dif[:,8]*6.0*time*DT
var_msd_elc = np.sqrt(elc[:,9])*6.0*time*DT
var_msd_elq = np.sqrt(elq[:,9])*6.0*time*DT
var_msd_dif = np.sqrt(dif[:,9])*6.0*time*DT

In [8]:

m1 = 0.059 # slope 1 from fit
q1 = 0.6   # intercept 1 from fit
m2 = 0.059 # slope 2 from fit
q2 = 1.8   # intercept 2 from fit
m3 = 0.0   # slope 3 from fit
q3 = 5.4   # intercept 3 from fit [x30]

In [9]:

plt.plot(time, 6.0*(m1*time + q1), '--', color='tab:blue', )
plt.plot(time, 6.0*(m2*time + q2), '--', color='tab:red')
plt.plot(time, 6.0*(m3*time + q3), '--', color='tab:green')
plt.plot(time, msd_elc)
plt.plot(time, msd_elq)
plt.plot(time, msd_dif*30)
plt.fill_between(time, msd_elc + var_msd_elc, msd_elc - var_msd_elc, alpha=0.4, label='$\mathbf{J^\prime}$')
plt.fill_between(time, msd_elq + var_msd_elq, msd_elq - var_msd_elq, alpha=0.4, label='$\mathbf{J}$')
plt.fill_between(time, 30.*(msd_dif + var_msd_dif), 30.*(msd_dif - var_msd_dif), alpha=0.4, label='$\mathbf{J^\prime} - \mathbf{J}$ [$\\times 30$]')
plt.xlim(0,300)
plt.ylim(0,120)
plt.xlabel('$\mathcal{T}$ [fs]')
plt.ylabel('$\\langle |\Delta \mathbf{\\mu}|^2 \\rangle$ [$\AA^2$]')
plt.tick_params('both', bottom=True, top=False, left=True, right=False, direction='in')
plt.legend(frameon=False)
plt.xticks([0,100,200,300])
plt.yticks([0,40,80,120])
plt.savefig('msdVStime.pdf', bbox_inches = 'tight', pad_inches = 0)

Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode

Graphical abstract¶

In [10]:

plt.plot(time, 6.0*(m1*time + q1), '--', color='tab:blue', )
plt.plot(time, 6.0*(m2*time + q2), '--', color='tab:red')
#plt.plot(time, 6.0*(0.0*time + 5.4), '--', color='tab:green')
plt.plot(time, msd_elc)
plt.plot(time, msd_elq)
#plt.plot(time, msd_dif*30)
plt.fill_between(time, msd_elc + var_msd_elc, msd_elc - var_msd_elc, alpha=0.4, 
                 label='$\mathbf{J^\prime}$')
plt.fill_between(time, msd_elq + var_msd_elq, msd_elq - var_msd_elq, alpha=0.4, 
                 label='$\mathbf{J} \\neq \mathbf{J}^\prime$')
#plt.fill_between(time, 30.*(msd_dif + var_msd_dif), 30.*(msd_dif - var_msd_dif), alpha=0.4, label='$\mathbf{J^\prime} - \mathbf{J}$ [$\\times 30$]')
plt.xlim(0,300)
plt.ylim(0,120)
plt.xlabel('$\mathcal{T}$ [fs]')
plt.ylabel('$\\langle |\Delta \mathbf{\\mu}|^2 \\rangle$ [$\AA^2$]')
plt.tick_params('both', bottom=True, top=False, left=True, right=False, direction='in')
plt.legend(frameon=False)
plt.xticks([0,100,200,300])
plt.yticks([0,40,80,120])

plt.text(180, 40, '$\\sigma^\prime = \\sigma$',
         fontsize=35)

plt.savefig('graphical.pdf', bbox_inches = 'tight', pad_inches = 0)

Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode

Nota bene: the variance given by the program "analisi" is the variance on the mean value, i.e. $1/B$ the variance, being $B$ the number of blocks. We have $\\ell=3$ components (columns of the original files of the time series of the fluxes).

In [11]:

blocchi = 40
cols = 3

plt.plot(np.arange(len(elc[:,9]))/(len(elc[:,9])-1), 1e3*4/cols*DT**2*elc[26,8]**2*(np.arange(len(elc[:,9]))/(len(elc[:,9])-1)), ':', color='black', 
         label='$\\frac{4}{l}\lambda^2 \,(\mathcal{T}\, / \,\mathcal{T}_{tot})$')

plt.plot(np.arange(len(elc[:,9]))/(len(elc[:,9])-1), 1e3*DT**2*(elc[:,7])*(blocchi), label='GK')

plt.plot((np.arange(len(elc[:,9])))/(len(elc[:,9])-1), 1e3*DT**2*(elc[:,9])*(blocchi), label='HE')

plt.plot((np.arange(len(elc[:,9])))/(len(elc[:,9])-1), 1e3*DT**2*(3*elc[:,9])*(blocchi), '--', color='tab:orange', label='3 HE')


plt.xlim(0,1)
plt.ylim(0,6)
plt.xlabel('$\mathcal{T}\,/\,\mathcal{T}_{tot}$')
plt.ylabel('var$({\lambda}_{X})$ [$10^{-3} \AA^4$/fs$^2$]')
plt.tick_params('both', bottom=True, top=False, left=True, right=False, direction='in')
plt.legend(frameon=False)
plt.yticks([0,2,4,6])
plt.xticks([0.5,1])
plt.savefig('varianze_app_40B.pdf', bbox_inches = 'tight', pad_inches = 0)

Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode

This is very noisy. Futhermore, the theoretical estimate was shown to be an upper bound for the variance on the GK integral (see Jones and Mandadapu, JCP 2012, cited ref.). With a larger number of blocks and thus more statistics the correct theoretical value is obtained.

In [12]:

elc = np.loadtxt('GKHE_class_120B.dat', skiprows=13)
elq = np.loadtxt('GKHE_quant_120B.dat', skiprows=13)
dif = np.loadtxt('GKHE_diff_120B.dat', skiprows=13)

In [13]:

blocchi = 120
cols = 3

plt.plot(np.arange(len(elc[:,9]))/(len(elc[:,9])-1), 1e3*4/cols*DT**2*elc[26,8]**2*(np.arange(len(elc[:,9]))/(len(elc[:,9])-1)), ':', color='black', 
         label='$\\frac{4}{l}\lambda^2 \,(\mathcal{T}\, / \,\mathcal{T}_{tot})$')

plt.plot(np.arange(len(elc[:,9]))/(len(elc[:,9])-1), 1e3*DT**2*(elc[:,7])*(blocchi), label='GK')

plt.plot((np.arange(len(elc[:,9])))/(len(elc[:,9])-1), 1e3*DT**2*(elc[:,9])*(blocchi), label='HE')

plt.plot((np.arange(len(elc[:,9])))/(len(elc[:,9])-1), 1e3*DT**2*(3*elc[:,9])*(blocchi), '--', color='tab:orange', label='3 HE')


plt.xlim(0,1)
plt.ylim(0,6)
plt.xlabel('$\mathcal{T}\,/\,\mathcal{T}_{tot}$')
plt.ylabel('var$({\lambda}_{X})$ [$10^{-3} \AA^4$/fs$^2$]')
plt.tick_params('both', bottom=True, top=False, left=True, right=False, direction='in')
plt.legend(frameon=False)
plt.yticks([0,2,4,6])
plt.xticks([0.5,1])
plt.savefig('varianze_app_40B.pdf', bbox_inches = 'tight', pad_inches = 0)

Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode
Substituting symbol T from STIXNonUnicode

Graphical abstract¶

Appendix¶