# Import Bokeh modules for interactive plotting
import bokeh.io
import bokeh.mpl
import bokeh.plotting
# Magic function to make matplotlib inline; other style specs must come AFTER
%matplotlib inline
#%matplotlib notebook
# This enables SVG graphics inline. There is a bug, so uncomment if it works.
%config InlineBackend.figure_formats = {'svg',}
# This enables high resolution PNGs. SVG is preferred, but has problems
# rendering vertical and horizontal lines
# %config InlineBackend.figure_formats = {'png', 'retina'}
# Set up Bokeh for inline viewing
bokeh.io.output_notebook()
We will load all calculated correlation energies obtained by G0W0, and than we will fit them with some analytic expression.
First, we vary the dielectric function $\varepsilon$. We load the data from the disc and display them.
from scipy import *
epsx = [1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0] # our clculated dielectric functions
ec=[]
for eps in epsx:
file = 'work/ec_eps_'+str(eps)+'_lam_0.0.dat'
dat = loadtxt(file).transpose()
rs = dat[0]
ec.append(dat[1])
ec = array(ec)
# Here we take the ratio ec/ec(lambda=1).
for i in range(1,len(ec)):
ec[i,:] *= 1/ec[0,:]
ec[0,:] *= 1/ec[0,:]
from pylab import *
%matplotlib inline
colors=['b','r','m','k','g','c','y','b','r']
for i in range(len(ec)):
plot(rs, ec[i,:], label='$\\varepsilon='+str(epsx[i])+'$')
legend(loc='best')
xlabel('$r_s$', fontsize='x-large')
ylabel('$E_c/E_c(\\varepsilon=1)$',fontsize='x-large')
grid()
The ratio of the correlation energies ($E_c$ is the correlation energy density, not to be confused with dielectric function $\varepsilon$) can be well fited by the function
\begin{eqnarray} \frac{E_c(\varepsilon,\lambda=0)}{E_c(\varepsilon=1,\lambda=0)} =\frac{1+b}{\varepsilon^a + b \varepsilon^d} \end{eqnarray}where
\begin{eqnarray} && a = a_0 + a_1 r_s + a_2 r_s^2\\ && b = 10^{-3} \times \frac{b_0 \sqrt{r_s}+ b_1 r_s}{1+ (b_2 r_s)^9}\\ && d = d_0 + d_1 r_s + d_2 r_s^2 + d_3 r_s^3 \end{eqnarray}and
\begin{eqnarray} && (a_0,a_1,a_2) = (1.74596971, -0.0892907, 0.00658866)\\ && (b_0, b_1, b_2) = (1.63289109, 1.15291480, 0.149402)\\ && (d_0,d_1,d_2,d_3) = (3.64370598, 0.03636027, -0.03886317, 0.00693599) \end{eqnarray}This is the expression coded below.
To get some rought approximation how $E_c$ depends on $\epsilon$, we take an average value (over $r_s$) of parameters $a$, $b$ and $d$ to obtain:
\begin{eqnarray} \frac{E_c}{E_c(\varepsilon=1)} =\frac{1.005}{\varepsilon^{1.6} + 0.005 \varepsilon^{3.6}} \end{eqnarray}which shows that $E_c$ falls off rughly $1/\varepsilon^{1.6}$, which is faster then exchange $1/\varepsilon$ but slower than the second order diagram ($1/\varepsilon^2$).
def EC_ratio_dielec(rs,eps):
# The coefficients a, b, d depend on rs, and here is their more precise fit:
ca = [ 1.74596971, -0.0892907, 0.00658866]
cb = [ 1.63289109, 1.15291480, 0.149402]
cd = [3.64370598, 0.03636027, -0.03886317, 0.00693599]
a = ca[0] + rs * ca[1] + rs**2 * ca[2]
b = 0.001*(cb[0]*sqrt(rs)+cb[1]*rs)/(1+ (cb[2]*rs)**9)
d = cd[0] + rs * cd[1] + rs**2 * cd[2] + rs**3 * cd[3]
return (1.+b)/(eps**a + b*eps**d)
Next we check how good is the fit. We plot the data with dots, and fit with lines.
%matplotlib inline
for i in range(len(rs)):
plot(epsx, ec[:,i], 'o'+colors[i], label='rs='+str(rs[i]))
eps_dense = linspace(epsx[0],epsx[-1],100)
efit = EC_ratio_dielec(rs[i],eps_dense)
plot(eps_dense, efit, '-'+colors[i])
legend(loc='best',numpoints=1)
xlabel('$\\varepsilon$',fontsize='x-large')
ylabel('$E_c/E_c(\\varepsilon=1)$',fontsize='x-large')
show()
for i in range(len(epsx)):
plot(rs, ec[i,:], 'o'+colors[i], label='$\\varepsilon='+str(rs[i])+'$')
efit = EC_ratio_dielec(rs,epsx[i])
plot(rs, efit, '-'+colors[i])
legend(loc='best',numpoints=1)
xlabel('$r_s$',fontsize='x-large')
ylabel('$E_c/E_c(\\varepsilon=1)$',fontsize='x-large')
show()
We next load the data, which show the dependence on yukawa screening length $\lambda$.
from scipy import *
lamx = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.75, 1.0, 1.25, 1.5, 2.0, 2.5, 3.0]
eps=1.0
ec=[]
for lam in lamx:
file = 'work/ec_eps_'+str(eps)+'_lam_'+str(lam)+'.dat'
dat = loadtxt(file).transpose()
rs = dat[0]
ec.append(dat[1])
ec = array(ec)
for i in range(1,len(ec)):
ec[i,:] *= 1/ec[0,:]
ec[0,:] *= 1/ec[0,:]
Next we fit the dependence on $\lambda$. We obtain the following fit:
\begin{eqnarray} \frac{E_c(\varepsilon=1,\lambda)}{E_c(\varepsilon=1,\lambda=0)} =\exp\left({a_0+a_1 r_s + a_2 r_s^2+a_3 r_s^3+a_4 r_s^4+a_5 r_s^5}\right)\; (1-\alpha+\beta r_s^2) + (\alpha-\beta r_s^2) \end{eqnarray}with
\begin{eqnarray} && \alpha=0.008\\ && \beta=0.00112\\ && a_n = \lambda C_{n,0} + \lambda^2 C_{n,1} + \lambda^3 C_{n,2} + \lambda^4 C_{n,3} + \lambda^5 C_{n,4}+\lambda^6 C_{n,5}+\lambda^7 C_{n,6} \end{eqnarray}and
\begin{eqnarray} C = \left(\begin{array}{ccccccc} 0.15805009& -0.77391602& 1.23971169& -1.04865383& 0.47809619& -0.11057964& 0.01016968\\ -0.306851& -0.77296572& 0.8791705& -0.69185034& 0.33779654& -0.08858483& 0.00935635\\ 0.13215843& -0.2776552& 0.45727548& -0.31469164& 0.10787374& -0.01661214& 0.0007591\\ -0.03086548& 0.0549528& -0.07252823& 0.04177618& -0.01084882& 0.00062192& 0.0001177\\ 0.00273230889& -0.00357007233& 0.00425309814& -0.00198811211& 0.000233761378& 0.000106803015& -2.50612307e-05\\ -9.28530649e-05& 8.09009085e-05& -9.43747991e-05& 3.89520548e-05& -3.10149723e-07& -4.23041605e-06& 8.02291467e-07 \end{array}\right) \end{eqnarray}Below we implement the yukawa screening and check how good is the fit.
def EC_ratio_yukawa(lamx,rs):
c=array([[0.15805009, -0.77391602, 1.23971169, -1.04865383, 0.47809619, -0.11057964, 0.01016968],
[-0.306851, -0.77296572, 0.8791705, -0.69185034, 0.33779654, -0.08858483, 0.00935635],
[0.13215843, -0.2776552, 0.45727548, -0.31469164, 0.10787374, -0.01661214, 0.0007591],
[-0.03086548, 0.0549528, -0.07252823, 0.04177618, -0.01084882, 0.00062192, 0.0001177],
[0.00273230889, -0.00357007233, 0.00425309814, -0.00198811211, 0.000233761378, 0.000106803015, -2.50612307e-05],
[-9.28530649e-05, 8.09009085e-05, -9.43747991e-05, 3.89520548e-05, -3.10149723e-07, -4.23041605e-06, 8.02291467e-07]])
x=array(lamx)
a = zeros((6,len(lamx)))
for n in range(6):
a[n,:] = x*(c[n,0] + x*(c[n,1] + x*(c[n,2] + x*(c[n,3] + x*(c[n,4] + x*(c[n,5] + c[n,6]*x))))))
def limit(rs):
return -rs**2*0.00112 + 0.008
def dlimit(rs):
return -2*0.00112*rs
rs = array(rs)
fn = zeros((len(lamx),len(rs)))
qn = zeros((len(lamx),len(rs)))
for i in range(len(lamx)):
xte = a[0,i] + rs*(a[1,i] + rs*(a[2,i] + rs*(a[3,i] + rs*(a[4,i] + rs*a[5,i]))))
te = exp(xte)
lmrs = limit(rs)
fn[i,:] = te[:]*(1-lmrs) + lmrs
das = rs*(a[1,i] + rs*(2*a[2,i] + rs*(3*a[3,i] + rs*(4*a[4,i] + rs*5*a[5,i]))))
qn[i,:] = -1./3.*te*(1-lmrs)*das - 1./3.*rs*(1-te)*dlimit(rs)
# Ec = Ev * fn
# Vc = Vc * fn + Ec * qn
return fn
from pylab import *
%matplotlib inline
colors=['b','r','m','k','g','c','y','b','r','m','k','g','c','y']
fn = EC_ratio_yukawa(lamx,rs)
for i in range(len(lamx)):
plot(rs, ec[i,:], 'o'+colors[i], label='$\\lambda='+str(lamx[i])+'$')
plot(rs, fn[i,:], '-'+colors[i])
legend(loc='best',numpoints=1)
grid()
xlabel('$r_s$',fontsize='x-large')
ylabel('$\\varepsilon_c/\\varepsilon_c(\\lambda=0)$',fontsize='x-large')
show()
for i in range(len(rs)):
plot(lamx, ec[:,i], 'o'+colors[i], label='$r_s='+str(rs[i])+'$')
plot(lamx, fn[:,i], '-'+colors[i])
legend(loc='best',numpoints=1)
grid()
xlabel('$\\lambda$',fontsize='x-large')
ylabel('$\\varepsilon_c/\\varepsilon_c(\\lambda=0)$',fontsize='x-large')
show()
In this work we used the total energy (using Migdal-Galitskii formula) to compute exchange correlation energy. We could use part of the Luttinger-Ward functional, which would lead to a slightly different fit, which was used in PRL 115, 196403 (2015). The old fit was:
def EC_ratio_old(rs,lmbda):
def f1(coef,x):
return (x*coef[1]+x**2*coef[2])/(1+x**2*coef[3]+x**4*coef[4]+x**6*coef[5])
def f2(coef,x):
return (x**2*coef[1]+x**3*coef[2])/(1+x**2*coef[3]+x**4*coef[4])
def f3(coef,x):
return (x**3*coef[1]+x**4*coef[2])/(1+coef[3]*x**2)
def f4(coef,x):
return x**4*(coef[1]+coef[2]*x**2)
coef1 = [0, 0.12238912, 0.73648662, 0.96044695, -0.07501634, 0.00207808]
coef2 = [0, 0.05839362, 0.11969474, 0.10156124, 0.01594125]
coef3 = [0, 0.00827519, 0.00557133, 0.01725079]
coef4 = [0, 5.29134419e-04, 4.49628225e-06]
cfs = zeros(5)
cfs[1] = exp(f1(coef1,lmbda))-1.
cfs[2] = exp(f2(coef2,lmbda))-1.
cfs[3] = exp(f3(coef3,lmbda))-1.
cfs[4] = exp(f4(coef4,lmbda))-1.
A = 1.0 + rs*(cfs[1]+rs*(cfs[2]+rs*(cfs[3]+rs*cfs[4])))
return 1/A
Here we compare the two different fits, to show that they mostly agree, in particular in the relavent range of small $r_s$.
%matplotlib inline
colors=['b','r','m','k','g','c','y','b','r','m','k','g','c','y']
fn = EC_ratio_yukawa(lamx,rs)
# This is obtained with different scheme of computing E_xc within G0W0
for i in range(len(lamx)):
fo = EC_ratio_old(rs,lamx[i])
plot(rs, fo, '--'+colors[i])
plot(rs, fn[i,:], '-'+colors[i],label='$\\lambda='+str(lamx[i])+'$')
legend(loc='best')
grid()
xlabel('$r_s$',fontsize='x-large')
ylabel('$\\varepsilon_c/\\varepsilon_c(\\lambda=0)$',fontsize='x-large')
show()
Finally, for completeness we also show how to compute the correlation potential, which is obtained from correlation density by functional derivative
\begin{eqnarray} V_c = \frac{\delta}{\delta\rho}\int d^3r\; E_c(r)\; \rho(r) = E_c(r) + \rho \frac{\delta E_c(r)}{\delta\rho(r)} \end{eqnarray}If the correlation energy is screened by
\begin{eqnarray} E_c = f(r_s) \; E_c^0 \end{eqnarray}we obtain for $V_c$ the following formula:
\begin{eqnarray} V_c = f(r_s) \left( E_c^0 + \rho \frac{\delta E^0_c(r)}{\delta\rho(r)}\right) + E_c^0 \rho \frac{\delta r_s}{\delta\rho(r)}\;\frac{d f(r_s)}{d r_s}= f(r_s)\; V_c^0 -\frac{1}{3}r_s \frac{d f(r_s)}{d r_s}\; E_c^0 \end{eqnarray}Next we define:
$$g(r_s)\equiv -\frac{1}{3}r_s \frac{d f(r_s)}{d r_s}$$so that we have
\begin{eqnarray} && E_c = f(r_s)\; E_c^0\\ && V_c = f(r_s)\;V_c^0 + g(r_s)\; E_c^0 \end{eqnarray}Below we implement both $f(r_s)$ and $g(r_s)$.
def EC_ratio_yukawa(lamx,rs):
c=array([[0.15805009, -0.77391602, 1.23971169, -1.04865383, 0.47809619, -0.11057964, 0.01016968],
[-0.306851, -0.77296572, 0.8791705, -0.69185034, 0.33779654, -0.08858483, 0.00935635],
[0.13215843, -0.2776552, 0.45727548, -0.31469164, 0.10787374, -0.01661214, 0.0007591],
[-0.03086548, 0.0549528, -0.07252823, 0.04177618, -0.01084882, 0.00062192, 0.0001177],
[0.00273230889, -0.00357007233, 0.00425309814, -0.00198811211, 0.000233761378, 0.000106803015, -2.50612307e-05],
[-9.28530649e-05, 8.09009085e-05, -9.43747991e-05, 3.89520548e-05, -3.10149723e-07, -4.23041605e-06, 8.02291467e-07]])
x=array(lamx)
a = zeros((6,len(lamx)))
for n in range(6):
a[n,:] = x*(c[n,0] + x*(c[n,1] + x*(c[n,2] + x*(c[n,3] + x*(c[n,4] + x*(c[n,5] + c[n,6]*x))))))
def limit(rs):
return -rs**2*0.00112 + 0.008
def dlimit(rs):
return -2*0.00112*rs
rs = array(rs)
fn = zeros((len(lamx),len(rs)))
qn = zeros((len(lamx),len(rs)))
for i in range(len(lamx)):
xte = a[0,i] + rs*(a[1,i] + rs*(a[2,i] + rs*(a[3,i] + rs*(a[4,i] + rs*a[5,i]))))
te = exp(xte)
lmrs = limit(rs)
fn[i,:] = te[:]*(1-lmrs) + lmrs
das = rs*(a[1,i] + rs*(2*a[2,i] + rs*(3*a[3,i] + rs*(4*a[4,i] + rs*5*a[5,i]))))
qn[i,:] = -1./3.*te*(1-lmrs)*das - 1./3.*rs*(1-te)*dlimit(rs)
# Ec = Ev * fn
# Vc = Vc * fn + Ec * qn
return (fn,qn)
def EC_ratio_dielec(rs,eps):
# The coefficients a, b, d depend on rs, and here is their more precise fit:
ca = [ 1.74596971, -0.0892907, 0.00658866]
cb = [ 1.63289109, 1.15291480, 0.149402]
cd = [3.64370598, 0.03636027, -0.03886317, 0.00693599]
a = ca[0] + rs * ca[1] + rs**2 * ca[2]
b = 0.001*(cb[0]*sqrt(rs)+cb[1]*rs)/(1+ (cb[2]*rs)**9)
d = cd[0] + rs * cd[1] + rs**2 * cd[2] + rs**3 * cd[3]
fn = (1.+b)/(eps**a + b*eps**d)
r_da_dr = rs * ca[1] + 2*rs**2 * ca[2]
r_dd_dr = rs * cd[1] + 2*rs**2 * cd[2] + 3*rs**3 * cd[3]
eps_a = eps**a
eps_d = eps**d
eps_abd = eps_a+b*eps_d
b2r_9 = (cb[2]*rs)**9
r_db_dr = 0.001*(cb[0]*sqrt(rs)*(0.5- 17./2.*b2r_9)+cb[1]*rs*(1-b2r_9*8))/(1+ b2r_9)**2
qn = -1./3.*(r_db_dr*(eps_a-eps_d)-(eps_a*r_da_dr + b*eps_d*r_dd_dr)*(1+b)*log(eps))/eps_abd**2
# Ec = Ev * fn
# Vc = Vc * fn + Ec * qn
return (fn,qn)