03. Structure Optimization with VASP

The structure that serves as input to transport calculations must be finalized through relaxation. This tutorial adopts as standard the workflow of using the plane-wave code VASP for the structure-determination step and transferring the result to SIESTA for cross-validation. VASP also reappears later as the hybrid-functional reference when correcting DFT level errors (preview of Chapter 11).

This cross-validation workflow is not limited to the carbon chain; it is the standard procedure reused as is for other systems later (e.g., the metal-electrode molecular junction). In summary:

1. Lattice and coordinate transfer — move the lattice vectors and coordinates from CONTCAR into fdf, changing only units and format (snippet below). Check the bond lengths of the converted result against the original.
2. Band-shape comparison — overlay the band dispersion and presence/absence of a gap from the two codes with the same functional.
3. Relative energy ordering comparison — compare the sign and magnitude of the energy difference between structures A and B. Comparing absolute total energies is forbidden — pseudopotentials and energy references differ between codes.

Learning objectives

• Understand VASP's role in the transport workflow (structure determination, reference electronic structure)
• Write the 4 main input files (INCAR/POSCAR/KPOINTS/POTCAR) and assemble the POTCAR
• Run the carbon chain relaxation and read OSZICAR/OUTCAR/CONTCAR
• Apply the force criterion ($|F| < 0.01$ eV/Å) for convergence
• Convert the relaxed structure to SIESTA fdf format and cross-validate the two codes

Background

Why VASP for structure optimization

Structure optimization is possible with SIESTA too. The reason this tutorial nevertheless uses VASP is methodological independence.

• A plane-wave basis can be converged systematically with a single parameter, ENCUT, and since the basis is not attached to atoms, it is free of the basis-optimization issues that require care in LCAO.
• If the structure and lattice constant computed with two codes using different bases agree, one gains confidence that the result is not an artifact of either basis. As long as the same PBE functional is used, the structures from the two codes should nearly coincide.
• Since VASP also handles the hybrid-functional calculations in later chapters (the reference for level corrections), using the same code from the structure stage onward keeps things consistent.

The 4 main input files

VASP reads 4 files with fixed names from the working directory.

File	Content
INCAR	Calculation parameters (tag = value)
POSCAR	Lattice vectors and atomic coordinates
KPOINTS	k-point sampling
POTCAR	PAW pseudopotential (licensed asset)

Assembling the POTCAR

The POTCAR is made by concatenating the per-element PAW files in exactly the element order of the POSCAR. Denoting the potpaw database path provided by the licensed institution as POTCAR_DIR:

cat $POTCAR_DIR/C/POTCAR > POTCAR

This example has only carbon, so one line suffices, but in a multi-element system, if the order in POSCAR line 6 (the element list) and the cat order disagree, the calculation silently proceeds with the wrong pseudopotential attached to the atoms. Make it a habit to always verify the element order with grep TITEL POTCAR after assembly.

POTCAR must not be distributed

The POTCAR is an asset bound to the VASP license. The example repository of this site does not include POTCAR, and uploading POTCAR to a public repository is a license violation. Assemble it yourself from your institution's potpaw database.

Input files

The example is the same 4-atom carbon chain cell as Chapters 01 and 02. The complete files are in the example repository at code/ch03-vasp.

One important design choice is built into the initial structure. In perfectly equidistant cumulene all atoms are symmetry-equivalent, so the forces are exactly zero, and a relaxation may terminate at the very first step. This is a symmetry-preserving stationary point that numerically cannot explore the unstable mode along the Peierls distortion. Therefore, a small perturbation along the Peierls distortion direction (bond lengths ±0.02 Å) is applied to the initial structure so that the relaxation can explore the lower-energy BLA structure.

POSCAR

C4 chain (cumulene + 0.02 Ang BLA perturbation)
1.0
   15.0000000000    0.0000000000    0.0000000000
    0.0000000000   15.0000000000    0.0000000000
    0.0000000000    0.0000000000    5.1600000000
C
4
Direct
  0.5000000000  0.5000000000  0.0000000000
  0.5000000000  0.5000000000  0.2461240310
  0.5000000000  0.5000000000  0.5000000000
  0.5000000000  0.5000000000  0.7461240310

