The design of scalable parallel simulation codes for complex phenomena is challenging. For simulations that rely on PDE solution on complex 3D geometries, all the components ranging from the initialization phases to the ultimate linear system solution must be efficiently parallelized. In this paper we address the solution of 3D harmonic Maxwell's equations using a parallel geometric full multigrid scheme where only a coarse mesh, that resolves the object geometry, has to be supplied by an external mesh generator. The electromagnetic problem is solved on a refined mesh that complies with the discretization constraint that links the mesh size and the wavelength. The mesh hierarchy is built in parallel using an automatic mesh refinement technique and parallel matrix-free calculations are performed to process all the finer meshes in the multigrid hierarchy. Only the linear system on the coarse mesh is assembled and solved by a parallel sparse direct solver. We illustrate the numerical robustness of the proposed scheme on large 3D complex geometries and assess the parallel scalability of the implementation on large problems with up to 1.3 billion unknowns on a few hundred cores.