Morten: I hope you had a nice flight and that you quickly adjust to the new day-night cycle!
I am unable to reproduce Gustavs result for the N=14, r_s = 1.0 case electron gas, and I have tried several approaches but end up with similar results.
My main code is written in C++ on GitHub, but the scripting below reproduces the results from the C++ code. Could you just read through this notebook and insert your comments where you see fit. Also, please feel free to alter any code as you like.
Mvh Audun
My first aim this semester is to set up the interactions for the electron gas and benchmark the reference energy of these expressions with G.Baardsens results. I will therefore first need the expressions for the one- and two body contributions to the energy. Thereafter, I need to initialize a given number of states up to a certain energy level, and finally use this states and their related expressions to calculate the Hartree-Fock energy (reference energy) for comparison with G.Baardsens results. This reference energy will not include the part containing the Madelung constant.
We will use atomic units, meaning that
$$\hbar = m_e = c = e = r_b = 1 $$The electron gas will consist of a periodic lattice of cubic cells, each with a volume
$$\Omega = L^3 = \frac{4}{3}\pi N r_s^3$$where N is the number of occupied orbitals and $r_s$ is a radius assosciated with each constituent particle.
An in-depth overview of this system is given in G.Baardsens dissertation.
For the one body expression (single particle energy), we have
$$ \langle p | h |q \rangle = \delta_{p,q}\frac{\hbar^2 }{2 m} \hat{k}_p^2 = \delta_{p,q}\frac{2 \hbar^2 }{m} (\frac{\pi}{L})^2 (n_x^2 + n_y^2 + n_z^2)_p $$where in atomic units $\hbar = m = c = e = 1$, $\hat{k}$ is the k-vector, $L^3 \equiv \Omega$ is the volume of the box, and the integers $n_i$ denotes the quantized energy levels along each axis.
The antisymmetrized interaction contribution to the hamiltonian is
$$ \langle p q || r s \rangle = \frac{4 \pi c^2}{\Omega} \delta_{k_p + k_q, k_r + k_s} [\delta_{ms_p,ms_r}\delta_{ms_q, ms_s} (1-\delta_{k_p,k_r})\frac{1}{|\hat{k}_r - \hat{k}_p|^2} - \delta_{ms_p,ms_s}\delta_{ms_q, ms_r} (1-\delta_{k_p,k_s})\frac{1}{|\hat{k}_s - \hat{k}_p|^2}] $$Question for Morten:
In Gustavs dissertation the expression for the two body interaction is given as above. In the description of my master thesis, there is however a parameter $\mu$:
$$ \langle p q || r s \rangle = \frac{4 \pi c^2}{\Omega} \delta_{k_p + k_q, k_r + k_s} [\delta_{ms_p,ms_r}\delta_{ms_q, ms_s} (1-\delta_{k_p,k_r})\frac{1}{\mu^2 + (\hat{k}_r - \hat{k}_p)^2} - \delta_{ms_p,ms_s}\delta_{ms_q, ms_r} (1-\delta_{k_p,k_s})\frac{1}{\mu^2 + (\hat{k}_s - \hat{k}_p)^2}] $$It does not say what this parameter is supposed to be and I cannot find it in other similar expressions. Is this a typo?
It should be straight forward to implement the above expressions and calculate the reference energy, but I have written multiple codes in both Python and C++ and fail to reproduce Gustav's results. Before I do the actual calculations I will try to gain some insight by considering some of the expressions by hand.
Since the basis is already minimized, there is no need for any HF-iterations. This means that the coefficient matrix for the HF-expansion is the identity matrix. We therefore only need to consider two-body interactions of the form $\langle i j || i j \rangle $.
This simplifies the two-body expression, since
$$ \langle i j || i j \rangle = \frac{4 \pi c^2}{\Omega} \delta_{k_i + k_j, k_i + k_j} [\delta_{ms_i,ms_i}\delta_{ms_j, ms_j} (1-\delta_{k_i,k_i})\frac{1}{|\hat{k}_i - \hat{k}_i|^2} - \delta_{ms_i,ms_j}\delta_{ms_j, ms_i} (1-\delta_{k_i,k_j})\frac{1}{|\hat{k}_j - \hat{k}_i|^2}] $$so that the direct term always cancel in the kroenecker delta $(1-\delta_{k_i,k_i})$:
$$ \langle i j || i j \rangle = - \delta_{ms_i,ms_j} \frac{4 \pi c^2}{\Omega}\frac{1}{|\hat{k}_j - \hat{k}_i|^2} $$In summary, we then expect (for the reference energy) that we only get contributions from the exchange part, where i and j are in parallel spin alignment. Morten: Is this true?
In [27]:
from numpy import *
class electronbasis():
def __init__(self, N, rs, Nparticles):
self.rs = rs
self.states = []
self.nstates = 0
self.nparticles = Nparticles
Nm = int(sqrt(N) + 1)
#Creating the basis
for x in range(-Nm, Nm):
for y in range(-Nm, Nm):
for z in range(-Nm,Nm):
e = x*x + y*y + z*z
if e <=N:
self.states.append([e, x,y,z, 1])
self.states.append([e, x,y,z,-1])
self.nstates += 2
self.states.sort() #Sorting the basis in increasing energy
self.L3 = (4*pi*self.nparticles*self.rs**3)/3.0
self.L2 = self.L3**(2/3.0)
self.L = pow(self.L3, 1/3.0)
for i in range(self.nstates):
self.states[i][0] *= 2*(pi**2)/self.L**2 #Multiplying in the missing factors in the single particle energy
self.states = array(self.states) #converting to array to utilize vectorized calculations
def hfenergy(self, nParticles):
#Calculating the HF-energy (reference energy)
e0 = 0.0
if nParticles<=self.nstates:
for i in range(nParticles):
e0 += self.h(i,i)
for j in range(nParticles):
if j != i:
e0 += .5*self.v(i,j,i,j)
else:
#Safety for cases where nParticles exceeds size of basis
print "Not enough basis states."
return e0
def h(self, p,q):
#Return single particle energy
return self.states[p,0]*(p==q)
def veval(self, p,q,r,s):
#A test for evaluating the two-body interaction
val = ""
if self.kdplus(p,q,r,s):
val+= "kdplus "
if self.kdspin(p,r):
val += "Direct[kdspin_pr "
if self.kdspin(q,s):
val += "kdspin_qs "
if self.kdwave(p,r) != 0:
val += "kdwave!=0 "
val += str(self.absdiff2(r,p))
val += "] "
if self.kdspin(p,s):
val += "Exchange[kdspin_pr "
if self.kdspin(q,r):
val += "kdspin_qs "
if self.kdwave(p,s) != 0:
val += "kdwave!=0 "
val += str(self.absdiff2(s,p))
val += "] "
return val
def vevalHF(self, N):
#Evaluation of all expressions of two-body contributions to the HF-energy
for i in range(N):
for j in range(N):
if i!= j:
print "<",i,j,"|",i,j,"> =",self.veval(i,j,i,j)
def v(self,p,q,r,s):
#Two body interaction
val = self.kdplus(p,q,r,s)*4*pi/self.L3
terms = 0.0
term1 = 0.0
term2 = 0.0
if self.kdspin(p,r)*self.kdspin(q,s)==1:
if self.kdwave(p,r) != 1.0:
term1=(4*self.absdiff2(r,p)*pi**2)/self.L2
terms += 1.0/term1
if self.kdspin(p,s)*self.kdspin(q,r)==1:
if self.kdwave(p,s) != 1.0:
term2=(4*self.absdiff2(s,p)*pi**2)/self.L2
terms -= 1.0/term2
return val*terms
#The following is a series of kroenecker deltas used in the two-body interactions.
#Run kd_integrity() to ensure that they work as intended.
def kdi(self,a,b):
#Kroenecker delta integer
return 1.0*(a==b)
def kda(self,a,b):
#Kroenecker delta array
d = 1.0
#print a,b,
for i in range(len(a)):
d*=(a[i]==b[i])
return d
def kdplus(self,p,q,r,s):
#Kroenecker delta wavenumber p+q,r+s
return self.kda(self.states[p][1:4]+self.states[q][1:4],self.states[r][1:4]+self.states[s][1:4])
def kdspin(self,p,q):
#Kroenecker delta spin
return self.kdi(self.states[p][4], self.states[q][4])
def kdwave(self,p,q):
#Kroenecker delta wavenumber
return self.kda(self.states[p][1:4],self.states[q][1:4])
def absdiff2(self,p,q):
val = 0.0
for i in range(1,4):
val += (self.states[p][i]-self.states[q][i])*(self.states[p][i]-self.states[q][i])
#if val == 0:
# print "div0"
return val
def kd_integrity(self):
#test integrity of kroenecker deltas
print "Array KD :", self.kda([0,1,2], [0,1,2]) == True
print "Integer KD :", self.kdi(1,1) == True
print "Opposite spin :", self.kdspin(0,1) == False
print "Equal spin :", self.kdspin(1,1) == True
print "Wavenumber equal :", self.kdwave(1,0) == True
print "Wavenumber not equal:", self.kdwave(1,2) == False
In [28]:
bs = electronbasis(3,1.0, 14) #max energy (n_x^2 + n_y^2 + n_z^2) <= 1, r_S = 1.0
print bs.states #in increasing order of energy, k_x, k_Y, k_z, m_s
The one body contribution to the HF-energy should now only be the sum of all the elements in the first column:
In [29]:
e_sp = sum(bs.states[0:14,0])
print "Single particle energy:", e_sp, "a.u."
The two body contribution is found by evaluating:
In [30]:
e_v = 0.0
for i in range(14):
for j in range(14):
if i!=j:
e_v += .5*bs.v(i,j,i,j)
print "Interaction energy :", e_v, "a.u."
print "Total energy :", e_v + e_sp, " a.u."
print "Total energy per particle:", 2*(e_v + e_sp)/14.0, " rydberg"
print "Gustav Baardsens results :", "1.93434 rydberg"
As can be seen from the above results, this deviates from Gustavs results considerably.
I get the same results in C++/FermionMingle, so this seems to be an error related to the way I have implemented the equations rather than a trivial bug.
Note (13.02.2015) The bug was found to be that an extra factor of 2 in the interaction energy. The above code now yields results very close to Gustavs results, and the C++ code (Www.github.com/audunsh/fermicc) is exactly the same.
In [31]:
print bs.v(0,2,0,2)!=0
print bs.v(0,1,0,1)==0
print bs.v(1,1,1,1)==0
print bs.v(1,3,1,3)!=0
print bs.v(3,1,1,3)==bs.v(1,3,3,1)
print bs.v(3,5,3,5)!=0
print bs.v(4,6,4,6)!=0
print bs.v(3,4,3,4)==0
bs.kd_integrity()
print bs.vevalHF(14)
In [ ]: