EXAMPLE · E3

Automating convergence tests

Automate cutoff, k-point and force convergence tests with shell scripts.

Goal

Learn the standard procedure of scripted convergence tests for ecutwfc, the k-grid, and forces, including the habit of judging in meV/atom and meV/Å. Background in Chapter 05.

New here

Item Role
Deriving inputs with sed One reference input generates the whole scan
grep '^!' Extracts only the converged total energies
meV/atom conversion ΔE × 13605.7 / nat

Input files

The reference input is the same as E1 (si.scf.in).

conv_ecut.sh · conv_kpts.sh · conv_force.sh

#!/bin/bash
# ecutwfc convergence scan; ecutrho follows at 8x (the PAW/US convention).
#
# For each cutoff E:
#   1. sed derives a scan input from si.scf.in, rewriting the two cutoff lines
#   2. pw.x runs it
#   3. grep '^!' picks the line "!    total energy = ..." (the converged value;
#      unmarked "total energy" lines are intermediate SCF iterations)
# Finally awk converts to meV/atom relative to the densest point:
#   dE [meV/atom] = (E - E_ref) * 13605.7 / nat     (1 Ry = 13605.7 meV)

NAT=2
printf "# ecutwfc(Ry)  E_total(Ry)   dE_vs_last(meV/atom)\n" > conv_ecut.dat
LAST=""
for E in 20 25 30 35 40 45 50 60 70 80; do
  sed -e "s/ecutwfc *=.*/ecutwfc      = $E/" \
      -e "s/ecutrho *=.*/ecutrho      = $((E*8))/" si.scf.in > tmp_e$E.in
  pw.x -in tmp_e$E.in > tmp_e$E.out
  EN=$(grep '^!' tmp_e$E.out | tail -1 | awk '{print $5}')
  echo "$E  $EN" >> conv_ecut.dat
  LAST=$EN
done
awk -v nat=$NAT -v ref="$LAST" '!/^#/{printf "%6s  %16s  %10.3f\n",$1,$2,($2-ref)*13605.7/nat}' conv_ecut.dat

The k-point scan (conv_kpts.sh) rewrites the K_POINTS line the same way. The force scan breaks the symmetry first and then reads the total force with head -1, never tail -1:

#!/bin/bash
# Force-based convergence scan (the criterion that matters for ML training
# data). On a perfectly symmetric structure every force is zero, so first
# break the symmetry: move the second Si from 0.25 to 0.26 along x.
sed 's/  Si  0.25  0.25  0.25/  Si  0.26  0.25  0.25/' si.scf.in > si_disp.in
for E in 30 40 50 60 70 80 90; do
  sed -e "s/ecutwfc *=.*/ecutwfc      = $E/" \
      -e "s/ecutrho *=.*/ecutrho      = $((E*8))/" si_disp.in > tmp_f$E.in
  pw.x -in tmp_f$E.in > tmp_f$E.out
  # The FIRST "atom 1 ... force" match after 'Forces acting on atoms' is the
  # total force. Later matches are the contribution breakdown (the last one
  # is the ~1e-6 SCF correction), so head -1 here, never tail -1.
  F=$(grep 'atom    1 type  1   force' tmp_f$E.out | head -1 | awk '{print $7}')
  echo "$E  $F"   # Ry/bohr -> eV/A: multiply by 25.7110
done

Run

bash conv_ecut.sh
bash conv_kpts.sh  > conv_kpts.dat
bash conv_force.sh > conv_force.dat

Output and figure

Si convergence: dE vs ecutwfc, dE vs k-grid, force vs ecutwfc
Measured silicon convergence (QE 7.5, PAW, ecutrho = 8×ecutwfc scanned together). Left: the cutoff scan is monotonic thanks to the variational principle. Center: the k-grid scan, where non-monotonic behavior would also be normal. Right: forces on the distorted structure; here 40 Ry already puts the force within 0.03 meV/Å.

The verdict in numbers (reference: the densest scan point):

ecutwfc ΔE (meV/atom) k-grid ΔE (meV/atom)
20 13.72 2³ 1279
25 5.03 4³ 94.2
30 1.76 6³ 12.0
40 0.91 8³ 1.95
50 0.27 10³ 0.36
60 0.14 12³ 0.073
  • Against a 1 meV/atom criterion, the passing line is around ecutwfc 40 Ry with a 10×10×10 grid. The 30 Ry / 8³ of E1 was a provisional teaching value.
  • The force on the distorted structure (Fx ≈ 0.0288 Ry/bohr = 0.741 eV/Å) moves by 0.5 meV/Å between 30 and 40 Ry, then stays within 0.03 meV/Å. For ML training data, judge with this force criterion (~1 meV/Å).

Exercises

  1. Repeat the scan with ecutrho pinned at 4x and watch what PAW does (the pitfall from Chapter 04).
  2. Compare the cutoff at which the energy converges with the cutoff at which the force converges. Which is higher?
  3. Write your own Python script to plot the scans (the one used for this guide is the plotting script included with the guide's source).
Common mistakes

Extracting forces with grep 'atom 1 ... force' | tail -1. The last match is not the total force; it is the SCF correction term (~10⁻⁶) from the contribution breakdown further down the output. The total force is the first match after Forces acting on atoms. I made exactly this mistake when first running this example; the script above is the corrected version.

05 Cutoff and k-point convergence · 04 Pseudopotentials