Python script for convergence test

Question

I want to do a convergence test of Ni using a Python script in Pyiron. Can someone please suggest how to write a script to test energy and k-mesh convergence?

Answer 1

There is an example in the pyiron documentation to calculate energy volume curves: https://pyiron.readthedocs.io/en/latest/source/notebooks/energy_volume_curve.html

And the corresponding jupyter notebook is available in the pyiron repository: https://github.com/pyiron/pyiron_atomistics/blob/master/notebooks/energy_volume_curve.ipynb

To adjust the energy cutoff and kpoint mesh, you can write an additional loop and update them using job.set_encut(300.0) for an energy cutoff of 300.0eV and job.set_kpoints([3,3,3]) for a 3x3x3 kpoint mesh. Finally while the example is written for GPAW the same code works for VASP by replacing pr.job_type.Gpaw with pr.job_type.VASP.

Going one step further, there is even a youtube tutorial for calculating an energy volume curve: https://www.youtube.com/watch?v=ZCVRLS4b5YU

Best,

Jan

Answer 2

Here is the simple bash script for testing kpoint energy convergence for VASP you can use same logic in python but you can also use this.

for i in `seq  1 1 5` # change the range needed
do
cat <<EOF >KPOINTS
MONK #header file
0
M
$i $i $i
0 0 0
EOF
mkdir $i
cp INCAR $i/
cp POSCAR $i/
cp POTCAR $i/
cp KPOINTS $i/
cd $i
vasp-5.4.4  # This is vasp run command add mpi if needed
# post process 
E=`grep 'F' OSZICAR|tail -n 1 | awk '{ print $5}'`;                
echo $i $E >>../Result.txt
cd ..
done

keep POSCAR, POTCAR, INCAR and this script in same folder
Execute bash script.sh check Result.txt for output

Answer 3

This is the script I wrote for convergence testing in Quantum ESPRESSO using Bash and Python.

According to Stefaan Cottenier, hydrostatic pressure on a unit cell is the property tested for convergence because it is quite sensitive to numerical precision. If this property is converged with respect to the k-mesh, many other less sensitive properties will be converged, too. Remember that you have to use the tstress keyword in your QE input.

The script generates a new QE input file with a new k-mesh when the convergence test fails. Run conv-test.sh to begin the convergence test.

(1) The Bash code: conv-test.sh

#!/bin/bash

# Generate a text file record.txt to be used for pressure bookkeeping
touch record.txt

# Run the QE input and pass the iteration number to the python script
for i in {1..15}
do
    pw.x -input input-$i.in > output-$i.out
    echo $i | python3 conv-test.py
done

(2) Python code: conv-test.py

# Open "record.txt", read all pressures stored in it from previous runs, and write them to the "stress" list  
record = open("record.txt", "r+")
stress = []
for line in record:
    stress.append(float(line.strip()))
record.close()

# Read the iteration number for the current run from the Bash script
I = input("Iteration number: ")
i = int(I)

# Open the output file with that current iteration number, read the pressure from it, and append it the "stress" list
F1 = open("output-%d.out" %i, "r+")
for line in F1:
    a0 = line.strip()
    a1 = " ".join(a0.split())
    if a1.find("total stress") != -1:
        a2 = a1.replace("total stress (Ry/bohr**3) (kbar) P= ", "")
        stress.append(float(a2))
F1.close()

# Add the value of pressure for this current iteration to "record.txt"
if len(stress) >= 1:
    record = open("record.txt", "a+")
    record.write(str(stress[len(stress)-1])+"\n")
    record.close()

try:
    # If this is the first run and "record.txt" contains a single stress value, increment the k-mesh anyway.
    if len(stress) <= 1:
        F1 = open("input-%d.in" %i, "r+")
        NewContent = ""
        for line in F1:
            NewLine = line.replace("%d %d %d" %(i, i, i), "%d %d %d" %(i+1, i+1, i+1))
            NewContent += NewLine + "\n"
        F1.close()
        F2 = open("input-%d.in" %(i+1), "w+")
        F2.write(NewContent)
        F2.close()

    # If "record.txt" contains at least two values of stress, compare the last two values ans check if they are equal within the required tolerance.
    elif round(stress[len(stress)-1], 2) - round(stress[len(stress)-2], 2) <= 0.001:
        F1 = open("input-%d.in" %i, "r+")
        NewContent = ""
        for line in F1:
            NewLine = line.replace("%d %d %d" %(i, i, i), "%d %d %d" %(i+1, i+1, i+1))
            NewContent += NewLine + "\n"
        F1.close()
        F2 = open("input-%d.in" %(i+1), "w+")
        F2.write(NewContent)
        F2.close()
    # If they are equal within the required tolerance, print "Convergence is achieved."
    else:
        print("Covergence is achieved.")
    print(stress)
except:
    print("Program is terminated.")
    exit()

Word of caution: So far, the script has the following drawbacks:

• It only tests convergence for cubic k-meshes.
• It doesn't terminates properly when convergence is achieved.
• It is designed for static (i.e., scf) calculations. With vc-relax calculations, you will have multiple pressures in the output file and the script will pick the first value while the last value is the one that should be picked. In this case, the script must be modified.

The script is experimental, so all suggestions of enhancement are welcomed.