Automating convergence tests

Goal

This is where QE skill forks. 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

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: measured

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

• 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 trap of 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 in the repository at .build/plot_qe.py).

Related chapters

05 Cutoff and k-point convergence · 04 Pseudopotentials