Electronic Structure: DOS, Magnetism, and Charge

Electronic structure analysis takes the question "What is the total energy of this structure" one step further: At which energies are there nearby electronic states? What atoms or orbitals do they come from? How do electronic states change with crystal momentum? The answers can help understand bonding, magnetism, and reactivity, but a DOS plot does not automatically give you catalytic activity or conductivity.

This chapter continues the silicon primitive cell from the previous chapter. finish first final static self-consistent evaluation , and then Complete copy of POSCAR, POTCAR and corresponding CHGCAR to the new properties directory. You can also practice the procedure on a fixed initial 5.43 Å structure, but the results must be labeled "calculations for this fixed geometry" and cannot be called optimized properties.

Occupation patterns and electron smearing affect points, energy, and the look and feel of DOS. The following table is a selection guide and is not a universal formula for all materials. In particular, just because a certain metal template uses ISMEAR=1, it cannot be directly attached to the insulator.

Smearing can make numerical iterations smooth, but excessive smearing will also drag the density of states into the band gap. Increasing NEDOS only makes the energy axis sample finer and will not make up for the lack of k-points. For metals, increasing smearing can usually reduce some numerical difficulties, but it cannot be a reason to obscure the details of the real Fermi surface. You are ultimately responsible for the physical quantities you want to report, not just for "the curve looks good."

Check basis: ISMEAR value · smearing method and scope of application · NEDOS

Preparation: New si-dos , copy the converged statically calculated POSCAR, POTCAR, CHGCAR. Keeping the same structure, PBE approximation and cutoff energy. used here ICHARG=1 , that is, reading the existing density and making a first guess, but continuing to be self-consistent on a denser grid; it is different from ICHARG=11, which fixes the charge density.

The following grid 16³, ENCUT=520 eV and NBANDS=16 are the starting point for this set of two-atom silicon learning inputs and do not mean that convergence has been done. If a higher cutoff energy has been selected in the previous chapter, it must be replaced simultaneously. POTCAR remains unchanged, especially CHGCAR corresponding to another Si potential cannot be mixed in.

SYSTEM = Si PBE dense-mesh DOS
GGA = PE
ENCUT = 520
PREC = Accurate
LREAL = .FALSE.
ISPIN = 1
ISTART = 0
ICHARG = 1
ALGO = Normal
EDIFF = 1E-7
NELM = 120
ISMEAR = -5
NBANDS = 16
LORBIT = 11
NEDOS = 2001
IBRION = -1
NSW = 0
LCHARG = .TRUE.
LWAVE = .TRUE.

Si dense uniform mesh for DOS
0
Gamma
16 16 16
0 0 0

Output: See DOSCAR, and vasprun.xml or vaspout.h5 when the corresponding output is enabled. When drawing a graph, the horizontal axis is usually changed to E−EF, or the top of the semiconductor valence band can be set to zero; no matter which benchmark is used, it must be stated in the legend. For systems with a band gap, the Fermi energy given by the program has a conventional dependence on the band gap. Do not directly interpret its small movements as charge transfer.

import py4vasp
calc = py4vasp.Calculation.from_path("si-dos")
calc.dos.plot()
# Requires LORBIT output and calculation results that support HDF5
calc.dos.plot("Si(s), Si(p)")

This code requires use in a Jupyter environment that has py4vasp-compatible installed and that the computational output supported by that tool actually exists. please si-dos Start the notebook in the parent directory of the notebook; if the notebook is already in the results directory si-dos inside, then put from_path("si-dos") Change to from_path(".") , to avoid looking for non-existent si-dos/si-dos . It reads the results of your runs and does not automatically execute VASP. If there is no HDF5 output, a tool that can read vasprun.xml or DOSCAR can be used, but the tool version, projection column meaning and energy zero point should be recorded, and the column order of different spin modes cannot be blindly guessed.

How to confirm that the results are available: SCF converges; DOS tends to be stable as the k-grid continues to be refined; unoccupied bands cover the energy range studied; the number of integrated states is compatible with the electron count of the potential used; the total DOS and projection weight changes can be explained; the upper and lower curves in the spin resolution diagram have a clear agreement. Spin-down DOS drawn on the negative semi-axis is often just a drawing convention and does not represent "negative density of states".

