We consider a molecule in solution allowing for more than one thermally accessible electronic state, as is common for paramagnetic molecules. However, the theory is fully general and is also applicable to diamagnetic molecules as a special case. Due to interactions with the environment in solution, e.g., collisions with other molecules, the molecule can undergo electronic transitions as well as whole-body rotations. Both processes are assumed to occur on a timescale much faster than the timescale of the nuclear spin dynamics, such that the dynamics of the magnetic nucleus can be described with a Hamiltonian that is a time average of the true Hamiltonian. Since the nuclear Zeeman Hamiltonian is linear in the magnetic field ⃗ B, this effect can be described by averaging ⃗ B over all possible thermally accessible electronic states and molecular orientations, weighted with their probability. We treat the electrons quantum-mechanically and the rotations classically, assuming that the molecule is otherwise rigid. This is similar to how paramagnetic molecules are usually modelled in semi-classical relaxation theory. We make use of two different types of coordinate systems: A laboratory frame, with orthonormal basis vectors ⃗ e X , ⃗ e Y , ⃗ e Z , in which the external field B lab 0 is constant, and a molecule-fixed frame, with orthonormal basis vectors ⃗ e x , ⃗ e y , ⃗ e z , which is rotating together with the molecule such that the coordinates of the atoms expressed in this frame are constant.