This tutorial targets beginners to intermediate users of SIESTA. It covers basic runs, k-points convergence, bulk and slab geometry optimization, surface reconstructions, density of states, and workfucntion. A second part deals with nudged elestic band (NEB) calculations. The tutorial utilizes the atomistic simulation enviroment (ASE) python wrapper for SIESTA. All cells are supposed to run separately, but in sequence, there are no global variables. The tutorial requires SIESTA 5.2.4 and ASE (3.29 not strict). A few python packages should be installed like matplotlib, pymatgen, etc. You will need a proper instalation and linking between ASE and SIESTA, described in the ASE manuals. You will need to download and adjust the path to the pseudopotential files, the tutorial used pseudodojo SIESTA PBE pseudopotentials. Basic understanding of python is required. Once you go through it, you will enjoy coding computational materials science!
About the author: Dr. Aleksandar Staykov, associate professor at Kyushu Universtity, Japan. E-mail: alex@i2cner.kyushu-u.ac.jp
for more information on theory, read my book which covers DFT fundamentals, wavefunctions, methods and applicatoins:
Ab Initio Simulations in Materials Science: Hands-On Introduction to Electronic Structure Modeling with VASP
https://www.amazon.com/Initio-Simulations-Materials-Science-Hands/dp/B0FGSNL5QR
SIESTA
SIESTA is a linear combination of atomic orbitals (LCAO) Density Functional Theory (DFT) code that employs norm‑conserving pseudopotentials and highly efficient numerical atomic basis sets. In this tutorial, all calculations were carried out using SIESTA 5.4.2 together with PseudoDojo-curated pseudopotentials to ensure consistency and high transferability.
For the electronic structure calculations, we used the standard double‑zeta polarized (DZP) basis set distributed with SIESTA, which provides a good balance between accuracy and computational cost for typical materials simulations.
The Atomistic Simulation Environment (ASE) served as the workflow manager, handling input generation, job execution, and post‑processing. Instead of relying on SIESTA’s built‑in geometry optimization routines, all structural relaxations were performed using ASE’s optimizers, which offer more flexibility and tighter integration with Python-based workflows.
Single‑point calculation performed with SIESTA through ASE introduces the essential structure of a SIESTA workflow: defining the atomic system, selecting the basis set and pseudopotentials, configuring the calculator, and executing the run to obtain the electronic ground‑state properties without performing any geometry optimization.
from ase import Atoms
from ase.build import bulk
from ase.calculators.siesta import Siesta
from ase.units import Ry
from ase.visualize import view
import os
#create a new directory for the calculation
os.mkdir("scf_siesta")
os.chdir("scf_siesta")
# use ASE to create a diamond structure of carbon
atoms = bulk('C', 'diamond', a=3.567, cubic=True)
# Minimal Siesta calculator
atoms.calc = Siesta(label='diamond', # system label
xc='PBE', # exchange-correlation functional
pseudo_path='/Users/staykov/pseudodojo', #path to pseudopotentials
pseudo_qualifier='gga', #if necessary, specify the qualifier for the pseudopotentials
symlink_pseudos=True,
mesh_cutoff=200 * Ry, # cutoff for the real-space grid
energy_shift=0.01 * Ry,
basis_set='DZP', # basis set
spin='non-polarized', # spin configuration
kpts=(4, 4, 4), # k-point grid
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 100}, # additional arguments for the SIESTA input file
)
energy = atoms.get_total_energy() # execute the calculation and get the total energy
os.chdir("../")
print("Energy:", energy)
view(atoms, viewer='x3d') # atomistic visualization of the structure using x3d viewer in IPython
Energy: -1308.986275
Job completed
Converging the k‑point mesh in SIESTA is a critical prerequisite for obtaining reliable and reproducible results. The required k‑point density depends on the Bravais lattice of the system and scales with the reciprocal‑space dimensions of the periodic simulation cell. In practice, this means that larger real‑space cells require fewer k‑points, while smaller cells demand a denser mesh. The appropriate k‑point grid must be determined specifically for the system under study. It cannot—and should not—be copied from another calculation, even if the material appears similar. Each structure has its own reciprocal‑space characteristics, and the optimal grid is therefore unique.A standard way to determine the correct mesh is to converge the total energy with respect to the k‑point sampling. By systematically increasing the grid density and monitoring how the total energy changes, one can identify the smallest mesh that yields stable, well‑converged results. This procedure ensures that all subsequent calculations—geometry optimizations, electronic structure analysis, or property evaluations—are built on a solid numerical foundation.
from ase import Atoms
from ase.calculators.siesta import Siesta
import os
from ase.build import bulk
from ase.units import Ry
import numpy as np
import matplotlib.pyplot as plt
os.mkdir("siesta_kpts")
os.chdir("siesta_kpts")
atoms = bulk('C', 'diamond', a=3.567, cubic=True)
kk_array = []
energy_array = []
for kk in range(1, 9, 1): # loop over k-points from 1x1x1 to 8x8x8
atoms.calc = Siesta(label='diamond',
xc='PBE',
pseudo_path='/Users/staykov/pseudodojo',
pseudo_qualifier='gga',
symlink_pseudos=True,
mesh_cutoff=200 * Ry,
energy_shift=0.01 * Ry,
basis_set='DZP',
spin='non-polarized',
kpts=(kk, kk, kk),
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 100},
)
energy = atoms.get_total_energy()
print("k-points:", kk)
kk_array.append(kk)
print("Energy:", energy)
energy_array.append(energy)
# Plot k-point energy convergence
plt.figure(figsize=(10, 6)) # Set the figure size (optional)
plt.plot(kk_array, energy_array, marker='o', linestyle='-', color='b') # Plot the data
plt.title('Energy vs k-points') # Add a title to the plot
plt.xlabel('k-points') # Label for the x-axis
plt.ylabel('Energy') # Label for the y-axis
plt.grid(True) # Add a grid (optional)
plt.show() # Display the plot
os.chdir("../")
Gamma-point calculation with interaction between periodic images Some features might not work optimally: e.g. DM initialization from atomic data Job completed
k-points: 1 Energy: -1293.182263
Job completed
k-points: 2 Energy: -1308.23107
Job completed
k-points: 3 Energy: -1308.933258
Job completed
k-points: 4 Energy: -1308.986276
Job completed
k-points: 5 Energy: -1308.994332
Job completed
k-points: 6 Energy: -1308.995062
Job completed
k-points: 7 Energy: -1308.995312 k-points: 8 Energy: -1308.995333
Job completed
Geometry and lattice optimization in this tutorial are performed using the ASE optimization algorithms, while SIESTA is used strictly as a force and energy calculator. This separation of roles is intentional. Many DFT packages in materials science—including SIESTA—provide basic geometry optimizers that are functional but often limited in robustness, flexibility, and convergence behavior. In contrast, external frameworks such as ASE offer modern, well‑tested optimization algorithms that integrate smoothly with Python workflows and provide more reliable convergence for both atomic positions and cell degrees of freedom. The optimizer used here is the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. BFGS is a quasi‑Newton method that combines gradient information with an evolving approximation of the inverse Hessian. This allows it to accelerate convergence by effectively estimating curvature without performing expensive second‑derivative calculations. For comparison, SIESTA’s built‑in Conjugate Gradient (CG) optimizer relies solely on first‑order information. CG constructs search directions that are conjugate with respect to the Hessian but never stores or approximates the Hessian itself. While efficient, CG typically converges more slowly and less smoothly than quasi‑Newton methods, especially for complex materials or systems with soft modes. To enable full lattice optimization, ASE provides the UnitCellFilter class, which exposes the cell degrees of freedom to the optimizer. When wrapped around the atomic configuration, UnitCellFilter ensures that both atomic positions and lattice vectors are updated consistently during the optimization process.
from ase import Atoms
from ase.calculators.siesta import Siesta
from ase.units import Ry
from ase.build import bulk
from ase.visualize import view
import os
from ase.optimize.bfgs import BFGS
from ase.filters import UnitCellFilter
from ase.io import write
os.mkdir("siesta_opt")
os.chdir("siesta_opt")
atoms = bulk('C', 'diamond', a=3.567, cubic=True)
print("Initial cell:")
print(atoms.cell) # print the initial cell parameters before the optimization
# Minimal Siesta calculator for testing
atoms.calc = Siesta(label='diamond',
xc='PBE',
pseudo_path='/Users/staykov/pseudodojo',
pseudo_qualifier='gga',
symlink_pseudos=True,
mesh_cutoff=230 * Ry,
energy_shift=0.01 * Ry,
basis_set='DZP',
spin='non-polarized',
kpts=(5, 5, 5),
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 100},
)
energy = atoms.get_total_energy()
ucf = UnitCellFilter(atoms) # set up a unit cell filter to allow optimization of the cell parameters
opt = BFGS(ucf, trajectory='cellopt.traj') # set up a BFGS optimizer to optimize the cell parameters of the structure
opt.run(fmax=0.01) # converge the optimization until the maximum force is below 0.01 eV/Å
# Run a single-point calculation
os.chdir("../")
print("Final cell:")
print(atoms.cell) # print the optimized cell parameters after the optimization is complete
write('cell_siesta.traj', atoms) # the optimization trajectory is saved as a set of geometry snapshots in a .traj file, which can be visualized later using ASE or other visualization tools
view(atoms, viewer='x3d')
Initial cell: Cell([3.567, 3.567, 3.567])
Job completed
Step Time Energy fmax BFGS: 0 00:44:32 -1308.994797 0.337870
Job completed
BFGS: 1 00:44:40 -1308.999467 0.307149
Job completed
BFGS: 2 00:44:51 -1309.022572 0.015012 BFGS: 3 00:45:00 -1309.022615 0.001113 Final cell: Cell([[3.591774264265943, 2.714633142111684e-20, 9.566565489278782e-21], [-3.6111392130239047e-20, 3.591774264265943, -2.8717733890076676e-20], [4.7555682922468305e-20, -2.9035486792415458e-21, 3.591774264265943]])
Job completed
The Density of States (DOS) is one of the fundamental electronic properties computed in Density Functional Theory. It describes how many electronic states are available at each energy level and provides direct insight into the material’s electronic behavior. By examining the DOS, we can identify the valence band, conduction band, and the band gap, which together determine whether a material behaves as an insulator, semiconductor, or metal. Because of this, DOS analysis is an essential step in interpreting and understanding the electronic structure obtained from DFT calculations.
from ase import Atoms
from ase.calculators.siesta import Siesta
from ase.units import Ry
import os
from ase.io import read
from ase.dft.dos import DOS # ASE calss for calculating the density of states (DOS)
import matplotlib.pyplot as plt
atoms = read('cell_siesta.traj')
atoms.pbc=True
os.mkdir("siesta_dos")
os.chdir("siesta_dos")
atoms.calc = Siesta(label='diamond',
xc='PBE',
pseudo_path='/Users/staykov/pseudodojo',
pseudo_qualifier='gga',
symlink_pseudos=True,
mesh_cutoff=230 * Ry,
energy_shift=0.01 * Ry,
basis_set='DZP',
spin='non-polarized',
kpts=(5, 5, 5),
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 100, 'WriteEigenvalues': True, 'SaveHS': True,},
)
energy = atoms.get_potential_energy()
dos = DOS(atoms.calc,
width=0.05,
npts=3000)
E = dos.get_energies()
D = dos.get_dos()
os.chdir("../")
# Select only -10 to +10 eV around Ef
mask = (E >= -10) & (E <= 10)
plt.plot(E[mask], D[mask])
plt.axvline(0, color='k', linestyle='--')
plt.xlim(-10, 10)
plt.xlabel(r'$E - E_F$ (eV)')
plt.ylabel('DOS (states/eV)')
plt.show()
Job completed
While bulk calculations are essential for understanding intrinsic material properties, many applications in catalysis, electronics, and nanotechnology require detailed knowledge of surface properties. In periodic DFT codes such as SIESTA, surfaces are modeled using slab geometries, where a finite number of atomic layers represent the surface and vacuum is added to isolate it from periodic images. ASE provides convenient tools for constructing such slabs. Using its built‑in ase.build.surface function, one can cleave a bulk crystal along a chosen set of Miller indices. In this tutorial, we generate the diamond (100) surface by slicing the bulk diamond structure along the (1 0 0) plane and adding sufficient vacuum to avoid interactions between repeated slabs. Once the slab is constructed, the geometry optimization is applied only to the atomic positions, not to the cell parameters. The lattice constants are optimized at the bulk level and then kept fixed during surface relaxation. This ensures that the slab inherits physically meaningful bulk geometry while allowing the surface layers to relax into their energetically preferred configuration.
from ase import Atoms
from ase.calculators.siesta import Siesta
from ase.units import Ry
from ase.build import surface
from ase.visualize import view
from ase.io import read
from ase.io import write
import os
from ase.optimize.bfgs import BFGS
atoms = read('cell_siesta.traj')
atoms.pbc=True
slab = surface(atoms, (1,0,0), 2, vacuum=10.0) # create a slab of the optimized structure with 4 layers and 15 Å of vacuum
slab.pbc=True
# center slab in vacuum
slab.center(axis=2)
os.mkdir("siesta_slab_opt")
os.chdir("siesta_slab_opt")
# Minimal Vasp calculator for testing
slab.calc = Siesta(label='diamond',
xc='PBE',
pseudo_path='/Users/staykov/pseudodojo',
pseudo_qualifier='gga',
symlink_pseudos=True,
mesh_cutoff=230 * Ry,
energy_shift=0.01 * Ry,
basis_set='DZP',
spin='non-polarized',
kpts=(5, 5, 1),
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 100},
)
energy = slab.get_total_energy()
opt = BFGS(slab, trajectory='cellopt.traj')
opt.run(fmax=0.03)
# Run a single-point calculation
os.chdir("../")
write('slab_siesta.traj', slab)
view(slab, viewer='x3d')
Job completed
Step Time Energy fmax BFGS: 0 00:49:41 -2603.623001 2.261051
Job completed
BFGS: 1 00:50:26 -2603.969982 0.840002
Job completed
BFGS: 2 00:51:33 -2604.021646 0.202029
Job completed
BFGS: 3 00:52:20 -2604.025396 0.139162
Job completed
BFGS: 4 00:52:59 -2604.027757 0.095486
Job completed
BFGS: 5 00:53:28 -2604.029964 0.127078
Job completed
BFGS: 6 00:54:01 -2604.031869 0.096370
Job completed
BFGS: 7 00:54:35 -2604.032464 0.031558 BFGS: 8 00:54:55 -2604.032525 0.003870
Job completed
The surface Density of States (DOS) provides essential insight into how a surface modifies the electronic structure of a material. Unlike the bulk, where all bonds are fully coordinated, a surface exposes atoms with unsatisfied valence—leading to electronic features that can differ dramatically from the bulk behavior. The DOS of the diamond (100) surface illustrates this clearly. Although bulk diamond is a wide‑band‑gap insulator, the (100) surface shows no band gap and exhibits metallic character. At first glance this result may seem counterintuitive. However, the explanation lies in the orbital configuration of the surface atoms. The top‑layer carbon atoms possess dangling bonds that are not compensated by neighboring atoms. These dangling‑bond states fall within the band gap of bulk diamond and create partially filled surface states, giving rise to metallic behavior. This metallicity is not just a computational artifact—it is consistent with experimental observations, where clean diamond surfaces often display surface conductivity due to these unsaturated bonds. Understanding this effect is crucial when studying surface chemistry, adsorption, catalysis, or electronic devices based on diamond surfaces.
from ase import Atoms
from ase.calculators.siesta import Siesta
from ase.units import Ry
import os
from ase.io import read
from ase.dft.dos import DOS
import matplotlib.pyplot as plt
atoms = read('slab_siesta.traj')
atoms.pbc=True
os.mkdir("siesta_dos_slab")
os.chdir("siesta_dos_slab")
atoms.calc = Siesta(label='diamond',
xc='PBE',
pseudo_path='/Users/staykov/pseudodojo',
pseudo_qualifier='gga',
symlink_pseudos=True,
mesh_cutoff=230 * Ry,
energy_shift=0.01 * Ry,
basis_set='DZP',
spin='non-polarized',
kpts=(5, 5, 1),
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 100, 'WriteEigenvalues': True, 'SaveHS': True,},
)
energy = atoms.get_potential_energy()
dos = DOS(atoms.calc,
width=0.05,
npts=3000)
E = dos.get_energies()
D = dos.get_dos()
os.chdir("../")
# Select only -10 to +10 eV around Ef
mask = (E >= -10) & (E <= 10)
plt.plot(E[mask], D[mask])
plt.axvline(0, color='k', linestyle='--')
plt.xlim(-10, 10)
plt.xlabel(r'$E - E_F$ (eV)')
plt.ylabel('DOS (states/eV)')
plt.show()
Job completed
In reality, a pristine diamond surface exists only under ultrahigh‑vacuum (UHV) conditions and even then only for short periods of time. Under ambient or experimental conditions, the surface rapidly becomes passivated, most commonly by hydrogen atoms. To model realistic surfaces, we therefore construct a hydrogen‑terminated diamond (100) slab and optimize its geometry again. ASE provides convenient tools for building such passivated surfaces. First, we identify the top and bottom surface carbon atoms, which are the ones carrying dangling bonds. Hydrogen atoms are then added to these sites to saturate the dangling bonds and restore the insulating character of the surface. To avoid artificial interactions between hydrogens, we arrange them in a pattern that maximizes the distance between neighboring H atoms, ensuring a physically meaningful termination. After adding the hydrogens, the slab is relaxed again, allowing both the surface carbon atoms and the hydrogen atoms to settle into their energetically preferred configuration. This passivated surface serves as a more realistic model for studying electronic properties, adsorption, and surface chemistry under typical experimental conditions.
import numpy as np
from ase.io import read
from ase.build import surface
from ase import Atom
from ase.visualize import view
atoms = read('cell_siesta.traj')
atoms.pbc = True
slab = surface(atoms, (1,0,0), 2, vacuum=10.0)
slab.center(axis=2)
# Highest carbon z coordinate
zmax = max(atom.position[2] for atom in slab if atom.symbol == 'C')
# Select top-layer carbon atoms
tol = 0.5 # Å tolerance
top_carbons = [i for i, atom in enumerate(slab)
if atom.symbol == 'C' and zmax - atom.position[2] < tol]
CH = 0.8
pos = slab[top_carbons[0]].position + np.array([-0.5, -0.5, CH])
slab.append(Atom('H', pos))
pos = slab[top_carbons[0]].position + np.array([0.5, 0.5, CH])
slab.append(Atom('H', pos))
pos = slab[top_carbons[1]].position + np.array([-0.5, 0.5, CH])
slab.append(Atom('H', pos))
pos = slab[top_carbons[1]].position + np.array([0.5, -0.5, CH])
slab.append(Atom('H', pos))
# Highest carbon z coordinate
zmin = min(atom.position[2] for atom in slab if atom.symbol == 'C')
top_carbons = [i for i, atom in enumerate(slab)
if atom.symbol == 'C' and atom.position[2] - zmin < tol]
pos = slab[top_carbons[0]].position - np.array([-0.5, -0.5, CH])
slab.append(Atom('H', pos))
pos = slab[top_carbons[0]].position - np.array([0.5, 0.5, CH])
slab.append(Atom('H', pos))
pos = slab[top_carbons[1]].position - np.array([-0.5, 0.5, CH])
slab.append(Atom('H', pos))
pos = slab[top_carbons[1]].position - np.array([0.5, -0.5, CH])
slab.append(Atom('H', pos))
view(slab, viewer='x3d')
Self‑consistent field (SCF) convergence and geometry relaxation in LCAO‑based DFT codes such as SIESTA depend strongly on the quality of the initial atomic configuration. When the starting geometry is far from the true minimum, the SCF cycle may become unstable, stall, or fail to converge entirely. This is exactly what happens for the H‑passivated diamond (100) surface: the initial guess contains significant strain and unfavorable bond orientations, making the SCF procedure overly sensitive. A practical and widely used strategy in SIESTA is to temporarily loosen the computational settings to help the optimizer explore a broader region of configuration space. Reducing the basis‑set size, lowering the k‑point density, and decreasing the real‑space mesh cutoff all make the SCF cycle more forgiving. Although these settings are not suitable for final production calculations, they allow the geometry optimizer to move the system closer to a physically meaningful structure where full‑accuracy parameters can later be restored. In this tutorial, we apply this approach by reducing the k‑point mesh to 3×3×3 and lowering the real‑space grid cutoff to 170 Ry. With these relaxed settings, the SCF converges reliably, enabling the geometry optimizer to correct the initial distortions. Once the structure is sufficiently improved, the calculation can be repeated with fully converged parameters to obtain accurate energies and electronic properties.
import numpy as np
from ase.io import read
from ase.io import write
from ase.build import surface
from ase import Atom
from ase.visualize import view
from ase.calculators.siesta import Siesta
from ase.units import Ry
from ase.optimize.bfgs import BFGS
import os
atoms = read('cell_siesta.traj')
atoms.pbc = True
slab = surface(atoms, (1,0,0), 2, vacuum=10.0)
slab.center(axis=2)
# Highest carbon z coordinate
zmax = max(atom.position[2] for atom in slab if atom.symbol == 'C')
# Select top-layer carbon atoms
tol = 0.5 # Å tolerance
top_carbons = [i for i, atom in enumerate(slab)
if atom.symbol == 'C' and zmax - atom.position[2] < tol]
# Add H atoms 1.09 Å above each top carbon
CH = 0.8
pos = slab[top_carbons[0]].position + np.array([-0.5, -0.5, CH])
slab.append(Atom('H', pos))
pos = slab[top_carbons[0]].position + np.array([0.5, 0.5, CH])
slab.append(Atom('H', pos))
pos = slab[top_carbons[1]].position + np.array([-0.5, 0.5, CH])
slab.append(Atom('H', pos))
pos = slab[top_carbons[1]].position + np.array([0.5, -0.5, CH])
slab.append(Atom('H', pos))
# Highest carbon z coordinate
zmin = min(atom.position[2] for atom in slab if atom.symbol == 'C')
top_carbons = [i for i, atom in enumerate(slab)
if atom.symbol == 'C' and atom.position[2] - zmin < tol]
pos = slab[top_carbons[0]].position - np.array([-0.5, -0.5, CH])
slab.append(Atom('H', pos))
pos = slab[top_carbons[0]].position - np.array([0.5, 0.5, CH])
slab.append(Atom('H', pos))
pos = slab[top_carbons[1]].position - np.array([-0.5, 0.5, CH])
slab.append(Atom('H', pos))
pos = slab[top_carbons[1]].position - np.array([0.5, -0.5, CH])
slab.append(Atom('H', pos))
os.mkdir("siesta_slab_h_opt")
os.chdir("siesta_slab_h_opt")
# Minimal Vasp calculator for testing
slab.calc = Siesta(label='diamond',
xc='PBE',
pseudo_path='/Users/staykov/pseudodojo',
pseudo_qualifier='gga',
symlink_pseudos=True,
mesh_cutoff=170 * Ry,
energy_shift=0.01 * Ry,
basis_set='SZP',
spin='non-polarized',
kpts=(3, 3, 1),
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 100},
)
energy = slab.get_total_energy()
opt = BFGS(slab, trajectory='cellopt.traj')
opt.run(fmax=0.03)
# Run a single-point calculation
os.chdir("../")
write('slab_h_siesta.traj', slab)
view(slab, viewer='x3d')
Job completed
Step Time Energy fmax BFGS: 0 01:01:06 -2718.438802 4.527510
Job completed
BFGS: 1 01:01:38 -2719.981279 1.780527
Job completed
BFGS: 2 01:01:59 -2720.481454 1.434775
Job completed
BFGS: 3 01:02:21 -2720.991475 1.375183
Job completed
BFGS: 4 01:02:40 -2721.247333 1.085346
Job completed
BFGS: 5 01:03:10 -2721.500232 0.657005
Job completed
BFGS: 6 01:03:48 -2721.619501 0.446299
Job completed
BFGS: 7 01:04:22 -2721.664135 0.329834
Job completed
BFGS: 8 01:04:53 -2721.676295 0.236303
Job completed
BFGS: 9 01:05:19 -2721.686622 0.193424
Job completed
BFGS: 10 01:05:52 -2721.703653 0.265693
Job completed
BFGS: 11 01:06:28 -2721.740113 0.258966
Job completed
BFGS: 12 01:07:02 -2721.406243 2.621486
Job completed
BFGS: 13 01:07:35 -2721.781388 0.208029
Job completed
BFGS: 14 01:07:53 -2721.792243 0.242205
Job completed
BFGS: 15 01:08:11 -2721.799845 0.237188
Job completed
BFGS: 16 01:08:35 -2721.809403 0.307038
Job completed
BFGS: 17 01:08:48 -2721.815537 0.279699
Job completed
BFGS: 18 01:09:09 -2721.826765 0.190260
Job completed
BFGS: 19 01:09:23 -2721.829630 0.177947
Job completed
BFGS: 20 01:09:43 -2721.835275 0.122043
Job completed
BFGS: 21 01:10:00 -2721.838567 0.070469
Job completed
BFGS: 22 01:10:10 -2721.840696 0.046812
Job completed
BFGS: 23 01:10:20 -2721.841084 0.046455
Job completed
BFGS: 24 01:10:26 -2721.841148 0.044807
Job completed
BFGS: 25 01:10:32 -2721.841227 0.040622
Job completed
BFGS: 26 01:10:42 -2721.841704 0.110535
Job completed
BFGS: 27 01:11:06 -2721.837149 0.453669
Job completed
BFGS: 28 01:11:31 -2721.845083 0.212256
Job completed
BFGS: 29 01:11:42 -2721.846595 0.270620
Job completed
BFGS: 30 01:11:58 -2721.848010 0.328499
Job completed
BFGS: 31 01:12:15 -2721.849624 0.399308
Job completed
BFGS: 32 01:12:28 -2721.851475 0.462431
Job completed
BFGS: 33 01:12:50 -2721.860179 0.877147
Job completed
BFGS: 34 01:13:20 -2721.878957 1.343878
Job completed
BFGS: 35 01:13:37 -2722.131107 2.006929
Job completed
BFGS: 36 01:13:52 -2722.643428 2.200675
Job completed
BFGS: 37 01:14:05 -2723.305788 4.960323
Job completed
BFGS: 38 01:14:19 -2723.969136 5.566324
Job completed
BFGS: 39 01:14:32 -2725.158783 3.405078
Job completed
BFGS: 40 01:14:44 -2724.906573 5.368603
Job completed
BFGS: 41 01:14:58 -2725.634297 4.582136
Job completed
BFGS: 42 01:15:12 -2726.042963 1.420036
Job completed
BFGS: 43 01:15:21 -2726.145033 0.815896
Job completed
BFGS: 44 01:15:32 -2726.271770 0.830943
Job completed
BFGS: 45 01:15:42 -2726.354121 0.691162
Job completed
BFGS: 46 01:15:50 -2726.413416 0.639350
Job completed
BFGS: 47 01:16:00 -2726.479910 0.574065
Job completed
BFGS: 48 01:16:12 -2726.547338 0.714338
Job completed
BFGS: 49 01:16:24 -2726.607379 0.665869
Job completed
BFGS: 50 01:16:36 -2726.666391 0.486329
Job completed
BFGS: 51 01:16:45 -2726.695615 0.305305
Job completed
BFGS: 52 01:16:54 -2726.716234 0.347288
Job completed
BFGS: 53 01:17:01 -2726.728957 0.310198
Job completed
BFGS: 54 01:17:10 -2726.743628 0.375634
Job completed
BFGS: 55 01:17:19 -2726.766859 0.436328
Job completed
BFGS: 56 01:17:29 -2726.814164 0.743877
Job completed
BFGS: 57 01:17:40 -2726.880010 1.112112
Job completed
BFGS: 58 01:17:50 -2726.984994 1.125572
Job completed
BFGS: 59 01:18:03 -2727.197017 1.024027
Job completed
BFGS: 60 01:18:13 -2727.260821 0.897828
Job completed
BFGS: 61 01:18:22 -2727.315533 0.568791
Job completed
BFGS: 62 01:18:34 -2727.326522 0.697510
Job completed
BFGS: 63 01:18:44 -2727.366079 0.414972
Job completed
BFGS: 64 01:18:52 -2727.389595 0.371964
Job completed
BFGS: 65 01:19:02 -2727.436705 0.332142
Job completed
BFGS: 66 01:19:11 -2727.455800 0.372752
Job completed
BFGS: 67 01:19:20 -2727.469173 0.330448
Job completed
BFGS: 68 01:19:28 -2727.476056 0.230206
Job completed
BFGS: 69 01:19:34 -2727.480233 0.243488
Job completed
BFGS: 70 01:19:41 -2727.484216 0.264698
Job completed
BFGS: 71 01:19:49 -2727.490884 0.330562
Job completed
BFGS: 72 01:19:56 -2727.502728 0.493365
Job completed
BFGS: 73 01:20:05 -2727.529346 0.863879
Job completed
BFGS: 74 01:20:15 -2727.568042 1.181201
Job completed
BFGS: 75 01:20:24 -2727.642074 1.151166
Job completed
BFGS: 76 01:20:37 -2727.811166 1.078242
Job completed
BFGS: 77 01:20:47 -2727.918579 0.751178
Job completed
BFGS: 78 01:21:00 -2727.988067 0.781713
Job completed
BFGS: 79 01:21:11 -2728.003279 0.887910
Job completed
BFGS: 80 01:21:22 -2728.071650 0.808238
Job completed
BFGS: 81 01:21:31 -2728.107872 0.420869
Job completed
BFGS: 82 01:21:40 -2728.152838 0.432755
Job completed
BFGS: 83 01:21:50 -2728.203863 0.422761
Job completed
BFGS: 84 01:21:59 -2728.228764 0.349727
Job completed
BFGS: 85 01:22:07 -2728.244182 0.303725
Job completed
BFGS: 86 01:22:16 -2728.252549 0.334123
Job completed
BFGS: 87 01:22:25 -2728.259661 0.336627
Job completed
BFGS: 88 01:22:32 -2728.264756 0.206908
Job completed
BFGS: 89 01:22:39 -2728.267364 0.199678
Job completed
BFGS: 90 01:22:47 -2728.268671 0.146715
Job completed
BFGS: 91 01:22:54 -2728.269545 0.139585
Job completed
BFGS: 92 01:23:01 -2728.270251 0.128328
Job completed
BFGS: 93 01:23:07 -2728.271134 0.117841
Job completed
BFGS: 94 01:23:14 -2728.272022 0.108062
Job completed
BFGS: 95 01:23:21 -2728.272921 0.105820
Job completed
BFGS: 96 01:23:27 -2728.273390 0.096876
Job completed
BFGS: 97 01:23:33 -2728.273624 0.088617
Job completed
BFGS: 98 01:23:39 -2728.273807 0.083671
Job completed
BFGS: 99 01:23:44 -2728.274103 0.080260
Job completed
BFGS: 100 01:23:52 -2728.274529 0.076242
Job completed
BFGS: 101 01:23:59 -2728.274973 0.068225
Job completed
BFGS: 102 01:24:05 -2728.275281 0.059409
Job completed
BFGS: 103 01:24:11 -2728.275438 0.056594
Job completed
BFGS: 104 01:24:16 -2728.275537 0.055067
Job completed
BFGS: 105 01:24:20 -2728.275617 0.055394
Job completed
BFGS: 106 01:24:26 -2728.275684 0.054729
Job completed
BFGS: 107 01:24:31 -2728.275737 0.054082
Job completed
BFGS: 108 01:24:38 -2728.275795 0.052638
Job completed
BFGS: 109 01:24:43 -2728.275858 0.050996
Job completed
BFGS: 110 01:24:48 -2728.275955 0.050856
Job completed
BFGS: 111 01:24:54 -2728.276063 0.050295
Job completed
BFGS: 112 01:25:01 -2728.276184 0.050496
Job completed
BFGS: 113 01:25:08 -2728.276303 0.054978
Job completed
BFGS: 114 01:25:14 -2728.276410 0.054426
Job completed
BFGS: 115 01:25:21 -2728.276524 0.048967
Job completed
BFGS: 116 01:25:27 -2728.276633 0.050109
Job completed
BFGS: 117 01:25:33 -2728.276714 0.047601
Job completed
BFGS: 118 01:25:38 -2728.276775 0.044655
Job completed
BFGS: 119 01:25:44 -2728.276817 0.041925
Job completed
BFGS: 120 01:25:49 -2728.276854 0.040035
Job completed
BFGS: 121 01:25:59 -2728.276900 0.038253
Job completed
BFGS: 122 01:26:04 -2728.276950 0.037064
Job completed
BFGS: 123 01:26:11 -2728.277008 0.038722
Job completed
BFGS: 124 01:26:17 -2728.277053 0.040113
Job completed
BFGS: 125 01:26:22 -2728.277093 0.041086
Job completed
BFGS: 126 01:26:28 -2728.277132 0.043760
Job completed
BFGS: 127 01:26:34 -2728.277198 0.047269
Job completed
BFGS: 128 01:26:40 -2728.277336 0.051881
Job completed
BFGS: 129 01:26:46 -2728.277589 0.063843
Job completed
BFGS: 130 01:26:53 -2728.277970 0.065306
Job completed
BFGS: 131 01:27:00 -2728.278301 0.063479
Job completed
BFGS: 132 01:27:07 -2728.278567 0.060112
Job completed
BFGS: 133 01:27:14 -2728.278685 0.058747
Job completed
BFGS: 134 01:27:19 -2728.278779 0.053158
Job completed
BFGS: 135 01:27:25 -2728.278889 0.054691
Job completed
BFGS: 136 01:27:32 -2728.278974 0.052589
Job completed
BFGS: 137 01:27:38 -2728.279057 0.050583
Job completed
BFGS: 138 01:27:45 -2728.279138 0.046574
Job completed
BFGS: 139 01:27:52 -2728.279221 0.041974
Job completed
BFGS: 140 01:27:58 -2728.279340 0.038935
Job completed
BFGS: 141 01:28:05 -2728.279394 0.037179
Job completed
BFGS: 142 01:28:12 -2728.279437 0.036684
Job completed
BFGS: 143 01:28:18 -2728.279458 0.035480
Job completed
BFGS: 144 01:28:24 -2728.279475 0.036196
Job completed
BFGS: 145 01:28:29 -2728.279491 0.037129
Job completed
BFGS: 146 01:28:35 -2728.279517 0.037168
Job completed
BFGS: 147 01:28:40 -2728.279540 0.035390
Job completed
BFGS: 148 01:28:47 -2728.279562 0.037623
Job completed
BFGS: 149 01:28:51 -2728.279574 0.037406
Job completed
BFGS: 150 01:28:57 -2728.279581 0.036522
Job completed
BFGS: 151 01:29:03 -2728.279587 0.036862
Job completed
BFGS: 152 01:29:09 -2728.279595 0.036729
Job completed
BFGS: 153 01:29:14 -2728.279606 0.036267
Job completed
BFGS: 154 01:29:21 -2728.279624 0.035423
Job completed
BFGS: 155 01:29:29 -2728.279653 0.033351
Job completed
BFGS: 156 01:29:35 -2728.279695 0.033305
Job completed
BFGS: 157 01:29:41 -2728.279747 0.032673
Job completed
BFGS: 158 01:29:48 -2728.279776 0.030369 BFGS: 159 01:29:54 -2728.279793 0.027660
Job completed
The fully optimized geometry differs substantially from our initial guess, confirming that the original configuration was far from the true minimum. After relaxing the structure with the loosened computational settings, we restored the high‑accuracy parameters. At this stage, the system was already close to equilibrium, and the SCF cycle was stable. As a result, the final high‑precision geometry optimization converged in only a few steps, demonstrating the effectiveness of the two‑stage relaxation strategy.
import numpy as np
from ase.io import read
from ase.io import write
from ase import Atom
from ase.visualize import view
from ase.calculators.siesta import Siesta
from ase.units import Ry
from ase.optimize.bfgs import BFGS
import os
atoms = read('slab_h_siesta.traj')
atoms.pbc = True
os.mkdir("siesta_slab_h_opt_high")
os.chdir("siesta_slab_h_opt_high")
# Minimal Vasp calculator for testing
slab.calc = Siesta(label='diamond',
xc='PBE',
pseudo_path='/Users/staykov/pseudodojo',
pseudo_qualifier='gga',
symlink_pseudos=True,
mesh_cutoff=230 * Ry,
energy_shift=0.01 * Ry,
basis_set='DZP',
spin='non-polarized',
kpts=(5, 5, 1),
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 100},
)
energy = slab.get_total_energy()
opt = BFGS(slab, trajectory='cellopt.traj')
opt.run(fmax=0.03)
# Run a single-point calculation
os.chdir("../")
write('slab_h_siesta_high.traj', slab)
view(slab, viewer='x3d')
Job completed
Step Time Energy fmax BFGS: 0 01:31:46 -2739.771834 0.701842
Job completed
BFGS: 1 01:32:20 -2739.812459 0.350674
Job completed
BFGS: 2 01:32:57 -2739.826412 0.268452
Job completed
BFGS: 3 01:33:23 -2739.834422 0.192907
Job completed
BFGS: 4 01:33:47 -2739.840794 0.165390
Job completed
BFGS: 5 01:34:10 -2739.846731 0.154870
Job completed
BFGS: 6 01:34:35 -2739.851928 0.130103
Job completed
BFGS: 7 01:34:58 -2739.854433 0.066471
Job completed
BFGS: 8 01:35:19 -2739.855451 0.059154
Job completed
BFGS: 9 01:35:36 -2739.856139 0.049899
Job completed
BFGS: 10 01:35:54 -2739.856772 0.051844
Job completed
BFGS: 11 01:36:12 -2739.857209 0.041879
Job completed
BFGS: 12 01:36:30 -2739.857536 0.030676
Job completed
BFGS: 13 01:36:47 -2739.857850 0.034641
Job completed
BFGS: 14 01:37:05 -2739.858110 0.033096 BFGS: 15 01:37:23 -2739.858247 0.020238
Job completed
The Density of States (DOS) of the hydrogen‑passivated diamond (100) surface shows a large band gap, closely matching that of bulk diamond. This behavior is expected: once the surface carbon atoms are saturated with hydrogen, the dangling‑bond states responsible for the metallic character of the pristine surface are removed. With these mid‑gap states eliminated, the electronic structure of the slab returns to that of an insulator, demonstrating that hydrogen termination effectively restores the bulk‑like electronic properties at the surface.
from ase import Atoms
from ase.calculators.siesta import Siesta
from ase.units import Ry
import os
from ase.io import read
from ase.dft.dos import DOS
import matplotlib.pyplot as plt
atoms = read('slab_h_siesta_high.traj')
atoms.pbc=True
os.mkdir("siesta_dos_slab_h")
os.chdir("siesta_dos_slab_h")
atoms.calc = Siesta(label='diamond',
xc='PBE',
pseudo_path='/Users/staykov/pseudodojo',
pseudo_qualifier='gga',
symlink_pseudos=True,
mesh_cutoff=230 * Ry,
energy_shift=0.01 * Ry,
basis_set='DZP',
spin='non-polarized',
kpts=(5, 5, 1),
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 100, 'WriteEigenvalues': True, 'SaveHS': True,},
)
energy = atoms.get_potential_energy()
dos = DOS(atoms.calc,
width=0.05,
npts=3000)
E = dos.get_energies()
D = dos.get_dos()
os.chdir("../")
# Select only -10 to +10 eV around Ef
mask = (E >= -10) & (E <= 10)
plt.plot(E[mask], D[mask])
plt.axvline(0, color='k', linestyle='--')
plt.xlim(-10, 10)
plt.xlabel(r'$E - E_F$ (eV)')
plt.ylabel('DOS (states/eV)')
plt.show()
Job completed
The Density of States (DOS) tells us how electronic states are distributed relative to the Fermi level, but it does not provide the absolute energy of the Fermi level. In periodic DFT calculations, all eigenvalues are referenced to an arbitrary internal zero, determined by the basis set, pseudopotentials, and numerical setup. As a result, absolute energy levels cannot be compared directly between two different materials or even between two separate calculations of the same material using different computational parameters. To obtain meaningful, comparable energy references, we compute the work function or, more precisely, the ionization potential of the slab. This quantity gives the energy of the valence‑band maximum (VBM) relative to the vacuum level, which is the same physical reference for every slab. By anchoring the electronic structure to the vacuum level, we can compare band edges, Fermi levels, and work functions across different systems in a consistent way. To perform work‑function or ionization‑potential calculations, the slab must include a sufficiently thick vacuum region. The electrostatic potential must reach a flat plateau in the vacuum, indicating that the potential is constant and free from interactions with the slab. Once this plateau is identified, we extract the vacuum electrostatic potential and reference the VBM (or Fermi level) to it. SIESTA outputs the electrostatic potential on a real‑space grid in the file ElectrostaticPotential.grid.nc, which contains the volumetric data needed to compute the averaged potential along the slab normal. To generate this file, several options must be enabled in the SIESTA input: 'WriteEigenvalues': True, 'SaveHS': True, 'SaveElectrostaticPotential': True, 'WriteVH': True,
import numpy as np
import matplotlib.pyplot as plt
import re
from ase import Atoms
from ase.calculators.siesta import Siesta
from ase.units import Ry
import os
from ase.io import read
slab = read('slab_h_siesta_high.traj')
slab.pbc=True
os.mkdir("siesta_wf_slab_h")
os.chdir("siesta_wf_slab_h")
slab.calc = Siesta(label='slab',
xc='PBE',
pseudo_path='/Users/staykov/pseudodojo',
pseudo_qualifier='gga',
symlink_pseudos=True,
mesh_cutoff=230 * Ry,
energy_shift=0.01 * Ry,
basis_set='DZP',
spin='non-polarized',
kpts=(5, 5, 1),
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 100, 'WriteEigenvalues': True, 'SaveHS': True,'SaveElectrostaticPotential': True, 'WriteVH': True,},
)
energy = slab.get_potential_energy()
import sisl
import numpy as np
import matplotlib.pyplot as plt
# --------------------------------------------------
# Read electrostatic potential
# --------------------------------------------------
grid = sisl.get_sile("ElectrostaticPotential.grid.nc").read_grid()
V = grid.grid
# planar average
V_z = V.mean(axis=(0,1))
# z coordinate
cell_z = grid.cell[2,2]
z = np.linspace(0, cell_z, len(V_z))
# --------------------------------------------------
# Vacuum level
# --------------------------------------------------
# for symmetric slab simply use the maximum value
V_vac = np.max(V_z)
# alternatively use averaging around plateau:
# vacuum_mask = z > (cell_z - 5.0)
# V_vac = np.mean(V_z[vacuum_mask])
print(f"Vacuum level = {V_vac:.6f} eV")
# --------------------------------------------------
# Read eigenvalues
# --------------------------------------------------
with open("slab.EIG") as f:
lines = f.readlines()
# first line contains Fermi level
E_F = float(lines[0])
# collect all eigenvalues
eigs = []
for line in lines[2:]:
parts = line.split()
# skip band index at beginning of each k-point block
if len(parts) > 0:
try:
int(parts[0])
parts = parts[1:]
except:
pass
for x in parts:
try:
eigs.append(float(x))
except:
pass
eigs = np.array(eigs)
# --------------------------------------------------
# Determine VBM and CBM
# --------------------------------------------------
occupied = eigs[eigs < E_F]
unoccupied = eigs[eigs > E_F]
VBM = occupied.max()
CBM = unoccupied.min()
print(f"Fermi level = {E_F:.6f} eV")
print(f"VBM = {VBM:.6f} eV")
print(f"CBM = {CBM:.6f} eV")
# --------------------------------------------------
# Ionization potential
# --------------------------------------------------
IP = V_vac - VBM
# Electron affinity if desired
EA = V_vac - CBM
print()
print(f"Vacuum - VBM (Ionization Potential) = {IP:.6f} eV")
print(f"Vacuum - CBM (Electron Affinity) = {EA:.6f} eV")
# --------------------------------------------------
# Plot electrostatic potential
# --------------------------------------------------
plt.figure(figsize=(8,5))
plt.plot(z, V_z, lw=2)
plt.axhline(
V_vac,
linestyle='--',
label=f'Vacuum level = {V_vac:.2f} eV'
)
plt.axhline(
VBM,
linestyle='--',
label=f'VBM = {VBM:.2f} eV'
)
plt.xlabel("z (Å)")
plt.ylabel("Electrostatic potential (eV)")
plt.title("Planar averaged electrostatic potential")
plt.legend()
plt.tight_layout()
plt.show()
os.chdir("../")
Job completed
Vacuum level = -0.908470 eV Fermi level = -1.742358 eV VBM = -3.487912 eV CBM = 0.377073 eV Vacuum - VBM (Ionization Potential) = 2.579443 eV Vacuum - CBM (Electron Affinity) = -1.285543 eV
The ionization potential of the hydrogen‑passivated diamond (100) slab is calculated to be 2.58 eV, whereas the pristine diamond (100) surface exhibits a much higher value of 6.66 eV. This dramatic difference is well‑known experimentally and highlights the crucial role of surface reconstruction and surface termination in determining electronic properties. Hydrogen termination fundamentally alters the surface electronic environment. By saturating the dangling bonds on the topmost carbon atoms, hydrogen atoms create a surface dipole layer that shifts the electrostatic potential downward. This dipole makes it energetically easier to remove an electron from the surface, resulting in a much lower ionization potential. In practical terms, the H‑terminated surface becomes far more favorable for electron emission, a property exploited in applications such as negative‑electron‑affinity (NEA) diamond devices. In contrast, the pristine diamond (100) surface lacks these stabilizing dipoles. Its unsaturated dangling bonds produce mid‑gap states and a higher surface potential, making electron extraction significantly more difficult. The high ionization potential reflects this unfavorable electronic environment. These results demonstrate how surface chemistry directly controls the electronic structure, and they emphasize the importance of accurate surface modeling when studying real materials or designing diamond‑based electronic and optoelectronic devices.
import numpy as np
import matplotlib.pyplot as plt
import re
from ase import Atoms
from ase.calculators.siesta import Siesta
from ase.units import Ry
import os
from ase.io import read
import sisl
import numpy as np
import matplotlib.pyplot as plt
slab = read('slab_siesta.traj')
slab.pbc=True
os.mkdir("csiesta_wf_slab_h_non")
os.chdir("csiesta_wf_slab_h_non")
slab.calc = Siesta(label='slab',
xc='PBE',
pseudo_path='/Users/staykov/pseudodojo',
pseudo_qualifier='gga',
symlink_pseudos=True,
mesh_cutoff=230 * Ry,
energy_shift=0.01 * Ry,
basis_set='DZP',
spin='non-polarized',
kpts=(5, 5, 1),
fdf_arguments={'DM.MixingWeight': 0.1, 'MaxSCFIterations': 100, 'WriteEigenvalues': True, 'SaveHS': True,'SaveElectrostaticPotential': True, 'WriteVH': True,},
)
energy = slab.get_potential_energy()
# --------------------------------------------------
# Read electrostatic potential
# --------------------------------------------------
grid = sisl.get_sile("ElectrostaticPotential.grid.nc").read_grid()
V = grid.grid
# planar average
V_z = V.mean(axis=(0,1))
# z coordinate
cell_z = grid.cell[2,2]
z = np.linspace(0, cell_z, len(V_z))
# --------------------------------------------------
# Vacuum level
# --------------------------------------------------
# for symmetric slab simply use the maximum value
V_vac = np.max(V_z)
# alternatively use averaging around plateau:
# vacuum_mask = z > (cell_z - 5.0)
# V_vac = np.mean(V_z[vacuum_mask])
print(f"Vacuum level = {V_vac:.6f} eV")
# --------------------------------------------------
# Read eigenvalues
# --------------------------------------------------
with open("slab.EIG") as f:
lines = f.readlines()
# first line contains Fermi level
E_F = float(lines[0])
# collect all eigenvalues
eigs = []
for line in lines[2:]:
parts = line.split()
# skip band index at beginning of each k-point block
if len(parts) > 0:
try:
int(parts[0])
parts = parts[1:]
except:
pass
for x in parts:
try:
eigs.append(float(x))
except:
pass
eigs = np.array(eigs)
# --------------------------------------------------
# Determine VBM and CBM
# --------------------------------------------------
occupied = eigs[eigs < E_F]
unoccupied = eigs[eigs > E_F]
VBM = occupied.max()
CBM = unoccupied.min()
print(f"Fermi level = {E_F:.6f} eV")
print(f"VBM = {VBM:.6f} eV")
print(f"CBM = {CBM:.6f} eV")
# --------------------------------------------------
# Ionization potential
# --------------------------------------------------
IP = V_vac - VBM
# Electron affinity if desired
EA = V_vac - CBM
print()
print(f"Vacuum - VBM (Ionization Potential) = {IP:.6f} eV")
print(f"Vacuum - CBM (Electron Affinity) = {EA:.6f} eV")
# --------------------------------------------------
# Plot electrostatic potential
# --------------------------------------------------
plt.figure(figsize=(8,5))
plt.plot(z, V_z, lw=2)
plt.axhline(
V_vac,
linestyle='--',
label=f'Vacuum level = {V_vac:.2f} eV'
)
plt.axhline(
VBM,
linestyle='--',
label=f'VBM = {VBM:.2f} eV'
)
plt.xlabel("z (Å)")
plt.ylabel("Electrostatic potential (eV)")
plt.title("Planar averaged electrostatic potential")
plt.legend()
plt.tight_layout()
plt.show()
os.chdir("../")
Vacuum level = 0.459204 eV Fermi level = -6.185412 eV VBM = -6.200983 eV CBM = -6.181915 eV Vacuum - VBM (Ionization Potential) = 6.660187 eV Vacuum - CBM (Electron Affinity) = 6.641119 eV
Job completed
