Exact exchange is a primordial ingredient in Kohn–Sham Density Functional Theory based Molecular Dynamics (MD) simulations whenever thermodynamic properties, kinetics, barrier heights or excitation energies have to be predicted with high accuracy. However, the cost of such calculations is often prohibitive, restricting the use of first principles MD to (semi-)local density functionals, in particular in a plane wave basis. We have recently proposed the use of coordinate-scaled orbitals to reduce the cost of the most expensive Fourier transforms during the calculation of the exact exchange potential of isolated systems. Here, we present the implementation and parallelisation of this coordinate scaling approach in the CPMD code and analyse its performance under different parallelisation schemes. We show that speedups increase with system size and that with an optimal configuration, speedups of up to one order of magnitude are possible with respect to conventional calculations. Simulations that have previously taken one week can therefore be finished within less than a day.