Common misreadings: The overlap of two peaks can only provide clues to orbital hybridization and is not sufficient to prove a chemical bond alone; high DOS near the Fermi energy does not automatically mean high conductivity, and also involves speed, scattering and carrier lifetime. Testing an explanation with additional evidence is more reliable than making up a story for each peak.

Check basis: ICHARG · LORBIT projection definition · unoccupied bands number NBANDS · py4vasp DOS

The uniform k-grid takes care of the electron density and Brillouin zone integrals, and the band structure path shows the eigenvalue changes in a specific direction. Taking many points along a few lines of high symmetry is no substitute for a self-consistent grid covering the entire Brillouin zone. So you can't start from an empty folder, just put a Line-mode KPOINTS and expect to get a trusted band structure.

1. Take uniform mesh SCF or DOS results for the same structure that has converged.
2. create si-bands , copy its POSCAR, POTCAR, CHGCAR.
3. Keeping PBE, ENCUT, magnetism and geometry consistent, instead use the following fixed density INCAR and path KPOINTS.
4. Use vasp_std ;After completion, read the coordinates of the eigenvalues, projections and matching paths and check the convergence of the band structure of interest.

SYSTEM = Si PBE bands from converged density
GGA = PE
ENCUT = 520
PREC = Accurate
LREAL = .FALSE.
ISPIN = 1
ISTART = 0
ICHARG = 11
ALGO = Normal
EDIFF = 1E-7
NELM = 120
ISMEAR = 0
SIGMA = 0.05
NBANDS = 16
LORBIT = 11
IBRION = -1
NSW = 0
LCHARG = .FALSE.
LWAVE = .FALSE.

ICHARG=11 By fixing the density input to CHGCAR, the electron solution still needs to converge, but it no longer searches for a new self-consistent density corresponding to the waypoint. ISTART=0 Avoid misreading old wavefunctions from another set of k-grids. LCHARG=.FALSE. and LWAVE=.FALSE. Reduce unnecessary output; original copies of input CHGCAR should be retained.

The following is the path of the silicon primitive cell in this tutorial: L–Γ–X–U and K–Γ two segments, U and K are disconnected. The primitive cell vectors must be arranged in the order of the previous chapter, that is, a₁=(0,a/2,a/2), a₂=(a/2,0,a/2), a₃=(a/2,a/2,0). The path coordinates are based on this set of three reciprocal lattice vectors; the overall translation of the origin does not change the path, but changing the unit cell base will change its fractional coordinates.

Diamond Si primitive: L-G-X-U | K-G
40
Line-mode
Reciprocal
0.500000 0.500000 0.500000 L
0.000000 0.000000 0.000000 G

0.000000 0.000000 0.000000 G
0.500000 0.000000 0.500000 X

0.500000 0.000000 0.500000 X
0.625000 0.250000 0.625000 U

0.375000 0.375000 0.750000 K
0.000000 0.000000 0.000000 G

40 is the number of samples per segment including the endpoint. Increasing it makes the curves denser but does not improve the integration accuracy of the previous CHGCAR. Here is a set of teaching display directions. It cannot be guaranteed that the extreme values of the band structure of all materials are on the high symmetry line. When reliable ground state band gaps are required, the valence band top and conduction band bottom should also be located in sufficiently dense three-dimensional sampling.

import py4vasp
scf = py4vasp.Calculation.from_path("si-static")
ef = scf.dos.read()["fermi_energy"]
calc = py4vasp.Calculation.from_path("si-bands")
calc.band.plot(fermi_energy=ef)

in si-static and si-bands Run this Jupyter code in the common parent directory of the fermi_energy Parameters compatible with py4vasp. Both directories must save the data generated by this actual calculation and can be read by this version. vaspout.h5 ;This requires VASP's HDF5 support. It is assumed here that the CHGCAR used by the band structure comes from si-static ;if from the previous section si-dos , and synchronize the reference directory of the code to si-dos . Structure, functional and potential must match. The code extracts EF from the self-consistent result of the uniform grid and then explicitly gives it to the band structure drawing tool. Therefore, this reference EF has been subtracted from the figure and cannot be subtracted manually again. When no HDF5 output is available, EFs should be recorded from the corresponding self-consistent OUTCAR and other plotting tools that support raw eigenvalues and explicit energy alignment should be used. This example cannot be executed directly.

