Time-harmonic electromagnetic scattering by inhomogeneous, three-dimensional structures within a free space environment can be described by electric- and magnetic field, volume integral equations involving the free space Green function. A comprehensive set of Galerkin projection formulations (also known as moment methods) for the numerical solution of these equations is presented, together with comparative numerical results. Such formulations are widely used for particle scattering analysis, optical near field calculation, etc. Results are obtained with higher-order, divergence- and curl-conforming basis functions on iso-parametric, tetrahedral meshes. The results demonstrate that all formulations converge with similar accuracy in the case of an analytically-solvable test problem. When modeling flux densities as solution variables, it is argued that solenoidal function spaces should be used, rather than the standard divergence-conforming function spaces; this assertion is supported by the results. Some of the formulations involve solving for curl-conforming fields; such fields can be discretized with fewer unknowns than divergence-conforming ones, implying lower computational costs. Additionally, some formulations yield system matrices which are approximately half-way sparse, meaning that computational costs will be down by a factor of 2 when iterative solvers are employed, which is the case for the widely-used fast methods.