Skip to main content

Convergence testing

Convergence with cutoff energy using PWTK

We can automate the previous self consistent calculation by varying a certain parameter. Say we want to check the total energy of the system for various values of ecutwfc. We can do that by using pwtk script.

# load the pw.x input from file

# open a file for writing resulting total energies
set fid [open etot_vs_ecutwfc.dat w]

# loop over different "ecut" values
foreach ecut { 12 16 20 24 28 32 } {

# name of I/O files: $ & $name.out
set name si_scf_ecutwfc-$ecut

# set the pw.x "ecutwfc" variable
SYSTEM "ecutwfc = $ecut"

# run the pw.x calculation
runPW $

# extract the "total energy" and write it to file
set Etot [::pwtk::pwo::totene $name.out]
puts $fid "$ecut $Etot"

close $fid

To run the above script:

pwtk si_scf_ecutoff.pwtk

Now we can plot the total energy with respect to ecutwfc. The data is in etot-vs-ecutwfc.dat

We will use matplotlib to make the plots. Here is the python code for plotting:

import matplotlib.pyplot as plt
from matplotlib import rcParamsDefault
import numpy as np
%matplotlib inline

x, y = np.loadtxt('../src/silicon/etot-vs-ecutwfc.dat', delimiter=' ', unpack=True)
plt.plot(x, y, "o-", markersize=5, label='Etot vs ecutwfc')
plt.xlabel('ecutwfc (Ry)')
plt.ylabel('Etot (Ry)')

Convergence test using UNIX shell script

We can do the convergence test with various parameters. We can calculate the total energy of the system by varying various parameters. We will use the shell script to automate the process with different cutoff energy values.


for CUTOFF in 10 15 20 25 30 35 40
cat > ${NAME}_${CUTOFF}.in << EOF
calculation = 'scf',
prefix = 'silicon'
outdir = './tmp/'
pseudo_dir = './pseudos/'
ibrav = 2,
celldm(1) = 10.0,
nat = 2,
ntyp = 1,
ecutwfc = $CUTOFF
mixing_beta = 0.6

Si 28.086 Si.pz-vbc.UPF

Si 0.0 0.0 0.0
Si 0.25 0.25 0.25

K_POINTS (automatic)
6 6 6 1 1 1

pw.x < ${NAME}_${CUTOFF}.in > ${NAME}_${CUTOFF}.out
echo ${NAME}_${CUTOFF}
grep ! ${NAME}_${CUTOFF}.out


Make sure the file has executable permission for the user:

chmod 700

Run the script file:

# or

We can plot the energy vs cutoff energy, and choose a reasonable value.


Initially, I had problem in running the script in macOS. The problem occurred because the script file format was set to DOS. The file format can be checked in following way:

Open the file in vi editor. vi Now in vi editor command mode (ESC key), type :set ff? This would tell you the file format. Now to change file format, use the command :set fileformat=unix

Convergence test against the number of k-points

We can run similar convergence test against another parameter, and choose the best value of that particular parameter. Here we will try to calculate the number of k-points in the Monkhorst-Pack mesh.


set fid [open etot-vs-kpoint.dat w]

foreach k { 2 4 6 8 } {

set name si_scf_kpoints-$k

K_POINTS automatic "$k $k $k 1 1 1"
runPW $

set Etot [::pwtk::pwo::totene $name.out]
puts $fid "$k $Etot"

close $fid

Run pwtk program:

pwtk si_scf_kpoints.pwtk
x, y = np.loadtxt('../src/silicon/etot-vs-kpoint.dat', delimiter=' ', unpack=True)
plt.plot(x, y, "o-", markersize=5, label='Etot vs kpoints')
plt.xlabel('# kpoints')
plt.ylabel('Etot (Ry)')

Convergence against lattice constant

Calculating total energy with respect to varying lattice constant.


# please uncomment & insert value as determined in the "ecutwfc" exercise
SYSTEM { ecutwfc = 30 }

# please uncomment & insert values as determined in the "kpoints" exercise
K_POINTS automatic { 6 6 6 1 1 1 }

set fid [open etot-vs-alat.dat w]

foreach alat { 9.7 9.8 9.9 10.0 10.1 10.2 10.3 10.4 10.5 10.6 10.7 } {

set name si_scf_alat-$alat

SYSTEM "celldm(1) = $alat"
runPW $

set Etot [::pwtk::pwo::totene $name.out]
puts $fid "$alat $Etot"

close $fid

Run the above code:

pwtk si_scf_alat.pwtk
x, y = np.loadtxt('../src/silicon/etot-vs-alat.dat', delimiter=' ', unpack=True)
plt.plot(x, y, "o-", markersize=5, label='Etot vs alat')
plt.xlabel('alat (Bohr)')
plt.ylabel('Etot (Ry)')

Note on CPU time

  • CPU time is proportional to the number of plane waves used for the calculation. Number of plane wave is proportional to the (ecutwfc)3/2

  • CPU time is proportional to the number if inequivalent k-points

  • CPU time increases as N3, where N is the number of atoms in the system.