1.1. Silicon Supercell

A calculation is said to be numerically converged when further refinement of these parameters produces changes in the target property smaller than a chosen tolerance (e.g., \(<\) 1 meV/atom in total energy, \(<\) 0.01 eV in a band gap). This is distinct from physical convergence, which refers to reaching a stationary state with respect to structural degrees of freedom (e.g., geometry optimization or relaxation). [PAYNE1992]

In this tutorial we’ll show how to perform a numerical convergence of real-space grid resolution and Brillouin-zone sampling density in a silicon supercell ( Fig. 1.1.1) using RESCU code.

Fig. 1.1.1 Conventional silicon \(2\times 2\times 2\) supercell visualization.

The calculation steps are given as follows:

1. Converge the real-space grid first

   • Fix k-mesh to a small value (e.g., \(2\times 2\times 2\)).

   • Test from coarse to fine domain.lowres values.

   • Select the coarsest grid meeting the \(\Delta E\) criterion (e.g., \(\leq 10^{-4}\) eV/atom).

2. Converge the k-points

   • Fix domain.lowres to the converged value.

   • Increase k-mesh density until \(\Delta E\leq 10^{-4}\) eV/atom.

3. Finalize SCF settings

   • Total energy tolerance for production runs;

   • Residual tolerance;

   • Adjust mixing parameters if oscillations occur.

1.1.1. Project Layout

To keep your workspace organized, we recommend creating a project folder with three subdirectories:

numerical_convergence/
├── generate_bulk_structures_for_convergence.py
├── xyz/
│   └── Si_2x2x2.xyz
├── lowres/
│   ├── scf_base.input
│   ├── run_lowres_convergence.py
│   └── plot_lowres_energy.py
└── kpoints/
    ├── scf_base.input
    ├── run_kpoint_convergence.py
    └── plot_kpoints_energy.py

Note

• After downloading, place the scripts into the subfolders shown above (generate_* at project root; run_*/plot_* and template *.input files into lowres/ and kpoints/ respectively).

• The plotting scripts write PNGs to the current directory (no dedicated figs/ folder is required on your machine).

1.1.2. SCF Convergence Criteria

The self-consistent field (SCF) procedure solves the Kohn–Sham equations iteratively, updating the electron density (or Kohn–Sham potential) until self-consistency is reached. The accuracy of the final result depends on the convergence criteria applied during this process. The total energy tolerance sets the maximum allowed change in the total Kohn–Sham energy between consecutive SCF iterations:

\[\Delta E_{\mathrm{tot}}^{(n)} = E_{\mathrm{tot}}^{(n)} - E_{\mathrm{tot}}^{(n-1)},\]

A typical production threshold for ground-state calculations is on the order of \(10^{-6}\) Ha (\(\approx 2.7\times 10^{-5}\) eV) per atom, which is sufficient for most energy differences to be meaningful at the meV scale. Looser thresholds may be acceptable for rough preliminary runs; tighter thresholds are required for sensitive quantities such as phonon frequencies or small energy differences between competing phases.

In RESCU, SCF convergence is controlled by a two-component vector: the charge variation per valence electron (dρ) and the total-energy variation per valence electron (dE). Convergence based on these quantities is often more robust than relying on total-energy change alone, since the latter can be small even when the density is not fully converged (e.g., in metallic systems). [PAYNE1992], [KRESSE1996]

1.1.3. Real-Space Grid Convergence

The real-space grid spacing (distance between grid points) used to represent wavefunctions and potentials is set in RESCU by the domain.lowres parameter. Smaller values correspond to finer grids (higher spatial resolution), which generally improve accuracy but increase memory usage and computational time. A well-converged domain.lowres ensures that total energies, forces, and band edges are not affected by the discretization of the real-space grid. In this tutorial we choose the largest domain.lowres for which \(\Delta E \le 10^{-4}\) eV/atom. General guidance on real-space grid (mesh) convergence can be found in standard documentation for real-space and localized-basis codes. [SIESTA_GRID]

To avoid coupled errors we’ll start converging the domain.lowres before k-points. [KRATZER2019]

Step-by-step procedure

1.1.3.1. Generate the bulk structure

To begin we’ll generate the silicon unit cell (\(2\times 2\times 2\)) using lattice parameter of a = 5.43 Å.