Vertical axis and band gap: Do not consider the re-estimated EF in non-self-consistent path calculations unconditionally as a reliable energy reference, since these path points are not normal volume integral meshes. One should adopt the same structure self-consistent calculated reference or explicit valence band top alignment strategy and verify that the post-processing tool has automatically subtracted EF. Repeatedly subtracting once will cause the overall energy level to become misaligned. When reporting the band gap, indicate functional, SOC or not, structure, sampling method and direct/indirect nature.

Semi-local functionals such as PBE often cannot give accurate experimental basic band gaps; ordinary Kohn–Sham band structures are also not complete optical excitation spectra. It is observed that the calculation is different from the experiment. First identify whether the comparison is the basic band gap, optical band gap or spectral peak at a certain temperature, and then choose a more appropriate method.

Checkpoints: The path labels match the actual unit cell; the eigenvalues of the adjacent segments are continuous at the same endpoint; U–K are not mistakenly connected into the calculated segment; the unoccupied bands are sufficient and the high-band structure of interest is converged; the path is regenerated after the structure undergoes symmetry changes. Defects Supercell will have band structure folding, and you cannot simply use the folded dense curves to correspond to primitive cells one by one.

Check basis: DFT band structure workflow · Official diamond silicon primitive cell and path · py4vasp band structure · hybrid functional band structure

Nonmagnetic calculations constrain the two spin channels to be equivalent; collinear spin calculations allow them to have different electron densities. Models containing Fe, Co, Ni, certain defects, or open-shell molecules cannot be considered non-magnetic just because the default settings work. There may be multiple local solutions to the self-consistent equation, and the initialization method will affect which one ends up.

Use the two-atom conventional bcc Fe cell below to explore how the initial spin affects the solution. The specified 2.87 Å geometry is a teaching model, not a structure optimized in this tutorial. The exercise compares candidate magnetic states at the same fixed geometry; it does not cover all possible magnetic orders. Use one licensed PAW-PBE Fe dataset. Revalidate the 520 eV cutoff and 12³ mesh for Fe rather than reusing convergence conclusions for Si.

bcc Fe conventional cell; teaching geometry
2.8700000000
1.0000000000 0.0000000000 0.0000000000
0.0000000000 1.0000000000 0.0000000000
0.0000000000 0.0000000000 1.0000000000
Fe
2
Direct
0.0000000000 0.0000000000 0.0000000000
0.5000000000 0.5000000000 0.5000000000

bcc Fe conventional cell uniform mesh
0
Gamma
12 12 12
0 0 0

SYSTEM = bcc Fe FM initial-state test
GGA = PE
ENCUT = 520
PREC = Accurate
LREAL = .FALSE.
LASPH = .TRUE.
ISPIN = 2
MAGMOM = 2.5 2.5
ISYM = 0
ISTART = 0
ICHARG = 2
ALGO = Normal
EDIFF = 1E-7
NELM = 160
ISMEAR = 1
SIGMA = 0.15
LORBIT = 11
IBRION = -1
NSW = 0
LCHARG = .TRUE.
LWAVE = .TRUE.

Create another one for the reverse initial state fe-afm directory, use exactly the same POSCAR, POTCAR, KPOINTS, change the magnetic moment line in INCAR to the following settings, still from ISTART=0; ICHARG=2 Start. Do not copy the CHGCAR or old WAVECAR containing the previous magnetization density at this time, otherwise the "changed initial magnetic moment" may not be re-initialized as you expected.

MAGMOM = 2.5 -2.5

How to compare: Confirm that each candidate state SCF has converged, record the final total magnetic moment, the local magnetic moment of each atom and the uniformly defined energy, and see which one has the lower energy under the same model and numerical accuracy. If two different initial states end up with the same solution, this fact should also be recorded. If we compare the true equilibrium magnetic states, we usually need to allow the corresponding magnetic states to optimize the structure separately, and then make the final static comparison with uniform accuracy.

The total magnetic moment approaching zero does not alone prove nonmagnetism: the two opposing local magnetic moments can cancel. To distinguish between nonmagnetic, ferromagnetic, antiferromagnetic and ferrimagnetic, one must look at the atomic resolution magnetic moment and spatial arrangement at the same time. The initial primitive cell must also accommodate the envisaged magnetic period; a unit cell that is too small cannot "grow" an antiferromagnetic mode that it does not allow via an INCAR parameter.

