[gpaw-users] gpaw-users Digest, Vol 175, Issue 17
DAMEGARIN SAMBIANI
sambianidamegarin at gmail.com
Mon Sep 2 14:34:06 CEST 2024
*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/4cb1ea8a/attachment-0001.htm>
More information about the gpaw-users
mailing list