[gpaw-users] Request for clarification on certain functionalities of GPAW
DAMEGARIN SAMBIANI
sambianidamegarin at gmail.com
Mon Sep 30 18:23:51 CEST 2024
Dear GPAW developer,
I would like to seek your assistance regarding the use of the GPAW Poisson
solver with a 2D density. I am currently working on a project that requires
this approach, and I would like to ensure that I am proceeding correctly.
Here is a proposal of how I plan to use the Poisson solver with a 2D
density:
#@pytest.mark.skip(reason='TODO')
def test_generalizedlaue():
from gpaw.poisson import FDPoissonSolver
from gpaw.poisson import PoissonSolver, FFTPoissonSolver
from gpaw.grid_descriptor import GridDescriptor
from gpaw.solvation.poisson import SolvationPoissonSolver,
WeightedFDPoissonSolver,PolarizationPoissonSolver
import numpy as np
ny, nz = 158, 158
nx = 158
gd = GridDescriptor((nx, ny, nz), (nx, ny, nz), pbc_c=(True, True,
True))
#1
poisson1 = FDPoissonSolver(nn=3, relax='J', eps=2e-10, maxiter=1000)
poisson1.set_grid_descriptor(
gd)
#2
#poisson2_solver= PoissonSolver(nn=3, eps=2e-10)
#poisson2_solver.set_grid_descriptor(gd)
#3
#poisson3_solver= FFTPoissonSolver()
#poisson3_solver.set_grid_descriptor(gd)
phi2_g = gd.zeros()
phi2_g[-1, 0, :] = 3 # Première ligne
phi2_g[-1, -1, :] = 3 # Dernière ligne
phi2_g[-1, :, 0] = 3 # Première colonne
phi2_g[-1, :, -1] = 3 # Dernière colonne
#print(phi2_g)
rho_3d = gd.zeros()
deltarho = np.ones((158,158))
rho_3d[-1, :, :] = deltarho
poisson1.solve(phi2_g, rho_3d)
print(phi2_g.shape)
print(phi2_g)
test_generalizedlaue()
SAMBIANI
PhD student at University of Lomé
Le ven. 27 sept. 2024 à 07:44, Jens Jørgen Mortensen <jjmo at dtu.dk> a écrit :
> 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 --------------
An HTML attachment was scrubbed...
URL: <http://listserv.fysik.dtu.dk/pipermail/gpaw-users/attachments/20240930/e1986e52/attachment-0001.htm>
More information about the gpaw-users
mailing list