For many years ab initio electronic structure calculations based upon density functional theory have been one of the main application areas in high performance computing (HPC). Typically, the Kohn-Sham equations are solved by minimisation of the total energy functional, using a plane wave basis set for valence electrons and pseudopotentials to obviate the representation of core states. One of the best known and widely used software for performing this type of calculation is the Vienna Ab initio Simulation Package, VASP, which currently offers a parallelisation strategy based on the distribution of bands and plane wave coefficients over the machine processors. We report here an improved parallelisation strategy that also distributes the k-point sampling workload over different processors, allowing much better scalability for massively parallel computers. As a result, some difficult problems requiring large k-point sampling become tractable in current computing facilities. We showcase three important applications: dielectric function of epitaxially strained indium oxide, solution energies of tetravalent dopants in metallic VO 2, and hydrogen on graphene.