[gpaw-users] Request for Assistance with density of state calculation
DAMEGARIN SAMBIANI
sambianidamegarin at gmail.com
Tue Sep 3 13:31:30 CEST 2024
Dear GPAW Developer,
I am reaching out again because I have not yet received a response to my
previous question. My query was about the difference in results when
calculating the density of states with GPAW using the *get_dos *function
compared to using the density of states from the *GreenFunction *class.
I also tried to implement the density of states using wave functions with
the following relation:
*\[\text{Im}\left[ \frac{1}{2\pi} \sum_{\alpha} \Psi_{\alpha E}(i, t)
\Psi_{\alpha E}^{\dagger}(i, t) \right]\]*
The results obtained seem to be closer to those provided by the *GreenFunction
*module. How can this be explained? You will find the complete code below.
Thank you in advance for your help and clarification.
#!/usr/bin/env python
"""Transport exersice
This file should do the same as pt_h2_lcao.py, but extracts the Hamiltonians
manually instead of using gpawtransport, which currently does not work
"""
#black real.py pour indenter le code python3
*import mathfrom ase.transport.tools import subdiagonalize, cutcouplingfrom
gpaw.lcao.pwf2 import get_lcao_xc, get_xc2# from diff_matrices import
Diff_mat_1D, Diff_mat_2Dfrom scipy.sparse.linalg import spsolveimport
scipy.sparse as spimport numpy as npimport sysfrom
ase.transport.greenfunction import GreenFunctionfrom
ase.transport.selfenergy import LeadSelfEnergy, BoxProbefrom numpy import
linalgfrom gpaw.fftw import check_fft_size ,create_plans,FFTPlans, FFTPlan,
FFTWPlan, FFTWPlans, NumpyFFTPlan, NumpyFFTPlans# from
gpaw.transport.greenfunction import GreenFunctionfrom gpaw.density import
Density, RealSpaceDensity# from gpaw.densities i*mport Densities
*import picklefrom ase.transport.calculators import
TransportCalculatorimport pickle as picklefrom gpaw.lcao.tools import (
remove_pbc, get_lcao_hamiltonian, get_lead_lcao_hamiltonian,
get_realspace_hs, lead_kspace2realspace,)from gpaw.coulomb import
get_vxcfrom gpaw.lcao.projected_wannier import get_bfsfrom gpaw.lcao.el_ph
import ElectronPhononCouplingMatrixfrom gpaw.lcao.tightbinding import
TightBindingfrom gpaw import GPAW, Mixer, FermiDiracfrom ase import
Atomsfrom scipy.linalg import solveimport matplotlib.pyplot as plt*
def fermi_dirac(E, Ef, KBT):
x = (E - Ef) / KBT
if x > 0:
exp_term = np.exp(-x)
return 1.0 / (1.0 + exp_term)
else:
exp_term = np.exp(x)
return exp_term / (1.0 + exp_term)
a = 2.41 # Pt binding length
b = 0.90 # H2 binding length
c = 1.70 # Pt-H binding length
L = 7.00 # width of unit cell
###################################################################################################################################
# *Scattering region #*
###################################################################################################################################
# Setup the Atoms for the scattering region.
atoms = Atoms("Pt5H2Pt5", pbc=(1, 0, 0), cell=[9 * a + b + 2 * c, L, L])
atoms.positions[:5, 0] = [i * a for i in range(5)]
atoms.positions[-5:, 0] = [i * a + b + 2 * c for i in range(4, 9)]
atoms.positions[5:7, 0] = [4 * a + c, 4 * a + c + b]
atoms.positions[:, 1:] = L / 2.0
# Attach a GPAW calculator
calc = GPAW(
h=0.3,
xc="PBE",
basis="szp(dzp)",
occupations=FermiDirac(width=0.1),
kpts=(2, 1, 1),
mode="lcao",
txt="pt_h2_lcao_scat.txt",
mixer=Mixer(0.1, 5, weight=100.0),
symmetry={"point_group": False, "time_reversal": False},
)
atoms.calc = calc
atoms.get_potential_energy() # Converge everything!
Ef = atoms.calc.get_fermi_level()
H_skMM, S_kMM = get_lcao_hamiltonian(calc)
e, dos = calc.get_dos(spin=0, npts=2001, width=0.2) #* DOS CALCULATION *
data = np.column_stack((e, dos))
np.savetxt("DOS.txt", data, fmt='%e', header='Energy (eV) DOS')
#np.savetxt("DOS.txt", e, dos)
plt.plot(e, dos)
plt.ylabel("DOS")
plt.savefig('DOS_Pt5-H2-Pt5.png', dpi=300)
#plt.show()
dos2 = calc.get_all_electron_density(gridrefinement=2)
print(dos2.shape) #(176, 48, 48)
#np.savetxt("ALL_electron_Density.txt", dos2)
################################################################################################################################
# Realspace Hamiltonien
Scattering region
################################################################################################################################
S1 = TightBinding(atoms, calc)*#USE TightBinding to extract hamiltonien in
real space*
H_NMM, S_NMM = S1.h_and_s() # (2,158,158)
H, S = H_NMM[0], S_NMM[0]
fmt = ["%16.10E" for _ in range(H.shape[1])]
np.savetxt("H.txt", H, fmt=fmt)
fmt = ["%16.10E" for _ in range(S.shape[1])]
np.savetxt("S.txt", H, fmt=fmt)
pickle.dump((H.astype(complex), S.astype(complex)), open("scat_hs.pickle",
"wb"), 2)
################################################################################################################################
# *STRUCTURE DE BANDE
DE Pt5-H2-Pt5 Scattering region*
################################################################################################################################
from ase.dft.kpoints import bandpath
calc.write('pt_h2_gs.gpw')
calc = GPAW('pt_h2_gs.gpw').fixed_density(
nbands=158,
symmetry='off',
kpts={'path': 'GXG', 'npoints': 60},
convergence={'bands': 8})
bs = calc.band_structure()
bs.plot(filename='bandP5H2Pt51.png', show=True,emin =-10)#, emax=30.0)
#-------------------
# Use four Pt atoms in the lead, so only take those from before
################################################################################################################################
# *Left principal layer #*
#################################################################################################################################
atoms = atoms[:4].copy()
atoms.set_cell([4 * a, L, L])
# Attach a GPAW calculator
calc = GPAW(
h=0.3,
xc="PBE",
basis="szp(dzp)",
occupations=FermiDirac(width=0.1),
kpts=(4, 1, 1), # More kpts needed as the x-direction is shorter
mode="lcao",
txt="pt_h2_lcao_llead.txt",
mixer=Mixer(0.1, 5, weight=100.0),
symmetry={"point_group": False, "time_reversal": True},
)
atoms.calc = calc
atoms.get_potential_energy() # Converge everything!
Ef = atoms.calc.get_fermi_level()
################################################################################################################################
# *Realspace Hamiltonien
Left principal layer*#
################################################################################################################################
S2 = TightBinding(atoms, calc) *#use TightBinding to extract hamiltonien in
real space*
H_NMM2, S_NMM2 = S2.h_and_s()
H, S = H_NMM2[0], S_NMM2[0]
fmt = ["%16.10E" for _ in range(H.shape[1])] #
np.savetxt("H1.txt", H, fmt=fmt)
fmt = ["%16.10E" for _ in range(S.shape[1])] #
np.savetxt("S1.txt", S, fmt=fmt)
#bzk_kc, weight_k, H_skMM, S_kMM = get_lead_lcao_hamiltonian(calc,
direction="x")
#ibzk_t_kc, weights_t_k, h_skii, s_kii = get_realspace_hs(H_skMM, S_kMM ,
bzk_kc, weight_k)
# Only use first kpt, spin, as there are
pickle.dump((H, S), open("lead1_hs.pickle", "wb"), 2)
################################################################################################################################
# * Right principal layer*
#
################################################################################################################################
pickle.dump((H, S), open("lead2_hs.pickle", "wb"), 2)
"""Transport"""
# Read in the hamiltonians
h, s = pickle.load(open("scat_hs.pickle", "rb"))
h1, s1 = pickle.load(open("lead1_hs.pickle", "rb"))
h2, s2 = pickle.load(open("lead2_hs.pickle", "rb"))
h_mm = h
s_mm = s
pl1 = len(h1) // 2
pl2 = len(h2) // 2
h1_ii = h1[:pl1, :pl1]
h1_ij = h1[:pl1, pl1 : 2 * pl1]
s1_ii = s1[:pl1, :pl1]
s1_ij = s1[:pl1, pl1 : 2 * pl1]
h2_ii = h2[:pl2, :pl2]
h2_ij = h2[pl2 : 2 * pl2, :pl2]
s2_ii = s2[:pl2, :pl2]
s2_ij = s2[pl2 : 2 * pl2, :pl2]
# if hc1 is None:
nbf = len(h_mm)
h1_im = np.zeros((pl1, nbf), complex)
s1_im = np.zeros((pl1, nbf), complex)
h1_im[:pl1, :pl1] = h1_ij
s1_im[:pl1, :pl1] = s1_ij
hc1 = h1_im
sc1 = s1_im
# if hc2 is None
h2_im = np.zeros((pl2, nbf), complex)
s2_im = np.zeros((pl2, nbf), complex)
h2_im[-pl2:, -pl2:] = h2_ij
s2_im[-pl2:, -pl2:] = s2_ij
hc2 = h2_im
sc2 = s2_im
"""Procedure of extract sigma matrice"""
eta = 1e-4
eta1 = 1e-5
eta2 = 1e-5
*selfenergies = [ LeadSelfEnergy((h1_ii, s1_ii), (h1_ij, s1_ij), (h1_im,
s1_im), eta1), LeadSelfEnergy((h2_ii, s2_ii), (h2_ij, s2_ij), (h2_im,
s2_im), eta2),]L1 = LeadSelfEnergy((h1_ii, s1_ii), (h1_ij, s1_ij), (h1_im,
s1_im), eta1)L2 = LeadSelfEnergy((h2_ii, s2_ii), (h2_ij, s2_ij), (h2_im,
s2_im), eta2)gf = GreenFunction(h, s, selfenergies, eta=1e-4)*
# Calcul de la DOS et PDOS à une énergie de fermi
"""Total density of states -1/pi Im(Tr(GS))"""
dos_value = gf.dos(Ef)
"""Projected density of states -1/pi Im(SGS/S)"""
pdos_values = gf.pdos(Ef)
fmt = "%16.10E"
np.savetxt("pdos", pdos_values, fmt=fmt)
*# transport part*
E_min = -60#15 # eV
E_max = 10 # eV
Ne = 150
dE = (E_max - E_min) / Ne
N = 158
nmode = 158
co = 0 + 0j
Ci = 1j
e_charge = 1.602176634e-19
KBT = 8.6173303e-4 #
Pour T=10K
hbar = 6.582119570e-16 # eV.s
psi_st = np.zeros((N, 1), dtype=np.complex128)
psi_x = np.zeros((N, 1), dtype=np.complex128)
tmp_R = np.zeros((N, 1, 1))
tmp_Q = np.zeros((N, 1, 1), dtype=np.complex128)
eigval = np.zeros((N, 1), dtype=np.complex128)
eigvec1 = np.zeros((N, N), dtype=np.complex128)
eigvec2 = np.zeros((N, N), dtype=np.complex128)
stock_psi = np.zeros((Ne , N, N), dtype=np.complex128)
rho = np.zeros((N, N), dtype=np.complex128)
rho_i = np.zeros((N, N), dtype=np.complex128)
sigma_mm1 = np.zeros((N, N), dtype=np.complex128)
sigma_mm2 = np.zeros((N, N), dtype=np.complex128)
sigma_tot = np.zeros((N, N), dtype=np.complex128)
gamma = np.zeros((N, N), dtype=np.complex128)
data_matrix1 = np.zeros((Ne, 2))
data_matrix2 = np.zeros((Ne, 2))
"""print('===== SIGMA =======')
sigma_mm1 = L1.retarded(Ef)
sigma_mm2 = L2.retarded(Ef)
sigma_tot = sigma_mm1 + sigma_mm2
fmt = ["(%16.10E,%16.10E)" for _ in range(sigma_tot.shape[1])]
np.savetxt("sigmatot",sigma_tot, fmt=fmt)
gamma = 1.0j * (sigma_tot - sigma_tot.T.conj())
eigval, eigvec1 = np.linalg.eig(gamma)
print('===== FIN SIGMA =======')"""
ldos_values = [ ]
print("================= Début du calcul stationnaire
=======================")
with open("DOS_st.txt", 'a') as file:
if file.tell() == 0:
file.write('#Energy (eV) DOS\n')
for i in range(Ne):
E = E_min + (i - 1) * dE
# eV
sigma_mm1 = L1.retarded(E)
sigma_mm2 = L2.retarded(E)
sigma_tot = sigma_mm1 + sigma_mm2
gamma = 1.0j * (sigma_tot - sigma_tot.T.conj())
eigval, eigvec1 = np.linalg.eig(gamma)
GR = E * s_mm - h_mm - sigma_tot
#GR = solve((E + 1j * eta) * s_mm - h_mm - sigma_tot,
np.identity(N))
GA = np.transpose(np.conj(GR))
# eV
rho_i+= dE * (GR @ sigma_tot @ GA) / np.pi
rho_mode = np.zeros((N, N), dtype=np.complex128)
DOS_ST = np.zeros((N, N), dtype=np.complex128)
for j in range(1, nmode):
vec1 = eigvec1[:, j].reshape((nmode, 1))
vec = vec1 * np.sqrt(eigval[j])
# eV
psi_st = solve(GR, vec)
psi_x = psi_st #np.fft.ifft(psi_st) * nmode
#plans = create_plans([N, 1, 1], float)
#tmp_Q[:, 0, 0] = psi_st.reshape((N))
#plans.tmp_Q[:, 0, 0]
#print(plans.tmp_R[:, 0, 0])
#plans.ifft()
#plans.tmp_R[:, 0, 0]
#print(plans.tmp_R[:, 0, 0])
#stock_psi[i, j, :] = plans.tmp_R[:, 0, 0].reshape((N))
#psi_x.reshape((N))
stock_psi[i, j, :] = psi_x.reshape((N))
rho_mode += (
Ci
* dE
# eV
* (1.0 / (2.0 * np.pi))
* fermi_dirac(E, Ef, KBT)
* np.matmul(psi_x, np.transpose(np.conj(psi_x)))
)
DOS_ST += Ci * dE * np.matmul(psi_x,
np.transpose(np.conj(psi_x)))
rho += rho_mode
ldos_value1 = gf.dos(E) *# DOS calculation 1*
data_matrix1[i, 0] = E
data_matrix1[i, 1] = ldos_value1
ldos_value2 = DOS_ST.imag.trace() / np.pi *# DOS calculation 2*
data_matrix2[i, 0] = E
data_matrix2[i, 1] = ldos_value2
# Sauvegarder l'énergie et le LDOS pour cette itération
file.write(f'{E:e} {ldos_value1:e}\n')
np.savetxt("DOS1.txt", data_matrix1) #* DOS 1 *
np.savetxt("DOS2.txt", data_matrix2) # *DOS 2*
The custom DOS code I'm using is developed by me. It employs the
wavefunction method to calculate the density of states using the following
formula:
*\[ \text{Im}\left[ \frac{1}{2\pi} \sum_{\alpha} \Psi_{\alpha E}
\Psi_{\alpha E}^{\dagger} \right] \]*
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://listserv.fysik.dtu.dk/pipermail/gpaw-users/attachments/20240903/8c43f723/attachment-0001.htm>
More information about the gpaw-users
mailing list