[gpaw-users] Request for clarification on certain functionalities of GPAW
DAMEGARIN SAMBIANI
sambianidamegarin at gmail.com
Fri Sep 13 12:08:45 CEST 2024
*Hi GPAW Developer,*
I am writing to you again regarding several questions I have asked about
the use of GPAW, for which I have not yet received a response. I understand
that you must be very busy, but I would greatly appreciate it if you could
take a moment to help me clarify these points.
*First, *I would like to know how to extract the Hamiltonian in real space
from an LCAO basis. I have proposed a method using the *TightBinding*
module with the *h_and_s* function, but I would like to get your opinion on
this approach or know if there is a more appropriate method.
*Second*, during my calculations of the density of states (DOS) using the
*get_dos* function, I observed different results compared to those obtained
with the GreenFunction module implemented in ASE(*gf = GreenFunction(h, s,
selfenergies, eta)*). I do not understand why there is this difference and
would be grateful for any clarification on this matter.
Thank you in advance for your valuable help, and I look forward to your
responses.
Please find my code below for reference.
Best regards,
SAMBIANI Damegarin
import math
from ase.transport.tools import subdiagonalize, cutcoupling
from gpaw.lcao.pwf2 import get_lcao_xc, get_xc2
from scipy.sparse.linalg import spsolve
import scipy.sparse as sp
import numpy as np
import sys
from ase.transport.greenfunction import GreenFunction
from ase.transport.selfenergy import LeadSelfEnergy, BoxProbe
from numpy import linalg
from gpaw.density import Density, RealSpaceDensity
import pickle
from ase.transport.calculators import TransportCalculator
import pickle as pickle
from gpaw.lcao.tools import (
remove_pbc,
get_lcao_hamiltonian,
get_lead_lcao_hamiltonian,
get_realspace_hs,
lead_kspace2realspace,
)
from gpaw.lcao.tightbinding import TightBinding
from gpaw import GPAW, Mixer, FermiDirac
from ase import Atoms
from scipy.linalg import solve
import matplotlib.pyplot as plt
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()
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()
################################################################################################################################
# Realspace Hamiltonien
Scattering region
################################################################################################################################
S1 = TightBinding(atoms, calc) #use module TightBinding to
extract hamiltonien in real space
H_NMM, S_NMM = S1.h_and_s() # (2,158,158) # hamiltonien in real space
H, S = H_NMM[0], S_NMM[0]
np.savetxt("H.txt", H)
np.savetxt("S.txt", S)
pickle.dump((H.astype(complex), S.astype(complex)), open("scat_hs.pickle",
"wb"), 2)
################################################################################################################################
# BANDE STRUCTURE OF
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)
H_NMM2, S_NMM2 = S2.h_and_s()
H, S = H_NMM2[0], S_NMM2[0]
np.savetxt("H1.txt", H)
np.savetxt("S1.txt", S)
pickle.dump((H, S), open("lead1_hs.pickle", "wb"), 2)
pickle.dump((H, S), open("lead2_hs.pickle", "wb"), 2)
################################################################################################################################
# Right principal layer #
################################################################################################################################
"""Transport"""
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"))
calc = TransportCalculator(
h=h,
h1=h1,
h2=h2, # hamiltonian matrices
s=s,
s1=s1,
s2=s2, # overlap matrices
align_bf=1,
)
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
align_bf = 1
diff = (h_mm[align_bf, align_bf] - h1_ii[align_bf, align_bf]) /
s_mm[align_bf, align_bf]
h_mm -= diff * s_mm
"""Procedure of extract sigma matrice"""
eta = 1e-5
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)
E_min = -1.5 #60 # 15
# eV
E_max = 1.5 #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 #V.femtseconde
psi_st = np.zeros((N, 1), dtype=np.complex128)
psi_x = np.zeros((N, 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)"""
ldos_values = []
print("===== FIN SIGMA =======")
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
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 gf 1
data_matrix1[i, 0] = E
data_matrix1[i, 1] = ldos_value1
ldos_value2 = DOS_ST.imag.trace() / np.pi
# DOS WF 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)
np.savetxt("DOS2.txt", data_matrix2)
np.savetxt("rho.txt", rho)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://listserv.fysik.dtu.dk/pipermail/gpaw-users/attachments/20240913/7c5aae16/attachment-0001.htm>
More information about the gpaw-users
mailing list