Non-perturbative calculation of orbital and spin effects in molecules subject to non-uniform magnetic fields

External non-uniform magnetic fields acting on molecules induce non-collinear spin densities and spin-symmetry breaking. This necessitates a general two-component Pauli spinor representation. In this paper, we report the implementation of a general Hartree-Fock method, without any spin constraints, for non-perturbative calculations with finite non-uniform fields. London atomic orbitals are used to ensure faster basis convergence as well as invariance under constant gauge shifts of the magnetic vector potential. The implementation has been applied to investigate the joint orbital and spin response to a field gradient-quantified through the anapole moments-of a set of small molecules. The relative contributions of orbital and spin-Zeeman interaction terms have been studied both theoretically and computationally. Spin effects are stronger and show a general paramagnetic behavior for closed shell molecules while orbital effects can have either direction. Basis set convergence and size effects of anapole susceptibility tensors have been reported. The relation of the mixed anapole susceptibility tensor to chirality is also demonstrated.