The chemical reactivity of radical pairs is strongly influenced by the interactions of electronic and nuclear spins. A detailed understanding of these effects requires a quantum description of the spin dynamics that considers spin-dependent reaction rates, interactions with external magnetic fields, spin-spin interactions, and the loss of spin coherence caused by coupling to a fluctuating environment. Modeling real chemical and biochemical systems, which frequently involve radicals with multinuclear spin systems, poses a severe computational challenge. Here, we implement a method based on the stochastic Schrödinger equation in the software package MolSpin. Large electron-nuclear spin systems can be simulated efficiently, with asymmetric spin-selective recombination reactions, anisotropic hyperfine interactions, and nonzero exchange and dipolar couplings. Spin-relaxation can be modeled using the stochastic time-dependence of spin interactions determined by molecular dynamics and quantum chemical calculations or by allowing rate coefficients to be explicitly time-dependent. The flexibility afforded by this approach opens new avenues for exploring the effects of complex molecular motions on the spin dynamics of chemical transformations.