• Line 1 is a comment; line 2 is the global scale factor.
• Lines 3–5 are the lattice vectors — the same cell as the fdf of Chapter 01 ($15 \times 15 \times 5.16$ Å).
• Lines 6–7 are the element list and counts. This order determines the POTCAR assembly order.
• Direct means fractional coordinates. The $z$ fractions 0.2461 and 0.7461 give an arrangement with bond lengths alternating at 1.27/1.31 Å (a ±0.02 Å perturbation from the equidistant 1.29 Å).

INCAR

SYSTEM = C4 chain relaxation

PREC   = Accurate
ENCUT  = 500
LREAL  = .FALSE.

ISMEAR = 0
SIGMA  = 0.05

EDIFF  = 1E-8
NELM   = 200

IBRION = 2
ISIF   = 2
NSW    = 100
EDIFFG = -0.01

INCAR commentary

Tag	Value	Meaning
PREC	Accurate	Precision preset for the FFT grid etc. — the foundation of good force calculations
ENCUT	500	Plane-wave cutoff (eV). The carbon PAW potential is hard, so take it generously
LREAL	.FALSE.	For small cells, reciprocal-space projection is more accurate
ISMEAR	0	Gaussian smearing. A safe choice for a system whose metallic/semiconducting character may change during the relaxation
SIGMA	0.05	Smearing width (eV). Kept small to reduce energy contamination from smearing
EDIFF	1E-8	Electronic SCF convergence criterion (eV). Force quality is directly tied to SCF convergence, so keep it tight
NELM	200	Maximum number of SCF iterations
IBRION	2	Conjugate-gradient ionic optimization
ISIF	2	Relax atomic positions only; the cell is fixed
NSW	100	Maximum number of ionic steps
EDIFFG	-0.01	Negative = force criterion. Stop when all atomic forces satisfy $

KPOINTS

Gamma-centered 1x1x32 (transport axis z only)
0
Gamma
 1 1 32
 0 0 0

The k-point principle is the same as with SIESTA — 1 along the vacuum directions ($x$, $y$), dense only along the periodic direction ($z$). For the 4-atom cell, 32 corresponds to 128 in primitive-cell terms, which is ample sampling.

If you also want to optimize the lattice constant — the ISIF 3 pitfall

ISIF = 3 relaxes the atoms and the cell together, but it is hard to use directly for this system. Since the stress acts along all three directions, the lattice along the vacuum directions ($x$, $y$) changes too, and the vacuum size is not a physically meaningful degree of freedom to optimize. Selective lattice relaxation of only chosen axes is not a stock feature of standard VASP. Hence there are two ways to determine the lattice constant $c$ of the 1D chain.

1. Manual $c$ scan (recommended) — take several points of $c$ in the range 5.00–5.30 Å, run an ISIF = 2 relaxation at each, and find the minimum of the $E(c)$ curve (Exercise 1).
2. Leave the cell stress as is, check the $z$-direction stress of the ISIF = 2 result, and use it only as a hint for the direction of the $c$ correction.

Also, in calculations where the cell changes, the basis set is defined relative to the cell, so a Pulay stress error appears — the standard prescription when doing cell optimization is to raise ENCUT further.

Running

Run once all 4 files including the POTCAR are in place. The executable name (vasp_std, etc.) and parallel settings depend on the institutional build.

mpirun -np 4 vasp_std > vasp.out

Redirect the standard output to vasp.out to keep a run log — unlike SIESTA, VASP does not automatically save its log to a file, so without redirection the run record is lost.

Analyzing the output

OSZICAR — energy per ionic step

       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.315293846529E+02   -0.31529E+02   ...
...
   1 F= -.31581462E+02 E0= -.31581523E+02  d E =-.315815E+02
DAV:   1    ...
...
   2 F= -.31581891E+02 E0= -.31581952E+02  d E =-.428999E-02
   ...

• The DAV: lines are electronic SCF iterations; an F= line prefixed with a step number marks the completion of one ionic step.
• F is the free energy, and E0 is the energy extrapolated to zero smearing. When using ISMEAR = 0, energy comparisons between structures are done with E0.
• Watch whether d E shrinks toward zero with each ionic step.

OUTCAR — forces and the convergence message

The per-atom forces at each ionic step are in the TOTAL-FORCE table of the OUTCAR.

grep -A 8 "TOTAL-FORCE" OUTCAR | tail -9

 POSITION                                       TOTAL-FORCE (eV/Angst)
 -----------------------------------------------------------------------------------
      7.50000      7.50000      0.00000       0.000000    0.000000    0.004213
      ...

