[gpaw-users] gpaw-users Digest, Vol 175, Issue 17

DAMEGARIN SAMBIANI sambianidamegarin at gmail.com
Mon Sep 2 12:00:36 CEST 2024


Dear GPAW Developers,

Hello,

I am reaching out to you because I have a concern regarding the calculation
of the density of states. GPAW offers several modules for this purpose,
which leaves me uncertain about which one to use.

In my case, I used the get_dos function to calculate the electronic density
of states. However, when I use the GreenFunction module, which also has a
DOS implementation, to perform the same calculation, I notice that the
resulting curves are different. This discrepancy is puzzling, and I am
unsure how to proceed.

My goal is to ensure that I am correctly extracting the Hamiltonians and
that, using the wavefunction method, I obtain the same density of states.

I have included a snippet of the code below for your reference.

Thank you in advance for your help.

*from gpaw.lcao.tightbinding import TightBinding*



*import sysfrom ase.transport.greenfunction import GreenFunctionfrom
ase.transport.selfenergy import LeadSelfEnergy, BoxProbe*

# 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()

#*Dos calculation*

e, dos = calc.get_dos(spin=0, npts=2001, width=0.2) *# calcul of DOS *
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)  #*Plot DOS*

*#EXTRACT REAL SPACE HAMILTONIEN WITH THE MODULE OF TIGHTBINDING from lcao*

S1 = TightBinding(atoms, calc)
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)

*#USE GREEN FUNCTION MODULE*

gf = GreenFunction(h, s, selfenergies, eta=1e-4)


E_min     = -15                                                    # 0.56
 # eV
E_max     =  10                                                     # 1.44
 # eV
Ne        = 150
dE        = (E_max - E_min) / Ne
N         = 158
nmode     = 158
co        = 0 + 0j
Ci        = 1j

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
        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

       * #First calculation of DOS*

        ldos_value1 = gf.dos(E)
        data_matrix1[i, 0] = E
        data_matrix1[i, 1] = ldos_value2

        *#second calculation of Dos*

        ldos_value2 = DOS_ST.imag.trace() / np.pi
        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_value:e}\n')

np.savetxt("DOS1.txt", data_matrix1)*# SAVE DOS 1*
np.savetxt("DOS2.txt", data_matrix2)#* SAVE DOS 2*

Le mer. 21 août 2024 à 10:27, <gpaw-users-request at listserv.fysik.dtu.dk> a
écrit :

> Send gpaw-users mailing list submissions to
>         gpaw-users at listserv.fysik.dtu.dk
>
> To subscribe or unsubscribe via the World Wide Web, visit
>         https://listserv.fysik.dtu.dk/mailman/listinfo/gpaw-users
> or, via email, send a message with subject or body 'help' to
>         gpaw-users-request at listserv.fysik.dtu.dk
>
> You can reach the person managing the list at
>         gpaw-users-owner at listserv.fysik.dtu.dk
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of gpaw-users digest..."
>
>
> Today's Topics:
>
>    1. Re:  Request for Assistance with Band Structure of H? Using
>       GPAW (Mikael Kuisma)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Tue, 20 Aug 2024 13:21:33 +0000
> From: Mikael Kuisma <kuisma at dtu.dk>
> To: DAMEGARIN SAMBIANI <sambianidamegarin at gmail.com>
> Cc: "gpaw-users at listserv.fysik.dtu.dk"
>         <gpaw-users at listserv.fysik.dtu.dk>
> Subject: Re: [gpaw-users]  Request for Assistance with Band Structure
>         of H? Using GPAW
> Message-ID:
>         <
> AS8P192MB189592471C371AF88A07294DB68D2 at AS8P192MB1895.EURP192.PROD.OUTLOOK.COM
> >
>
> Content-Type: text/plain; charset="utf-8"
>
> Hi!
>
> In the first part, where you plot the DOS, you have already plotted the
> closest thing that a molecule has to a band structure: density of states.
> Whatever you try to do next is ill-defined, as molecules are localized
> particles, thus non-periodic, and thus do not have any band structure. If
> you were to calculate the band structure of a molecule in a large enough
> cell, you would get that each band is completely flat, and it would anyway
> reduce to all of the information available in the density of states plot.
>
> I suggest to read basics of band structure theory, and then start from the
> GPAW tutorial:
>
> https://wiki.fysik.dtu.dk/gpaw/tutorialsexercises/electronic/band_structure/bands.html
>
> https://wiki.fysik.dtu.dk/gpaw/tutorialsexercises/electronic/bandstructures/bandstructures.html#bandstructures
>
> BR.
> Mikael
>
> ________________________________
> From: DAMEGARIN SAMBIANI <sambianidamegarin at gmail.com>
> Sent: Tuesday, August 20, 2024 11:02 AM
> To: Mikael Kuisma <kuisma at dtu.dk>
> Subject: Request for Assistance with Band Structure of H? Using GPAW
>
>
> Dear GPAW Developers,
>
> I am writing to request your assistance in plotting the band structure of
> the H? molecule using GPAW. Despite my efforts, I am still facing
> difficulties in obtaining a correct plot of the band structure.
>
> Attached, you will find the code I used for this simulation. I would
> greatly appreciate any guidance you can provide to help resolve this issue.
>
> Thank you in advance for your assistance and expertise.
>
> from ase.build import molecule
> from gpaw import GPAW, PW, FermiDirac
> import matplotlib.pyplot as plt
>
> # Cr?ez la mol?cule H2
> h2 = molecule('H2')
> h2.center(vacuum=3.0)
>
> # Calcul de l'?tat fondamental
> calc = GPAW(mode=PW(200),
>             xc='PBE',
>             kpts=(1, 1, 1),
>             random=True,
>             occupations=FermiDirac(0.01),
>             txt='H2_gs.txt')
> h2.calc = calc
> h2.get_potential_energy()
> ef = calc.get_fermi_level()
> print(ef)
> calc.write('H2_gs.gpw')
> e, dos = calc.get_dos()
> plt.plot(e, dos)
> plt.ylabel('DOS')
> plt.show()
>
> # Bande
> from ase import Atoms
> from ase.build import molecule
> from ase.visualize import view
> from gpaw import GPAW
> from ase.dft.kpoints import bandpath
> import matplotlib.pyplot as plt
> from gpaw.spinorbit import soc_eigenstates
>
> calc = GPAW('H2_gs.gpw').fixed_density(
>     nbands=158,
>     symmetry='off',
>     kpts={'path': 'GX', 'npoints': 60},
>     convergence={'bands': 8})
> bs = calc.band_structure()
> bs.plot(filename='bandH2.png', show=True, emax=10.0)
>
> Sincerely,
>
> SAMBIANI
>
> PhD student
> -------------- next part --------------
> An HTML attachment was scrubbed...
> URL: <
> http://listserv.fysik.dtu.dk/pipermail/gpaw-users/attachments/20240820/7442d19a/attachment-0001.htm
> >
>
> ------------------------------
>
> Subject: Digest Footer
>
> _______________________________________________
> gpaw-users mailing list
> gpaw-users at listserv.fysik.dtu.dk
> https://listserv.fysik.dtu.dk/mailman/listinfo/gpaw-users
>
>
> ------------------------------
>
> End of gpaw-users Digest, Vol 175, Issue 17
> *******************************************
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://listserv.fysik.dtu.dk/pipermail/gpaw-users/attachments/20240902/2253ec26/attachment-0001.htm>


More information about the gpaw-users mailing list