We investigate in detail the initial susceptibility, magnetization curves, and microstructure of ferrofluids in various concentration and particle dipole moment ranges by means of molecular dynamics simulations. We use the Ewald summation for the long-range dipolar interactions, take explicitly into account the translational and rotational degrees of freedom, coupled to a Langevin thermostat. When the dipolar interaction energy is comparable with the thermal energy, the simulation results on the magnetization properties agree with the theoretical predictions very well.For stronger dipolar couplings, however, we find systematic deviations from the theoretical curves.We analyze in detail the observed microstructure of the fluids under different conditions. The formation of clusters is found to enhance the magnetization at weak fields and thus leads to a larger initial susceptibility. The influence of the particle aggregation is isolated by studying ferrosolids, which consist of magnetic dipoles frozen in at random locations but which are free to rotate.Due to the artificial suppression of clusters in ferro-solids the observed susceptibility is considerably lowered when compared to ferrofluids.