Download: generate_bulk_structures_for_convergence.py and place it at the project root (numerical_convergence/), then run:

cd numerical_convergence
python generate_bulk_structures_for_convergence.py

This script generates the following file: xyz/Si_2x2x2.xyz

Download: Si_2x2x2.xyz

1.1.3.2. Prepare base input file

1. Create the folder

   numerical_convergence/lowres/
   

2. Add the file (exact name)

   scf_base.input
   

3. Use these placeholders (ALL CAPS) so the script can replace them for each run

   ``LOWRES``, numeric resolution parameter we'll sweep
   
   ``TAG``, short label for the case/run.
   

   SCF input file

   Listing 1.1.1 lowres/scf_base.input (real-space grid sweep)

   units.length = 'Angstrom';
   LCAO.status = true;
   info.calculationType = 'self-consistent';
   info.savepath = './results/scf_TAG';
   info.outfile = './results/scf_TAG.out';
   
   % Lattice vectors (conventional cubic cell for Si)
   a = 5.43;  % Experimental lattice constant for Si in Angstrom
   domain.latvec = [a*2, 0, 0; 0, a*2, 0; 0, 0, a*2];
   units.domain.latvec = 'Angstrom';
   
   % Atom positions (diamond structure for Si)
   atom.xyz = '../xyz/Si_2x2x2.xyz';
   units.atom.xyz = 'Angstrom';
   
   % Pseudopotential for Si
   element(1).species = 'Si';
   element(1).path = './Si_DZP_PBE.mat';
   
   % Functional (GGA PBE)
   functional.list = {'XC_GGA_X_PBE', 'XC_GGA_C_PBE'};
   functional.libxc = true;
   
   % k-points
   kpoint.gridn = [2,2,2];
   kpoint.shift = [0.5, 0.5, 0.5];
   kpoint.sampling = 'gauss';
   
   % Numerical settings
   domain.lowres = LOWRES;   % <-- Placeholder for lowres value
   eigensolver.emptyBand = 8;
   
   % SCF settings
   option.maxSCFiteration = 200;
   spin.type = 'degenerate';
   

   The keywords specify the following:

   • LCAO.status If LCAO.status = true, RESCU performs a calculation using the atomic orbital basis. This is generally faster than performing a purely real space calculation, but the atomic orbital basis must be validated for correctness.

   • info.calculationType This input parameter specifies the type of the calculation that RESCU should perform.

   • info.savepath The path and name prefix of the file where the results will be written. RESCU appends the rank number to the file name prefix, which for serial calculations is always 0.

   • info.outfile File in which writes the simulation information.

   • domain.latvec Lattice vectors. The matrix must contain three row-vectors.

   • units.domain.latvec and units.atom.xyz Specify units as ‘Angstrom’ for readability.

   • atom.xyz Is an array of row-vectors containing the position of each atom.

   • element(1).species and element(1).path Tag and path to the pseudopotential file for each element.

   • functional.list List of exchange-correlation functionals used.

   • functional.libxc Set to true to use LibXC for more functionals.

   • kpoint.gridn Sets the number of k-points used to sample the first Brillouin zone.

   • kpoint.shift Shift the k-point grid uniformly. Monkhorst-Pack grids are gamma-centered for odd grids and shifted by [0.5,0.5,0.5] for even grids.

   • domain.lowres Controls the real space resolution (finer grid = smaller value = more accurate but slower).

   • eigensolver.emptyBand Number of bands included in the calculation in addition to the number of valence bands.

   • option.maxSCFiteration Maximal number of self-consistent steps.

   • spin.type Set to ‘degenerate’ for non-spin-polarized calculations.

4. Now, just run the SCF calculation by using the following script. (Note: please simplify the script just for SCF test)

   Download: run_lowres_convergence.py. Place it in numerical_convergence/lowres/ directory and run:

   cd numerical_convergence/lowres
   python run_lowres_convergence.py
   

1.1.3.3. Analyze energy results

After the calculations, the script below could be used to extract the following information:

• A printed table of ΔE (relative to the finest grid) for each domain.lowres.

• Plot a figure (energy_vs_lowres.png) of total energy vs. grid spacing (saved in the current directory);

To extract the energy results Download: plot_lowres_energy.py. Place it in numerical_convergence/lowres/ and run:

cd numerical_convergence/lowres
python plot_lowres_energy.py

From the printed table and plot, you can check how the convergence of the total energy per atom changes as the real-space grid spacing (lowres) is varied.

Example energy output:

Convergence check relative to finest grid = 0.20 Å

  lowres=0.20 ⇒ ΔE = 0.00e+00 eV/atom  → converged (tol=1.0e-04)
  lowres=0.25 ⇒ ΔE = -4.55e-07 eV/atom  → converged (tol=1.0e-04)
  lowres=0.30 ⇒ ΔE = -3.74e-06 eV/atom  → converged (tol=1.0e-04)
  lowres=0.40 ⇒ ΔE = -8.92e-06 eV/atom  → converged (tol=1.0e-04)
  lowres=0.50 ⇒ ΔE = -1.36e-04 eV/atom  → NOT converged (tol=1.0e-04)
  lowres=0.60 ⇒ ΔE = -2.34e-03 eV/atom  → NOT converged (tol=1.0e-04)

Fig. 1.1.2 Example: total energy convergence vs. domain.lowres for silicon unit cell \(2\times 2\times 2\) .

The finest grid is 0.20 Å. Every other row reports how much the energy changes relative to this reference. For lowres = 0.25 … 0.40, the ΔE \(< 10^{-4}\) eV/atom. Hence these rows are converged, while 0.50 and 0.60 are not. You can pick a more economical grid (e.g., 0.30–0.40) and still stay within the tolerance.

1.1.4. k-Point Convergence

The k-point mesh controls the sampling of the Brillouin zone in periodic systems. Insufficient k-point sampling can cause significant errors in total energies, forces, and band structures, especially in metals and systems with small Brillouin zones.

Once the real-space grid is converged, fix domain.lowres with the value within the tolerance and test \(k_{1}, k_{2}, k_{3}\) Monkhorst-Pack meshes. Here we’ll use k = 1, …,8.

Tip

For even k, set kpoint.shift = [0.5, 0.5, 0.5] to avoid :math:`Gamma-`point bias. For HSE calculations set the Γ-center grid is mandatory; for PBE it is optional.

Step-by-step procedure

1.1.4.1. Bulk structure

1. Ensure bulk structure is available

   Reuse the xyz/ directory from the real-space grid step or generate the structure using the script,

generate_bulk_structures_for_convergence.py, and run at project root.

1.1.4.2. Prepare base input file

1. Create the folder

   numerical_convergence/kpoints/
   

2. Add the file (exact name)

   scf_base.input
   

3. Use these placeholders (ALL CAPS) so the script can replace them for each run.

   ``KGRID``, numeric *k-*points parameter we'll sweep;
   
   ``SHIFTLINE``, use kpoint.shift = [0.5,0.5,0.5] for even grids, otherwise remove the line
   
   ``TAG``, short label for the case/run.
   

   SCF input file

   Listing 1.1.2 kpoints/scf_base.input (k-point sweep)

   units.length = 'Angstrom';
   LCAO.status = true;
   info.calculationType = 'self-consistent';
   info.savepath = './results/scf_TAG';  % This will be replaced below
   info.outfile = './results/scf_TAG.out';
   
   % Lattice vectors (conventional cubic cell for Si)
   a = 5.43;  % Experimental lattice constant for Si in Angstrom
   domain.latvec = [a*2, 0, 0; 0, a*2, 0; 0, 0, a*2];
   units.domain.latvec = 'Angstrom';
   
   % Atom positions (diamond structure for Si)
   atom.xyz = '../xyz/Si_2x2x2.xyz';
   units.atom.xyz = 'Angstrom';
   
   % Pseudopotential for Si
   element(1).species = 'Si';
   element(1).path = './Si_DZP_PBE.mat';
   
   % Functional (GGA PBE)
   functional.list = {'XC_GGA_X_PBE', 'XC_GGA_C_PBE'};
   functional.libxc = true;
   
   % k-points (using a fixed grid as already tested)
   kpoint.gridn = [KGRID,KGRID,KGRID];
   SHIFTLINE % Γ-centered grid: mandatory for HSE; optional for PBE (kept consistent here)
   
   % Numerical settings
   domain.lowres = 0.2;
   eigensolver.emptyBand = 8;
   
   % SCF settings
   option.maxSCFiteration = 200;
   spin.type = 'degenerate';
   