When restarting the magnetization task, even though MAGMOM is no longer responsible for resetting the magnetization guess, it still participates in the symmetry determination. Arbitrary deletions may change the allowed symmetries. When switching initial state studies, a new directory, initialization from scratch, and complete recording are more reliable than changing only one row of magnetic moments. Magnetic anisotropy, non-collinearity and SOC are taken to the next stage and numerical accuracy is improved to the level required for target energy differences.

Check basis: ISPIN · MAGMOM and restart · The difference between ISYM · LASPH

The localized d/f electrons of some transition metals or rare earth compounds are underdescribed by ordinary semilocal functionals. DFT+U adds corrections to selected local orbital subspaces, but its effect depends on the subspace, potential file, functional, oxidation state, and U/J definition used. Simply moving the U on someone else's material to the new system will not automatically lead to more accurate results.

The following is Parameter template to be configured : Assume that you have been provided with a unit cell of Ni₂O₂ capable of accommodating the magnetic order under study, the POSCAR elemental order is Ni O, the number is 2 2, and the first two atoms are Ni. It provides no structure and does not fill in U; so it cannot be run directly. Determine U first based on a scheme such as the original paper or consistent linear response, save its source, and then replace the placeholders with numbers.

# TEMPLATE ONLY: supply a validated Ni2O2 magnetic cell and choose U_eff.
# POSCAR species must be Ni O with counts 2 2 for the MAGMOM below.
SYSTEM = Ni2O2 PBE+U magnetic template
GGA = PE
ENCUT = YOUR_CONVERGED_ENCUT
PREC = Accurate
LREAL = .FALSE.
LASPH = .TRUE.
ISPIN = 2
MAGMOM = 2.0 -2.0 0.0 0.0
ISYM = 0
ISTART = 0
ICHARG = 2
ALGO = Normal
EDIFF = 1E-7
NELM = 160
ISMEAR = 0
SIGMA = 0.05
LDAU = .TRUE.
LDAUTYPE = 2
LDAUL = 2 -1
LDAUU = YOUR_U_EFF 0.0
LDAUJ = 0.0 0.0
LMAXMIX = 4
LORBIT = 11
IBRION = -1
NSW = 0
LCHARG = .TRUE.
LWAVE = .TRUE.

The array length must be equal to the number of element types: two numbers for this Ni O example, and four initial values for four atoms for MAGMOM. If O Ni order is used, LDAUL, LDAUU, LDAUJ, and POTCAR are all put into the new order together. Different variants of the same element in the potential file have different valence electron treatments and cannot be switched at will in the same set of controls.

When receiving DOS/band structure: First obtain a self-consistent CHGCAR with the complete DFT+U model. Set the appropriate LMAXMIX in both the initial SCF calculation and the subsequent fixed-density calculation: commonly 4 for d elements and 6 for f elements. Keep the same LDAU parameters, magnetic state, and structure, then choose ICHARG=1 for self-consistent DOS or, where applicable, ICHARG=11 for a band-path calculation. Adding the missing setting only during post-processing cannot recover occupation-matrix information that was never saved.

How to justify your choice: Report the U/J form, numerical values, orbitals, sources and applicable chemical environments used, and examine whether the core conclusions are sensitive to reasonable parameter intervals. If a change in U results in a change in magnetic state, bonding, or stability ordering, this is a result that needs discussion. Just because a certain U has the lowest calculated absolute total energy, it cannot be said to be the "best U"; different Us correspond to different energy functionals and cannot be compared like this.

DFT+U is often accompanied by local minima of multiple magnetic states and orbital occupancies. It is necessary to test different reasonable initial states and check the final local magnetic moment, occupation, energy gap and structure; numerical convergence only means that a certain solution has been found, and does not guarantee that it is the most relevant physical solution. For metals, covalent systems and the coexistence of multiple valence states, more careful arguments are needed instead of treating U as a unified band gap corrector.

Check basis: DFT+U form · LDAUL · LDAUU · LDAUJ · LMAXMIX · LASPH

