In this work, we focus on an implementation involving PW basis sets, whose particular advantage in MD simulations is that they significantly simplify the calculation of nuclear gradients. This is due to the fact that they are defined without any explicit dependence on nuclear positions, which leads to vanishing Pulay forces.