Dear all,
I tried to compare a specific matrix element computed in the rf part of abinit with a finite difference approach.
First I choose to look at the contribution of the L q-point on the gamma k point of Diamond. And specifically at the contribution of the LO mode on the first electronic band at gamma. I choose this because of it is a non degenerate mode and a non degenerate band. The corresponding matrix element is:
<n'k+q | dV/dR | nk> = <1L| H(1) | 1G>
1) RF using abinit
File for generating DDB and GKK
Code: Select all
ndtset 10
#
# DATASET 1 : make ground state wavefunctions and density
#
tolwfr1 1.0d-18
nline1 8 # This is to expedite the convergence of higher-lying bands
rfphon1 0 # for DS1 do _not_ do perturbation
nqpt1 0 # for DS1 do _not_ do perturbation
prtwf1 1 # need GS wavefunctions for further runs
getwfk1 0
# enforce calculation of forces at each SCF step
optforces 1
qpt2 0.00000000E+00 0.00000000E+00 0.00000000E+00
qpt3 2.50000000E-01 0.00000000E+00 0.00000000E+00
qpt4 5.00000000E-01 0.00000000E+00 0.00000000E+00
qpt5 2.50000000E-01 2.50000000E-01 0.00000000E+00
qpt6 5.00000000E-01 2.50000000E-01 0.00000000E+00
qpt7 -2.50000000E-01 2.50000000E-01 0.00000000E+00
qpt8 5.00000000E-01 5.00000000E-01 0.00000000E+00
qpt9 -2.50000000E-01 5.00000000E-01 2.50000000E-01
# DS5 get DDK perturbation (stored in GKK matrix element files)
tolwfr10 1.0d-18
qpt10 0 0 0
rfphon10 0
rfelfd10 2
prtwf10 0
# General data for all phonon calculations:
rfatpol 1 2 # Treat displacement of the two atoms.
rfdir 1 1 1 # Do all direction
rfphon 1 # Do phonon response
prtwf 0 # for response function runs,
# do not keep first order wavefunctions
tolvrs 1.0e-12
getwfk 1
nqpt 1
prepgkk 1 # force all perturbations to be calculated for q-points considered
prtgkk 1 # print out GKK files containing electron-phonon coupling
# Common data for GS and perturbation runs
ngkpt 4 4 4
kptopt 3
# use a centered grid for the kpoints: obligatory for el-ph for the moment
nshiftk 1
shiftk 0.0 0.0 0.0
ecut 20
acell 3*6.6709983131
rprim
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
nband 10
# Cell dependent parameters
occopt 1
natom 2
typat 1 1
xred 0.00 0.00 0.00
0.25 0.25 0.25
nstep 300
ntypat 1
znucl 6
diemac 6
enunit 2The merge input files:
Code: Select all
telphon_2.ddb.out
Total ddb for CaS FCC system
8
telphon_1o_DS2_DDB
telphon_1o_DS3_DDB
telphon_1o_DS4_DDB
telphon_1o_DS5_DDB
telphon_1o_DS6_DDB
telphon_1o_DS7_DDB
telphon_1o_DS8_DDB
telphon_1o_DS9_DDBCode: Select all
telphon_3o_GKK.bin # Name of output file
0 # binary (0) or ascii (1) output
telphon_1o_DS1_WFK # GS wavefunction file
0 48 48 # number of 1WF files, of GKK files, and of perturbations in the GKK files
telphon_1o_DS2_GKK1 # names of the 1WF then the (eventual) GKK files
telphon_1o_DS2_GKK2
telphon_1o_DS2_GKK3
telphon_1o_DS2_GKK4
telphon_1o_DS2_GKK5
telphon_1o_DS2_GKK6
telphon_1o_DS3_GKK1
telphon_1o_DS3_GKK2
telphon_1o_DS3_GKK3
telphon_1o_DS3_GKK4
telphon_1o_DS3_GKK5
telphon_1o_DS3_GKK6
telphon_1o_DS4_GKK1
telphon_1o_DS4_GKK2
telphon_1o_DS4_GKK3
telphon_1o_DS4_GKK4
telphon_1o_DS4_GKK5
telphon_1o_DS4_GKK6
telphon_1o_DS5_GKK1
telphon_1o_DS5_GKK2
telphon_1o_DS5_GKK3
telphon_1o_DS5_GKK4
telphon_1o_DS5_GKK5
telphon_1o_DS5_GKK6
telphon_1o_DS6_GKK1
telphon_1o_DS6_GKK2
telphon_1o_DS6_GKK3
telphon_1o_DS6_GKK4
telphon_1o_DS6_GKK5
telphon_1o_DS6_GKK6
telphon_1o_DS7_GKK1
telphon_1o_DS7_GKK2
telphon_1o_DS7_GKK3
telphon_1o_DS7_GKK4
telphon_1o_DS7_GKK5
telphon_1o_DS7_GKK6
telphon_1o_DS8_GKK1
telphon_1o_DS8_GKK2
telphon_1o_DS8_GKK3
telphon_1o_DS8_GKK4
telphon_1o_DS8_GKK5
telphon_1o_DS8_GKK6
telphon_1o_DS9_GKK1
telphon_1o_DS9_GKK2
telphon_1o_DS9_GKK3
telphon_1o_DS9_GKK4
telphon_1o_DS9_GKK5
telphon_1o_DS9_GKK6The anaddb file:
Code: Select all
# turn on calculation of the electron-phonon quantities
elphflag 1
nqpath 1
qpath
0.0 0.0 0.0
ep_prt_yambo 1
# impose acoustic sum rule in a symmetric way
asr 2
dipdip 0 # Dipole-dipole interaction treatment
# bravais lattice necessary
brav 1 # Bravais lattic simple cubic
ngqpt 4 4 4 # Monkhorst Pack indices
nqshft 1 # Number of q-points in repeted basic q-cell
q1shft 0.0 0.0 0.0
ifcflag 1
ifcana 1
# print dielectric matrix with freq dependence
dieflag 0
# print out eigenvectors and symmetrize dyn matrix
eivec 1
#Wavevector list number 1 (Reduced coordinates and normalization factor)
nph1l 0
# This line added when defaults were changed (v5.3) to keep the previous, old behaviour
symdynmat 0All this give me a yambo_elphon_gkk_bymode file that contain the matrix element I'm interested in:
Code: Select all
reduced coord qpoint no 3 0.5000000000E+00 0.0000000000E+00 0.0000000000E+00
matrix elements of all phonon modes for this q-point
kpoint 1
el bands n,np 1 1
mat el for phonon mode num = 1 0.6148462851E-12 0.1910195268E-12
mat el for phonon mode num = 2 0.3322124045E-12 0.6254075303E-12
mat el for phonon mode num = 3 -0.1993910522E-11 0.5519944168E-12
mat el for phonon mode num = 4 0.1166239653E-02 -0.1471911143E-12
mat el for phonon mode num = 5 0.4120494136E-13 0.2338180477E-12
mat el for phonon mode num = 6 -0.4321510669E-11 -0.3949574423E-12We can see that the LO phonon mode give a matrix element of 0.1166239653E-02-0.1471911143E-12j
2) Finite difference using abinit
For the finite difference calculation I used the exact same psp, ecut, acell, conv. criteria etc..
As I'm considering the L phonon mode I made a 2x1x1 supercell with 4 atoms (2 per primitive cell). The k grid is then 2x4x4.
As a first check I computed the phonon frequency by making finite difference on the total energy and obtain very good result:
For the LO phonon I have 0.0589 Ha by finite difference (the convergence over a small displacement parameter (lambda or h) is made for finite difference calculations). This result compare well with the DFPT phonon freq. of LO phonon of 0.005890477 Ha.
Here is my python script used to determine the displacement that I have to made in order to make the finite difference calculation. This script is based on solving the eigenvalue problem based on the knowledge of the dynamical matrix obtained by DFPT.
Code: Select all
import numpy as np
# ------------------------------------------------------------------
# Find the eigenvector and eigenenergies from the dynamical matrix
# computed using DFPT in ABINIT
# ------------------------------------------------------------------
# Location of the Dynamical matrix file
L_file = '/home/Samuel/Documents/WorkDiam/AHC_comparison_with_yambo/0/Dynamical-mat-L.dat'
# Carbon mass in a.u.
Mc=21894.16693
def read_dynamical_matrix(datafile):
with open(datafile,'r') as f:
raw = list()
for line in f:
if line[0]=='#': continue
parts=line.split()
if not parts: continue
val=float(parts[4])
if np.abs(val) < 0.000000001:
val=0.0
raw.append(val)
dynmat = np.array(raw)
dynmat = dynmat.reshape(6,6)
return dynmat
dynmatL=read_dynamical_matrix(L_file)
[eigvalL,eigvectL]=np.linalg.eig(dynmatL)
eigvectL=np.transpose(eigvectL)
print "eigenvalues at L :"
for n in eigvalL:
a = np.sqrt(n/Mc)
print a
print "eigenvectors at L :"
for n in eigvectL: print n
# -------------------------------------------------------------------
# Create Abinit input file for finite difference on eigen energies and other
# -------------------------------------------------------------------
# Value and number of the displacement (always 1/2 of the previous)
h = np.array([0.1,0.05,0.025,0.0125,0.00625,0.003125,0.0015625])
# Cell parameters
acell = np.array([13.336533,6.6682665177,6.6682665177])
rprim = np.array([[0.0,0.5,0.5],[0.5,0.0,0.5],[0.5,0.5,0.0]])
cell = np.mat([rprim[0]*acell[0],rprim[1]*acell[1],rprim[2]*acell[2]])
xred = np.mat([[0.0,0.0,0.0],[0.125,0.25,0.25],[0.5,0.0,0.0],\
[0.625,0.25,0.25]])
# In the specific case of the L point, h = -h
# Indice that record the number of xcart in the input file
j = 1
# Equilibrium position
print " "
print "Contribution of the L point on the electronic eigenenergies at Gamma"
print "and the first part of the NDDW."
print "Using finite difference."
print " "
print " xcart%.0f %.15f %.15f %.15f "%(j,(xred[0]*cell)[0,0],\
(xred[0]*cell)[0,1],(xred[0]*cell)[0,2])
print " %.15f %.15f %.15f "%((xred[1]*cell)[0,0],\
(xred[1]*cell)[0,1],(xred[1]*cell)[0,2])
print " %.15f %.15f %.15f "%((xred[2]*cell)[0,0],\
(xred[2]*cell)[0,1],(xred[2]*cell)[0,2])
print " %.15f %.15f %.15f "%((xred[3]*cell)[0,0],\
(xred[3]*cell)[0,1],(xred[3]*cell)[0,2])
print " "
j += 1
for disp in h:
# LO mode
i = 1
print " xcart%.0f %.15f %.15f %.15f "%(j,(xred[0]*cell)[0,0]+\
eigvectL[i][0]*disp,(xred[0]*cell)[0,1]+eigvectL[i][1]*disp,\
(xred[0]*cell)[0,2]+eigvectL[i][2]*disp)
print " %.15f %.15f %.15f "%((xred[1]*cell)[0,0]+\
eigvectL[i][3]*disp,(xred[1]*cell)[0,1]+eigvectL[i][4]*disp,\
(xred[1]*cell)[0,2]+eigvectL[i][5]*disp)
print " %.15f %.15f %.15f "%((xred[2]*cell)[0,0]-\
eigvectL[i][0]*disp,(xred[2]*cell)[0,1]-eigvectL[i][1]*disp,\
(xred[2]*cell)[0,2]-eigvectL[i][2]*disp)
print " %.15f %.15f %.15f "%((xred[3]*cell)[0,0]-\
eigvectL[i][3]*disp,(xred[3]*cell)[0,1]-eigvectL[i][4]*disp,\
(xred[3]*cell)[0,2]-eigvectL[i][5]*disp)
print " "
j += 1
As we have a supercell there is a folding that will happen. The matrix element we want to compute is actually:
<n'k+q | dV/dR | nk> = <1L| H(1) | 1G> = <2G|H(1)|1G>
It can be shown (by making a taylor expansion of |psi_1(lambda)> and other quantities that :
<2G|H(1)|1G> = (e_1(lambda) - e_2(0)) <psi_2(0)|psi_1(1)> + O(lambda^2)
with indices indicating the band number and (0) an equilibrium position (unperturbed quantity) and with psi_1(1) = (psi_1(lambda)-psi_1(0))/lambda the first order perturbation obtained by finite difference also.
Using the following script (using the cut3d tool of abinit):
Code: Select all
# ----------------------------------------------------------------
# This python script is used to compute the matrix element
# <n,k+q|H(1)|n',k>
# Here we specfically take n=1, n'=2 at gamma
# This correspond to k=gamma and q=L
# ---------------------------------------------------------------
import os
import numpy as np
def read_eigenenergies(datafile):
flag = False
with open(datafile,'r') as f:
for line in f:
if flag:
eigenenergies = map(float,line.split())
break
if "kpt= 0.0000 0.0000 0.0000" in line:
flag = True
return eigenenergies
def use_cut3d(datafile,out,gamma,band):
if datafile[len(datafile)-3:len(datafile)] == "POT":
os.system("rm cut3d.files")
with open("cut3d.files","a") as files:
files.write(str(datafile)+"\n")
files.write("1 \n")
files.write("5 \n") #3D formatted data
#(output the bare 3D data - one column)
files.write(str(out)+"\n")
files.write("0 \n")
os.system("cut3d < cut3d.files ")
if datafile[len(datafile)-3:len(datafile)] == "WFK":
os.system("rm wfk.files")
with open("wfk.files","a") as files:
files.write(str(datafile)+"\n")
files.write("1 \n")
files.write("0 \n")
files.write(str(gamma)+" \n") # Gamma point
files.write(str(band)+" \n") # Band nb 6
files.write("0 \n")
files.write("0 \n")
files.write("1 \n") # real 3D data one column
files.write(str(out)+" \n")
files.write("0 \n")
os.system("cut3d < wfk.files ")
def first_order_wf(wfk0,wfk1,lam):
S=0.0; N=0.0
psi0 = np.loadtxt(wfk0).view(complex)
psi1 = np.loadtxt(wfk1).view(complex)
S = (psi1-psi0)/lam
# np.savetxt('first_order_wfc.txt', S)
return S
def integrationwf(wfk1,wfk2):
S=0.0; N=0.0
psi1 = np.loadtxt(wfk1).view(complex)
if type(wfk2) != str:
psi2 = wfk2
S = psi1.conj()*psi2
N1 = psi1.conj()*psi1
return S.sum()/N1.sum()
else:
psi2 = np.loadtxt(wfk2).view(complex)
S = psi1.conj()*psi2
N1 = psi1.conj()*psi1
return S.sum()/N1.sum()
# The base name for the supercell files are assumed to be named "supercell"
# and the primitive as assumed to be named "primitive"
# Value and number of the displacement
h = np.array([0.1,0.05,0.025,0.0125,0.00625,0.003125,0.0015625])
# Initialisation values
k = 0
# Frequency (computed in DFPT using ABINIT)
omega = np.array([0.00480084,0.005785943,0.002453274,0.005671174])
# Carbon mass in a.u.
Mc = 21894.16693
# --------------------------------------------------
# Compute <psi|dVscf|psi> = de_1G
# --------------------------------------------------
# Second derivative using finite centred differences
# The second dimension is the number of h in geometric progression
# h,h/2,h/4,h/8 etc..
# The third index correspond to modes (here 4)
# The first index correspond to the Richardson extrapolation
D = np.zeros([len(h),len(h),4])
# Prepare equilibrium files
# k=gamma
# q = L
use_cut3d("supercello_DS1_WFK","wfkEQ",1,2) # psi_2^{(0)}
use_cut3d("supercello_DS1_WFK","wfkEQ",1,1) # psi_2^{(0)}
energy2 = read_eigenenergies("supercello_DS1_EIG")
while k < len(h):
use_cut3d("supercello_DS"+str(k+2)+"_WFK","wfk",1,1) # psi_1(lambda)
psi1 = first_order_wf("wfkEQ_k1_b1_s1","wfk_k1_b1_s1",h[k])
print "psi1", psi1
LO = integrationwf("wfkEQ_k1_b2_s1",psi1)
print "LO:", LO
energy1 = read_eigenenergies("supercello_DS"+str(k+2)+"_EIG")
# Define the phase
# phase = <psi1(0)|psi1(lambda)>/| <psi1(0)|psi1(lambda)>|
phase_int = integrationwf("wfkEQ_k1_b1_s1","wfk_k1_b1_s1")
phase_mod = np.abs(phase_int)
phase = phase_int/phase_mod
# print "phase",phase
elph=(energy1[0]-energy2[1])*LO/phase
print "elph",elph
# g12k=elph*(np.sqrt(1/(2*Mc*omega[1])))
# print "g12k", g12k
# LA mode
D[0,k,0] = elph
k += 1
l = 0
while l < len(h):
k = 0
while k < len(h):
if (l-1) >= 0 and (k-1) >= 0:
# LA mode
D[l,k,0] = (D[l-1,k,0]-(D[l-1,k-1,0]*(1.0/(2**(2*l)))))/\
(1-(1.0/(2**(2*l))))
if D[l-1,k-1,0] == 0.0:
D[l,k,0] = 0.0
# LO mode
k += 1
l += 1
print "Richardson extrapolation of <psi_2^{(0)}|dVscf|psi_1^{(0)} for LO"
print DI obtain a value (after convergence of the finite difference using a Richardson extrapolation) of 0.05747 + 1E-12j
This value compare poorly with the one obtained by the 77_ddb/prt_gkk_yambo routine.
Am I missing something? Is there additional stuff in the yambo_elphon_gkk_bymode file? (not the yambo_elphon_gkksqtw_bymode)
Thank you.
Samuel.
PS: If you want to run the calculations, it take less than 10 min on a normal desktop computer.