We determine the magnetic susceptibility of thermal QCD matter by means of first principles lattice simulations using staggered quarks with physical masses. A novel method is employed that only requires simulations at zero background field, thereby circumventing problems related to magnetic flux quantization. After a careful continuum limit extrapolation, diamagnetic behavior (negative susceptibility) is found at low temperatures and strong paramagnetism (positive susceptibility) at high temperatures. We revisit the decomposition of the magnetic susceptibility into spin- and orbital angular momentum- related contributions. The spin term - related to the normalization of the photon lightcone distribution amplitude at zero temperature - is calculated non-perturbatively and extrapolated to the continuum limit. Having access to both the full magnetic susceptibility and the spin term, we calculate the orbital angular momentum contribution for the first time. The results reveal the opposite of what might be expected based on a free fermion picture. We provide a simple parametrization of the temperature- and magnetic field-dependence of the QCD equation of state that can be used in phenomenological studies.