Judge convergence in two ways — ① are all force components of the last step smaller than 0.01 eV/Å, and ② does vasp.out (or OUTCAR) contain the reached required accuracy message? A calculation that stopped by exhausting NSW has not converged.

CONTCAR — the final structure

The CONTCAR is the final (or last-step) structure written in the same format as the POSCAR. What to verify physically in this example is the final BLA — computing the bond lengths from the fractional coordinates, PBE leaves a weak bond alternation (PBE's underestimation of BLA is a known tendency, and hybrid functionals give larger BLA — we meet this again in Chapter 11).

Transferring the structure to SIESTA

Take the relaxed structure into the transport workflow (SIESTA). Unit and format conversion is the crux — CONTCAR uses fractional (Direct) coordinates multiplied by a scale factor, while our fdf convention is Cartesian coordinates in Å.

import numpy as np

with open("CONTCAR") as f:
    lines = f.readlines()

scale = float(lines[1].split()[0])
cell = scale * np.array(
    [[float(x) for x in lines[i].split()[:3]] for i in (2, 3, 4)]
)
symbols = lines[5].split()
counts = [int(x) for x in lines[6].split()]
natoms = sum(counts)

i = 7
if lines[i].strip().lower().startswith("s"):      # skip the Selective dynamics line
    i += 1
cartesian = lines[i].strip().lower()[0] in "ck"   # distinguish Cartesian/Direct
pos = np.array(
    [[float(x) for x in lines[i + 1 + n].split()[:3]] for n in range(natoms)]
)
if not cartesian:
    pos = pos @ cell                              # fractional -> Angstrom Cartesian

species = []
for s, c in zip(symbols, counts):
    species += [s] * c
spmap = {s: j + 1 for j, s in enumerate(dict.fromkeys(species))}

print("LatticeConstant 1.0 Ang")
print("%block LatticeVectors")
for v in cell:
    print(f"  {v[0]:14.8f} {v[1]:14.8f} {v[2]:14.8f}")
print("%endblock LatticeVectors")
print()
print("AtomicCoordinatesFormat Ang")
print("%block AtomicCoordinatesAndAtomicSpecies")
for r, s in zip(pos, species):
    print(f"  {r[0]:14.8f} {r[1]:14.8f} {r[2]:14.8f}  {spmap[s]}")
print("%endblock AtomicCoordinatesAndAtomicSpecies")

Replace the corresponding blocks in the fdf with the printed blocks. The same conversion is possible with sisl (read with sisl.get_sile("CONTCAR").read_geometry() and save as fdf). Whichever route you take, do not skip the verification of checking the bond lengths of the converted result against the original CONTCAR.

Cross-validation

Rerun the SIESTA SCF (Chapter 01) and bands (Chapter 02) with the transferred structure and compare the two codes.

• What may be compared — equilibrium lattice constant and bond lengths, band shape (dispersion, presence/absence of a gap), relative energies (the energy difference between structures A and B). With the same PBE, the results should be close despite the basis difference (LCAO vs plane-wave); a large discrepancy points to the convergence parameters of one side (basis, ENCUT, k-points).
• What must not be compared — absolute total energies. Pseudopotentials and energy references differ between codes, making absolute comparisons meaningless.

Exercises

1. Lattice constant $c$ scan — with the equidistant structure, vary $c$ over 5.00, 5.08, 5.16, 5.24, 5.32 Å (in each POSCAR keep lines 3–5 and the fractional coordinates, modifying only $c$), obtain the single-point energy E0, fit a parabola to $E(c)$, and find the equilibrium $c$. Compare with the tutorial's common value of 5.16 Å.
2. Initial-perturbation dependence — run relaxations with initial BLA perturbations of 0.00, 0.02, 0.05 Å and compare the final BLA. Explain why the relaxation stays exactly at the equidistant structure with a 0.00 perturbation (the symmetry argument), and check whether the final BLA of the perturbed cases converges to the same value regardless of the perturbation size.
3. Cross-validation practice — transfer the relaxed CONTCAR to fdf with the snippet above, compute the SIESTA bands, and overlay them with the equidistant cumulene bands of Chapter 02. Check whether a gap has opened at the Fermi crossing near Γ due to the BLA, and how large the gap is.

With the structure finalized, we are ready to move on to transport. Starting from the next chapter, the NEGF formalism and the electrode calculation begin.