The simulation of condensed matter in first principles Molecular Dynamics (FPMD) heavily relies on Kohn-Sham Density Functional Theory (KS-DFT) calculations. The accuracy of such simulations is governed by the reliability of the underlying potential energy surface (PES) and the accessible time scale; both factors are crucially influenced by the choice of exchange-correlation (xc) functional. Hybrid xc functionals, which include a fraction of exact exchange, are among the most accurate xc approximations available. However, in particular in combination with the plane wave basis commonly used in FPMD, the computational overhead due to the evaluation of exchange integrals significantly limits the time scale of the simulation, thus hampering convergence of observables. This work aims at improving the availability of hybrid xc functionals in a plane wave basis while reducing their computational overhead. Part I provides an introduction to the field of FPMD. Starting from the quantum mechanical partition function, we show the simplifications behind the single-state Born-Oppenheimer approximation, followed by a discussion of the theoretical foundations of KS-DFT, an overview of practical approximations, and their implementation in a plane wave/pseudopotential formalism. In Part II, we outline the implementation and performance of two popular families of xc functionals. We first derive the plane wave expressions for the Coulomb-Attenuation Method (CAM) and compare its performance in a plane wave basis to common atom-centred basis sets. A new xc functional, CAM-O3LYP, is reported, which outperforms the commonly employed CAM-B3LYP in some notorious cases. We then present the implementation of the M05 to M11 families of Minnesota functionals in a plane wave basis. A detailed convergence analysis suggests that a fine enough integration mesh has to be provided in order to obtain accurate energy differences, and a comparison of reaction enthalpies obtained in plane waves and atom-centred bases reveals an unusual sensitivity of certain Minnesota functionals to the choice of basis. We provide a rationale for the observed differences by comparing electron densities between different functionals and basis sets demonstrating that plane waves make it possible to conveniently obtain values for basis set limits for functionals that can be difficult to converge in atom-centred bases. Part III introduces a coordinate-scaling scheme for the calculation of exact exchange integrals for isolated systems in a plane wave basis. By applying the scaling relations of the exact exchange functional to the Hartree-Fock term, we show that it is possible to effectively reduce the cost of the Fast Fourier Transforms (FFT) that constitute the bottleneck of the calculation while retaining good accuracy. We then present the practical implementation and parallelisation of the method in the CPMD code and show that speedups can reach values of up to one order of magnitude, thereby considerably improving the computational footprint of hybrid functionals. In Part IV, we offer a perspective on further work by discussing recent additions to the library of dispersion-corrected atom-centred potentials (DCACP) and by introducing a method for structural refinement based on electron density difference maps, density-difference driven molecular dynamics (d3MD). We conclude by providing a brief overview of the results obtained in this work and by discussing future perspectives.