01. First SIESTA Calculation — 1D Carbon Chain SCF

As the first hands-on exercise, we perform a self-consistent field (SCF) calculation of a 1D carbon chain (cumulene). The 4-atom cell built in this chapter is the example system that runs through the entire tutorial and becomes the raw material for the electrode and device in later chapters. We go through the fdf input format, structure definition, basis set and grid, k-point sampling, and SCF convergence criteria one by one.

Learning objectives

• Learn the basic syntax of the fdf input format (keyword–value pairs, %block, %include, unit notation)
• Define the lattice, atomic coordinates, and chemical species in fdf
• Understand LCAO basis (DZP) and MeshCutoff as two distinct concepts
• Apply the k-point sampling principle for 1D systems (dense only along the transport axis)
• Read the SCF output (dDmax, Eharris, E_KS, Fermi energy) and judge convergence

Background

The fdf input format

SIESTA's input is a text file in fdf (flexible data format) format. The rules are simple.

• Keyword–value pairs: written as keyword value on one line. Example: MeshCutoff 300. Ry
• Keyword matching is lenient: case-insensitive, and the characters ., -, _ are ignored. MeshCutoff, Mesh.Cutoff, and mesh-cutoff are all the same keyword.
• Physical quantities carry units: specify the unit after the value, as in 300. Ry, 1.29 Ang, 0.01 eV/Ang.
• Blocks: multi-line data are wrapped in %block name ... %endblock name. Lattice vectors, atomic coordinates, and the k-grid are typical examples.
• File splitting: %include filename pulls in another fdf file. The convention of separating the structure into struct.fdf and including it from the main input is used in later chapters (electrode/device). In this chapter we use a single file.
• Comments: everything after # is ignored.
• Logical values: T/F (or true/false).

Example system — the cumulene carbon chain

A one-dimensional chain of carbon atoms with equal spacing $d = 1.29$ Å is called cumulene (all C–C bonds have equal double-bond character). The common specification for this tutorial is as follows.

Item	Value
Atomic spacing	$d = 1.29$ Å (equidistant)
Unit cell	4 C atoms, $c = 4d = 5.16$ Å
Transport axis	$z$
Vacuum	$x = y = 15$ Å

The reason for placing 15 Å of vacuum in the $x$ and $y$ directions is to cut off artificial interactions between neighboring image chains under periodic boundary conditions. Since LCAO basis orbitals are exactly zero beyond the atomic cutoff radius, vacuum just large enough that orbitals do not overlap is sufficient.

LCAO basis and MeshCutoff — two different concepts

SIESTA uses local orbitals attached to atoms (LCAO), not plane waves, as its basis. The basis size is set by PAO.BasisSize, and this tutorial uses DZP (double-zeta polarized) as the standard. DZP consists of 2 radial functions per valence orbital plus polarization orbitals, and is a safe choice for most production calculations.

MeshCutoff, independent of the basis, sets the fineness of the real-space grid on which the electron density and potentials are represented. Its name resembles ENCUT of plane-wave codes and is easily confused with it, but the roles differ — in SIESTA, basis quality is governed by PAO.BasisSize and grid integration accuracy by MeshCutoff. This tutorial uses 300 Ry as the standard.

k-point sampling — for 1D, only the transport axis

This system is periodic only along $z$; $x$ and $y$ are vacuum. Vacuum directions have no band dispersion, so a single k-point (Γ) suffices, and only the $z$ direction is sampled densely. A metallic chain is sensitive to sampling near the Fermi surface (a Fermi point in 1D), so we use a generous $1 \times 1 \times 64$.

SCF convergence criterion — why $10^{-8}$

Convergence of the SCF loop is judged by SCF.DM.Tolerance, the maximum change of the density matrix (DM) between iterations. This tutorial uses 1.0d-8 for all calculations. Even a loose, default-level value ($10^{-4}$) lets the SCF "converge", but the total energy in that state fluctuates in its last digits, so energy differences between structures cannot be trusted. This tutorial deals with fine energy scales — e.g., comparing the energy difference between cumulene and polyyne (at the meV-per-atom level) — so make a tight criterion your habit from the start.

Input file

