[gpaw-users] Request for clarification on certain functionalities of GPAW
Jens Jørgen Mortensen
jjmo at dtu.dk
Fri Sep 27 09:43:59 CEST 2024
Here is an example of how to extract the LCAO-hamiltonian in real space:
import numpy as np
from ase import Atoms
from gpaw import GPAW
from gpaw.lcao.tools import get_lcao_hamiltonian
def real_space(weight_k, kpt_kc, A_kMM, disp_c):
B_MM = np.zeros(A_kMM.shape[1:])
for w, kpt_c, A_MM in zip(weight_k, kpt_kc, A_kMM):
phase = np.exp(-2j * np.pi * kpt_c @ disp_c)
B_MM += w * (phase * A_MM).real
return B_MM
atoms = Atoms('H',
cell=[1.1, 5, 5],
pbc=[True, False, False])
atoms.center(axis=(1, 2))
atoms.calc = GPAW(
mode='lcao',
basis='dzp',
kpts=[9, 1, 1],
txt='h.txt')
atoms.get_potential_energy()
H_skMM, S_kMM = get_lcao_hamiltonian(atoms.calc)
weight_k = atoms.calc.get_k_point_weights()
kpt_kc = atoms.calc.get_ibz_k_points()
for d in range(-3, 4):
S_MM = real_space(weight_k, kpt_kc, S_kMM, [d, 0, 0])
H_MM = real_space(weight_k, kpt_kc, H_skMM[0], [d, 0, 0])
print(d, S_MM)
For this to work, you need to have enough k-points and you should not
use any symmeries to reduce the k-points (time reversal is OK).
Jens Jørgen
On 9/13/24 12:08, DAMEGARIN SAMBIANI via gpaw-users wrote:
> *
> *
>
> *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)
>
>
>
> _______________________________________________
> gpaw-users mailing list
> gpaw-users at listserv.fysik.dtu.dk
> https://listserv.fysik.dtu.dk/mailman/listinfo/gpaw-users
-------------- next part --------------
import numpy as np
from ase import Atoms
from gpaw import GPAW
from gpaw.lcao.tools import get_lcao_hamiltonian
def real_space(weight_k, kpt_kc, A_kMM, disp_c):
B_MM = np.zeros(A_kMM.shape[1:])
for w, kpt_c, A_MM in zip(weight_k, kpt_kc, A_kMM):
phase = np.exp(-2j * np.pi * kpt_c @ disp_c)
B_MM += w * (phase * A_MM).real
return B_MM
atoms = Atoms('H',
cell=[1.1, 5, 5],
pbc=[True, False, False])
atoms.center(axis=(1, 2))
atoms.calc = GPAW(
mode='lcao',
basis='dzp',
kpts=[9, 1, 1],
txt='h.txt')
atoms.get_potential_energy()
H_skMM, S_kMM = get_lcao_hamiltonian(atoms.calc)
weight_k = atoms.calc.get_k_point_weights()
kpt_kc = atoms.calc.get_ibz_k_points()
for d in range(-3, 4):
S_MM = real_space(weight_k, kpt_kc, S_kMM, [d, 0, 0])
H_MM = real_space(weight_k, kpt_kc, H_skMM[0], [d, 0, 0])
print(d, S_MM)
More information about the gpaw-users
mailing list