We check the convergence with $N_\text{kpt}$ for the calculation of the vacancy Green function for FCC and HCP structures. In particular, we will look at:
with increasing k-point density.
In [1]:
import sys
sys.path.extend(['../'])
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
%matplotlib inline
import onsager.crystal as crystal
import onsager.GFcalc as GFcalc
Create an FCC and HCP lattice.
In [2]:
a0 = 1.
FCC, HCP = crystal.Crystal.FCC(a0, "fcc"), crystal.Crystal.HCP(a0, chemistry="hcp")
print(FCC)
print(HCP)
We will put together our vectors for consideration:
and our sitelists and jumpnetworks.
In [3]:
FCCR = np.array([0,2.,2.])
HCPR1, HCPR2 = np.array([4.,0.,0.]), np.array([2.,0.,2*np.sqrt(8/3)])
In [4]:
FCCsite, FCCjn = FCC.sitelist(0), FCC.jumpnetwork(0, 0.75)
HCPsite, HCPjn = HCP.sitelist(0), HCP.jumpnetwork(0, 1.01)
We use $N_\text{max}$ parameter, which controls the automated generation of k-points to iterate through successively denser k-point meshes.
In [5]:
FCCdata = {pmaxerror:[] for pmaxerror in range(-16,0)}
print('kpt\tNkpt\tG(0)\tG(R)\tG diff')
for Nmax in range(1,13):
GFFCC = GFcalc.GFCrystalcalc(FCC, 0, FCCsite, FCCjn, Nmax=Nmax)
Nreduce, Nkpt, kpt = GFFCC.Nkpt, np.prod(GFFCC.kptgrid), GFFCC.kptgrid
for pmax in sorted(FCCdata.keys(), reverse=True):
GFFCC.SetRates(np.ones(1), np.zeros(1), np.ones(1)/12, np.zeros(1), 10**(pmax))
g0,gR = GFFCC(0,0,np.zeros(3)), GFFCC(0,0,FCCR)
FCCdata[pmax].append((Nkpt, g0, gR))
Nkpt,g0,gR = FCCdata[-8][-1] # print the 10^-8 values
print("{k[0]}x{k[1]}x{k[2]}\t".format(k=kpt) +
" {:5d} ({})\t{:.12f}\t{:.12f}\t{:.12f}".format(Nkpt, Nreduce,
g0, gR,g0-gR))
In [6]:
HCPdata = []
print('kpt\tNkpt\tG(0)\tG(R1)\tG(R2)\tG(R1)-G(0)\tG(R2)-G0')
for Nmax in range(1,13):
GFHCP = GFcalc.GFCrystalcalc(HCP, 0, HCPsite, HCPjn, Nmax=Nmax)
GFHCP.SetRates(np.ones(1), np.zeros(1), np.ones(2)/12, np.zeros(2), 1e-8)
g0,gR1,gR2 = GFHCP(0,0,np.zeros(3)), GFHCP(0,0,HCPR1), GFHCP(0,0,HCPR2)
Nreduce, Nkpt, kpt = GFHCP.Nkpt, np.prod(GFHCP.kptgrid), GFHCP.kptgrid
HCPdata.append((Nkpt, g0, gR1, gR2))
print("{k[0]}x{k[1]}x{k[2]}\t".format(k=kpt) +
"{:5d} ({})\t{:.12f}\t{:.12f}\t{:.12f}\t{:.12f}\t{:.12f}".format(Nkpt, Nreduce,
g0, gR1, gR2,
g0-gR1, g0-gR2))
First, look at the behavior of the error with $p_\text{max}$(error) parameter. The k-point integration error scales as $N_\text{kpt}^{5/3}$, and we see the $p_\text{max}$ error is approximately $10^{-8}$.
In [7]:
print('pmax\tGinf\talpha (Nkpt^-5/3 prefactor)')
Ginflist=[]
for pmax in sorted(FCCdata.keys(), reverse=True):
data = FCCdata[pmax]
Nk53 = np.array([N**(5/3) for (N,g0,gR) in data])
gval = np.array([g0 for (N,g0,gR) in data])
N10,N5 = np.average(Nk53*Nk53),np.average(Nk53)
g10,g5 = np.average(gval*Nk53*Nk53),np.average(gval*Nk53)
denom = N10-N5**2
Ginf,alpha = (g10-g5*N5)/denom, (g10*N5-g5*N10)/denom
Ginflist.append(Ginf)
print('{}\t{}\t{}'.format(pmax, Ginf, alpha))
Plot the error in the Green function for FCC (at 0, maximum R, and difference between those GF). We extract the infinite value by fitting the error to $N_{\mathrm{kpt}}^{-5/3}$, which empirically matches the numerical error.
In [8]:
# plot the errors from pmax = 10^-8
data = FCCdata[-8]
Nk = np.array([N for (N,g0,gR) in data])
g0val = np.array([g0 for (N,g0,gR) in data])
gRval = np.array([gR for (N,g0,gR) in data])
gplot = []
Nk53 = np.array([N**(5/3) for (N,g0,gR) in data])
for gdata, start in zip((g0val, gRval, g0val-gRval), (0,1,2)):
N10,N5 = np.average(Nk53[start:]*Nk53[start:]),np.average(Nk53[start:])
denom = N10-N5**2
g10 = np.average(gdata[start:]*Nk53[start:]*Nk53[start:])
g5 = np.average(gdata[start:]*Nk53[start:])
Ginf,alpha = (g10-g5*N5)/denom, (g10*N5-g5*N10)/denom
gplot.append(np.abs(gdata-Ginf))
fig, ax1 = plt.subplots()
ax1.plot(Nk, gplot[0], 'k', label='$G(\mathbf{0})$ error $\sim N_{\mathrm{kpt}}^{-5/3}$')
ax1.plot(Nk, gplot[1], 'b', label='$G(\mathbf{R})$ error $\sim N_{\mathrm{kpt}}^{-5/3}$')
ax1.plot(Nk, gplot[2], 'b--', label='$G(\mathbf{0})-G(\mathbf{R})$ error')
ax1.set_xlim((1e2,2e5))
ax1.set_ylim((1e-11,1))
ax1.set_xscale('log')
ax1.set_yscale('log')
ax1.set_xlabel('$N_{\mathrm{kpt}}$', fontsize='x-large')
ax1.set_ylabel('integration error $G-G^\infty$', fontsize='x-large')
ax1.legend(bbox_to_anchor=(0.6,0.6,0.4,0.4), ncol=1,
shadow=True, frameon=True, fontsize='x-large')
ax2 = ax1.twiny()
ax2.set_xscale('log')
ax2.set_xlim(ax1.get_xlim())
ax2.set_xticks([n for n in Nk])
ax2.set_xticklabels(["${:.0f}^3$".format(n**(1/3)) for n in Nk])
ax2.set_xlabel('k-point grid', fontsize='x-large')
ax2.grid(False)
ax2.tick_params(axis='x', top='on', direction='in', length=6)
plt.show()
# plt.savefig('FCC-GFerror.pdf', transparent=True, format='pdf')
Plot the error in Green function for HCP.
In [9]:
# plot the errors from pmax = 10^-8
data = HCPdata
Nk = np.array([N for (N,g0,gR1,gR2) in data])
g0val = np.array([g0 for (N,g0,gR1,gR2) in data])
gR1val = np.array([gR1 for (N,g0,gR1,gR2) in data])
gR2val = np.array([gR2 for (N,g0,gR1,gR2) in data])
gplot = []
Nk53 = np.array([N**(5/3) for (N,g0,gR1,gR2) in data])
for gdata, start in zip((g0val, gR1val, gR2val, g0val-gR1val, g0val-gR2val), (3,3,3,3,3)):
N10,N5 = np.average(Nk53[start:]*Nk53[start:]),np.average(Nk53[start:])
denom = N10-N5**2
g10 = np.average(gdata[start:]*Nk53[start:]*Nk53[start:])
g5 = np.average(gdata[start:]*Nk53[start:])
Ginf,alpha = (g10-g5*N5)/denom, (g10*N5-g5*N10)/denom
gplot.append(np.abs(gdata-Ginf))
fig, ax1 = plt.subplots()
ax1.plot(Nk, gplot[0], 'k', label='$G(\mathbf{0})$ error $\sim N_{\mathrm{kpt}}^{-5/3}$')
ax1.plot(Nk, gplot[1], 'b', label='$G(\mathbf{R}_1)$ error $\sim N_{\mathrm{kpt}}^{-5/3}$')
ax1.plot(Nk, gplot[2], 'r', label='$G(\mathbf{R}_2)$ error $\sim N_{\mathrm{kpt}}^{-5/3}$')
ax1.plot(Nk, gplot[3], 'b--', label='$G(\mathbf{0})-G(\mathbf{R}_1)$ error')
ax1.plot(Nk, gplot[4], 'r--', label='$G(\mathbf{0})-G(\mathbf{R}_2)$ error')
ax1.set_xlim((1e2,2e5))
ax1.set_ylim((1e-11,1))
ax1.set_xscale('log')
ax1.set_yscale('log')
ax1.set_xlabel('$N_{\mathrm{kpt}}$', fontsize='x-large')
ax1.set_ylabel('integration error $G-G^\infty$', fontsize='x-large')
ax1.legend(bbox_to_anchor=(0.6,0.6,0.4,0.4), ncol=1,
shadow=True, frameon=True, fontsize='medium')
ax2 = ax1.twiny()
ax2.set_xscale('log')
ax2.set_xlim(ax1.get_xlim())
ax2.set_xticks([n for n in Nk])
# ax2.set_xticklabels(["${:.0f}$".format((n*1.875)**(1/3)) for n in Nk])
ax2.set_xticklabels(['6','10','16','20','26','30','36','40','46','50','56','60'])
ax2.set_xlabel('k-point divisions (basal)', fontsize='x-large')
ax2.grid(False)
ax2.tick_params(axis='x', top='on', direction='in', length=6)
plt.show()
# plt.savefig('HCP-GFerror.pdf', transparent=True, format='pdf')