Create a working directory and save the following as input.fdf. The complete file is also available in the example repository at code/ch01-siesta-scf.

input.fdf

# 01. 1D carbon chain (cumulene) — first SCF
# 4-atom cell, d(C-C) = 1.29 Ang, c = 5.16 Ang, 15 Ang vacuum
# Required file: C.psml (PseudoDojo, PBE, scalar-relativistic)

SystemName   Cumulene C4 chain
SystemLabel  siesta

NumberOfAtoms    4
NumberOfSpecies  1
%block ChemicalSpeciesLabel
  1  6  C
%endblock ChemicalSpeciesLabel

LatticeConstant 1.0 Ang
%block LatticeVectors
  15.000000   0.000000   0.000000
   0.000000  15.000000   0.000000
   0.000000   0.000000   5.160000
%endblock LatticeVectors

AtomicCoordinatesFormat Ang
%block AtomicCoordinatesAndAtomicSpecies
   7.500000   7.500000   0.000000  1
   7.500000   7.500000   1.290000  1
   7.500000   7.500000   2.580000  1
   7.500000   7.500000   3.870000  1
%endblock AtomicCoordinatesAndAtomicSpecies

# --- basis / grid ---
PAO.BasisSize  DZP
MeshCutoff     300. Ry

# --- exchange-correlation ---
XC.Functional  GGA
XC.Authors     PBE

# --- k-points: 1D chain, transport axis z ---
%block kgrid_Monkhorst_Pack
   1   0   0    0.0
   0   1   0    0.0
   0   0  64    0.0
%endblock kgrid_Monkhorst_Pack

# --- SCF ---
MaxSCFIterations   300
SCF.DM.Tolerance   1.0d-8
SCF.Mixer.Weight   0.2
SCF.Mixer.History  15
DM.UseSaveDM       T

# --- output ---
WriteForces        T

Line-by-line commentary

System definition

• SystemName — a human-readable description string. It is echoed in the output and has no effect on the calculation.
• SystemLabel — the prefix of all output files. Here it is siesta, so files such as siesta.DM, siesta.XV, and siesta.EIG are generated.
• %block ChemicalSpeciesLabel — chemical species definition. Each line is species index, atomic number, label. The label C is tied to the pseudopotential file name (C.psml).

Structure definition

• LatticeConstant — the scale factor for the lattice vectors. Setting it to 1.0 Ang makes the numbers in LatticeVectors directly Å.
• %block LatticeVectors — the 3 lattice vectors, one per row. Vacuum $x = y = 15$ Å, and $c = 5.16$ Å along $z$.
• AtomicCoordinatesFormat Ang — declares atomic coordinates in Å-unit Cartesian. Fractional (fractional coordinates with respect to the lattice vectors) is also common.
• %block AtomicCoordinatesAndAtomicSpecies — each line is x y z species-index. The 4 atoms are placed equidistantly at $z = 0, 1.29, 2.58, 3.87$ Å, centered in the cell ($x = y = 7.5$ Å).

Basis, grid, XC, k-points — exactly the values explained in the background section. The k-grid block consists of a $3 \times 3$ integer matrix plus a shift for each row (last column); the diagonal entries are the divisions along each direction.

SCF settings

• MaxSCFIterations 300 — the maximum number of SCF iterations. If this is reached before convergence, the calculation stops, so give it headroom.
• SCF.DM.Tolerance 1.0d-8 — the convergence criterion (see background section).
• SCF.Mixer.Weight 0.2, SCF.Mixer.History 15 — density mixing (Pulay mixing) parameters. Mix in 20% of the new density and use the history of the previous 15 steps. If convergence oscillates, reducing the weight is the first remedy.
• DM.UseSaveDM T — if siesta.DM exists in the directory, reuse it as the initial density. Repeated calculations of the same system get a large SCF speedup (the basic technique of restarting).

SystemLabel determines the output file names

The prefix of every output file is the value of SystemLabel. If this value is scrambled from calculation to calculation, post-processing scripts repeatedly fail to find files. This tutorial fixes the convention — standalone calculations use siesta, electrode calculations use Electrode (Chapter 05), and device calculations use trans (Chapter 06 onward). All TBtrans and sisl post-processing in later chapters assumes this convention.

