We present the Gaussian and plane waves (GPW) method and its implementation in Quickstep which is part of the freely available program package CP2K. The GPW method allows for accurate density functional calculations in gas and condensed phases and can be effectively used for molecular dynamics simulations. We show how derivatives of the GPW energy functional, namely ionic forces and the Kohn–Sham matrix, can be computed in a consistent way. The computational cost of computing the total energy and the Kohn–Sham matrix is scaling linearly with the system size, even for condensed phase systems of just a few tens of atoms. The efficiency of the method allows for the use of large Gaussian basis sets for systems up to 3000 atoms, and we illustrate the accuracy of the method for various basis sets in gas and condensed phases. Agreement with basis set free calculations for single molecules and plane wave based calculations in the condensed phase is excellent. Wave function optimisation with the orbital transformation technique leads to good parallel performance, and outperforms traditional diagonalisation methods. Energy conserving Born–Oppenheimer dynamics can be performed, and a highly efficient scheme is obtained using an extrapolation of the density matrix. We illustrate these findings with calculations using commodity PCs as well as supercomputers.