Spin–orbit coupling links the spin of an electron to its orbital motion. For heavy elements, band structure splitting, certain topological and magnetic anisotropy problems, it may change key conclusions. However, "Open LSORBIT" only opens the model terms and is far from completing physical interpretation or accuracy verification.

This example continues to use Nonmagnetic primitive cell of two Si atoms . create si-soc , copying the previously statically calculated POSCAR, POTCAR, CHGCAR, and sufficiently converged uniform KPOINTS for non-magnetic PBE, use instead vasp_ncl . Here ICHARG=1 re-self-consistent SOC at the existing density, instead of just reading a few energies at the old density.

If you copy the Slurm script of the introductory chapter, you must also command -v vasp_std with srun vasp_std The program name in the two lines is changed to vasp_ncl ; Changing only INCAR will not change the actually started program. Save the verified script in si-soc , enter the directory and then submit.

SYSTEM = Si PBE SOC self-consistent example
GGA = PE
ENCUT = 520
PREC = Accurate
LREAL = .FALSE.
LASPH = .TRUE.
GGA_COMPAT = .FALSE.
ISTART = 0
ICHARG = 1
ALGO = Normal
EDIFF = 1E-7
NELM = 160
ISMEAR = 0
SIGMA = 0.05
LSORBIT = .TRUE.
LNONCOLLINEAR = .TRUE.
SAXIS = 0 0 1
MAGMOM = 0 0 0  0 0 0
ISYM = -1
NBANDS = 32
LORBIT = 11
IBRION = -1
NSW = 0
LCHARG = .TRUE.
LWAVE = .TRUE.

Both atoms start with zero magnetic moments in this nonmagnetic Si SOC example. For a magnetic material, specify three components per atom and interpret their directions in the spin coordinate system defined by SAXIS. Components reported for an arbitrary SAXIS are not automatically laboratory-frame Cartesian components. When restarting from an existing magnetization density, changing MAGMOM alone does not arbitrarily rotate the stored moments on selected atoms.

SOC band structure How to proceed: First validate the SOC SCF result on the uniform mesh. Then use its CHGCAR and the same SOC settings for an applicable fixed-density PBE band calculation. Retain vasp_ncl, three-component MAGMOM, SAXIS, symmetry settings, and projection conventions. Reading a CHGCAR obtained without SOC does not make the resulting calculation a self-consistent SOC ground state.

Success check: Confirm that the ncl executable is used and that SOC and noncollinear settings are active in OUTCAR. Check electronic convergence, the stability of the splitting with denser meshes and more unoccupied bands, and consistent structure and energy alignment relative to the calculation without SOC. For precise energy differences, rerun the no-SOC static reference self-consistently with the same LASPH and GGA_COMPAT settings so that numerical changes are not mistaken for a pure SOC effect. Magnetic anisotropy energies are often far smaller than ordinary bond energies and need stricter mesh, electronic-threshold, and symmetry tests. The teaching precision on this page is not sufficient evidence for quantitative conclusions.

Opening the SOC and seeing two curves separate does not alone prove that the topology is nontrivial; related judgments may also require symmetry, topological invariants, surface states, and the stability of the approximation used. First, accurately describe "what band structure changes occur under this model", and then discuss its physical meaning.

Check basis: LSORBIT · LNONCOLLINEAR and Version Notes · SAXIS coordinate system · NBANDS and Spinors

Numerical accuracy supplement: GGA_COMPAT official definition .

The overall density of states and band structure describe electrons in energy and momentum space, while charge density puts the problem back into real space. First make it clear whether you want to answer "where are the electrons rearranged" or "how many electrons are allocated to each atom according to a certain spatial division". The two are related, but they are not the same quantity, nor can a certain partition charge be directly equated to an integer oxidation state.

On the confirmed Si static structure, create si-charge , following the same POSCAR, POTCAR and convergence KPOINTS, using the following full INCAR. This step does static self-consistent from scratch and outputs reconstructed densities, without using unvalidated densities from the end of structural optimization or molecular dynamics.

SYSTEM = Si static charge analysis example
GGA = PE
ENCUT = 520
PREC = Accurate
LREAL = .FALSE.
ISPIN = 1
ISTART = 0
ICHARG = 2
ALGO = Normal
EDIFF = 1E-7
NELM = 120
ISMEAR = 0
SIGMA = 0.05
IBRION = -1
NSW = 0
LCHARG = .TRUE.
LWAVE = .FALSE.
LAECHG = .TRUE.

