The ability to quantitatively predict and analyze the rate of electron spin relaxation of open-shell systems is important for electron paramagnetic resonance and paramagnetic nuclear magnetic resonance spectroscopies. We present a combined molecular dynamics (MD), quantum chemistry (QC), and spin dynamics simulation method for calculating such spin relaxation rates. The method is based on the sampling of a MD trajectory by QC calculations, to produce instantaneous parameters of the spin Hamiltonian used, in turn, to numerically solve the Liouville-von Neumann equation for the time evolution of the spin density matrix. We demonstrate the approach by simulating the relaxation of electron spin in an aqueous solution of Ni(2+) ion. The spin-lattice (T1) and spin-spin (T2) relaxation rates are extracted directly from the simulations of the time dependence of the longitudinal and transverse magnetization, respectively. Good agreement with the available, indirectly obtained experimental data is obtained by our method.