Placing the pseudopotential

Copy the C.psml downloaded in the installation chapter into the working directory (the same location as input.fdf). SIESTA looks in the working directory for a pseudopotential file with the same name as the label in ChemicalSpeciesLabel. If the file is missing, the run terminates immediately with an error, so verify that these two files exist before running.

ls
# C.psml  input.fdf

Running

Serial execution:

siesta < input.fdf > siesta.out

With an MPI-parallel build:

mpirun -np 4 siesta < input.fdf > siesta.out

A calculation of this size (4 atoms, DZP) finishes within minutes on a laptop. You can watch the progress in another terminal with tail -f siesta.out.

Analyzing the output

siesta.out is organized roughly as header (version and input echo) → basis generation log → SCF cycle → final energies and forces.

The SCF cycle

In the middle of the file, find a table of the following form (the numbers and column layout vary slightly with version and build — below is an example excerpt).

        iscf     Eharris(eV)        E_KS(eV)     FreeEng(eV)     dDmax    Ef(eV)  dHmax(eV)
   scf:    1     -621.383914     -598.130561     -598.152804  1.240718 -8.309820  4.879274
   scf:    2     -600.560423     -611.238379     -611.260622  0.421559 -5.410181  1.680392
   ...
   scf:   34     -616.283719     -616.283719     -616.305962  0.000000 -4.591428  0.000001

SCF cycle converged after 34 iterations

Meaning of each column:

• Eharris — the Harris functional energy. During convergence it approaches along a different path from E_KS, and the two coincide when converged. The difference between the two values is an intuitive gauge of SCF convergence.
• E_KS — the Kohn–Sham total energy.
• FreeEng — the free energy including the entropy term from electronic smearing.
• dDmax — the maximum change of density matrix elements between iterations. Convergence is declared when this drops below SCF.DM.Tolerance ($10^{-8}$).
• Ef — the Fermi energy (eV) at that iteration.
• dHmax — the maximum change of Hamiltonian elements. Used as an auxiliary convergence indicator alongside the DM criterion.

Two things to check for convergence — ① is the SCF cycle converged message present, and ② is the dDmax of the last iteration at or below the tolerance you set? A calculation that stopped by hitting MaxSCFIterations has not converged, and its results must not be used.

Total energy and Fermi energy

In the final summary toward the end of the file:

siesta: E_KS(eV) =             -616.2837

siesta: Final energy (eV):
...
siesta:         Total =    -616.283719
siesta:         Fermi =      -4.591428

• The total energy (Total) is referenced to the pseudopotential, so its absolute value is meaningless by itself; only differences between calculations under the same conditions (basis, grid, k-points) carry physical meaning.
• The Fermi energy is used throughout later chapters as the energy reference ($E - E_F$) for bands and transmission.

Generated files

File	Content
siesta.DM	Converged density matrix — reused as the initial density on recalculation (DM.UseSaveDM)
siesta.XV	Final structure (coordinates + velocities)
siesta.EIG	Eigenvalues at the sampled k-points — the first line is the Fermi energy
siesta.FA	Per-atom forces (WriteForces T)

In equidistant cumulene all atoms are symmetry-equivalent, so it is normal for the forces in siesta.FA to come out close to zero.

Exercises

1. k-point convergence test — vary $k_z$ over 16, 32, 64, 128 and tabulate the total energy. Find the point at which the energy difference between two adjacent settings drops below 1 meV per atom. Explain, from the viewpoint of Fermi-point sampling, why k-convergence is slow in a metallic 1D system.
2. MeshCutoff convergence test — vary MeshCutoff over 200, 250, 300, 400 Ry and check the change in total energy. When the grid becomes coarse, the eggbox effect appears — the energy fluctuates with atomic position. Confirm numerically that 300 Ry is a reasonable compromise.
3. Basis size comparison — change PAO.BasisSize to SZ, DZ, DZP and compare the total energy and Fermi energy. Think about why the energy decreases as the basis grows (the variational principle).

In the next chapter, we add a band path and a PDOS block to this calculation and verify why cumulene is a metal — Chapter 02 — Band Structure and DOS.