4. Now, just run the SCF calculation by using the following script.

   Download: run_kpoint_convergence.py. Place it in numerical_convergence/kpoints/ directory and run:

   cd numerical_convergence/kpoints
   python run_kpoint_convergence.py
   

The script runs SCF for each k-point mesh, applying a shift for even grids.

1.1.4.3. Analyze k-point results

After the calculations, the script below could be used to extract the following information:

• A printed table of ΔE (relative to the largest k-mesh) for each k-point.

• Plot a figure (energy_vs_kpoints.png) of total energy vs. k-mesh size (saved in the current directory)

To extract the energy results, download: plot_kpoints_energy.py. Place it in numerical_convergence/kpoints/ and run:

cd numerical_convergence/kpoints
python plot_kpoints_energy.py

The script will compare total energies per atom against the reference (largest k-point mesh) and report whether the calculation is converged within the specified tolerance.

Example output:

Convergence check relative to largest mesh = 8×8×8
  k= 1 ⇒ ΔE = 1.08e-01 eV/atom  → NOT converged (tol=1.0e-04)
  k= 2 ⇒ ΔE = 2.65e-03 eV/atom  → NOT converged (tol=1.0e-04)
  k= 3 ⇒ ΔE = 1.16e-04 eV/atom  → NOT converged (tol=1.0e-04)
  k= 4 ⇒ ΔE = 7.26e-06 eV/atom  → converged (tol=1.0e-04)
  k= 5 ⇒ ΔE = 6.17e-07 eV/atom  → converged (tol=1.0e-04)
  k= 6 ⇒ ΔE = 6.93e-08 eV/atom  → converged (tol=1.0e-04)
  k= 7 ⇒ ΔE = 8.65e-09 eV/atom  → converged (tol=1.0e-04)
  k= 8 ⇒ ΔE = 0.00e+00 eV/atom  → converged (tol=1.0e-04)

Fig. 1.1.3 Example: total energy convergence vs. k-point mesh for Si.

The largest mesh is \(8\times 8\times 8\) (reference). For \(k\ge 4\), the energy differences satisfy \(|\Delta E| < 10^{-4}\) eV/atom; e.g., \(5\times 5\times 5\) is already within tolerance and may be chosen as an economical mesh for production. [MONKHORST1976]

Tip

Always converge k-points after grid spacing. Doing both at once can couple errors and mislead threshold selection. [KRATZER2019]

For metals, combine a finer k-mesh with a finite electronic smearing (e.g., Methfessel–Paxton or Marzari–Vanderbilt “cold” smearing) and converge both with respect to your property of interest. [METHFESSEL1989], [MARZARI1999]

1.1.5. References

[MARTIN]

R. M. Martin. Electronic Structure: Basic Theory and Practical Methods, 2nd ed. Cambridge University Press (2020).

[KRATZER2019] (1,2)

P. Kratzer. The Basics of Electronic Structure Theory for Periodic Systems Frontiers in Chemistry 7, 106 (2019).

[PAYNE1992] (1,2)

M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias, and J. D. Joannopoulos. Iterative minimization techniques for ab initio total-energy calculations: molecular dynamics and conjugate gradients. Rev. Mod. Phys. 64, 1045 (1992).

[LEJAEGHERE2016]

K. Lejaeghere et al.. Reproducibility in density functional theory calculations of solids. Science 351, aad3000 (2016).

[KRESSE1996]

G. Kresse and J. Furthmüller. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B 54, 11169 (1996).

[MONKHORST1976]

H. J. Monkhorst and J. D. Pack. Special points for Brillouin-zone integrations. Phys. Rev. B 13, 5188 (1976).

[SIESTA_GRID]

The real-space grid – convergence of MeshCutoff. SIESTA Documentation (accessed 2025).

[METHFESSEL1989]

M. Methfessel and A. T. Paxton. High-precision sampling for Brillouin-zone integration in metals. Phys. Rev. B 40, 3616 (1989).

[MARZARI1999]

N. Marzari, D. Vanderbilt, A. De Vita, and M. C. Payne. Thermal Contraction and Disordering of the Al(110) Surface. Phys. Rev. Lett. 82, 3296 (1999).