[gpaw-users] gpaw-users Digest, Vol 175, Issue 17
Mikael Kuisma
kuisma at dtu.dk
Mon Sep 2 14:25:46 CEST 2024
Hi!
It is a bit hard to reproduce, since the script you provide doesn't import atoms, or define many of the variables. You are also assigning to stock_psi variable which is not defined. Could you provide the entire script, which runs. Where did you get this custom DOS code you are using?
BR,
Mikael
________________________________
From: gpaw-users <gpaw-users-bounces at listserv.fysik.dtu.dk> on behalf of DAMEGARIN SAMBIANI via gpaw-users <gpaw-users at listserv.fysik.dtu.dk>
Sent: Monday, September 2, 2024 12:00 PM
To: gpaw-users at listserv.fysik.dtu.dk <gpaw-users at listserv.fysik.dtu.dk>
Subject: Re: [gpaw-users] gpaw-users Digest, Vol 175, Issue 17
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 sys
from ase.transport.greenfunction import GreenFunction
from 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<mailto:gpaw-users-request at listserv.fysik.dtu.dk>> a écrit :
Send gpaw-users mailing list submissions to
gpaw-users at listserv.fysik.dtu.dk<mailto: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<mailto: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<mailto: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<mailto:kuisma at dtu.dk>>
To: DAMEGARIN SAMBIANI <sambianidamegarin at gmail.com<mailto:sambianidamegarin at gmail.com>>
Cc: "gpaw-users at listserv.fysik.dtu.dk<mailto:gpaw-users at listserv.fysik.dtu.dk>"
<gpaw-users at listserv.fysik.dtu.dk<mailto: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<mailto: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<mailto:sambianidamegarin at gmail.com>>
Sent: Tuesday, August 20, 2024 11:02 AM
To: Mikael Kuisma <kuisma at dtu.dk<mailto: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<mailto: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/0592aa66/attachment-0001.htm>
More information about the gpaw-users
mailing list