[gpaw-users] gpaw-users Digest, Vol 175, Issue 17
DAMEGARIN SAMBIANI
sambianidamegarin at gmail.com
Mon Sep 2 14:40:46 CEST 2024
Hi,
Thank you for your feedback.
The custom DOS code I’m using is developed by me. It employs the
wavefunction method to calculate the density of states using the following
formula:
*Im[ 1 / (2π) ∑**_{α**} ΨαE * ΨαE†]*
I hope this clarifies the origin of the code. Please let me know if you
need any additional information or the complete working script.
Best regards,
Le lun. 2 sept. 2024 à 12:34, DAMEGARIN SAMBIANI <
sambianidamegarin at gmail.com> a écrit :
>
>
> *Hi!*
>
> *Thank you for your feedback.*
>
> I apologize for the oversight. You will find the complete script below,
> including the necessary imports and variable definitions.
> #!/usr/bin/env python
> """Transport exersice
>
> This file should do the same as pt_h2_lcao.py, but extracts the
> Hamiltonians
> manually instead of using gpawtransport, which currently does not work
> """
> #black real.py pour indenter le code python3
> 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.transport.greenfunction import GreenFunction
> from gpaw.density import Density, RealSpaceDensity
>
> # from gpaw.densities import Densities
> 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.coulomb import get_vxc
> from gpaw.lcao.projected_wannier import get_bfs
> from gpaw.lcao.el_ph import ElectronPhononCouplingMatrix
> 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
>
>
>
>
> def fermi_dirac(E, Ef, KBT):
> x = (E - Ef) / KBT
> if x > 0:
> exp_term = np.exp(-x)
> return 1.0 / (1.0 + exp_term)
> else:
> exp_term = np.exp(x)
> return exp_term / (1.0 + exp_term)
>
> 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()
> H_skMM, S_kMM = get_lcao_hamiltonian(calc)
> e, dos = calc.get_dos(spin=0, npts=2001, width=0.2)
> 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()
> dos2 = calc.get_all_electron_density(gridrefinement=2)
> print(dos2.shape) #(176, 48, 48)
> #np.savetxt("ALL_electron_Density.txt", dos2)
>
> ################################################################################################################################
> # Realspace Hamiltonien
> Scattering region
>
> ################################################################################################################################
> 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)
>
>
>
> ################################################################################################################################
> # elph coupling
>
> ################################################################################################################################
> #from gpaw.elph.electronphonon import ElectronPhononCoupling
> from gpaw.mpi import world
> #from gpaw.elph import DisplacementRunner
> #elph = DisplacementRunner(atoms=atoms, calc=calc,
> # supercell=(3, 3, 3), name='elph',
> # calculate_forces=True)
> #supercell=(3, 3, 3)
> #elph = ElectronPhononCoupling(atoms, calc, supercell=supercell,
> # name='elph', calculate_forces=True)
> #elph.run()
> #q = [[0., 0., 0.], [1. / 11., 1. / 11., 1. / 11.]]
> #M = elph.get_M(modes, q)
>
>
> ################################################################################################################################
> # STRUCTURE DE BANDE
> DE 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)
> #-------------------
> # base = get_bfs(calc)
>
> ################################################################################################################################
> # Exchange correlation
> potentiel
>
> ################################################################################################################################
> """Vxc = get_vxc(paw, spin=0, U=None)"""
>
> # 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]
> fmt = ["%16.10E" for _ in range(H.shape[1])] #
> np.savetxt("H1.txt", H, fmt=fmt)
> fmt = ["%16.10E" for _ in range(S.shape[1])] #
> np.savetxt("S1.txt", S, fmt=fmt)
> #bzk_kc, weight_k, H_skMM, S_kMM = get_lead_lcao_hamiltonian(calc,
> direction="x")
> #ibzk_t_kc, weights_t_k, h_skii, s_kii = get_realspace_hs(H_skMM, S_kMM ,
> bzk_kc, weight_k)
> # Only use first kpt, spin, as there are
> pickle.dump((H, S), open("lead1_hs.pickle", "wb"), 2)
>
> ################################################################################################################################
> # Right principal layer
> #
>
> ################################################################################################################################
> pickle.dump((H, S), open("lead2_hs.pickle", "wb"), 2)
> """Transport"""
> # Read in the hamiltonians
> 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"))
>
> #tcalc = 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
> """Procedure of extract sigma matrice"""
> eta = 1e-4
> 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=1e-4)
> # Calcul de la DOS et PDOS à une énergie de fermi
> """Total density of states -1/pi Im(Tr(GS))"""
> dos_value = gf.dos(Ef)
>
> #exit(1)
> """Projected density of states -1/pi Im(SGS/S)"""
> pdos_values = gf.pdos(Ef)
> fmt = "%16.10E"
> np.savetxt("pdos", pdos_values, fmt=fmt)
>
>
> E_min = -60#15 # eV
> E_max = 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 = 0.65 #0.65 #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)
> data_matrix1[i, 0] = E
> data_matrix1[i, 1] = ldos_value1
>
> 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_value1:e}\n')
>
> np.savetxt("DOS1.txt", data_matrix1)
> np.savetxt("DOS2.txt", data_matrix2)
>
> Le lun. 2 sept. 2024 à 12:25, Mikael Kuisma <kuisma at dtu.dk> a écrit :
>
>> 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>
>> 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/85a5a4dd/attachment-0001.htm>
More information about the gpaw-users
mailing list