The Random-Phase approximation (RPA) provides an appealing framework for semi-local density functional theory. In its current formulation, it is cost-effective and has a better scaling behaviour compared to other wavefunction based correlation methods. To broaden the application field for RPA, it is necessary to have first order properties available. RPA nuclear gradients allow for structure optimizations and data sampling for machine learning applications. We report on an efficient implementation of RPA nuclear gradients for massively parallel computers. We apply the implementation to two polymorphs of the benzene crystal obtaining very good cohesive and relative energies. Different correction and extrapolation schemes are investigated for further improvement of the results and in order to estimate error bars.