# 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.

src/silicon/si_scf_ecutoff.pwtk
# load the pw.x input from fileload_fromPWI pw.scf.silicon.in# open a file for writing resulting total energiesset fid [open etot_vs_ecutwfc.dat w]# loop over different "ecut" valuesforeach ecut { 12 16 20 24 28 32 } {    # name of I/O files: $name.in &$name.out    set name si_scf_ecutwfc-$ecut # set the pw.x "ecutwfc" variable SYSTEM "ecutwfc =$ecut"    # run the pw.x calculation    runPW $name.in # 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:

notebooks/silicon-scf.ipynb
import matplotlib.pyplot as pltfrom matplotlib import rcParamsDefaultimport numpy as np%matplotlib inlineplt.rcParams["figure.dpi"]=150plt.rcParams["figure.facecolor"]="white"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)')plt.legend(frameon=False)plt.show()

## 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.

src/silicon/si_script.sh
#!/bin/shNAME="ecut"for CUTOFF in  10 15 20 25 30 35 40docat > ${NAME}_${CUTOFF}.in << EOF &control    calculation = 'scf',    prefix = 'silicon'    outdir = './tmp/'    pseudo_dir = './pseudos/' / &system    ibrav =  2,    celldm(1) = 10.0,    nat =  2,    ntyp = 1,    ecutwfc = $CUTOFF / &electrons mixing_beta = 0.6 /ATOMIC_SPECIES Si 28.086 Si.pz-vbc.UPFATOMIC_POSITIONS (alat) Si 0.0 0.0 0.0 Si 0.25 0.25 0.25K_POINTS (automatic) 6 6 6 1 1 1EOFpw.x <${NAME}_${CUTOFF}.in >${NAME}_${CUTOFF}.outecho${NAME}_${CUTOFF}grep !${NAME}_${CUTOFF}.outdone Make sure the file has executable permission for the user: chmod 700 si_script.sh Run the script file: ./si_script.sh# orsh si_script.sh We can plot the energy vs cutoff energy, and choose a reasonable value. caution 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 si_script.sh 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. src/silicon/si_scf_kpoints.pwtk load_fromPWI pw.scf.silicon.inset 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$name.in    set Etot [::pwtk::pwo::totene $name.out] puts$fid "$k$Etot"}close $fid Run pwtk program: pwtk si_scf_kpoints.pwtk notebooks/silicon-scf.ipynb 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)')plt.legend(frameon=False)plt.show() ## Convergence against lattice constant​ Calculating total energy with respect to varying lattice constant. src/silicon/si_scf_alat.pwtk load_fromPWI pw.scf.silicon.in# please uncomment & insert value as determined in the "ecutwfc" exerciseSYSTEM { ecutwfc = 30 }# please uncomment & insert values as determined in the "kpoints" exerciseK_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$name.in    set Etot [::pwtk::pwo::totene $name.out] puts$fid "$alat$Etot"}close \$fid

Run the above code:

pwtk si_scf_alat.pwtk
notebooks/silicon-scf.ipynb
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)')plt.legend(frameon=False)plt.show()

## 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.