LAECHG=.TRUE. Generate AECCAR0, AECCAR1, and AECCAR2, which involve nuclear density, atomic superposition valence density, and self-consistent valence density respectively. What is used for common Bader reference partitioning is the sum of AECCAR0 and AECCAR2, not the sum of all three files. You must first confirm that self-consistent is complete, that these files belong to the same final static job, and perform a convergence check on the fine FFT mesh.

After installing and documenting versions of the Henkelman Group Bader tool and its accompanying scripts, you can run common processes provided by the authors in this results directory:

chgsum.pl AECCAR0 AECCAR2
bader CHGCAR -ref CHGCAR_sum

This process uses the reconstructed reference density to determine spatial partitioning and then integrates the valence electron density given by CHGCAR. The results are usually in ACF.dat. If the number of valence electrons assigned to an atom is recorded as N _Bader , a convention for net atomic charge is:

Under this sign convention, positive values represent missing electrons relative to the potential file valence electron count. ZVAL must be taken from the actual POTCAR, and the atomic number cannot be substituted casually; if another definition of "getting electrons as positive" is used, the inverse sign must be clearly written. Check whether the sum of the number of electrons of all atoms, the charge in the vacuum area, and the total number of electrons are self-consistent, and gradually increase NGXF/NGYF/NGZF to see whether the partition number of each atom is stable.

Ideal symmetrical crystals of the same element can be used as a consistency check, but systems with interfaces, vacancies, magnetism, and vacuum are more complex. The non-integer charge given by Bader depends on the partitioning and numerical grid; it supports electron rearrangement analysis and cannot be simply interpreted as "this atom is strictly +1.37 valence".

Charge-density differences require compatible calculations

For the combination of adsorbate A and surface B, a commonly used definition is Δρ=ρ _A+B−ρ_A−ρ_B . The reference calculations for A and B retain their positions in the combined structure, removing another portion of the atoms but not relaxing them separately again; otherwise density differences due to geometric changes are mixed into the plot. The total system and the fragment system must use the same unit cell, grid origin, fine FFT size, compatibility potential and functional, the same cutoff energy, and the number of electrons and charge states must meet the reference process you define.

It is recommended that the real-space grid be read by a post-processor that understands the VASP file format and charge conservation checked with full-space integration, rather than line-by-line subtraction of two CHGCAR files. For a neutral segment with a consistent definition of the total number of electrons, the whole cell integral of Δρ should be close to zero within the numerical error; the charged reference needs to be re-judged according to its definition. When drawing, check several isosurface thresholds at the same time to avoid choosing only the one graph that best supports the preset conclusion.

Complete the exercise: Write a legend for a DOS, band structure, or differential charge diagram that includes input sources, coordinate/energy bases, numerical convergence, and interpretation boundaries. Then list at least one additional check that would overturn the current interpretation. enter Advanced tasks such as adsorption, phonon, and molecular dynamics , maintain the same standard of evidence.

Check basis: LAECHG output · Bader tool author's note · CHGCAR data meaning · Fine FFT mesh NGXF · Number of electrons NELECT

Read the file description and corresponding chapters first, and then complete the model, POTCAR and pre-calculation. Files with the suffix .template contain parameters to be changed and cannot be submitted as they are.

The main text is an independently organized study instruction. For label meanings and version differences, refer to the official documentation corresponding to your program version.

1. ISMEAR
2. Smearing technique
3. ICHARG
4. LORBIT
5. NEDOS
6. NBANDS
7. py4vasp DOS
8. DFT band structure process
9. Official Diamond Silicon Tutorial
10. py4vasp band structure
11. hybrid functional band structure
12. ISPIN
13. MAGMOM
14. ISYM
15. LASPH
16. LDAUTYPE
17. LDAUL
18. LDAUU
19. LDAUJ
20. LMAXMIX
21. LSORBIT
22. LNONCOLLINEAR
23. SAXIS
24. LAECHG
25. Bader tool author's note
26. CHGCAR
27. NGXF
28. NELECT
29. GGA_COMPAT