When you run the shell script you should get a series of directories, cutoff_${cutoff}Ry. Run the input files in each directory (you may want to setup a script to do this).

At the end you should have a set of output files that contain the total energy of the system and the forces on each atom.

TASK

Extract and plot the total energy of the system as a function of cutoff

Extract and plot the total force on the system as a function of cutoff (search for 'SUM OF ATOMIC FORCES')

Extract and plot the force on some chosen atoms from the system as a function of cutoff

Analysis

What is converged?

Compare the convergence of forces to the default convergence criteria for geometry optimization.

What sets the required cutoff? It is the basis set (which is dictated by the pseudopotentials). You will need to be able to represent the Gaussian with largest exponent well on the realspace grids. Oxygen, being very electronegative (on the right of the period table with many protons) has very contracted 2s states. You can see in the output

we see that the largest exponent is only 2.7 Bohr-2, so can be represented on a much coarser grid.

Task
If you like, have a look at the BASIS_MOLOPT file (in the data directory, or online here) and see how the exponents change across the periodic table

The convergence is largely dominated by the calculation of the gradient terms in a GGA functional (compare a simulation with LDA to the PBE used here). The evaluation of these terms on the grids are demanding, and very dependent on the functional.